Language comparison benchmark: Monte Carlo

Monte Carlo

The area under the curve y = sin(x²) and above the x-axis for x values between 0 and √π should be determined. This should be done with a Monte Carlo integration. The program should check 10000000 random points inside a rectangle where x is in the range [0, 2) and y is in the range [0. 1). A random point (x, y) is inside the region if y < sin(x²) holds.

C

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

int main(int argc, char *argv[]) {
  struct timeval time;

  gettimeofday(&time,NULL);
  srand48((unsigned int) time.tv_usec);

  int i, N = 10000000, num = 0;
  double x, y;

  for (i = 0; i < N; i++) {
    x = 2 * drand48();
    y = drand48();
    if (y < sin(x*x)) num++;
  }

  printf("%f\n", 2.0 * num / N);
  return 0;
}

Measurement:

prompt> gcc -lm -O2 sinx2area.c -o sinx2area
prompt> time ./sinx2area
0.894774

real    0m0.431s
user    0m0.431s
sys     0m0.000s

C++

#include <cstdio>
#include <cstdint>
#include <random>
#include <chrono>

using namespace std;
using randgen_t = minstd_rand;
using int32 = int_fast32_t;

static constexpr const int32 N = 10'000'000;

int main (int argc, char *argv[]) {
	minstd_rand::result_type default_seed =
		std::chrono::system_clock::now().time_since_epoch().count();
	randgen_t rand_gen{default_seed};

	int32 n = 0;
	for (int32 i = 0; i < N; ++i) {
		double x = (2.0 * rand_gen()) / randgen_t::max();
		double y = (1.0 * rand_gen()) / randgen_t::max();
		if (y < sin(x*x)) {
			n++;
		}
	}
	
	printf("%f\n", (2.0 * n) / N);
}

Measurement:

prompt> g++ -O3 -flto sinx2area.cpp -lm -o sinx2area
prompt> time ./sinx2area
0.894844

real    0m0.290s
user    0m0.286s
sys     0m0.004s

Go

package main

import (
   "fmt"
   "math"
   "math/rand"
   "time"
)

const N = 10000000
var x, y float64
var num int = 0

func main() {
   rand.Seed(time.Now().UnixNano())
   for i := 0; i < N; i++ {
      x := 2.0*rand.Float64()
      y := rand.Float64()
      if y < math.Sin(x*x) {
         num++
      }
   }
   fmt.Printf("%.7f\n", 2.0*float64(num)/float64(N))
}

Measurement:

prompt> go build sinx2area.go
prompt> time ./sinx2area
0.8946618

real    0m0.608s
user    0m0.601s
sys     0m0.008s

Java

public class sinx2area {
   public static void main (String[] args) {
      int N = 10000000;
      int num = 0;
      for (int i = 0; i < N; i++) {
         double x = 2.0 * Math.random();
         double y = Math.random();
         if (y < Math.sin(x*x)) num += 1;
      }
      System.out.println(2.0 * num / N);
    }
}

Measurement:

prompt> time java sinx2area
0.8944406

real    0m0.526s
user    0m0.539s
sys     0m0.016s

JavaScript

(function(){
	var N = 10000000;
	var num = 0;
	for (var i = 0; i < N; i ++) {
		var x = 2 * Math.random();
		var y = Math.random();
		if (y < Math.sin(x * x))
			num++;
	}
	console.log(2 * num / N);
})();

Measurement:

prompt> time node sinx2area.js
0.8952678

real    0m0.475s
user    0m0.453s
sys     0m0.024s

Perl

#!/usr/bin/perl

$n=10000000;
$tot=0;

for($i=1;$i<=$n;$i++){
    $x=2*rand();
    $y=rand();
    
    if($y<sin($x*$x)){
        $tot++;
    }
}
$area=(2.0*$tot)/$n;
print "$area\n";

Measurement:

prompt> time perl sinx2area.pl
0.8940136

real    0m3.046s
user    0m3.038s
sys     0m0.000s

Python 2

#!/usr/bin/python
import random
import math

max=10000000
tot=0

for i in range(1,max+1):
    x=2*random.random()
    y=random.random()
    if y < math.sin(x*x):
        tot=tot+1

area=(2.0*tot/max)
print "%0.5f" %(area)

Measurement:

prompt> time python sinx2area.py 
0.89431

real    0m4.480s
user    0m4.402s
sys     0m0.076s

Ruby

#!/usr/bin/ruby

N = 10000000
num = 0
for i in 1 .. N do
    x = 2 * rand()
    y = rand()
    num += 1 if (y < Math.sin(x**2))
end
area = (2.0 * num) / N
puts(area)

Measurement:

prompt> time ruby sinx2area.rb
0.8949598

real    0m2.960s
user    0m2.947s
sys     0m0.013s

Seed7

$ include "seed7_05.s7i";
  include "float.s7i";
  include "math.s7i";

const proc: main is func
  local
    const integer: N is 10000000;
    var integer: num is 0;
    var float: x is 0.0;
    var float: y is 0.0;
  begin
    for N do
      x := rand(0.0, 2.0);
      y := rand(0.0, 1.0);
      num +:= ord(y < sin(x*x));
    end for;
    writeln(2.0 * flt(num) / flt(N) digits 7);
  end func;

Measurement:

prompt> s7c -O2 -oc3 sinx2area
prompt> time ./sinx2area
0.8947734

real    0m0.223s
user    0m0.223s
sys     0m0.000s