#/***********************************************************
#    random.rb -- カイ2乗分布
#              -- ガンマ分布
#              -- 三角分布
#              -- 指数分布
#              -- 正規分布
#              -- ベータ分布
#              -- 累乗分布
#              -- ロジスティック分布
#              -- ワイブル分布
#              -- Cauchy (コーシー) 分布
#              -- F分布
#              -- t分布
#***********************************************************/

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

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

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

PI = 3.141592653589793238
E =  2.718281828459045235

def triangular_rnd  # 三角分布 
    return rnd() - rnd()
end

def power_rnd1(n)  # 累乗分布 1 
    p = rnd()
    for i in 0...n
        if ((r = rnd()) > p); p = r; end
    end
    return p
end

def power_rnd(n)  # 累乗分布 2 
    return rnd() ** (1.0 / (n + 1).to_f)
end

def exp_rnd  # 指数分布 
    return -Math::log(1 - rnd())
end

def Cauchy_rnd1  # Cauchy (コーシー) 分布 1 
    return Math::tan(PI * rnd())
end

def Cauchy_rnd  # Cauchy (コーシー) 分布 2 
    begin 
        x = 1 - rnd();  y = 2 * rnd() - 1
    end while (x * x + y * y > 1)
    return y / x.to_f
end

def nrnd1  # 正規分布 1 
    return rnd() + rnd() + rnd() + rnd() + rnd() + rnd()\
         + rnd() + rnd() + rnd() + rnd() + rnd() + rnd() - 6
end

$sw_rnd2 = 0
def nrnd2  # 正規分布 2 
    if ($sw_rnd2 == 0) 
        $sw_rnd2 = 1
        $t = Math::sqrt(-2 * Math::log(1 - rnd()));  $u = 2 * Math::PI * rnd()
        return $t * Math::cos($u)
    else 
        $sw_rnd2 = 0
        return $t * Math::sin($u)
    end
end

$sw = 0; $r1 = 0; $r2 = 0; $s = 0
def nrnd  # 正規分布 3 
    if ($sw == 0) 
        $sw = 1
        begin 
            $r1 = 2 * rnd() - 1
            $r2 = 2 * rnd() - 1
            $s = $r1 * $r1 + $r2 * $r2
        end while ($s > 1 || $s == 0)
        $s = Math::sqrt(-2 * Math::log($s) / $s.to_f)
        return $r1 * $s
    else 
        $sw = 0
        return $r2 * $s
    end
end

def gamma_rnd1(two_a)  # ガンマ分布 
    x = 1
    (two_a / 2).downto(1) { x *= rnd() }
    x = -Math::log(x)
    if ((two_a & 1) != 0)   # if two_a is odd 
        r = nrnd();  x += 0.5 * r * r
    end
    return x
end

def gamma_rnd(a)  # ガンマ分布, a > 0 
    if (a > 1) 
        t = Math::sqrt(2 * a - 1)
        begin 
            begin 
                # 次の４行は y = tan(PI * rnd()) と同値 
                begin 
                    x = 1 - rnd();  y = 2 * rnd() - 1
                end while (x * x + y * y > 1)
                y /= x.to_f
                x = t * y + a - 1
            end while (x <= 0)
            u = (a - 1) * Math::log(x / (a - 1).to_f) - t * y
        end while (u < -50 || rnd() > (1 + y * y) * Math::exp(u))
    else 
        t = E / (a + E)
        begin 
            if (rnd() < t) 
                x = rnd() ** 1 / a.to_f;  y = Math::exp(-x)
            else 
                x = 1 - Math::log(rnd());  y = x ** (a - 1)
            end
        end while (rnd() >= y)
    end
    return x
end

def beta_rnd1(a,  b)  # ベータ分布 1 
    begin 
        x = rnd() ** (1 / a.to_f);  y = rnd() ** (1 / b.to_f)
    end while (x + y > 1)
    return x / (x + y).to_f
end

def beta_rnd(a,  b)  # ベータ分布 2 
    temp = gamma_rnd(a)
    return temp / (temp + gamma_rnd(b)).to_f
end

def chisq_rnd1( n)  # カイ2乗分布 
    s = 0
    for i in 0...n; t = nrnd(); s += t * t; end
    return s
end

def chisq_rnd( n)  # カイ2乗分布 
    return 2 * gamma_rnd(0.5 * n)
end

def F_rnd( n1,  n2)  # F分布 
    return (chisq_rnd(n1) * n2) / (chisq_rnd(n2) * n1)
end

def t_rnd(n)  # t分布 
    if (n <= 2); return nrnd() / Math::sqrt(chisq_rnd(n) / n.to_f); end
    begin 
        a = nrnd()
        b = a * a / (n - 2).to_f
        c = Math::log(1 - rnd()) / (1 - 0.5 * n).to_f
    end while (Math::exp(-b-c) > 1 - b)
    return a / Math::sqrt((1 - 2.0 / n.to_f) * (1 - b))
end

def logistic_rnd  # ロジスティック分布 
    r = rnd()
    return Math::log(r / (1 - r).to_f)
end

def Weibull_rnd(alpha)  # Weibull (ワイブル) 分布 
    return (-Math::log(1 - rnd())) ** (1 / alpha.to_f)
end

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

histo = []

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

printf("***** メニュー *****\n")
printf("  1: 三角分布\n")
printf("  2: 累乗分布 1\n")
printf("  3: 累乗分布 2\n")
printf("  4: 指数分布\n")
printf("  5: コーシー分布 1\n")
printf("  6: コーシー分布 2\n")
printf("  7: 正規分布 1\n")
printf("  8: 正規分布 2\n")
printf("  9: 正規分布 3\n")
printf(" 10: ガンマ分布 1\n")
printf(" 11: ガンマ分布 2\n")
printf(" 12: ベータ分布 1\n")
printf(" 13: ベータ分布 2\n")
printf(" 14: カイ２乗分布 1\n")
printf(" 15: カイ２乗分布 2\n")
printf(" 16: F分布\n")
printf(" 17: t分布\n")
printf(" 18: ロジスティック分布\n")
printf(" 19: ワイブル分布\n")
printf("? ");  choice = gets.to_i

case choice 
when 2,3,10,11,14,15,17,19
    printf("引数 (1個)? ");  a = gets
when 12,13,16
    printf("引数 (2個)? ")
    a, b = gets.scan(/\S+/)
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 = triangular_rnd(); puts "triangular"
    when 2; x = power_rnd1(a.to_i)
    when 3; x = power_rnd(a.to_f)
    when 4; x = exp_rnd()
    when 5; x = Cauchy_rnd1()
    when 6; x = Cauchy_rnd()
    when 7; x = nrnd1()
    when 8; x = nrnd2()
    when 9; x = nrnd()
    when 10; x = gamma_rnd1(a.to_i)
    when 11; x = gamma_rnd(a.to_f)
    when 12; x = beta_rnd1(a.to_f, b.to_f)
    when 13; x = beta_rnd(a.to_f, b.to_f)
    when 14; x = chisq_rnd1(a.to_i)
    when 15; x = chisq_rnd(a.to_f)
    when 16; x = F_rnd(a.to_f, b.to_f)
    when 17; x = t_rnd(a.to_f)
    when 18; x = logistic_rnd()
    when 19; x = Weibull_rnd(a.to_f)
    end
    ix =  (2 * x).floor + 10
    if (ix >= 0 && ix < 20); histo[ix] += 1; end
    s1 += x;  s2 += x * x
end

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

exit 0