円周率の計算 (1) - シンセ・アンプラグド

今回からは、シンセ関係から離れて、円周率の計算について触れたいと思います。
それと言うのも、20 個ほどの LSTTL を使った「RETROF-8」という 8 ビット・コンピュータを製作されている方がいて、Web サイト上で進行状況を公開されています。

TTLのみでコンピュータを自作する

ハードウェアは、ほぼ完成し、現在、円周率の計算プログラムの作成が進行中です。
8 ビット CPU の乏しいリソースでの計算は、もちろん、計算桁数世界一を目指すような世界とは全く異なりますが、「8 ビット縛り」という点に興味を覚え、私も限界を追求してみたくなりました。
8 ビット CPU のターゲットとして、ATmega シリーズを想定し、結果をシリアルから出力する Arduino 用のスケッチとして実現しました。　バイナリ・スケッチ・サイズは 5324 バイトです。
スケッチは長くなるので、次回、リストを掲載します。　PC 上の gcc でコンパイルしても実行できます。
(8 月 17 日追記: atan255 の公式の 44*atan(1/2072) の符号が誤っていたのを修正しました。)
「8 ビット縛り」の条件として、次のようなものを考えました。

計算に使う「多倍長変数」として、uint8_t つまり符号なし 8 ビット整数の配列を使い、その配列要素数は 8 ビット・インデクスでアクセス可能な 255 を上限とする。
加減算が 8 ビットなのはもちろん、乗算は、被乗数 8 ビット x 乗数 8 ビット = 積 16 ビット、除算は、被除数 16 ビット ÷ 除数 8 ビット = 商 8 ビット、余り 8 ビットの演算に限る。

RETROF-8 では、4 K バイトのアドレス空間の中にデータおよびプログラムを収めなければなりませんが、Arduino 版ではプログラムやデータサイズについては制限を設けません。
また、計算時間についても、減らすための努力は特にしていません。
1. で配列要素数の上限が 256 ではなく、255 なのは、配列要素数の数値自体も 8 ビット幅におさめるためで、数値の整数部に 1 バイト使うので、数値の小数部を表現する要素数は最大 254 になります。
C 言語の仕様では、処理系のデフォルト int サイズより小さなサイズの整数型の変数の値は、int 型のサイズまで「プロモーション」されてから演算されるので、AVR-GCC では、8 ビット整数の演算でも 16 ビット・サイズにプロモーションされて演算ルーチンが呼ばれます。
そのため、C 言語の中だけの記述では、8 ビット演算のみ使用するという 2. での縛りは確実とはいえません。
そこで、乗算および除算の部分はインライン・アセンブラ機能を使って、アセンブリ言語で書いた8 ビット演算ルーチンを使用しています。
また、AVR の (ハードウェア) 乗算命令やシフト/ローテイト命令を使わず、加減算命令だけを使って計算する方法も選択できるようにしました。
計算方式としては、リソースの乏しいコンピュータ用としては定番の Machin (マチン) の公式などの π を表す arctan の公式を使って、arctan 自体の値は級数展開の公式を使って求める方法で実行しています。
また、数値は、精度を稼ぐために、多倍長整数にバイナリの表現で格納しています。
そのため、バイナリの表現で π が求まったあと、10 進数として表示するためにも演算が必要になります。
後で示す実行結果の経過時間には、この表示のための演算に要する時間は含んでいません。
254 バイト、つまり 254 x 8 = 2032 ビットの 2 進小数部で表現できる精度を 10 進数の桁数に換算するには、桁数を d として、
$\qquad\qquad 10^d\, =\, 2^{254\cdot 8}=2^{2032}$
の関係から、両辺の常用対数を取って、
$\qquad\qquad \log_{10}(10^d)\, =\, \log_{10}(2^{2032})$
$\qquad\qquad d\, \cdot\, \log_{10}(10)\, =\, 2032\, \cdot\, \log_{10}(2)$
$\qquad\qquad \begin{eqnarray} d &=& 2032 \,\cdot\, \log_{10}(2) \,=\, 2032\,\cdot\, 0.301 \\ &=& 611.7\end{eqnarray}$
となり、10 進 611 桁程度が精度の限界となります。
加減算命令に限った場合の実行結果を下に示します。

Machin (digits = 182), Elapsed time = 1458 ms
pi = 3.
1415926535 8979323846 2643383279 5028841971 6939937510 
5820974944 5923078164 0628620899 8628034825 3421170679 
8214808651 3282306647 0938446095 5058223172 5359408128 
4811174502 8410270193 8521105559 51

