#/***********************************************************
#    numint.rb -- 数値積分
#***********************************************************/

def f(x)  # 被積分関数 $f(x)$ 
    return 4 / (1 + x * x);  # 例 $4 / (1 + x^2)$ 
end

nmax = 32

a = 0
b = 1

printf("    n  台形           中点           Simpson\n")
h = b - a;  trapezoid = h * (f(a) + f(b)) / 2.0
n = 1; while (n <= nmax)
    midpoint = 0
    for i in 1..n; midpoint += f(a + h * (i - 0.5)); end
    midpoint *= h
    simpson = (trapezoid + 2 * midpoint) / 3.0
    printf("%5d % -14g % -14g % -14g\n",\
        n, trapezoid, midpoint, simpson)
    h /= 2.0;  trapezoid = (trapezoid + midpoint) / 2.0
    n *= 2 
end

exit 0