Pi - Computing Pi in C

Computing Pi in C

Dik T. Winter wrote a 160-byte C program to compute the first 800 digits of pi. We analyze his code here. The original code

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

can be rewritten as:

#include <stdio.h>

int main() {
    int r[2800 + 1];
    int i, k;
    int b, d;
    int c = 0;

    for (i = 0; i < 2800; i++) {
        r[i] = 2000;
    }

    for (k = 2800; k > 0; k -= 14) {
        d = 0;

        i = k;
        for (;;) {
            d += r[i] * 10000;
            b = 2 * i - 1;

            r[i] = d % b;
            d /= b;
            i--;
            if (i == 0) break;
            d *= i;
        }
        printf("%.4d", c + d / 10000);
        c = d % 10000;
    }

    return 0;
}

Then we see that during the first iteration of the $k$ loop, the code is essentially computing the digits of

2 (1 + \frac{1}{3} (1 + \frac{2}{5} (1 + . . . (1 + \frac{2799}{2 \cdot 2799 + 1} (1 + 0)) . . .)))

By Beeler et al. 1972, Item 120, this is an approximation of $π$ .

We note each term in the approximation gives an additional bit of precision (see above link) thus 14 terms give 4 decimal digits of precision each time (since $2^{14} > 10^{4}$ ). This is why $2800 = 14 \times 200$ terms are used.

To get a better idea of what the program does, let us write

\frac{a_{0}}{b_{0}} = 1, \frac{a_{1}}{b_{1}} = \frac{1}{3}, \frac{a_{2}}{b_{2}} = \frac{2}{5}, . . .

where each $a_{i}, b_{i}$ are coprime positive integers.

Set $n = 2799$ . Then our goal is to compute the digits of

P = 2 (1 + \frac{a_{1}}{b_{1}} (1 + \frac{a_{2}}{b_{2}} (1 + . . . (1 + \frac{a_{n}}{b_{n}} (1)) . . .)

Let $x = 2 \cdot 10^{7}$ . In the first iteration, we compute

P_{0} = x + a_{1} ⌊ \frac{1}{b_{1}} (x + . . . a_{n - 1} ⌊ \frac{1}{b_{n - 1}} (x + a_{n} ⌊ \frac{1}{b_{n}} x ⌋) ⌋ . . .) ⌋

Let $x = q_{n} b_{n} + r_{n}$ and $q_{i} b_{i} + r_{i} = x + q_{i + 1} b_{i + 1}$ for $i < n$ where the $q_{i}, r_{i}$ are positive integers. Then $P_{0} = q_{0}$ , and we have

10^{7} \cdot P = q_{0} + r_{0} + \frac{a_{1}}{b_{1}} (r_{1} + \frac{a_{2}}{b_{2}} (r_{2} + . . . (r_{n - 1} + \frac{a_{n}}{b_{n}} (r_{n})) . . .) .

The program then prints out $q_{0} \mod 10^{4}$ , which should be the first four digits of an eight-digit number.

The first four digits are subtracted away. Then on the next iteration, we roughly compute $10^{4}$ times the error term

r_{0} + \frac{a_{1}}{b_{1}} (r_{1} + \frac{a_{2}}{b_{2}} (r_{2} + . . . (r_{n - 1} + \frac{a_{n}}{b_{n}} (r_{n})) . . .)

add any possible carry to the unprinted digits from the last iteration, and print the first four digits. This process is repeated until 200 digits are being printed. At each step, we can forget about 14 terms because from the above bounds we know they cannot possibly affect any other digits.

However there is one thing I don't understand. Why is it guaranteed only four digits will be printed each time? For example if we add a carry to 9999 then it will take five digits. In fact, 9999 appears in the first 200 digits, but fortunately there is no carry immediately after that.