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

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



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

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)

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

2011年10月26日水曜日

平方剰余

MathworldのQuadratic Residueを見ると, 10月1日のブログUlam Spiralのような, なんやらフラクタル風の図がある. これもやはり自分でも描くべしと, 調べてみた.

Quadratic Residue, つまり平方剰余に興味を持つ人には時々出会う.

まずaがpの平方剰余であるとは, 0<x<pのxについて, x2=a (mod p)となるxがあることだ. 従って, この範囲のについて計算すると, たとえばp=10として

(define p 10)
(map (lambda (x) (modulo (* x x) p)) (a2b 1 p))
=>(1 4 9 6 5 6 9 4 1)

従って10に対しては, 1,4,5,6,9が平方剰余であり, ここに現れない2,3,7,8が平方非剰余(Quadratic Nonresidue)である. 上の計算で見る通り, 1からp-1のxに対して, 現れる平方剰余は対称なので, 真ん中まで計算すれば良い.

従って

(define (quadratic-residue p)
(map (lambda (x) (modulo (* x x) p))
(a2b 1 (+ (quotient p 2) 1))))

(quadratic-residue 10) => (1 4 9 6 5)
(quadratic-residue 11) => (1 4 9 5 3)
(quadratic-residue 12) => (1 4 9 4 1 0)
(quadratic-residue 13) => (1 4 9 3 12 10)

0は平方剰余と言わないらしいから, 0を除き, nubを使って重複を省き, 大きさの順にするには,

(define (quadratic-residue p)
(sort (nub (filter (lambda (x) (> x 0))
(map (lambda (x) (modulo (* x x) p))
(a2b 1 (+ (quotient p 2) 1)))
)) <))

(map (lambda (p) (quadratic-residue p)) (a2b 10 14))
=> ((1 4 5 6 9) (1 3 4 5 9) (1 4 9) (1 3 4 9 10 12))

これを見るとどれにも1,4,9があるが, xが1,2,3の時のもので当然だ.

そこで, p=2から40までの平方剰余を計算してみると

(for-each (lambda (p) (display p)
(display (quadratic-residue p)) (newline))
(a2b 2 41))
=>
2(1)
3(1)
4(1)
5(1 4)
6(1 3 4)
7(1 2 4)
8(1 4)
9(1 4 7)
10(1 4 5 6 9)
11(1 3 4 5 9)
12(1 4 9)
13(1 3 4 9 10 12)
14(1 2 4 7 8 9 11)
15(1 4 6 9 10)
16(1 4 9)
17(1 2 4 8 9 13 15 16)
18(1 4 7 9 10 13 16)
19(1 4 5 6 7 9 11 16 17)
20(1 4 5 9 16)
21(1 4 7 9 15 16 18)
22(1 3 4 5 9 11 12 14 15 16 20)
23(1 2 3 4 6 8 9 12 13 16 18)
24(1 4 9 12 16)
25(1 4 6 9 11 14 16 19 21 24)
26(1 3 4 9 10 12 13 14 16 17 22 23 25)
27(1 4 7 9 10 13 16 19 22 25)
28(1 4 8 9 16 21 25)
29(1 4 5 6 7 9 13 16 20 22 23 24 25 28)
30(1 4 6 9 10 15 16 19 21 24 25)
31(1 2 4 5 7 8 9 10 14 16 18 19 20 25 28)
32(1 4 9 16 17 25)
33(1 3 4 9 12 15 16 22 25 27 31)
34(1 2 4 8 9 13 15 16 17 18 19 21 25 26 30 32 33)
35(1 4 9 11 14 15 16 21 25 29 30)
36(1 4 9 13 16 25 28)
37(1 3 4 7 9 10 11 12 16 21 25 26 27 28 30 33 34 36)
38(1 4 5 6 7 9 11 16 17 19 20 23 24 25 26 28 30 35 36)
39(1 3 4 9 10 12 13 16 22 25 27 30 36)
40(1 4 9 16 20 24 25 36)

これを使って描いたのが次の図である.


p=200までを描くと


これは, 最初に述べた, Mathworldにあった図と同じである.

PostScriptのプログラムは

40 40 translate
/d 2 def
/dot {2 dict begin /b exch def /a exch def a d mul
b d mul moveto d 0 rlineto 0 d rlineto d neg 0 rlineto
closepath fill end} def

1 1 200{/p exch def
1 1 p 1 sub{/x exch def
/a x x mul p mod def
a 0 gt {p a dot} if} for} for

のようになっている.

高さ1,4,9のような横線の他に, 右下がりの線も目立つ. つまり座標でいうと(5,4)(6,3)(7,2)(8,1)の線; (9,7)(10,6)(11,5)(12,4)(13,3)(14,2)(15,1)の線などだ. 最初の線はxとyの和が9, 次のでは和が16なのに気づく.

最初の(5,4)は, 5を法として, 4は平方剰余であるということだが, 4に法の5を足してみると9になり, 3掛ける3を5で割った余りが4なのである. その次の(6,3)は同じ9は6を法として3であるということで, この線は9を9未満の数で割った剰余の線. 次は16を16未満の数で割った剰余の線であった.

2011年10月1日土曜日

Ulam Spiral

WikipediaにUlam Spiralという項目があった.

1から図のように螺旋状に自然数を配置し, それが素数ならその場所に点を打つというものである.





私は100万くらいまでの素数のビット表を持っているから, こういう図を再現するのは何でもない.

Wikipediaには縦横200ドットの図があるので, 同じものを書いてみた.

まったく同じ図が出来て安心する.




PostScriptのプログラムはこんな具合いだ.


