-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Data.Polynomial.Interpolation.Lagrange
-- Copyright   :  (c) Masahiro Sakai 2012
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  portable
--
-- References:
--
-- * Lagrange polynomial <http://en.wikipedia.org/wiki/Lagrange_polynomial>
--
-----------------------------------------------------------------------------
module ToySolver.Data.Polynomial.Interpolation.Lagrange
  ( interpolate
  ) where

import ToySolver.Data.Polynomial (UPolynomial, X (..))
import qualified ToySolver.Data.Polynomial as P

interpolate :: (Eq k, Fractional k) => [(k,k)] -> UPolynomial k
interpolate :: [(k, k)] -> UPolynomial k
interpolate [(k, k)]
zs = [UPolynomial k] -> UPolynomial k
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([UPolynomial k] -> UPolynomial k)
-> [UPolynomial k] -> UPolynomial k
forall a b. (a -> b) -> a -> b
$ do
  (k
xj,k
yj) <- [(k, k)]
zs
  let lj :: Polynomial k v -> Polynomial k v
lj Polynomial k v
x = [Polynomial k v] -> Polynomial k v
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [k -> Polynomial k v
forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant (k
1 k -> k -> k
forall a. Fractional a => a -> a -> a
/ (k
xj k -> k -> k
forall a. Num a => a -> a -> a
- k
xm)) Polynomial k v -> Polynomial k v -> Polynomial k v
forall a. Num a => a -> a -> a
* (Polynomial k v
x Polynomial k v -> Polynomial k v -> Polynomial k v
forall a. Num a => a -> a -> a
- k -> Polynomial k v
forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant k
xm) | (k
xm,k
_) <- [(k, k)]
zs, k
xj k -> k -> Bool
forall a. Eq a => a -> a -> Bool
/= k
xm]
  let x :: UPolynomial k
x = X -> UPolynomial k
forall a v. Var a v => v -> a
P.var X
X
  UPolynomial k -> [UPolynomial k]
forall (m :: * -> *) a. Monad m => a -> m a
return (UPolynomial k -> [UPolynomial k])
-> UPolynomial k -> [UPolynomial k]
forall a b. (a -> b) -> a -> b
$ k -> UPolynomial k
forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant k
yj UPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
* UPolynomial k -> UPolynomial k
forall v. Ord v => Polynomial k v -> Polynomial k v
lj UPolynomial k
x