Spigot アルゴリズム
Spigot アルゴリズムは上の桁から順々に出力する方法で， 上手く利用すれば短いプログラムを作れたり， 手計算で大きな苦労無く π を計算できたりする． （ただし，正確に出力する桁数 $n$ に比例する容量と， $O(n^2)$ の計算時間が必要．）
158Byte で 2400 桁[JB02]
まずは[JB02]に掲載されていたプログラムを紹介する． わずか 158Byte で 2400 桁（小数点以下 2399 桁分）の $\pi$ を求めることができる．
int a=10000,b,c=8400,d,e,f[8401],g;main(){for(;b-c;)
f[b++]=a/5;for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),
e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}
しかしながらこの形では分かりにくいので見やすい（と思われる）形へ変形する． さらに，意味を取りやすいよう変数名を変更したものを右に書いた．
#include <stdio.h>
int a = 10000;
int c = 8400;
int b;
int d;
int e;
int g;
int f[8401];
int main(void) {
for (b = 0; b < c; b++) {
f[b] = a / 5;
}
e = 0;
for (c = 8400; c > 0; c -= 14) {
d = 0;
for (b = c - 1; b > 0; b--) {
g = 2 * b - 1;
d = d * b + f[b] * a;
f[b] = d % g;
d /= g;
}
printf("%04d", e + d / a);
e = d % a;
}
return 0;
}
#include <stdio.h>
int base = 10000; // 基底
int n = 8400; // 計算項数
int i; // ループ変数
int temp; // 一時変数/繰り上がり
int out; // 出力値
int denom; // 分母
int numerator[8401]; // 分子
int main(void) {
for (i = 0; i < n; i++) {
numerator[i] = base / 5;
}
out = 0;
for (n = 8400; n > 0; n -= 14) {
temp = 0;
for (i = n - 1; i > 0; i--) {
denom = 2 * i - 1;
temp = temp * i + numerator[i] * base;
numerator[i] = temp % denom;
temp /= denom;
}
printf("%04d", out + temp / base);
out = temp % base;
}
return 0;
}
解説
計算アルゴリズム
このプログラムでは Euler の arctan 展開 に $x=1$ を代入し，以下のように変形したものを利用している．
\[ \pi = 2\left(1+\frac13\left(1+\frac25\left(1+\frac37\left(1+\cdots\left(1+\frac{k}{2k+1}\left(1+\cdots\right) \cdots \right.\right.\right.\right.\right) \] \[ = 2+\frac13\left(2+\frac25\left(2+\frac37\left(2+\cdots\left(2+\frac{k}{2k+1}\left(2+\cdots\right) \cdots \right.\right.\right.\right) \]
c のループが進むごとに 1000 倍し，
整数部を出力する，という作業を繰り返している．
1000 (==a) 倍しているのは
f[b] * a の部分である．
その中にある b のループでは，
各項の繰り上がりを計算している．
式とプログラムとの対応は
という感じになっている．
プログラムでは分かりづらいが，d は
g+2 == 2*b+1 で割った後の商なので
$\dfrac{1}{2k+1}$ に当たる部分は無い．
計算項数
1 項進めるごとに足される項が $\dfrac{k}{2k+1}$ 倍となるので， $d$ 桁の計算に必要な計算項数 $n$ は
\[ \prod_{k=1}^n \frac{k}{2k+1} \lt 2^{-n} \lt 10^{-d} \] \[ n \gt (\log_210)d \gt 3.3219 d \]だが，上記の実装では 1 回ループを進むごとに 4 桁ずつ表示するため $d=4D$ である必要がある．このとき
\[ n \gt \log_210^{4D} \gt 4\times3.3219D = 13.2876D \]から（本来は上から抑える証明が必要だが） 1 回計算が進むごとに 13.316 より少し大きな 14 項を削って計算している． 無論全項計算しても良いのだが，$14D$ 番目以降の項は計算しても結果に違いを与えないため 0 とみなして計算している．
133 Byte で 15000 桁[FB02]
こちらも基本的には上記と同じ方式である．
グローバル変数での宣言では 0 が代入されること，
0 の値は条件分岐で「偽」と判断されることを利用することで
上記プログラムの 2 (==a/5)
を代入するループ部分を計算用ループに混合させ，
プログラムのサイズを削減している．
a[52514],b,c=52514,d,e,f=1e4,g,h;
main(){for(;b=c-=14;h=printf("%04d",e+d/f))
for(e=d%=f;g=--b*2;d/=g)d=d*b+f*(h?a[b]:f/5),a[b]=d%--g;}
#include <stdio.h>
int a[52514];
int b;
int c;
int d = 0;
int e;
int f = 10000;
int g;
int h = 0;
int main() {
for (c = 52500; c > 0; c -= 14) {
d %= f;
e = d;
for (b = c - 1; b > 0; --b) {
g = 2 * b - 1;
d = d * b + f * (h ? a[b] : (f / 5));
a[b] = d % g;
d /= g;
}
h = printf("%04d", e + d / f);
}
return 0;
}
#include <stdio.h>
int nume[52514];
int i;
int n;
int carry = 0;
int digit;
int base = 10000;
int denom;
int first = 0;
int main() {
for (n = 52500; n > 0; n -= 14) {
carry %= base;
digit = carry;
for (i = n - 1; i > 0; --i) {
denom = 2 * i - 1;
carry = carry * i + base * (first ? nume[i] : (base / 5));
nume[i] = carry % denom;
carry /= denom;
}
first = printf("%04d", digit + carry / base);
}
return 0;
}