Pi by Monte-Carlo

Question

I was inspired by this SO post to investigate a good Java8 way to calculate Pi based on simulation.

I used a similar task to learn about parallel programming on both CUDA, and Intel Xeon Phi processors. Those systems are more geared for parallel programming, but I felt inclined to apply it to 'regular' Java anyway.

Wikipedia has an small section showing this, and includes the following graphic:

The following code performs the above simulation:

 import java.util.Arrays;
 import java.util.Random;

 /**
  * Approximate the value of Pi by using a Monte-Carlo simulation for the area of a circle of radius 1.
  *
  * @author rolf
  *
  */
 public class PiGuess {


     private static final ThreadLocal<Random> LOCAL_RANDOM = new ThreadLocal<Random>() {
         protected Random initialValue() {
             return new Random();
         };
     };

     private static final int CPU_COUNT = Runtime.getRuntime().availableProcessors();

     /**
      * Split a number of samples as evenly as possible over the number of available processors.
      * @param samples the samples to split
      * @return an array containing the number of samples to process on each processor.
      */
     private static final long[] apportion(final long samples) {
         int core = CPU_COUNT;
         final long[] portions = new long[core];
         long remaining = samples;

         while (core > 0) {
             final long part = (remaining - 1 + core) / core;
             core--;
             portions[core] = part;
             remaining -= part;
         }
         return portions;
     }

     /**
      * Calculate the approximate area of a circle (radius 1.0) based on a sample system on a single quadrant of the circle.
      * A parallel mechanism is used to improve performance.
      * @param samples the number of samples to take
      * @return the area of the circle.
      */
     public static final double sampleCircleArea(final long samples) {

         /*
         Monte-Carlo simulation for the area of a circle.

         A circle of radius 1 just fits in a square of sides 2.
         In one quadrant of the square (area 1 by 1) we have a quarter circle
         If we put the center of the circle at origin 0,0, and then randomly sample points
         in that quadrant, we can tell whether that point is in the circle if the ray from the
         origin is shorter than the radius of the circle.
         If the point is at (x,y), then the ray is the 'hypotenuse' (using Pythagoras).

         We know the area of the quadrant, we can sample millions of points in the quadrant,
         and we can calculate a ratio of the quadrant's area that is inside the circle.

         This sampled area, multiplied by 4, gives the area of the circle.

         */

         // how many samples to process in each thread.
         long[] counts = apportion(samples);

         // add up how many samples appear in the circle
         long inside = Arrays.stream(counts).parallel().map(s -> samplePortion(s)).sum();

         // convert the quadrant area back to the circle area.
         return (4.0 * inside) / samples;
     }

     /**
      * Internal sampling method that counts the number of input samples that are inside the circle too.
      * @param samples the samples to calculate
      * @return the count of samples that are inside the circle.
      */
     private static final long samplePortion(final long samples) {
         final Random rand = LOCAL_RANDOM.get();
         long inside = 0;
         for (int i = 0; i < samples; i++) {
             if (isInside(rand)) {
                 inside++;
             }
         }
         return inside;
     }

     /**
      * The core test for each sample, does a random point in the quadrant lie inside the circle.
      * @param rand the source for the random circle.
      * @return true if the random sample is inside the circle.
      */
     private static final boolean isInside(final Random rand) {
         final double x = rand.nextDouble();
         final double y = rand.nextDouble();
         return x * x + y * y <= 1.0;
     }


     public static void main(String[] args) {
         double[] calculations = new double[100];
         for (int i = 0; i < calculations.length; i++) {
             double calc = sampleCircleArea(100000);
             calculations[i] = calc;
             System.out.printf("Loop %d guesses Pi at %.5f%n", i, calc);
         }
         System.out.printf("Overall calculation is %.5f%n", Arrays.stream(calculations).average().getAsDouble());
     }

 }

When I run it with 1,000,000 samples, I end with the following output:

Loop 97 guesses Pi at 3.14128
Loop 98 guesses Pi at 3.14102
Loop 99 guesses Pi at 3.14142
Overall calculation is 3.14154

My questions of particular interest are:

1. Can performance be improved
2. Are Java 8 mechanisms used appropriately? Are there any simplifications I missed?
3. Are there any other observations?

Answer 1

The code and documentation look fine in general, so this review will focus on minor optimizations on a per-method basis.

Use ThreadLocal.withInitial

You can set your LOCAL_RANDOM with the following:

private static final ThreadLocal<Random> LOCAL_RANDOM = ThreadLocal.withInitial(Random::new);

You need to give a Supplier<Random> as argument, and that is exactly what Random::new is.

Reduce logic involved in apportion

In my opinion too much is going on is this method, there needs to be a way to write it in a cleaner way. An example is the following, though it still does not eliminate all increment and decrement operations:

/**
 * Split a number of samples as evenly as possible over the number of available processors.
 * 
 * This method will for example split 77 over 4 as [20, 19, 19, 19]
 * 
 * @param samples the samples to split
 * @return an array containing the number of samples to process on each processor.
 */
private static final long[] apportion(final long samples) {
    int cores = CPU_COUNT;
    long[] portions = new long[cores];

    //calculate minimum amount of samples per core
    long minSamples = samples / cores;
    Arrays.setAll(portions, index -> minSamples);

    //calculate remaining samples after evenly dividing the minimum amount over all cores
    long remaining = samples - (minSamples * cores);

    //split the remaining samples over the cores        
    for (int coreId = 0; remaining > 0; coreId++) {
        portions[coreId]++;
        remaining--;
    }

    return portions;
}

Streamify samplePortion

This suggestion can be a matter of taste, but I believe a piece of code written in a functional programming variant is easier to understand than some calculations involving loops and more importantly it is less error-prone as you cannot possibly implement the count() method incorrectly!

/**
 * Internal sampling method that counts the number of input samples that are inside the circle too.
 * @param samples the samples to calculate
 * @return the count of samples that are inside the circle.
 */
private static final long samplePortion(final long samples) {
    final Random random = LOCAL_RANDOM.get();
    return LongStream.range(0, samples)
        .filter(sample -> isInside(random))
        .count();
}

This solution also reduces the lines of code used, hence it should in theory be easier to maintain.

Use underscores in numeric literals

For clarify you might write the following in the main method:

double calc = sampleCircleArea(100_000);

Note the 100_000 over 100000.

Answer 2

Overall everything looks very appropriate for the given task and I agree with skiwi's points.

If your class PiGuess is not getting instantiated make it final and add a private constructor.

public final class PiGuess {
    private PiGuess() {}
}

Which brings me to the next point: Your code is not object-oriented. Which it doesn't have to be if your goal was to implement a prototypical application.

If you also wanted to exercise OOP, forget about my first comment and consider the following: Remove the static keyword from the methods, add a constructor to PiGuess which takes the sample size and the number of calculations as parameter and add a public method that basicly does what your main method does now (no computing in the constructor, that's not what they are for).

Apart from that, double[] calculations = new double[100]; could be final.