Can a Monte Carlo pi calculation be used for a world record?

Question

I have this random function to calculate pi Monte Carlo style:

max=10000000;
format long;

in = 0;
tic
for k=1:max
    x = rand();
    y = rand();
    if sqrt(x^2 + y^2) < 1
        in = in + 1;
    end
end

toc
calc_pi = 4*in/max 
epsilon = abs(pi - calc_pi)

calcPi(100000000);

If I could iterate this 10e100 times, could this algorithm compete with the world record? If so, how can I find the number of iteration that will give the Nth digit?

Answer 1

This is a nice exercise for calculating pi, but it is probably a very inefficient one. Some remarks:

• My statistics are rusty before I had my coffee, but I guess the error scales with 1 / sqrt(n_guess). To get N digits, you need an error of 10^(-N), so you would need roughly (10^N)^2 random guesses. If you would do 1e100 guesses, as you proposed, you would only get on the order of 50 digits of pi! The number of iteration is thus some exponentional function of the number of required digits, which is horribly slow. A good algorithm is maybe linear in the number of digits you want.

• With the large number of guesses required, you have to start questioning the quality of your random number generator.

• Your algorithm will be limited by floating-point errors to 1e-16 or so. Calculating digits of pi requires some sort of arbitrary precision number format.

• To speed up your algorithm, you can leave out the sqrt().

• Don't use a variable called max, this overwrites an existing function. Use n_guess or so.

Quick and dirty test to prove my theory (after coffee):

pie = @(n) 4 * nnz(rand(n,1).^2 + rand(n, 1).^2 < 1) / n;
ntrial = round(logspace(1, 8, 100));
pies = arrayfun(pie, ntrial);

loglog(ntrial, abs(pies - pi), '.k', ntrial, ntrial.^-.5, '--r')
xlabel('ntrials')
ylabel('epsilon')
legend('Monte Carlo', '1 / sqrt(ntrial)')