５．連分数法（Continued Fraction Method）

ちょっと横道にそれる。
パズル等でよく使われる、12345679 という８桁の数がある（８が抜けていることに注意）。
この数は、

12345679× 9 = 111111111
12345679×18 = 222222222
12345679×27 = 333333333
12345679×36 = 444444444
12345679×45 = 555555555
12345679×54 = 666666666
12345679×63 = 777777777
12345679×72 = 888888888
12345679×81 = 999999999

という性質がある（何故こうなるかは、１行目より明らかと思う）。
この 12345679 は、素数か否か？

prmdiv()を使えば一発でわかるが、これは 37 という約数を持つ。

12345679 = 37×333667

で、この 333667 も素数であり、これで素因数分解は完了である。
37程度であれば、筆算でも見つけられる範囲である。

似たような数で、123456789 という 9桁の数はどうか。
これは、9 という自明な約数を持つ。

123456789 = 9×13717421

この 13717421 は素数か？
これも prmdiv() を使えば一発でわかるが、 3607 という約数を持つ。

13717421 = 3607×3803

であり、さすがにこの 3607 を筆算で求めるのは、かなり困難であろう。
しかし、これを筆算で見つける方法がある。

中学以来おなじみの、以下のような因数分解の公式がある。

x² - y² = (x + y) (x - y)

さて、ある数 n が上のように２つの２乗数の差で表すことができたとしよう。
このとき、上の因数分解の公式をつかって、n を x+y と x-y に素因数分解することができる。

x の初期値として、

x = int(sqrt(n)) + 1

から始める。y の初期値は、x² - n の平方根の整数部、すなわち、

y = int(sqrt(x² - n))

とする。この初期値から以下の処理を繰り返す。

1. w = x² - n - y² を求める
2. w = 0 なら、x+y, x-y が求める素因数である
3. w > 0 なら、y = y + 1 とする
4. w < 0 なら、x = x + 1 とする
5. 1. へ

プログラムは以下のとおり。

 5   ' Fermat method
10   input "n=";N
20   X=isqrt(N)+1:Y=isqrt(X^2-N):print X,Y
30   W=X^2-N-Y^2
40   if W=0 then print X+Y,X-Y,X,Y:end
50   if W>0 then inc Y:goto 30
60   if W<0 then inc X:goto 30

y に比べて x の増え方は大きいので、x が変わるごとに、y を計算し直すようにした方がよい。
このプログラムを実行して、13717421 を代入すると、x の初期値は 3704。
次の 3705 で、y=98 となり、因数 3803、3607 が見つかる。
この方法は、素因数が n の平方根に近いとき、一瞬にして見つけられる。この方法を Fermat の方法 と呼ぶ。

さて、n = x² - y² とならないまでも、

x² - y²≡0 (mod n)

であれば、x-y と n の gcd をとって、n の素因数が求められる場合がある。
このような x, y を求める方法としては、Pell方程式の解法を利用する。

Pell方程式の解を求める過程で、

p_i² - nq_i² = (-1)ⁱ * b_i

という式が得られる。
この右辺が２乗数であれば、

p_i² - (-1)ⁱ * b_i≡0 (mod n)

であるから、話は終わり。
もし２乗数でなくても、いくつかの (-1)ⁱ * b_i を求めて、
その積が２乗数になるように構成してやればよい。
この方法を 連分数法 (Continued Fraction Method) と呼ぶ（このアルゴリズムは Legendre による）。

Morrison と Brillhart は、この方法により、７番目の Fermat 数の素因数分解、

F₇ = 2²⁷ + 1 = 2¹²⁸ + 1
  = 59649589127497217×5704689200685129054721

を得た（1970）。

このアルゴリズムは途中の連分数を求めていく過程で大きな数となり処理時間がかかるため、
現在では、複数多項式二次ふるい法にとって変わられている。