2011年10月28日金曜日

素因数探し

大きい数が素数かどうか知りたいことがある. 素数と分かればそれでよし. 素数でなければ, 素因数が知りたくなるのが人情だ. 素数であるかは素数性のテストがあるので, それによればよい. 素因数探しはそれに較べ困難である.

Knuth先生のTAOCP第2巻の4.5.4項は素因数に分解する話題である. 最初のアルゴリズム4.5.4Aは2,3,5,...と順に割ってみる方法である. 素数を全部覚えているわけにもいかぬから, 2,3,5のあとは4,2,4,2,...と足して疑似素数を発生する. この辺はもう少し凝れるが, とりあえずはここまで.

その後にあるアルゴリズム4.5.4Cはちょっと面白いので, 今回はその話にしよう. TAOCPによると, この方法は1643年にFermatが使ったものらしい.

分解したい奇数nが与えられた時, このアルゴリズムは下の図の左上のように, 正の整数aとbについて, a2の正方形からb2の正方形を引いた面積をnにするのである. 灰色の部分がnになる.

そういうaとbは, 右上の図のように, aとbの差が1の時, からなず存在する. その隙間の狭い面積をnにするのである. nは奇数だから出来るわけだ. つまりa=(n+1)/2, b=(n-1)/2とすると, a2-b2=(n2+2n+1)/2-(n2-2n+1)/2=4n/4=nである.

このようなaとbが分かれば, n=a2-b2=(a+b)(a-b)だから, a+bとa-bが素因数である. 2つの素因数の差は2b.



素因数がいくつもあると, aとbの組はいくつもあったりする. n=5005なら下の2つの図のように, 712-62も732-182も5005になる.

求め方はこうだ. 最初a=floor(√n), b=0とする. つまりa2がnに等しいか, ぎりぎりに近いが少し小さ目にする. そしてr=a2-b2を計算し, r<nならaを1増やす. r>nならbを1増やす. r=nならそのaとbが求めるものだ. 素因数を求めるのに割り算をしていない.

たとえばn=9なら, a=3, b=0で決まる. 素因数は3.

n=15ならa=3, b=0, r=9から始め, r<nだからaを4にする. するとr=16になり, r>nだから今度はbを1にする. するとr=15iになり, a=4, b=1に決まる. 素因数は5と3.

n=21ならa=4, b=0, r=16から始める. そろそろ面倒になってきたから, プログラムを書こう.

(define (fermat n)
(define (try a b)
(let ((r (- (* a a) (* b b))))
(display (list a b r)) (newline) ;a,b,rを出力
(cond ((< r n) (try (+ a 1) b))
((= r n) (cons (+ a b) (- a b))) ;素因数が決まる
((> r n) (try a (+ b 1))))))
(try (inexact->exact (floor (sqrt n))) 0))

実行してみると

(fermat 21)
(4 0 16)
(5 0 25)
(5 1 24)
(5 2 21)
=> (7 . 3)

上の図の下の左は

(fermat 5005)
(70 0 4900)
(71 0 5041)
(71 1 5040)
(71 2 5037)
(71 3 5032)
(71 4 5025)
(71 5 5016)
(71 6 5005)
=> (77 . 65)

これで分かるように, このアルゴリズムはaとbの差の大きい方の解を得る.

TAOCPのアルゴリズム4.5.4Cでは, rの計算を加減算だけで出来るように, うえのaとbの代りにx=2a+1, y=2b+1を使い, rもnと比較するのでなく, r-nにして0と比較する.

こういうアルゴリズムだ.

C1 x←2(floor(√n)), y←1, r←floor(√n)2-n.
C2 if r=0,終了 n=((x-y)/2)((x+y-2)/2).
C3 r←r+x, x←x+2.
C4 r←r-y, y←y+2.
C5 if r>0 →C4, else →C2.

(TAOCP風の記述では, C2, C3のような各ステップの終わりに, 行き先(→)の指定がなければ, 次のステップへ進むことが了解されている.)

このプログラムの意外なのは, C2でr=0でなければxとyを増やしてしまう点だ. xを増やした結果はr=0にならなず, r>0になるから, yも同時に増やしている. n=19(素数)でトレースしてみる. 赤字のrはnより小さく, aを増やす時を示す.

(fermat 19)
(4 0 16)
(5 0 25)
(5 1 24)
(5 2 21)
(5 3 16)
(6 3 27)
(6 4 20)
(6 5 11)
(7 5 24)
(7 6 13)
(8 6 28)
(8 7 15)
(9 7 32)
(9 8 17)
(10 8 36)
(10 9 19)
=> (19 . 1)

nが平方数なら, n=9の例のように一発で決る. そうでないなら, aは√nより小さいからr<nになり, aを増やす. その後bを増やすとrは減って, nに等しくなるか(つまり停止するか), r<nになりaを増やす. 先ほどはr>nだったrからb2を引いてr<nになったところへ, bより大きいa2を足すのだから, b2を引くまえのrより大きくなり, r=nとなるはずはないのである. (上の結果の赤字の上下のrの値を見較べると, 下の方が大きいのが分かる.)

前のプログラムも

(define (fermat n)
(define (try a b)
(let ((r (- (* a a) (* b b))))
(display (list a b r)) (newline) ;a,b,rを出力
(cond ((< r n) (try (+ a 1) (+ b 1)))
((= r n) (cons (+ a b) (- a b))) ;素因数が決まる
((> r n) (try a (+ b 1))))))
(try (inexact->exact (floor (sqrt n))) 0))

と改良できて,

(fermat 19)
(4 0 16)
(5 1 24)
(5 2 21)
(5 3 16)
(6 4 20)
(6 5 11)
(7 6 13)
(8 7 15)
(9 8 17)
(10 9 19)
=> (19 . 1)

たしかにこの方がスマートだ.

0 件のコメント: