#/***********************************************************
#    gamma.rb -- ガンマ関数
#             -- ベータ関数
#***********************************************************/

DBL_DIG = 15                   # <float.h>
PI = 3.14159265358979324       # $\pi$ 
LOG_2PI = 1.83787706640934548  # $\log 2\pi$ 
N = 8

B0 = 1                  # 以下はBernoulli数 
B1 = (-1.0 / 2.0)
B2 = ( 1.0 / 6.0)
B4 = (-1.0 / 30.0)
B6 = ( 1.0 / 42.0)
B8 = (-1.0 / 30.0)
B10 = ( 5.0 / 66.0)
B12 = (-691.0 / 2730.0)
B14 = ( 7.0 / 6.0)
B16 = (-3617.0 / 510.0)

def loggamma(x)  # ガンマ関数の対数 
    v = 1
    while (x < N);  v *= x;  x += 1;  end
    w = 1 / (x * x).to_f
    return ((((((((B16 / (16 * 15))  * w + (B14 / (14 * 13))) * w\
                + (B12 / (12 * 11))) * w + (B10 / (10 *  9))) * w\
                + (B8  / ( 8 *  7))) * w + (B6  / ( 6 *  5))) * w\
                + (B4  / ( 4 *  3))) * w + (B2  / ( 2 *  1))) / x\
                + 0.5 * LOG_2PI - Math::log(v) - x + (x - 0.5) * Math::log(x)
end

def gamma( x)  # ガンマ関数 
    if (x < 0)
        return PI / (Math::sin(PI * x) * Math::exp(loggamma(1 - x))).to_f
    end
    return Math::exp(loggamma(x))
end

def beta(x,  y)  # ベータ関数 
    return Math::exp(loggamma(x) + loggamma(y) - loggamma(x + y))
end

printf("  x       Gamma(x)\n")
(-5.5).step(0.51, 1.0) do |x|
    printf("%4.1f  % .*g\n", x, DBL_DIG, gamma(x))
end
1.step(5.1, 1.0) do |x|
    printf("%4.1f  % .*g\n", x, DBL_DIG, gamma(x))
end
(10).step(30.1, 5.0) do |x|
    printf("%4.1f  % .*g\n", x, DBL_DIG, gamma(x))
end

exit 0