{-# LANGUAGE ParallelListComp #-}
module Math.Polynomial.Interpolation where
import Math.Polynomial
import Math.Polynomial.Lagrange
import Data.List
polyInterp :: Fractional a => [(a,a)] -> a -> a
polyInterp :: [(a, a)] -> a -> a
polyInterp xys :: [(a, a)]
xys = [a] -> a
forall a. [a] -> a
head ([a] -> a) -> (a -> [a]) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[a]] -> [a]
forall a. [a] -> a
last ([[a]] -> [a]) -> (a -> [[a]]) -> a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(a, a)] -> a -> [[a]]
forall a. Fractional a => [(a, a)] -> a -> [[a]]
neville [(a, a)]
xys
neville :: Fractional a => [(a,a)] -> a -> [[a]]
neville :: [(a, a)] -> a -> [[a]]
neville xys :: [(a, a)]
xys x :: a
x = [[a]]
table
where
(xs :: [a]
xs,ys :: [a]
ys) = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
xys
table :: [[a]]
table = [a]
ys [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
:
[ [ ((a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
x_j) a -> a -> a
forall a. Num a => a -> a -> a
* a
p1 a -> a -> a
forall a. Num a => a -> a -> a
+ (a
x_i a -> a -> a
forall a. Num a => a -> a -> a
- a
x) a -> a -> a
forall a. Num a => a -> a -> a
* a
p0) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
x_i a -> a -> a
forall a. Num a => a -> a -> a
- a
x_j)
| p0 :: a
p0:p1 :: a
p1:_ <- [a] -> [[a]]
forall a. [a] -> [[a]]
tails [a]
row
| a
x_j <- [a]
xs
| a
x_i <- [a]
x_is
]
| [a]
row <- [[a]]
table
| [a]
x_is <- [[a]] -> [[a]]
forall a. [a] -> [a]
tail ([a] -> [[a]]
forall a. [a] -> [[a]]
tails [a]
xs)
, Bool -> Bool
not ([a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [a]
x_is)
]
nevilleDiffs :: Fractional a => [(a,a)] -> a -> [[(a,a)]]
nevilleDiffs :: [(a, a)] -> a -> [[(a, a)]]
nevilleDiffs xys :: [(a, a)]
xys x :: a
x = [[(a, a)]]
table
where
(xs :: [a]
xs,ys :: [a]
ys) = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
xys
table :: [[(a, a)]]
table = [a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
ys [a]
ys [(a, a)] -> [[(a, a)]] -> [[(a, a)]]
forall a. a -> [a] -> [a]
:
[ [ ( (a
x_j a -> a -> a
forall a. Num a => a -> a -> a
- a
x) a -> a -> a
forall a. Num a => a -> a -> a
* (a
c1 a -> a -> a
forall a. Num a => a -> a -> a
- a
d0) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
x_j a -> a -> a
forall a. Num a => a -> a -> a
- a
x_i)
, (a
x_i a -> a -> a
forall a. Num a => a -> a -> a
- a
x) a -> a -> a
forall a. Num a => a -> a -> a
* (a
c1 a -> a -> a
forall a. Num a => a -> a -> a
- a
d0) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
x_j a -> a -> a
forall a. Num a => a -> a -> a
- a
x_i)
)
| (_c0 :: a
_c0,d0 :: a
d0):(c1 :: a
c1,_d1 :: a
_d1):_ <- [(a, a)] -> [[(a, a)]]
forall a. [a] -> [[a]]
tails [(a, a)]
row
| a
x_j <- [a]
xs
| a
x_i <- [a]
x_is
]
| [(a, a)]
row <- [[(a, a)]]
table
| [a]
x_is <- [[a]] -> [[a]]
forall a. [a] -> [a]
tail ([a] -> [[a]]
forall a. [a] -> [[a]]
tails [a]
xs)
, Bool -> Bool
not ([a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [a]
x_is)
]
iterativePolyFit :: (Fractional a, Eq a) => [(a,a)] -> Poly a
iterativePolyFit :: [(a, a)] -> Poly a
iterativePolyFit = Endianness -> [a] -> Poly a
forall a. (Num a, Eq a) => Endianness -> [a] -> Poly a
poly Endianness
LE ([a] -> Poly a) -> ([(a, a)] -> [a]) -> [(a, a)] -> Poly a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(a, a)] -> [a]
forall a. Fractional a => [(a, a)] -> [a]
loop
where
loop :: [(a, a)] -> [a]
loop [] = []
loop xys :: [(a, a)]
xys = a
c0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [(a, a)] -> [a]
loop (Int -> [(a, a)] -> [(a, a)]
forall a. Int -> [a] -> [a]
drop 1 [(a, a)]
xys')
where
c0 :: a
c0 = [(a, a)] -> a -> a
forall a. Fractional a => [(a, a)] -> a -> a
polyInterp [(a, a)]
xys 0
xys' :: [(a, a)]
xys' =
[ (a
x,(a
y a -> a -> a
forall a. Num a => a -> a -> a
- a
c0) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
x)
| (x,y) <- [(a, a)]
xys
]
lagrangePolyFit :: (Fractional a, Eq a) => [(a,a)] -> Poly a
lagrangePolyFit :: [(a, a)] -> Poly a
lagrangePolyFit xys :: [(a, a)]
xys = [Poly a] -> Poly a
forall a. (Num a, Eq a) => [Poly a] -> Poly a
sumPolys
[ a -> Poly a -> Poly a
forall a. (Num a, Eq a) => a -> Poly a -> Poly a
scalePoly a
f ((Poly a, a) -> Poly a
forall a b. (a, b) -> a
fst (Poly a -> a -> (Poly a, a)
forall a. (Num a, Eq a) => Poly a -> a -> (Poly a, a)
contractPoly Poly a
p a
x))
| a
f <- (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) [a]
ys [a]
phis
| a
x <- [a]
xs
]
where
(xs :: [a]
xs,ys :: [a]
ys) = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
xys
p :: Poly a
p = [a] -> Poly a
forall a. (Num a, Eq a) => [a] -> Poly a
lagrange [a]
xs
phis :: [a]
phis = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ((a, a) -> a
forall a b. (a, b) -> b
snd ((a, a) -> a) -> (a -> (a, a)) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Poly a -> a -> (a, a)
forall a. (Num a, Eq a) => Poly a -> a -> (a, a)
evalPolyDeriv Poly a
p) [a]
xs