import Data.List (intersperse, zip4, zipWith5, zipWith6)
import Text.Printf (printf)

main :: IO()
main = do
  let points = [-5.0, -3.0, -1.0, 1.0, 3.0]
      ys = map f points

  putStrLn $ unlines $ (interpolate points ys) ++ ["          : 0"]

-- function defined by a variant
f :: Double -> Double
f x = (sin x) + cubeRoot (2*x)

cubeRoot :: Double -> Double
cubeRoot x = signum x * (abs x) ** (1/3)

interpolate :: [Double] -> [Double] -> [String]
interpolate points values = zipWith splineToString [0..(length $ tail values)] $
    computeSplines $ back $ forth $ prepareCoefficients points values
  where
    -- h is a length between two successive points
    hs = zipWith (-) (tail points) points

    -- in:  list of x and y coordinates of knots
    -- out: list of coefficients for three-diagonal matrix (a, c, b, f)
    prepareCoefficients :: [Double] -> [Double] -> [(Double, Double, Double, Double)]
    prepareCoefficients _ ys = zip4 as cs bs fs
      where
        as = 0 : hs
        cs = map (2*) $ zipWith (+) (tail hs) hs -- c = 2(h_i-1 + h_i)
        bs = (tail hs) ++ [0] -- bs starts from second element of hs
        fs = zipWith6 (\a b c d e g -> 6*(((a-b)/c) - (d-e)/g))
                      (drop 2 ys)
                      (tail ys)
                      (tail hs)
                      (tail ys)
                      ys
                      hs

    -- in:  list of (a, c, b, f)
    -- out: list of (alpha, beta) (in reverse order)
    --      first pair in out would be (0, x), where x is x_n - part of a solution
    forth :: [(Double, Double, Double, Double)] -> [(Double, Double)]
    forth coeffs = compAlphaBeta coeffs [(0, 0)]
      where
        compAlphaBeta :: [(Double, Double, Double, Double)] -> [(Double, Double)] -> [(Double, Double)]
        compAlphaBeta [] ac = ac
        compAlphaBeta ((a, c, b, f):cs) ac@((x, y):_) = compAlphaBeta cs ((alpha, beta) : ac)
          where
            -- x and y are alpha and beta for i-1
            alpha = -b/(a * x + c)
            beta  = (f - a * y)/(a * x + c)

    -- in:  list of (alpha, beta) (in reverse order)
    --      first pair in this list would be (_, x), where x is x_n; we should just extract it
    -- out: list of c
    back :: [(Double, Double)] -> [Double]
    back ((_, x_nth):pairs) = compResults pairs [x_nth]
      where
        compResults :: [(Double, Double)] -> [Double] -> [Double]
        compResults [] ac = 0 : ac ++ [0]
        compResults ((a, b):ps) ac@(x_prev:_) = compResults ps ((a * x_prev + b) : ac)

    -- in:  list of c
    -- out: list of (a, b, c, d) -- spline coefficients
    computeSplines :: [Double] -> [(Double, Double, Double, Double)]
    computeSplines s = zip4 as bs cs ds
      where
         as = zipWith3 (\x y z -> (x-y)/(6*z))
                (tail s)
                s
                hs
         bs = map (/2) s
         cs = zipWith5 (\a b c d e -> ((a-b)/c) - ((2*c*d + c*e)/6))
                (tail values)
                values
                hs
                s
                (tail s)
         ds = values

    -- in:  i and (a, b, c, d)
    -- out: string representation of given spline function
    splineToString :: Int -> (Double, Double, Double, Double) -> String
    splineToString i (a, b, c, d)
      | i == 0    = "spline(x) = " ++ eq ++ " \\"
      | otherwise = "          : " ++ eq ++ " \\"
      where
        eq = "x >= " ++ printf "%.0f" (points !! i) ++ " && " ++
             "x < "  ++ printf "%.0f" (points !! (i+1)) ++ " ? "
             ++ smartPrintDouble a ++ "*(x" ++ x_i ++ ")**3"
             ++ smartPrintDouble b ++ "*(x" ++ x_i ++ ")**2"
             ++ smartPrintDouble c ++ "*(x" ++ x_i ++ ")"
             ++ printf "%+.6f" d
        x_i = smartPrintInt $ -(points !! (i+1))


smartPrint :: String -> Double -> String
smartPrint format x
  | x > 0     = " + " ++ printf format x
  | otherwise = " - " ++ (printf format $ abs x)

smartPrintDouble, smartPrintInt :: Double -> String
smartPrintDouble = smartPrint "%.6f"
smartPrintInt = smartPrint "%.0f"