while 1
    read x$
    if x$ ="end" then print "**Over**": end

    read a, b, N, knownValue

    print " Function y ="; x$; " from "; a; " to "; b; " in "; N; " steps"
    print " Known exact value ="; knownValue

    areaLR = IntegralByLeftRectangle(   x$, a, b, N)
    areaRR = IntegralByRightRectangle(  x$, a, b, N)
    areaMR = IntegralByMiddleRectangle( x$, a, b, N)
    areaTr = IntegralByTrapezium(       x$, a, b, N)
    areaSi = IntegralBySimpsonRule(     x$, a, b, N)

    print "Left rectangle method   "; using( "##########.##########", areaLR); " diff "; knownValue-areaLR; tab(70); (knownValue-areaLR)/knownValue*100;" %"
    print "Right rectangle method  "; using( "##########.##########", areaRR); " diff "; knownValue-areaRR; tab(70); (knownValue-areaRR)/knownValue*100;" %"
    print "Middle rectangle method "; using( "##########.##########", areaMR); " diff "; knownValue-areaMR; tab(70); (knownValue-areaMR)/knownValue*100;" %"
    print "Trapezium  method       "; using( "##########.##########", areaTr); " diff "; knownValue-areaTr; tab(70); (knownValue-areaTr)/knownValue*100;" %"
    print "Simpson's Rule          "; using( "##########.##########", areaSi); " diff "; knownValue-areaSi; tab(70); (knownValue-areaSi)/knownValue*100;" %"

    print

wend

end

'------------------------------------------------------
    'we have N sizes, that gives us N+1 points
    'point 0 is a
    'point N is b
    'point i is xi =a +i *h
    'Often, precision is (sharper?) then single step area
    'So there should be EXACT number of steps, hence loop by integer i.

function IntegralByLeftRectangle( x$, a, b, N)
    h = ( b -a) /N
    s = 0
    for i = 0 to N -1
        x = a +i *h
        s = s + h *eval( x$)
    next
    IntegralByLeftRectangle = s
end function

function IntegralByRightRectangle( x$, a, b, N)
    h =( b -a) /N
    s = 0
    for i =1 to N
        x = a +i *h
        s = s + h *eval( x$)
    next
    IntegralByRightRectangle = s
end function

function IntegralByMiddleRectangle( x$, a, b, N)
    h =( b -a) /N
    s = 0
    for i =0 to N -1
        x = a +i *h +h /2
        s = s + h *eval( x$)
    next
    IntegralByMiddleRectangle = s
end function

function IntegralByTrapezium( x$, a, b, N)
'Formula is h*((f(a)+f(b))/2 + sum_{i=1}^{N-1} (f(x_i)))
    h  =( b -a) /N
    x  = a
    fa =eval( x$)
    x  =b
    fb =eval( x$)
    s = h *( fa +fb) /2
    for i =1 to N -1
        x = a +i *h
        s = s + h *eval( x$)
    next
    IntegralByTrapezium = s
end function

function IntegralBySimpsonRule( x$, a, b, N)
    'Simpson
    'N should be even.
    if N mod 2 then N =N +1
    'It really doesn't look right to double number of points from N to 2N -
    ' - this method is most accurate of all presented!
    'So we use NN as N/2, and N will be 2NN
    'Formula is h/6*( f(a)+f(b) + 4*(f(x_1)+f(x_3)+...+f(x_{2NN-1})+ 2*(f(x_2)+f(x_4)+...+f(x_{2NN-2})) )
    'Somehow I messed up h/6, h/3 and what is h, regarding "n=number of double intervals of size 2h"
    NN =N /2

    h  =( b -a) /N
    x  =a
    fa =eval (x$)
    x  =b
    fb =eval( x$)
    s = h /3 *( fa +fb)
    for i =1 to 2 *NN -1 step 2
        x = a +i *h
        s = s + h /3 *4 *eval( x$)  'odd points
    next
    for i =2 to 2 *NN -2 step 2
        x = a +i *h
        s = s + h /3 *2 *eval( x$)  'even points
    next

    IntegralBySimpsonRule = s
end function

'=======================================================
data "x^3",  0,    1,     100,          0.25
data "x^-1", 1,  100,    1000,          4.605170
data "x",    0, 5000,    1000,   12500000.0   '   should use 5 000 000 steps
data "x",    0, 6000,    1000,   18000000.0   '   should use 6 000 000 steps
data "end"

end