## Sunday, 19 July 2015

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

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");

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

sw.Stop();

Console.WriteLine("Press Enter Key");

}

{
double piApproximation = 0;
int inCircle = 0;
double x, y = 0;

int[] localCounters = new int[numberOfCores];

for (int i = 0; i < numberOfCores; i++)
{
int procIndex = i; //closure capture
{
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;

});
}

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++;
}

});

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"));

}

}
}



1. Programming is combination of intelligent and creative work. Programmers can do anything with code. The entire Programming tutorials that you mention here on this blog are awesome. Beginners Heap also provides latest tutorials of Programming from beginning to advance level.
Be with us to learn programming in new and creative way.

1. Hi, Great.. Tutorial is just awesome..It is really helpful for a newbie like me.. I am a regular follower of your blog. Really very informative post you shared here. Kindly keep blogging. If anyone wants to become a .Net developer learn from Dot Net Online Training from India. or learn thru ASP.NET Essential Training Online . Nowadays Dot Net has tons of job opportunities on various vertical industry.
JavaScript Online Training from India

2. Hi, i have another faster method implemented:

{

int[][] l_oResult = new int[numberOfCores][];
int chunk = (p_nTotal - 1) / numberOfCores;
for (int i = 0; i < numberOfCores; i++)
{
l_oResult[i] = new int[2] { i * chunk, (i + 1) * chunk };
}
l_oResult[numberOfCores - 1][1] += ((p_nTotal - 1) % numberOfCores);
return l_oResult;
}
static int SInglePIApprox( int min, int max)
{
Random rnd = new Random();
int inCircle = 0;
double x, y = 0;
for (int i = min; i < max; i++)
{
x = rnd.NextDouble();
y = rnd.NextDouble();

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

private static void MyParallelMonteCarloPiApproximationSimulation()
{

// Parallelized
int[] summ = new int[numberOfCores];
Parallel.For(
0,
numberOfCores,
delegate(int i)
{
summ[i] = SInglePIApprox(StartEnd[i][0], StartEnd[i][1]);
});

int inCircle = summ.Sum();

double piApproximation = 4 * ((double)inCircle / (double)iterations);
Console.WriteLine();
Console.WriteLine("Approximated Pi = {0}", piApproximation.ToString("F8"));

}

3. If you need your ex-girlfriend or ex-boyfriend to come crawling back to you on their knees (no matter why you broke up) you need to watch this video
right away...

(VIDEO) Why your ex will NEVER get back...