module Data.YAP.Polynomial (
Polynomial,
polynomial,
coefficients,
degree,
evaluate,
pretty,
approximations,
compose,
differentiate,
integrate
) where
import Prelude.YAP
import Data.YAP.Algebra
newtype Polynomial a = P [a]
instance Functor Polynomial where
fmap f (P as) = P (map f as)
instance (Eq a, AbelianGroup a) => Eq (Polynomial a) where
P as == P bs = eq as bs
where eq [] ys = allZero ys
eq xs [] = allZero xs
eq (x:xs) (y:ys) = x == y && eq xs ys
allZero :: (Eq a, AbelianGroup a) => [a] -> Bool
allZero = all (== zero)
instance (Eq a, Show a, AbelianGroup a) => Show (Polynomial a) where
showsPrec p x =
showParen (p > 10) $ showString "polynomial " . shows (coefficients x)
instance (AbelianGroup a) => AbelianGroup (Polynomial a) where
zero = P []
P xs + P ys = P (add xs ys)
negate (P xs) = P (map negate xs)
add :: (AbelianGroup a) => [a] -> [a] -> [a]
add [] ys = ys
add xs [] = xs
add (x:xs) (y:ys) = x+y : add xs ys
instance (Ring a) => Ring (Polynomial a) where
P xs * P ys = P (mul xs ys)
fromInteger i = P [fromInteger i]
mul :: (Ring a) => [a] -> [a] -> [a]
mul [] _ys = []
mul _xs [] = []
mul (x:xs) ys = add (map (x*) ys) (zero:mul xs ys)
instance (Eq a, Field a) => EuclideanDomain (Polynomial a) where
divMod a b
| null bs = error "division by zero"
| otherwise = (P (reverse d), P (reverse m))
where (d, m) = divModAux n as bs
as = rev_coefficients a
bs = rev_coefficients b
n = length as length bs
associate (P as)
| allZero as = P []
| otherwise = P (map (/ last as) (init as) ++ [1])
unit (P as)
| allZero as = P [1]
| otherwise = P [last as]
divModAux :: (Eq a, Field a) => Int -> [a] -> [a] -> ([a],[a])
divModAux n as _
| n < 0 = ([], as)
divModAux n (0:as) bs = (0:d, m)
where (d, m) = divModAux (n1) as bs
divModAux n (a:as) bs@(b:bs') = (c:d, m)
where c = a/b
(d, m) = divModAux (n1) (sub as (map (c*) bs')) bs
sub :: (AbelianGroup a) => [a] -> [a] -> [a]
sub xs [] = xs
sub (x:xs) (y:ys) = (xy) : sub xs ys
sub [] _ = error "can't happen"
polynomial :: [a] -> Polynomial a
polynomial = P
coefficients :: (Eq a, AbelianGroup a) => Polynomial a -> [a]
coefficients (P as) = reverse (dropWhile (== zero) (reverse as))
rev_coefficients :: (Eq a, AbelianGroup a) => Polynomial a -> [a]
rev_coefficients (P as) = dropWhile (== zero) (reverse as)
degree :: (Eq a, AbelianGroup a) => Polynomial a -> Int
degree a = length (rev_coefficients a)
evaluate :: (Ring a) => Polynomial a -> a -> a
evaluate (P as) x = foldr (\ a v -> a + x*v) zero as
approximations :: (Ring a) => Polynomial a -> a -> [a]
approximations (P as) x = scanr1 (+) (zipWith (*) as (iterate (*x) 1))
pretty :: (Ord a, Show a, Ring a) => Polynomial a -> String -> String
pretty (P as) x = case dropWhile ((== zero) . fst) (reverse (zip as terms)) of
[] -> "0"
(a,t):ats -> showFirst a ++ t ++ showRest ats
where terms = "" : x : [x ++ "^" ++ show n | n <- [2..]::[Int]]
showFirst a
| a < 0 = '-':show (negate a)
| a == 1 = ""
| otherwise = show a
showRest [] = ""
showRest [(a,t)]
| a < 0 = " - " ++ show (negate a) ++ t
| a > 0 = " + " ++ show a ++ t
| otherwise = ""
showRest ((a,t):ats)
| a < 0 = " - " ++ show (negate a) ++ t ++ showRest ats
| a == 1 = " + " ++ t ++ showRest ats
| a > 0 = " + " ++ show a ++ t ++ showRest ats
| otherwise = showRest ats
compose :: (Ring a) => Polynomial a -> Polynomial a -> Polynomial a
compose (P as) q = evaluate (P (map constant as)) q
constant :: a -> Polynomial a
constant a = P [a]
differentiate :: (Ring a) => Polynomial a -> Polynomial a
differentiate (P as) = P (zipWith (*) (tail as) (map fromInteger [1..]))
integrate :: (Field a) => Polynomial a -> Polynomial a
integrate (P as) = P (zero : zipWith (/) as (map fromInteger [1..]))