2
$\begingroup$

I have been trying to write a function that duplicates PowerMod[a, b, n], computing a^b mod n. I am currently testing using 3^x mod 353 and varying x. I have found that my results match that of the built-in function until I hit x = 10^(308). At this point, I get a recursion limit error. Is there a reason this is occurring, considering Mathematica's built in function still works at these values? 

pmod[a_, b_, mod_] := 
  Module[{l, z, binarylist = IntegerDigits[b, 2], val = 1},
    l = Length[binarylist];
    Clear[z];
    z[1] = a;
    z[j_] := z[j] = Mod[z[j - 1]^2, mod];
    z[l];
    Do[
      If[binarylist[[j]] == 1,
        val * = z[l - j + 1]; val = Mod[val, 
        mod]],
      {j, 1, l}]; 
    val]

I use l - j + 1 because I want when j = 1, if binarylist[[j]] = 1; val *= z[j], when j = 2; want val*=z[j - 1], ..., j = l, val *=z[1]. This is a consequence of Mathematica, lists starting at 1, not 0.

pmod[3, 10^305, 353]

140

PowerMod[3, 10^305, 353]

140

pmod[3, 10^308, 353]

\$RecursionLimit::reclim2: Recursion depth of 1024 exceeded during evaluation of Mod[z$50178[4-1]^2,353].

185

PowerMod[3, 10^308, 353]

58

Edit

I thought this may be due to 10^308 exceeding 2^1024, but my math shows that happens at 10^309. If I should be using 2^1023 ( I don't see why i would be, but I may just be overthinking this), then that explains the error. ( It isn't GREATER than 2^1024, but including 0 in the array of the digits gives it 1024 elements. See my answer to my question below.)

Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ There are better ways to implement the Russian peasant algorithm, but as a starting point: z = Reverse[NestList[Mod[#^2, mod] &, a, l - 1]]; and then do val *= z[[j]] within your loop. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Mar 8, 2019 at 4:44
  • $\begingroup$ I know this isn't the best way, I was just curious why it was breaking down. I added some bits at the end explaining what my first thought on the issue may have been. I will implement your suggestions, and see if it stills breaks. Thanks :) $\endgroup$
    – Shinaolord
    Commented Mar 8, 2019 at 4:46
  • 1
    $\begingroup$ The warning message is pretty informative; your implementation has a recursive computation of z that is hitting the current limit. If you want to do an experiment, try Block[{$RecursionLimit = 2048}, pmod[3, 10^308, 353]]. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Mar 8, 2019 at 4:48
  • $\begingroup$ That is where My original thought at the error came from, considering 10^308 may have been (i believe, with the inclusion of 2^0=1 in binary, it is the cause of the error), and hence the list binarylist was greater than or equal to 1024, exceeding the recursion limit. I'll try modifying that limit now, as you suggested. $\endgroup$
    – Shinaolord
    Commented Mar 8, 2019 at 4:50
  • $\begingroup$ @J.M.iscomputer-less That was the problem. I can either post an answer to my own question, or allow you to point out via your own answer that I was just exceeding $RecursionLimit by having a number that when expressed in binary, had >= 1024 digits. Either way, thanks, I probably could have figured this out eventually, though who knows how long. (if you are computer less, do you just remember all this stuff? I know, bad joke.) $\endgroup$
    – Shinaolord
    Commented Mar 8, 2019 at 4:52

1 Answer 1

Reset to default
3
$\begingroup$

The issue lies in that Mathematica has a built in limit to the depth of it's recursion, which is 1024. So we cannot, without changing that limit calculate the values of z[1024],z[1025], etc. This can be changed by inserting at the beginning of the code within Module, a recursion limit greater than 1024 (I chose $RecursionLimit=2*1024). The reason this occurs at 10^308 is because the binary expansion of 10^308 has 1024 digits, the limit for recursion within Mathematica. A working code(up to 2048 binary digits, of course), is provided below:

pmod[a_, b_, mod_] := Module[{l, z, binarylist = IntegerDigits[b, 2], val = 1},
  $RecursionLimit=1024*2;
  l = Length[binarylist];
  Clear[z];
  z[1] = a;
  z[j_]:=z[j]= Mod[z[j - 1]^2, mod];
  z[l];
  Do[
    If[
        binarylist[[j]] == 1,
           val *= z[l-j+1]; val = Mod[val, mod]],
  {j, 1, l}];  
    val]

There are, however, more efficient way to implement this function, which I suspect will show it self as the number of digits increases to very large amounts, as those are likely to be implemented by the inherent PowerMod function built-in to Mathematica.

Improve this answer
$\endgroup$
7
  • 2
    $\begingroup$ It's often not a good idea to set $RecursionLimit too high globally; if you must do so, localize it in a Block[], so you have something like Block[{$RecursionLimit = 1024*2}, Module[(* stuff *)]]. Even then, all you have done is to delay the problem to b == 10^616 or so. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Mar 8, 2019 at 5:06
  • 2
    $\begingroup$ Just for your reference, here's is one way to (compactly) do the Russian peasant method: pmod[a_, b_, mod_] := Fold[Mod[If[#2 == 1, a #, #] #, mod] &, 1, IntegerDigits[b, 2]]. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Mar 8, 2019 at 5:08
  • $\begingroup$ I was indeed slightly worried about globally setting $RecursionLimit too high. I forgot that module does not automatically localize all the settings changes you make (I've only had to deal with changing $Option values probably half a dozen times). I realize all this does it delay the problem. I'll have to find a way to parse how Mathematica itself has this defined, as if I recall correctly, almost all numbers(maybe even all?) begin to exhibit a pattern that can be utilized to shorten this process. . I would expect that to greatly increase the utility of such a function. Maybe something $\endgroup$
    – Shinaolord
    Commented Mar 8, 2019 at 5:11
  • $\begingroup$ along the lines of if a pattern is noticed ( such as z[k]=z[h], z[k+1]=z[h+1], etc), to just repeatedly add that to the end of the list until the length is reached. Interesting things to think about, nonetheless. I will try to parse your compact, code, I generally have a hard time parsing such things quickly as I don't use Mathematica too much, nor the compact syntax used within it. $\endgroup$
    – Shinaolord
    Commented Mar 8, 2019 at 5:11
  • 4
    $\begingroup$ In general, in a situation like this, instead of a bare z[l] use Do[z[k], {k,l}] which does only one recursion step at a time. I do this a lot in my own code. $\endgroup$
    – Somos
    Commented Mar 8, 2019 at 5:31

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.