モンテカルロ法による各言語の速度比較

#はじめに

モンテカルロ法
モンテカルロ法を使用して各言語で円周率を求めることにより各言語の速度比較を行います。
モンテカルロ法とは、1.0>=x>=0, 1.0>=y>=0に任意の（乱数の）点を落とし、その落とされた点が原点(0,0）からの円内に落ちたものと、円外に落ちたものの比を求め、その比が円周率になるというものです。
もちろん、計算で求められた円周率に近似するものの、正確な値はでませんが、プログラミング言語で実装することには意味があると思っています。
理由としては、
・浮動小数点計算の速度が比較できる。
・Loop速度が比較できる。
・メモリの使用は少量であり、ページングは（おそらく）発生しない。
・DiskIOが発生しないため、外部記憶装置の影響を受けにくい。
・正解の円周率は既知の情報であり、デバッグが楽。
などがあります。

2020年6月28日 rubyを追加しました。
2020年6月28日 pythonの計算方法を変更しました。
2020年6月28日 rustのループ回数指定に誤りがあったので修正しました。
2020年6月28日 rustを@tatsuya6502 さんご提供のXoroshiro版のソースを追加測定
2020年7月1日 Juliaを追加しました。
2020年7月18日 Node.jsを追加しました。
2021年1月30日 D言語を追加しました。

以下、各言語ソース
#perl

perlは今まで試した中で非常に高速のイメージがあったため、ベンチマーク（≒基準）として、用意しました。
ソースもシンプルですので、オブジェクト指向等によるオーバーヘッドもなく、高速に動作するはずです。

montecarlo.pl

use Time::HiRes qw( gettimeofday tv_interval);

local $ave = 0;

local $time_start, $time_end, $time_split = 0;

local $max_times    = 10000000;
local $test_times   = 100;

local $cnt = 0;

for ( $i = 1 ; $i <= $test_times ; $i++ ) {

    print($i."\n");

    $time_start = [gettimeofday];
    $cnt = 0;
    for ( $j = 0 ; $j < $max_times ; $j++) {
        $x = rand();
        $y = rand();

        if ( $x * $x + $y * $y <= 1) {
            $cnt++;
        }
    }

    $time_split = tv_interval($time_start);
    print("pi ".(4 * $cnt / $max_times)."\n");
    print("split:".$time_split."\n");

    $ave += $time_split;
}

print("AVE:".($ave / $test_times)."\n");
~                                                                                                                    ~                                                                                                                    

#c言語
とりあえず、「Cなら速いんじゃね？」ということで用意しました。
厳密にチューニングすればもっと早くなるのかもしれませんが、
コーディングのしやすさでこのくらいで抑えてみました。

montecarlo.c

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

float montecarlo();

void main() {

    float ave = 0;
    int times = 100;

    for (int i = 1 ; i <= times ; i++){
        printf("%03d\n", i );
        ave += montecarlo();
    }

    printf("AVE %f\n\n", ave / times);

}

float montecarlo() {

    clock_t time_start, time_end;
    double  time_split;
    int times = 10000000;
    float x, y;

    time_start = clock();

    int cnt = 0;
    for ( int i = 0; i < times ; i++) {
        x = (float)rand() / (float)RAND_MAX;
        y = (float)rand() / (float)RAND_MAX;

        if ( (x * x + y * y) <= 1) {
            cnt++;
        }
    }

    time_end = clock();

    time_split = (double)(time_end - time_start) / 1000000;
    printf("%f\n", 4.0 * cnt / times);
    printf("time_split : %lf\n", time_split);

    return(time_split);

}

#go言語
goは前評判も良かったことから、perlと異なるベンチマークとして用意しました。
コンパイル言語ですので、c言語にどれだけ肉薄できるか？という意味でも興味があります。

montecarlo.go

package main

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

var i, j, cnt int
var x, y float64
var time_start, time_end int64
var ave, time_split float64

const MAX_TIMES = 10000000
const TEST_TIMES = 100

func main() {
        //fmt.Printf("%03d\n", 10)

        //fmt.Printf("%d\n", time.Now())
        //fmt.Printf("%f\n", float64(time.Now().UnixNano()) / 1000000000.0)
        //fmt.Printf("%f\n", float64(time.Now().UnixNano()) / 1000000000.0)
        //fmt.Printf("%f\n", float64(time.Now().UnixNano()) / 1000000000.0)

        rand.Seed(time.Now().UnixNano())

        ave := 0.0
        for i := 1 ; i <= TEST_TIMES ; i++ {

                fmt.Printf("%03d times\n", i)

                time_start      := time.Now().UnixNano()
                cnt := 0
                for j := 0 ; j <  MAX_TIMES ; j++ {

                        x := rand.Float64()
                        y := rand.Float64()

                        if x * x + y * y <= 1 {
                                cnt++
                        }

                }
                fmt.Printf("pi : %f\n", 4 * float64(cnt) / float64(MAX_TIMES))

                time_end        := time.Now().UnixNano()
                time_split      := float64(time_end - time_start) / 1000000000
                fmt.Printf("time : %f\n", time_split)
                ave += time_split
                //ave := time_split

        }

        fmt.Printf("AVE : %f\n", ave / TEST_TIMES)

}

#Java

コンパイル言語としては現状で鉄板に位置づけられると思います。（各々ご意見はあるでしょうが・・・。）
VMが間に入ることによりC言語よりは遅いことが予想されますが、過去の経験では致命的なほど遅いわけではないという印象です。

montecarlo.java

import java.util.Random;
import java.util.*;


class Montecarlo {
    public static void main(String[] args) {

        //System.out.println( getNowTime() );

                // 円の半径を指定する
                double iR = 1.0;

                // 描画回数を指定する
                int MAX_POINT = 10000000;

        // テスト回数を指定する
        int TEST_TIMES = 100;

                // Random オブジェクトの定義
                Random rand = new Random();

        // 開始、終了時刻
        long time_start, time_end = 0;

        double time_split;
        double ave = 0.0;

                // Random値の待避用変数の定義
                double ranValue = 0.0;
                double POS_X = 0.0;
                double POS_Y = 0.0;

        for ( int j = 1 ; j <= TEST_TIMES ; j ++) {

            //System.out.println("%03d", j);
            System.out.printf("%03d\n", j);

            // 開始時刻
            time_start = System.nanoTime();
                //System.out.println( "start : " + time_start );

            // 円弧状に落ちた件数をカウントする
            int HIT_COUNT = 0;

            // 出力用ワーク変数の初期化
            for ( int i = 0 ; i < MAX_POINT ; i++ ) {

                POS_X = rand.nextDouble();
                POS_Y = rand.nextDouble();

                if ( iR >= POS_X * POS_X + POS_Y * POS_Y )  {
                    HIT_COUNT++;
                }

            }
            // System.out.println(HIT_COUNT );
            System.out.println( (  4 * HIT_COUNT  / (double)MAX_POINT ) );

            // 終了時刻
            time_end = System.nanoTime();

                //System.out.println( "end   : " + time_end );

            time_split = (double)(time_end - time_start) / 1000000000.0;

                //System.out.println( "split : " + (time_end - time_start) );
                System.out.println( "split : " + time_split );

            ave += time_split;

        }

        //System.out.println( getNowTime() );
        System.out.println( "AVE : " + ave / TEST_TIMES );

    }

        /// 現在時刻をyyyy/mm/dd hh:mm:ss(jpn) 形式で取得する
        public static String getNowTime() {

                return(String.format("%1$tF %1$tT" ,Calendar.getInstance()));

        }

}

#lua
script言語としては軽いと評判ですので加えてみました。
言語仕様としてはOSに近しい命令を使用するため早いというイメージです。

montecarlo.lua

local socket = require("socket")

math.randomseed(os.time())

MAX_LOOP = 10000000
TEST_TIMES = 100

ave = 0.0
for i = 1, TEST_TIMES do
    print(string.format("%03d times", i))
    time_start  = socket.gettime()

    cnt = 0
    for j = 1 , MAX_LOOP do
        x = math.random()
        y = math.random()

        if x * x + y * y <= 1 then
            cnt = cnt + 1
        end

    end

    time_end    = socket.gettime()

    time_split  = time_end - time_start

    print(string.format("pi    : %f", 4 * cnt / MAX_LOOP))
    print(string.format("split : %f", time_split))

    ave = ave + time_split

end

print(string.format("ave : %f", ave / TEST_TIMES))

#python
言わずと知れたpythonです。
今回は2系、3系で比較するとともに、pypyでの比較を行います。

montecarlo.py

import random
import time
import math

_times = 10000000

def pi():
    _time_start = time.time()

    _cnt = 0
    for i in range(_times):

        _x = random.random()
        _y = random.random()

        #if  _x ** 2 + _y **2 <= 1:
        if  _x * _x + _y * _y <= 1:

            _cnt += 1

    print( 4.0 * _cnt / _times )

    _time_end = time.time()

    _time_split = _time_end - _time_start
    print(_time_split)

    return _time_split

_test = 100

_ave = 0.0
for _i in range(_test):
    print('{:03d} times'.format(_i+1))
    _ave += pi()

print('ave: {} s'.format(_ave/_test))

#Rust
RustはAI,機械学習で注目されている言語とのことです。
pythonよりは早く、Cよりは開発しやすい言語とのことです。

Cargo.toml

[package]
name = "montecarlo"
version = "0.1.0"
authors = ["ubuntu"]
edition = "2018"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[dependencies]
rand = "0.6"
floating-duration="0.1.2"

src/main.rs

fn main() {

    use rand::Rng;                  // 乱数用
    use std::time::Instant;         // 時刻取得用
    use floating_duration::TimeAsFloat; // 時刻→floatへの変換用

    println!("Hello, world!");

    let loop_max = 10000000;        // 打点の回数
    let test_max = 100;             // 平均を出すためのループ

    let mut x: f32;
    let mut y: f32;

    let mut rng = rand::thread_rng();//randseedの設定

    let mut split_time: f32 = 0.0;
    //for _test_count in 1..test_max {
    for _test_count in 1..=test_max {
        println!("times:{:?}", _test_count );
        let start_time = Instant::now();

        let mut count = 0;
        //for _loot_count in 1..loop_max {
        for _loot_count in 1..=loop_max {
            x = rng.gen_range(0., 1.);  //乱数0.0~1.0
            y = rng.gen_range(0., 1.);

            if x * x + y * y <= 1.0 {
                count = count + 1;
            }

        }
        println!("pi:{}", 4.0 * count as f32 / loop_max as f32); // as f32は型キャスト
        let end_time = start_time.elapsed().as_fractional_secs();// 時間を秒に変換
        println!("time:{:?}", end_time );

        split_time = split_time + end_time as f32;

    }

    println!("AVE:{}", split_time / test_max as f32);

}

@tatsuya6502 さんご提供のXoroshiro版ソース

Cargo.toml

[package]
name = "montecarlo"
version = "0.1.0"
authors = ["ubuntu"]
edition = "2018"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[dependencies]
rand = "0.7"
rand_xoshiro = "0.4"

src/main.rs

// 乱数用
use rand::Rng;
use rand_xoshiro::{rand_core::SeedableRng, Xoroshiro128Plus};

// 時刻取得用
use std::time::{Instant, SystemTime};

const LOOP_MAX: usize = 10_000_000; // 打点の回数
const TEST_MAX: usize = 100; // 平均を出すためのループ

fn main() -> Result<(), Box<dyn std::error::Error>> {
    // random seedを設定
    let now = SystemTime::now().duration_since(SystemTime::UNIX_EPOCH)?;
    let mut rng = Xoroshiro128Plus::seed_from_u64(now.as_millis() as u64);

    let mut split_time: f32 = 0.0;
    for _test_count in 1..=TEST_MAX {
        println!("times:{}", _test_count);
        let start_time = Instant::now();

        let mut count = 0;
        for _loot_count in 0..LOOP_MAX {
            // 乱数を生成 [0.0, 1.0) ・・ 0.0以上 1.0未満
            let x: f32 = rng.gen_range(0.0, 1.0);
            let y: f32 = rng.gen_range(0.0, 1.0);

            if x * x + y * y <= 1.0 {
                count += 1;
            }
        }
        println!("pi:{}", 4.0 * count as f32 / LOOP_MAX as f32); // as f32は型キャスト
        let end_time = start_time.elapsed().as_secs_f32(); // 時間を秒に変換
        println!("time:{}", end_time);

        split_time += end_time as f32;
    }

    println!("AVE:{}", split_time / TEST_MAX as f32);
    Ok(())
}

#Nim
比較用に・・・

montecarlo.nim

import random
import times

randomize()

const TEST_MAX = 10000000
const LOOP_MAX = 100

var x, y = 0.0

var ave = 0.0

for loop in 1..LOOP_MAX:
    echo $loop & " times"

    var time_start = cpuTime()
    var count = 0
    for test in 1..TEST_MAX:

        x = rand(0.0..1.0)
        y = rand(0.0..1.0)

        # echo $x
        # echo $y

        if ( x * x + y * y <= 1):
            count += 1

    echo "pi:" &  repr(float(4.0 * float(count) / float(TEST_MAX)))

    var time_split = cpuTime() - time_start

    echo "split" & repr(time_split)
    ave += time_split

echo "AVE:" & repr(ave / float(LOOP_MAX))
                                                                                                                    ~                                                                                                                    

#Ruby
@scivola さんにソースを提供していただきましたのでrubyを追加しました。

montecarlo.rb

TIMES = 10000000

def pi
  time_start = Time.now

  cnt = 0
  TIMES.times do
    x = rand
    y = rand

    if  x * x + y * y <= 1
      cnt += 1
    end
  end

  puts 4.0 * cnt / TIMES

  time_end = Time.now

  time_split = time_end - time_start
  puts time_split

   time_split
end

test = 100

ave = 0.0

test.times do |i|
    puts "%03d times" % (i + 1)
    ave += pi
end

puts "ave: #{ ave / test } s"

#Julia
もう一つ伸びそうな気がするのですが、チューニングのポイントがわかりません。

montecarlo.jl

using Printf: @printf, @sprintf
using Random


function montecarlo()
    local MAX_TIMES::Int32
    local TEST_TIMES::Int32

    local ave::Float32

    MAX_TIMES = 100
    TEST_TIMES = 10000000

    ave = 0.0
    time_sum = @elapsed for i = 1:MAX_TIMES

        @printf("times %03d\n", i )
        cnt = 0
        time_split = @elapsed for j = 1:TEST_TIMES

             x = rand(Float32)
             y = rand(Float32)

             if x * x + y * y <=1
                 cnt = cnt + 1
             end
        end

        @printf("pi: %f\n", 4 * cnt / TEST_TIMES)
        @printf("split : %f\n", time_split )

    end

    @printf("\n")
    return time_sum  / MAX_TIMES
end

@printf("AVE : %f\n", montecarlo())

#Node.js
@Kogia_sima さんよりソースを提供していただきましたので、Node.jsの計測を追加しました。

montecarlo.js

const performance = require('perf_hooks').performance;

const LOOP_MAX = 10000000;
const TEST_MAX = 100;

function pi() {
    const start = performance.now();

    let cnt = 0;

    for (let i = 0; i < LOOP_MAX|0; ++i) {
        const x = Math.random();
        const y = Math.random();

        if (x * x + y * y <= 1.0) {
            cnt += 1;
        }
    }

    const end = performance.now();
    const elapsed = (end - start) / 1000.0;

    console.log(4.0 * cnt / LOOP_MAX);
    console.log(elapsed);
    return elapsed;
}

if (require.main === module) {
    let ave = 0.0;

    for (let i = 1; i <= TEST_MAX; ++i) {
        console.log(`${i} times`);
        ave += pi();
    }

    console.log(`AVE: ${ave / TEST_MAX}`)
}

#D言語
@devmynoteさんにソース提供していただきましたのでD言語の計測を追加しました。

monte1.d

// ldc2 -m64 -release -O monte1.d
import std.algorithm;
import std.datetime;
import std.random;
import std.range;
import std.stdio;

void main()
{
    double    ave   = 0.0;
    immutable times = 100;
    for ( int i = 1; i <= times; i++ ){
        writefln("times : %03d", i);
        ave += montecarlo();
    }
    writefln("AVE   : %.8f\n", ave / times);
}

double montecarlo()
{
    immutable times = 10_000_000;
    long time_start = Clock.currStdTime();                  // 戻り値の単位hnsecs
    auto rnd = Random(cast(uint)time_start);
    long cnt =
        generate!(() => [uniform01(rnd), uniform01(rnd)])   // 乱数を2つ生成
        .takeExactly(times)                                 // times回繰り返す
        .map!(a => a[0] * a[0] + a[1] * a[1])               // それぞれ2乗して足して
        .count!("a <= b")(1.0);                             // 1.0以下となる場合をカウント
    long time_end = Clock.currStdTime();
    double time_split = cast(double)(time_end - time_start) / 10_000_000;   // 時間を秒に変換
    writefln("pi    : %.8f", 4.0 * cnt / times);
    writefln("split : %.8f", time_split);

    return( time_split );
}

AVE : 7.62902263

#結果

perl 5.30
	AVE:8.93247789s

java 14.0.1
	AVE: 1.4161036380199996s

c
	gcc 9.3.0

	gcc montecarlo.c -o pi

	gcc 最適化なし
	AVE:1.501472s

	gcc 実行速度の最適化 -O3
	AVE:1.425400s

python 2.7.18
	AVE:14.216013s

pypy 7.3.1 with GCC 9.3.0
	ave: 0.983056025505

python 3.8.2
	AVE:12.485190s
pypy3 7.3.1 with GCC 9.3.0
	ave: 1.1475282835960388 s

go 1.13.8
	AVE:2.132493

lua 5.3
	AVE:6.715329

rust (cargo 1.44.0)
	cargo run --release
	AVE 0.507196783
	
	cargo run
	AVE 16.449156

Nim 1.0.6
    nim c -r montecarlo.nim
    AVE 7.2063023393
    
    nim c -r -d:release montecarlo.nim
    AVE:0.3236947073200001

#2020年6月28日 追記

ruby を追加しました。
　Ruby 2.7.0dev (2019-12-25 master e1e1d92277) [aarch64-linux]
    ave: 7.5495061204099985 s

pythonのソースを修正しました。
 if  _x ** 2 + _y **2 <= 1:    ①
　↓
 if  _x * _x  + _y * _y <= 1:  ②

　※python 3.8.2の計測
　①再取得 
    ave: 11.132939925193787 s
　②
    ave: 9.153075976371765 s
　※pypy3の計測
　①再取得
    ave: 1.0760090994834899 s
　②
    ave: 1.3655683732032775 s
 あれ？

rustのループを修正しました。
    for _test_count in 1..test_max {
        for _loot_count in 1..loop_max {　①
     ↓
    for _test_count in 1..=test_max {
        for _loot_count in 1..=loop_max {　②

 cargo run --release の計測
　①再取得 
    AVE:0.5022585
　②
    AVE:0.51400733

rustで@tatsuya6502 さんご提供のXoroshiro版ソースについて計測しました。

　AVE:0.21163946

#2020年7月1日 追記

Juliaを追加しました。
  
  julia montecarlo.jl
    AVE : 0.616752

#2020年7月18日 追記

Node.jsを追加しました。
  node v10.21.0

  node montecarlo.js
    AVE: 1.0223767548998444

結論の順番が誤っていたため修正

正：pypy3 > c = java > go
誤：c = java >= pypy3

#2021年1月30日 追記

D言語を追加しました。
  binary    /usr/bin/ldc2
  version   1.12.0 (DMD v2.082.1, LLVM 6.0.1)
  config    /etc/ldc2.conf (armv6-unknown-linux-gnueabihf)
  OVERVIEW: LDC - the LLVM D compiler

  ldc2 -m32 -release -O monte1.d
  ./monte1
    AVE   : 7.15935545

#結論
このレベルであれは Nim rustが最速との評価になりました。
速度的には、

Nim > rust > pypy > C = java > pypy3 > go > lua > ruby > perl > python2系 > python3系
rust(Xoroshiro版) > Nim > rust > Julia > pypy > Node.js > pypy3 > C = java > go > lua > D言語 > ruby > perl > python2系 > python3系

となりました。
順当にコンパイル言語>スクリプト言語となっています。
もちろん、私のコーディングが悪く、「もっとこうすれば早くなる」というのはあるかもしれません。
しかし、多少工夫をしたところで1段階変わるだけで、さほど順位の入れ替わりがあるとは考えにくいです。

特筆すべきはpythonで、pypyを使用することによりC言語に近しい速度で動作する（可能性）が確認できました（もちろん、そのままのソースがpypyで動かない場合もありますし、ライブラリがpypyに対応していない場合もあります）。
しかしながらpythonの柔軟性のある記述で高速化できるのであれば試してみる価値があるとおもいます。

2020年6月28日追記
rubyが予想以上に速くなっており驚きを隠せません。以前試したときはc:ruby=1:20くらいでしたので、今回もそのくらいと予想していましたが見事に裏切られました。
@scivola さんに提案をいただかなければ無視していました。ありがとうございます。

pythonは今回のソース修正でpython3 で起動した場合は早くなるものの、pypy3で起動した場合は遅くなるという謎の結果になりました。（識者のコメント希望）
ちなみに本稿のソースではコメントして違いを表していますが、実測のソースでは置き換えていますので、コメントの評価で遅くなったということはありません。（昔のBasicとかでそういうのがあった気がしますが・・・）

rustは疑似乱数生成器の違いが如実に現れました。
@tatsuya6502 さんご提供のXoroshiro版ソースにより疑似乱数生成器も適材適所であることが分かりました。
ご指摘ありがとうございました。

2020年7月18日追記
@Kogia_simaさんよりソース提供を受け、Node.jsを追加しました。
Node.jsってこんなに早かったんだとビックリしています。
@Kogia_simaさんソース提供ありがとうございました。

（なぜjavascriptではなくNode.jsとしているかというと、今回はNode.jsから呼び出して計測しているからです。）
なお、この計測結果はUbuntuではなくrasbianでの計測結果となります。

2021年1月30日追記
@devmynoteさんよりソース提供を受け、D言語を追加しました。
@devmynoteさんソース提供ありがとうございました。
結果は、AVE : 7.15935545と、感覚的には非常に遅い結果となってしまいました。
（感覚的には、遅くてもgo並みと思っていました。）
@devmynoteさんにお聞きしたところ、筆者の環境がバージョンが古いことが分かりましたので、後日バージョンを最新化して再計測する予定です。

本稿はたまーに更新します。
よさそうな言語がありましたら追加していこうと思います。