Help us understand the problem. What is going on with this article?

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

はじめに

モンテカルロ法
モンテカルロ法を使用して各言語で円周率を求めることにより各言語の速度比較を行います。
モンテカルロ法とは、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を追加しました。

以下、各言語ソース

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())

結果

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

結論

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

Nim > rust > pypy > C = java > pypy3 > go > lua > ruby > perl > python2系 > python3系
rust(Xoroshiro版) > Nim > rust > Julia > pypy > C = java >= pypy3 > go > lua > 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版ソースにより疑似乱数生成器も適材適所であることが分かりました。
ご指摘ありがとうございました。

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

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした