２．ρ method

素因数分解したい数 n について、次のような数列を考える。

x₀ = 2
x_i+1 ≡ x_i²+1 (mod n)   (i≥0)

ある数を n で割った余りは、

0, 1, 2, ...... , n-1

の n 通りしかあり得ない。
よって、上の計算をずっと続けていくと、いつかは、

x_i≡x_j (mod n)

となるようなｉ、ｊが見つかるはずである。
この時、

g = gcd(n, x_j - x_i)

を計算してやると、1 < g < n となり、n の素因数が見つかる可能性がある。
この方法を Pollard の ρ method と呼ぶ。
（この名前の由来は、上の j 以降は、j-i の周期で必ず同じ値をとるので、
　そのループをギリシャ文字の ρ に見立てて、ρ method と呼ばれている。）

上の条件を満たす時、

x_2m≡x_m (mod n)

を満たすようなｍが存在する。
よって、

x_2m - x_m (m = 1, 2, ...)

を順次計算していき、n との最大公約数をとればよい。
プログラムは以下のとおり。

10   ' rho method
20   input "n=";N
30   P=prmdiv(N):if P>1 then print P;"*";:N=N\P:goto 30
40   A=2:B=A
50   '
60   A=(A^2+1)@N
70   B=((B^2@N+1)^2+1)@N
80   G=gcd(B-A,N):if G>1 then print:print G:N=N\G:goto 80
90   goto 60

実行して、10³⁷ - 1 を入力すると、ちゃんと 2028119 という素因数をはじき出した。

上のプログラムで、A、B の添え字について着目すると、

A : ０   →   １   →   ２   →   ３   →   ４   ……
B : ０→(１)→２→(３)→４→(５)→６→(７)→８→ ……

B の（ ）に囲まれた部分は、B-A の計算に寄与していない。
というよりも、計算に必要な値を得るために、２つに１つは捨てている。
これは、もったいない。

そこで、Brent は以下のように改良型アルゴリズムを考え出した。
x_j - x_i として、以下の順に計算していく。添え字のみ抜き出すと以下のとおり。

2-1

3-2
4-2

5-4
6-4
7-4
8-4

9-8
10-8
...
16-8

17-16
...

これなら、毎回計算する x_i+1 を全て利用することができて、無駄がない。
かつ、A の値の計算回数も抑えられている。
プログラムは以下のとおり。

 10   ' rho method : Brent version
 20   input "n=";N
 30   P=prmdiv(N):if P>1 then print P;"*";:N=N\P:goto 30
 40   A=2:B=(A^2+1)@N
 50   '
 60   G=gcd(B-A,N):if G>1 then print:print G:N=N\G:goto 60
 70   I=1
 80   J=2^I:A=B
 90   for K=1 to J
100     B=B^2@N+1
110     G=gcd(B-A,N):if G>1 then print:print G:N=N\G:goto 110
120   next K
130   inc I:goto 80

素因数はそう簡単には見つからないので、何も毎回 gcd をとる必要はない。
x_j - x_iを 100個ぐらいためておいてから、gcd をとればよい。
このように改良したプログラムは以下のとおり。

 10   ' rho method : Brent version
 20   input "n=";N
 30   P=prmdiv(N):if P>1 then print P;"*";:N=N\P:goto 30
 40   A=2:B=(A^2+1)@N
 50   '
 60   G=gcd(B-A,N):if G>1 then print:print G:N=N\G:goto 60
 70   I=1:C=0:W=1
 80   J=2^I:A=B
 90   for K=1 to J
100     B=B^2@N+1:inc C
110     W=W*(B-A)@N:if C@100=0 then gosub 140:C=0
120   next K
130   inc I:goto 80
140   G=gcd(W,N):if G>1 then print:print G:N=N\G:goto 140
150   return

さらに云うなら、初期値は２である必要はない。
また、x_i+1 を求める計算式も、x²+1 である必要はない。 例えば、x³ - 1 であってもよい。
これらの修正は容易なので、プログラム例は示さない。