-- |
-- Module: Symtegration.Polynomial
-- Description: Polynomials for Symtegration.
-- Copyright: Copyright 2024 Yoo Chung
-- License: Apache-2.0
-- Maintainer: dev@chungyc.org
--
-- This modules defines a type class that concrete types representing polynomials
-- should be an instance of.  It includes important algorithms operating on
-- polynomials.  In particular, algorithms for polynomial division and
-- the extended Euclidean algorithm are included.
module Symtegration.Polynomial
  ( -- * Polynomials
    Polynomial (..),
    monic,
    mapCoefficients,

    -- * Algorithms
    divide,
    pseudoDivide,
    extendedEuclidean,
    diophantineEuclidean,
    greatestCommonDivisor,
    subresultant,
    differentiate,
    integrate,
    squarefree,
  )
where

import Data.Monoid (Sum (..))

-- $setup
-- >>> import Data.Ratio ((%))
-- >>> import Symtegration.Symbolic
-- >>> import Symtegration.Symbolic.Simplify
-- >>> import Symtegration.Polynomial.Indexed

-- | Polynomials must support the operations specified in this type class.
-- All powers must be non-negative.
class (Integral e, Num c) => Polynomial p e c where
  -- | Returns the degree of a given polynomial.
  --
  -- The following returns 9 for the highest term in \(3x^9 + 2x^4 + x\):
  --
  -- >>> degree (3 * power 9 + 2 * power 4 + power 1 :: IndexedPolynomial)
  -- 9
  degree :: p e c -> e

  -- | Returns the coefficient for the term with the given power.
  --
  -- The following returns 4 from the \(4x^3\) term in \(x^4 + 4x^3 + 3\):
  --
  -- >>> coefficient (power 4 + 4 * power 3 + 3 :: IndexedPolynomial) 3
  -- 4 % 1
  coefficient :: p e c -> e -> c

  -- | Returns the leading coefficient.
  --
  -- The following returns 6 from the \(6x^3\) term in \(6x^3 + 2x^2\):
  --
  -- >>> leadingCoefficient (6 * power 3 + 2 * power 2 :: IndexedPolynomial)
  -- 6 % 1
  --
  -- The leading coefficient is never zero unless the polynomial itself is zero.
  leadingCoefficient :: p e c -> c

  -- | Returns the polynomial without the leading term.
  --
  -- >>> deleteLeadingTerm (2 * power 3 + power 1 + 2 :: IndexedPolynomial)
  -- x + 2
  deleteLeadingTerm :: p e c -> p e c

  -- | Fold the terms, i.e., the powers and coefficients, using the given monoid.
  -- Only terms with non-zero coefficients will be folded.
  -- Folding is ordered from lower to higher terms.
  --
  -- For example with \(3x^5 - 2x + 7\),
  --
  -- >>> foldTerms (\e c -> show (e, c)) (3 * power 5 - 2 * power 1 + 7 :: IndexedPolynomial)
  -- "(0,7 % 1)(1,(-2) % 1)(5,3 % 1)"
  foldTerms :: (Monoid m) => (e -> c -> m) -> p e c -> m

  -- | Multiplies a polynomial by a scalar.
  --
  -- The following divides \(6x + 2\) by 2:
  --
  -- >>> scale (1 % 2) (6 * power 1 + 2 :: IndexedPolynomial)
  -- 3x + 1
  scale :: c -> p e c -> p e c

  -- | Returns a single term with the variable raised to the given power.
  --
  -- The following is equivalent to \(x^5\):
  --
  -- >>> power 5 :: IndexedPolynomial
  -- x^5
  power :: e -> p e c

-- | Scale the polynomial so that its leading coefficient is one.
--
-- >>> monic $ 4 * power 2 + 4 * power 1 + 4 :: IndexedPolynomial
-- x^2 + x + 1
--
-- The exception is when the polynomial is zero.
--
-- >>> monic 0 :: IndexedPolynomial
-- 0
monic :: (Polynomial p e c, Eq c, Fractional c) => p e c -> p e c
monic :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic p e c
p
  | p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e c
p c -> c -> Bool
forall a. Eq a => a -> a -> Bool
== c
0 = p e c
p
  | Bool
