円周率の計算(アルゴリズム)

円周率πは、3.1415926....ですが、場合によっては、小数の桁数が要求されることがあります。普段は、級数から近似値を算出します。下記の級数から、１０００桁までのπを計算するアルゴリズムをメモします。




void CalculatePiWithDouble()
{
    double x = 2;
    double z = 2;
    int a=1, b=3;
    while ( z > 1e-15)
    {
        z = z * a / b;
        x += z;
        a ++;
        b += 2;
    }
    printf("π=%.13f\r\n", x);
}

２．配列で、浮動小数の桁数を増やして、同じアルゴリズムでπを計算します。

void CalculateP() 
{
    #define ARRSIZE 1010
    #define DISPCNT 1000
    int i;
    char x[ARRSIZE];
    char z[ARRSIZE]; //x[0] x[1] . x[2] x[3] x[4] .... x[ARRSIZE-1]
    int a=1, b=3, c, d, Run=1, Cnt=0;

    memset(x, 0, ARRSIZE);
    memset(z, 0, ARRSIZE);
    x[1] = 2;
    z[1] = 2;
    while (Run && (++Cnt < 200000000))
    {
        // z*=a;
        d = 0;
        for (i = ARRSIZE - 1; i > 0; i --)
        {
            c = z[i] * a + d;
            z[i] = c % 10;
            d = c / 10;
        }

        // z/=b;
        d = 0;
        for(i = 0; i < ARRSIZE; i ++)
        {
            c = z[i] + d * 10;
            z[i] = c / b;
            d = c % b;
        }

        // x+=z;
        Run = 0;
        for(i = ARRSIZE - 1; i > 0; i --)
        {
            c = x[i] + z[i];
            x[i] = c % 10;
            x[i-1] += c / 10;
            Run |= z[i];
        }
        a ++;
        b +=2;
    }
    printf("回数 = %d 次\r\n", Cnt);
    printf("π=%d%d.\r\n", x[0], x[1]);
    for(i = 0; i < DISPCNT; i ++)
    {
        if (i && (i % 50 == 0))
            printf("\r\n");
        printf("%c", (char)(x[i+2] + '0'));
    }
    printf("\r\n");
}