#/***********************************************************
#    monte.rb -- モンテカルロ法
#***********************************************************/

#**** 一様乱数 (線形合同法) *****************************
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

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

def monte1(n)
    hit = 0
    for i in 0...n
        x = rnd();  y = rnd()
        if (x * x + y * y < 1); hit += 1; end
    end
    p =  hit / n.to_f
    printf("pi = %6.4f +- %6.4f\n",
        4 * p, 4 * Math::sqrt(p * (1 - p) / (n - 1).to_f))
end

def monte2(n)
    sum = sumsq = 0
    for i in 0...n
        x = rnd()
        y = Math::sqrt(1 - x * x)
        sum += y;  sumsq += y * y
    end
    mean = sum / n.to_f
    sd = Math::sqrt((sumsq / n.to_f - mean * mean) / (n - 1).to_f)
    printf("pi = %6.4f +- %6.4f\n", 4 * mean, 4 * sd)
end

def monte3(n)
    a = (Math::sqrt(5) - 1) / 2.0

    x = sum = 0
    for i in 0...n
        if ((x += a) >= 1); x -= 1; end
        sum += Math::sqrt(1 - x * x)
    end
    printf("pi = %6.4f\n", 4 * sum / n.to_f)
end

init_rnd(123456789)
monte1(10000)

init_rnd(123456789)
monte2(10000)

monte3(10000)

exit 0