2018年8月22日水曜日
e1000桁の計算
60年くらい前によく行われていた, Napier's constant eの値を1000桁計算するパラメトロン計算機PC-1のプログラムについてのファイルがあったので, 私のブログに載せておくことにした.
自然対数の底 e の値を小数点以下1000桁まで計算することは,1960年代の計算機で流行した. Taylor展開
e = 1 + 1/1! + 1/2! + ... + 1/n!
で計算するのだが, n = 460くらいで1000桁の精度がある.
実際には上の式で計算するのではない.
e = 1 + 1/1(1 + 1/2(1 + 1/3(1 + ... +1/459(1 + 1/460)...)))
を使う. さらに例えば4項ずつまとめて計算することを考える.
Qa-4 = 1 + 1/(a - 3)(1 + 1/(a - 2)(1 + 1/(a - 1)(1 + 1/a・Qa)))
a = 460, 456, 452, ... 0
Q460 = 1
e = Q0
これはさらに簡単になって
Qa-4 = {a(1 + (a - 1)(1 + (a - 2)(1 + (a - 3)))) + Qa}
/{a(a - 1)(a - 2)(a - 3)}
しかし 460×459×458×457 = 44192868840 となり, 36ビットのPC-1のレジスタの容量235 = 34359738368には収まらない. 連続する4個の整数の積は4で割れるから, 2ビット右シフトして計算すれば, この問題は解決する.以下のプログラムはそのようにしてeの値を計算している.
次は10進数での出力である. 小数を10進法で出力するには, 全体を10倍し,整数部分にでた数を印字する. これでは乗算に時間をかけて1桁しか印字できないので,以下のプログラムでは, 小数部分に一度に1010をかけ,10桁分の整数部を作り出す. また多倍長の乗算を10区間に分割し,各区間毎に10桁の整数を1桁ずつ印字するようになっている.
作業場所は10進で1000桁なので, 3300ビット必要. PC-1の36ビットの長語を100語使う
(300〜498 300番地の方が小数点側).
作業場所は10進で1000桁なので, 3300ビット必要. PC-1の36ビットの長語を100語使う (300~498 300番地の方が小数点側).
0h=80,0i=100,0p=180,0n=272,0m=300,
70:
0 it
1 jl h 作業場所クリアへ
2 it
3 jl i 計算プログラムへ
4 it
5 jl p 印字プログラムへ
6 0
0h:
0 a 8n ┐
1 x 13r ┘ plant link
2 p 2n
3 x 7r
4 s 3n
5 z 11r
6 l 35 アキュムレータを0に
7 tl( ) 作業場所をクリア
8 p 7r
9 a 8n
10 jl 3r
11 pl 16n Q460を1に
12 tl m
13 jl( )
0i:
0 a 8n ┐ e1000桁計算ルーチン
1 x 3r ┘ plant link
2 pl 18n ┐
3 zl( ) ┘ a = 0なら帰る
4 sl 16n a - 1
5 tl 20n
6 sl 16n a - 2
7 tl 22n
8 sl 16n a - 3
9 tl 24n
10 vl 22n ×(a - 2)
11 ll 34
12 vl 20n ×(a - 1)
13 ll 35
14 vl 18n ×a
15 ll 34 a(a-1)(a-2)(a-3)
16 tl 26n を26nLへ
17 p 2n
18 x 22r
19 x 25r
20 s 3n
21 z 29r
22 ql( ) ┐
23 r 33 │ Qを4で割る
24 ll 33 │
25 tl( ) ┘
26 p 22r
27 a 8n
28 jl 18r
29 pl 24n a - 3
30 al 16n + 1
31 vl 22n x(a - 2)
32 ll 35
33 al 16n + 1
34 vl 20n x(a - 1)
35 ll 35
36 al 16n + 1
37 vl 18n x a, aは4の倍数なので
38 ll 33 最後の2ビットは0
39 al m Rレジスタは0
40 tl m
41 p 2n
42 x 46r
43 x 48r
44 s 3n
45 z 52r
46 ql( ) ┐ 多倍精度/単精度
47 dl 26n │ a(a-1)..(a-3)で割る
48 tl( ) ┘
49 p 46r
50 a 8n
51 jl 42r
52 pl 24n
53 sl 16n
54 tl 18n
55 jl 3r
0p:
0 a 8n ┐ 出力ルーチン
1 x 68r ┘ plant link
2 q 16n
3 p 4n
4 x 8r
5 x 12r
6 s 5n
7 z 16r
8 pl( ) ┐
9 wl 12n │
10 tl 28 │ 1回目の1010倍
11 ql 28 │
12 tl( ) ┘
13 p 8r
14 s 8n
15 jl 4r
16 p 4n ┐
17 tllt │ 2. と印字
18 s 70r │
19 tllt ┘
20 p 6n 20行用カウンタ設定
21 x 9n
22 p 7n 5ブロック用カウンタ設定
23 x 10n
24 z 65r
25 pl 28
26 tl 68
27 p n 1010倍10回カウンタ
28 x 11n 初回は0に設定
29 q 16n
30 p 4n
31 x 39r
32 x 43r
33 s 5n
34 z 62r
35 p 11n
36 z 47r 10回カウンタを調べる
37 s 17n
38 x 11n
39 pl( ) ┐
40 wl 12n │
41 tl 28 │ 1010倍
42 ql 28 │
43 tl( ) ┘
44 p 39r
45 s 8n
46 jl 31r
47 ql 68 1010倍10回終り
48 tl 28 1文字 印字へ
49 p n
50 dl 14n 109で割る
51 l 18
52 a 56r
53 x 54r
54 p ( ) 文字コードを取る
55 tllt 印字
56 q 2n アドレス部分は52で使う
57 vl n
58 ql 28
59 tl 68
60 p 1n
61 jl 37r
62 p 10n
63 s 17n
64 jl 23r
65 p 9n
66 s 17n
67 z ( ) link
68 x 9n
69 p 8n
70 l 12 CR
71 tllt
72 l 2 LF
73 tllt
74 r 1 SP
75 tllt
76 tllt
77 jl 22r
0n:
0 0 ┐
1 10 ┘ 10(長語)
2 229676 (0) 300 以下文字コード
3 233972 (1) 500 と定数
4 197106 (2) 498
5 172332 (3) 300
6 135188 (4) 20
7 217093 (5) 5
8 245762 (6) 2
9 180224 (7) ( ) カウンタ
10 143360 (8) ( ) カウンタ
11 184320 (9) ( ) カウンタ
12 38146 ┐
13 254976 ┘ 1010
14 3814 ┐
15 182784 ┘ 109
16 0 ┐
17 1 ┘ 1 (長語)
18 0 ┐
19 460 ┘ a (長語)
この後27nまで作業場所
jl70.
登録:
コメントの投稿 (Atom)
0 件のコメント:
コメントを投稿