除算サブルーチン

パラメトロン計算機 PC-1 が完成し稼働し始めたのは, 1958年3月26日ということになっているが, この時は乗算も除算も出来なかった. その後, それらの回路を設計し, やがて乗算が付き, 夏頃には除算も出来るようになった.

イニシアルオーダーは十進二進変換もするが, その際の10倍に乗算を使わず, シフトしたり足したりしているのは, 当時, 乗算命令が未実装だったからである.

従って, 始めの頃は, 乗算や除算のサブルーチンが存在した. 私は数値計算には興味がないので, 除算は必要なく, 高橋先生が作られたはずの除算のサブルーチンがどうなっていたかは知らない. 先生のライブラリのファイルが手元に来たので, 除算サブルーチンを探したが, 見つからなかった.

ところで, PC-1がお手本にした, Cambridge大学のEDSACには, 最後まで除算はなく, サブルーチンを使っていた. そういえば, EDSACでは除算はどうやったか知りたくなり, 例のEDSACの本を見ると, 思いの他簡単であった.

a₀←被除数,
c₀←除数-1
と初期化し, 漸化式
a_k+1←-a_kc_k+a_k,
c_k+1←-c_k².
で繰り返すと, c_k=0の時, a_kが商になる.

これをそのままSchemeでプログラムにする.

(define (div n d) ;n 被除数 d 除数
 (define (dv a c)
  (display (list a c)) ;途中の様子を見る
  (if (< (abs c) 0.0001) a
   (dv (- a (* a c)) (- (* c c)))))
 (dv n (- d 1)))

PC-1もEDSACも当時の計算機は, -1≤n<1の固定小数点数しか扱わない. 利用者は, 計算がこの範囲になるよう, 絶えず位取りに注意する必要があった.

そこで, 上のサブルーチンを(div 0.4 0.8), (div 0.4 0.5)でやってみると,

> (div 0.4 0.8)
(.4 -.19999999999999996)
(.48 -.03999999999999998)
(.4992 -1.5999999999999983e-3)
(.49999871999999995 -2.5599999999999945e-6)
=> .49999871999999995

> (div 0.4 0.5)
(.4 -.5)
(.6000000000000001 -.25)
(.7500000000000001 -.0625)
(.7968750000000001 -.00390625)
(.7999877929687501 -.0000152587890625)
=> .7999877929687501

正しい商は, 0.5と0.8である. 目がちらちらするから, 途中のa_k, c_kを整数でシミュレートしてから同様な表にすると,

0.4/0.8
k a          c           delta
0 0.4        -0.2        0.1
1 0.48       -0.04       0.02
2 0.4992     -0.0016     0.0008
3 0.4999872  -0.00000256 0.0000128
4 0.49999999 -0.00000000 0.00000000
    99967232    00065536   00032768

0.4/0.5
k a          c           delta
0 0.4        -0.5        0.4
1 0.60       -0.25       0.2
2 0.7500     -0.0625     0.05
3 0.796875   -0.00390625 0.003125
4 0.79998779 -0.00001525 0.00001220
    296875      87890625   703125

左端のkはaとcの添字, 右端のdeltaは正確なn/dとa_kとの差である. 最後のk=4の行は長いので, 少数点以下8桁目で折り返してある.

こういうのの性質をしらべるには, 各回におけるinvariantを発見するのが常道である. ためつすがめつするまでもなく, この例ではcとdeltaには比例関係があることが分かる. つまり, 上の計算ではdelta/cは-1/2, 下のではdelta/cは-4/5, つまり除算の商の符号を変えたものであった.

すなわち, n/d=a_k-(n/d)*c_kということだ.

ここまで分かると, 次は帰納法で証明するだけである.

k=0の時

a₀-(n/d)*c₀=n-(n/d)*(d-1)=n/d.

k=kで

a_k-(n/d)*c_k=n/d

が成立しているとする. その時k+1の式は

a_k+1-(a/d)*c_k+1
=-a_kc_k+a_k-(n/d)*(-c_k²)
=-a_kc_k+a_k+(n/d)*c_k²
=c_k(a_k-(n/d)*c_k)+a_k
=a_k-(n/d)*c_k)=n/d.

やはり, n/dになる. 一方 c_kはどんどん小さくなっていくから, c_k=0の時, a_k→n/d.

