Pi (π) stands as one of the most fascinating mathematical constants, classified as a transcendental number beginning with 3.14159265358.... The quest to calculate pi with ever-increasing accuracy has captivated mathematicians throughout history.

Intriguingly, mathematicians hypothesize that pi may be a normal number—if proven true, this would mean every possible finite sequence of digits appears somewhere within its infinite decimal expansion. Under this assumption, personal identifiers like names, birth dates, and phone numbers could theoretically be located within pi's digits. This conjecture has inspired numerous enthusiasts to search for meaningful patterns with in the constant (tools like Pi-Search enable such investigations).

Computational Milestones

From ancient geometric approaches to modern computational algorithms, humanity's ability to determine pi's value has grown dramatically. As of March 2024, Solidigm (a storage company in California) announced completion of pi calculations to approximately 105 trillion digits, requiring 75 days of computation.

This article explores efficient algorithms for computing pi and their practical implementation in code.

Quick Access via Libraries

For rapid access to many digits, Python's mpmath library provides a straightforward solution:

from mpmath import mp
mp.dps = 1000
print(mp.pi)

This yields 1000 digits almost instantly, leveraging C implementations internally.

The Wallis Product Formula

Among the simplest rapidly-converging series for pi is the Wallis product. The mathematical formulation follows:

π/2 = 1 + 1/3 + (1/3)×(2/5) + (1/3)×(2/5)×(3/7) + ...

Reorganized for computation:

π = 2 + 2/3 + (2/3)(2/5) + (2/3)(2/5)(3/7) + ...

Standard Double Implementation

#include <stdio.h>

int main() {
    double result = 2.0, term = 2.0;
    int numerator = 1;
    int denominator = 3;
    
    while (term > 0) {
        term *= (double)numerator / denominator;
        result += term;
        numerator++;
        denominator += 2;
    }
    
    printf("π=%.14f\n", result);
    return 0;
}

Output: π=3.1415926535898

This approach yields 14 correct decimal places, limited by IEEE 754 double precision (8 bytes).

High-Precision Calculation with GMP

For extended precision, the GNU Multiple Precision Arithmetic Library provides arbitrary-precision floating-point types:

#include <stdio.h>
#include <gmp.h>

int main() {
    mpf_set_default_prec(256);
    
    mpf_t total, term, num, den, threshold;
    mpf_init(total);
    mpf_init(term);
    mpf_init(num);
    mpf_init(den);
    mpf_init(threshold);
    
    mpf_set_ui(total, 2);
    mpf_set_ui(term, 2);
    mpf_set_ui(num, 1);
    mpf_set_ui(den, 3);
    mpf_set_str(threshold, "0.00000000000000000001", 20);
    
    int iterations = 0;
    while (mpf_cmp(term, threshold) > 0) {
        mpf_mul(term, term, num);
        mpf_div(term, term, den);
        mpf_add(total, total, term);
        mpf_add_ui(num, num, 1);
        mpf_add_ui(den, den, 2);
        gmp_printf("[%d]%.50Ff\n", ++iterations, total);
    }
    
    gmp_printf("[Final]π=%.26Ff\n", total);
    
    mpf_clear(total);
    mpf_clear(term);
    mpf_clear(num);
    mpf_clear(den);
    mpf_clear(threshold);
    
    return 0;
}

Compile with: gcc program.c -lgmp

Arbitrary Precision with Arrays

For complete control over digit count, implementing big integer arithmetic using arrays offers maximum flexibility:

#include <stdio.h>
#include <string.h>

