{-# LANGUAGE ParallelListComp #-}
{-# OPTIONS_GHC -O3 #-}
import Data.List (zip5)
import Graphics.Gnuplot.Simple
import Numeric.FastMath
import Text.Printf (printf)
f_sum :: Num a => [(t -> a)] -> t -> a
f_sum xs = \x -> sum [f x | f <- xs]
f_product :: Num a => [(t -> a)] -> t -> a
f_product xs = \x -> product [f x | f <- xs]
delete :: Int -> [a] -> [a]
delete k xs = let (ys, zs) = splitAt k xs in ys ++ tail zs
fact :: (Enum a, Num a) => a -> a
fact n = product [1..n]
lagrange :: Fractional a => [a] -> [a] -> a -> a
lagrange xs ys
| n /= length ys = error "Incorrect input size"
| otherwise = f_sum [\x -> f x * (y / a) | f <- l1 | a <- l2 | y <- ys]
where n = length xs
l1 = map f_product [[\x -> x - a | a <- delete k xs] | k <- [0..n-1]]
l2 = map product [let xk = xs !! k in [xk - a | a <- delete k xs] | k <- [0..n-1]]
step :: Fractional a => a -> a -> Int -> a
step a b n = abs (b - a) / (fromIntegral $ n - 1)
memo2 :: (Enum a, Num a) => (a -> a -> b) -> [[b]]
memo2 f = map (\x -> map (f x) [0..]) [0..]
newtonK :: Fractional a => Int -> a -> a -> [a] -> a -> a
newtonK k x0 h deltas x = q k x / (fromIntegral . fact $ k) * (deltas !! k)
where q k = f_product [\x -> (x - x0) / h - (fromIntegral a) | a <- [0..k-1]]
newton :: Fractional a => a -> a -> [a] -> a -> a
newton a b ys = f_sum [newtonK k a h deltas | k <- [0..n-1]]
where h = step a b n
n = length ys
deltas = map (\(a, b) -> delta' a b) [(i, 0) | i <- [0..n-1]]
delta 0 n = ys !! n
delta 1 n = ys !! (n + 1) - ys !! n
delta m n = delta' (m - 1) (n + 1) - delta' (m - 1) n
deltastore = memo2 delta
delta' i j = deltastore !! i !! j
splitting :: (Enum t, Fractional t) => t -> t -> Int -> [t]
splitting a b n = [a, a + h .. b]
where h = step a b n
optimalSplitting :: Floating t => t -> t -> Int -> [t]
optimalSplitting a b n = [0.5 * (b - a) * c i + (b + a) | i <- [n,n-1..0]]
where c i = cos $ pi * (fromIntegral $ 2 * i + 1) / (fromIntegral $ 2 * n + 2)
printTable :: [Double] -> [Double] -> [Double] -> [Double] -> [Double] -> IO ()
printTable s1 s2 s3 s4 s5 = sequence_ [printRow s | s <- (zip5 s1 s2 s3 s4 s5)]
where printRow (c1, c2, c3, c4, c5) = printf "%13.7f | %13.10f | %13.10f | %13.10f | %13.10f\n" c1 c2 c3 c4 c5
main :: IO ()
main = sequence_ $
plotFuncs options xs [f1, f2, g] :
printf "%13s | %13s | %13s | %13s | %13s\n" "x" "y" "lagrange" "lagrange opt" "newton" :
printTable xs1 ys1 fs1 fs2 gs :
[]
where a = 1
b = 11
n = 21
h = \x -> 10 `logBase` x + 7 / (2 * x + 6)
xs = splitting a b (n * 10)
xs1 = splitting a b n
xs2 = optimalSplitting a b n
ys1 = map h xs1
ys2 = map h xs2
fs1 = map f1 xs1
fs2 = map f2 xs1
gs = map g xs1
f1 = lagrange xs1 ys1 :: Double -> Double
f2 = lagrange xs2 ys2 :: Double -> Double
g = newton a b ys1 :: Double -> Double
options = [Title "Interpolation", XRange (a - 1, b + 1)]