そういう風にEDSACの除算サブルーチンは出来ていたのであった.

世界一周プログラム

1月15日は私の恩師, 高橋秀俊先生の誕生日で, OBが先生のお宅に集まった. 先生が他界されたのは, 1985年だから, それから27年が経った. 成人の日は1月の第2月曜に変ったから, 15日ではなくなったが, 成人の日には今でもOBが某所に集まっている.

今回, その会で思わぬものを渡された. 高橋先生が東大にいらした頃の, パラメトロン計算機 PC-1 のライブラリのファイルである.



先生は1975年に東大を定年退官され, 慶応義塾大学へ移られた. 東大の教授室を後片付けした筑波大学の長谷部君が, そのファイルを見つけ, 自分の部屋で保管していたが, その長谷部君も定年になり, 筑波大学の教授室を片付けた. そしてそのファイルを私に委ねてきたのであった.

PC-1のライブラリを整備していた頃, 紙の質はまだ悪く, 日焼けで茶色になり, 縁もぼろぼろになってきた. 私が東大工学部にいた時, このままでは持たないと思って全部コピーしなおし, 製本した. その折, もう1部を作って高橋先生に差し上げたら大変喜ばれ, 「これを見ると, PC-1をやっていた頃が昨日のように思い出されます」という, 葉書の礼状を下さった.

先生がPC-1のライブラリファイルを, 東大に残していかれたとは, この正月まで知らなかった.

こうして手元にめぐってきたライブラリを見ると, 私の書いたイニシアルオーダーの前に, 先生が書かれたイニシアルオーダーがあった. 確かに1958年3月26日稼働開始の直後は先生のイニシアルオーダーが使われていたが, 長いし遅いので, 連休中に私が書き直したのが, その後1964年まで, 長く使われた.

先生のファイルには, "self reproducing program"もあった. これは計算機のメモリーをテストするため, 短いプログラムを走らせると, そのプログラムが自分を, アドレス部を修正しながら, 少し先のメモリーにコピーし, コピーが終わると, 新しい部分が走るようなものである.

Illinois大学のIlliac Iで使っていたと聞いていて, 高橋先生も書いてみられたのだろう.

Illiacのは, アドレス部がそれぞれその場所になっている, 2つのプログラムがあり, その差から新しいコピーを作っていると聞いたが, PC-1のは違っていた. (Illiacではこのプログラムをleapfrog(カエル跳び)と呼んでいたらしい.)

先生のスケッチは途中らしく, 正しくなさそうなので, 実際に使われていたプログラムを紹介しよう.



やはり, PC-1の命令から説明しなければなるまい. 1語は18ビットで, そのうち左端の6ビットが命令コード, 上のプログラムで, 各命令の先頭の"o"や"p"や"j"が入る. その右の1ビットが, l/sビットで, オペランドが長語(36ビット)か, 短語(18ビット)を区別し, 長語なら, 命令に"l"をつけ, そのビットが1になる. 残りの11ビットがアドレス部だが, PC-1は512語しかなかったので, 9ビットだけが必要, 上の2ビットは無視する.



0番地は"o"命令, outputの意味で, アキュムレータの左端の6ビットをテレタイプコードとして送り出す. アドレス部は, 命令の置いてある番地とすることになっているから, "0"である.

1番地は"p 16". "p"はpositive loadで, 16番地の短語をアキュムレータに取り出す. 従ってアキュムレータは, 16番地の命令"jl 17"になる.

2番地"jl"は, jump命令で, "j"の場合, lは長短ではなく, この場合は無条件ジャンプである. 従って, 3番地へ飛ぶ.

3番地"x 8"は, アキュムレータのアドレス部だけ, 8番地のアドレス部へ入れる命令. したがって, 8番地は "p 17"になる. 8番地のアドレスにかっこがついているのは, このアドレスは変更されるという意味である.

4番地の"a"はadd, 5番地の"s"はsubtractで, アキュムレータの"p 17"に, 16番地の"jl 17"を足し, 0番地の"o 0"を引く. 命令コード部も変るが今は関係なく, 17に17を足し0を引くから, アドレス部は17増えて34になり, 6番地, 7番地の"x"命令で, 9番地と12番地のアドレス部に入れられる.

