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.

0 件のコメント: