# 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;
}
```

22 August 2023