このように, 定数を足すのにも, アドレス部は絶えず変るので, 定数を持っているわけにはいかず, 差を使う必要があった.

8番地に来て, "p 17", つまり"jl 3"をアキュムレータに持ってきて, それを9番地の"t", transfer命令で34番地にいれる. 10, 11番地は4, 5番地と同様で, アキュムレータのアドレス部を17増やし, 20にしてそれを34番地のアドレス部に入れる. 最終的に34番地は"jl 20"になる.

13番地は"p 8"だから, アキュムレータは"p 17"になる. 14, 15番地は, 14を足し, 15を引くから, 1引くことになり, アドレス部は16になり, 17番地へ飛び, さらに3番地へ飛ぶと, 今度は8番地, 9番地, 12番地のアドレス部がそれぞれ1減り, 16番地の"jl 17"が"jl 34"になって33番地に入る.

これを繰り返していると, やがて1番地の"p 16"が18番地に"p 33"として入り, 次のサイクルでは, "o 0"が"o 17"になって7番地の"jl 3"に上書きされる. コピーサイクルが17番地に到達した時点からコピーされたプログラムが走り出し, 同じことを繰り返すのである.

制御が新しいコピーに移ると, そこに出力命令があるので, プリンターが1文字打出す. つまりプリンターが文字を打ち続けている限り, メモリーは正常であるわけだ. 512/17=30だから, 30回のコピーでメモリーの始めに戻る. PC-1では, これを世界とみて, このプログラムを世界一周と呼んだ. 80日世界一周という映画があったから, ○○秒世界一周ともいったが, それが何秒だったかは, 残念ながら思い出せない.

素因数探し

昔は計算機はなくても, 計算用の道具はいろいろあった. 12月3日のブログに書いたfactor stencilもそういうものの1つである. 竹内端三先生はfactor stencilを「因数型紙」と訳された.

Derrick Norman Lehmerがfactor stencilを作ったのは1925年頃のことである. それ以前, 1914年にLehmerは10006721までの素数表を完成し出版した. といってみるだけは簡単であるが, 665000個の素数があるわけで, (Legendreの素数定理で計算すると(define (pi x) (/ x (- (log x) 1.08366))) (pi 10006721) => 665556.99...) 1ページに5000個ずつ, 133ページに収められている.

その第1ページ目には48593までの素数が縦100行, 横50列に印刷してある. factor stencilはこの第1ページの素数の, それぞれのRに対応する位置に穴を開けて作られた. Lehmerの当時の寄書に, "The device for cutting the stencils is already constructed and..."とあるから, 何かの器具を工夫したのであろうが, 詳しいことは分からない.

かくして作られた初版のfactor stencilには, しかし誤りが多数あったようで, 1937年にMichigan大学のJ.D.Elderがそれらを指摘するとともに, そういうものを作るならHollerithカードで作るのはどうかと提案した. Hollerithカードは, 今から50年くらい前, FortranプログラムをパンチしたいわゆるIBMカードで, 20世紀の初めの頃のアメリカの国勢調査の集計用にHerman Hollerithが開発した. 統計機械で読むため, 穴の位置は正確であり, パンチ用の機械も沢山あったはずだ.

Lermerもその案に賛成し, Elderに多くのノウハウを提供した. そういう次第でHollerithカード版のfactor stencilが完成出版されたらしい. 1枚のHollerithカードの, 縦10行横80列の位置を使い, Lehmerの5000個の代りに, 7枚の組で5600個の素数を収めたという. これが完成する直前にLehmerは他界した.

Hollerithカードのfactor stencilがどのようなものであったかに私は興味を持った. Lehmer流に素数を1から始めると1枚目のカードの最後の800番目は6131であり, そこまでの素数について, R=-1とR=2(12月3日のブログの表の1行目と2行目

    3 5 7111317192329313741434753596167717379838997
-1    X     X X     X   X X     X   X     X     X X
 2      X     X   X   X   X   X         X X X   X X

)について, 多分このようであったろうと想像して描いてみたのが下の図だ.




IBMカードは昔のプログラマにはお馴染だが, 最近は見たこともない人が多かろう. この図のように, 横187ミリ, 縦83ミリ程度のカードで, 向きが分かるように, 隅に1ヶ所にコーナーカットがある. 下の拡大図で分かるように, 0から9の行番号は大きい数字で, 0から79の列番号は小さい数字で示す. 本来のIBMカードの列番号は1から80であり, 行は0の上にXとYとがある.





