Psst.. new poll here.

[email protected] web/email now available. Want one? Go here.

Cannot use outlook/hotmail/live here to register as they blocking our mail servers. #microsoftdeez

Obey the Epel!

## Paste

Pasted as Haskell by Minoru ( 14 years ago )

```
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"
```

## Revise this Paste