add ToySolver.Data.Polynomial.Interpolation.Hermite

src/ToySolver/Data/Polynomial/Interpolation/Hermite.hs

@@ -0,0 +1,78 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  ToySolver.Data.Polynomial.Interpolation.Hermite
+-- Copyright   :  (c) Masahiro Sakai 2020
+-- License     :  BSD-style
+--
+-- Maintainer  :  [email protected]
+-- Stability   :  provisional
+-- Portability :  non-portable
+--
+-- References:
+--
+-- * Lagrange polynomial <https://en.wikipedia.org/wiki/Hermite_interpolation>
+--
+-----------------------------------------------------------------------------
+module ToySolver.Data.Polynomial.Interpolation.Hermite
+  ( interpolate
+  ) where
+
+import Data.List (unfoldr)
+
+import ToySolver.Data.Polynomial (UPolynomial, X (..))
+import qualified ToySolver.Data.Polynomial as P
+
+
+-- | Compute a hermite Hermite interpolation from a list
+-- \([(x_0, [y_0, y'_0, \ldots y^{(m)}_0]), (x_1, [y_1, y'_1, \ldots y^{(m)}_1]), \ldots]\).
+interpolate :: (Eq k, Fractional k) => [(k,[k])] -> UPolynomial k
+interpolate xyss = sum $ zipWith (\c p -> P.constant c * p) cs ps
+  where
+    x = P.var X
+    ps = scanl (*) (P.constant 1) [(x - P.constant x') | (x', ys') <- xyss, _ <- ys']
+    cs = map head $ dividedDifferenceTable xyss
+
+
+type D a = Either (a, Int, [a]) ((a, a), a)
+
+
+dividedDifferenceTable :: forall a. Fractional a => [(a, [a])] -> [[a]]
+dividedDifferenceTable xyss = unfoldr f xyss'
+  where
+    xyss' :: [D a]
+    xyss' =
+      [ Left (x, len, zipWith (\num den -> num / fromInteger den) ys (scanl (*) 1 [1..]))
+      | (x, ys) <- xyss
+      , let len = length ys
+      ]
+
+    d :: D a -> [a]
+    d (Left (_x, n, ys)) = replicate n (head ys)
+    d (Right (_, y)) = [y]
+
+    gx1, gx2, gy :: D a -> a
+    gx1 (Left (x, _n, _ys)) = x
+    gx1 (Right ((x1, _x2), _y)) = x1
+    gx2 (Left (x, _n, _ys)) = x
+    gx2 (Right ((_x1, x2), _y)) = x2
+    gy (Left (_x, _n, ys)) = head ys
+    gy (Right (_, y)) = y
+
+    f :: [D a] -> Maybe ([a], [D a])
+    f [] = Nothing
+    f xs = Just (concatMap d xs, h xs)
+      where
+        h :: [D a] -> [D a]
+        h (z1 : zzs) =
+          (case z1 of
+             Left (x, n, _ : ys@(_ : _)) -> [Left (x, n-1, ys)]
+             _ -> [])
+          ++
+          (case zzs of
+             z2 : _ -> [Right ((gx1 z1, gx2 z2), (gy z2 - gy z1) / (gx2 z2 - gx1 z1))]
+             [] -> [])
+          ++
+          h zzs
+        h [] = []

test/TestPolynomial.hs

@@ -23,6 +23,7 @@ import qualified ToySolver.Data.Polynomial.GroebnerBasis as GB
 import qualified ToySolver.Data.Polynomial.Factorization.FiniteField as FactorFF
 import qualified ToySolver.Data.Polynomial.Factorization.Hensel.Internal as Hensel
 import qualified ToySolver.Data.Polynomial.Factorization.Zassenhaus as Zassenhaus
+import qualified ToySolver.Data.Polynomial.Interpolation.Hermite as HermiteInterpolation
 import qualified ToySolver.Data.Polynomial.Interpolation.Lagrange as LagrangeInterpolation
 
 import qualified Data.Interval as Interval
@@ -910,6 +911,30 @@ case_Lagrange_interpolation_2 = p @?= q
         ]
     q = 6*x^2 - 11*x + 6
 
+-- https://en.wikipedia.org/wiki/Hermite_interpolation
+case_Hermite_interpolation = p @?= q
+  where
+    x :: UPolynomial Rational
+    x = P.var X
+    p = HermiteInterpolation.interpolate
+        [ (-1, [2, -8, 56])
+        , (0, [1, 0, 0])
+        , (1, [2, 8, 56])
+        ]
+    q = x^8 + 1
+
+prop_Hermite_interpolation_random =
+  forAll upolynomials $ \p ->
+    forAll (choose (0, 2)) $ \m ->
+    let d = P.deg p
+        -- m = 2
+        n = (d + 1 + m) `div` (m+1)
+        -- d <= n (m + 1) - 1
+        xs = genericTake n [-1, 0 ..]
+        ds = [(x', genericTake (m+1) [P.eval (\_ -> x') q | q <- iterate (\q -> P.deriv q X) p]) | x' <- xs]
+        p' = HermiteInterpolation.interpolate ds
+    in counterexample (show (p, ds, p')) $ p == p'
+
 -- ---------------------------------------------------------------------
 
 -- http://www.math.tamu.edu/~geller/factoring.pdf

toysolver.cabal

@@ -351,6 +351,7 @@ Library
      ToySolver.Data.Polynomial.Factorization.SquareFree
      ToySolver.Data.Polynomial.Factorization.Zassenhaus
      ToySolver.Data.Polynomial.GroebnerBasis
+     ToySolver.Data.Polynomial.Interpolation.Hermite
      ToySolver.Data.Polynomial.Interpolation.Lagrange
      ToySolver.FileFormat
      ToySolver.FileFormat.Base