-- | Lagrange interpolation
--
-- See <https://en.wikipedia.org/wiki/Lagrange_polynomial>

{-# LANGUAGE BangPatterns, DataKinds, KindSignatures, GeneralizedNewtypeDeriving, TypeFamilies #-}
module Math.Algebra.Polynomial.Univariate.Lagrange where

--------------------------------------------------------------------------------

import GHC.TypeLits

{-
import Data.Array ( Array , (!) , listArray , assocs ) 
import Data.List

import Data.Proxy

import Math.Algebra.Polynomial.Misc
import Math.Algebra.Polynomial.Pretty
-}

import Math.Algebra.Polynomial.Class

import qualified Math.Algebra.Polynomial.FreeModule as ZMod
import Math.Algebra.Polynomial.FreeModule ( FreeMod , ZMod , QMod )

import Math.Algebra.Polynomial.Univariate

--------------------------------------------------------------------------------
-- * Lagrange interpolation

-- | The Lagrange interpolation for a list of @(x_i,y_i)@ pairs: the resulting polynomial P
-- is the one with minimal degree for which @P(x_i) = y_i@
lagrangeInterp :: KnownSymbol var => [(Rational,Rational)] -> Univariate Rational var 
lagrangeInterp :: [(Rational, Rational)] -> Univariate Rational var
lagrangeInterp [(Rational, Rational)]
xys = Univariate Rational var
final where
  final :: Univariate Rational var
final = [Univariate Rational var] -> Univariate Rational var
forall p. AlmostPolynomial p => [p] -> p
sumP [ CoeffP (Univariate Rational var)
-> Univariate Rational var -> Univariate Rational var
forall p. AlmostPolynomial p => CoeffP p -> p -> p
scaleP ([Rational]
ys[Rational] -> Int -> Rational
forall a. [a] -> Int -> a
!!Int
j) ([Rational] -> Int -> Univariate Rational var
forall (var :: Symbol).
KnownSymbol var =>
[Rational] -> Int -> Univariate Rational var
lagrangePoly [Rational]
xs Int
j) | Int
j<-[Int
0..Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1] ] where
  m :: Int
m = [(Rational, Rational)] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Rational, Rational)]
xys
  ([Rational]
xs,[Rational]
ys) = [(Rational, Rational)] -> ([Rational], [Rational])
forall a b. [(a, b)] -> ([a], [b])
unzip [(Rational, Rational)]
xys

-- | The minimal degree polynomial which evaluates to zero at the given inputs, except a single one
-- on which it evaluates to 1
lagrangePoly :: KnownSymbol var => [Rational] -> Int -> Univariate Rational var
lagrangePoly :: [Rational] -> Int -> Univariate Rational var
lagrangePoly [Rational]
xs Int
j = FreeMod Rational (U var) -> Univariate Rational var
forall coeff (var :: Symbol).
FreeMod coeff (U var) -> Univariate coeff var
Uni (FreeMod Rational (U var) -> Univariate Rational var)
-> FreeMod Rational (U var) -> Univariate Rational var
forall a b. (a -> b) -> a -> b
$ Rational -> FreeMod Rational (U var) -> FreeMod Rational (U var)
forall b c. (Ord b, Eq c, Num c) => c -> FreeMod c b -> FreeMod c b
ZMod.scale (Rational
1Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
denom) FreeMod Rational (U var)
numer where
  numer :: FreeMod Rational (U var)
numer  = [FreeMod Rational (U var)] -> FreeMod Rational (U var)
forall b c.
(Ord b, Monoid b, Eq c, Num c) =>
[FreeMod c b] -> FreeMod c b
ZMod.product [ Int -> FreeMod Rational (U var)
term Int
i    | Int
i<-[Int
0..Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1] , Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
j ]
  denom :: Rational
denom  = [Rational] -> Rational
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product      [ Int -> Rational
x Int
j Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Int -> Rational
x Int
i | Int
i<-[Int
0..Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1] , Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
j ]
  m :: Int
m      = [Rational] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Rational]
xs
  x :: Int -> Rational
x Int
i    = [Rational]
xs [Rational] -> Int -> Rational
forall a. [a] -> Int -> a
!! Int
i
  term :: Int -> FreeMod Rational (U var)
term Int
i = [(U var, Rational)] -> FreeMod Rational (U var)
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList 
    [ (Int -> U var
forall (var :: Symbol). Int -> U var
U Int
1 ,     Rational
1 )
    , (Int -> U var
forall (var :: Symbol). Int -> U var
U Int
0 , - Int -> Rational
x Int
i )
    ]

--------------------------------------------------------------------------------