列0, 1の各行に対する素数は

     0  1  2  3  4  5  6  7  8  9
列0  1  2  3  5  7 11 13 17 19 23
列1 29 31 37 41 43 47 53 59 61 67

である. 列2以降に対しても, 素数表があれば, 穴の位置と素数は対応づけられる. 上の表と図の穴の位置とを較べてみて欲しい.

こういう絵を描くと, 実際に穴のあちあカードを作ってみたくなる. 鎌倉のFabLabの田中君と相談し, craftroboを使って試作したR=2のカードが下だ. まだ調整の必要があるらしく, うまく切れていない穴もあるが, 感触は得られた.



プログラムのカードと違い, 穴の数が多くてバイナリカードのようであり, 丁寧に扱わないとすぐに切れそうなのが心配だ. Edlerの作ったカードのセットはどこかに残っていないかしら.

中央値関数の排他的論理和

TAOCPには面白い話が満載である. 今回は多数決関数関連の話だ. この本では, 多数決関数のことを中央値(median)関数というので, そう書くことにしよう.

中央値関数は, 奇数個の引数を大きさの順に並べ, その中央の値を返す(平均値でないことに注意).

    x≤y≤zの時 <xyz>=y.

あるグループで流行っている最中限ゲームみたいなものである.

論理関数としては, x,yを偽か真とする時, 常に反対するか常に賛成する偽か真のサクラの入力を加え, <偽xy>がxとyの論理積, <真xy>が論理和である.

実数値については, <-∞xy>はxとyの最小値, <+∞xy>は最大値になる.

こういう話題がいちおう終わってから,

m>0について,
x₀⊕x₁⊕...⊕x_2m=<¬x₀s₁s₁...s_2m>
ここで, s_j=<x₀x_jx_j+1...x_j+m-1¬x_j+m¬x_j+m+1...¬x_j+2m-1>,
ただし, k≥1について, x_2m+k=x_kとする;

と書いてあった.

つまり, 何変数であっても, その排他的論理和がサクラ(定数)なしの2段の中央値関数で得られるのである. 美しい.

そもそも排他的論理和は単調関数ではないから, 中央値関数1段だけでは得られない. 3変数の排他的論理和, つまり全加算器もパラメトロンでは2段必要であった.

上の見事な式は, TAOCPによると, 量子化学が専門の佐々木不可止(ふかし)君が考えたという.

佐々木君は東大物理の学生時代に, 高橋研によく出入りしていて, 多数決論理に関心があったのは知っていたが, こういうことも考えたのかと驚いた. 北海道大学の化学科の教授を定年退職し, 苫小牧駒沢大学に移られたが, そこも退職の後, 惜しいことに2007年に他界された.

あの長々とした式は眼がくらくらするので, TAOCPにあったようにm=1とm=2について書き直すと

x₀⊕x₁⊕x₂=<¬x₀<x₀x₁¬x₂><x₀x₂¬x₁>>

x₀⊕...⊕x₄=<¬x₀<x₀x₁x₂¬x₃¬x₄><x₀x₂x₃¬x₄¬x₁><x₀x₃x₄¬x₁¬x₂><x₀x₄x₁¬x₂¬x₃>>

となる.

この式をテストしてみよう. ます多数決関数mを定義する. nは否定をとる関数.

(define (m . xs)
 (define (n x) (- 1 x))
 (let ((l (length xs)))
  (if (> (apply + xs) (/ (- l 1) 2)) 1 0)))
