5
\$\begingroup\$

I've implemented my own 'power-mod' procedure which will compute a^b 'mod' n for some large b and even larger n. I think the way I've done it is a common procedure, and it does seem well optimised, but I really need to squeeze the maximum performance out of it as it is currently the bottleneck in my program (~75% of the CPU time).

First I will explain in words the procedure (ironically, in an imperative style), then display the code (which is quite short).


In Words

Calculating a^b mod n:

  1. Obtain the binary expansion of b in Big Endian form. Call it bE.

For example, b = 10 becomes b = [1,0,1,0] which is the reverse of its conventional binary representation 0101. This part is well optimised (~5% of the CPU time).

  1. Drop the first element of bE (essentially because this will always be 1-- I don't care about b=0 as this is trivial).

  2. Start with result = 1. (i.e. a^0)

  3. If bE is empty, then return result.

  4. Take the next element from bE to be k.

    • If k is 0, then do result = result^2 'mod' n
    • If k is 1, then do result = result^2 * a 'mod' n
  5. Go to step 4.


The Code

-- This part computes the Big Endian binary representation list of an Integer `k`. It's not a bottleneck, so I'm not so concerned about this.
binarizeE k = binarizeETab k []
  where
    binarizeETab 0 xs  = xs
    binarizeETab k xs = binarizeETab (fst kdivmod) ((snd kdivmod):xs)
      where
        kdivmod = k `divMod` 2

-- The part that I want to optimise to oblivion!
-- Doing `rem` and then ending with a final `mod` may make an incredibly minute improvement for extremely large numbers as compared with using `mod` each and every step.
largepowmod num pow modbase = (operate num opList) `mod` modbase
  where
    opList = drop 1 $ binarizeE pow
    operate k    []   = k
    operate k (0:ops) = operate (k^2         `rem` modbase) ops
    operate k (1:ops) = operate ((k^2 * num) `rem` modbase) ops

My Optimistic Goal

The background here is that I've implemented the Baillie–PSW Primality Test and am trying to improve its performance so that it is comparable with professional software like Mathematica. This bottleneck is happening during the Miller-Rabin Test procedure. Since Haskell is compiled, and Mathematica on my machine utilizes the same CPU % as my program, I don't see why I couldn't achieve something close. At the moment, they are indeed close for some tests (1.1s versus 1.4s), but for fun I'd like to see if maybe I can even surpass Mathematica!


Additional Information / Benchmark

The base a is in my case always 2. Maybe this lends itself to some neat trick. The exponent b and modulo base n vary greatly, but a typical example would be something like n an odd number somewhere in the vicinity of 10^1000 and b ~= n/2.

Example test comparison of my current implementation and Mathematica's:

Set b(n) = (n-1)/2. Then evaluating 2^b mod n for every odd n between 10^1000 and 10^1000 + 10000 takes ~110s for my program, and ~92s for Mathematica.

Testing code below for reference.

Haskell:

import PrimeStuff
import PrimeStuff.PQTrials
import System.Environment
import Control.DeepSeq

main :: IO()
main = do
        let trialNs   = [10^1000 + 1, 10^1000 + 3.. 10^1000 + 9999]
        let modTest n = largepowermod 2 ((n-1) `div` 2) n
        let test      = map modTest trialNs

        let sol   = deepseq test "Done."
        putStrLn(sol)

Mathematica:

nList = 10^1000 + (2*Range[5000] - 1);
Timing[results = PowerMod[2, (# - 1)/2, #] & /@ nList;][[1]]

Ideas

I am always using a=2 as the base, so in big endian form, a^b looks like [1,0,0,0.....0] with b zeros. Is there perhaps a way of taking the modulo of this with respect to n by manipulating the 0s and 1s directly? Would such a thing be faster than what the computer/compiler is already doing with my current code?

\$\endgroup\$
  • \$\begingroup\$ Although Haskell is compiled, there is a chance that Mathematica has written the primality test in C (which is usually faster than Haskell). \$\endgroup\$ – Dair Sep 25 '16 at 5:52
  • \$\begingroup\$ Did you also leave the types off in your code? \$\endgroup\$ – Zeta Sep 26 '16 at 9:04
  • \$\begingroup\$ @Zeta - Yes, does this affect performance? \$\endgroup\$ – Myridium Sep 26 '16 at 9:05
  • \$\begingroup\$ @Myridium It may very well, but that depends on the actual number size. You said your numbers are huge, but it's not clear how huge they actually are. Could you add some examples (e.g. main = print $ largepowmod (10^3) (10^6 + 3) (10 ^ 9), the timings you got and comparison timings from Mathematica? \$\endgroup\$ – Zeta Sep 26 '16 at 9:12
  • \$\begingroup\$ @Zeta - Alright, I'll add types (can also change binary expansion to type [Int] instead of [Integer]). The base of the exponent is always 2 actually-- perhaps this lends itself to a trick. The exponent and modular base are both on the order of ~10^100. I shall add the extra info soon. \$\endgroup\$ – Myridium Sep 26 '16 at 9:21
4
\$\begingroup\$

TL;DR

Use strictness annotations at the right point to get rid of thunks:

largepowmod :: Integer -> Integer -> Integer -> Integer
largepowmod num pow modbase = (operate num opList) `mod` modbase
  where
    opList = drop 1 $ binarizeE pow
    operate !k    []   = k
    operate !k (0:ops) = operate (k^2         `rem` modbase) ops
    operate !k (1:ops) = operate ((k^2 * num) `rem` modbase) ops
--          ^

Note that this is equivalent to using foldl'.

Proper benchmarks

First of all, let's get reliable benchmarks. We can use Criterion for this:

import Criterion.Main
import Text.Printf (printf)

main :: IO()
main = defaultMain [
  bgroup "largepowmod" $ 
    flip map ([1,3..5] ++ [9999]) $ \i ->
      bench (printf "10^1000 + %d" i) $   
        let n = 10^1000 + i :: Integer
        in whnf (largepowmod 2 ((n - 1) `div` 2)) n
  ]

On my machine, I get about 16.4ms per call with -O2:

benchmarking largepowmod/10^1000 + 1
time                 16.14 ms   (16.07 ms .. 16.22 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 16.18 ms   (16.13 ms .. 16.23 ms)
std dev              124.6 μs   (99.83 μs .. 164.8 μs)

benchmarking largepowmod/10^1000 + 3
time                 16.44 ms   (16.26 ms .. 16.63 ms)
                     1.000 R²   (0.999 R² .. 1.000 R²)
mean                 16.35 ms   (16.29 ms .. 16.42 ms)
std dev              168.3 μs   (135.2 μs .. 223.2 μs)

benchmarking largepowmod/10^1000 + 5
time                 16.24 ms   (16.19 ms .. 16.30 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 16.35 ms   (16.30 ms .. 16.41 ms)
std dev              138.4 μs   (89.59 μs .. 182.9 μs)

benchmarking largepowmod/10^1000 + 9999
time                 16.45 ms   (16.34 ms .. 16.58 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 16.37 ms   (16.30 ms .. 16.43 ms)
std dev              165.2 μs   (129.3 μs .. 235.4 μs)

But let's not stop here. Let's have a look at the profile:

$ ghc --make Test.hs -O2 -prof -auto-all && ./Test +RTS -s -p && head -n20 Test.prof

        Tue Sep 27 13:32 2016 Time and Allocation Profiling Report  (Final)

           Test +RTS -s -p -RTS --output benchmark.html

        total time  =       19.65 secs   (19652 ticks @ 1000 us, 1 processor)
        total alloc = 8,760,057,056 bytes  (excludes profiling overheads)

COST CENTRE         MODULE                %time %alloc

largepowmod.operate Main                   91.9   83.5
getGCStats          Criterion.Measurement   5.1   10.6
getOverhead         Criterion.Monad         2.2    1.5
applyGCStats        Criterion.Measurement   0.1    1.5


                                                                                                        individual     inherited
COST CENTRE                                  MODULE                                   no.     entries  %time %alloc   %time %alloc

MAIN                                         MAIN                                     273           0    0.0    0.0   100.0  100.0
 main                                        Main                                     550           0    0.0    0.0   100.0   99.9
  main.\                                     Main                                     745           0    0.0    0.0    91.9   83.5
   largepowmod                               Main                                     746        1208    0.0    0.0    91.9   83.5
    largepowmod.operate                      Main                                     752     4011768   91.9   83.5    91.9   83.5
    largepowmod.opList                       Main                                     751        1208    0.0    0.0     0.0    0.0

So most of the time is used in operate. It's allocating most of the data. However, that's a little bit of a red herring, since binarizeE is inlined at this point, and poor operate is merely allocating the list that will be generated lazily.

Getting stricter

However, we're not strict enough. In operate we create a large chunk of intermediate data:

operate k    []   = k
operate k (0:ops) = operate (k^2         `rem` modbase) ops
operate k (1:ops) = operate ((k^2 * num) `rem` modbase) ops
--                           ^^^^^^^^^^^^^^^^^^^^^^^^^

Haskell is lazy, therefore (k^2 * num) `rem` modbase isn't evaluated until it's actually needed. Therefore we allocate \$\mathcal O(\log_2 p)\$ additional terms. This isn't necessary. We can provide a small fix to make sure that k has been fully evaluated:

largepowmod' :: Integer -> Integer -> Integer -> Integer
largepowmod' num pow modbase = (operate num opList) `mod` modbase
  where
    opList = drop 1 $ binarizeE pow
    operate !k    []   = k
    operate !k (0:ops) = operate (k^2         `rem` modbase) ops
    operate !k (1:ops) = operate ((k^2 * num) `rem` modbase) ops

This needs -XBangPatterns. How does this fare? We change our main:

main :: IO()
main = defaultMain [
  bgroup "largepowmod" $ 
    flip map ([1,3..5] ++ [9999]) $ \i ->
      bench (printf "10^1000 + %d" i) $   
        let n = 10^1000 + i :: Integer
        in nf (largepowmod 2 ((n - 1) `div` 2)) n
 , bgroup "largepowmod'" $ 
    flip map ([1,3..5] ++ [9999]) $ \i ->
      bench (printf "10^1000 + %d" i) $   
        let n = 10^1000 + i :: Integer
        in nf (largepowmod' 2 ((n - 1) `div` 2)) n
  ]

And run it again (without -prof):

benchmarking largepowmod/10^1000 + 1
time                 16.14 ms   (16.07 ms .. 16.22 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 16.18 ms   (16.13 ms .. 16.23 ms)
std dev              124.6 μs   (99.83 μs .. 164.8 μs)

benchmarking largepowmod/10^1000 + 3
time                 16.44 ms   (16.26 ms .. 16.63 ms)
                     1.000 R²   (0.999 R² .. 1.000 R²)
mean                 16.35 ms   (16.29 ms .. 16.42 ms)
std dev              168.3 μs   (135.2 μs .. 223.2 μs)

benchmarking largepowmod/10^1000 + 5
time                 16.24 ms   (16.19 ms .. 16.30 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 16.35 ms   (16.30 ms .. 16.41 ms)
std dev              138.4 μs   (89.59 μs .. 182.9 μs)

benchmarking largepowmod/10^1000 + 9999
time                 16.45 ms   (16.34 ms .. 16.58 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 16.37 ms   (16.30 ms .. 16.43 ms)
std dev              165.2 μs   (129.3 μs .. 235.4 μs)

benchmarking largepowmod'/10^1000 + 1
time                 14.54 ms   (14.49 ms .. 14.60 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 14.51 ms   (14.46 ms .. 14.54 ms)
std dev              101.0 μs   (64.28 μs .. 177.1 μs)

benchmarking largepowmod'/10^1000 + 3
time                 14.53 ms   (14.49 ms .. 14.58 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 14.53 ms   (14.50 ms .. 14.56 ms)
std dev              78.32 μs   (62.67 μs .. 95.66 μs)

benchmarking largepowmod'/10^1000 + 5
time                 14.42 ms   (14.36 ms .. 14.48 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 14.56 ms   (14.51 ms .. 14.65 ms)
std dev              180.4 μs   (92.90 μs .. 298.2 μs)

benchmarking largepowmod'/10^1000 + 9999
time                 14.60 ms   (14.49 ms .. 14.85 ms)
                     0.998 R²   (0.993 R² .. 1.000 R²)
mean                 14.61 ms   (14.53 ms .. 14.87 ms)
std dev              314.4 μs   (69.97 μs .. 618.3 μs)

We went from 16.14ms to 14.77ms (0.89 of original time), or if we extrapolate this to your 110s to 98s, which is closer to Mathematica's time.

\$\endgroup\$
  • \$\begingroup\$ I probably should have mentioned that my benchmark was with foldl' implemented as per @Gurkenglas' suggestion. Also, the timing of even a long test like this (~2m) varies wildly on my machine. This is probably because it doesn't have a fixed CPU rate. Consequently, the difference of 16.14->14.77 is well within the fluctuations I get between identical tests! Is there a reliable way of profiling which can eliminate some of this noise? i.e. the .prof file says the number of ticks taken, but this too fluctuates wildly. \$\endgroup\$ – Myridium Sep 27 '16 at 13:08
  • \$\begingroup\$ Nonetheless, I appreciate your helpful sequel to the profiling 101 you gave me on my previous question. (I've been referring to that quite a bit). \$\endgroup\$ – Myridium Sep 27 '16 at 13:09
  • \$\begingroup\$ @Myridium: Criterion uses multiple passes to get reliable information. If the CPU fluctuates very much, Criterion will report outliers and a large std dev (and a non-1.0 R²). \$\endgroup\$ – Zeta Sep 27 '16 at 13:29
  • \$\begingroup\$ Oh! Right, I can just check the mean! The depth of profiling data excites me... I very much appreciate your help; it's because of your previous answer that I was able to implement this Primality Test efficiently in the first place. \$\endgroup\$ – Myridium Sep 27 '16 at 13:31
  • \$\begingroup\$ @Myridium Ah, the prime sieve. I completely forgot about that one. Either way, if you use Criterion, make sure to use --output <filename>.html to get a nice report. For more information, see serpentine.com/criterion. \$\endgroup\$ – Zeta Sep 27 '16 at 13:33
3
\$\begingroup\$

An equivalent rewrite of largepowmod:

largepowmod num pow modbase = foldl foo num opList `mod` modbase
  where
    opList = drop 1 $ binarizeE pow
    foo k 0 = k^2       `rem` modbase
    foo k 1 = k^2 * num `rem` modbase

Let's assume this one doesn't hit a bottleneck:

largepowmod num pow modbase = foldl foo 1 (binarizeE pow) `mod` modbase
  where
    foo k 0 = k^2       `rem` modbase
    foo k 1 = k^2 * num `rem` modbase

Maybe some strictness can help? I've heard that foldl is pretty much never what you want compared to foldl':

largepowmod num pow modbase = foldl' foo 1 (binarizeE pow) `mod` modbase
  where
    foo k 0 = k^2       `rem` modbase
    foo k 1 = k^2 * num `rem` modbase

Another thing that comes to mind is that k*k might be faster than k^2 for obvious reasons. looks it up Oh hey maybe you want to look at http://hackage.haskell.org/package/base-4.9.0.0/docs/src/GHC.Real.html#%5E

\$\endgroup\$
  • \$\begingroup\$ Unfortunately, neither k*k or foldl' make a discernible difference to the computation time. Nonetheless, I'll take them on board. \$\endgroup\$ – Myridium Sep 26 '16 at 0:34
  • \$\begingroup\$ I still don't understand foldl', but the source you provide says it has a special case for small exponents anyway, so I think a^2 would be perfectly equivalent to a*a. \$\endgroup\$ – Myridium Sep 26 '16 at 10:47
  • \$\begingroup\$ Concerning your k^2 remark, you want to look at the documentation, too, right? There is an optimization rule for ^2 ;) \$\endgroup\$ – Zeta Sep 27 '16 at 12:43
  • \$\begingroup\$ Ah, I expected any documentation to be above the relevant sourcecode, not below. ^2 is fine then, I'd even say more idiomatic. \$\endgroup\$ – Gurkenglas Sep 27 '16 at 16:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.