素因数探し

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

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について, a²の正方形からb²の正方形を引いた面積をnにするのである. 灰色の部分がnになる.

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

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



素因数がいくつもあると, aとbの組はいくつもあったりする. n=5005なら下の2つの図のように, 71²-6²も73²-18²も5005になる.

求め方はこうだ. 最初a=floor(√n), b=0とする. つまりa²がnに等しいか, ぎりぎりに近いが少し小さ目にする. そしてr=a²-b²を計算し, 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)²-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からb²を引いてr<nになったところへ, bより大きいa²を足すのだから, b²を引くまえの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)

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

平方剰余

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

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

まずaがpの平方剰余であるとは, 0<x<pのxについて, x²=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未満の数で割った剰余の線であった.

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

再帰曲線

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

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

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

下の図を見て欲しい.



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

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

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

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

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

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

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

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

そこで, 'で+-の反転を表わすとすると,
X₃=X₂'+Y₂'.

結局
X₀=+, Y₀=-,
X_n=X_n-1'+Y_n-1',
Y_n=X_n-1'-Y_n-1'
となることが分かる.

Schemeでプログラムしてみる. この引数のnは上の漸化式の2^n-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))
=> (- + + - - - + + - + + + - - +)

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

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



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

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

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

整数立方根

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

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

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



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

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

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




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

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



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



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

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

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

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

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

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

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

多面体描画道楽

Archimedes多面体13種のうち, rhomb(斜方とか菱形)という形容詞のつくものが2つある. 斜方立方八面体と斜方二十十二面体である. 何ゆえに斜方といわれるのか.





これらの図に示すように, この2種類の多面体には, それぞれ3種類の面がある.

斜方立方八面体では, 青は元々の立方体の面, 緑は正八面体の面である. 斜方二十十二面体では, 青は元々の正十二面体の面, 緑は正二十面体の面である.

では, それぞれの赤の面はなにか.

赤の面を延長して, 緑の面の上方での交点を決め, 赤の正方形に外接する菱形を考えることが出来る. 赤は元々そういう面の多面体の面である.

この多面体がどちらも菱形になるので, 斜方といわれる所以である. つまり, 斜方(赤)立方(青)八(緑)面体, 斜方(赤)二十(緑)十二(青)面体という命名であったわけだ.

ではその斜方多面体を描いてみよう.


上は斜方立方八面体である. ｘ, y, zの対称軸も描いてある. z軸(青)のまわりに少し回転し, y軸(緑)のまわりにも僅かに回転したようになっている.


そのz軸, y軸まわりの回転をやめ, x軸方向から眺めたのが, この上の図だ. そこで, 最初の赤い正方形を含む面を考えてみると, この図の赤線で示す菱形が得られる. その座標が分かれば, 次の図が描けるわけだ.



これは, 対角線の比が1:√2の菱形が12個で構成されているので, 斜方十二面体(rhombic dodecahedron)という.


斜方二十十二面体の方は, 面が多いせいか, 多少手ごわい.



これが元の斜方二十十二面体で, 上と同様にz軸, y軸について回転してある.



回転角を0とし, x方向から眺めると, このように見えるはずだ. 赤線はそこにかぶせた菱形である. この方は, 菱形30枚で構成され, 斜方三十面体(triacontahedron)という.



菱形の対角線の比は1:φ(黄金比)である.

同一の菱形だけで出来る多面体はこの2つだけらしい. またKeplerはこの2つの多面体のあることを知っていたといわれる.

3シリンダ機関車

ずいぶん前のこのブログに3シリンダ機関車のことを書いた. そこで引用しているURLが, 最近なくなっているのに気づいたので, もう一度3シリンダ機関車の弁装置について書くことにする.

今でも蒸気機関車は人気があり, 各地で復活されているが, その構造はあまり理解されていないに違いない.

簡単にいってしまえば, 機関車の左右の前方につけた2個のシリンダ(気筒)の中にピストンを入れ, ピストンの両側から交互に蒸気を補給してピストンを前後に動かす. ピストンの往復運動を主連棒で繋ぎ, 動輪を回転して機関車を進める. ピストンがどちらかの端に寄っていると, 蒸気の補給が出来ず, 死点といって力が出せないが, 反対側のピストンが最も力の出る中央にいるように, 左右のピストンの位相を90度づらせて配置してある. このことは, 鉄道博物館で実物の機関車を見ると確認出来る.

ピストンの前後に交互に蒸気の補給をしているのが, シリンダと並んだ蒸気室で, その中に滑り弁があり, これが蒸気室の中を前後に滑りながら, 蒸気をピストンの一方に送り込む.

Walschaertの走り装置の図はここにあるので参照されたい.

滑り弁は, 動輪に動力を伝える主連棒と90度づれたエキセン棒から, 逆転装置を経由した弁棒で押し引きされて運動する. この辺の装置は大変に込み入っているから, 説明は省略!

ところで, 表題の3シリンダ機関車は, 2つのシリンダでは力不足の時, 第3のシリンダを両シリンダの中央に搭載し, 合わせて3本の主連棒で動輪を駆動する方式である. 動輪軸の中央はクランク状になっていて, そこに中央の主連棒が接続されている.

問題は, 左右の動輪のエキセン棒のようなものが, 中央のクランクからも出ているかということだが, そういう複雑な走り装置にはなっておらず, 中央のシリンダの弁の制御は両側の弁の制御から構成していたというのが前回のブログであった.

下の図を見て欲しい. Cylinder AとCylinder Bと書いたのが, 両端のシリンダの蒸気室で, 左(機関車の後方)から弁棒が入り, 途中に滑り弁があり, 弁棒の延長が右方に突き出している. 機関車の先頭にAOD, BDCのような梃子があり, 支点AはCylinder Aの弁棒の先に, 支点BはCylinder Bの弁棒の先に固定されている. 梃子だから, 図で左右に揺れると支点の上下移動もあるはずだが, 今は各支点は左右にだけ動くとしている.

支点Oは固定点である. またOD=OA/2, つまりAが動くとDはAの半分の振幅で逆位相に動く. そのDを支点として, DB=DCの梃子があり, CがCylinder Cの弁棒になっている. そうすると, CはAとBが120度の位相差で動くと, さらに120度の位相差で動き, 3番目の滑り弁が望み通りに動くのである.

このアニメーションがここにある.

アニメの下でくるくる回る3芒星は, 120度の位相で回転することを示す.