5
$\begingroup$

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}]

enter image description here

Update 3

I've added timings for Mont version 0 based on @Simon-Woods answer. Great timings ...

$\endgroup$
9
  • $\begingroup$ Look into Compile with CompilationTarget-> "C". $\endgroup$
    – 0xFE
    Commented Apr 23, 2013 at 14:13
  • 1
    $\begingroup$ You have to be careful though because Compile expects machine-sized integers in most places. You'd need to use an alternate representation. $\endgroup$
    – 0xFE
    Commented Apr 23, 2013 at 14:42
  • 1
    $\begingroup$ I tied it and because of the "machine-size" integers it's useless for this scenario. $\endgroup$ Commented Apr 23, 2013 at 15:24
  • $\begingroup$ You could probably do some trick like representing it as a list of integers to get it to compile, but I don't know if that overhead would be more or less than the benefit you would get from compiling. $\endgroup$
    – 0xFE
    Commented Apr 24, 2013 at 2:33
  • 2
    $\begingroup$ I think your questions are good ones, but unfortunately your usual goal of speeding up number-theoretical functions by rewriting them at the top level is in general going to be somewhere between extremely difficult and impossible. Might I suggest that you instead consider doing this with GMP in a MathLink program? $\endgroup$ Commented Apr 24, 2013 at 10:49

1 Answer 1

7
+50
$\begingroup$

For the example problem I get about a factor of 4 speedup over PowerMod by memoizing Mont. This of course means that Mont should not contain any global variables so I rewrote the code slightly:

MontExp[b_, e_, n_] := Module[
  {RLength, R, RM1, RInverse, NPrime, M, Result},
  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];
  Do[
   Result = Mont[Result, Result, n, NPrime, RM1, RLength];
   If[expBit == 1, Result = Mont[Result, M, n, NPrime, RM1, RLength]]
   , {expBit, IntegerDigits[e, 2]}];
  Mont[Result, 1, n, NPrime, RM1, RLength]];

mem : Mont[u_, v_, n_, NPrime_, RM1_, RLength_] :=
 mem = Module[{t = u v}, BitShiftRight[(t + BitAnd[t*NPrime, RM1] n), RLength]]

Testing:

p = 2^20000 + 1;

Timing[res1 = PowerMod[2, p, p];]
(* {3.546, Null} *)

Timing[res2 = MontExp[2, p, p];]
(* {0.765, Null} *)

res1 == res2
(* True *)
$\endgroup$
4
  • 2
    $\begingroup$ You might consider using the memoization notation I proposed here. I think with exposure it will be easy to read and also more concise and less error-prone. (+1) $\endgroup$
    – Mr.Wizard
    Commented Apr 30, 2013 at 12:09
  • 1
    $\begingroup$ @Mr.Wizard, good idea, thanks. It makes for much easier reading when the pattern is long. $\endgroup$ Commented Apr 30, 2013 at 12:17
  • $\begingroup$ @SimonWoods, thanks for memoizing trick, take a look at my updates $\endgroup$ Commented Apr 30, 2013 at 15:19
  • $\begingroup$ @MohsenAfshin, I should read the question more carefully - I didn't notice that Mont version 0 was the faster one. Interesting how jagged the timings are for MontExp compared to PowerMod. $\endgroup$ Commented Apr 30, 2013 at 16:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.