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:


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


Drop the first element of bE (essentially because this  will always be 1-- I don't care about b=0 as this is trivial).
Start with result = 1. (i.e. a^0)
If bE is empty, then return result. 
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

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, #] &amp; /@ 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?
TL;DR

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

largepowmod :: Integer -&gt; Integer -&gt; Integer -&gt; 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 -&gt;
      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 &amp;&amp; ./Test +RTS -s -p &amp;&amp; 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 -&gt; Integer -&gt; Integer -&gt; 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 -&gt;
      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 -&gt;
      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.