#/***********************************************************
#    irandom.rb -- 幾何分布
#               -- 2項分布
#               -- Poisson (ポアソン) 分布
#***********************************************************/

#**** 一様乱数 (線形合同法) *****************************

ULONG_MAX = 4294967295
ULONG_MAX1 = 4294967296
$seed = 1;  # 任意 

def init_rnd(x)
    $seed = x
end

def irnd
    $seed = ($seed * 1566083941) % ULONG_MAX1 + 1
    return $seed
end

def rnd  # 0 <= rnd() < 1 
    return (1.0 / (ULONG_MAX + 1.0)) * irnd()
end

#********************************************************

#include <math.h>

def geometric_rnd1(p)  # 幾何分布 1 
    n = 1
    while (rnd() > p); n += 1; end
    return n
end

def geometric_rnd(p)  # 幾何分布 2 
    return (Math::log(1 - rnd()) / Math::log(1 - p)).ceil
end

def binomial_rnd(n,  p)  # 2項分布 
    r = 0
    for i in 0...n
        if (rnd() < p); r += 1; end
    end
    return r
end

def Poisson_rnd(lambda)  # Poisson (ポアソン) 分布 
    lambda = Math::exp(lambda) * rnd()
    k = 0
    while (lambda > 1) 
        lambda *= rnd();  k += 1
    end
    return k
end

#********************************************************

histo = []

init_rnd(Time.new.to_i)  # 初期化 

printf("***** メニュー *****\n")
printf("  1: 幾何分布 1\n")
printf("  2: 幾何分布 2\n")
printf("  3: ２項分布\n")
printf("  4: ポアソン分布\n")
printf("? ");  choice = gets.to_i

case choice 
when 1, 2, 4
    printf("引数 (1個)? ");  a = gets.to_f
when 3
    printf("引数 (2個)? ")
    a, b = gets.scan(/\S+/)
    a = a.to_f; b = b.to_f
end
printf("個数? ");  n = gets.to_i

for i in 0...20; histo[i] = 0; end
s1 = s2 = 0
for i in 0...n
    case choice 
    when 1
        x = geometric_rnd1(a)
    when 2
        x = geometric_rnd(a)
    when 3
        x = binomial_rnd(a, b)
    when 4
        x = Poisson_rnd(a)
    end
    if (x >= 0 && x < 20); histo[x] += 1; end
    s1 += x;  s2 += x * x
end

for i in 0...20
    printf("%4d: %5.1f%%\n", i, 100.0 * histo[i] / n)
end
s1 /= n;  s2 = Math::sqrt((s2 - n * s1 * s1) / (n - 1))
printf("平均 %g  標準偏差 %g\n", s1, s2)

exit 0