RUSTで円周率をモンテカルロ法で求めた

はじめに

RUSTが速くて安全で、とても良い言語だという評判なので、
試しに、RUSTでモンテカルロ法で円周率を求めるプログラムを作ってみた。
アルゴリズム自体は、よくある"hit-or-miss Monte Carlo"法である。
(0, 0), (1, 0), (0, 1) (1, 1)で囲まれた正方形の中にランダムに点を生成した際に、
原点から半径1以下にある確率は、PI/4となることを利用し、円周率PIを推定する。

乱数には、メルセンヌツイスターを利用。利用したクレートは、mt19937。
C++のBoostのような定番ライブラリーはまだ無い印象。（調査不足かもしれませんが。）

---

ソースコード

use std::env;

fn main() {
    let mut argv = Vec::<String>::new();
    for argument in env::args() {
        argv.push(argument);
    }
    let argc = argv.iter().len();

    if argc != 3 {
        println!("need exactly two arguments");
        println!("usage: /path/to/get_pi num_rand seed");
        std::process::exit(1);
    }

    let num_rand: u32 = argv[1].parse().unwrap();
    let seed: u32 = argv[2].parse().unwrap();

    let mut m = mt19937::MT19937::new_with_slice_seed(&[seed]);

    let mut count: u32 = 0;
    for _i in 0..num_rand {
        let x = mt19937::gen_res53(&mut m);
        let y = mt19937::gen_res53(&mut m);
        count += unity_if_in_circle(x, y);
    }

    let pi = count as f64 / num_rand as f64 * 4.0;
    println!("*****************************");
    println!(
        "With {} samples and random number seed of {}, ",
        num_rand, seed
    );
    println!("pi is estimated to be {}", pi);
    println!("*****************************");

    std::process::exit(0);
}

fn unity_if_in_circle(x: f64, y: f64) -> u32 {
    let r2 = x.powi(2) + y.powi(2);
    if r2 < 1.0 {
        1
    } else {
        0
    }
}

---

実行

ターミナルで./path/to/get_pi 100000000 1
と打てば、

*****************************
With 100000000 samples and random number seed of 1,
pi is estimated to be 3.14175024
*****************************

と結果が出力される。
seedを変えると、

*****************************
With 100000000 samples and random number seed of 2,
pi is estimated to be 3.14167144
*****************************

や

*****************************
With 100000000 samples and random number seed of 3,
pi is estimated to be 3.14157348
*****************************

と結果が得られる。
サンプル数の平方が4桁なので、おおよそ4桁程度の精度が出ている。
一桁精度を上げるには、100倍のサンプルが必要。

---

検討事項やコメントなど

env::argsの利用方法がよく分からない。
「clapを使え」というのはありそうだけれども。
たかだか2個の整数を読むだけなので、今回は素朴にenv::argsを利用してみた。
しかし、env::args（で得られるArgs）の取り回しがよく分からず、また、argcに対応するものも見つけられず。
最終的に、各引数は文字列のベクトルに格納して、そのベクトル長を調べてargcを定義。
もう少し良いやり方がある気がするのだが。

let pi = count as f64 / num_rand as f64 * 4.0という文は綺麗（そのまま英語で読める）。
C++だと(double)count/rand*4.0といった所かな。

関数のunity_if_in_circleは切り分けなくても、
count += if x.powi(2) + y.powi(2) < 1.0 { 1 } else { 0 } ;
で十分だったかな。C++で書くにしても、三項演算子で済ませたような気もします。