2011年11月30日水曜日

平方根の連分数展開

TAOCPの第4.5.3項の演習問題12は, quadratic irrationalityという話題である. よく分からぬが, どうやら平方数でない整数の平方根の連分数展開をしているらしい.

そのアルゴリズムを, さっそくSchemeに書き直した.

(define (cf d n)
(let* ((x (sqrt d)) (v 1)
(a (inexact->exact (floor x))) (u a))
(define (loop)
(display (list 'v v 'a a 'u u)) (newline)
(set! v (/ (- d (* u u)) v))
(set! a (inexact->exact (floor (/ (+ x u) v))))
(set! u (- (* a v) u))
(set! n (- n 1))
(if (> n 0) (loop) 'ok))
(loop)))

関数cfは, dの平方根の連分数の部分商をn回計算するものである. 31の平方根は

(cf 31 13)
(v 1 a 5 u 5)
(v 6 a 1 u 1)
(v 5 a 1 u 4)
(v 3 a 3 u 5)
(v 2 a 5 u 5)
(v 3 a 3 u 4)
(v 5 a 1 u 1)
(v 6 a 1 u 5)
(v 1 a 10 u 5)
(v 6 a 1 u 1)
(v 5 a 1 u 4)
(v 3 a 3 u 5)
(v 2 a 5 u 5)

となる. このaの値の5,1,1,3,5,3,1,1,10,1,1,3,5が, √31の連分数の部分商である. √31=5+1/(1+1/(1+1/(3+1/(5+....)))).

アルゴリズムに戻って, 変数はv,a,uの3個(nは繰返し数). let*で設定するvの初期値は1, aは√dの整数部分, つまり最初の部分商, uもaであって,
vn+1=(d-un2)/vn,
an+1=floor((√d+un)/vn+1),
un+1=an+1vn+1-un
を繰り返す. aに√dがあるが, 整数部分をとるので下の方はいらない. vもuも整数のままで, 誤差が入らないから, いつまでも計算が続く.

この演習問題の次の設問は, ある回数から先は, 0<un<√d, 0<vn<2√dを証明せよというのだ. 私は証明しようとも思わないが, そういうことはこの計算が周期的になるわけで, 平方根は無理数なのに, 連分数にすると繰り返すというのは面白い.

さて, このアルゴリズムは何をやっていたか調べてみる. √dをx0, 最初の部分商をa0とすると, (1)のように書ける. (TEXのプログラムでは, D,A,U,Vは大文字になっている.) ただし, u-1, v0は(0)のようだ. uの添字がvのより1だけ小さいのに注意. x=(√d+u)/vを各xの標準形とする.



(1)の添字をnだけ増やすと, 一般形xnが(2)のように得られる. anを左辺に移すと, 1/xn+1が(3)である. 上下を引っくり返せば(4)のxn+1になるが, 分母に平方根があるので, 有理化し, 左上に√dが来るように, vnで割る. これを(4)の最後のような標準形だと思うと, 対応する項から(5)の漸化式を得る. unを書いたら, 添字を1増やすとun+1が出来る. vn+1も作れるが, unの式を入れると, vの漸化式になり, たしかにschemeで書いたプログラムが得られる.

というわけで, 今回もアルゴリズムが解明出来た.

2011年11月19日土曜日

素因数探し

TAOCPの素因数探しの最初のアルゴリズム4.5.4Aは, 2,3,...と素数で次々を割ってみるやつだ.

A1. t←0, k←0, n←N.
A2. if n=1, 終了.
A3. q←floor(n/dk), r←n mod dk.
A4. if r≠0 →A6.
A5. t←t+1, pt←dk, n←q, →A2.
A6. if q>dk, (k←k+1, →A3).
A7. t←t+1, pt←n, 終了.

ここでd0=2,d1=3, ....は割ってみるべき, 次々の素数で, d2=5の次は, 2,4,2,4,...と足していくという方法があるとこのアルゴリズムには書いてある. つまり, d=2,3,5,7,11,13,17,19,23,25,29,31,35,...とする. これをみると, 25,35以外は素数で, 一見うまく行きそうだが, 要するに2と3で割れない数を仮に素数としている. 多くの素数を記憶するわけにはいかぬので, 疑似素数列を生成するわけだ.

上のプログラムをSchemeでコーディングしてみる. ptは素因数の列の配列なので, Schemeではリストpにする. A5とA7でtを1増やすのは, 配列pの添字を進めるので, リストではいらない. t←t+1, pt←n(set! p (cons n p))となる. kは約数の配列の添字だ. A6でkを増やすのは次の約数にするのだから, ここは(nextd)として, 後で考える.

(define (a2) (if (= n 1) p (a3)))
(define (a3) (set! q (quotient n d))
(set! r (modulo n d)) (a4))
(define (a4) (if (not (= r 0)) (a6) (a5)))
(define (a5) (set! p (cons d p)) (set! n q) (a2))
(define (a6) (if (> q d) (begin (nextd) (a3)) (a7)))
(define (a7) (cons n p))

a4の(not (= r 0))(> r 0)でよいが, その辺にはこだわらぬ. dは初期値を2とし, 1,2,2,4と順に足せば3,5,7,11が得られる. そこでこのリスト(1 2 2 4)の最後を(...2 ^ 2...)の間の ^ のところに繋げれば無限リスト(1 2 2 4 2 4 2 4 ...)が作れて, その後も2,4,2,4と足せるはずだ. このリストをddといおう. 次のようにして作る.

(define dd '(1 2 2 4))
(set-cdr! (list-tail dd 3) (cddr dd))

そうすれば, (nextd)

(define (nextd)
(set! d (+ d (car dd))) (set! dd (cdr dd)))
でよい.
その上と下に

(define (algorithm454a n)
(let ((p '()) (d 2) (q 0) (r 0))

(a2)))

をつければ完成で, (algorithm454a 3628800)でよべば, (7 5 5 3 3 3 3 2 2 2 2 2 2 2 2)が得られる.

このプログラムでは, 1をfactorizeしようとすると, a2でいきなり止り, 空リストが返る. TAOCPを読み直すと, 「every positive integer n」は pkを素数, p1 ≤ p2 ≤...≤ ptとして,

n=p1p2...pt

のように一意に表わせる. n=1の時はt=0で成り立つ とあるので, 空リストは当然だ. n=0ならどうか. 0をpositive integerとするかどうかだが, 上のプログラムは止らない.

ところで, 50年も前, パラメトロン計算機でこういう計算をしたころ, 疑似素数列に現れる非素数を, 我々は臨時素数と呼んでいた.

2,4,2,4と足すということは, 6の幅の間に2つの数をテストするから, 疑似素数の全整数に対する頻度は1/3=0.333...である.

2と3と5で割れない疑似素数列なら, その頻度はどうなるか. TAOCPは, 計算時間は20%節約になると説明する.

0から2,3,5のLCM, 30の前までで, 29までで, 2でも3でも5でも割れない数は, 1,7,11,13,17,19,23,29の8個なので, 8/30=0.266... . これが1/3の何%かを計算すると, (8/30)/(1/3)=8/10だから, 計算時間は80%になり, 20%の節約になる.

ついでに, TAOCPに7までにすると, さらに14%節約とあるのを確認しよう.

2*3*5*7=210個のうち, 2でも3でも5でも7でも割れないものは,

(filter (lambda (n) (> (* (modulo n 2) (modulo n 3)
(modulo n 5) (modulo n 7)) 0)) (a2b 0 (* 2 3 5 7))) =>
(1 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79
83 89 97 101 103 107 109 113 121 127 131 137 139 143 149
151 157 163 167 169 173 179 181 187 191 193 197 199 209)

このリストのlengthは48なので, 頻度は48/210. (48/210)/(8/30)=6/7=0.857... これを0.86と思えばTAOCPのいうように節約は14%である.

ところで, 2でも3でも5でも割れない数は, 0から29までに8個あると上に述べた. 計算方法はこうだ. 0から29までの数を書き, 2で割れるもの, 3で割れるもの, 5でわれるものに印をつけると図のようになる.



2で割れる15個の数には赤線を引いた. 3の10個は緑で, 5の6個は青である. これらを30から引くと, 30 - 15 - 10 - 6=-1. 6や10のように2重線のものは, 2回引いたから, 実は引きすぎである. これは1回ずつ戻さなければならない. それは6で割れる0,6,12,18,24の(30/6=)5個, 10で割れる0,10,20の(30/10=)3個, 15で割れる0と15の(30/15=)2個で, -1+5+3+2=9になる. しかし, 0は今度は3回戻されてしまった. ここは本来は1回引く場所であるが, -3+3をやってしまった. そこで(30/30=)1回を引いて, 結局四角で囲んだ8個が残る.

こういう計算法をprinciple of inclusion and exclusionという. 日本語では包除原理とか和積の原理とかいうらしい.

この原理により, 2,3,5,7,11までの頻度を計算してみよう.

(define t (* 2 3 5 7 11))
t=>2310

(- t
(+ (/ t 2) (/ t 3) (/ t 5) (/ t 7) (/ t 11))
(- (+ (/ t 2 3) (/ t 2 5) (/ t 2 7) (/ t 2 11) (/ t 3 5)
(/ t 3 7) (/ t 3 11) (/ t 5 7) (/ t 5 11) (/ t 7 11)))
(+ (/ t 2 3 5) (/ t 2 3 7) (/ t 2 3 11) (/ t 2 5 7)
(/ t 2 5 11) (/ t 2 7 11) (/ t 3 5 7) (/ t 3 5 11)
(/ t 3 7 11) (/ t 5 7 11))
(- (+ (/ t 2 3 5 7) (/ t 2 3 5 11) (/ t 2 3 7 11)
(/ t 2 5 7 11) (/ t 3 5 7 11)))
(+ (/ t 2 3 5 7 11)))
=> 480

(/ 480 2310)=>16/17=.20779...

大体1/5である.

上に書いたSchemeのプログラムは, dやpの列をリストにした以外はTAOCPのプログラムの焼き直しである. もう少しSchemeらしくしたのが, 次のプログラムだ.

(define (factor n)
(define dd '(1 2 2 4))
(set-cdr! (list-tail dd 3) (cddr dd))
(define (loop n d p)
(define (nextd)
(set! d (+ d (car dd))) (set! dd (cdr dd)))
(define (next)
(let ((q (quotient n d)))
(cond ((= (modulo n d) 0) (loop q d (cons d p)))
((> q d) (nextd) (next))
(else (cons n p)))))
(if (= n 1) p (next)))
(loop n 2 '()))

(map factor (a2b 1 30))
=> (() (2) (3) (2 2) (5) (3 2) (7) (2 2 2) (3 3) (5 2) (11)
(3 2 2) (13) (7 2) (5 3) (2 2 2 2) (17) (3 3 2) (19) (5 2 2)
(7 3) (11 2) (23) (3 2 2 2) (5 5) (13 2) (3 3 3) (7 2 2) (29))

プログラムは本質的にループだから, loopというプログラムを末尾再帰で呼ぶことにする. ある約数dで割れなければ, 次の約数を(nextd)で準備し, nextを呼ぶ. 素因数分解のたびにddを作り直すのも癪だが, そう長くはないから我慢する. さらに多くの素数の倍数をスキップするddの作り方については, またいつか書こう.

2011年11月7日月曜日

素因数探し

昨日のブログ(篩を使った素因数探し)は, TAOCPのアルゴリズムを理解するのが目的であったので, 同書のアルゴリズムの特徴であるgoto文をそのまま反映していた.

しかし, もっとSchemeらしくするにはどうするか.

前のアルゴリズムでは, 絶えずmoduloを取っているのが問題であった. あのアルゴリズムでは, 篩が2重の配列になっているので, 添字を配列の範囲に収めるためにmoduloを取るのである.

しかし, Scheme風にすると, 配列はリストになり, リストとなれば, 自転車のチェーンのように無限リストが作れる.

すなわち, (define foo '(0 1 2 3)) と設定し(図の上), (set-cdr! (list-tail foo 3) foo)
とすると, (list-tail foo 3)でfooのcdrを3回とり, 3のセルに達する. そのcdrのnilをfooに書き換えると(図の下), 3の次が0になり, 無限ループが出来る.



第2の改良点は, 篩の要素を0と1にせず, #fと#tにする. そうするとandが一発でとれる. TAOCPでは[述語]というIverson blacketを多用していて, これは述語が真のとき1, 偽のとき0になるものなので, 前回のプログラムもそうなっていた. これを(1 1 1... 1)とequal?で真偽を判定した.



lispのandは先頭からみて, 偽をみつけると途端に終了するから, この方が早いのである.

第3は, 無限リストを次々とcdr downするのに, いちいち代入するのではなく, 引数として末尾再帰で渡すことである.

このようにして書直したのが, 次のプログラムである.


(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) #t #f)) x)))
(set-cdr! (list-tail s (- (length s) 1)) s)
s))
(define s3 (makesieve 3 n))
(define s5 (makesieve 5 n))
(define s7 (makesieve 7 n))
(define s11 (makesieve 11 n))
(define s13 (makesieve 13 n))
(define s17 (makesieve 17 n))
(define s19 (makesieve 19 n))
(call-with-current-continuation
(lambda (throw)
(define (next x s3 s5 s7 s11 s13 s17 s19)
(if (and (car s3) (car s5) (car s7) (car s11)
(car s13) (car s17) (car s19))
(let ((y (sqrt (- (* x x) n))))
(if (integer? y) (throw (cons (+ x y) (- x y))))))
(next (+ x 1) (cdr s3) (cdr s5) (cdr s7) (cdr s11)
(cdr s13) (cdr s17) (cdr s19)))
(let ((x (inexact->exact (ceiling (sqrt n)))))
(next x
(list-tail s3 (modulo x 3))
(list-tail s5 (modulo x 5))
(list-tail s7 (modulo x 7))
(list-tail s11 (modulo x 11))
(list-tail s13 (modulo x 13))
(list-tail s17 (modulo x 17))
(list-tail s19 (modulo x 19)))))))


解が見つかったとき, 脱出するのにcall-with-current-continuationを使っているが, これが結局一番簡単なようである.

引数をぞろぞろ引き摺っていくのは, 素数の個数を変更するのに困るわけだが, とりあえずはこれでさくさく動く.

篩全体をリストにするプログラムも書いてはみたが, mapをとったりするので, 上のプログラムより遅かった. プログラムを動的に生成するという考えもあるが, 分かり難くもなり, 今はためらっている.

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へ来る. 結果を調べ, 必要なら再起動するようになっていた.

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



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