円周率計算プログラム

昔書いた円周率のプログラムに少し手を入れたもの。

言語は C です。
上の桁から計算していって、計算したものから順に表示します。
計算する様子をひたすら眺めていたい人向けです。
一応、10万桁程度なら数分で計算できます。
速さや桁数を追求したい人は他をあたってください。

#include <stdio.h>
#define M 167000
#define N 10
unsigned a[M], b[M], c[N], i, j, k, m, s, x;
main()
{
    for (i = -N, s = x = 3; j < M - 3; i++, s = 0, m = j) {
        for (j = k = 1; j < m || x; j++) {
            x = 100 * a[j] + x * k++;
            a[j] = x % (4 * k);
            x /= 4 * k++;
            b[j] = 100 * b[j] + x;
            s += b[j] / k;
            b[j] %= k;
        }
        for (k = 0; k < N; k++) {
            s += 100 * c[k];
            c[k] = s % 100;
            s /= 100;
        }
        if (i > M) continue;
        printf("%02d", s);
        if (i % 500 == 0) printf("\n[%d]", 2 * i);
        if (i % 25 == 0) printf("\n");
        else if (i % 5 == 0) printf(" ");
    }
}

・n 桁求めたいなら、プログラム2行目の M を 1.67n 程度にする。
　　現在の値は 167000 だから10万桁計算できる。
　　(漸近的には だが、少し余裕が必要なので 1.67n にしている)
・計算時間は O(n^2) なので、桁数が10倍になれば100倍の時間がかかる。
・M が 2^32/100*(3/32) ≒ 400万 を超えると32ビットの符号なし整数がオーバーフローして無意味な計算になる。
　　だから求められる桁数の上限は 400万/1.67 = 240万桁 程度。
・使ってる公式は arctan 系ではないが、比較的有名なもの。
・最初の最初に 0 が出力されるのは手抜きではなく仕様。

(プログラムは、短くしようとして変な書き方になってるけど、普段はもう少し行儀がいい)

追記 : ideone.com にも置いてみた http://ideone.com/ToWR1j