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).
a^b mod n:
- Obtain the binary expansion of
bin Big Endian form. Call it
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=0as this is trivial).
result = 1. (i.e.
bEis empty, then return
Take the next element from
kis 0, then do
result = result^2 'mod' n
kis 1, then do
result = result^2 * a 'mod' n
Go to step 4.
-- 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
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
b ~= n/2.
Example test comparison of my current implementation and Mathematica's:
b(n) = (n-1)/2. Then evaluating
2^b mod n for every odd
10^1000 + 10000 takes ~110s for my program, and ~92s for Mathematica.
Testing code below for reference.
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)
nList = 10^1000 + (2*Range - 1); Timing[results = PowerMod[2, (# - 1)/2, #] & /@ nList;][]
I am always using
a=2 as the base, so in big endian form,
a^b looks like
b zeros. Is there perhaps a way of taking the modulo of this with respect to
n by manipulating the
1s directly? Would such a thing be faster than what the computer/compiler is already doing with my current code?