{-# LANGUAGE NoImplicitPrelude #-}
{- |
This module implements polynomial functions on plain lists.
We use such functions in order to implement methods of other datatypes.

The module organization differs from that of @ResidueClass@:
Here the @Polynomial@ module exports the type
that fits to the NumericPrelude type classes,
whereas in @ResidueClass@ the sub-modules export various flavors of them.
-}
module MathObj.Polynomial.Core (
   horner, hornerCoeffVector, hornerArgVector,
   normalize,
   shift, unShift,
   equal,
   add, sub, negate,
   scale, collinear,
   tensorProduct, tensorProductAlt,
   mul, mulShear, mulShearTranspose,
   divMod, divModRev,
   stdUnit,
   progression, differentiate, integrate, integrateInt,
   mulLinearFactor,
   alternate,
   ) where

import qualified Algebra.Module               as Module
import qualified Algebra.Field                as Field
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 Data.List as List
import NumericPrelude.List (zipWithOverlap, )
import Data.Tuple.HT (mapPair, mapFst, forcePair, )
import Data.List.HT
          (dropWhileRev, switchL, shear, shearTranspose, outerProduct, )

import qualified NumericPrelude.Base as P
import qualified NumericPrelude.Numeric as NP

import NumericPrelude.Base    hiding (const, reverse, )
import NumericPrelude.Numeric hiding (divMod, negate, stdUnit, )


{- |
Horner's scheme for evaluating a polynomial in a ring.
-}
{-# INLINE horner #-}
horner :: Ring.C a => a -> [a] -> a
horner x = foldr (\c val -> c+x*val) zero

{- |
Horner's scheme for evaluating a polynomial in a module.
-}
{-# INLINE hornerCoeffVector #-}
hornerCoeffVector :: Module.C a v => a -> [v] -> v
hornerCoeffVector x = foldr (\c val -> c+x*>val) zero

{-# INLINE hornerArgVector #-}
hornerArgVector :: (Module.C a v, Ring.C v) => v -> [a] -> v
hornerArgVector x = foldr (\c val -> c*>one+val*x) zero


{- |
It's also helpful to put a polynomial in canonical form.
'normalize' strips leading coefficients that are zero.
-}
{-# INLINE normalize #-}
normalize :: (ZeroTestable.C a) => [a] -> [a]
normalize = dropWhileRev isZero

{- |
Multiply by the variable, used internally.
-}
{-# INLINE shift #-}
shift :: (Additive.C a) => [a] -> [a]
shift [] = []
shift l  = zero : l

{-# INLINE unShift #-}
unShift :: [a] -> [a]
unShift []     = []
unShift (_:xs) = xs

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


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

{-# INLINE negate #-}
negate :: (Additive.C a) => [a] -> [a]
negate = map NP.negate


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


collinear :: (Eq a, Ring.C a) => [a] -> [a] -> Bool
collinear (x:xs) (y:ys) =
   if x==zero && y==zero
     then collinear xs ys
     else scale x ys == scale y xs
-- here at least one of xs and ys is empty
collinear xs ys =
   all (==zero) xs && all (==zero) ys


{-# INLINE tensorProduct #-}
tensorProduct :: Ring.C a => [a] -> [a] -> [[a]]
tensorProduct = outerProduct (*)

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


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

{-# INLINE mul #-}
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)) []

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

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


divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])
divMod x y =
   mapPair (List.reverse, List.reverse) $
   divModRev (List.reverse x) (List.reverse y)

{-
snd $ Poly.divMod (repeat (1::Double)) [1,1]
-}
divModRev :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])
divModRev x y =
   let (y0:ys) = dropWhile isZero y
       -- the second parameter represents lazily (length x - length y)
       aux xs' =
         forcePair .
         switchL
           ([], xs')
           (P.const $
              let (x0:xs) = xs'
                  q0      = x0/y0
              in  mapFst (q0:) . aux (sub xs (scale q0 ys)))
   in  if isZero y
         then error "MathObj.Polynomial: division by zero"
         else aux x (drop (length y - 1) x)

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


{-# INLINE progression #-}
progression :: Ring.C a => [a]
progression = iterate (one+) one

{-# INLINE differentiate #-}
differentiate :: (Ring.C a) => [a] -> [a]
differentiate = zipWith (*) progression . drop 1

{-# INLINE integrate #-}
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.
-}
{-# INLINE integrateInt #-}
integrateInt :: (ZeroTestable.C a, Integral.C a) => a -> [a] -> [a]
integrateInt c x =
   c : zipWith Integral.divChecked x progression


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

{-# INLINE alternate #-}
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 (divChecked (degree*(degree-1)) 2)
                  (resultant xs (differentiate xs))
          `divChecked` last xs
-}