はじめに
モンテカルロ法
モンテカルロ法を使用して各言語で円周率を求めることにより各言語の速度比較を行います。
モンテカルロ法とは、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さんにお聞きしたところ、筆者の環境がバージョンが古いことが分かりましたので、後日バージョンを最新化して再計測する予定です。
本稿はたまーに更新します。
よさそうな言語がありましたら追加していこうと思います。