/* from https://math.stackexchange.com/questions/630134/calculating-an-pmod-m-in-the-general-case */ /* calculates (b^^i) mod m */ debugprint = 0; reducetower(b,m,i) = { local(n,rz,log2p); /* these are the recursion exit conditions */ if ((b%m)==0, return(0)); /* "0" if b divisible by m, or m=1 */ if ((i==2), return(power2n(b,m,b))); /* if i==2, calculate (b^b) mod m */ /* boundary cases for m=3n*(3^^3), just detect m>3^^3, similary for 2^^4 */ /* and similarily for m=2n*(2^^3), just detect m>2^^3 */ if ((b==2) && (i==3) && (m>16), return (16)); if ((b==2) && (i==4) && (m>65536), return (65536)); if ((b==3) && (i==3) && (m>7625597484987), return (7625597484987)); /* theoretically, there is also (b==4) && (i==3) && (m>4^^3) */ log2p=0; n=eulerphi(m); /* n replaces m in recursion, reducetower(b,n,i-1); */ if (n+1<m, /* if m is not a prime */ log2p=floor(log(m)/log(2)); /* repeating pattern length n guarenteed */ ); rz = reducetower(b,n,i-1); /* recursion with i-1, and n replacing m */ /* we know there are repeats beginning at or before log2p, so that */ /* b^log2p = b^(log2p+n) */ /* b^b^z can be replaced with b^(b^z mod n) as long as (b^z mod n)>log2p */ /* if b^z<log2p, keep adding (b^z)+n until >= log2p */ while ((log2p>0) && (rz<log2p), rz=rz+n); rz = power2n(b,m,rz); /* (b^rz) mod m */ if (debugprint, print("("b"^^"i")%"m"=("b"^(("b"^^"i-1")%"n"))%"m"="rz); ); return(rz); } power2n(z,m,i) = { local(t,y); if (i==0,return(1)); y=z; t=1; /* express i in base 2, multiply the subparts of i as they are calculated */ while ((y<>1) && (i>0), if ((i%2)==1, t=(t*y)%m;i=i-1); y=(y*y)%m; i=i/2; ); return(t); }