(*Recursive Trapezium Rule using functional concepts*)
(*Here is demonstrating the construction of a function to approximate the integral of an (R->R) function in [a,b]
  by means of the trapezium rule*)
(*This is done by first taking the mathematical definition, then abstracting to get an ML definition*)
(*Let statements are then used to hold functions that will later on be curried to build the function*)

(*I like this because it uses map and fold for good reasons to complete the task.
  And I learned a bit about local vs let and how to chain results together efficiently*)

(*Mathematical def -> ML def*)
trap f (a,b) n
= h/2*(f(a) + 2f(a+h) + 2f(a+2*h) + .. + f(b))
= h/2 * (foldr (op +) [f(a), 2f(a+h), 2f(a+2h) .. , f(b)])
= h/2 * (foldr (op +) f(a)::(map (fn x => 2*x) [f(a+h), f(a+2h), .. , f(a+(n-1)h)])::f(b)::[]
= h/2 * (foldr (op +) f(a)::(map (fn x => 2*x) (map f [a+h, a+2h, ..., a+(n-1)h])::f(b)::[]
= h/2 * (foldr (op +) f(a)::(map (fn x => 2*x) (map f (map (fn i => a+i*h) [1, 2, .., (n-1)])::f(b)::[]
= h/2 * (foldr (op +) f(a)::(map (fn x => 2*x) (map f (map (fn i => a+i*h) (list_gen (1,n-1))))::f(b)::[]

(*Building of function*)
fun list_gen (a, b) = if (a = b) then b::[]
              else a::list_gen ((a+1), b);

fun trap f (a,b) n =
let val h = (b-a)/real(n)
    val quad_indexes = list_gen
    val quad_points = map (fn i => a+real(i)*h)
    val quad_vals = map f
    val double = map (fn x => 2.0*x)
    val sum = foldr (op+) 0.0
in
    (h/2.0)*(sum ((f a)::(double (quad_vals (quad_points (list_gen (1, n-1)))))@[f b]))
end;


(*This could possibly be done in one handwritten recursive function, let's try*)
(*Attempt at "all in one" definition*)
fun trap f (a,b) n =
    let fun trap_helper (a,b) h =
        if (a=b) then (f a) + (f b)
        else 2.0*(f a) + (trap_helper f (a+h, b) h)
    in
        ((b-a)/(2.0*n)) * (trap_helper (a,b) ((b-a)/(2.0*n)))
    end

(*This does not work as reals are not equality types and therefore cannot be compared*)