2017年4月26日水曜日

したたり算法

この話題の前回のブログでは, 円周率を表すmixed radixの式を示し, 右の項から 正規化しながら十進の小数点以下の数字を求める話を書いた.

通常の計算ならそうするのだが, したたり算法では左から計算したい. それも 可能であって, 次のように計算が進行する.

下のtex出力の図で, 左上の式を見てほしい. 左端のvとuにすべての情報を 記憶させながら, 左側(外側)からかっこをはずす. 話の発端となった円周率の式 は, 左側の2行目のものである. 2+で始まっていたが, 外から1を掛けて標準形に してある.


最初の仕事は左かっこの外にあるvを, かっこ内の2つの項に掛けることだ. 今は 1だから左の3行目のように, vが消えた以外はなにも変らない.

次は先頭に来たuに相当する2を, 1/3の右のかっこの中にいれることである. それには この2に1/3の逆数, 3/1を掛ける. そうして出来た2×(3/1)+2をまとめて8 にする. 下から3行目のように, 1/3(8+ ... が出来て, 標準形であるが, かっこが 外側から1組減った.

また同じように1/3をかっこ内の8と2/5に掛ける. そうして出来た8/3を 2/15の逆数を掛けてかっこ内に入れる. すると右側の上の行のようになる. 後は 一々書かぬが, 同様に進むのである.

これを手でやるのは面倒なので, 以後は計算機にやってもらう. かっこ内は左からu とvと, 標準形の式にあったiの値である. 最後の項が5に相当するところ(4行目)までが 上の出力に見られる.

(8 1/3 2)
(22 2/15 3)
(160/3 2/35 4)
(122 8/315 5)
(1352/5 8/693 6)
(8818/15 16/3003 7)
(8832/7 16/6435 8)
(18782/7 128/109395 9)
(356984/63 128/230945 10)
...
これを更に進めると
(864779799556983658/40970475 67108864/450883717216034179 31)
のようになり, このuとvを掛けると3.1415926533011596になる. 小数点以下9 桁目まで合っている.

これが3.1415...なのは計算機が十進法にしてくれたからで, したたり算法では 自分で計算しなければならない. その過程が次の出力である.


最初の行は, 先程の長々とした分数のuとv, それとiから計算した残りの項 (31/63(2+...))である. u*vの整数部分を計算してみると
(floor->exact
 (*  67108864/450883717216034179
   864779799556983658/40970475)) => 3
と3になるので, 次の2行にわたる式のように, 先頭に3を括りだす. そしてかっこの中から3にvの逆数を掛けたものを引いて辻褄を合せる. かっこ 内を計算したものが最後の行である.

この最後の行のvとuと十進にしたいので10を掛けると, 整数部分に1が得られる. 3.14...の1のことだ.
(floor->exact
(* 10 67108864/450883717216034179
  2615629766097082765349437/2749482034790400)) => 1
ここでもTexで書くのは止め, 計算機にやってもらう. 次々のuとvと出力される 数である.
(2615629766097082765349437/2749482034790400
   67108864/450883717216034179 3)
(1536675519372845944725869/5498964069580800
   671088640/450883717216034179 1)
(58841914244318110336667/5498964069580800
   6710886400/450883717216034179 4)
(437921482322098289538739/109979281391616000
   67108864000/450883717216034179 1)
(136926162079932661882877/219958562783232000
   671088640000/450883717216034179 5)
(196056880918257839392441/10997928139161600000
   6710886400000/450883717216034179 9)
(60341900506756319941901/13747410173952000000
   67108864000000/450883717216034179 2)
(196925612577461046092237/549896406958080000000
   671088640000000/450883717216034179 6)
(48785647745580267174347/2199585627832320000000
   6710886400000000/450883717216034179 5)
(222531979586221607133547/109979281391616000000000
   67108864000000000/450883717216034179 3)
(8569388169424319751667/1099792813916160000000000
   671088640000000000/450883717216034179 3)
今日の計算はまずかっこを取り払い, つまりかっこの次の2を消し, この処理をi=31に なるまで進め, 今度は十進の各桁と順次に括りだすものであった. 大体は 左から右へ処理できたが, 本来のしたたり算法では, 出力できる状態に なると同時に出力するわけで, 次回はその話をしたい.

2017年4月13日木曜日

したたり算法

ことし(2017年)3月にオックスフォード大学のGibbons君が研究所のセミナーで話をしてくれることになった. 彼の話は面白いものもあるが, 難解なものも多いので, 少し昔の研究だが, Spigot Algorithmで円周率を計算する話をして貰った.

Spigot Algorithmとは, 円周率πや自然対数の底 e のような長い小数を, 頭の方から次々と出力する計算法である. これはちょうど水道の栓から水が一滴一滴したたり落ちる様子に似ているので, 水栓という意味のSpigotを使い, Spigot Algorithmという. 水栓算法というのが直訳だが, 殺風景なので, ここでは「したたり算法」といおう.

日本で電子計算機が使えるようになった1960年頃, あまり速くもなく, 記憶装置も少なかった当時の機械で, eの値を計算することが各センター間で流行した. 東大高橋研のパラメトロン計算機PC-1の記憶装置(18ビット512語)と計算速度(加算400マイクロ秒)では, eは小数点以下1000桁くらいがほどほど計算できる限界であり, eの多項式展開の460の階乗分の1くらいまで使えば1000桁の精度が得られると分ってそういう計算をしていた. まず十進1000桁に相当する値を二進で計算し, その後二進十進変換して2.7182...を出力した.

このように通常は何桁計算するかを初めに決めるのである. これに対してしたたり算法では普通には得ようとする桁を前もっては決めない. eの計算では, 長い時間かけて計算が進行し, 暫くは出力が現れないが, したたり算法ではいきなり出力が始まって見ていて楽しい. しかしそうなにもかも上手くいく筈はなく, 段々と疲れてきたが如く出力が遅くなる傾向がある.

前口上はこのくらいにして, Gibbons流の算法を何回かに分け, 少しずつ説明したい.

eの計算には e=1+1/1!+1/2!+1/3!+... を使うのが常識である. 解析概論でもこれを使ってeの値を計算してみせる(1953年版 p.75).

一方円周率には公式が山のようにあって, どれを使うか迷うが, その昔ENIACで2035桁まで計算した時はMachinの公式

π/4=4 arctan 1/5-arctan 1/239
arctan x=∑n=0(-1)nx2n+1/ (2n+1)

を使ったといわれている. (城,牧之内「計算機械」p.105.)

さて式が徐々に複雑になるので, 以後はTexで書いたものを見ていただく. Gibbonsの計算の式は下のTex出力の一番上のものである. これはeの式のように簡単には得られない. Wolframのページを見ると, その下のarctan xの式が見付かる. この式のxに1を代入すると, arctan 1 は45度つまりπ/4になる. xが1なら, x2n+1は1になり, (1+x2)n+1 は 2n+1になる.

πにするにはこの式を4倍にするわけで, それには22nを22n+2にする. 一方分母には2n+1があるので, 結局2の羃は2n+1になる. それが途中のhttp:...の行の下の式だ. このnに0, 1, 2, 3を代入して書いてみるとその下のようになり, 各項の共通の因子を前に出すと, その下の括弧の式になる.

ここでn=4の項を書くと, 一番下のようになる. 左にある分数は階乗の展開で, これをみると, 前に現れる共通因子 1/3, 2/5, 3/7は赤線で消したように消える. また青線のように分母の8, 6, 4, 2は分子の4, 3, 2, 1と右の5個の2のうちの4個とキャンセルし, 結局右辺にある4/9×2が残り, 一番上の式になった.

また次のTexの出力を見てほしい. 一番上は円周率π=3.1415...の十進法を強調して描いたものである. 1桁下る度の1/10の重みで足されることを示す.

その次は前の式で, 1桁下る度のその桁の値は常に2であることを示す. ただし重みは10ではなく, 3, 5/2, 7/3, ...と変っていく. 十進法をdecimal radixというが, 各桁で重みが変るのはmixed radixという.

mixed radixというと意外かも知れないが, 1日は24時間, 1時間は60分, 1分は60秒, 1秒は1000ミリ秒とか, 1里は36町, 1町は60間, 1間は6尺, 1尺は10寸とかなど身近になくはない.

πの場合はmixed radixではあるが, i番目の重みが(2i+1)/iと決っているのが有難く, その各桁の値が2に固定されているのが味噌である.

従ってもうπの値は分ったようなものだが, やはり十進でないと気がすまないから, mixed radixからdecimal radixへと基数変換(radix conversion)する必要がある.

上の出力で, Mixed radix to decimal conversionとある下が基数変換である. 最初はπの展開式を第5項までとったものだ. それでまずこれを10倍する. つまり 各桁の値を10倍する. それが次の行である.

十進の場合の各桁が0から9までの値しか取れないように, このmixed radixでも, 第i桁は2iまでしかとることができない. 従って下の桁から, 2iを越えていれば, 正規化する必要がある. (20)の下の方に1,9と見えるのは, 20を分母の11で除したときの商が1, 剰余が9であることを示す. つまりこの桁は正規化すると, 1が繰上がって9になるわけだ.

繰上がった1は, 分子の5を掛けて, 5として繰上がり, その桁の20と足して25になる. 次はその25を分母の9で除し, 商2, 剰余7になる. 商2に分子の4を掛け, 8を20に足し, 次の桁では28を正規化する.

これを繰り返すと, 第1桁では32を正規化し, 商10に分子1を掛けた10が繰上がり, 20と足して30になる. つまり円周率の10倍の整数部は30ということが分った.

各桁の値は剰余で置き換え, それを10倍すると, 300から始まる行のようになる.

これは最初の取った項数が少いし, 10倍と正規化も1回しかしていないので, 様子が分り難いが, 100項くらい取って, この操作を何回も繰り返したのが, 最後の数列である. この先頭がπを何回も10倍したもので, 3141592653までは正しいπの値が得られている.

今回の説明はこのくらいで.