```{-# LANGUAGE RebindableSyntax #-}

-----------------------------------------------------------------------------
-- |
-- Module      :  Data.YAP.Polynomial
-- Copyright   :  (c) Ross Paterson 2011
-- License     :  BSD-style (see the file libraries/base/LICENSE)
--
-- Maintainer  :  ross@soi.city.ac.uk
-- Stability   :  provisional
-- Portability :  portable
--
-- An example instance of the new classes: polynomials.
-- Some of these functions work with infinite polynomials,
-- i.e. formal power series.
--
-----------------------------------------------------------------------------

module Data.YAP.Polynomial (
-- * Polynomials
Polynomial,
polynomial,
-- * Queries
-- ** Finite polynomials
coefficients,
degree,
evaluate,
pretty,
-- ** Formal power series
approximations,
-- * Operations
compose,
differentiate,
integrate
) where

import Prelude.YAP
import Data.YAP.Algebra

-- Coefficients, least significant first.
-- There may be trailing zeroes, but the interface hides them.
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)

-- | If @b@ is non-zero, @'mod' a b@ has a smaller degree than @b@.
-- If @a@ is non-zero, @'associate' a@ has a leading coefficient of @1@.
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 (n-1) as bs
divModAux n (a:as) bs@(b:bs') = (c:d, m)
where c = a/b
(d, m) = divModAux (n-1) (sub as (map (c*) bs')) bs

sub :: (AbelianGroup a) => [a] -> [a] -> [a]
sub xs [] = xs
sub (x:xs) (y:ys) = (x-y) : sub xs ys
sub [] _ = error "can't happen"

-- | Construct a polynomial from a list of coefficients,
-- least significant first.
polynomial :: [a] -> Polynomial a
polynomial = P

-- | The coefficients of a finite polynomial, least significant first
-- and with no trailing zeroes.
coefficients :: (Eq a, AbelianGroup a) => Polynomial a -> [a]
coefficients (P as) = reverse (dropWhile (== zero) (reverse as))

-- | The coefficients of a finite polynomial, starting with the most
-- significant non-zero coefficient.
rev_coefficients :: (Eq a, AbelianGroup a) => Polynomial a -> [a]
rev_coefficients (P as) = dropWhile (== zero) (reverse as)

-- | The degree of a finite polynomial.
--
-- @'degree' p = length (coefficients p)@
degree :: (Eq a, AbelianGroup a) => Polynomial a -> Int
degree a = length (rev_coefficients a)

-- | Evaluate a polynomial for a given value of @x@.
--
-- @'evaluate' a x = 'zipWith' (*) ('coefficients' a) ('iterate' (*x) 1)@
--
-- (The implementation uses Horner's rule.)
evaluate :: (Ring a) => Polynomial a -> a -> a
evaluate (P as) x = foldr (\ a v -> a + x*v) zero as

-- | The infinite list of evaluations of truncations of the polynomial
-- or power series.
approximations :: (Ring a) => Polynomial a -> a -> [a]
approximations (P as) x = scanr1 (+) (zipWith (*) as (iterate (*x) 1))

-- other functions

-- | Pretty-print a polynomial, e.g.
--
-- @pretty (polynomial [3, 4, 0, 1, 5]) \"x\" = \"5x^4 + x^3 + 4x + 3\"@
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

-- | Composition of polynomials:
--
-- @'evaluate' ('compose' a b) = 'evaluate' a . 'evaluate' b@
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]

-- | Symbolic differentiation of polynomials.
differentiate :: (Ring a) => Polynomial a -> Polynomial a
differentiate (P as) = P (zipWith (*) (tail as) (map fromInteger [1..]))

-- | Symbolic integration of polynomials.
integrate :: (Field a) => Polynomial a -> Polynomial a
integrate (P as) = P (zero : zipWith (/) as (map fromInteger [1..]))
```