otherwise = c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (c
1 c -> c -> c
forall a. Fractional a => a -> a -> a
/ p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e c
p) p e c
p

-- | Maps the coefficients in a polynomial to form another polynomial.
--
-- For example, it can be used to convert a polynomial with 'Rational' coefficients
-- into a polynomial with 'Expression' coefficients.
--
-- >>> let p = 2 * power 1 + 1 :: IndexedPolynomial
-- >>> let q = mapCoefficients fromRational p :: IndexedSymbolicPolynomial
-- >>> simplify $ coefficient q 1
-- Number 2
mapCoefficients ::
  (Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c')) =>
  (c -> c') ->
  p e c ->
  p e c'
mapCoefficients :: forall (p :: * -> * -> *) e c c'.
(Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c')) =>
(c -> c') -> p e c -> p e c'
mapCoefficients c -> c'
f p e c
p = Sum (p e c') -> p e c'
forall a. Sum a -> a
getSum (Sum (p e c') -> p e c') -> Sum (p e c') -> p e c'
forall a b. (a -> b) -> a -> b
$ (e -> c -> Sum (p e c')) -> p e c -> Sum (p e c')
forall m. Monoid m => (e -> c -> m) -> p e c -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms e -> c -> Sum (p e c')
forall {p :: * -> * -> *} {e}.
Polynomial p e c' =>
e -> c -> Sum (p e c')
convertTerm p e c
p
  where
    convertTerm :: e -> c -> Sum (p e c')
convertTerm e
e c
c = p e c' -> Sum (p e c')
forall a. a -> Sum a
Sum (p e c' -> Sum (p e c')) -> p e c' -> Sum (p e c')
forall a b. (a -> b) -> a -> b
$ c' -> p e c' -> p e c'
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (c -> c'
f c
c) (e -> p e c'
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power e
e)

-- | Polynomial division.  It returns the quotient polynomial and the remainder polynomial.
--
-- For example, dividing \(p = x^3-12x^2-42\) by \(q = x^2 - 2x + 1\)
-- returns \(x-10\) as the quotient and \(-21x-32\) as the remainder,
-- since \(p = (x-10)q -21x - 32\):
--
-- >>> let p = power 3 - 12 * power 2 - 42 :: IndexedPolynomial
-- >>> let q = power 2 - 2 * power 1 + 1 :: IndexedPolynomial
-- >>> divide p q
-- (x + (-10),(-21)x + (-32))
divide ::
  (Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
  -- | Dividend polynomial being divided.
  p e c ->
  -- | Divisor polynomial dividing the dividend.
  p e c ->
  -- | Quotient and remainder.
  (p e c, p e c)
divide :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
divide p e c
p p e c
q = p e c -> p e c -> (p e c, p e c)
go p e c
0 p e c
p
  where
    go :: p e c -> p e c -> (p e c, p e c)
go p e c
quotient p e c
remainder
      | p e c
remainder p e c -> p e c -> Bool
forall a. Eq a => a -> a -> Bool
/= p e c
0, e
delta e -> e -> Bool
forall a. Ord a => a -> a -> Bool
>= e
0 = p e c -> p e c -> (p e c, p e c)
go (p e c
quotient p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
+ p e c
t) (p e c
remainder' p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
- p e c
qt')
      | Bool
otherwise = (p e c
quotient, p e c
remainder)
      where
        delta :: e
delta = p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
remainder e -> e -> e
forall a. Num a => a -> a -> a
- p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
q
        t :: p e c
t = c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e c
remainder c -> c -> c
forall a. Fractional a => a -> a -> a
/ p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e c
q) (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ e -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power e
delta
        -- remainder and q * t will have the same leading coefficients.
        -- Subtract them without the leading terms.
        -- Not necessary for purely numeric coefficients,
        -- but guarantees the cancellation of the leading terms when coefficients are symbolic.
        remainder' :: p e c
remainder' = p e c -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> p e c
deleteLeadingTerm p e c
remainder
        qt' :: p e c
qt' = p e c -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> p e c
deleteLeadingTerm (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ p e c
q p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
t

-- | Polynomial pseudo-division.  It returns the pseudo-quotient and pseudo-remainder polynomials.
--
-- Equivalent to \(b^{\delta+1} p\) divided by \(q\),
-- where \(p\) and \(q\) are polynomials with integer coefficients,
-- \(b\) is the leading coefficient of \(q\) and \(\delta=\max(-1, \deg(p) - \deg(q))\).
-- This guarantees the pseudo-quotient and pseudo-remainder exist,
-- even when the quotient and remainder do not when only integer coefficients are allowed.
--
-- For example, with \(p = 3x^3 + x^2 + x + 5\) and \(q = 5x^2 - 3x + 1\),
-- it is the case that \(5^2p = (15x + 14)q + (52x + 111)\):
--
-- >>> let p = 3 * power 3 + power 2 + power 1 + 5 :: IndexedPolynomial
-- >>> let q = 5 * power 2 - 3 * power 1 + 1 :: IndexedPolynomial
-- >>> pseudoDivide p q
-- (15x + 14,52x + 111)
pseudoDivide ::
  (Polynomial p e c, Eq (p e c), Num (p e c), Num c) =>
  -- | Dividend polynomial being pseudo-divided.
  p e c ->
  -- | Divisor polynomial pseudo-dividing the dividend.
  p e c ->
  -- | Pseudo-quotient and pseudo-remainder.
  (p e c, p e c)
pseudoDivide :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Num c) =>
p e c -> p e c -> (p e c, p e c)
pseudoDivide p e c
p p e c
q
  | p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
p e -> e -> Bool
forall a. Ord a => a -> a -> Bool
< p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
q = (p e c
0, p e c
p)
  | Bool
otherwise = e -> p e c -> p e c -> (p e c, p e c)
forall {b}. Integral b => b -> p e c -> p e c -> (p e c, p e c)
go (e
1 e -> e -> e
forall a. Num a => a -> a -> a
+ p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
p e -> e -> e
forall a. Num a => a -> a -> a
- p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
q) p e c
0 p e c
p
  where
    b :: c
b = p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e c
q
    go :: b -> p e c -> p e c -> (p e c, p e c)
go b
n p e c
quotient p e c
remainder
      | p e c
remainder p e c -> p e c -> Bool
forall a. Eq a => a -> a -> Bool
/= p e c
0, e
delta e -> e -> Bool
forall a. Ord a => a -> a -> Bool
>= e
0 = b -> p e c -> p e c -> (p e c, p e c)
go (b
n b -> b -> b
forall a. Num a => a -> a -> a
- b
1) p e c
quotient' p e c
remainder'
      | Bool
otherwise = (c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (c
b c -> b -> c
forall a b. (Num a, Integral b) => a -> b -> a
^ b
n) p e c
quotient, c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (c
b c -> b -> c
forall a b. (Num a, Integral b) => a -> b -> a
^ b
n) p e c
remainder)
      where
        delta :: e
delta = p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
remainder e -> e -> e
forall a. Num a => a -> a -> a
- p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
q
        t :: p e c
t = c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e c
remainder) (e -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power e
delta)
        quotient' :: p e c
quotient' = c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale c
b p e c
quotient p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
+ p e c
t
        -- Subtract with the leading terms deleted.
        -- The leading terms cancel out numerically,
        -- but guarantee cancellation when the coefficients are symbolic.
        remainder' :: p e c
remainder' = p e c -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> p e c
deleteLeadingTerm (c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale c
b p e c
remainder) p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
- p e c -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> p e c
deleteLeadingTerm (p e c
t p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
q)

-- | The extended Euclidean algorithm.  For polynomials \(p\) and \(q\),
-- it returns the greatest common divisor between \(p\) and \(q\).
-- It also returns \(s\) and \(t\) such that \(sp+tq = \gcd(p,q)\).
--
-- For example, for \(p=2x^5-2x\) and \(q=x^4-2x^2+1\), it is the case
-- that \(\gcd(p,q)=-x^2+1\) and \((-\frac{1}{4}x) p + (\frac{1}{2}x^2 + 1) q = -x^2+1\):
--
-- >>> let p = 2 * power 5 - 2 * power 1 :: IndexedPolynomial
-- >>> let q = power 4 - 2 * power 2 + 1 :: IndexedPolynomial
-- >>> extendedEuclidean p q
-- (((-1) % 4)x,(1 % 2)x^2 + 1,(-1)x^2 + 1)
extendedEuclidean ::
  (Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
  -- | Polynomial \(p\).
  p e c ->
  -- | Polynomial \(q\).
  p e c ->
  -- | \(s\), \(t\), and \(\gcd(p,q)\).
  (p e c, p e c, p e c)
extendedEuclidean :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c, p e c)
extendedEuclidean p e c
u p e c
v = p e c
-> p e c
-> p e c
-> p e c
-> p e c
-> p e c
-> (p e c, p e c, p e c)
forall {p :: * -> * -> *} {e} {c}.
(Polynomial p e c, Fractional c, Eq (p e c), Num (p e c)) =>
p e c
-> p e c
-> p e c
-> p e c
-> p e c
-> p e c
-> (p e c, p e c, p e c)
descend p e c
u p e c
v p e c
1 p e c
0 p e c
0 p e c
1
  where
    descend :: p e c
-> p e c
-> p e c
-> p e c
-> p e c
-> p e c
-> (p e c, p e c, p e c)
descend p e c
g p e c
0 p e c
s p e c
t p e c
_ p e c
_ = (p e c
s, p e c
t, p e c
g)
    descend p e c
a p e c
b p e c
a1 p e c
a2 p e c
b1 p e c
b2 = p e c
-> p e c
-> p e c
-> p e c
-> p e c
-> p e c
-> (p e c, p e c, p e c)
descend p e c
b p e c
r p e c
b1 p e c
b2 p e c
r1 p e c
r2
      where
        (p e c
q, p e c
r) = p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
divide p e c
a p e c
b
        r1 :: p e c
r1 = p e c
a1 p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
- p e c
q p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
b1
        r2 :: p e c
r2 = p e c
a2 p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
- p e c
q p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
b2

-- | Solves \(sa + tb = c\) for given polynomials \(a\), \(b\), and \(c\).
-- It will be the case that either \(s=0\) or
-- the degree of \(s\) will be less than the degree of \(b\).
--
-- >>> let a = power 4 - 2 * power 3 - 6 * power 2 + 12 * power 1 + 15 :: IndexedPolynomial
-- >>> let b = power 3 + power 2 - 4 * power 1 - 4 :: IndexedPolynomial
-- >>> let c = power 2 - 1 :: IndexedPolynomial
-- >>> diophantineEuclidean a b c
-- Just (((-1) % 5)x^2 + (4 % 5)x + ((-3) % 5),(1 % 5)x^3 + ((-7) % 5)x^2 + (16 % 5)x + (-2))
--
-- If there is no such \((s,t)\), then 'Nothing' is returned.
diophantineEuclidean ::
  (Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
  -- | Polynomial \(a\).
  p e c ->
  -- | Polynomial \(b\).
  p e c ->
  -- | Polynomial \(c\).
  p e c ->
  -- | \((s,t)\) such that \(sa + tb = c\).
  Maybe (p e c, p e c)
diophantineEuclidean :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> p e c -> Maybe (p e c, p e c)
diophantineEuclidean p e c
a p e c
b p e c
c
  | p e c
r p e c -> p e c -> Bool
forall a. Eq a => a -> a -> Bool
/= p e c
0 = Maybe (p e c, p e c)
forall a. Maybe a
Nothing
  | p e c
s' p e c -> p e c -> Bool
forall a. Eq a => a -> a -> Bool
/= p e c
0, p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
s' e -> e -> Bool
forall a. Ord a => a -> a -> Bool
>= p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
b = (p e c, p e c) -> Maybe (p e c, p e c)
forall a. a -> Maybe a
Just (p e c
r', p e c
t' p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
+ p e c
q' p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
a)
  | Bool
otherwise = (p e c, p e c) -> Maybe (p e c, p e c)
forall a. a -> Maybe a
Just (p e c
s', p e c
t')
  where
    (p e c
s, p e c
t, p e c
g) = p e c -> p e c -> (p e c, p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c, p e c)
extendedEuclidean p e c
a p e c
b
    (p e c
q, p e c
r) = p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
divide p e c
c p e c
g
    s' :: p e c
s' = p e c
q p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
s
    t' :: p e c
t' = p e c
q p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
* p e c
t
    (p e c
q', p e c
r') = p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
divide p e c
s' p e c
b

-- | Returns the greatest common divisor btween two polynomials.
--
-- Convenient wrapper over 'extendedEuclidean' which only returns the greatest common divisor.
greatestCommonDivisor ::
  (Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
  -- | Polynomial \(p\).
  p e c ->
  -- | Polynomial \(q\).
  p e c ->
  -- | \(\gcd(p,q)\).
  p e c
greatestCommonDivisor :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> p e c
greatestCommonDivisor p e c
p p e c
q = p e c
g
  where
    (p e c
_, p e c
_, p e c
g) = p e c -> p e c -> (p e c, p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c, p e c)
extendedEuclidean p e c
p p e c
q

-- | Returns the resultant and the subresultant polynomial remainder sequence for the given polynomials.
--
-- >>> subresultant (power 2 + 1) (power 2 - 1 :: IndexedPolynomial)
-- (4 % 1,[x^2 + 1,x^2 + (-1),(-2),0])
-- >>> subresultant (2 * power 2 - 3 * power 1 + 1) (5 * power 2 + power 1 - 6 :: IndexedPolynomial)
-- (0 % 1,[2x^2 + (-3)x + 1,5x^2 + x + (-6),17x + (-17),0])
-- >>> subresultant (power 3 + 2 * power 2 + 3 * power 1 + 4) (5 * power 2 + 6 * power 1 + 7 :: IndexedPolynomial)
-- (832 % 1,[x^3 + 2x^2 + 3x + 4,5x^2 + 6x + 7,16x + 72,832,0])
--
-- === __Reference__
--
-- See sections 1.4 and 1.5 in
-- [/Symbolic Integration I: Transcendental Functions/](https://doi.org/10.1007/b138171)
-- by Manuel Bronstein for the definition of resultants, subresultants,
-- polynomial remainder sequences, and subresultant polynomial remainder sequences.
subresultant ::
  (Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>
  -- | First element in the remainder sequence.
  p e c ->
  -- | Second element in the remainder sequence.
  p e c ->
  -- | The resultant and the subresultant polynomial remainder sequence.
  (c, [p e c])
subresultant :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>
p e c -> p e c -> (c, [p e c])
subresultant p e c
p p e c
q
  | p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
p e -> e -> Bool
forall a. Ord a => a -> a -> Bool
>= p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
q = ([p e c] -> [c] -> c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>
[p e c] -> [c] -> c
resultantFromSequence [p e c]
rs [c]
betas, [p e c]
rs)
  | Bool
otherwise = ((-c
1) c -> e -> c
forall a b. (Num a, Integral b) => a -> b -> a
^ (p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
q e -> e -> e
forall a. Num a => a -> a -> a
* p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
p) c -> c -> c
forall a. Num a => a -> a -> a
* c
resultant, [p e c]
prs)
  where
    ([p e c]
rs, [c]
betas) = (p e c, p e c) -> c -> c -> ([p e c], [c])
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>
(p e c, p e c) -> c -> c -> ([p e c], [c])
subresultantRemainderSequence (p e c
p, p e c
q) c
gamma c
beta
    gamma :: c
gamma = -c
1
    beta :: c
beta = (-c
1) c -> e -> c
forall a b. (Num a, Integral b) => a -> b -> a
^ (e
1 e -> e -> e
forall a. Num a => a -> a -> a
+ e
delta)
    delta :: e
delta = p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
p e -> e -> e
forall a. Num a => a -> a -> a
- p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
q

    (c
resultant, [p e c]
prs) = p e c -> p e c -> (c, [p e c])
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>
p e c -> p e c -> (c, [p e c])
subresultant p e c
q p e c
p

-- | Derives the subresultant polynomial remainder sequence for 'subresultant'.
-- Constructs \(\gamma_i\), \(\beta_i\), and the remainder sequence as it goes along.
-- Returns the remainder sequence and the sequence of \(\beta_i\).
subresultantRemainderSequence ::
  (Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>
  -- | The previous and current remainders in the sequence.
  (p e c, p e c) ->
  -- | \(\gamma_i\) as defined for the subresultant PRS.
  c ->
  -- | \(\beta_i\) as defined for the subresultant PRS.
  c ->
  -- | Polynomial remainder sequence and sequence of \(\beta_i\).
  ([p e c], [c])
subresultantRemainderSequence :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>
(p e c, p e c) -> c -> c -> ([p e c], [c])
subresultantRemainderSequence (p e c
rprev, p e c
rcurr) c
gamma c
beta
  | p e c
rcurr p e c -> p e c -> Bool
forall a. Eq a => a -> a -> Bool
/= p e c
0 = (p e c
rprev p e c -> [p e c] -> [p e c]
forall a. a -> [a] -> [a]
: [p e c]
rs, c
beta c -> [c] -> [c]
forall a. a -> [a] -> [a]
: [c]
betas)
  | Bool
otherwise = ([p e c
rprev, p e c
rcurr], [c
beta])
  where
    ([p e c]
rs, [c]
betas) = (p e c, p e c) -> c -> c -> ([p e c], [c])
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>
(p e c, p e c) -> c -> c -> ([p e c], [c])
subresultantRemainderSequence (p e c
rcurr, p e c
rnext) c
gamma' c
beta'
    (p e c
_, p e c
r) = p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Num c) =>
p e c -> p e c -> (p e c, p e c)
pseudoDivide p e c
rprev p e c
rcurr
    rnext :: p e c
rnext = c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (c
1 c -> c -> c
forall a. Fractional a => a -> a -> a
/ c
beta) p e c
r
    lc :: c
lc = p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e c
rcurr
    delta :: e
delta = p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
rprev e -> e -> e
forall a. Num a => a -> a -> a
- p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
rcurr
    delta' :: e
delta' = p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
rcurr e -> e -> e
forall a. Num a => a -> a -> a
- p e c -> e
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p e c
rnext
    gamma' :: c
gamma' = ((-c
lc) c -> e -> c
forall a b. (Num a, Integral b) => a -> b -> a
^ e
delta) c -> c -> c
forall a. Num a => a -> a -> a
* (c
gamma c -> e -> c
forall a b. (Fractional a, Integral b) => a -> b -> a
^^ (e
1 e -> e -> e
forall a. Num a => a -> a -> a
- e
delta))
    beta' :: c
beta' = (-c
lc) c -> c -> c
forall a. Num a => a -> a -> a
* (c
gamma' c -> e -> c
forall a b. (Num a, Integral b) => a -> b -> a
^ e
delta')

-- | Constructs the resultant based on the subresultant polynomial remainder sequence
-- and the sequence of \(\beta_i\) used to construct the subresultant PRS.
resultantFromSequence ::
  (Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>
  -- | Subresultant polynomial remainder sequence.
  [p e c] ->
  -- | Sequence of \(\beta_i\) used for deriving the subresultant PRS.
  [c] ->
  -- | Resultant.
  c
resultantFromSequence :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>
[p e c] -> [c] -> c
resultantFromSequence [p e c]
rs [c]
betas = [p e c] -> [c] -> c -> c -> c
forall {t} {b} {p :: * -> * -> *}.
(Fractional t, Polynomial p b t) =>
[p b t] -> [t] -> t -> t -> t
go [p e c]
rs [c]
betas c
1 c
1
  where
    go :: [p b t] -> [t] -> t -> t -> t
go (p b t
r : p b t
r' : p b t
r'' : [p b t]
rs') (t
beta : [t]
betas') t
c t
s
      | [] <- [p b t]
rs', p b t -> b
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p b t
r' b -> b -> Bool
forall a. Ord a => a -> a -> Bool
> b
0 = t
0
      | [] <- [p b t]
rs', p b t -> b
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p b t
r b -> b -> Bool
forall a. Eq a => a -> a -> Bool
== b
1 = p b t -> t
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p b t
r'
      | [] <- [p b t]
rs' = t
s t -> t -> t
forall a. Num a => a -> a -> a
* t
c t -> t -> t
forall a. Num a => a -> a -> a
* p b t -> t
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p b t
r' t -> b -> t
forall a b. (Num a, Integral b) => a -> b -> a
^ p b t -> b
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p b t
r
      | Bool
otherwise = [p b t] -> [t] -> t -> t -> t
go (p b t
r' p b t -> [p b t] -> [p b t]
forall a. a -> [a] -> [a]
: p b t
r'' p b t -> [p b t] -> [p b t]
forall a. a -> [a] -> [a]
: [p b t]
rs') [t]
betas' t
c' t
s'
      where
        s' :: t
s' | b -> Bool
forall a. Integral a => a -> Bool
odd (p b t -> b
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p b t
r), b -> Bool
forall a. Integral a => a -> Bool
odd (p b t -> b
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p b t
r') = -t
s | Bool
otherwise = t
s
        c' :: t
c' = t
c t -> t -> t
forall a. Num a => a -> a -> a
* ((t
beta t -> t -> t
forall a. Fractional a => a -> a -> a
/ (t
lc t -> b -> t
forall a b. (Num a, Integral b) => a -> b -> a
^ (b
1 b -> b -> b
forall a. Num a => a -> a -> a
+ b
delta))) t -> b -> t
forall a b. (Num a, Integral b) => a -> b -> a
^ p b t -> b
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p b t
r') t -> t -> t
forall a. Num a => a -> a -> a
* (t
lc t -> b -> t
forall a b. (Num a, Integral b) => a -> b -> a
^ (p b t -> b
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p b t
r b -> b -> b
forall a. Num a => a -> a -> a
- p b t -> b
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p b t
r''))
        lc :: t
lc = p b t -> t
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p b t
r'
        delta :: b
delta = p b t -> b
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p b t
r b -> b -> b
forall a. Num a => a -> a -> a
- p b t -> b
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> e
degree p b t
r'
    go [p b t]
_ [t]
_ t
_ t
_ = t
0

-- | Returns the derivative of the given polynomial.
--
-- >>> differentiate (power 2 + power 1 :: IndexedPolynomial)
-- 2x + 1
differentiate :: (Polynomial p e c, Num (p e c), Num c) => p e c -> p e c
differentiate :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate p e c
p = Sum (p e c) -> p e c
forall a. Sum a -> a
getSum (Sum (p e c) -> p e c) -> Sum (p e c) -> p e c
forall a b. (a -> b) -> a -> b
$ (e -> c -> Sum (p e c)) -> p e c -> Sum (p e c)
forall m. Monoid m => (e -> c -> m) -> p e c -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms e -> c -> Sum (p e c)
forall {e} {p :: * -> * -> *} {c}.
(Polynomial p e c, Num (p e c)) =>
e -> c -> Sum (p e c)
diffTerm p e c
p
  where
    diffTerm :: e -> c -> Sum (p e c)
diffTerm e
0 c
_ = p e c -> Sum (p e c)
forall a. a -> Sum a
Sum p e c
0
    diffTerm e
e c
c = p e c -> Sum (p e c)
forall a. a -> Sum a
Sum (p e c -> Sum (p e c)) -> p e c -> Sum (p e c)
forall a b. (a -> b) -> a -> b
$ c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (e -> c
forall a b. (Integral a, Num b) => a -> b
fromIntegral e
e c -> c -> c
forall a. Num a => a -> a -> a
* c
c) (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ e -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power (e
e e -> e -> e
forall a. Num a => a -> a -> a
- e
1)

-- | Returns the integral of the given polynomial.
--
-- >>> integrate (power 2 + power 1 :: IndexedPolynomial)
-- (1 % 3)x^3 + (1 % 2)x^2
integrate :: (Polynomial p e c, Num (p e c), Fractional c) => p e c -> p e c
integrate :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Fractional c) =>
p e c -> p e c
integrate p e c
p = Sum (p e c) -> p e c
forall a. Sum a -> a
getSum (Sum (p e c) -> p e c) -> Sum (p e c) -> p e c
forall a b. (a -> b) -> a -> b
$ (e -> c -> Sum (p e c)) -> p e c -> Sum (p e c)
forall m. Monoid m => (e -> c -> m) -> p e c -> m
forall (p :: * -> * -> *) e c m.
(Polynomial p e c, Monoid m) =>
(e -> c -> m) -> p e c -> m
foldTerms e -> c -> Sum (p e c)
forall {p :: * -> * -> *} {e} {c}.
(Polynomial p e c, Fractional c) =>
e -> c -> Sum (p e c)
integrateTerm p e c
p
  where
    integrateTerm :: e -> c -> Sum (p e c)
integrateTerm e
e c
c = p e c -> Sum (p e c)
forall a. a -> Sum a
Sum (p e c -> Sum (p e c)) -> p e c -> Sum (p e c)
forall a b. (a -> b) -> a -> b
$ c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (c
c c -> c -> c
forall a. Fractional a => a -> a -> a
/ (c
1 c -> c -> c
forall a. Num a => a -> a -> a
+ e -> c
forall a b. (Integral a, Num b) => a -> b
fromIntegral e
e)) (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ e -> p e c
forall (p :: * -> * -> *) e c. Polynomial p e c => e -> p e c
power (e
e e -> e -> e
forall a. Num a => a -> a -> a
+ e
1)

-- | Returns the squarefree factorization of the given polynomial.
--
-- Specifically, for a polynomial \(p\), find \([p_1, p_2, \ldots, p_n]\) such that
--
-- \[ p = \sum_{k=1}^n p_k^k \]
--
-- where all \(p_k\) are squarefree, i.e., there is no polynomial \(q\) such that \(q^2 = p_k\).
--
-- For example, the squarefree factorization of \(x^8 + 6x^6 + 12x^4 + 8x^2\)
-- is \(x^2 (x^2 + 2)^3\):
--
-- >>> squarefree (power 8 + 6 * power 6 + 12 * power 4 + 8 * power 2 :: IndexedPolynomial)
-- [1,x,x^2 + 2]
squarefree :: (Polynomial p e c, Eq (p e c), Num (p e c), Eq c, Fractional c) => p e c -> [p e c]
squarefree :: forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Eq c, Fractional c) =>
p e c -> [p e c]
squarefree p e c
0 = [p e c
0]
squarefree p e c
p
  | (p e c
x : [p e c]
xs) <- p e c -> p e c -> [p e c]
forall {p :: * -> * -> *} {e} {c}.
(Polynomial p e c, Fractional c, Eq c, Eq (p e c), Num (p e c)) =>
p e c -> p e c -> [p e c]
factor p e c
u p e c
v = c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale c
c p e c
x p e c -> [p e c] -> [p e c]
forall a. a -> [a] -> [a]
: [p e c]
xs
  | Bool
otherwise = [c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale c
c p e c
1]
  where
    c :: c
c = p e c -> c
forall (p :: * -> * -> *) e c. Polynomial p e c => p e c -> c
leadingCoefficient p e c
p
    q :: p e c
q = c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
Polynomial p e c =>
c -> p e c -> p e c
scale (c
1 c -> c -> c
forall a. Fractional a => a -> a -> a
/ c
c) p e c
p
    q' :: p e c
q' = p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate p e c
q
    g :: p e c
g = p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ p e c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> p e c
greatestCommonDivisor p e c
q p e c
q'
    (p e c
u, p e c
_) = p e c
q p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
`divide` p e c
g
    (p e c
v, p e c
_) = p e c
q' p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
`divide` p e c
g
    factor :: p e c -> p e c -> [p e c]
factor p e c
s p e c
y
      | p e c
z p e c -> p e c -> Bool
forall a. Eq a => a -> a -> Bool
== p e c
0 = [p e c
s]
      | Bool
otherwise = p e c
f p e c -> [p e c] -> [p e c]
forall a. a -> [a] -> [a]
: p e c -> p e c -> [p e c]
factor p e c
s' p e c
y'
      where
        z :: p e c
z = p e c
y p e c -> p e c -> p e c
forall a. Num a => a -> a -> a
- p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Num (p e c), Num c) =>
p e c -> p e c
differentiate p e c
s
        f :: p e c
f = p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq c, Fractional c) =>
p e c -> p e c
monic (p e c -> p e c) -> p e c -> p e c
forall a b. (a -> b) -> a -> b
$ p e c -> p e c -> p e c
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> p e c
greatestCommonDivisor p e c
s p e c
z
        (p e c
s', p e c
_) = p e c
s p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
`divide` p e c
f
        (p e c
y', p e c
_) = p e c
z p e c -> p e c -> (p e c, p e c)
forall (p :: * -> * -> *) e c.
(Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>
p e c -> p e c -> (p e c, p e c)
`divide` p e c
f