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