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
bin Big Endian form. Call itbE.
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 be1-- I don't care aboutb=0as this is trivial).Start with
result = 1. (i.e.a^0)If
bEis empty, then returnresult.Take the next element from
bEto bek.- If
kis 0, then doresult = result^2 'mod' n - If
kis 1, then doresult = result^2 * a 'mod' n
- If
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?
main = print $ largepowmod (10^3) (10^6 + 3) (10 ^ 9), the timings you got and comparison timings from Mathematica? \$\endgroup\$[Int]instead of[Integer]). The base of the exponent is always2actually-- 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\$