# L-13 MCS 572 Wed 8 Feb 2023 : traprulerec.jl

# recursive application of the Trapezoidal rule

using Printf

"""
    traprule(f::Function,a::Float64,b::Float64,n::Int64)

Applies the composite Trapezoidal rule to approximate the integral
of f over [a,b] with n+1 function evaluations.
"""
function traprule(f::Function,a::Float64,b::Float64,n::Int64)
    h = (b-a)/n
    y = (f(a) + f(b))/2
    x = a
    for i=1:n-1
        x = x + h
        y = y + f(x)
    end
    return h*y
end

"""
    rectraprule(level::Int64,depth::Int64,
                f::Function,a::Float64,b::Float64,n::Int64)

Recursive application of the composite Trapezoidal rule to f on [a,b]
with n+1 function evaluations, at the leaves of the recursion,
when level equals the given depth.
"""
function rectraprule(level::Int64,depth::Int64,
                     f::Function,a::Float64,b::Float64,n::Int64)
    if level == depth
        return traprule(f,a,b,n)
    else
        middle = (b-a)/2

        return (rectraprule(level+1,depth,f,a,middle,n)
              + rectraprule(level+1,depth,f,middle,b,n))
    end
end

"""
Does a basic check, using exp(x) for x in [0,1].
"""
function test()
    exact = exp(1.0) - 1.0
    nbr = 10^8
    approx1 = traprule(exp,0.0,1.0,nbr)
    approx2 = rectraprule(0,3,exp,0.0,1.0,nbr)
    err1 =  abs(exact - approx1)
    strerr1 = @sprintf("%.2e", err1)
    err2 =  abs(exact - approx2)
    strerr2 = @sprintf("%.2e", err2)
    println(@sprintf("%.16e", exact))
    print(@sprintf("%.16e",approx1))
    println("   error : $strerr1")
    print(@sprintf("%.16e",approx2))
    println("   error : $strerr2")
end

test()