Converting Conor McBride's Polynomial datatype to raising/falling factorial base.

Poly.hs

{-# LANGUAGE FlexibleInstances, DeriveFunctor #-}

import Control.Applicative
import Data.Unfolder
import Data.Unfoldable
import Data.Ratio
import Test.QuickCheck hiding (choose)


class Polynomial p where
  eval :: (Eq a, Enum a, Num a) => p a -> a -> a
  maxDegree :: p a -> Int
  
isEq :: (Polynomial p, Polynomial q) => p (Ratio Integer) -> q (Ratio Integer) -> Bool
isEq p q = all (\d -> let i = toInteger d % 1 in eval p i == eval q i) [0 .. max (maxDegree p) (maxDegree q)]


newtype Raising a = Raising [a] deriving (Functor, Show)

instance Polynomial Raising where
  eval (Raising cs) x = sum (zipWith (*) cs (scanl (*) 1 [x, x + 1..]))
  maxDegree (Raising []) = 0
  maxDegree (Raising cs) = length cs - 1

zipWithR :: Num a => (a -> a -> a) -> Raising a -> Raising a -> Raising a
zipWithR g (Raising a) (Raising b) = Raising (zipWith' g a b)


newtype Falling a = Falling [a] deriving (Functor, Show)

instance Polynomial Falling where
  eval (Falling cs) x = sum (zipWith (*) cs (scanl (*) 1 [x, x - 1..]))
  maxDegree (Falling []) = 0
  maxDegree (Falling cs) = length cs - 1

zipWithF :: Num a => (a -> a -> a) -> Falling a -> Falling a -> Falling a
zipWithF g (Falling a) (Falling b) = Falling (zipWith' g a b)


zipWith' :: Num a => (a -> a -> a) -> [a] -> [a] -> [a]
zipWith' _ [] [] = []
zipWith' f [] (y:ys) = f 0 y : zipWith' f [] ys
zipWith' f (x:xs) [] = f x 0 : zipWith' f xs []
zipWith' f (x:xs) (y:ys) = f x y : zipWith' f xs ys


data Poly a = P0 | P1 | Poly a :+: Poly a | a :* Poly a | Sh (Poly a) | X (Poly a) | Sum (Poly a)
  deriving (Functor, Show)

instance Polynomial Poly where
  eval P0 _ = 0
  eval P1 _ = 1
  eval (p :+: q) x = eval p x + eval q x
  eval (k :* p) x = k * eval p x
  eval (Sh p) x = eval p (x + 1)
  eval (X p) x = x * eval p x
  eval (Sum p) x = sm (eval p) x
  maxDegree P0 = 0
  maxDegree P1 = 0
  maxDegree (p :+: q) = max (maxDegree p) (maxDegree q)
  maxDegree (_ :* p) = maxDegree p
  maxDegree (Sh p) = maxDegree p
  maxDegree (X p) = 1 + maxDegree p
  maxDegree (Sum p) = 1 + maxDegree p

sm :: (Eq a, Num a) => (a -> a) -> (a -> a)
sm _ 0 = 0
sm f n = sm f (n - 1) + f (n - 1)

diff :: Poly a -> Poly a
diff P0 = P0
diff P1 = P0
diff (p :+: q) = diff p :+: diff q
diff (k :* p) = k :* diff p
diff (Sh p) = Sh (diff p)
diff (X p) = Sh p :+: X (diff p)
diff (Sum p) = p

lemma :: Poly (Ratio Integer) -> Bool
lemma p = isEq (Sh p) (p :+: diff p)


poly2Raising :: Poly (Ratio Integer) -> Raising (Ratio Integer)
poly2Raising P0 = Raising []
poly2Raising P1 = Raising [1]
poly2Raising (p :+: q) = zipWithR (+) (poly2Raising p) (poly2Raising q)
poly2Raising (k :* p) = (k *) `fmap` poly2Raising p
poly2Raising (Sh p) = let Raising cs = poly2Raising p in Raising $ h cs 0 (\_ t -> t)
  where
    h [] _ k = k 0 []
    h (c:cs) n k = h cs (n + 1) (\d t -> k (n * (c + d)) (c + d : t)) 
poly2Raising (X p) = let Raising cs = poly2Raising p in Raising $ zipWith (-) (0:cs) $ zipWith (*) cs [0..] ++ [0]
poly2Raising (Sum p) = let Raising cs = poly2Raising p in Raising $ 0 : zipWith (/) cs [1..] -- Wrong
  -- But works if sm f n = sm f (n - 1) + f n

poly2Falling :: Poly (Ratio Integer) -> Falling (Ratio Integer)
poly2Falling P0 = Falling []
poly2Falling P1 = Falling [1]
poly2Falling (p :+: q) = zipWithF (+) (poly2Falling p) (poly2Falling q)
poly2Falling (k :* p) = (k *) `fmap` poly2Falling p
poly2Falling (Sh p) = let Falling cs = poly2Falling p in Falling . zipWith (+) cs $ zipWith (*) (tail cs) [1..] ++ [0]
poly2Falling (X p) = let Falling cs = poly2Falling p in Falling . (0 :) . zipWith (+) cs $ zipWith (*) (tail cs) [1..] ++ [0]
poly2Falling (Sum p) = let Falling cs = poly2Falling p in Falling $ 0 : zipWith (/) cs [1..]


instance Unfoldable Poly where
  unfold fa = choose
    [ pure P0
    , pure P1
    , (:+:) <$> unfold fa <*> unfold fa
    , (:*) <$> fa <*> unfold fa
    , Sh <$> unfold fa
    , X <$> unfold fa
    , Sum <$> unfold fa
    ]

instance Arbitrary a => Arbitrary (Poly a) where
  arbitrary = arbitraryDefault

main :: IO ()
main = do
  quickCheck (\p -> isEq p (poly2Falling p))
  quickCheck lemma