func main()
    println("Compute Intergal_a^b f(x)dx using Simpson formula");
    print("a = ");
    a = double(scan());
    print("b = ");
    b = double(scan());
    print("n = ");
    n = int(scan());

    s0 = simpson(sin, a, b, n);
    println("integral_a^b sin(x) dx = ", s0);

    s1 = simpson(cube, a, b, n);
    println("integral_a^b x^3 dx = ", s1);
endfunc

function cube(x)
    return x^3;
endfunc

function simpson(f, a, b, n)
    k = 2*n;
    h = (b - a)/k;
    s = f(a);
    x = a;
    while (0 == 0)
        x = x + h;
        s = s + 4.*f(x);
        x = x + h;
        if (x >= b - h/2.)
            break;
        endif
        s = s + 2.*f(x);
    endwhile
    s = s + f(b);
    return s*h/3.;
endfunc

# Starter
main();