はじめに
Crystal言語で正規乱数が出したくなった。しかし、Crystalは数値計算を目的としたプログラミング言語ではないので標準ライブラリにはなかった。
Box-Muller法
ボックス=ミュラー法というものを使えば、一様分布に従う確率変数から、標準正規分布に従う確率変数を生成できるそうだ。
実装
Seedを設定して毎回同じ乱数が出るようにしたかった。
- 標準ライブラリのRandomをオープンクラスで拡張する
- 標準ライブラリのRandomを継承してRandom2を作る
- 新しくRandクラスを作り、知らないメソッドが呼ばれたらRandomに委譲する
の3つのパターンが考えられて、一番最後のものを採用することにした。
class Rand
def initialize
r = Random.new
initialize(random: r)
end
def initialize(seed : UInt64)
r = Random.new(seed)
initialize(random: r)
end
def initialize(random : Random)
@random = random
@iset = 0.to_i32
@gset = 0.to_f64
end
forward_missing_to @random
# Generates a random number from a standard normal distribution.
def randn : Float64
if @iset == 0
rsq, v1, v2, fac = 0.0, 0.0, 0.0, 0.0
loop do
v1 = 2.0 * @random.next_float - 1.0
v2 = 2.0 * @random.next_float - 1.0
rsq = v1 * v1 + v2 * v2
break if rsq < 1.0 && rsq != 0.0
end
fac = Math.sqrt(-2.0 * Math.log(rsq) / rsq)
@gset = v1 * fac
@iset = 1
return v2 * fac
else
@iset = 0
return @gset
end
end
# Generates a random number from a normal distribution with the specified mean and standard deviation.
def randn(mean, std_dev) : Float64
mean = mean.to_f
std_dev = std_dev.to_f
raise "std_dev must be positive" if std_dev.negative?
mean + std_dev * randn()
end
end
Cのコードなどを見ながら見よう見まねで書いている。何か間違っている部分を発見したらコメント欄で指摘してほしい。
ライブラリにした
このコードは再利用する機会がありそうなのでShardにした。Crystalは発展途上の言語なので、自分で書いたコードを再利用しやすいように積極的にライブラリとして公開していきたい。
この記事は以上です。