function rect_left(integer rid, atom x, atom h)
    if atom(h) then end if  -- suppress warning
    return call_func(rid,{x})
end function

function rect_mid(integer rid, atom x, atom h)
    return call_func(rid,{x+h/2})
end function

function rect_right(integer rid, atom x, atom h)
    return call_func(rid,{x+h})
end function

function trapezium(integer rid, atom x, atom h)
    return (call_func(rid,{x})+call_func(rid,{x+h}))/2
end function

function simpson(integer rid, atom x, atom h)
    return (call_func(rid,{x})+4*call_func(rid,{x+h/2})+call_func(rid,{x+h}))/6
end function

function cubed(atom x)
    return power(x,3)
end function

function recip(atom x)
    return 1/x
end function

function ident(atom x)
    return x
end function

function integrate(integer m_id, integer f_id, atom a, atom b, integer steps)
atom accum = 0,
     h = (b-a)/steps
    for i=0 to steps-1 do
        accum += call_func(m_id,{f_id,a+h*i,h})
    end for
    return h*accum
end function

function smartp(atom N)
string res
    if N=floor(N) then return sprintf("%d",N) end if
    res = sprintf("%12f",round(N,1000000))
    if find('.',res) then
        res = trim_tail(res,"0")
        res = trim_tail(res,".")
    end if
    return res
end function

procedure test(sequence tests)
string name
atom a, b, steps, rid
    printf(1,"Function     Range     Iterations       L-Rect       M-Rect       R-Rect      Trapeze      Simpson\n")
    for i=1 to length(tests) do
        {name,a,b,steps,rid} = tests[i]
        printf(1,"  %-5s %6d - %-5d %10d  %12s %12s %12s %12s %12s\n",{name,a,b,steps,
                            smartp(integrate(routine_id("rect_left"),   rid,a,b,steps)),
                            smartp(integrate(routine_id("rect_mid"),    rid,a,b,steps)),
                            smartp(integrate(routine_id("rect_right"),  rid,a,b,steps)),
                            smartp(integrate(routine_id("trapezium"),   rid,a,b,steps)),
                            smartp(integrate(routine_id("simpson"),     rid,a,b,steps))})
    end for
end procedure

constant tests = {{"x^3", 0,    1,     100, routine_id("cubed")},
                  {"1/x", 1,  100,    1000, routine_id("recip")},
                  {"x",   0, 5000, 5000000, routine_id("ident")},
                  {"x",   0, 6000, 6000000, routine_id("ident")}}

test(tests)