Calculating PI in C# using Monte-Carlo simulation

The following code sample shows numeric compuation of the number PI using Monte-Carlo simulation. First, a sequential approach is used. Then the Parallel.For construct in TPL is used. In the end, we use Tasks in TPL.


To compute PI we use the same approach. We consider the unit circle inscribed in a square around origo with corners at coordinates (-1,-1), (-1, 1), (1,1) and (1,-1). The number PI can be defined as generating random numbers and looking at the ratio of the numbers inside the circle M divided upon the total numbers generated N. We know that the following can then be expected:

(a) M / N = PI / 4

Why? Because the square has got an area equal to four, remember that the unit square got sides equal to the number 2 and its area is therefore 2 * 2 = 4. The unit circle got a radius of 1, hence its area is PI * 1^2 = PI. The ratio to be expected between the areas of the unit circle and unit rectangle therefore gives the formula above. We can further compute the approximated numeric value of PI equal to:

(b) PI = 4 * (M / N)

This expression (b) is directly from the previous expression (a)
Let's move on the code sample, review the code. I have included a screen shot at the end. The conclusion I got after testing showed after several runs shows that the sequential version runs in about 3.5 seconds on my eight core system with about half the time, about 1.8 seconds using Parallel.For - The last version using Tasks and Tasks.WaitAll give about 1.7 seconds and the quickest compuation, about twice as fast. The iterations I used in the demo was 80 million.
Here is the code written in C#:

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Threading;
using System.Threading.Tasks;

namespace MonteCarloPiApproximation
{
    class Program
    {
        static int numberOfCores = Environment.ProcessorCount; 
        static int iterations = 10000000 * numberOfCores;

        static void Main(string[] args)
        {
            Console.WriteLine("Monte Carlo numeric simulation of PI");
            Console.WriteLine("Iteration limit: " + iterations);
            Console.WriteLine("Number of processor cores on system: " + Environment.ProcessorCount);
            var sw = new Stopwatch();
            sw.Start();

            Console.WriteLine("\nMONTE CARLO SIMULATION");

            MonteCarloPiApproximationSerialSimulation();
            sw.Stop();

            Console.WriteLine("Serial simulation: (ms)" + sw.ElapsedMilliseconds);
            Console.WriteLine();

            sw.Restart();

            Console.WriteLine("\nMONTE CARLO SIMULATION");
            MonteCarloPiApproximationParallellForSimulation(); 

            sw.Stop();

            Console.WriteLine("Parallell simulation using Parallel.For: (ms)" + sw.ElapsedMilliseconds);
            Console.WriteLine();

            sw.Restart();

            Console.WriteLine("\nMONTE CARLO SIMULATION");

            MonteCarloPiApproximationParallelTasksSimulation(); 

            Console.WriteLine("Parallell simulation using parallell Tasks: (ms)" + sw.ElapsedMilliseconds);
            Console.WriteLine();

            sw.Stop();           

            Console.WriteLine("Press Enter Key");



            Console.ReadKey(); 
        }

        private static void MonteCarloPiApproximationParallelTasksSimulation()
        {
            double piApproximation = 0;
            int inCircle = 0;
            double x, y = 0;

            int[] localCounters = new int[numberOfCores];
            Task[] tasks = new Task[numberOfCores];

            for (int i = 0; i < numberOfCores; i++)
            {
                int procIndex = i; //closure capture 
                tasks[procIndex] = Task.Factory.StartNew(() =>
                {
                    int localCounterInside = 0;

                    Random rnd = new Random();

                    for (int j = 0; j < iterations / numberOfCores; j++)
                    {
                        x = rnd.NextDouble();
                        y = rnd.NextDouble();
                        if (Math.Sqrt(x * x + y * y) <= 1.0)
                            localCounterInside++;
                    } 
                    localCounters[procIndex] = localCounterInside;

                });               
            }

            Task.WaitAll(tasks);
            inCircle = localCounters.Sum(); 

            piApproximation = 4 * ((double)inCircle / (double)iterations);

            Console.WriteLine();
            Console.WriteLine("Approximated Pi = {0}", piApproximation.ToString("F8"));
           
        }      

        private static void MonteCarloPiApproximationParallellForSimulation()
        {
            double piApproximation = 0;
            int inCircle = 0;
            double x, y = 0;
                   
            Parallel.For(0, numberOfCores, new ParallelOptions{ MaxDegreeOfParallelism = numberOfCores }, i =>
            {
              
                int localCounterInside = 0;

                Random rnd = new Random(); 

                for (int j = 0; j < iterations / numberOfCores; j++)
                {
                    x = rnd.NextDouble();
                    y = rnd.NextDouble();
                    if (Math.Sqrt(x*x+y*y) <= 1.0)
                        localCounterInside++;                                                        
                }

                Interlocked.Add(ref inCircle, localCounterInside); 
                            
            }); 

            piApproximation = 4 * ((double)inCircle / (double)iterations);

            Console.WriteLine();
            Console.WriteLine("Approximated Pi = {0}", piApproximation.ToString("F8"));
            
        }

        private static void MonteCarloPiApproximationSerialSimulation()
        {
            double piApproximation = 0;
            int total = 0;
            int inCircle = 0; 
            double x,y = 0;
            Random rnd = new Random(); 

            while (total < iterations)
            {
                x = rnd.NextDouble(); 
                y = rnd.NextDouble();

                if ((Math.Sqrt(x*x+y*y) <= 1.0))
                    inCircle++;

                total++;                
                piApproximation =  4 * ((double)inCircle / (double)total); 
            } //while 


            Console.WriteLine();
            Console.WriteLine("Approximated Pi = {0}", piApproximation.ToString("F8"));

        }




    }
}