How can I make my program that calculates π faster?

Question

For my first Raku program, I thought it might be a fun challange to port over a program I'd written in C# about two years ago, and I was right: it was.

There's just one problem: while the C# program runs in 5 seconds, the Raku program takes nearly 5 minutes. I think the Raku program is doing more or less the same as the C# code, and so it should be equally fast (or equally slow, depending on how you look at it).

Here's my code:

use v6;

class Point {
    has Num $.x;
    has Num $.y;

    method is-inside-unit-circle(--> Bool) {
        $!x²+$!y² <= 1;
    }
}

sub random-points(Int $count --> Seq) {
    gather for 0..^$count {
        take Point.new(x => 1.rand.Num, y => 1.rand.Num);
    }
}

sub compute-pi(Int $batch --> Seq) {
    my Num $total = 0e0;
    my Num $count = 0e0;

    gather for 0..^∞ {
        my Seq $pointsInside = random-points($batch).grep: *.is-inside-unit-circle;

        $total += $batch;
        $count += $pointsInside.elems;
        my Num $ratio = $count / $total;

        take $ratio * 4;
    }
}

{ say "π ≈ $_" } for compute-pi(100_000).head(500);

And here's part of the output:

% time raku estimation.raku
π ≈ 3.12912
π ≈ 3.12874
π ≈ 3.1326533333333333
…
π ≈ 3.1419357429718877
π ≈ 3.1419323446893785
π ≈ 3.14193696
raku estimation.raku  272,98s user 1,78s system 100% cpu 4:34,76 total

My two questions are:

1. How can I make it faster?
2. Is my code idiomatic, and how could I make it more idiomatic?

---

For completeness' sake, here's the C# code:

using System;
using System.Collections.Generic;
using System.Linq;

internal class Program
{
    internal readonly struct Point
    {
        private readonly double _x;
        private readonly double _y;

        internal Point(double x, double y)
        {
            _x = x;
            _y = y;
        }

        internal bool IsInsideUnitCircle()
        {
            return Math.Pow(_x, 2) + Math.Pow(_y, 2) <= 1;
        }
    }

    internal class PointGenerator
    {
        private readonly Random _random;

        public PointGenerator() => _random = new Random();

        public IEnumerable<Point> GeneratePoints(int count)
        {
            for (var i = 0; i < count; i++)
            {
                yield return new Point(_random.NextDouble(), _random.NextDouble());
            }
        }
    }
    
    private static IEnumerable<double> ComputePi(int batch)
    {
        var pointGenerator = new PointGenerator();

        var total = 0.0;
        var count = 0.0;

        while (true)
        {
            var pointsInside = pointGenerator.GeneratePoints(batch)
                .Where(point => point.IsInsideUnitCircle());

            total += batch;
            count += pointsInside.Count();
            var ratio = count / total;

            yield return ratio * 4;
        }
    }

    private static void Main()
    {
        foreach (var estimate in ComputePi(100_000).Take(500))
        {
            Console.WriteLine($"π ≈ {estimate}");
        }
    }
}

And part of its output:

% time dotnet run -c Release
π ≈ 3,13948
π ≈ 3,13904
π ≈ 3,138786666666667
…
π ≈ 3,1415373493975904
π ≈ 3,14154501002004
π ≈ 3,14155416
dotnet run -c Release  5,43s user 0,30s system 103% cpu 5,555 total

Answer 1

You random-points sub has a few things which can be removed or simplified.

• 1.rand returns a Num, so there is no reason to call the method .Num on its result.
• Rather than gather for take you can use the xx operator.
• You only want a positive number for the count, so declare that with UInt

sub random-points(UInt $count --> Seq) {
    Point.new(x => 1.rand, y => 1.rand) xx $count
}

xx seems to be on the order of 10× faster than gather for take.
The reason is mostly that take works by throwing a CONTROL and gather captures that message.
On my computer that change alone takes it from nearly 9 minutes (8m49.881s) to around 2 and a half minutes (2m27.582s). (The other changes have no meaningful effect on the run time.)

