# My version of PowerMod breaks down around 10^308

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

• 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. – J. M.'s ennui Mar 8 '19 at 4:44
• 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 :) – Shinaolord Mar 8 '19 at 4:46
• 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]]. – J. M.'s ennui Mar 8 '19 at 4:48 • 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. – Shinaolord Mar 8 '19 at 4:50 • @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.) – Shinaolord Mar 8 '19 at 4:52

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.

• 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. – J. M.'s ennui Mar 8 '19 at 5:06
• 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]]. – J. M.'s ennui Mar 8 '19 at 5:08
• 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 – Shinaolord Mar 8 '19 at 5:11
• 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. – Shinaolord Mar 8 '19 at 5:11
• 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. – Somos Mar 8 '19 at 5:31