1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
use crate::{Normal, set_state, update_and_uniform};
use std::cell::Cell;
impl Normal {
/// コンストラクタ
/// * `_seed_1` - 乱数の種
/// * `_seed_2` - 乱数の種。`_seed_1`と同じ値の場合、コンストラクタ側で変更する。
pub fn new(_seed_1: u32, _seed_2: u32) -> Self {
let _seed_other = if _seed_1 != _seed_2 { _seed_2 } else { (_seed_1 as u64 + 1192u64) as u32};
Self {
xyzw_1: set_state(_seed_1),
xyzw_2: set_state(_seed_other),
even_flag: Cell::<bool>::new(false),
even_result: Cell::<f64>::new(0f64),
}
}
/// 標準正規分布に従う乱数を返す
/// * 平均値 0
/// * 標準偏差 1
pub fn sample(&self) -> f64 {
// アルゴリズム 3.2
// step 1 & 5: 偶数回目の乱数は、奇数回目で計算したもう一つの値を返す
if self.even_flag.get() {
self.even_flag.set(false);
self.even_result.get()
}
else {
loop {
// step 2: 独立な一様乱数を2個生成する
let u1: f64 = update_and_uniform(&self.xyzw_1);
let u2: f64 = update_and_uniform(&self.xyzw_2);
// step 3: 中間変数を生成する
let v1 = 2f64 * u1 - 1f64;
let v2 = 2f64 * u2 - 1f64;
let v = v1.powi(2) + v2.powi(2);
// step 4: 0 < v < 1 のとき、戻り値計算に移る
if 0f64 < v && v < 1f64 {
let w: f64 = (-2f64 * v.ln() / v).sqrt();
// step 5: 計算した乱数を返す
self.even_result.set(v2 * w); // y2
self.even_flag.set(true);
return v1 * w; // y1
}
}
}
}
}
#[macro_export]
/// 正規分布のインスタンスを生成するマクロ
/// * `() =>` - 乱数の種は自動生成
/// * `($seed_1: expr, $seed_2: expr) =>` - 乱数の種を指定する
macro_rules! create_normal {
() => {{
let seeds: (u32, u32) = $crate::create_seeds();
$crate::Normal::new(seeds.0, seeds.1)
}};
($seed_1: expr, $seed_2: expr) => {
$crate::Normal::new($seed_1 as u32, $seed_2 as u32)
};
}