Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
    
Look into Compile with CompilationTarget-> "C". –  0xFE Apr 23 '13 at 14:13
1  
You have to be careful though because Compile expects machine-sized integers in most places. You'd need to use an alternate representation. –  0xFE Apr 23 '13 at 14:42
1  
I tied it and because of the "machine-size" integers it's useless for this scenario. –  Mohsen Afshin Apr 23 '13 at 15:24
    
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. –  0xFE Apr 24 '13 at 2:33
2  
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? –  Oleksandr R. Apr 24 '13 at 10:49
show 4 more comments

1 Answer

up vote 5 down vote accepted
+50

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 *)
share|improve this answer
1  
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) –  Mr.Wizard Apr 30 '13 at 12:09
1  
@Mr.Wizard, good idea, thanks. It makes for much easier reading when the pattern is long. –  Simon Woods Apr 30 '13 at 12:17
    
@SimonWoods, thanks for memoizing trick, take a look at my updates –  Mohsen Afshin Apr 30 '13 at 15:19
    
@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. –  Simon Woods Apr 30 '13 at 16:18
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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