int main() {
    const int BUFFER_SIZE = 1010;
    const int DIGIT_COUNT = 1000;
    
    int result[BUFFER_SIZE];
    int working[BUFFER_SIZE];
    int p = 1, q = 3, carry, temp, position;
    int active = 1, iterations = 0;
    
    memset(result, 0, sizeof(result));
    memset(working, 0, sizeof(working));
    result[1] = 2;
    working[1] = 2;
    
    while (active && ++iterations < 100000) {
        // Multiply: working = working * p
        carry = 0;
        for (position = BUFFER_SIZE - 1; position >= 0; position--) {
            temp = working[position] * p + carry;
            working[position] = temp % 10;
            carry = temp / 10;
        }
        
        // Divide: working = working / q
        carry = 0;
        for (position = 0; position < BUFFER_SIZE; position++) {
            temp = working[position] + carry * 10;
            working[position] = temp / q;
            carry = temp % q;
        }
        
        // Accumulate: result += working
        active = 0;
        for (position = BUFFER_SIZE - 1; position > 0; position--) {
            temp = result[position] + working[position];
            result[position] = temp % 10;
            result[position - 1] += temp / 10;
            active |= working[position];
        }
        
        p++;
        q += 2;
    }
    
    printf("Completed in %d iterations\n", iterations);
    printf("Pi = %d.\n", result[1]);
    
    for (position = 0; position < DIGIT_COUNT; position++) {
        if (position > 0 && position % 100 == 0)
            printf("\n");
        printf("%d", result[position + 2]);
    }
    
    return 0;
}

Execution reveals that after 3,343 iterations, the first 999 decimal places are accurate.

Scaling to 10,000 Digits

Adjusting constants to:

const int BUFFER_SIZE = 10100;
const int DIGIT_COUNT = 10000;

yields 10,000 correct digits after 33,538 iterations. The ratio confirsm approximately 3.3 iterations per accurate digit.

Extending to Other Constants

The same digit-by-digit methodology applies to calculating other mathematical constants. For Euler's number e = 2.71828...:

#include <stdio.h>

int main() {
    int digits = 100;
    int factorial[101];
    int n, remainder, position;
    
    for (position = 0; position <= digits; position++)
        factorial[position] = 10;
    
    printf("e=2.");
    remainder = 0;
    
    for (position = 0; position < digits; position++) {
        for (n = digits; n > 0; n--) {
            remainder = remainder + factorial[n] * 10;
            factorial[n] = remainder % n;
            remainder = remainder / n;
        }
        printf("%d", (remainder / 10) % 10);
        remainder = remainder % 10;
    }
    
    printf("\n\n");
    return 0;
}

Verification via mpmath:

import mpmath
mpmath.mp.dps = 1000
print(mpmath.e)

Historical Context: Ramanujan's Formulas

In 1914, mathematician Srinivasa Ramanujan published 14 formulas for pi in his paper. His famous series achieves 8 correct decimal digits per term. In 1985, Gosper used this formula to compute 17,500,000 digits. The Cuhdnovsky brothers refined Ramanujan's approach in 1989, creating the Chudnovsky formula yielding 15 digits per term, and calculated 4,044,000,000 digits by 1994.

Code Golf: Computing Mersenne Prime Digits

Bellard's IOCCC 2000 entry calculates the Mersenne prime 2^6972593-1 (2,098,960 decimal digits) using Number Theoretic Transform techniques:

#include <stdio.h>

int m = 754974721, N, t[1 << 22], a, *p, i, e = 1 << 22, j, s, b, c, U;
f(d) {
    for (s = 1 << 23; s; s /= 2, d = d * 1LL * d % m)
        if (s < N)
            for (p = t; p < t + N; p += s)
                for (i = s, c = 1; i; i--)
                    b = *p + p[s], p[s] = (m + *p - p[s]) * 1LL * c % m,
                    *p++ = b % m, c = c * 1LL * d % m;
    for (j = 0; i < N - 1;) {
        for (s = N / 2; !((j ^= s) & s); s /= 2);
        if (++i < j)
            a = t[i], t[i] = t[j], t[j] = a;
    }
}

main() {
    *t = 2;
    U = N = 1;
    while (e /= 2) {
        N *= 2;
        U = U * 1LL * (m + 1) / 2 % m;
        f(362);
        for (p = t; p < t + N;)
            *p++ = (*p * 1LL * *p % m) * U % m;
        f(415027540);
        for (a = 0, p = t; p < t + N;)
            a += (6972593 & e ? 2 : 1) * *p,
            *p++ = a % 10, a /= 10;
    }
    while (!*--p);
    t[0]--;
    while (p >= t)
        printf("%d", *p--);
}

Compact Implementation

An elegant IOCCC submission computing 800 digits:

#include <stdio.h>
long 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);
}

Conclusion

From textbook formulas to optimized libraries, computing pi demonstrates the beautiful intersection of mathematical theory and computational practice. The array-based approach, while conceptually simple, remains remarkably effective—proving that sometimes the most fundamental techniques yield the highest precision.