# Large Power Mod - Extreme Efficiency

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? • Although Haskell is compiled, there is a chance that Mathematica has written the primality test in C (which is usually faster than Haskell). – Dair Commented Sep 25, 2016 at 5:52 • Did you also leave the types off in your code? – Zeta Commented Sep 26, 2016 at 9:04 • @Zeta - Yes, does this affect performance? Commented Sep 26, 2016 at 9:05 • @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?
– Zeta
Commented Sep 26, 2016 at 9:12
• @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. Commented Sep 26, 2016 at 9:21

# 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.

• 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. Commented Sep 27, 2016 at 13:08
• 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). Commented Sep 27, 2016 at 13:09
• @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²).
– Zeta
Commented Sep 27, 2016 at 13:29
• 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. Commented Sep 27, 2016 at 13:31
• @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.
– Zeta
Commented Sep 27, 2016 at 13:33

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

• Unfortunately, neither k*k or foldl' make a discernible difference to the computation time. Nonetheless, I'll take them on board. Commented Sep 26, 2016 at 0:34
• 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. Commented Sep 26, 2016 at 10:47
• Concerning your k^2 remark, you want to look at the documentation, too, right? There is an optimization rule for ^2 ;)
– Zeta
Commented Sep 27, 2016 at 12:43
• Ah, I expected any documentation to be above the relevant sourcecode, not below. ^2 is fine then, I'd even say more idiomatic. Commented Sep 27, 2016 at 16:58