円周率.jp > プログラム > Spigot アルゴリズム

Spigot アルゴリズム

Spigot アルゴリズムは上の桁から順々に出力する方法で, 上手く利用すれば短いプログラムを作れたり, 手計算で大きな苦労無く π を計算できたりする. (ただし,正確に出力する桁数 n に比例する容量と, n2 に比例する計算時間が必要.)

158Byte で 2400 桁[JB02]

まずは[JB02]に掲載されていたプログラムを紹介する. わずか 158Byte で 2400 桁(小数点以下 2399 桁分)の π を求めることができる.

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 を代入し, 以下のように変形したものを利用している.

π = 2 \Bigl(1+\frac{1}{3} \Bigl(1+\frac{2}{5} \Bigl(1+\frac{3}{7} \Bigl(1+… \Bigl(1+\frac{k}{2k+1} \Bigl(1+… \Bigr) \Bigr) \Bigr)\Bigr) \Bigr)
= 2+\frac{1}{3} \Bigl(2+\frac{2}{5} \Bigl(2+\frac{3}{7} \Bigl(2+… \Bigl(2+\frac{k}{2k+1} \Bigl(2+… \Bigr) \Bigr) \Bigr)\Bigr) \Bigr)

c のループが進むごとに 1000 倍し, 整数部を出力する,という作業を繰り返している. 1000 (==a) 倍しているのは f[b] * a の部分である. その中にある b のループでは, 各項の繰り上がりを計算している. 式とプログラムとの対応は

\frac{{\tt 1}}{{\tt g}} \Bigl( f[b]bd

という感じになっている. プログラムでは分かりづらいが,dg+2 == 2*b+1 で割った後の商なので \frac{1}{2k+1} に当たる部分は無い.

計算項数

1 項進めるごとに足される項が \frac{k}{2k+1} 倍となるので, d 桁の計算に必要な計算項数 n

\prod_{k=1}%5en \frac{k}{2k+1} < 2-n < 10-d
n > (log210)d > 3.3219d

だが,上記の実装では 1 回ループを進むごとに 4 桁ずつ表示するため d=4D である必要がある.このとき

n > log2104D > 4×3.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;
}