I'm trying to write a Montgomery exponentiation based on this which can compete with Mathematica PowerMod
. We know that PowerMod
uses square and multiply technique.
The speedup must be obtained by replacing modular operation mod n with modular operation mod 2^x. Can we accomplish this in Mathematica such that it take over from PowerMod
?
This is my implementation:
Global variables:
RLength = 0; R = 0; RM1=0; RInverse = 0; NPrime = 0; verbose = False;
and the MontExp (b^e mod n):
MontExp[b_, e_, n_] := (RLength = BitLength[n]; R = 2^RLength; RM1=R-1;
RInverse = PowerMod[R, -1, n]; NPrime = PowerMod[-n, -1, R];
M = Mod[b*R, n]; Result = Mod[R, n];
If[verbose,
Print["MontParams: R=", R, ", RInverse=", RInverse, " ,NPrime=",
NPrime, " ,M=", M]];
Do[Result = Mont[Result, Result, n];
If[expBit == 1, Result = Mont[Result, M, n]], {expBit,
IntegerDigits[e, 2]}]; Result = Mont[Result, 1, n]; Return[Result])
Mont function version 0:
Mont[u_, v_, n_] := (z = Mod[u*v*RInverse, n];
If[verbose, Print["Monto ", u, " * ", v, " => ", z]]; Return[z]);
Mont function version 1:
Mont[u_, v_, n_] := (t = u*v;
z = BitShiftRight[(t + Mod[t*NPrime, R] n), RLength];
If[verbose, Print["Monto ", u, " * ", v, " => ", z]]; Return[z]);
Mont function version 2:
Mont[u_, v_, n_] := (t = u*v;
z = BitShiftRight[(t + BitAnd[t*NPrime, RM1] n), RLength];
If[verbose, Print["Monto ", u, " * ", v, " => ", z]]; Return[z]);
and Timings:
p = 2^20000 + 1;
Mathematica PowerMod:
Timing[PowerMod[2, p, p] == 2]
{1.529, False}
Mont v0
{3.432, False}
Mont v1
{7.332, False}
Mont v2
{3.541, False}
As you can see as I tried to improve it with binary shifts instead of modular operations, it had negative impact on the speed. That's probably because of non-native implementation in Mathematica. Any idea to improve it?
--
Update
I've learnt that Mod[b, 2^n] == BitAnd[b, 2^n-1]
so I changed the Version 2 to use BitAnd
, but yet no gain in compare to original PowerMod
...
Update 2
It seems that because of its reliance on shifts, the speedups are only for 2^k+1 numbers. However I saw an amazing result from @Simoon-Woods answer:
list = (2^# + 1) & /@ Range[5000, 100000, 5000];
PowerModTimings = First[Timing[PowerMod[2, #, #]]] & /@ list
{0.047,0.265,0.717,1.529,2.871,4.336,6.506,9.173,11.934,16.879,20.623,25.772,30.373,37.3,45.49,55.131,63.274,73.788,85.114,96.112}
MontExpTimingsV2 = First[Timing[MontExp[2, #, #]]] & /@ list (*Mont version 2*)
{0.063,0.156,0.312,0.483,5.258,0.92,4.711,1.56,8.081,18.642,2.949,3.51,16.13,18.268,36.91,5.569,15.413,28.236,106.143,60.388}
MontExpTimingsV0 = First[Timing[MontExp[2, #, #]]] & /@ list (*Mont version 0*)
{0.047,0.093,0.188,0.234,2.418,0.53,2.246,0.858,3.417,8.58,1.451,1.669,7.051,8.003,18.221,2.824,6.91,13.245,51.231,29.047}
And plotting the result:
ListLinePlot[{PowerModTimings, MontExpTimingsV2, MontExpTimingsV0}]
Update 3
I've added timings for Mont
version 0 based on @Simon-Woods answer. Great timings ...
Compile
withCompilationTarget-> "C"
. $\endgroup$ – 0xFE Apr 23 '13 at 14:13Compile
expects machine-sized integers in most places. You'd need to use an alternate representation. $\endgroup$ – 0xFE Apr 23 '13 at 14:42