{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-}

{- |
Polynomials and rational functions in a single indeterminate.
Polynomials are represented by a list of coefficients.
All non-zero coefficients are listed, but there may be extra '0's at the end.

Usage:
Say you have the ring of 'Integer' numbers
and you want to add a transcendental element @x@,
that is an element, which does not allow for simplifications.
More precisely, for all positive integer exponents @n@
the power @x^n@ cannot be rewritten as a sum of powers with smaller exponents.
The element @x@ must be represented by the polynomial @[0,1]@.

In principle, you can have more than one transcendental element
by using polynomials whose coefficients are polynomials as well.
However, most algorithms on multi-variate polynomials
prefer a different (sparse) representation,
where the ordering of elements is not so fixed.

If you want division, you need "Number.Ratio"s
of polynomials with coefficients from a "Algebra.Field".

You can also compute with an algebraic element,
that is an element which satisfies an algebraic equation like
@x^3-x-1==0@.
Actually, powers of @x@ with exponents above @3@ can be simplified,
since it holds @x^3==x+1@.
You can perform these computations with "Number.ResidueClass" of polynomials,
where the divisor is the polynomial equation that determines @x@.
If the polynomial is irreducible
(in our case @x^3-x-1@ cannot be written as a non-trivial product)
then the residue classes also allow unrestricted division
(except by zero, of course).
That is, using residue classes of polynomials
you can work with roots of polynomial equations
without representing them by radicals
(powers with fractional exponents).
It is well-known, that roots of polynomials of degree above 4
may not be representable by radicals.
-}

module MathObj.Polynomial(T(..), fromCoeffs, showsExpressionPrec, const,
                  eval, compose, equal, add, sub, negate,
                  shift, unShift,
                  mul, scale, divMod,
                  tensorProduct, tensorProduct',
                  mulShear, mulShearTranspose,
                  horner, horner',
                  progression, differentiate, integrate, integrateInt,
                  fromRoots, alternate)
where

import qualified Algebra.Differential         as Differential
import qualified Algebra.VectorSpace          as VectorSpace
import qualified Algebra.Module               as Module
import qualified Algebra.Vector               as Vector
import qualified Algebra.Field                as Field
import qualified Algebra.PrincipalIdealDomain as PID
import qualified Algebra.Units                as Units
import qualified Algebra.IntegralDomain       as Integral
import qualified Algebra.Ring                 as Ring
import qualified Algebra.Additive             as Additive
import qualified Algebra.ZeroTestable         as ZeroTestable
import qualified Algebra.Indexable            as Indexable

import Algebra.Module((*>))
import Algebra.ZeroTestable(isZero)

import Control.Monad (liftM)
import qualified Data.List as List
import NumericPrelude.List
   (zipWithOverlap, dropWhileRev, shear, shearTranspose, outerProduct)

import Test.QuickCheck (Arbitrary(arbitrary,coarbitrary))

import qualified Prelude     as P98
import qualified PreludeBase as P
import qualified NumericPrelude as NP

import PreludeBase    hiding (const)
import NumericPrelude hiding (divMod, negate, stdUnit)

newtype T a = Cons {coeffs :: [a]}

fromCoeffs :: [a] -> T a
fromCoeffs = lift0

lift0 :: [a] -> T a
lift0 = Cons

lift1 :: ([a] -> [a]) -> (T a -> T a)
lift1 f (Cons x0) = Cons (f x0)

lift2 :: ([a] -> [a] -> [a]) -> (T a -> T a -> T a)
lift2 f (Cons x0) (Cons x1) = Cons (f x0 x1)

{-
Functor instance is e.g. useful for showing polynomials in residue rings.
@fmap (ResidueClass.concrete 7) (polynomial [1,4,4::ResidueClass.T Integer] * polynomial [1,5,6])@
-}

instance Functor T where
  fmap f (Cons xs) = Cons (map f xs)

plusPrec, appPrec :: Int
plusPrec = 6
appPrec  = 10

instance (Show a) => Show (T a) where
  showsPrec p (Cons xs) =
    showParen (p >= appPrec) (showString "Polynomial.fromCoeffs " . shows xs)

