(* @param f float -> float
 * @param a, b a < b
 * @return result
 * @example
 * let pi = 3.1415926535;;
 * integral sin 0.0 pi;;
 *)
let integral f a b =
    let n = 1e5 in
    let delta = (b -. a) /. n in
    let x i = a +. float_of_int(i) *. delta in
    let y i = f (x i) in
    let delta_area i =
        (y (i - 1) +. y i) *. delta /. 2.0
    in
    let rec loop i acc =
        if x i > b then acc
        else loop (i + 1) (delta_area i +. acc)
    in loop 1 0.0