1970年頃のことだ. HITAC 10, FACOM R, MACC 7など, 1語が16ビットの4K
語のミニコンが多く使われていた. これらのミニコンは, 正面に押しボタ
ンのようなスイッチが並んでいて, 電源を入れた後, スイッチを操作し,
数語の「初期入力プログラム」を手動で入力する. プログラム全体が入力
し終わると初期入力プログラムが動き出し, 紙テープのプログラムを読み
込めた.
ところがこのメーカー製プログラムの手の入力が面倒で, 誰もがもっと短
いプログラムに出来ないか工夫を凝らし, 徐々に短いプログラムが出現し
た. これが一番短いと思っていると, 他から更に短いのが出来たという話
が来る. するとまた知恵を絞り, その語数のプログラムが作れるのであっ
た.
つまり「何語で出来る」という情報が重要で, それを聞くとその語数のプ
ログラムに向って注力できる.
最近ツイッターを眺めていたら, Automatic Computer@nt_a_c さんの
Automatic Computer@nt_a_c·3月18日
PC-1で円周率、400桁は出せた。もう数日早くできればちょうど良かったのだが。後はどのくらいコードを削れるか。
Automatic Computer@nt_a_c·3月24日
円周率、できた。1000桁、印刷と合わせて3回読み込み、合わせて40分くらい。https://automaticcomputer.github.io/PC1onUnity/pi....
というツイートがあった. 60年も前の計算機 PC-1でプログラムを書く奇
特な人がまだいるかと驚く一方, どうすればあの容量のメモリーでπが1000桁計算
できるか私も考え始めた.
πの1000桁については, このブログでも2012年11月5日に書いている. 簡単に計算するには
やはりMachinの式でやることになろう.
π/4 = 4 tan
-1(1/5)- tan
-1(1/239)
tan
-1(1/5)=1/5-1/(3·5
3)+1/(5·5
5)
-1/(7·5
7)+...
tan
-1(1/239)=1/239-1/(3·239
3)+1/(5·239
5)-1/(7·239
7)+...
PC-1は1語が18ビット. 計算は2語をつなげた長語36ビットを使う. 十進10
桁で1長語を使うから, 1000桁計算には100長語が必要で, eの1000桁では
300番地から500番地の100長語を作業場所に使っている. 512短語のPC-1で
は100長語はどうしても2組しかとれないから,それで上の計算をしなけれ
ばならない.
普通に考えると16×tan
-1(1/5)を計算するには, atanを置く場所tsを0にし,
多倍長の16/5を作る. これをssとする. それを1で割り(その商をqsとする), qsをtsに足す.
ssを25で割りssにし, そのssを3でわりqsとし, tsからqsを引く. ...
これをやるとtan
-1(1/5)を計算するだけで100長語が3組必要になりすでに
破綻する.
tsはtan
-1(1/5)とtan
-1(1/239)で共用出来る. ssは
1/5が済んだ後1/239用に使えばよい.
残る手段はqsを作ってから, 下の桁からtsに足したりtsから引いたりするのを
なんとかすることだ.
我々が小学校以来習った加減算は下の桁から始めるものであった. それに対して
上級になって習った算盤では上の桁から始めた. 小学生のころ, 上の桁から加減算
が出来るのが不思議であったが, そのうち繰上げをローカルに処理しながら下へ進んで
いるからだと納得した. さらに一度繰上げが通過した桁に, その後から再び繰上げ
が来ることはないことも証明できた. 基数が10から2^35と大きくなると,
繰上げが何桁も伝搬する確率は小さくなるが, 無視できない. い
要するにqsの結果はいらない. 各々の桁の商はqsに足すか引くかし, 剰余は次の
桁の除算に使うと作業場所は要らなくなり, なんとかなりそうだ.
そこでSchemeを使ってシミュレートしたのが次のプログラムである.
lnは長語数. b35, d10は二進35桁, 十進10桁の定数. tsはatanの値を蓄積する
場所, ssはatan(1/z)の計算で1/z, 1/z^3, 1/z^5,...を保持する場所である.
(zは5と239)
ldivはssをd = z^2で次々割る. ldiv+- はssをn(1,3,5,..)で割りながら,
部分商をtsに足すか引くかする.
(define ln 100)
(define b35 (expt 2 35))
(define d10 (expt 10 10))
(define (ldiv ss d)
(let ((r 0) (n 0))
(do ((i 0 (+ i 1))) ((= i ln))
(set! n (+ (fix:lsh r 35) (list-ref ss i)))
(list-set! ss i (quotient n d))
(set! r (modulo n d)))))
(define (ldiv+- ts ss ff d)
(let ((r 0) (n 0) (s 0) (c 0))
(do ((i 0 (+ i 1))) ((= i ln))
(set! n (+ (fix:lsh r 35) (list-ref ss i)))
(set! r (modulo n d))
(set! s ((if ff + -) (list-ref ts i)
(quotient n d)))
(list-set! ts i (fix:and s (- b35 1)))
(set! c (fix:and s b35))
(do ((j (- i 1) (- j 1))) ((= c 0))
(set! s ((if ff + -) (list-ref ts j) 1))
(list-set! ts j (fix:and s (- b35 1)))
(set! c (fix:and s b35))))))
(define (lmul ns m)
(let ((e 0) (p 0))
(do ((i (- ln 1) (- i 1))) ((= i 0) e)
(set! p (+ (* (list-ref ns i) m) e))
(list-set! ns i (modulo p b35))
(set! e (quotient p b35)))))
(define ts (make-list ln 0))
(define ss (make-list ln 0)) (list-set! ss 0 80)
(do ((n 1 (+ n 2)) (ff #t (not ff))) ((> n 1433))
(ldiv ss (* 5 5))
(ldiv+- ts ss ff n))
(define ss (make-list ln 0))(list-set! ss 0 (* 4 239))
(do ((n 1 (+ n 2)) (ff #f (not ff))) ((> n 423))
(ldiv ss (* 239 239))
(ldiv+- ts ss ff n))
(do ((n 0 (+ n 1))) ((= n 100))
(display (list (lmul ts d10))))
この出力の始めの方は次のようだ.
1415926535 8979323846 2643383279 5028841971
6939937510 5820974944 5923078164 628620899
8628034825 3421170679 8214808651 3282306647
938446095 5058223172 5359408128 4811174502
8410270193 8521105559 6446229489 5493038196
...
10^10進法の各桁の後に空白があるのは, 十進の一番上が0なのが分るためだ.
これで一応は計算出来るが, 考えてみるに, ssは徐々に小くなっていく.
その途中にずーっと先頭から除算を始めるのはやはり無駄である. つまり
除算を開始する位置を制御すれば計算は早く終わる. 下はそのように修正
したプログラムで, 変数bが除算開始位置を覚えている.
(define b 0)
(define (ldiv ss d)
(let ((r 0) (n 0) (q 0))
(do ((i b (+ i 1))) ((= i ln))
(set! n (+ (fix:lsh r 35) (list-ref ss i)))
(set! q (quotient n d))
(list-set! ss i q)
(if (= q 0) (set! b (+ i 1)))
(set! r (modulo n d)))))
(define (ldiv+- ts ss ff d)
(let ((r 0) (n 0) (s 0) (c 0))
(do ((i b (+ i 1))) ((= i ln))
(set! n (+ (fix:lsh r 35) (list-ref ss i)))
(set! r (modulo n d))
(set! s ((if ff + -)
(list-ref ts i) (quotient n d)))
(list-set! ts i (fix:and s (- b35 1)))
(set! c (fix:and s b35))
(do ((j (- i 1) (- j 1))) ((= c 0))
(set! s ((if ff + -) (list-ref ts j) 1))
(list-set! ts j (fix:and s (- b35 1)))
(set! c (fix:and s b35))))))
(define (lmul ns m)
(let ((e 0) (p 0))
(do ((i (- ln 1) (- i 1))) ((= i 0) e)
(set! p (+ (* (list-ref ns i) m) e))
(list-set! ns i (modulo p b35))
(set! e (quotient p b35)))))
(define start (real-time-clock))
(define ts (make-list ln 0))
(define ss (make-list ln 0)) (list-set! ss 0 80)
(set! b 0)
(do ((n 1 (+ n 2)) (ff #t (not ff))) ((> n 1433))
(ldiv ss (* 5 5)) (ldiv+- ts ss ff n))
(define ss (make-list ln 0)) (list-set! ss 0 (* 4 239))
(set! b 0)
(do ((n 1 (+ n 2)) (ff #f (not ff))) ((> n 423))
(ldiv ss (* 239 239)) (ldiv+- ts ss ff n))
(do ((n 0 (+ n 1))) ((= n 100))
(display (lmul ts d10)) (display " "))
上のプログラムをPC-1のプログラムにしたのが次だ. 作業番地tsはeの1000桁と同様,
300番地以降を使う. プログラムの最初にあるインタールードでクリアしておく.
ssは0番地から200番地までを使う. つまりイニシアルオーダーR0は壊す.
計算だけで記憶場所を使うので, 結果tsの出力は別のプログラムにする.
まだコメントが少ないが, 徐々に追加するつもりだ.
0p=200,
0 0p=196,0p:p2, インタールード 0をaccへ
1 tl(300), 作業領域ts(300から498)をクリア
2 p1r,
3 a40, 2を足す
4 x1r,
5 s8r, 終了か
6 z40, R0へ戻る
7 jlr, ループする
8 500,jlp.
0 0p=196,0n=276,0m=496,0p:a7n,
1 x16r, atan サブルーチン 帰り番地を入れる
2 p4n, 作業場所ssをクリア
3 x7r,
4 s19n, 終了?
5 z11r,
6 p4n, 0をaccに
7 tl( ), 作業場所へ
8 p7r, アドレスを増やす
9 a7n,
10 jl3r,
11 pl10n, ssの最初の語に数を入れる
12 tl0,
13 pl4n, atan 計算ループの回数
14 tl2n,
15 sl8n, その回数が0になる迄
16 zl( ), 帰り番地
17 q4n, z2の除算ループ. 最初の剰余をRregに
18 p22n,
19 x23r, ssの語のアドレスを順に作る
20 x25r,
21 s19n, 終り?
22 z34r,
23 ql( ), 被除数をaccとRregに
24 dl12n, z^2による除算
25 tl( ), 商を格納
26 sl4n, 商が零かテスト 1を引き零なら負にする
27 kl31r,
28 p22n, 除算開始の位置b
29 a7n, bを2増やす
30 x22n,
31 p23r,
32 a7n, ssのアドレスを増やす
33 jl19r,
34 q4n,
35 p22n, ssをnで除し, tsに加減する
36 x43r,
37 a20n,
38 x46r,
39 x49r,
40 x52r,
41 s21n, ループの終か
42 z70r,
43 ql( ), 被除数を作り
44 dl2n, nで割る.
45 tln, 商を一旦しまう
46 pl( ), ts
47 aln, aln<->sln 商を足すか引く
48 k51r, sign bitが立っていたら繰上げ処理へ
49 tl( ), 商をしまってこの語は終り
50 jl67r,
51 cl14n, sign bitを消して
52 tl( ), しまう
53 p46r, 繰上げ処理
54 s7n,
55 x58r, 繰上げ先の語のアドレス
56 x61r,
57 x64r,
58 pl( ),
59 al4n, aln<->sln 繰上げを足すか引く
60 k63r, sign bitがあれば繰上げ
61 tl( ), 繰上げがなければ処理は終り
62 jl67r,
63 cl14n, sign bitを消す
64 tl( ), しまう
65 p58r, 繰上げ先のアドレスを進める
66 jl54r,
67 p43r, ssの語のアドレスを進める
68 a7n,
69 jl36r,
70 p79r, バイナリスイッチ切替え
71 b47r,
72 t47r,
73 p79r,
74 b59r,
75 t59r,
76 pl2n, nを増やす
77 al6n,
78 jl14r,
79 i, a xor s = i
0 0n:,, データ領域
2 ,,
4 ,1,
6 ,2,
8 ,1435,
10 ,80, =5*16 atan(1/5)で使う
12 ,25,
14 131071,262143, マスクパターン
16 425,956, 956=239*4
18 57121,194, 57121=239*239
20 300,494,
22 0,
0 0m:it, 主ルーチン
1 jlp, atan(1/5)
2 p16n, atan(1/239)用に
3 t9n, データを入替え
4 p17n,
5 t11n,
6 p18n,
7 t13n,
8 p4n,
9 t22n,
10 p4n,
11 t22n,
12 it,
13 jlp, atan(1/239)
14 0, 計算終り 停止する
jlm.
このプログラムは私のPC-1シミュレータに追加してある. シミュレータを開き, R0
を読み込み, Place Tapeをクリック. メニューの下方にあるpi1000cmpを
クリックすると計算のプログラムが設定される. そこでFree Runをクリック;
さらにInit Startをクリックするとプログラムが走り出す.
計算は私のMacBook Pro環境では約3分半で終り, シミュレータは停止する.
結果の出力には作業場所
を同じにしてあるので, eの1000桁の出力プログラムを多少修正して使う.
計算プログラムがR0領域を破壊したから, 出力プログラムの読込み前にもう一度
R0をロードする. 後はPlace TapeからPi1000prtを選び, そのプログラムを
読込み起動すればπの1000桁が出力される.
久し振りにPC-1のプログラムを書いたが, インタプリタで必要な情報を取り出し
たり, Schemeのプログラムで予備テストが出来たり, 便利な世の中である.