showsExpressionPrec :: (Show a, ZeroTestable.C a, Additive.C a) =>
   Int -> String -> T a -> String -> String
showsExpressionPrec p var poly =
    if isZero poly
      then showString "0"
      else
        let terms = filter (not . isZero . fst)
                       (zip (coeffs poly) monomials)
            monomials = id :
                        showString "*" . showString var :
                        map (\k -> showString "*" . showString var
                                 . showString "^" . shows k)
                            [(2::Int)..]
            showsTerm x showsMon = showsPrec (plusPrec+1) x . showsMon
        in showParen (p > plusPrec)
           (foldl (.) id $ List.intersperse (showString " + ") $
            map (uncurry showsTerm) terms)

{- |
Horner's scheme for evaluating an polynomial
-}

horner' :: Ring.C a => a -> [a] -> a
horner' x = foldr (\c val -> c+x*val) zero

{-
***** Module types are more general,
but most times this flexibility is not needed and
let type inference fail.
-}

horner :: Module.C a b => a -> [b] -> b
horner x = foldr (\c val -> c+x*>val) zero

eval :: Module.C a b => T b -> a -> b
eval (Cons y) x = horner x y

{- |
'compose' is the functional composition of polynomials.

It fulfills
  @ eval x . eval y == eval (compose x y) @
-}

-- compose :: Module.C a b => T b -> T a -> T a
-- compose (Cons x) y = horner y (map const x)
compose :: (Ring.C a) => T a -> T a -> T a
compose (Cons x) y = horner' y (map const x)

{- |
It's also helpful to put a polynomial in canonical form.
'normalize' strips leading coefficients that are zero.
-}

normalize :: (ZeroTestable.C a) => [a] -> [a]
normalize = dropWhileRev isZero

{- |
Multiply by the variable, used internally.
-}

shift :: (Additive.C a) => [a] -> [a]
shift [] = []
shift l  = zero : l

unShift :: [a] -> [a]
unShift []     = []
unShift (_:xs) = xs

const :: a -> T a
const x = lift0 [x]

equal :: (Eq a, ZeroTestable.C a) => [a] -> [a] -> Bool
equal x y = and (zipWithOverlap isZero isZero (==) x y)

instance (Eq a, ZeroTestable.C a) => Eq (T a) where
  (Cons x) == (Cons y) = equal x y

instance (Indexable.C a, ZeroTestable.C a) => Indexable.C (T a) where
  compare = Indexable.liftCompare coeffs

instance (ZeroTestable.C a) => ZeroTestable.C (T a) where
  isZero (Cons x) = isZero x


add, sub :: (Additive.C a) => [a] -> [a] -> [a]
add = (+)
sub = (-)

negate :: (Additive.C a) => [a] -> [a]
negate = map NP.negate

instance (Additive.C a) => Additive.C (T a) where
  (+)    = lift2 add
  (-)    = lift2 sub
  zero   = lift0 []
  negate = lift1 negate


scale :: Ring.C a => a -> [a] -> [a]
scale s = map (s*)


instance Vector.C T where
   zero  = zero
   (<+>) = (+)
   (*>)  = Vector.functorScale

instance (Module.C a b) => Module.C a (T b) where
   (*>) x = lift1 (x *>)

instance (Field.C a, Module.C a b) => VectorSpace.C a (T b)


tensorProduct :: Ring.C a => [a] -> [a] -> [[a]]
tensorProduct = outerProduct (*)

tensorProduct' :: Ring.C a => [a] -> [a] -> [[a]]
tensorProduct' xs ys = map (flip scale ys) xs

{- |
'mul' is fast if the second argument is a short polynomial,
'MathObj.PowerSeries.**' relies on that fact.
-}

mul :: Ring.C a => [a] -> [a] -> [a]
{- prevent from generation of many zeros
   if the first operand is the empty list -}
mul [] = P.const []
mul xs = foldr (\y zs -> let (v:vs) = scale y xs in v : add vs zs) []
-- this one fails on infinite lists
--    mul xs = foldr (\y zs -> add (scale y xs) (shift zs)) []

mulShear :: Ring.C a => [a] -> [a] -> [a]
mulShear xs ys = map sum (shear (tensorProduct xs ys))

