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분 미만으로 계산할 수 있었다. 하드웨어의 진보는 대단해 !

또, PC에서의 계산 기록에 대해서는, 와카마츠등지 이츠키씨가 1990년 2월에, 후지쯔 FM-TOWNS(메모리 2 MB)를 이용해, 184시간 45분을 걸쳐π를 100만 자리수까지 계산했다.

arctan를에 의한 고전적인 공식으로 그 밖에 이하의 것이 있다.
Gauss의 공식
　　　π/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배로 끝난다.

Gauss·르잘돌의 공식
　　　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의 각 치는, 요구하려고 하는 정밀도 이상으로 계산해 둘 필요가 있다.