(m 1 0 0) => 0
(m 1 0 1) => 1
(m (n 1）(n 1) 1) => 0

すると5変数(m=2)の排他的論理和はこう書ける.

(define (ex x0 x1 x2 x3 x4)
 (let ((s1 (m x0 x1 x2 (n x3) (n x4)))
       (s2 (m x0 x2 x3 (n x4) (n x1)))
       (s3 (m x0 x3 x4 (n x1) (n x2)))
       (s4 (m x0 x4 x1 (n x2) (n x3))))
  (m (n x0) s1 s2 s3 s4)))
(map (lambda (x) (apply ex (n2bs x 5))) (a2b 0 32)) =>
(0 1 1 0 1 0 0 1 1 0 0 1 0 1 1 0 1 0 0 1 0 1 1 0 0 1 1 0 1 0 0 1)

でどうやらよさそうだ.

パラメトロンは通常3入力止りだが, 5入力が許されるとしてこの回路はこう描ける.



これをさらに納得すべく, s₁,...s₄や, その計算の元になっているx₁+x₂-x₃-x₄
などを, x₀=0である前半の16個について出力してみよう. 自己相対関数なので, 半分やればいいことになっている.

xの否定は1-xと計算し, (l-1)/2の2と比較したから, 閾値と較べる値は, 本来は
    x₁+x₂+(1-x₃)+(1-x₄)=2+x₁+x₂-x₃-x₄.

従って, 先頭の2を外して,
    x₁+x₂-x₃-x₄
の和sumについてs₁はsum≤0なら偽(x₀=0を除外したから, 等号が要る), sum>0なら真である.

s₁と添字がmだけ違うs₃は, 足す項と引く項がちょうど反対だから, 和も正負逆で同じ値になる. こういう対を相棒といおう. ところで, 真が奇数個だと, 和が0になることはないので, 相棒同士の真理値は反対になり, x₀がキャスティングボートをとって, exは真になる. すなわち排他的論理和である.

ところが, 真が偶数個だと, どこかに相棒同士で和が0になるものがあり, x₀は0だから, sは2つとも偽になる. 先頭の¬x₀が1を主張しても, 和が0でない他の相棒同士は打ち消し合い, この2個(あるいは4個)の0連合には勝てず, exは0になる.

   + s1 + s2 + s3 + s4 ex
0( 0 0  0 0  0 0  0 0  0)
1(-1 0 -1 0  1 1  1 1  1)
2(-1 0  1 1  1 1 -1 0  1)
3(-2 0  0 0  2 1  0 0  0)
4( 1 1  1 1 -1 0 -1 0  1)
5( 0 0  0 0  0 0  0 0  0)
6( 0 0  2 1  0 0 -2 0  0)
7(-1 0  1 1  1 1 -1 0  1)
8( 1 1 -1 0 -1 0  1 1  1)
9( 0 0 -2 0  0 0  2 1  0)
a( 0 0  0 0  0 0  0 0  0)
b(-1 0 -1 0  1 1  1 1  1)
c( 2 1  0 0 -2 0  0 0  0)
d( 1 1 -1 0 -1 0  1 1  1)
e( 1 1  1 1 -1 0 -1 0  1)
f( 0 0  0 0  0 0  0 0  0)

左端は行数である. 各行のかっこ内の9個の値は, 左からx₁+x₂-x₃-x₄, s₁, x₂+x₃-x₄-x₁, s₂, x₃+x₄-x₁-x₂, s₃, x₄+x₁-x₂-x₃, s₄, exである. 偶数個の場合, どの対で和が0になるかを調べたのが下だ. 赤の-は, x₁とx₂の和がx₃とx₄の和に等しい(つまりs₁s₃が0になる)場合, 青の-は, x₂とx₃の和がx₄とx₁の和に等しい(つまりs₂s₄が0になる)場合である.


0(0 0-0-0-0-)
1(0 0 0 0 1 )
2(0 0 0 1 0 )
3(0 0 0-1 1-)
4(0 0 1 0 0 )
5(0 0-1-0-1-)
6(0 0-1 1-0 )
7(0 0 1 1 1 )
8(0 1 0 0 0 )
9(0 1-0 0-1 )
a(0 1-0-1-0-)
b(0 1 0 1 1 )
c(0 1 1-0 0-)
d(0 1 1 0 1 )
e(0 1 1 1 0 )
f(0 1-1-1-1-)


和が0になる対は本当にあるか. こういう図を描いてみた.



いま真というか1の数は偶数である. s₁とs_m+1の範囲で1の数が同じであれば, 一発で和が同じところは見つかったわけでハッピーだ.

一方, s₁がオール0でs_m+1がオール1なら, 左の密度は0, 右の密度は1. そこで図の下へs₂, s₃のように加算の範囲をずらしていけば, どこかで2つの範囲の1の数が一致する, 同じ密度になる場所があり, そこで赤や青の-が引ける対が見つかるはずである(平均値の定理を習ったのは高校生の時だっけ).

こういうわけで, 偶数の場合はexが偽になるのである.

「なるほど」であった.

実はサクラを使うと, どんな関数でも2段の中央値関数で作れるが, そのことはいずれ.

平方根の連分数展開

Norman Lehmerの論文に, factor stencilはNの連分数展開でも使えると書いてあった.

その始めの方. まずx₁=√N=q₁+1/x₂とおく. 但し, q₁はNの平方根の最大整数. 同様にx₂=(√N+P₂)/Q₂=q₂+1/x₃のようにする. (ははぁ. TAOCPの連分数展開のUとVがPとQになっているわけだ.)

順に進んで, x_n=(√N+P_n)/Q_n=q_n+1/x_n+1になるというまでは明快だが, その後に
「従って
P_n=q_n-1Q_n-1-P_n-1,
Q_n=q_n-1(P_n-1-P_n)+Q_n-2
がexpediouslyに得られる」
と書いてある.

上の方の式は11月30日のブログの(4)と(5)の間にあるU_n=A_nV_n-U_n-1と同じだが, 下の方の式はどう導くか.

今回のブログはその計算法である.

(0)はx_nの標準形である. (1)はx_n-1で, そのx_nに(0)を代入し, 前回のように計算を続けると(1)の最後になる. (1)の最初と分母を等しいとしたのが(2)で, ただちに(3)と書ける. 分子の√Nの後に続く項を等しいとしたのが(4)で, 途中で(3)を利用する. P_nとP_n-1を交換したのが上の式(5)だ.

一方, (3)の添字を1減らしたのが(6)で, (3)-(6)が(7),(8)である. 1/Q_n-1は, (5)を移項した(9)から(10)のように得られ, これを(8)に代入すると(11)を経て, 下の方の式(12)になるのが分かる.

漸化式はP₁=0, Q₁=1から始まるが, 問題は, Q_n-2があるので, Q₀が必要なことだ.

実はQ₀=Nなのだが, この計算は次の通り.

これで, 連分数のP,Q,qが次々と求まることになる. なんだかまた久しぶりにこのような計算をした.

素因数探し

11月6日の本ブログにある自転車チェーンの篩を作ったのはHenry Lehmerだ. その父親のNorman LehmerもCalifornia大学Berkeley校の数学の教授であった. 父のLehmerが素因数探しの道具として, Factor Stencilを考案したという話がある.

それがどういうものかが分かってきた. 今回はその話をしよう.

aが素数pの平方剰余であるとは, x²=a mod pとなる何かのxがあることなのは周知のとおり. こういうa(≠0)の時, Legendreの記法で(a/p)=+1とする. (htmlなので, 分子分母が横並びだが, 通常は上下に書く.)

xがないときは, aを平方非剰余といい, (a/p)=-1とする.

p=7の時, 剰余は(0,1,2,4), 非剰余は(3,5,6)である. ここで, 0を除外し, 非剰余同士を掛けると, 3*3=2 mod 7, 3*5=1 mod 7, 5*5=4 mod 7のように剰余になり, 剰余と非剰余を掛けると2*3=6 mod 7のように非剰余のまま, 剰余同士を掛けると4*4=2 mod 7で剰余であり, +1と-1で掛算になる.

ところで, ある整数Nを素因数分解する場合, 任意の整数aについて, a²-N=Rは, Nがpで割れるなら, mod pで考えると, Nの項はないのと同じだから, Rはこのpの平方剰余になる.

aを√Nの程度にとると, Rは小さくなり, 平方剰余の乗算に規則でもっと分解すると, さらに小さい剰余か非剰余が得られ, それから素因数の形が決るという.

たとえば, 平方剰余に-1があるなら, すなわちp-1なら, 素因数は4n+1, 2があるなら, 素因数は8n±1の形, というふうに, 素因数はlinear formになっている. それを組み合せて, 素因数の探索の範囲を狭めたい. ところが, 平方剰余が段々大きくなってくると, linear formの形も複雑になり, 式の形を決めにくくなる. Normanの論文によると, Rが113と199の場合, formの形は11088にもなるそうだ.

そこで, 各平方剰余について, 候補になる素数の表をあらかじめ作っておこうというのが, この発想である.

小さい範囲で, テストしてみよう. 下がその剰余(縦方向)と素数(横方向)の表である. 10月26日のブログの図の素数のところを抜き出したものと思えばよい.

    3 5 7111317192329313741434753596167717379838997
-1    X     X X     X   X X     X   X     X     X X
 2      X     X   X   X   X   X         X X X   X X
 3        X X     X     X     X   X X   X X   X   X
 5        X     X   X X   X       X X   X   X   X
 7  X           X   X X X     X X X           X
11    X X       X       X   X   X           X X X X
13  X         X   X X       X   X   X       X
17          X   X           X X X X   X       X X
19  X X       X       X           X X X X X X

Normanの素数表には1があるというので, この左に1と2の列があるかも知れないが, 除数としては不要なので, とりあえず3からの表だ. 7の下を見ると, 2と11にXがあり, 2と11=4 mod 7が平方剰余であることを示す. 他の素数についても同じだ. たとえば,

(quadres 97)=>(0 1 2 3 4 6 8 9 11 12 16 18 22 24 25 27 31 32
 33 35 36 43 44 47 48 49 50 53 54 61 62 64 65 66 70 72 73 75
 79 81 85 86 88 89 91 93 94 95 96)

のうち, 19までの素数, 2,3,11と, 96(= -1 mod 97)があるから-1とにXがある.

また, 横向きに眺めると, -1の行は, 5,13,17,29,37,41,53,61,73,89,97で, すべて4n+1. 2の行の7,23,31,47,71,79は8n-1, 17,41,73,89,97は8n+1である.

さて, これを使い, N=2279(=43x53) の素因数を探してみる. (sqrt N) => 47.738873049120045 だから, aを48から順に増やしながら, a²-Nを計算し, その素因数分解もしておく.

 (map (lambda (a) (factorize (- (* a a) n))) (a2b 48 80))
=> ((5 5) (61 2) (17 13) (23 7 2) (17 5 5) (53 5 2) (13 7 7)
 (373 2) (857) (97 5 2) (31 7 5) (601 2) (1321) (103 7 2)
 (313 5) (13 13 5 2) (79 23) (139 7 2) (67 31) (17 13 5 2)
 (67 7 5) (73 17 2) (2621) (1381 2) (83 7 5) (61 5 5 2)
 (139 23) (239 7 2) (269 13) (73 5 5 2) (761 5) (283 7 2))

ここに並んでいるのは, 掛けると+1になるもの同士である. (5 5)からは, 5が+1かも, -1かも知れないということが分かるだけだ. (61 2)では, 61が上の表の剰余にないから, 無視. (17 13)は, とりあえず13と17は+1と+1か, -1と-1が分かっただけ. その次の(17 5 5)はラッキーだ. 5と5は打ち消すから, 17が+1と判明した. したがって, (17 13)から, 13も+1と決る.

(13 7 7)からも13が剰余と確信出来る. さらに進むと, (13 13 5 2)や(17 13 5 2)から, 2と5は+1同士か-1同士か, いずれにしろ同じ仲間ということが分かるが, 役にたつかどうかは不明.

しかし, 13と17が得られたので, 13と17の行にXのある素数を探すと, ブラボー! 43と53が見つかった. ついでだが, 43と53の列では, 2と5は空白で, 同じ組なことが示せる. いまは1組だけが見つかったが, 通常は候補がたくさん得られ, 実際に割ってみる必要がある.

Norman Lehmerが1914年に出版したFactor Stencilは, 剰余が±238まで, 素数が48593までの素数表であって, 各剰余ごとに1枚で, Xの代りに素数の位置に孔をあけてある. 今の例では13と17の紙を取り出してきて重ねると, 条件に合う素数のところから光が洩れる仕掛けであった. 当時としては大変な労作であった. いまでは, MacBookで簡単に実験が出来, ありがたい時代だ.

平方根の連分数展開

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であって,
v_n+1=(d-u_n²)/v_n,
a_n+1=floor((√d+u_n)/v_n+1),
u_n+1=a_n+1v_n+1-u_n
を繰り返す. aに√dがあるが, 整数部分をとるので下の方はいらない. vもuも整数のままで, 誤差が入らないから, いつまでも計算が続く.

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

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



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

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