# L-13 MCS 572 Wed 8 Feb 2023 : traprulerecmt.jl
"""
Recursive application of the Trapezoidal rule with tasking.

Run typing
time JULIA_NUM_THREADS=8 julia traprulerecmt.jl 4
at the command prompt for the depth of recursion equal to 4,
with tasks mapped to 8 threads, or alternatively, type
julia -t 8 traprulerecmt.jl 4
at the command prompt in a terminal window.
"""

using Printf
import Base.Threads.@spawn

"""
    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

        t = @spawn rectraprule(level+1,depth,f,a,middle,n)
        return rectraprule(level+1,depth,f,middle,b,n) + fetch(t)
    end
end

"""
Runs with the given depth,
on using exp(x) for x in [0,1].
"""
function test(depth::Int64)
    exact = exp(1.0) - 1.0
    nbr = 10^8
    approx = rectraprule(0,depth,exp,0.0,1.0,nbr)
    err =  abs(exact - approx)
    println(@sprintf("%.16e", exact))
    strerr = @sprintf("%.2e", err)
    print(@sprintf("%.16e",approx))
    println("   error : $strerr")
end

if length(ARGS) > 0
    nbr = parse(Int64, ARGS[1])
    test(nbr)
else
    test(3)
end