Spigot アルゴリズム

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

158Byte で 2400 桁[JB02]

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

  1. int a=10000,b,c=8400,d,e,f[8401],g;main(){for(;b-c;)
  2. f[b++]=a/5;for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),
  3. e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}

しかしながらこの形では分かりにくいので見やすい（と思われる）形へ変形する． さらに，意味を取りやすいよう変数名を変更したものを右に書いた． 

  1. #include <stdio.h>
  2.  
  3. int a = 10000;
  4. int c = 8400;
  5. int b;
  6. int d;
  7. int e;
  8. int g;
  9. int f[8401];
  10.  
  11. int main(void) {
  12.   for (b = 0; b < c; b++) {
  13.     f[b] = a / 5;
  14.   }
  15.  
  16.   e = 0;
  17.   for (c = 8400; c > 0; c -= 14) {
  18.     d = 0;
  19.     for (b = c - 1; b > 0; b--) {
  20.       g = 2 * b - 1;
  21.       d = d * b + f[b] * a;
  22.       f[b] = d % g;
  23.       d /= g;
  24.     }
  25.     printf("%04d", e + d / a);
  26.     e = d % a;
  27.   }
  28.   return 0;
  29. }
  1. #include <stdio.h>
  2.  
  3. int base = 10000;     // 基底
  4. int n = 8400;         // 計算項数
  5. int i;                // ループ変数
  6. int temp;             // 一時変数/繰り上がり
  7. int out;              // 出力値
  8. int denom;            // 分母
  9. int numerator[8401];  // 分子
  10.  
  11. int main(void) {
  12.   for (i = 0; i < n; i++) {
  13.     numerator[i] = base / 5;
  14.   }
  15.  
  16.   out = 0;
  17.   for (n = 8400; n > 0; n -= 14) {
  18.     temp = 0;
  19.     for (i = n - 1; i > 0; i--) {
  20.       denom = 2 * i - 1;
  21.       temp = temp * i + numerator[i] * base;
  22.       numerator[i] = temp % denom;
  23.       temp /= denom;
  24.     }
  25.     printf("%04d", out + temp / base);
  26.     out = temp % base;
  27.   }
  28.   return 0;
  29. }

解説

計算アルゴリズム

このプログラムでは Euler の arctan 展開 に x=1 を代入し，以下のように変形したものを利用している．

π=2(1+13(1+25(1+37(1+⋯(1+k2k+1(1+⋯)⋯) =2+13(2+25(2+37(2+⋯(2+k2k+1(2+⋯)⋯)

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

⋯1g(f[b]+b⋅d⋯

という感じになっている． プログラムでは分かりづらいが，d は g+2 == 2*b+1 で割った後の商なので 12k+1 に当たる部分は無い．

計算項数

1 項進めるごとに足される項が k2k+1 倍となるので， d 桁の計算に必要な計算項数 n は

n∏k=1k2k+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) を代入するループ部分を計算用ループに混合させ， プログラムのサイズを削減している． 

  1. a[52514],b,c=52514,d,e,f=1e4,g,h;
  2. main(){for(;b=c-=14;h=printf("%04d",e+d/f))
  3. for(e=d%=f;g=--b*2;d/=g)d=d*b+f*(h?a[b]:f/5),a[b]=d%--g;}
  1. #include <stdio.h>
  2.  
  3. int a[52514];
  4. int b;
  5. int c;
  6. int d = 0;
  7. int e;
  8. int f = 10000;
  9. int g;
  10. int h = 0;
  11.  
  12. int main() {
  13.   for (c = 52500; c > 0; c -= 14) {
  14.     d %= f;
  15.     e = d;
  16.     for (b = c - 1; b > 0; --b) {
  17.       g = 2 * b - 1;
  18.       d = d * b + f * (h ? a[b] : (f / 5));
  19.       a[b] = d % g;
  20.       d /= g;
  21.     }
  22.     h = printf("%04d", e + d / f);
  23.   }
  24.   return 0;
  25. }
  1. #include <stdio.h>
  2.  
  3. int nume[52514];
  4. int i;
  5. int n;
  6. int carry = 0;
  7. int digit;
  8. int base = 10000;
  9. int denom;
  10. int first = 0;
  11.  
  12. int main() {
  13.   for (n = 52500; n > 0; n -= 14) {
  14.     carry %= base;
  15.     digit = carry;
  16.     for (i = n - 1; i > 0; --i) {
  17.       denom = 2 * i - 1;
  18.       carry = carry * i + base * (first ? nume[i] : (base / 5));
  19.       nume[i] = carry % denom;
  20.       carry /= denom;
  21.     }
  22.     first = printf("%04d", digit + carry / base);
  23.   }
  24.   return 0;
  25. }