mulShearTranspose :: Ring.C a => [a] -> [a] -> [a]
mulShearTranspose xs ys = map sum (shearTranspose (tensorProduct xs ys))

instance (Ring.C a) => Ring.C (T a) where
  one         = const one
  fromInteger = const . fromInteger
  (*)         = lift2 mul


divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])
divMod x y =
    let (y0:ys) = dropWhile isZero (reverse y)
        aux l xs' =
          if l < 0
            then ([], xs')
            else
              let (x0:xs) = xs'
                  q0      = x0/y0
                  (d',m') = aux (l-1) (sub xs (scale q0 ys))
              in  (q0:d',m')
        (d, m) = aux (length x - length y) (reverse x)
    in  if isZero y
          then error "MathObj.Polynomial: division by zero"
          else (reverse d, reverse m)

instance (ZeroTestable.C a, Field.C a) => Integral.C (T a) where
  divMod (Cons x) (Cons y) =
     let (d,m) = divMod x y
     in  (Cons d, Cons m)

stdUnit :: (ZeroTestable.C a, Ring.C a) => [a] -> a
stdUnit x = case normalize x of
    [] -> one
    l  -> last l

instance (ZeroTestable.C a, Field.C a) => Units.C (T a) where
  isUnit (Cons []) = False
  isUnit (Cons (x0:xs)) = not (isZero x0) && all isZero xs
  stdUnit    (Cons x) = const        (stdUnit x)
  stdUnitInv (Cons x) = const (recip (stdUnit x))

{-
Polynomials are a Euclidean domain, so no instance is necessary
(although it might be faster).
-}

instance (ZeroTestable.C a, Field.C a) => PID.C (T a)

progression :: Ring.C a => [a]
progression = iterate (one+) one

differentiate :: (Ring.C a) => [a] -> [a]
differentiate = zipWith (*) progression . tail

integrate :: (Field.C a) => a -> [a] -> [a]
integrate c x = c : zipWith (/) x progression

{- |
Integrates if it is possible to represent the integrated polynomial
in the given ring.
Otherwise undefined coefficients occur.
-}
integrateInt :: (ZeroTestable.C a, Integral.C a) => a -> [a] -> [a]
integrateInt c x =
   c : zipWith Integral.safeDiv x progression


instance (Ring.C a) => Differential.C (T a) where
  differentiate = lift1 differentiate


fromRoots :: (Ring.C a) => [a] -> T a
fromRoots = Cons . foldl (flip mulLinearFactor) [1]

mulLinearFactor :: Ring.C a => a -> [a] -> [a]
mulLinearFactor x yt@(y:ys) = Additive.negate (x*y) : yt - scale x ys
mulLinearFactor _ [] = []

alternate :: Additive.C a => [a] -> [a]
alternate = zipWith ($) (cycle [id, Additive.negate])

{-
see htam: Wavelet/DyadicResultant

resultant :: Ring.C a => [a] -> [a] -> [a]
resultant xs ys =

discriminant :: Ring.C a => [a] -> [a]
discriminant xs =
   let degree = genericLength xs
   in  parityFlip (safeDiv (degree*(degree-1)) 2)
                  (resultant xs (differentiate xs))
          `safeDiv` last xs
-}

instance (Arbitrary a, ZeroTestable.C a) => Arbitrary (T a) where
   arbitrary = liftM (fromCoeffs . normalize) arbitrary
   coarbitrary = undefined


{- * legacy instances -}

{- |
It is disputable whether polynomials shall be represented by number literals or not.
An advantage is, that one can write
let x = polynomial [0,1]
in  (x^2+x+1)*(x-1)
However the output looks much different.
-}
legacyInstance :: a
legacyInstance =
   error "legacy Ring.C instance for simple input of numeric literals"

instance (Ring.C a, Eq a, Show a, ZeroTestable.C a) => P98.Num (T a) where
   fromInteger = const . fromInteger
   negate = Additive.negate -- for unary minus
   (+)    = legacyInstance
   (*)    = legacyInstance
   abs    = legacyInstance
   signum = legacyInstance

instance (Field.C a, Eq a, Show a, ZeroTestable.C a) => P98.Fractional (T a) where
   fromRational = const . fromRational
   (/) = legacyInstance