Gauss (digits = 325), Elapsed time = 4985 ms
pi = 3.
1415926535 8979323846 2643383279 5028841971 6939937510 
5820974944 5923078164 0628620899 8628034825 3421170679 
8214808651 3282306647 0938446095 5058223172 5359408128 
4811174502 8410270193 8521105559 6446229489 5493038196 
4428810975 6659334461 2847564823 3786783165 2712019091 
4564856692 3460348610 4543266482 1339360726 0249141273 
7245870066 0631558817 48768

Stoermer (digits = 453), Elapsed time = 10898 ms
pi = 3.
1415926535 8979323846 2643383279 5028841971 6939937510 
5820974944 5923078164 0628620899 8628034825 3421170679 
8214808651 3282306647 0938446095 5058223172 5359408128 
4811174502 8410270193 8521105559 6446229489 5493038196 
4428810975 6659334461 2847564823 3786783165 2712019091 
4564856692 3460348610 4543266482 1339360726 0249141273 
7245870066 0631558817 4881520920 9628292540 9171536436 
7892590360 0113305305 4882046652 1384146951 9415116094 
3305727036 5759591953 0921861173 8193261179 3105118548 
036

atan255 (digits = 610), Elapsed time = 29842 ms
pi = 3.
1415926535 8979323846 2643383279 5028841971 6939937510 
5820974944 5923078164 0628620899 8628034825 3421170679 
8214808651 3282306647 0938446095 5058223172 5359408128 
4811174502 8410270193 8521105559 6446229489 5493038196 
4428810975 6659334461 2847564823 3786783165 2712019091 
4564856692 3460348610 4543266482 1339360726 0249141273 
7245870066 0631558817 4881520920 9628292540 9171536436 
7892590360 0113305305 4882046652 1384146951 9415116094 
3305727036 5759591953 0921861173 8193261179 3105118548 
0744623799 6274956735 1885752724 8912279381 8301194912 
9833673362 4406566430 8602139494 6395224737 1907021798 
6094370277 0539217176 2931767523 8467481846 7669405132 
0005681261

ここで、(digits = ) の数字は、単に表示桁数を表しているだけです。
表示桁数は実際に正しい値となる桁数 + 2 〜 3 桁を選んであり、正しい桁数と、実行時間は下の表にまとめました。

公式	有効桁数	実行時間 (s) (加減算のみ)	実行時間 (s) (シフトあり)
Machin	180	1.5	0.3
Gauss	322	5.0	0.9
Stoermer	451	10.9	1.8
atan255	608	29.8	4.8

加減算のみの場合と、シフト/ローテイト命令を使った場合との実行時間の比率は 5 程度になっています。
表の公式の名称が「atan255」となっているのは、初項が arctan(1/255) となる公式で、名前が付けられているかどうか分からないので、仮の名称を「atan255」としているものです。
atan255 を使った場合には、限界に近い有効桁数 608 桁が得られ、加減算のみでも実行時間は 30 秒ほどです。
最後に、各公式を掲載しておきます。
Machin (マチン) の公式

$\qquad \frac{\,\pi\,}{4}\, =\, 4 \,\cdot\, \arctan\left( \frac{1}{\,5\,} \right) \, -\, \arctan\left( \frac{1}{\,239\,} \right)$

Gauss (ガウス) の公式

$\qquad \frac{\,\pi\,}{4}\, =\, 12 \,\cdot\, \arctan\left( \frac{1}{\,18\,} \right) \, +\, 8\,\cdot\, \arctan\left( \frac{1}{\,57\,} \right) \, -\, 5\,\cdot\, \arctan\left( \frac{1}{\,239\,} \right)$

Störmer (シュテルマー) の公式

$\qquad \begin{eqnarray}\frac{\,\pi\,}{4}& =& 44 \,\cdot\, \arctan\left( \frac{1}{\,57\,} \right) \, +\, 7\,\cdot\, \arctan\left( \frac{1}{\,239\,} \right) \\& -& 12\,\cdot\, \arctan\left( \frac{1}{\,682\,} \right) \, +\, 24\,\cdot\, \arctan\left( \frac{1}{\,12943\,} \right)\end{eqnarray}$

atan255 (仮) の公式

$\qquad \begin{eqnarray}\frac{\,\pi\,}{4}& =& 227 \,\cdot\, \arctan\left( \frac{1}{\,255\,} \right) \, -\, 100\,\cdot\, \arctan\left( \frac{1}{\,682\,} \right) \\& +& 44\,\cdot\, \arctan\left( \frac{1}{\,2072\,} \right) \, +\, 51\,\cdot\, \arctan\left( \frac{1}{\,2943\,} \right) \\& -& 27\,\cdot\, \arctan\left( \frac{1}{\,12943\,} \right) \, +\, 88\,\cdot\, \arctan\left( \frac{1}{\,16432\,} \right)\end{eqnarray}$