import "Prelude"
include "Vector"


@private
@inline
chooseInterval :: Vector Double -> Double -> (<e> a) -> (Integer -> <e> a) -> (<e> a) -> <e> a
chooseInterval table x down inRange up
    | x < table!0        = down
    | x >= table!(len-1) = up
    | otherwise          = search 0 (len-1)
  where
    len = length table
    search low high | high <= low+1 = inRange low
                    | otherwise = do
       middle = (low + high) `div` 2
       if x < table!middle
       then search low middle
       else search middle high


interpolateNearest :: Vector Double -> Vector Double -> Double -> Double
interpolateNearest xs ys x = chooseInterval xs x
    (ys!0)
    (\i -> if x - xs!i < xs!(i+1) - x then ys!i else ys!(i+1))
    (ys!(length xs-1))
    
interpolateLinear :: Vector Double -> Vector Double -> Double -> Double
interpolateLinear xs ys x = chooseInterval xs x
    (ys!0)
    (\i -> lin (xs!i) (xs!(i+1)) (ys!i) (ys!(i+1)))
    (ys!(length ys-1))
  where
    lin x0 x1 y0 y1 = y0 + (y1-y0) * (x-x0) / (x1-x0)

interpolateCubic :: Vector Double -> Vector Double -> Double -> Double
interpolateCubic xs ys x = chooseInterval xs x
    (ys!0)
    (\i -> let a = max 0 (i-1) ; b = min (i+2) (length xs-1) in
        cubic (xs!a) (xs!i) (xs!(i+1)) (xs!b)
              (ys!a) (ys!i) (ys!(i+1)) (ys!b))
    (ys!(length ys-1))
  where
    cubic x0 x1 x2 x3 y0 y1 y2 y3 = let
        l = x2-x1
        dx1 = (y2-y0) * l / (x2-x0) 
        dx2 = (y3-y1) * l / (x3-x1)
        t = (x-x1) / l
      in 
        y1 + t * (dx1 + t * ((-3*y1 + 3*y2 - 2*dx1 - dx2) + t * (2*y1 - 2*y2 + dx1 + dx2)))