/n 1 def /x 0 def /y 0 def /d 2.5 def /p 1 def
/q {n primep {x y dot}if /n n 1 add def} def
/x+ {q /x x d add def} def
/x- {q /x x d sub def} def
/y+ {q /y y d add def} def
/y- {q /y y d sub def} def
100 {
p {x+} repeat p {y+} repeat /p p 1 add def
p {x-} repeat p {y-} repeat /p p 1 add def} repeat


primepは引数nが素数なら真, 合成数なら偽を返す. dotはx, yに点を打つ.

2011年9月25日日曜日

再帰曲線

このブログにdragon曲線のことを書いたのは, 2年以上も前であった.

dragon曲線は, 紙を同じ向きに半分半分と折り, 折り目を90度になるように広げたものであった.

最近, この折る向きを毎回逆にしたらどうなるかと思った. 紙を折ってやってみるのは, やはり3,4回が限度であり, 手元の計算機の威力を借りたくなった.

下の図を見て欲しい.



左端のAは折る前である. 下の黒丸は出発点を示す. その上の方に右向きの矢印が示すように, 最初は右に折る. するとBのようになる. 紙の長さが倍にになっているが, 今は折り方が問題なので, 分かりやすく描いてある.

出発点から途中右折したので, その角には+が付けてある.

Bを今度は左へ折る. 最初の折り目の+は左の下へ移動する. Cの出発点から辿ると, 最初は左折, 次がBで出来た右折, さらに右折する.

Cを次は右へ折るとDになる. 上の方の+や-は, 今回出来たのもで, 下の方のは左がCの下にあった折り目, 右がCの上にあったものだ.

出発点からの右折左折を抜き出すと

B +
C - + +
D + - - + + + -
この伝でつづけると次は
E - + + - - - + + - + + + - - +
となる.

一段上の列の要素を間に+と-を交互に入れたものになっているが, dragon曲線の場合と同じで, このパターンに着目する.

EをX3, DをX3とすると, X3の左半分はX2の+-を反対にしたものだ. 右半分は中央だけが違って右半分を殆んど同じである.

そこで, 'で+-の反転を表わすとすると,
X3=X2'+Y2'.

結局
X0=+, Y0=-,
Xn=Xn-1'+Yn-1',
Yn=Xn-1'-Yn-1'
となることが分かる.

Schemeでプログラムしてみる. この引数のnは上の漸化式の2n-1になっている.

(define (r s) (if (eq? s '+) '- '+))

(define (x n i) (cond ((= i n) '+)
((< i n) (r (x (/ n 2) i)))
((> i n) (r (y (/ n 2) (- i n))))))
(define (y n i) (cond ((= i n) '-)
((< i n) (r (x (/ n 2) i)))
((> i n) (r (y (/ n 2) (- i n))))))

(map (lambda (i) (x 8 i)) (a2b 1 16))
=> (- + + - - - + + - + + + - - +)


これで右折左折の情報が得られたから, いよいよ交互折りの曲線を描くことする.

そして出来たのが次の図である. X1からX8までが, 色を変えながら描いてある. 左中ほど上の, 直角に下の曲がった青の線がB(X1)である.



左から右向きに出発するのが, 出発点である. その下の緑の線がC(X2)に対応する.

折角描いてみたが, 竜にも鳳凰にもならず, 単に三角形が出来ただけであった. 思うにdragon曲線を考えた人も, これもやっては見たが, つまらない結果だったので, このことは書いて置かなかったのかもしれない.

つまらない絵しか描けないことが分かっただけでも, 一応の知見が得られたというべきか.

2011年9月2日金曜日

整数立方根

平方根に較べ, 立方根はおおいに疎外されている感があるが, 立方根を計算したいこともある.

その場合, 実用的には
(expt 9270 (/ 1 3)) => 21.00680051861007
なことで済ませてしまう.

さて, Warrenの「ハッカーの楽しみ」を眺めていたら, 整数立方根という話題があった(訳書の226ページ). こういうプログラムが書いてある(32ビット用だ).



同書の巻頭の「推薦の辞」に私が書いたように, 本書にはそのアルゴリズムでよい理由が殆んど書いてない.

このアルゴリズムは下の説明のように出来ているのだ.

例によってSchemeに書き直すと




最初のsの設定が30でなく18なのは, そんなに大きい数でテストすることもないからだ. また最後にyだけでなく, xも返すのは, 余りも欲しいからである.

途中の経過を見るdisplayが3行ある. 9270の立方根を計算してみる.



すなわち, 立方根は21で, 剰余が9. 整数の範囲で開立をやめ, 剰余が得られるから, 整数立方根というわけである. 同じものを十進法で書くと,



つまりbの始めの方は, xより小さくなるまで218, 215, 212,... と減っていき, 4096がxより初めて小さくなった時点でyを1とする. yは最後は二進で10101だが, その最初の1(=24)が決ったのである.

xの方はそれを引いて, 9270-4096=5174になった.

次はyのその下のビットを1にするか, 0にするかを決めることである. y=16の時, その下のビットは8だから(16+8)3=163+3*162*8+3*16*82+83=13824

しかし, 最初の4096はすでに引いてあるので, 比較するbは, 上の4段目の13824-4096=9728である. こうするよりは, y=2, (y+1)3-y3=3y(y+1)+1を作り, 左シフトする方が楽というのが, プログラムの趣旨であろう.

8のビットは立たなかった. 次は4のビットで, (16+4)3-163=3904である. これはxより小さいからyは1増えて101になった.

このようにしてyのビットを上から次々を立てていき, 1ビットごとにきめていく. 最後に残ったxが開立の剰余である.

このプログラムの場合は, 方針は最初から見え見えであったが.