{-# LANGUAGE RebindableSyntax #-}
{- |
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, dilate, shrink,
   ) 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.Reverse.StrictSpine as Rev
import qualified Data.List as List
import NumericPrelude.List (zipWithOverlap, )
import Data.Tuple.HT (mapPair, mapFst, forcePair, )
import Data.List.HT (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, )


{- $setup
>>> import qualified MathObj.Polynomial.Core as PolyCore
>>> import qualified MathObj.Polynomial as Poly
>>> import qualified Data.List as List
>>> import qualified Test.QuickCheck as QC
>>> import Test.QuickCheck ((==>))
>>> import Data.Tuple.HT (mapPair, mapSnd)
>>> import NumericPrelude.Numeric
>>> import NumericPrelude.Base
>>> import Prelude ()
>>>
>>> intPoly :: [Integer] -> [Integer]
>>> intPoly = id
>>>
>>> ratioPoly :: [Rational] -> [Rational]
>>> ratioPoly = id
-}


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

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

{-# INLINE hornerArgVector #-}
hornerArgVector :: (Module.C a v, Ring.C v) => v -> [a] -> v
hornerArgVector :: v -> [a] -> v
hornerArgVector v
x = (a -> v -> v) -> v -> [a] -> v
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\a
c v
val -> a
ca -> v -> v
forall a v. C a v => a -> v -> v
*>v
forall a. C a => a
onev -> v -> v
forall a. C a => a -> a -> a
+v
valv -> v -> v
forall a. C a => a -> a -> a
*v
x) v
forall a. C a => a
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 :: [a] -> [a]
normalize = (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
Rev.dropWhile a -> Bool
forall a. C a => a -> Bool
isZero

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

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

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


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

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


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


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


{- |
prop> \(QC.NonEmpty xs) (QC.NonEmpty ys) -> PolyCore.tensorProduct xs ys == List.transpose (PolyCore.tensorProduct ys (intPoly xs))
-}
{-# INLINE tensorProduct #-}
tensorProduct :: Ring.C a => [a] -> [a] -> [[a]]
tensorProduct :: [a] -> [a] -> [[a]]
tensorProduct = (a -> a -> a) -> [a] -> [a] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [[c]]
outerProduct a -> a -> a
forall a. C a => a -> a -> a
(*)

tensorProductAlt :: Ring.C a => [a] -> [a] -> [[a]]
tensorProductAlt :: [a] -> [a] -> [[a]]
tensorProductAlt [a]
xs [a]
ys = (a -> [a]) -> [a] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map ((a -> [a] -> [a]) -> [a] -> a -> [a]
forall a b c. (a -> b -> c) -> b -> a -> c
flip a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
scale [a]
ys) [a]
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 :: [a] -> [a] -> [a]
mul [] = [a] -> [a] -> [a]
forall a b. a -> b -> a
P.const []
mul [a]
xs = (a -> [a] -> [a]) -> [a] -> [a] -> [a]
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\a
y [a]
zs -> let (a
v:[a]
vs) = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
scale a
y [a]
xs in a
v a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
add [a]
vs [a]
zs) []
-- this one fails on infinite lists
--    mul xs = foldr (\y zs -> add (scale y xs) (shift zs)) []

{- |
prop> \xs ys  ->  PolyCore.equal (intPoly $ PolyCore.mul xs ys) (PolyCore.mulShear xs ys)
-}
{-# INLINE mulShear #-}
mulShear :: Ring.C a => [a] -> [a] -> [a]
mulShear :: [a] -> [a] -> [a]
mulShear [a]
xs [a]
ys = ([a] -> a) -> [[a]] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> a
forall a. C a => [a] -> a
sum ([[a]] -> [[a]]
forall a. [[a]] -> [[a]]
shear ([a] -> [a] -> [[a]]
forall a. C a => [a] -> [a] -> [[a]]
tensorProduct [a]
xs [a]
ys))

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


{- |
prop> \x y -> case (PolyCore.normalize x, PolyCore.normalize y) of (nx, ny) -> not (null (ratioPoly ny)) ==> mapSnd PolyCore.normalize (PolyCore.divMod nx ny) == mapPair (PolyCore.normalize, PolyCore.normalize) (PolyCore.divMod x y)
prop> \x y -> not (isZero (ratioPoly y)) ==> let z = fst $ PolyCore.divMod (Poly.coeffs x) y in  PolyCore.normalize z == z
prop> \x y -> case PolyCore.normalize $ ratioPoly y of ny -> not (null ny) ==> List.length (snd $ PolyCore.divMod x y) < List.length ny
-}
divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])
divMod :: [a] -> [a] -> ([a], [a])
divMod [a]
x [a]
y =
   ([a] -> [a], [a] -> [a]) -> ([a], [a]) -> ([a], [a])
forall a c b d. (a -> c, b -> d) -> (a, b) -> (c, d)
mapPair ([a] -> [a]
forall a. [a] -> [a]
List.reverse, [a] -> [a]
forall a. [a] -> [a]
List.reverse) (([a], [a]) -> ([a], [a])) -> ([a], [a]) -> ([a], [a])
forall a b. (a -> b) -> a -> b
$
   [a] -> [a] -> ([a], [a])
forall a. (C a, C a) => [a] -> [a] -> ([a], [a])
divModRev ([a] -> [a]
forall a. [a] -> [a]
List.reverse [a]
x) ([a] -> [a]
forall a. [a] -> [a]
List.reverse [a]
y)

{-
snd $ Poly.divMod (repeat (1::Double)) [1,1]
-}
{- |
The modulus will always have one element less than the divisor.
This means that the modulus will be denormalized in some cases,
e.g. @mod [2,1,1] [1,1,1] == [1,0]@ instead of @[1]@.
-}
divModRev :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])
divModRev :: [a] -> [a] -> ([a], [a])
divModRev [a]
x [a]
y =
   case (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile a -> Bool
forall a. C a => a -> Bool
isZero [a]
y of
      [] -> [Char] -> ([a], [a])
forall a. HasCallStack => [Char] -> a
error [Char]
"MathObj.Polynomial: division by zero"
      a
y0:[a]
ys ->
         let -- the second parameter represents lazily (length x - length (normalize y))
             aux :: [a] -> [b] -> ([a], [a])
aux [a]
xs' =
               ([a], [a]) -> ([a], [a])
forall a b. (a, b) -> (a, b)
forcePair (([a], [a]) -> ([a], [a]))
-> ([b] -> ([a], [a])) -> [b] -> ([a], [a])
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
               ([a], [a]) -> (b -> [b] -> ([a], [a])) -> [b] -> ([a], [a])
forall b a. b -> (a -> [a] -> b) -> [a] -> b
switchL
                 ([], [a]
xs')
                 (([b] -> ([a], [a])) -> b -> [b] -> ([a], [a])
forall a b. a -> b -> a
P.const (([b] -> ([a], [a])) -> b -> [b] -> ([a], [a]))
-> ([b] -> ([a], [a])) -> b -> [b] -> ([a], [a])
forall a b. (a -> b) -> a -> b
$
                    let (a
x0:[a]
xs) = [a]
xs'
                        q0 :: a
q0      = a
x0a -> a -> a
forall a. C a => a -> a -> a
/a
y0
                    in  ([a] -> [a]) -> ([a], [a]) -> ([a], [a])
forall a c b. (a -> c) -> (a, b) -> (c, b)
mapFst (a
q0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:) (([a], [a]) -> ([a], [a]))
-> ([b] -> ([a], [a])) -> [b] -> ([a], [a])
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [b] -> ([a], [a])
aux ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
sub [a]
xs (a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
scale a
q0 [a]
ys)))
         in  [a] -> [a] -> ([a], [a])
forall b. [a] -> [b] -> ([a], [a])
aux [a]
x (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop ([a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
ys) [a]
x)

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


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

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

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


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

{-# INLINE alternate #-}
alternate :: Additive.C a => [a] -> [a]
alternate :: [a] -> [a]
alternate = ((a -> a) -> a -> a) -> [a -> a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
($) ([a -> a] -> [a -> a]
forall a. [a] -> [a]
cycle [a -> a
forall a. a -> a
id, a -> a
forall a. C a => a -> a
Additive.negate])

{-# INLINE shrink #-}
shrink :: Ring.C a => a -> [a] -> [a]
shrink :: a -> [a] -> [a]
shrink a
k = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. C a => a -> a -> a
(*) ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a
ka -> a -> a
forall a. C a => a -> a -> a
*) a
forall a. C a => a
one)

{-# INLINE dilate #-}
dilate :: Field.C a => a -> [a] -> [a]
dilate :: a -> [a] -> [a]
dilate = a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
shrink (a -> [a] -> [a]) -> (a -> a) -> a -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a
forall a. C a => a -> a
Field.recip


{-
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
-}