Add module Data.Poly.Interpolation

Author: konsumlamm

poly.cabal

@@ -31,6 +31,7 @@ library
   hs-source-dirs: src
   exposed-modules:
     Data.Poly
+    Data.Poly.Interpolation
     Data.Poly.Laurent
     Data.Poly.Semiring
     Data.Poly.Orthogonal
@@ -89,6 +90,7 @@ test-suite poly-tests
     Dense
     DenseLaurent
     DFT
+    Interpolation
     Orthogonal
     Quaternion
     TestUtils

src/Data/Poly/Internal/Dense.hs

@@ -27,6 +27,7 @@ module Data.Poly.Internal.Dense
   , scale
   , pattern X
   , eval
+  , evalk
   , subst
   , deriv
   , integral
@@ -460,6 +461,14 @@ substitute' f (Poly cs) x = fst' $
   G.foldl' (\(acc :*: xn) cn -> acc `plus` f cn xn :*: x `times` xn) (zero :*: one) cs
 {-# INLINE substitute' #-}
 
+-- | Evaluate the kth derivative of the polynomial at a given point.
+evalk :: (G.Vector v a, Num a) => Int -> Poly v a -> a -> a
+evalk k (Poly cs) x = fst' $
+  G.ifoldl' (\(acc :*: xn) i cn -> acc + kth i * cn * xn :*: x * xn) (0 :*: 1) (G.drop k cs)
+  where
+    kth i = fromIntegral $ product [(i + 1)..(i + k)]
+{-# INLINE evalk #-}
+
 -- | Take the derivative of the polynomial.
 --
 -- >>> deriv (X^3 + 3 * X) :: UPoly Int

src/Data/Poly/Interpolation.hs

@@ -0,0 +1,39 @@
+{-# LANGUAGE CPP #-}
+
+module Data.Poly.Interpolation
+  ( lagrange
+  , hermite
+  ) where
+
+#if __GLASGOW_HASKELL__ < 910
+import Data.Foldable (foldl')
+#endif
+
+import Data.Poly.Internal.Dense
+import qualified Data.Vector.Generic as G
+
+-- | Compute the [Lagrange interpolating polynomial](https://en.wikipedia.org/wiki/Lagrange_polynomial).
+--
+-- This is the (unique) polynomial of minimal degree interpolating the given points.
+-- The values are given as @(x, y)@ pairs where @y@ is the value at @x@. The @x@ values must be distinct.
+lagrange :: (G.Vector v a, Eq a, Fractional a) => [(a, a)] -> Poly v a
+lagrange = fst . foldl' f (0, 1)
+  where
+    f (p, w) (x, y) =
+      let a = (y - eval p x) / eval w x
+      in (p + scale 0 a w, scale 1 1 w - scale 0 x w) -- (p + a * w, w * (X - x))
+{-# INLINABLE lagrange #-}
+
+-- | Compute the [Hermite interpolating polynomial](https://en.wikipedia.org/wiki/Hermite_interpolation).
+--
+-- This is the (unique) polynomial of minimal degree interpolating the given points and derivatives.
+-- The values are given as @(x, ys)@ pairs where @ys !! k@ is the k-th derivative at @x@. The @x@ values must be distinct.
+hermite :: (G.Vector v a, Eq a, Fractional a) => [(a, [a])] -> Poly v a
+hermite = fst . foldl' f (0, 1)
+  where
+    f (p, w) (x, ys) = let (_, p', w') = foldl' g (0, p, w) ys in (p', w')
+      where
+        g (k, p', w') y =
+          let a = (y - evalk k p' x) / evalk k w' x
+          in (k + 1, p' + scale 0 a w', scale 1 1 w' - scale 0 x w')
+{-# INLINABLE hermite #-}

test/Interpolation.hs

@@ -0,0 +1,40 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Interpolation (testSuite) where
+
+import Data.Function (on)
+import Data.List (nubBy)
+import Data.Poly hiding (scale)
+import Data.Poly.Interpolation
+import Test.Tasty
+import Test.Tasty.QuickCheck
+
+import TestUtils ()
+
+testSuite :: TestTree
+testSuite = localOption (QuickCheckMaxSize 10) $ testGroup "Interpolation"
+  [ testProperty "lagrange interpolates" $ \xys -> prop_lagrange (nubBy ((==) `on` fst) xys)
+  , testProperty "hermite interpolates"  $ \xys -> prop_hermite (nubBy ((==) `on` (\(x, _, _, _, _) -> x)) xys)
+  , testProperty "lagrange == hermite"   $ \xys -> prop_lagrange_hermite (nubBy ((==) `on` fst) xys)
+  ]
+
+prop_lagrange :: [(Rational, Rational)] -> Property
+prop_lagrange xys =
+  let p = lagrange xys :: VPoly Rational
+  in conjoin $ map (\(x, y) -> eval p x === y) xys
+
+prop_hermite :: [(Rational, Rational, Rational, Rational, Rational)] -> Property
+prop_hermite xys =
+  let
+    p = hermite (map (\(x, y, y', y'', y''') -> (x, [y, y', y'', y'''])) xys) :: VPoly Rational
+    p' = deriv p
+    p'' = deriv p'
+    p''' = deriv p''
+  in conjoin $ map (\(x, y, y', y'', y''') -> eval p x === y .&&. eval p' x === y' .&&. eval p'' x === y'' .&&. eval p''' x === y''') xys
+
+prop_lagrange_hermite :: [(Rational, Rational)] -> Property
+prop_lagrange_hermite xys =
+  let
+    p = lagrange xys :: VPoly Rational
+    q = hermite (map (\(x, y) -> (x, [y])) xys) :: VPoly Rational
+  in p === q

test/Main.hs

@@ -8,6 +8,7 @@ import qualified Dense
 import qualified DenseLaurent
 import qualified DFT
 import qualified Orthogonal
+import qualified Interpolation
 #ifdef SupportSparse
 import qualified Multi
 import qualified MultiLaurent
@@ -20,6 +21,7 @@ main = defaultMain $ testGroup "All"
     [ Dense.testSuite
     , DenseLaurent.testSuite
     , DFT.testSuite
+    , Interpolation.testSuite
     , Orthogonal.testSuite
 #ifdef SupportSparse
     , Sparse.testSuite