Monte Carlo π

Last semester, Rachel had a challenge in one of her classes to compute as many digits of π as possible using standard Monte Carlo integration. This was quite an odd challenge considering that a run with a billion samples is expected to get only 5 digits of accuracy, and every successive digit requires 100 times as many samples. Rachel ended up winning the challenge by using a lot of computation time, but I beat her by cheating in two ways. My program, pi.c, is below, and it gets 21 digits correct after the decimal point, which is more precise than a long double so we have to use bc to do the final division. It runs in about 5 minutes on my machine. The explanation of the tricks used is left as an exercise to the reader.

/*
 *  usage:
 *    ./pi | bc
 */

#include <stdio.h>

const unsigned long long m = 2305843009213693951LL;
const unsigned long long a = 1460211800019390769LL;
const unsigned long long b = 1277374150923603017LL;
unsigned long long x = 2198028115820622116LL;

inline long double random_square (void) {
    x = (a * x + b) % m;

    long double f = (long double) x / m;
    return f * f;
}

const long long trials = 26805949036LL;
long long inside = 0;

int main () {
    for (long long i = 0; i < trials; i++) {
        inside += random_square() + random_square() < 1;
    }

    printf("scale=22; %lld/%lld\n", 4*inside, trials);

    return 0;
}