PROGRAM integ
! compute integral of f(x) from x = a to x = b
! using rectangular approximation
CALL initial(a,b,sum,dx,delta,n)
DO
   CALL sumf(a,b,sum,dx,delta)
   CALL output(sum,dx,n)
   CALL double(dx,delta,n)
LOOP until key input
GET KEY k
END

SUB initial(a,b,sum,dx,delta,n)
    DECLARE DEF f
    INPUT prompt "lower limit a = ": a
    INPUT prompt "upper limit b = ": b
    LET n = 2                     ! initial number of intervals
    LET dx = (b - a)/n
    LET delta = dx                ! delta = dx only for n = 2
    LET sum = f(a)   ! -> (f(a) + f(b))/2 for trapezoidal rule
    PRINT
    PRINT "    n","rectangular approximation"
    PRINT
END SUB

SUB sumf(a,b,sum,dx,delta)
    DECLARE DEF f
    LET x = a + dx                ! dx is integration interval
    DO while x < b
       LET sum = sum + f(x)
       LET x = x + delta
    LOOP
END SUB

SUB double(dx,delta,n)
    ! double number of intervals
    LET delta = dx
    LET n = 2*n
    LET dx = 0.5*dx               ! dx does not equal delta for n > 2
END SUB

SUB output(sum,dx,n)
    PRINT using "#######": n;
    PRINT,
    PRINT using "####.#######": sum*dx
END SUB

DEF f(x) = cos(x)