Optimized Version of  Pollard's  r  Method

home   |   index
Rho Factorization Method:
 
The optimized version below
is roughly 3 times faster than
its simpler counterpart at left.

In our general discussion of factorization methods, we're giving the backbone of Pollard's rho method as the following terse UBASIC  implementation  (which doesn't check the primality of the factors it prints out).

input "Integer to be factored";N

U=1 : REM U is the value at index K   (K is a power of 2)
K=1
V=U : REM V is the value at index K+Q (the "running" value)
loop
   for Q=1 to K
      V=(1+V^2)@N : REM  x@y  stands for "x modulo y"
      D=gcd(V-U,N) : REM This is a time-consuming step !
      if D>1 then print D : N=N\D : if N=1 then end
   next Q
   U=V
   K=2*K
endloop

It's an overkill to compute at each iteration a  greatest common factor  (GCD)  which is almost always equal to 1...  Instead, we could "accumulate", as a product modulo N, many consecutive values  (say, 100 of them)  and compute instead the greatest common factor of  N  and  that  (which is still almost always equal to 1 when large numbers are involved).  This essentially replaces a GCD operation by a multiplication modulo N  (except once in 100 iterations).  Doing just that would  usually discover the same large factors.  However, we may even afford the luxury to design a procedure which retrieves  exactly  the same factors by retracing the latest steps of the above procedure whenever (at least) one of them is bound to be successful...  This is what's done in the program listed below.

 10   M=99
 20   dim A(M)
 30   '
 40   input "Integer to be factored";N
 60   '
 70   U=1:' Value at index K (K is a power of 2)
 80   K=1
 90   V=U:' Value at index K+Q
100   J=0:D=1
110   loop
120   Q=0:while Q<K:inc Q
130   V=(1+V*V)@N
140   A(J)=V-U
150   D=(A(J)*D)@N
160   inc J
170   if J<=M then 290
180   D=gcd(D,N)
190   if D=1 then 280
200   for J=0 to M
210   D=gcd(A(J),N)
220   if D=1 then 270
230   print D;"(";K+Q-M+J;")"
240   N=N\D
250   if N=1 then print:cancel for:goto 40
260   D=1:U=U@N
270   next J
280   J=0
290   wend
300   U=V
310   K=2*K
320   endloop

One reason the idiom of line 120 replaces the more straightforward "for Q=1 to K" (used in the didactic version) is that the latter UBASIC construct would cause an error when K reaches  4294967296,  as would happen in a couple of days with a 3 GHz microcomputer.


visits since June 3, 2006 Valid HTML