円周率.jp (http://xn--w6q13e505b.jp/program/spigot.html)

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;
}

解説

計算アルゴリズム

このプログラムでは Eular の 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 のループでは, 各項の繰り上がりを計算している. 式とプログラムとの対応は

\[ \cdots \frac{\tt 1}{\tt g} \left({\tt f[b]} + {\tt b}\cdot{\tt d} \right. \cdots \]

という感じになっている. プログラムでは分かりづらいが,dg+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;
}