# L-10 MCS 572 Wed 18 Sep 2024 : mtcomptrap.jl

# Multithreaded version of the composite trapezoidal rule.
# At the command prompt type for example:
# time JULIA_NUM_THREADS=8 julia mtcomptrap.jl
# to run the code with 8 threads.

using Printf
using Base.Threads

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

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

nt = nthreads()
println("The number of threads : $nt")
subapprox = zeros(nt)

f(x) = sqrt(1 - x^2)
dx = 1/nt
bounds = [i for i=0:dx:1]

# println(traprule(f, 0.0, 1.0, 1000))

timestart = time()
@threads for i=1:nt
    subapprox[i] = traprule(f, bounds[i], bounds[i+1], 1_000_000) 
end
approxpi = 4*sum(subapprox)
elapsed = time() - timestart

println("The approximation for Pi : $approxpi")
err = @sprintf("%.3e", pi - approxpi)
println("with error : $err")
println("The elapsed time : $elapsed seconds")