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

Improve this question
$\endgroup$
9
  • $\begingroup$ Look into Compile with CompilationTarget-> "C". $\endgroup$ 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$ 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$ 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

Reset to default
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 *)
Improve this answer
$\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$ 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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.