<< NumberTheory`NumberTheoryFunctions`

\(F[{y_, r_, s_}, p_] := 
      Module[{}, 
If[
          0 < y ≤ 
            p\/3, 
				Return[{	Mod[g*y, p], 	
              Mod[r + 1, 
                p - 1], 	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ s\ }]]; 
If[
          p\/3 < y <= \(2  p\)\/3, 
				Return[{	
              Mod[y*y, p], 	Mod[2*r, p - 1], 	
              Mod[2*s, p - 1]}]]; 
If[\(2  p\)\/3 < 
            y ≤ 
            p, 
				Return[{	
              Mod[h*y, p], 				\ \ \ \ r, 	
              Mod[s + 1, p - 1]}]];
]\)

PollardRho[g_, h_, p_] := 
      Module[{r0, y, k, t, done, w, ergs, a, b, d, a2, b2, 
          x}, 

r0  =  
          Random[Integer, {1, p - 1}]; 
y  =  {PowerMod[
              g, r0, p], r0, 
            0}; 

k  =  
          0; 
t  =  1; 
done  =  
          False; 
While[
          Not[done], 
w  =  y; 
While[
            Not[done] && 
              k < t, 
k++; 
y  =  
              F[y, p]; 
done  =  (w[[1]] == 
                  y[[1]]);
]; 
t  =  2*t;
]; 

ergs  =  {w[[2]], w[[3]], y[[2]], 
            y[[3]]}; 
a  =  
          Mod[ergs[[1]] - ergs[[3]], 
            p - 1]; 
b  =  
          Mod[ergs[[4]] - ergs[[2]], 
            p - 1]; 

d  =  
          GCD[b, p - 1]; 
a2  =  
          a/d; 
b2  =  
          b/d; 
p2  =  (p - 1)/
            d; 

x  =  
          Mod[a2*PowerMod[b2, -1, p2], p2]; 
k  =  
          0; 
While[
          PowerMod[g, x + k*p2, p] ≠ 
            h, k++]; 

Return[
          x + k*p2];
]

\(PollardRhoVerbose[g_, h_, p_] := 
      Module[{r0, y, k, t, done, w, ergs, a, b, d, a2, b2, 
          x}, 
r0\  = \ 
          Random[Integer, {1, p - 1}]; 
y = {PowerMod[g, 
              r0, p], r0, 0}; 
Print["Starte bei y = ", 
          y]; 
k = 0; 
t = 
          1; 
done\  = \ 
          False; 
Print["Berechne Folge:"]; 
While[
          Not[done], 
w\  = \ y; 
Print[
            k, ":  ", w]; 
While[
            Not[done] && 
              k < t, 
\(k++\); 
y = 
              F[y, p]; 
done = \((w[\([1]\)] == 
                  y[\([1]\)])\);
]; 
t\  = \ 2*t;
]; 
Print[k, ":  ", 
          y]; 
ergs = {w[\([2]\)], w[\([3]\)], y[\([2]\)], 
            y[\([3]\)]}; 
a = 
          Mod[ergs[\([1]\)] - ergs[\([3]\)], p - 1]; 
b = 
          Mod[ergs[\([4]\)] - ergs[\([2]\)], 
            p - 1]; 
Print["Also ist"]; 
Print["         g ^ ( ", 
          a, " ) ≡ h ^ ( ", 
          b, " ) ≡ g ^ ( x * ", b, " )   mod ", 
          p]; 
Print["Löse also"]; 
Print["         x * ", b, " ≡ ", 
          a, "   mod ", p - 1]; 
d\  = \ 
          GCD[b, p - 
              1]; 
Print["Der größte gemeinsame Teiler von ", b, " und ", p - 1, " ist ", 
          d]; 
Print["Reduziere die Gleichung also um den Faktor ", d]; 
a2\  = \ a/d; 
b2\  = \ 
          b/d; 
p2\  = \ \((p - 1)\)/
            d; 
Print["         x * ", 
          b2, " ≡ ", a2, "   mod ", 
          p2]; 
x\  = \ 
          Mod[a2*PowerMod[b2, \(-1\), p2], 
            p2]; 
Print["Also ist x = ", 
          x]; 
Print["Lösungen sind also ", 
          x, " + k * ", p2, \ " wobei 0 ≤ k < ", 
          d]; 
k = 0; 
While[
          PowerMod[g, x + k*p2, p] ≠ h, 
          Print["Teste:   k = ", k, "   ->   x = ", 
            x + k*p2, \ "  ->  Kein Erfolg"]; \(k++\)]; 
Print["Teste:   k = ", k, "   ->   ", 
          x + k*p2, \ "  ->  Erfolg"]; 
Return[
          x + k*p2];
]\)

p  =  2017;
g  =  5;
h  =  3;
{{g, h, p}, MultiplicativeOrder[g, p, h]}

PollardRho[g, h, p]

PollardRhoVerbose[g, h, p]

p  =  
        Random[Integer, {10^6, 10^7 - 1}];
While[! PrimeQ[p], p++];
g  =  PrimitiveRoot[p];
lsg  =  Random[Integer, {1, p - 1}];
h  =  PowerMod[g, lsg, p];
{{g, h, p}, {lsg}}

Timing[PollardRho[g, h, p]]

Timing[MultiplicativeOrder[g, p, h]]