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.