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に点を打つ.