Estimating Pi with Monte Carlo Simulation

_monteCarlo.md

Estimating Pi with Monte Carlo Simulation

Notes:

• I decided to use Random() for the random number generator, but I seed it with a RNGCryptoServiceProvider. That should be "random enough" for what we're trying to achieve (the first 5-7 digits of pi).
• Program utilizes Parallel.For to run simulations in parallel. It will literally thrash your CPU to 100%.
• Make sure to pick a radius that's a power of 2, otherwise the modulo won't distribute your random values evenly.

monteCarlo.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Security.Cryptography;
using System.Text;
using System.Threading.Tasks;

namespace MonteCarloPi
{
    class Program
    {
        static object lockObject = new object();

        static void Main(string[] args)
        {
            ulong radius = 2147483648; // 2^31
            ulong inside = 0;
            ulong total = 0;
            ulong batchSize = 200000;
            var seeder = new RNGCryptoServiceProvider();
            var bytes = new byte[4];
            var parallelFactor = Environment.ProcessorCount;

            while(true)
            {
                Parallel.For(0, parallelFactor, (s) =>
                {
                    seeder.GetBytes(bytes);
                    var seed = BitConverter.ToInt32(bytes, 0);
                    var res = Calc(radius, batchSize, seed);

                    lock (lockObject) // update the total count
                    {
                        inside += res;
                        total += batchSize;
                    }
                });

                var pi = (4 * inside * 1000000000) / (total);
                var piString = pi.ToString().Insert(1, ".");
                Console.WriteLine("Total: {0} million, Pi: {1}", total / 1000000, piString);
            }
        }

        static ulong Calc(ulong radius, ulong total, int seed)
        {
            var rand = new Random(seed);
            byte[] bytes = new byte[8];
            ulong i = 0;
            ulong inside = 0;

            while (i < total)
            {
                rand.NextBytes(bytes);
                var x = BitConverter.ToUInt64(bytes, 0) % radius;
                rand.NextBytes(bytes);
                var y = BitConverter.ToUInt64(bytes, 0) % radius;

                i++;
                if ((x * x + y * y) <= radius * radius)
                    inside++;
            }

            return inside;
        }
    }
}