{-# LANGUAGE ParallelListComp #-}
module Math.Polynomial.Interpolation where

import Math.Polynomial
import Math.Polynomial.Lagrange
import Data.List

-- |Evaluate a polynomial passing through the specified set of points.  The
-- order of the interpolating polynomial will (at most) be one less than 
-- the number of points given.
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

-- |Computes the tableau generated by Neville's algorithm.  Each successive
-- row of the table is a list of interpolants one order higher than the previous,
-- using a range of input points starting at the same position in the input
-- list as the interpolant's position in the output list.
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)
            ]

-- |Computes the tableau generated by a modified form of Neville's algorithm
-- described in Numerical Recipes, Ch. 3, Sec. 2, which records the differences
-- between interpolants at each level.  Each pair (c,d) is the amount to add
-- to the previous level's interpolant at either the same or the subsequent
-- position (respectively) in order to obtain the new level's interpolant.
-- Mathematically, either sum yields the same value, but due to numerical
-- errors they may differ slightly, and some \"paths\" through the table
-- may yield more accurate final results than others.
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]
:
            [ [ ( {-c-} (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)
                , {-d-} (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)
            ]

-- |Fit a polynomial to a set of points by iteratively evaluating the 
-- interpolated polynomial (using 'polyInterp') at 0 to establish the
-- constant coefficient and reducing the polynomial by subtracting that
-- coefficient from all y's and dividing by their corresponding x's.
-- 
-- Slower than 'lagrangePolyFit' but stable under different sets of 
-- conditions.
-- 
-- Note that computing the coefficients of a fitting polynomial is an 
-- inherently ill-conditioned problem.  In most cases it is both faster and 
-- more accurate to use 'polyInterp' or 'nevilleDiffs' instead of evaluating
-- a fitted polynomial.
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
                    ]

-- |Fit a polynomial to a set of points using barycentric Lagrange polynomials.
-- 
-- Note that computing the coefficients of a fitting polynomial is an 
-- inherently ill-conditioned problem.  In most cases it is both faster and 
-- more accurate to use 'polyInterp' or 'nevilleDiffs' instead of evaluating
-- a fitted polynomial.
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