pi  円周率πの値を多桁数で算出する

int pi(int *data, int keta);

data  （出力）円周率の値が10進4桁ずつ入った整数の配列
keta  （入力）小数点以下の桁数の指定

正常に計算できたときに0、エラーのとき-1。

#define LOG10_5    0.699
#define LOG10_239  2.378

int pi(int *data, int keta)
{
    int i;

    for (i = 0; i < (keta-1)/4+2; i++) data[i] = 0;
    if (series(data, keta, (int)(keta/(2*LOG10_5)+2), 4, 5, 1) != 0)
        return -1;
    if (series(data, keta, (int)(keta/(2*LOG10_239)+2), 1, 239, 0) != 0)
        return -1;
    return 0;
}

int series(int *data, int keta, int maxItem,
           unsigned bunsi, unsigned bunpo, unsigned plusMinus)
{
    unsigned bunpo2;
    unsigned *item, *itemDiv;
    unsigned *imp,  *ikp;
    int      *dap;
    int      len;
    int      i, j, k;
    long     x;

    len = (keta-1)/4+2 +1;
    if ((item = (unsigned *)calloc(len, sizeof(unsigned))) == NULL)
        return -1;
    for (imp = item, i = len; i--; ) *imp++ = 0;

    if ((itemDiv = (unsigned *)calloc(len, sizeof(unsigned))) == NULL) {
        free(item);
        return -1;
    }

    *item = bunsi * bunpo * 4;
    bunpo2 = bunpo * bunpo;

    for (j = 1; maxItem--; j += 2, plusMinus++) {
        x = 0;
        for (i = len, imp = item; i--; ) {
            x = x * 10000 + *imp;
            *imp++ = x / bunpo2;
            x %= bunpo2;
        }
        
        imp = item;
        ikp = itemDiv;
        if (j == 1) for (i = len; i--; ) *ikp++ = *imp++;
        else {
            x = 0;
            for (i = len; i--; ) {
                x = x * 10000 + *imp++;
                *ikp++ = x / j;
                x %= j;
            }
        }

        k = 0;
        dap = data + len;
        ikp = itemDiv + len;
        if (plusMinus & 1) {
            for (i = len; i--; ) {
                *--dap += *--ikp + k;
                if (*dap < 10000) k = 0;
                else {
                    *dap -= 10000;
                    k = 1;
                }
            }
        } else {
            for (i = len; i--; ) {
                *--dap -= *--ikp + k;
                if (*dap >= 0) k = 0;
                else {
                    *dap += 10000;
                    k = 1;
                }
            }
        }
    }
    free(itemDiv);
    free(item);
    return 0;
}

ここで使う式は、最もよく知られている、J. Machinの式
　　　π/4 = 4 arctan (1/5) - arctan (1/239)
である。彼は1706年にこの公式でπを小数点以下100桁を求めた。 同じ公式で、1949年にはENIACにより2037桁が約70時間で、1958 年にはIBM704により1万桁が約100分で求められている。

皆さんの使うコンピュータでどれぐらいの時間 ?

ちなみに、わたしの使っているPentium 120のDOS/Vマシンでは 1万桁は1分弱で計算できた。ハードウェアの進歩ってすごい !

また、パソコンでの計算記録については、若松登志樹氏が1990年 2月に、富士通FM-TOWNS（メモリ2MB）を用いて、184時間45分を かけ、πを100万桁まで計算した。

arctanをによる古典的な公式にほかに以下のものがある。
ガウスの公式
　　　π/4 = 12 arctan (1/18) + 8 arctan (1/57) - 5 arctan (1/239)
クリンジェンシェルナの公式
　　　π/4 = 8 arctan (1/10) - arctan (1/239) - 4 arctan (1/515)
シュテルマーの公式
　　　π/4 = 6 arctan (1/8) + 2 arctan (1/57) + arctan (1/239)
高野喜久雄の公式
　　　π/4 = 12 arctan (1/49) + 32 arctan (1/57) - 5 arctan (1/239) + 12 arctan (1/110443)
柴田昭彦の公式
　　　π/4 = 17 arctan (1/22) + 3 arctan (1/172) - 2 arctan (1/682) - 7 arctan (1/5357)

シュテルマーの公式に8の割り算を工夫すれば高速化期待できるし、 高野の公式を利用して並列処理をすれば良い成績が期待できよう。

arctanを使った公式でπを計算すると、桁数を2倍にするのに計算時間 は4倍の増加になるが、以下の方法では2.1〜2.2倍で済む。

ガウス・ルジャンドルの公式
　　　A = 1、B = (1/2)^1/2、T = 1/4、X = 1
として、AとBの差が必要とする精度より大きい間、次の計算を繰り返し 実行する。
　　　Y = A、A = (A + B)/2、B = (B * Y)^1/2、
　　　T = T - X * (Y-A)²、X = 2 * X
すると、πの値は(A + B)²/(4 * T) となる。ただし、A、B、T、 Yの各値は、求めようとする精度以上で計算しておく必要がある。