From 686a047d5ac1d5ef2fddb23496d52ed1baffd1bb Mon Sep 17 00:00:00 2001
From: konsumlamm <konsumlamm@gmail.com>
Date: Mon, 7 Apr 2025 14:31:32 +0200
Subject: [PATCH 1/3] Add module `Data.Poly.Interpolation`

---
 poly.cabal                      |  2 ++
 src/Data/Poly/Internal/Dense.hs |  9 ++++++++
 src/Data/Poly/Interpolation.hs  | 39 ++++++++++++++++++++++++++++++++
 test/Interpolation.hs           | 40 +++++++++++++++++++++++++++++++++
 test/Main.hs                    |  2 ++
 5 files changed, 92 insertions(+)
 create mode 100644 src/Data/Poly/Interpolation.hs
 create mode 100644 test/Interpolation.hs

diff --git a/poly.cabal b/poly.cabal
index 126530b..aeb6326 100644
--- a/poly.cabal
+++ b/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
diff --git a/src/Data/Poly/Internal/Dense.hs b/src/Data/Poly/Internal/Dense.hs
index 04e0636..3f6f8ab 100644
--- a/src/Data/Poly/Internal/Dense.hs
+++ b/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
diff --git a/src/Data/Poly/Interpolation.hs b/src/Data/Poly/Interpolation.hs
new file mode 100644
index 0000000..2b4f257
--- /dev/null
+++ b/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 #-}
diff --git a/test/Interpolation.hs b/test/Interpolation.hs
new file mode 100644
index 0000000..338fee1
--- /dev/null
+++ b/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
diff --git a/test/Main.hs b/test/Main.hs
index 4b45cfd..af756d9 100644
--- a/test/Main.hs
+++ b/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

From 44679ab440919454dea769a915c7ab2e9874010b Mon Sep 17 00:00:00 2001
From: konsumlamm <konsumlamm@gmail.com>
Date: Tue, 8 Apr 2025 12:15:36 +0200
Subject: [PATCH 2/3] Avoid CPP

---
 src/Data/Poly/Interpolation.hs | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/src/Data/Poly/Interpolation.hs b/src/Data/Poly/Interpolation.hs
index 2b4f257..83dda90 100644
--- a/src/Data/Poly/Interpolation.hs
+++ b/src/Data/Poly/Interpolation.hs
@@ -1,13 +1,10 @@
-{-# LANGUAGE CPP #-}
-
 module Data.Poly.Interpolation
   ( lagrange
   , hermite
   ) where
 
-#if __GLASGOW_HASKELL__ < 910
-import Data.Foldable (foldl')
-#endif
+import Prelude hiding (Foldable(..))
+import Data.Foldable
 
 import Data.Poly.Internal.Dense
 import qualified Data.Vector.Generic as G

From 0bea925d379d5c055568c9277c51267d4271434a Mon Sep 17 00:00:00 2001
From: konsumlamm <konsumlamm@gmail.com>
Date: Tue, 8 Apr 2025 16:39:11 +0200
Subject: [PATCH 3/3] Use `Map` to represent points

---
 poly.cabal                     |  2 ++
 src/Data/Poly/Interpolation.hs | 20 +++++++++++---------
 test/Interpolation.hs          | 23 +++++++++++------------
 3 files changed, 24 insertions(+), 21 deletions(-)

diff --git a/poly.cabal b/poly.cabal
index aeb6326..b9e5bb8 100644
--- a/poly.cabal
+++ b/poly.cabal
@@ -65,6 +65,7 @@ library
 
   build-depends:
     base >= 4.12 && < 5,
+    containers >= 0.5.4,
     deepseq >= 1.1 && < 1.6,
     primitive >= 0.6,
     semirings >= 0.5.2,
@@ -102,6 +103,7 @@ test-suite poly-tests
       SparseLaurent
   build-depends:
     base >=4.10 && <5,
+    containers,
     mod >=0.1.2,
     poly,
     QuickCheck >=2.12 && <2.14.3,
diff --git a/src/Data/Poly/Interpolation.hs b/src/Data/Poly/Interpolation.hs
index 83dda90..548fd30 100644
--- a/src/Data/Poly/Interpolation.hs
+++ b/src/Data/Poly/Interpolation.hs
@@ -4,7 +4,9 @@ module Data.Poly.Interpolation
   ) where
 
 import Prelude hiding (Foldable(..))
-import Data.Foldable
+import qualified Data.Foldable as F
+
+import Data.Map
 
 import Data.Poly.Internal.Dense
 import qualified Data.Vector.Generic as G
@@ -12,11 +14,11 @@ 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)
+-- The keys are the @x@ values and the associated @y@ is the value at @x@.
+lagrange :: (G.Vector v a, Eq a, Fractional a) => Map a a -> Poly v a
+lagrange = fst . foldlWithKey' f (0, 1)
   where
-    f (p, w) (x, y) =
+    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 #-}
@@ -24,11 +26,11 @@ lagrange = fst . foldl' f (0, 1)
 -- | 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)
+-- The keys are the @x@ values and the associated @ys@ are the values and derivatives at @x@, where @ys !! k@ is the k-th derivative.
+hermite :: (G.Vector v a, Eq a, Fractional a) => Map a [a] -> Poly v a
+hermite = fst . foldlWithKey' f (0, 1)
   where
-    f (p, w) (x, ys) = let (_, p', w') = foldl' g (0, p, w) ys in (p', w')
+    f (p, w) x ys = let (_, p', w') = F.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
diff --git a/test/Interpolation.hs b/test/Interpolation.hs
index 338fee1..967d32a 100644
--- a/test/Interpolation.hs
+++ b/test/Interpolation.hs
@@ -2,8 +2,7 @@
 
 module Interpolation (testSuite) where
 
-import Data.Function (on)
-import Data.List (nubBy)
+import Data.Map
 import Data.Poly hiding (scale)
 import Data.Poly.Interpolation
 import Test.Tasty
@@ -13,28 +12,28 @@ 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)
+  [ testProperty "lagrange interpolates" prop_lagrange
+  , testProperty "hermite interpolates"  prop_hermite
+  , testProperty "lagrange == hermite"   prop_lagrange_hermite
   ]
 
-prop_lagrange :: [(Rational, Rational)] -> Property
+prop_lagrange :: Map Rational Rational -> Property
 prop_lagrange xys =
   let p = lagrange xys :: VPoly Rational
-  in conjoin $ map (\(x, y) -> eval p x === y) xys
+  in conjoin $ fmap (\(x, y) -> eval p x === y) (toList xys)
 
-prop_hermite :: [(Rational, Rational, Rational, Rational, Rational)] -> Property
+prop_hermite :: Map 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 = hermite (fmap (\(y, y', y'', y''') -> [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
+  in conjoin $ fmap (\(x, (y, y', y'', y''')) -> eval p x === y .&&. eval p' x === y' .&&. eval p'' x === y'' .&&. eval p''' x === y''') (toList xys)
 
-prop_lagrange_hermite :: [(Rational, Rational)] -> Property
+prop_lagrange_hermite :: Map Rational Rational -> Property
 prop_lagrange_hermite xys =
   let
     p = lagrange xys :: VPoly Rational
-    q = hermite (map (\(x, y) -> (x, [y])) xys) :: VPoly Rational
+    q = hermite (fmap (\y -> [y]) xys) :: VPoly Rational
   in p === q