Function nextpp finds the next prime power after its argument.
nextpp[k_] := Module[{x},
x = k+1;
While[Length[FactorInteger[x]] > 1,
x++];
Return[x]
]
And here is the Pollard function. Instead of building a list of prime powers, which is somewhat slow, it loops and generates the next prime power using nextpp. The limit 26000 is an arbitrary one - you have to stop somewhere.
pollard[n_] :=
Module[{ b, q, g},
b=3;
q=2;
g = 1;
For[i = 2, (i < 26000) && (g == 1), i = nextpp[i],
b = PowerMod[b, First[First[FactorInteger[i]]], n];
q = Mod[q(b-1),n];
g = GCD[q, n];
Print[i, " ", g]
];
Return[g]
]
OK, here's a sample output (we could suppress the prints
to make it go faster):
In[82]:=
pollard[2^160+1]
2 1
3 1
4 1
5 1
7 1
8 1
9 1
11 1
13 1
16 1
17 1
19 1
23 1
25 1
27 1
29 1
31 1
32 1
37 1
41 1
43 1
47 1
49 1
53 1
59 1
61 1
64 1
67 1
71 1
73 1
79 1
81 1
83 1
89 1
97 1
101 1
103 1
107 1
109 1
113 1
121 1
125 1
127 1
128 641
Out[82]:=
641