2011年11月6日日曜日

素因数探し

前回のブログはTAOCPのAlgorithm4.5.4Cが話題であった.

TAOCPのその次はAlgorithm4.5.4Dで, 篩を使うものである. 今回はその説明をしたい.

以下の例で, 素因数を探す数Nは, 23番のMersenne数M23=223-1
=8388607=178481×47である. また, この方法では複数の素数を利用する. それに3,5,7,11を使おう.

TAOCPの説明の通りに進めると, まず下のような表を作る. 一番左が法にする素数mである. 次にその法に現れる数x, つまり0からm-1を書く. 更にその右は, xの自乗のmod mである. つまりmを法として, 自乗の数にはこれしか現れないことを確認する.



最後は0〜m-1のxについて, (x2-N)mod mを書く. この中で, 隣りの自乗の表にあるものだけが, 考慮に値するので, それを赤で示す. 例えばm=3だと, 2,0,0のうち, 自乗の表には0と1しかないので, 2は黒のまま, 0は赤にする. するとその2つの赤なので, xの1と2の場合に対応する.

m=5だと, 赤の位置, x=1か4ならよい.

上の表の最後の数列は,

(define n (- (expt 2 23) 1))
(map (lambda (m)
(let* ((x (a2b 0 m))
(x2 (map (lambda (x) (modulo (* x x) m)) x))
(s (map (lambda (x) (modulo (- (* x x) n) m)) x)))
(display (list x x2 s)))) '(3 5 7 11))

のように計算した.

このようにして, あるxについて, xがすべての素数の法で赤の位置に対応したとき, つまり, この表で篩われた時に, x2-Nがあるyの自乗かどうかを調べるのである.

TAOCPのアルゴリズムはごちゃごちゃしているが, わたし流にSchemeで書き直すとこうなる.

(define (algorithm454d n)
(define (makesieve m n)
(let* ((x (a2b 0 m))
(x2 (map (lambda (a) (modulo (* a a) m)) x))
(s (map (lambda (b)
(if (member (modulo (- (* b b) n) m) x2) 1 0)) x)))
s))
(define ms '(3 5 7 11 13 17 19))
(define m1 (map (lambda (x) 1) ms)) ;ms length 1's
(define ss (map (lambda (m) (makesieve m n)) ms))
(display ss)
(define x 0) (define ks '())

(define (d1)
(set! x (inexact->exact (ceiling (sqrt n)))) (display x)
(set! ks (map (lambda (m) (modulo x m)) ms))
(display ks) (d2))

(define (d2)
(if (equal? (map (lambda (s k) (list-ref s k)) ss ks) m1)
(d4) (d3)))

(define (d3)
(set! x (+ x 1))
(set! ks (map (lambda (k m) (modulo (+ k 1) m)) ks ms))
(d2))

(define (d4)
(let* ((d (- (* x x) n)))
(if (integer? (sqrt d)) (let ((y (sqrt d)))
(list (+ x y) (- x y)))
(d3))))
(d1))

もう少しScheme風にも改良できそうだが, 今はこの辺で. 途中でss, xとksの初期値を出力している. それらは次の通りである.

ss=((0 1 1) (0 1 0 0 1) (1 0 1 0 0 1 0)
(1 1 0 0 1 0 0 1 0 0 1) (0 0 0 1 1 0 1 1 0 1 1 0 0)
(1 0 1 1 1 1 0 0 0 0 0 0 1 1 1 1 0)
(1 0 1 1 1 0 1 0 0 0 0 0 0 1 0 1 1 1 0))
x=2897
kx=(2 2 6 4 11 7 9)

TAOCPのアルゴリズムでは, kを(-x)mod mと計算しているが, 上の表から分かるように, 篩の値は, 0を除いて対称的なので, マイナスにする理由がなく, 私のアルゴリズムでは, x mod mで計算する.

素数をたくさん使うと, d4に来る回数が少なくなるのは当然である. ms=(3 5 7 11 13 17 19)だと523回だが, (3 5 7 11 13 17 19 23 29 31 37 41 43 47)では11回であった.

ところで, 1926年に, Henry Lehmerが自転車のチェーンをつかって篩を作った話は有名である. Mountain ViewのComputer History Museumにはその複製が展示されている.



小さい素数はチェーンが短いので, 何回か繰り返て実装されている. 篩われる数がチェーンの輪の上端に来ると, 設定してあるピンで回路が切れ, 回転が止る. つまり上のアルゴリズムのstep D4へ来る. 結果を調べ, 必要なら再起動するようになっていた.

上の例で, ピンがどうなっているかを示したのが次の図である.



篩のピンが対称的なのがよく分かるではないか.

0 件のコメント: