Linear, cubic and Lagrange function interpolation and extrapolation

FunctionInterpolation

let lagrange xys =
  let n = # xys in
  let xys = Array.sortBy [x, _ → x] xys in
  let lagrange j x =
    for 0 1 (n-1) 1 [y, i →
      if i = j then y else
        let xi, _ = Array.get xys i in
        let xj, _ = Array.get xys j in
        y * (x - xi) / (xj - xi)] in
  [ x →
      for 0 1 (n-1) 0 [y, j →
        let _, yj = Array.get xys j in
        y + yj * lagrange j x ] ]

let search xys x =
  let i, _ =
    while (0, # xys) [i, k → k-i > 1] [i, k →
      let j = (i + k) / 2 in
      let xj, _ = Array.get xys j in
      if x < xj then i, j else j, k] in
  i

let linear =
  [ {} → [x → nan]
  | {_, y} → [_ → y]
  | xys →
      let n = # xys in
      let xys = Array.sortBy [x, _ → x] xys in
      [ x →
          let i = search xys x in
          let i = min (n-2) i in
          let x0, y0 = Array.get xys i in
          let x1, y1 = Array.get xys (i+1) in
          let d = (x - x0) / (x1 - x0) in
          (1 - d) * y0 + d * y1 ] ]

(** Precompute second derivatives for cubic interpolation. *)
let spline xys =
  let get = Array.get in
  let set = Array.Unsafe.set in
  let n = # xys in
  let u = Array.init (n-1) [_ → 0] in
  let d2y = Array.init n [_ → 0] in
  let () =
    for 1 1 (n-2) () [(), i →
      let x0, y0 = get xys (i-1) in
      let x1, y1 = get xys i in
      let x2, y2 = get xys (i+1) in
      let s = (x1 - x0) / (x2 - x0) in
      let p = s * get d2y (i-1) + 2 in
      let () = (s - 1) / p @ set d2y i in
      let () = (y2 - y1) / (x2 - x1) - (y1 - y0) / (x1 - x0) @ set u i in
      (6 * get u i / (x2 - x0) - s * get u (i-1)) / p
      @ set u i ] in
  let () = set d2y (n-1) 0 in
  let () =
    for (n-2) (-1) 0 () [(), k →
      get d2y k * get d2y (k+1) + get u k @ set d2y k ] in
  d2y

let cubic xys =
  let get = Array.get in
  let set = Array.Unsafe.set in
  let xys = Array.sortBy [x, _ → x] xys in
  let d2y = spline xys in
  let n = # xys in
  [ x →
      let i = search xys x in
      let i = min (n-2) i in
      let x0, y0 = get xys i in
      let x1, y1 = get xys (i+1) in
      let h = x1 - x0 in
      let a = (x1 - x) / h in
      let b = (x - x0) / h in
      a*y0 + b*y1 + ((a*a*a - a)*get d2y i +
        (b*b*b - b)*get d2y (i+1))*(h*h)/6 ]