The sweet spot for the gather take feature is where there isn't already a built-in for what you are doing, and the speed doesn't matter as much as the ease of writing it.

You could go to the work of creating your own Iterator class which could speed it up. That is why xx is so much faster, there is already an Iterator class that was created specifically for it.

(Note that I did actually test to see if writing a new Iterator class would help. It didn't help, in fact it was slower.)

---

There are similar issues in compute-pi.

To start off, there are improvements to be made to gather for 0..^∞.

1. Don't mix Int and Num endpoints for Range objects if you can help it. 0e0..∞
2. Don't use for for infinite loops, just use loop

gather loop {

    ...

    take $ratio * 4;
}

I would not store the filtered random-points Seq in a variable, because the only thing you need from it is the count.

my Int $points-inside = +random-points($batch).grep: *.is-inside-unit-circle;

The prefix + operator just coerces to a Numeric value. Obviously the only number that makes sense for a Seq is the same as .elems, the count of elements.
.elems could potentially return a Num, specifically Inf. Considering that would break the rest of the code anyway, I think it's fine to store it in an Int scalar.

---

Ok now that we are done with the more obvious cleanups we can try to make it faster.

Thankfully changing from gather for take in random-points to xx helped a lot. (No other change had any meaningful change on the amount of time taken.)

The first thing I did at this point was to run the profiler

raku --profile test.raku

There were no obvious slow parts of our code. In fact a good ⅔ of the run time was in the Raku runtime, not our code.

So let's add parallelism to see if we can reduce the time by using more CPU cores.

We could just make the sequence returned from random-points a race sequence to see if it helps.

my Int $pointsInside = +random-points($batch).race.grep(*.is-inside-unit-circle);

After that change it now takes 4 minutes, almost twice the time it took from before that change.
I did mess around with the batch and degree arguments to race, but no combination seemed to improve it meaningfully. At best I made it match the original timing.

So now we know that was a dead-end. random-points.grep was already about as fast as possible.

---

So let's add the asynchrony to compute-pi instead.

One obvious thing we could do is have a bunch of threads creating the sequences of random-points.grep.

Of course we don't know how many we need to create. So let's give that as an argument.

sub compute-pi(UInt $batch-size, UInt $batch-count --> Seq) {

Now we want to have an existing sequence that we can call race on.

my \batches = ^$batch-count;

Add race to it.
(We don't care about the order of random-points.grep values we get.)

my \batches = (^$batch-count).race(batch => 1);

The reason for the batch argument of 1 is that we know each random-points.grep takes a long time, so we don't want it to batch them together. We also don't want to wait for 64 of them to be calculated before we get any results.

Now to turn that into the sequence of values that we want by using map.

my \batches = (^$batch-count).race(batch => 1).map: {
    +random-points($batch-size).grep(*.is-inside-unit-circle)
}

Rather than have an infinite loop, we will loop over those batches

gather for batches -> Int $pointsInside {

So the last part of the file now looks like this:

sub compute-pi(UInt $batch-size, UInt $batch-count --> Seq) {

    my \batches = (^$batch-count).race(batch => 1).map: {
        +random-points($batch-size).grep(*.is-inside-unit-circle)
    }

    my Num $total = 0e0;
    my Num $count = 0e0;

    gather for batches -> Int $pointsInside {

        $total += $batch;
        $count += $pointsInside;
        my Num $ratio = $count / $total;

        take $ratio * 4;
    }
}

{ say "π ≈ $_" } for compute-pi(100_000, 500);

It now takes just over a minute (1m13.905s) with a cpu time of over 5 minutes (5m19.237s).
Which means we have more than 400% utilization. Which is fairly good on a 4 core CPU, with hyperthreading.

---

We can gain a bit more if we use native num variables instead of boxed Num variables.

You have to be careful with this as you can actually make it take more time if the value ends up needing to be boxed for an operation.

class Point {
    has num $.x;
    has num $.y;
    ...
}

...

sub compute-pi(UInt $batch-size, UInt $batch-count --> Seq) {

    my \batches = (^$batch-count).race(batch => 1).map: {
        random-points($batch-size).grep(*.is-inside-unit-circle).Num
    }

    my num $total = 0e0;
    my num $count = 0e0;

    gather for batches -> num $pointsInside {

        $total += $batch;
        $count += $pointsInside;
        my num $ratio = $count / $total;

        take $ratio * 4;
    }
}

That final change brings it down to about 42 seconds (0m42.771s) with 4 minutes of CPU time (4m0.929s).
Which also seems to increase the CPU utilization up over 500%.

---

The biggest improvement was going from gather for take to just xx. Which I'd argue was warranted just for how much it improved the clarity of the code. The other changes were of minimal benefit.

I think to get any significant improvement would require improving the compiler, runtime, or the VM.

Really what we ended up with is fairly fast considering how much of higher level language Raku is than C#, and how few people are working on it.

Answer 2

My original run time with your sample code is a consistent 5 minutes, 10 seconds.

You seem to have stumbled across a performance regression or something with exponentiation, because if I change $!x²+$!y² <= 1 to (($!x * $!x) + ($!y * $!y)) <= 1 it cuts 40 seconds off down to 4 minutes, 30 seconds. I'll see about filing an issue over that.

Your function compute-pi creates a Seq, and inside it calls random-points to make a Seq. So your ending compute-pi(100_000).head(500) will run the 'gather' construct inside random-points 500 times, generating a Seq 100_000 units long each time.

If instead you just have it return a number representing the number of randomly generated points that fit into a circle, on my machine that cuts execution time from just over 4 minutes, 30 seconds to 1 minute, 40 seconds:

sub inside-points(Int $count --> Num) {
  my Num $tot = 0e0;
  for 0..^$count {
    if Point.new(x => 1e0.rand, y => 1e0.rand).is-inside-unit-circle() {
      $tot++;
    }
  }   
  $tot;
}
# and then further down, inside compute-pi:
my Num $pointsInside = inside-points($batch);
$total += $batch;
$count += $pointsInside;

The next swap to make for performance is to swap out Raku's "Num" - similar to Double in Java or C# - for num64, which is a native 'double' in C# or Java. You can just do a find/replace on the code and switch Num to num64. That cuts the execution time for me from 100 seconds to 50. I'm not going to post that change here.

The final change is to eliminate object creation. You can remove the Point type entirely, and change the function random-points to:

sub inside-points(Int $count --> num64) {
  my num64 $tot = 0e0;
  my num64 $x; 
  my num64 $y; 
  for 0..^$count {
      $x = 1e0.rand;
      $y = 1e0.rand;
      # again, avoid the performance hit of $x²+$y²
      if (($x * $x) + ($y * $y)) <= 1e0 {
        $tot++;
      }
  }
  $tot;
}

That cuts the execution time to 6 seconds on my machine.

If performance is an issue, you can write Raku code in this kind of imperative style and it seems to work pretty well. But that said, I like your original implementation better. Your use of a Point type and sequences makes the intent easier to follow than in mine.

Edit: bonus round, I decided to try to see if I could speed it further by using concurrency. This version runs in about 4 seconds on my machine. Uncommenting the lines with the $totalcalcs variable will use it to prove that the program is still running properly, meaning it doesn't skip any loops.

I use asynchronous promises to breakup the calculation of random values into batches. I tried different batch sizes and 8_000 seems to give the best result.

# my atomicint $totalcalcs = 0;
sub inside-points-wrapped(Int $count --> num64) {
    my num64 $tot = 0e0;
    my num64 $x;
    my num64 $y; 
    for 0..^$count {
        $x = 1e0.rand;
        $y = 1e0.rand;
        if (($x * $x) + ($y * $y)) <= 1e0 {
          $tot++;
        }
        # $totalcalcs⚛++;
    }
    $tot;
}

my \batch_size = 8_000;
sub inside-points(int $count --> num64) {
    my $loops = $count div batch_size;
    my $rem = $count mod batch_size;
    my @inner_batches = batch_size xx $loops;
    @inner_batches.push($rem);
    my Promise @promises;
    for @inner_batches -> $batch {
      @promises.push(start { inside-points-wrapped($batch); });
    }
    my @result = await @promises; 
    [+] @result;
}

sub compute-pi(Int $batch --> Seq) {
    my num64 $total = 0e0;
    my num64 $count = 0e0;

    gather for 0..^∞ {
        my num64 $pointsInside = inside-points($batch);
        $total += $batch;
        $count += $pointsInside;
        my num64 $ratio = $count / $total;
        take $ratio * 4e0;
    }
}
{ say "π ≈ $_" } for compute-pi(100_000).head(500);
# say $totalcalcs;

Answer 3

Though the main question is how to have Raku perform faster, I am going to address performance with the C# code. Do not mark this an the accepted answer unless you can apply similar changes to your Raku code.

You added the C# code for completeness. I will not offer a CR of syntax or style but will address performance. On my particular PC, your code runs in 4 seconds but with the changes below it runs consistently less than 1.3 seconds.

Use simple squaring instead of Math.Pow

internal bool IsInsideUnitCircle() => (_x * _x) + (_y * _y) <= 1;

PointGenerator could be a static class

You are calling this 500 times, which means _random is created 500 times. Creating a new Random() is not a super fast operation. If you intentionally wanted a new random every 200_000 times (that is 100_000 iterations calling NextDouble() twice per Point). If you think you can live with seeding _random only once, you can make this class and its methods & properties static too.

Some LINQ is sluggish

In the ComputerPi method, I will avoid 2 LINQ calls: Where() and Count(). In my example, I am using a static PointGenerator class. You can change this to non-static as you originally had.

private static IEnumerable<double> ComputePi(int count)
{
    int total = 0;  // not a double like the original post
    int inside = 0; 
    while (true)
    {
        foreach (var point in PointGenerator.GeneratePoints(count))
        {
            // Avoid using the LINQ Where() and Count() methods
            if (point.IsInsideUnitCircle())
            {
                inside++;
            }
        }

        total += count;
        yield return 4.0 * (double)inside / (double)total;
    }
}

Do you want to time your entire app or just the computing PI part?

You are currently measuring how long your entire application runs. If you wanted to measure just the part where PI is being estimated, you could alter Main:

// I will only time how long it takes to estimate PI, not how long it takes to write results to the console.
var sw = Stopwatch.StartNew();
var estimates = ComputePi(100_000).Take(500).ToList();
sw.Stop();

for (var i = 0; i < estimates.Count; i++)
{
    Console.WriteLine($"trial = {i}, π ≈ {estimates[i]}, delta = {(Math.PI - estimates[i])}");
}

Console.WriteLine($"Elapsed time: {sw.Elapsed}");

Note I display a little something extra to the console window. This does have a blocking issue in that ALL calculations will have to completely run BEFORE you see anything at the console. This could impart a false sense that the application is frozen.

Again, altering Main is purely optional but it did shave 0.1+ seconds off my timings (relative to my PC). You could try this a few times on your PC to measure the difference, and then return to your original Main so that you have the benefit of scrolling messages as they occur.

Little things add up

These are all little things but LOTS of little things done wrong can add up to slower performance just like LOTS of little things done right can add to performance savings. In my case, I was able to run your Monte Carlo simulation in 1/3 the time.

Answer 4

You can speed it up as follows

use v6;

class Point {
    has Num $.x;
    has Num $.y;

    method is-inside-unit-circle(--> Bool) {
        $.x²+$.y² <= 1;
    }
}

sub random-points(Int $count --> RaceSeq) {
    [^$count].race.map({ Point.new(x => 1.rand, y => 1.rand) });
}

sub compute-pi(Int $batch --> Seq) {
    my Num $total = 0e0;
    my Num $count = 0e0;

    gather while True {
        my RaceSeq $pointsInside = random-points($batch).race.grep: *.is-inside-unit-circle;

        $total += $batch;
        $count += $pointsInside.elems;
        my Num $ratio = $count / $total;

        take $ratio * 4;
    }
}

say "π ≈ $_" for compute-pi(100_000).head(500);

Output:

real    2m53,695s
user    9m16,476s
sys     0m6,780s