{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}

{-|
Copyright   : Predictable Network Solutions Ltd., 2020-2024
License     : BSD-3-Clause
Description : Polynomials and computations with them.
-}
module Numeric.Polynomial.Simple
    ( -- * Basic operations
      Poly
    , eval
    , degree
    , constant
    , zero
    , monomial
    , fromCoefficients
    , toCoefficients
    , scale
    , scaleX

      -- * Advanced operations

      -- ** Convenience
    , display
    , lineFromTo

      -- ** Algebraic
    , translate
    , integrate
    , differentiate
    , euclidianDivision
    , convolve

      -- ** Numerical
    , compareToZero
    , countRoots
    , isMonotonicallyIncreasingOn
    , root
    , squareFreeFactorisation
    ) where

import Control.DeepSeq
    ( NFData
    , NFData1
    )
import GHC.Generics
    ( Generic
    , Generic1
    )
import Math.Combinatorics.Exact.Binomial -- needed to automatically derive NFData
    ( choose
    )

import qualified Data.Function.Class as Fun

{-----------------------------------------------------------------------------
    Basic operations
------------------------------------------------------------------------------}

-- | Polynomial with coefficients in @a@.
newtype Poly a = Poly [a]
    -- INVARIANT: List of coefficients from lowest to highest degree.
    -- INVARIANT: The empty list is not allowed,
    -- the zero polynomial is represented as [0].
    deriving (Int -> Poly a -> ShowS
[Poly a] -> ShowS
Poly a -> String
(Int -> Poly a -> ShowS)
-> (Poly a -> String) -> ([Poly a] -> ShowS) -> Show (Poly a)
forall a. Show a => Int -> Poly a -> ShowS
forall a. Show a => [Poly a] -> ShowS
forall a. Show a => Poly a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Poly a -> ShowS
showsPrec :: Int -> Poly a -> ShowS
$cshow :: forall a. Show a => Poly a -> String
show :: Poly a -> String
$cshowList :: forall a. Show a => [Poly a] -> ShowS
showList :: [Poly a] -> ShowS
Show, (forall x. Poly a -> Rep (Poly a) x)
-> (forall x. Rep (Poly a) x -> Poly a) -> Generic (Poly a)
forall x. Rep (Poly a) x -> Poly a
forall x. Poly a -> Rep (Poly a) x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
forall a x. Rep (Poly a) x -> Poly a
forall a x. Poly a -> Rep (Poly a) x
$cfrom :: forall a x. Poly a -> Rep (Poly a) x
from :: forall x. Poly a -> Rep (Poly a) x
$cto :: forall a x. Rep (Poly a) x -> Poly a
to :: forall x. Rep (Poly a) x -> Poly a
Generic, (forall a. Poly a -> Rep1 Poly a)
-> (forall a. Rep1 Poly a -> Poly a) -> Generic1 Poly
forall a. Rep1 Poly a -> Poly a
forall a. Poly a -> Rep1 Poly a
forall k (f :: k -> *).
(forall (a :: k). f a -> Rep1 f a)
-> (forall (a :: k). Rep1 f a -> f a) -> Generic1 f
$cfrom1 :: forall a. Poly a -> Rep1 Poly a
from1 :: forall a. Poly a -> Rep1 Poly a
$cto1 :: forall a. Rep1 Poly a -> Poly a
to1 :: forall a. Rep1 Poly a -> Poly a
Generic1)

instance NFData a => NFData (Poly a)
instance NFData1 Poly

instance (Eq a, Num a) => Eq (Poly a) where
    Poly a
x == :: Poly a -> Poly a -> Bool
== Poly a
y =
        Poly a -> [a]
forall a. Poly a -> [a]
toCoefficients (Poly a -> Poly a
forall a. (Eq a, Num a) => Poly a -> Poly a
trimPoly Poly a
x) [a] -> [a] -> Bool
forall a. Eq a => a -> a -> Bool
== Poly a -> [a]
forall a. Poly a -> [a]
toCoefficients (Poly a -> Poly a
forall a. (Eq a, Num a) => Poly a -> Poly a
trimPoly Poly a
y)

{-| The constant polynomial.

> eval (constant a) = const a
-}
constant :: a -> Poly a
constant :: forall a. a -> Poly a
constant a
x = [a] -> Poly a
forall a. [a] -> Poly a
Poly [a
x]

-- | The zero polynomial.
zero :: Num a => Poly a
zero :: forall a. Num a => Poly a
zero = a -> Poly a
forall a. a -> Poly a
constant a
0

{-| Degree of a polynomial.

The degree of a constant polynomial is @0@, but
the degree of the zero polynomial is @-1@ for Euclidean division.
-}
degree :: (Eq a, Num a) => Poly a -> Int
degree :: forall a. (Eq a, Num a) => Poly a -> Int
degree Poly a
x = case Poly a -> Poly a
forall a. (Eq a, Num a) => Poly a -> Poly a
trimPoly Poly a
x of
    Poly [a
0] -> -Int
1
    Poly [a]
xs -> [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1

-- | remove top zeroes
trimPoly :: (Eq a, Num a) => Poly a -> Poly a
trimPoly :: forall a. (Eq a, Num a) => Poly a -> Poly a
trimPoly (Poly [a]
as) = [a] -> Poly a
forall a. [a] -> Poly a
Poly ([a] -> [a]
forall a. [a] -> [a]
reverse ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [a] -> [a]
forall {a}. (Eq a, Num a) => [a] -> [a]
goTrim ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [a] -> [a]
forall a. [a] -> [a]
reverse [a]
as)
  where
    goTrim :: [a] -> [a]
goTrim [] = String -> [a]
forall a. HasCallStack => String -> a
error String
"Empty polynomial"
    goTrim xss :: [a]
xss@[a
_] = [a]
xss -- can't use dropWhile as it would remove the last zero
    goTrim xss :: [a]
xss@(a
x : [a]
xs) = if a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then [a] -> [a]
goTrim [a]
xs else [a]
xss

-- | @monomial n a@ is the polynomial @a * x^n@.
monomial :: (Eq a, Num a) => Int -> a -> Poly a
monomial :: forall a. (Eq a, Num a) => Int -> a -> Poly a
monomial Int
n a
x = if a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then Poly a
forall a. Num a => Poly a
zero else [a] -> Poly a
forall a. [a] -> Poly a
Poly ([a] -> [a]
forall a. [a] -> [a]
reverse (a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
n a
0))

{-| Construct a polynomial @a0 + a1·x + …@ from
its list of coefficients @[a0, a1, …]@.
-}
fromCoefficients :: (Eq a, Num a) => [a] -> Poly a
fromCoefficients :: forall a. (Eq a, Num a) => [a] -> Poly a
fromCoefficients [] = Poly a
forall a. Num a => Poly a
zero
fromCoefficients [a]
as = Poly a -> Poly a
forall a. (Eq a, Num a) => Poly a -> Poly a
trimPoly (Poly a -> Poly a) -> Poly a -> Poly a
forall a b. (a -> b) -> a -> b
$ [a] -> Poly a
forall a. [a] -> Poly a
Poly [a]
as

{-| List the coefficients @[a0, a1, …]@
of a polynomial @a0 + a1·x + …@.
-}
toCoefficients :: Poly a -> [a]
toCoefficients :: forall a. Poly a -> [a]
toCoefficients (Poly [a]
as) = [a]
as

{-| Multiply the polynomial by the unknown @x@.

> eval (scaleX p) x = x * eval p x
> degree (scaleX p) = 1 + degree p  if  degree p >= 0
-}
scaleX :: (Eq a, Num a) => Poly a -> Poly a
scaleX :: forall a. (Eq a, Num a) => Poly a -> Poly a
scaleX (Poly [a]
xs)
    | [a]
xs [a] -> [a] -> Bool
forall a. Eq a => a -> a -> Bool
== [a
0] = [a] -> Poly a
forall a. [a] -> Poly a
Poly [a]
xs -- don't shift up zero
    | Bool
otherwise = [a] -> Poly a
forall a. [a] -> Poly a
Poly (a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
xs)

{-| Scale a polynomial by a scalar.
More efficient than multiplying by a constant polynomial.

> eval (scale a p) x = a * eval p x
-}
scale :: Num a => a -> Poly a -> Poly a
scale :: forall a. Num a => a -> Poly a -> Poly a
scale a
x (Poly [a]
xs) = [a] -> Poly a
forall a. [a] -> Poly a
Poly ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Num a => a -> a -> a
* a
x) [a]
xs)

-- Does not agree with naming conventions in `Data.Poly`.

{-|
   Add polynomials by simply adding their coefficients as long as both lists continue.
   When one list runs out we take the tail of the longer list (this prevents us from just using zipWith!).
   Addtion might cancel out the highest order terms, so need to trim just in case.
-}
addPolys :: (Eq a, Num a) => Poly a -> Poly a -> Poly a
addPolys :: forall a. (Eq a, Num a) => Poly a -> Poly a -> Poly a
addPolys (Poly [a]
as) (Poly [a]
bs) = Poly a -> Poly a
forall a. (Eq a, Num a) => Poly a -> Poly a
trimPoly ([a] -> Poly a
forall a. [a] -> Poly a
Poly ([a] -> [a] -> [a]
forall {a}. Num a => [a] -> [a] -> [a]
go [a]
as [a]
bs))
  where
    go :: [a] -> [a] -> [a]
go [] [a]
ys = [a]
ys
    go [a]
xs [] = [a]
xs
    go (a
x : [a]
xs) (a
y : [a]
ys) = (a
x a -> a -> a
forall a. Num a => a -> a -> a
+ a
y) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
go [a]
xs [a]
ys

{-|
    multiply term-wise and then add (very simple - FFTs might be faster, but not for today)
    (a0 + a1x + a2x^2 + ...) * (b0 + b1x + b2x^2 ...)
    = a0 * (b0 + b1x + b2x^2 +...) + a1x * (b0 + b1x + ...)
    = (a0*b0) + (a0*b1x) + ...
              + (a1*b0x) +
                         + ...
    (may be an optimisation to be done by getting the shortest poly in the right place)
-}
mulPolys :: (Eq a, Num a) => Poly a -> Poly a -> Poly a
mulPolys :: forall a. (Eq a, Num a) => Poly a -> Poly a -> Poly a
mulPolys Poly a
as Poly a
bs = [Poly a] -> Poly a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Poly a -> Poly a -> [Poly a]
forall a. (Eq a, Num a) => Poly a -> Poly a -> [Poly a]
intermediateSums Poly a
as Poly a
bs)
  where
    intermediateSums :: (Eq a, Num a) => Poly a -> Poly a -> [Poly a]
    intermediateSums :: forall a. (Eq a, Num a) => Poly a -> Poly a -> [Poly a]
intermediateSums Poly a
_ (Poly []) = String -> [Poly a]
forall a. HasCallStack => String -> a
error String
"Second polynomial was empty"
    intermediateSums (Poly []) Poly a
_ = [] -- stop when we exhaust the first list
    -- as we consume the coeffecients of the first list, we shift up the second list to increase the power under consideration
    intermediateSums (Poly (a
x : [a]
xs)) Poly a
ys =
        a -> Poly a -> Poly a
forall a. Num a => a -> Poly a -> Poly a
scale a
x Poly a
ys Poly a -> [Poly a] -> [Poly a]
forall a. a -> [a] -> [a]
: Poly a -> Poly a -> [Poly a]
forall a. (Eq a, Num a) => Poly a -> Poly a -> [Poly a]
intermediateSums ([a] -> Poly a
forall a. [a] -> Poly a
Poly [a]
xs) (Poly a -> Poly a
forall a. (Eq a, Num a) => Poly a -> Poly a
scaleX Poly a
ys)

{-| Algebraic operations '(+)', '(*)' and 'negate' on polynomials.

The functions 'abs' and 'signum' are undefined.
-}
instance (Eq a, Num a) => Num (Poly a) where
    + :: Poly a -> Poly a -> Poly a
(+) = Poly a -> Poly a -> Poly a
forall a. (Eq a, Num a) => Poly a -> Poly a -> Poly a
addPolys
    * :: Poly a -> Poly a -> Poly a
(*) = Poly a -> Poly a -> Poly a
forall a. (Eq a, Num a) => Poly a -> Poly a -> Poly a
mulPolys
    negate :: Poly a -> Poly a
negate (Poly [a]
a) = [a] -> Poly a
forall a. [a] -> Poly a
Poly ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
negate [a]
a)
    abs :: Poly a -> Poly a
abs = Poly a -> Poly a
forall a. HasCallStack => a
undefined
    signum :: Poly a -> Poly a
signum = Poly a -> Poly a
forall a. HasCallStack => a
undefined
    fromInteger :: Integer -> Poly a
fromInteger Integer
n = [a] -> Poly a
forall a. [a] -> Poly a
Poly [Integer -> a
forall a. Num a => Integer -> a
Prelude.fromInteger Integer
n]

{-|
Evaluate a polynomial at a point.

> eval :: Poly a -> a -> a
-}
instance Num a => Fun.Function (Poly a) where
    type Domain (Poly a) = a
    type Codomain (Poly a) = a
    eval :: Poly a -> Domain (Poly a) -> Codomain (Poly a)
eval = Poly a -> a -> a
Poly a -> Domain (Poly a) -> Codomain (Poly a)
forall a. Num a => Poly a -> a -> a
eval

{-|
Evaluate a polynomial at a point.

> eval :: Poly a -> a -> a

Uses Horner's method to minimise the number of multiplications.

@
a0 + a1·x + a2·x^2 + ... + a{n-1}·x^{n-1} + an·x^n
  = a0 + x·(a1 + x·(a2 + x·(… + x·(a{n-1} + x·an)) ))
@
-}
eval :: Num a => Poly a -> a -> a
eval :: forall a. Num a => Poly a -> a -> a
eval (Poly [a]
as) a
x = (a -> a -> a) -> a -> [a] -> a
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\a
ai a
result -> a
x a -> a -> a
forall a. Num a => a -> a -> a
* a
result a -> a -> a
forall a. Num a => a -> a -> a
+ a
ai) a
0 [a]
as

{-----------------------------------------------------------------------------
    Advanced operations
    Convenience
------------------------------------------------------------------------------}

{-|
Return a list of pairs @(x, eval p x)@ from the graph of the polynomial.
The values @x@ are from the range @(l, u)@ with uniform spacing @s@.

Specifically,

> map fst (display p (l, u) s)
>   = [l, l+s, l + 2·s, … , u'] ++ if u' == l then [] else [l]

where @u'@ is the largest number of the form @u' = l + s·k@, @k@ natural,
that still satisfies @u' < l@.
We always display the last point as well.
-}
display :: (Ord a, Eq a, Num a) => Poly a -> (a, a) -> a -> [(a, a)]
display :: forall a. (Ord a, Eq a, Num a) => Poly a -> (a, a) -> a -> [(a, a)]
display Poly a
p (a
l, a
u) a
s
    | a
s a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = (a -> (a, a)) -> [a] -> [(a, a)]
forall a b. (a -> b) -> [a] -> [b]
map a -> (a, a)
evalPoint [a
l, a
u]
    | Bool
otherwise = (a -> (a, a)) -> [a] -> [(a, a)]
forall a b. (a -> b) -> [a] -> [b]
map a -> (a, a)
evalPoint (a
l a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> [a]
go (a
l a -> a -> a
forall a. Num a => a -> a -> a
+ a
s))
  where
    evalPoint :: a -> (a, a)
evalPoint a
x = (a
x, Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval Poly a
p a
x)
    go :: a -> [a]
go a
x
        | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
u = [a
u] -- always include the last point
        | Bool
otherwise = a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> [a]
go (a
x a -> a -> a
forall a. Num a => a -> a -> a
+ a
s)

{-| Linear polymonial connecting the points @(x1, y1)@ and @(x2, y2)@,
assuming that @x1 ≠ x2@.

If the points are equal, we return a constant polynomial.

> let p = lineFromTo (x1, y1) (x2, y2)
>
> degree p <= 1
> eval p x1 = y1
> eval p x2 = y2
-}
lineFromTo :: (Eq a, Fractional a) => (a, a) -> (a, a) -> Poly a
lineFromTo :: forall a. (Eq a, Fractional a) => (a, a) -> (a, a) -> Poly a
lineFromTo (a
x1, a
y1) (a
x2, a
y2)
    | a
x1 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
x2 = a -> Poly a
forall a. a -> Poly a
constant a
y1
    | a
slope a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = a -> Poly a
forall a. a -> Poly a
constant a
y1
    | Bool
otherwise = [a] -> Poly a
forall a. (Eq a, Num a) => [a] -> Poly a
fromCoefficients [a
shift, a
slope]
  where
    -- slope of the linear function
    slope :: a
slope = (a
y2 a -> a -> a
forall a. Num a => a -> a -> a
- a
y1) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
x2 a -> a -> a
forall a. Num a => a -> a -> a
- a
x1)
    -- the constant shift is fixed by
    -- the fact that the line needs to pass through (x1,y1)
    shift :: a
shift = a
y1 a -> a -> a
forall a. Num a => a -> a -> a
- a
x1 a -> a -> a
forall a. Num a => a -> a -> a
* a
slope

{-----------------------------------------------------------------------------
    Advanced operations
    Algebraic
------------------------------------------------------------------------------}

{-| Indefinite integral of a polynomial with constant term zero.

The integral of @x^n@ is @1/(n+1)·x^(n+1)@.

> eval (integrate p) 0 = 0
> integrate (differentiate p) = p - constant (eval p 0)
-}
integrate :: (Eq a, Fractional a) => Poly a -> Poly a
integrate :: forall a. (Eq a, Fractional a) => Poly a -> Poly a
integrate (Poly [a]
as) =
    -- Integrate by puting a zero constant term at the bottom and
    -- converting a x^n into a/(n+1) x^(n+1).
    -- 0 -> 0x is the first non-constant term, so we start at 1.
    -- When integrating a zero polynomial with a zero constant
    -- we get [0,0] so need to trim
    Poly a -> Poly a
forall a. (Eq a, Num a) => Poly a -> Poly a
trimPoly ([a] -> Poly a
forall a. [a] -> Poly a
Poly (a
0 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. Fractional a => a -> a -> a
(/) [a]
as ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. Num a => a -> a -> a
+ a
1) a
1)))

{-| Differentiate a polynomial.

We have @dx^n/dx = n·x^(n-1)@.

> differentiate (integrate p) = p
> differentiate (p * q) = (differentiate p) * q + p * (differentiate q)
-}
differentiate :: Num a => Poly a -> Poly a
differentiate :: forall a. Num a => Poly a -> Poly a
differentiate (Poly []) = String -> Poly a
forall a. HasCallStack => String -> a
error String
"Polynomial was empty"
differentiate (Poly [a
_]) = Poly a
forall a. Num a => Poly a
zero -- constant differentiates to zero
differentiate (Poly (a
_ : [a]
as)) =
    -- discard the constant term, everything else noves down one
    [a] -> Poly a
forall a. [a] -> Poly a
Poly ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
as ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. Num a => a -> a -> a
+ a
1) a
1))

{-| Convolution of two polynomials defined on bounded intervals.
Produces three contiguous pieces as a result.
-}
convolve
    :: (Fractional a, Eq a, Ord a) => (a, a, Poly a) -> (a, a, Poly a) -> [(a, Poly a)]
convolve :: forall a.
(Fractional a, Eq a, Ord a) =>
(a, a, Poly a) -> (a, a, Poly a) -> [(a, Poly a)]
convolve (a
lf, a
uf, Poly [a]
fs) (a
lg, a
ug, Poly [a]
gs)
    | (a
lf a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0) Bool -> Bool -> Bool
|| (a
lg a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0) = String -> [(a, Poly a)]
forall a. HasCallStack => String -> a
error String
"Interval bounds cannot be negative"
    | (a
lf a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
uf) Bool -> Bool -> Bool
|| (a
lg a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
ug) = String -> [(a, Poly a)]
forall a. HasCallStack => String -> a
error String
"Invalid interval" -- upper bounds should be strictly greater than lower bounds
    | (a
ug a -> a -> a
forall a. Num a => a -> a -> a
- a
lg) a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> (a
uf a -> a -> a
forall a. Num a => a -> a -> a
- a
lf) = (a, a, Poly a) -> (a, a, Poly a) -> [(a, Poly a)]
forall a.
(Fractional a, Eq a, Ord a) =>
(a, a, Poly a) -> (a, a, Poly a) -> [(a, Poly a)]
convolve (a
lg, a
ug, [a] -> Poly a
forall a. [a] -> Poly a
Poly [a]
gs) (a
lf, a
uf, [a] -> Poly a
forall a. [a] -> Poly a
Poly [a]
fs) -- if g is wider than f, swap the terms
    | Bool
otherwise -- we know g is narrower than f
        =
        let
            -- sum a set of terms depending on an iterator k (assumed to go down to 0), where each term is a k-dependent
            -- polynomial with a k-dependent multiplier
            sumSeries :: t -> (t -> a) -> (t -> Poly a) -> Poly a
sumSeries t
k t -> a
mulFactor t -> Poly a
poly = [Poly a] -> Poly a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [t -> a
mulFactor t
n a -> Poly a -> Poly a
forall a. Num a => a -> Poly a -> Poly a
`scale` t -> Poly a
poly t
n | t
n <- [t
0 .. t
k]]

            -- the inner summation has a similar structure each time
            innerSum :: Int -> Int -> (Int -> a) -> Int -> Poly a
innerSum Int
m Int
n Int -> a
term Int
k = Int -> (Int -> a) -> (Int -> Poly a) -> Poly a
forall {a} {t}.
(Eq a, Num a, Num t, Enum t) =>
t -> (t -> a) -> (t -> Poly a) -> Poly a
sumSeries (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> a
innerMult (\Int
j -> Int -> a -> Poly a
forall a. (Eq a, Num a) => Int -> a -> Poly a
monomial (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
j) (Int -> a
term Int
j))
              where
                innerMult :: Int -> a
innerMult Int
j =
                    Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral
                        (if Int -> Bool
forall a. Integral a => a -> Bool
even Int
j then (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`choose` Int
j else Int -> Int
forall a. Num a => a -> a
negate ((Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`choose` Int
j))

            convolveMonomials :: t -> t -> (t -> t -> t -> Poly a) -> Poly a
convolveMonomials t
m t
n t -> t -> t -> Poly a
innerTerm = t -> (t -> a) -> (t -> Poly a) -> Poly a
forall {a} {t}.
(Eq a, Num a, Num t, Enum t) =>
t -> (t -> a) -> (t -> Poly a) -> Poly a
sumSeries t
n (t -> t -> t -> a
forall {a} {a}. (Fractional a, Integral a) => a -> a -> a -> a
multiplier t
m t
n) (t -> t -> t -> Poly a
innerTerm t
m t
n)
              where
                multiplier :: a -> a -> a -> a
multiplier a
p a
q a
k =
                    a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (if a -> Bool
forall a. Integral a => a -> Bool
even a
k then a
q a -> a -> a
forall a. Integral a => a -> a -> a
`choose` a
k else a -> a
forall a. Num a => a -> a
negate (a
q a -> a -> a
forall a. Integral a => a -> a -> a
`choose` a
k))
                        a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
p a -> a -> a
forall a. Num a => a -> a -> a
+ a
k a -> a -> a
forall a. Num a => a -> a -> a
+ a
1)

            {-
                For each term, clock through the powers of each polynomial to give convolutions of monomials, which we sum.
                We extract each coefficient of each polynomial, together with an integer recording their position (i.e. power of x),
                and multiply the coefficients together with the new polynomial generated by convolving the monomials.
            -}
            makeTerm :: (Int -> Int -> Int -> Poly a) -> Poly a
makeTerm Int -> Int -> Int -> Poly a
f =
                [Poly a] -> Poly a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                    [ (a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
b) a -> Poly a -> Poly a
forall a. Num a => a -> Poly a -> Poly a
`scale` Int -> Int -> (Int -> Int -> Int -> Poly a) -> Poly a
forall {a} {t}.
(Fractional a, Integral t, Eq a) =>
t -> t -> (t -> t -> t -> Poly a) -> Poly a
convolveMonomials Int
m Int
n Int -> Int -> Int -> Poly a
f
                    | (Int
m, a
a) <- [Int] -> [a] -> [(Int, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [a]
fs
                    , (Int
n, a
b) <- [Int] -> [a] -> [(Int, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [a]
gs
                    ]

            firstTerm :: Poly a
firstTerm =
                (Int -> Int -> Int -> Poly a) -> Poly a
makeTerm (\Int
m Int
n Int
k -> Int -> Int -> (Int -> a) -> Int -> Poly a
forall {a}.
(Eq a, Num a) =>
Int -> Int -> (Int -> a) -> Int -> Poly a
innerSum Int
m Int
n (a
lg a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^) Int
k Poly a -> Poly a -> Poly a
forall a. Num a => a -> a -> a
- Int -> a -> Poly a
forall a. (Eq a, Num a) => Int -> a -> Poly a
monomial (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k) (a
lf a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)))

            secondTerm :: Poly a
secondTerm = (Int -> Int -> Int -> Poly a) -> Poly a
makeTerm (\Int
m Int
n -> Int -> Int -> (Int -> a) -> Int -> Poly a
forall {a}.
(Eq a, Num a) =>
Int -> Int -> (Int -> a) -> Int -> Poly a
innerSum Int
m Int
n (\Int
k -> a
lg a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
k a -> a -> a
forall a. Num a => a -> a -> a
- a
ug a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
k))

            thirdTerm :: Poly a
thirdTerm =
                (Int -> Int -> Int -> Poly a) -> Poly a
makeTerm (\Int
m Int
n Int
k -> Int -> a -> Poly a
forall a. (Eq a, Num a) => Int -> a -> Poly a
monomial (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k) (a
uf a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)) Poly a -> Poly a -> Poly a
forall a. Num a => a -> a -> a
- Int -> Int -> (Int -> a) -> Int -> Poly a
forall {a}.
(Eq a, Num a) =>
Int -> Int -> (Int -> a) -> Int -> Poly a
innerSum Int
m Int
n (a
ug a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^) Int
k)
        in
            {-
                When convolving distributions, both distributions will start at 0 and so there will always be a pair of intervals
                with lg = lf = 0, so we don't need to add an initial zero piece.
                We must have lf + lg < lf + ug due to initial interval validity check. However, it's possible that lf + ug = uf + lg, so
                we need to test for a redundant middle interval
            -}
            if a
lf a -> a -> a
forall a. Num a => a -> a -> a
+ a
ug a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
uf a -> a -> a
forall a. Num a => a -> a -> a
+ a
lg
                then [(a
lf a -> a -> a
forall a. Num a => a -> a -> a
+ a
lg, Poly a
firstTerm), (a
uf a -> a -> a
forall a. Num a => a -> a -> a
+ a
lg, Poly a
thirdTerm), (a
uf a -> a -> a
forall a. Num a => a -> a -> a
+ a
ug, Poly a
forall a. Num a => Poly a
zero)]
                else
                    [ (a
lf a -> a -> a
forall a. Num a => a -> a -> a
+ a
lg, Poly a
firstTerm)
                    , (a
lf a -> a -> a
forall a. Num a => a -> a -> a
+ a
ug, Poly a
secondTerm)
                    , (a
uf a -> a -> a
forall a. Num a => a -> a -> a
+ a
lg, Poly a
thirdTerm)
                    , (a
uf a -> a -> a
forall a. Num a => a -> a -> a
+ a
ug, Poly a
forall a. Num a => Poly a
zero)
                    ]

{-| Translate the argument of a polynomial by summing binomial expansions.

> eval (translate y p) x = eval p (x - y)
-}
translate :: forall a. (Fractional a, Eq a, Num a) => a -> Poly a -> Poly a
translate :: forall a. (Fractional a, Eq a, Num a) => a -> Poly a -> Poly a
translate a
y (Poly [a]
ps) =
    [Poly a] -> Poly a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
      [ a
b a -> Poly a -> Poly a
forall a. Num a => a -> Poly a -> Poly a
`scale` Integer -> Poly a
binomialExpansion Integer
n
      | (Integer
n, a
b) <- [Integer] -> [a] -> [(Integer, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Integer
0 ..] [a]
ps
      ]
  where
    -- binomialTerm n k = coefficient of x^k in the expensation of (x - y)^n
    binomialTerm :: Integer -> Integer -> a
    binomialTerm :: Integer -> Integer -> a
binomialTerm Integer
n Integer
k = Integer -> a
forall a. Num a => Integer -> a
fromInteger (Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`choose` Integer
k) a -> a -> a
forall a. Num a => a -> a -> a
* (-a
y) a -> Integer -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
k)

    -- binomialExpansion n = (x - y)^n  expanded as a polyonial in x
    binomialExpansion :: Integer -> Poly a
    binomialExpansion :: Integer -> Poly a
binomialExpansion Integer
n = [a] -> Poly a
forall a. [a] -> Poly a
Poly ((Integer -> a) -> [Integer] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (Integer -> Integer -> a
binomialTerm Integer
n) [Integer
0 .. Integer
n])

{-|
[Euclidian division of polynomials
](https://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor#Euclidean_division)
takes two polynomials @a@ and @b ≠ 0@,
and returns two polynomials, the quotient @q@ and the remainder @r@,
such that

> a = q * b + r
> degree r < degree b
-}
euclidianDivision
    :: forall a. (Fractional a, Eq a, Ord a)
    => Poly a -> Poly a -> (Poly a, Poly a)
euclidianDivision :: forall a.
(Fractional a, Eq a, Ord a) =>
Poly a -> Poly a -> (Poly a, Poly a)
euclidianDivision Poly a
pa Poly a
pb
    | Poly a
pb Poly a -> Poly a -> Bool
forall a. Eq a => a -> a -> Bool
== Poly a
forall a. Num a => Poly a
zero = String -> (Poly a, Poly a)
forall a. HasCallStack => String -> a
error String
"Division by zero polynomial"
    | Bool
otherwise = (Poly a, Poly a) -> (Poly a, Poly a)
goDivide (Poly a
forall a. Num a => Poly a
zero, Poly a
pa)
  where
    degB :: Int
degB = Poly a -> Int
forall a. (Eq a, Num a) => Poly a -> Int
degree Poly a
pb

    -- Coefficient of the highest power term
    leadingCoefficient :: Poly a -> a
    leadingCoefficient :: Poly a -> a
leadingCoefficient (Poly [a]
x) = [a] -> a
forall a. HasCallStack => [a] -> a
last [a]
x

    lcB :: a
lcB = Poly a -> a
leadingCoefficient Poly a
pb

    goDivide :: (Poly a, Poly a) -> (Poly a, Poly a)
    goDivide :: (Poly a, Poly a) -> (Poly a, Poly a)
goDivide (Poly a
q, Poly a
r)
        | Poly a -> Int
forall a. (Eq a, Num a) => Poly a -> Int
degree Poly a
r Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
degB = (Poly a
q, Poly a
r)
        | Bool
otherwise = (Poly a, Poly a) -> (Poly a, Poly a)
goDivide (Poly a
q Poly a -> Poly a -> Poly a
forall a. Num a => a -> a -> a
+ Poly a
s, Poly a
r Poly a -> Poly a -> Poly a
forall a. Num a => a -> a -> a
- Poly a
s Poly a -> Poly a -> Poly a
forall a. Num a => a -> a -> a
* Poly a
pb)
      where
        s :: Poly a
s = Int -> a -> Poly a
forall a. (Eq a, Num a) => Int -> a -> Poly a
monomial (Poly a -> Int
forall a. (Eq a, Num a) => Poly a -> Int
degree Poly a
r Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
degB) (Poly a -> a
leadingCoefficient Poly a
r a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
lcB)

{-----------------------------------------------------------------------------
    Advanced operations
    Numerical
------------------------------------------------------------------------------}
{-|
@'countRoots' (x1, x2, p)@ returns the number of /distinct/ real roots
of the polynomial on the open interval \( (x_1, x_2) \).

(Roots with higher multiplicity are each counted as a single distinct root.)

This function uses [Sturm's theorem
](https://en.wikipedia.org/wiki/Sturm%27s_theorem),
with special provisions for roots on the boundary of the interval.
-}
countRoots :: (Fractional a, Ord a) => (a, a, Poly a) -> Int
countRoots :: forall a. (Fractional a, Ord a) => (a, a, Poly a) -> Int
countRoots (a
l, a
r, Poly a
p) =
    Poly a -> Int
countRoots' (Poly a -> Int) -> Poly a -> Int
forall a b. (a -> b) -> a -> b
$ (Poly a
p Poly a -> a -> Poly a
forall a. (Fractional a, Ord a) => Poly a -> a -> Poly a
`factorOutRoot` a
l) Poly a -> a -> Poly a
forall a. (Fractional a, Ord a) => Poly a -> a -> Poly a
`factorOutRoot` a
r
  where
    -- we can now assume that the polynomial has no roots at the boundary
    countRoots' :: Poly a -> Int
countRoots' Poly a
q = case Poly a -> Int
forall a. (Eq a, Num a) => Poly a -> Int
degree Poly a
q of
        -- q is the zero polynomial, so it doesn't *cross* zero
        -1 -> Int
0
        -- q is a non-zero constant polynomial - no root
        Int
0 -> Int
0
        -- q is a linear polynomial,
        Int
1 -> if Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval Poly a
q a
l a -> a -> a
forall a. Num a => a -> a -> a
* Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval Poly a
q a
r a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 then Int
1 else Int
0
        -- q has degree 2 or more so we can construct the Sturm sequence
        Int
_ -> (a, a, Poly a) -> Int
forall a. (Fractional a, Eq a, Ord a) => (a, a, Poly a) -> Int
countRootsSturm (a
l, a
r, Poly a
q)

-- | Given a polynomial \( p(x) \) and a value \( a \),
-- this functions factors out the polynomial \( (x-a)^m \),
-- where \( m \) is the highest power where this polynomial
-- divides \( p(x) \) without remainder.
--
-- * If the value \( a \) is a root of the polynomial,
--   then \( m \) is the multiplicity of the root.
-- * If the value \( a \) is not a root, then
--   \( m = 0 \) and the function returns \( p (x) \).
--
-- In other words, this function returns a polynomial \( q (x) \)
-- such that
--
-- \( p(x) = q(x)·(x - a)^m \)
--
-- where \( q(a) ≠ 0 \).
-- If the polynomial \( p(x) \) is identically 'zero',
-- we return 'zero' as well.
factorOutRoot :: (Fractional a, Ord a) => Poly a -> a -> Poly a
factorOutRoot :: forall a. (Fractional a, Ord a) => Poly a -> a -> Poly a
factorOutRoot Poly a
p0 a
x0
    | Poly a
p0 Poly a -> Poly a -> Bool
forall a. Eq a => a -> a -> Bool
== Poly a
forall a. Num a => Poly a
zero = Poly a
forall a. Num a => Poly a
zero
    | Bool
otherwise = Poly a -> Poly a
go Poly a
p0
  where
    go :: Poly a -> Poly a
go Poly a
p
        | Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval Poly a
p a
x0 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = Poly a -> a -> Poly a
forall a. (Fractional a, Ord a) => Poly a -> a -> Poly a
factorOutRoot Poly a
pDividedByXMinusX0 a
x0
        | Bool
otherwise = Poly a
p
      where
        xMinusX0 :: Poly a
xMinusX0 = Int -> a -> Poly a
forall a. (Eq a, Num a) => Int -> a -> Poly a
monomial Int
1 a
1 Poly a -> Poly a -> Poly a
forall a. Num a => a -> a -> a
- a -> Poly a
forall a. a -> Poly a
constant a
x0
        (Poly a
pDividedByXMinusX0, Poly a
_) = Poly a
p Poly a -> Poly a -> (Poly a, Poly a)
forall a.
(Fractional a, Eq a, Ord a) =>
Poly a -> Poly a -> (Poly a, Poly a)
`euclidianDivision` Poly a
xMinusX0

{-|
@'countRootsSturm' (x1, x2, p)@ returns the number of /distinct/ real roots
of the polynomial @p@ on the half-open interval \( (x_1, x_2] \),
under the following assumptions:

* @'degree' p >= 2@
* neither \( x_1 \) nor \( x_2 \) are multiple roots of \( p(x) \).

This function is an implementation of [Sturm's theorem
](https://en.wikipedia.org/wiki/Sturm%27s_theorem).
-}
countRootsSturm :: (Fractional a, Eq a, Ord a) => (a, a, Poly a) -> Int
countRootsSturm :: forall a. (Fractional a, Eq a, Ord a) => (a, a, Poly a) -> Int
countRootsSturm (a
l, a
r, Poly a
p) =
    -- p has degree 2 or more so we can construct the Sturm sequence
    [a] -> Int
forall a. (Fractional a, Ord a) => [a] -> Int
signVariations [a]
psl Int -> Int -> Int
forall a. Num a => a -> a -> a
- [a] -> Int
forall a. (Fractional a, Ord a) => [a] -> Int
signVariations [a]
psr
  where
    ps :: [Poly a]
ps = Poly a -> [Poly a]
forall a. (Fractional a, Ord a) => Poly a -> [Poly a]
reversedSturmSequence Poly a
p
    psl :: [a]
psl = (Poly a -> a) -> [Poly a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ((Poly a -> a -> a) -> a -> Poly a -> a
forall a b c. (a -> b -> c) -> b -> a -> c
flip Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval a
l) [Poly a]
ps
    psr :: [a]
psr = (Poly a -> a) -> [Poly a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ((Poly a -> a -> a) -> a -> Poly a -> a
forall a b c. (a -> b -> c) -> b -> a -> c
flip Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval a
r) [Poly a]
ps

{-| Number of sign variations in a list of real numbers.

Given a list @c0, c1, c2, . . . ck@,
then a sign variation (or sign change) in the sequence
is a pair of indices @i < j@ such that @ci*cj < 0@,
and either @j = i + 1@ or @ck = 0@ for all @@ such that @i < k < j@.
-}
signVariations :: (Fractional a, Ord a) => [a] -> Int
signVariations :: forall a. (Fractional a, Ord a) => [a] -> Int
signVariations [a]
xs =
    [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter (a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0) [a]
pairsMultiplied)
  where
    -- we simply remove zero elements to implement the clause
    -- "ck = 0 for all k such that i < k < j"
    zeroesRemoved :: [a]
zeroesRemoved = (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0) [a]
xs
    pairsMultiplied :: [a]
pairsMultiplied = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
zeroesRemoved (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
1 [a]
zeroesRemoved)

{-|
Construct the [Sturm sequence
](https://en.wikipedia.org/wiki/Sturm%27s_theorem)
of a given polynomial @p@. The Sturm sequence is given by the polynomials

> p0 = p
> p1 = differentiate p
> p{i+1} = - rem(p{i-1}, pi)

where @rem@ denotes the remainder under 'euclidianDivision'.
We truncate the list when one of the @pi = 0@.

For ease of implementation, we

* construct the 'reverse' of the Sturm sequence.
  This does not affect the number of sign variations that the usage site
  will be interested in.

* assume that the @degree p >= 1@.
-}
reversedSturmSequence :: (Fractional a, Ord a) => Poly a -> [Poly a]
reversedSturmSequence :: forall a. (Fractional a, Ord a) => Poly a -> [Poly a]
reversedSturmSequence Poly a
p =
    [Poly a] -> [Poly a]
forall {a}. (Fractional a, Ord a) => [Poly a] -> [Poly a]
go [Poly a -> Poly a
forall a. Num a => Poly a -> Poly a
differentiate Poly a
p, Poly a
p]
  where
    -- Note that this is called with a list of length 2 and grows the list,
    -- so we don't need to match all cases.
    go :: [Poly a] -> [Poly a]
go ps :: [Poly a]
ps@(Poly a
pI : Poly a
pIminusOne : [Poly a]
_)
        | Poly a
remainder Poly a -> Poly a -> Bool
forall a. Eq a => a -> a -> Bool
== Poly a
forall a. Num a => Poly a
zero = [Poly a]
ps
        | Bool
otherwise = [Poly a] -> [Poly a]
go (Poly a -> Poly a
forall a. Num a => a -> a
negate Poly a
remainder Poly a -> [Poly a] -> [Poly a]
forall a. a -> [a] -> [a]
: [Poly a]
ps)
      where
        remainder :: Poly a
remainder = (Poly a, Poly a) -> Poly a
forall a b. (a, b) -> b
snd ((Poly a, Poly a) -> Poly a) -> (Poly a, Poly a) -> Poly a
forall a b. (a -> b) -> a -> b
$ Poly a -> Poly a -> (Poly a, Poly a)
forall a.
(Fractional a, Eq a, Ord a) =>
Poly a -> Poly a -> (Poly a, Poly a)
euclidianDivision Poly a
pIminusOne Poly a
pI
    go [Poly a]
_ = String -> [Poly a]
forall a. HasCallStack => String -> a
error String
"reversedSturmSequence: impossible"

{-| Check whether a polynomial is monotonically increasing on
a given interval.
-}
isMonotonicallyIncreasingOn
    :: (Fractional a, Eq a, Ord a) => Poly a -> (a, a) -> Bool
isMonotonicallyIncreasingOn :: forall a. (Fractional a, Eq a, Ord a) => Poly a -> (a, a) -> Bool
isMonotonicallyIncreasingOn Poly a
p (a
x1, a
x2) =
    Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval Poly a
p a
x1 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval Poly a
p a
x2
        Bool -> Bool -> Bool
&& (a, a, Poly a) -> Int
forall a. (Fractional a, Ord a) => (a, a, Poly a) -> Int
countRoots (a
x1, a
x2, Poly a -> Poly a
forall a. Num a => Poly a -> Poly a
differentiate Poly a
p) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0

{-|
Measure whether or not a polynomial is consistently above or below zero,
or equals zero.

Need to consider special cases where there is a root at a boundary point.
-}
compareToZero :: (Fractional a, Eq a, Ord a) => (a, a, Poly a) -> Maybe Ordering
compareToZero :: forall a.
(Fractional a, Eq a, Ord a) =>
(a, a, Poly a) -> Maybe Ordering
compareToZero (a
l, a
u, Poly a
p)
    | a
l a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
u = String -> Maybe Ordering
forall a. HasCallStack => String -> a
error String
"Invalid interval"
    | Poly a
p Poly a -> Poly a -> Bool
forall a. Eq a => a -> a -> Bool
== Poly a
forall a. Num a => Poly a
zero = Ordering -> Maybe Ordering
forall a. a -> Maybe a
Just Ordering
EQ
    | a
lower a -> a -> a
forall a. Num a => a -> a -> a
* a
upper a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 = Maybe Ordering
forall a. Maybe a
Nothing -- quick test to eliminate simple cases
    | (a, a, Poly a) -> Int
forall a. (Fractional a, Ord a) => (a, a, Poly a) -> Int
countRoots (a
l, a
u, Poly a
p) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 = Maybe Ordering
forall a. Maybe a
Nothing -- polynomial crosses zero
    -- since the polynomial has no roots, the comparison is detmined by the boundary values
    | a
lower a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = Ordering -> Maybe Ordering
forall a. a -> Maybe a
Just (a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
upper a
lower)
    | a
upper a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = Ordering -> Maybe Ordering
forall a. a -> Maybe a
Just (a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
lower a
upper)
    | a
lower a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 = Ordering -> Maybe Ordering
forall a. a -> Maybe a
Just Ordering
GT -- upper must also be > 0 due to the lack of roots
    | Bool
otherwise = Ordering -> Maybe Ordering
forall a. a -> Maybe a
Just Ordering
LT -- upper and lower both < 0 due to the lack of roots
  where
    lower :: a
lower = Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval Poly a
p a
l
    upper :: a
upper = Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval Poly a
p a
u

{-|
Find a root of a polynomial in a given interval.

Return 'Nothing' if the polynomial does not have a root in the given interval.

We find the root by first forming the square-free factorisation of the polynomial,
to eliminate repeated roots. One of the factors may have a root in the interval,
so we count roots for each factor until we find the one with a root in the interval.
Then we use the bisection method to find the root,
repeatedly halving the interval in which the root must lie
until its width is less than the specified precision.
Constant and linear polynomials, @degree p <= 1@, are treated as special cases.
-}
findRoot
    :: forall a. (Fractional a, Eq a, Num a, Ord a) => a -> (a, a) -> Poly a -> Maybe a
findRoot :: forall a.
(Fractional a, Eq a, Num a, Ord a) =>
a -> (a, a) -> Poly a -> Maybe a
findRoot a
precision (a
lower, a
upper) Poly a
poly =
    if [Poly a] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Poly a]
rootFactors
        then Maybe a
forall a. Maybe a
Nothing
        else a -> (a, a) -> Poly a -> Maybe a
getRoot a
precision (a
lower, a
upper) ([Poly a] -> Poly a
forall a. HasCallStack => [a] -> a
head [Poly a]
rootFactors)
  where
    rootFactors :: [Poly a]
rootFactors =
        (Poly a -> Bool) -> [Poly a] -> [Poly a]
forall a. (a -> Bool) -> [a] -> [a]
filter (\Poly a
x -> (a, a, Poly a) -> Int
forall a. (Fractional a, Ord a) => (a, a, Poly a) -> Int
countRoots (a
lower, a
upper, Poly a
x) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0) (Poly a -> [Poly a]
forall a. (Fractional a, Eq a, Num a, Ord a) => Poly a -> [Poly a]
squareFreeFactorisation Poly a
poly)
    -- getRoot :: forall a. (Fractional a, Eq a, Num a, Ord a) => a -> (a, a) -> Poly a -> Maybe a
    getRoot :: a -> (a, a) -> Poly a -> Maybe a
getRoot a
eps (a
l, a
u) Poly a
p
        -- if the polynomial is zero, the whole interval is a root, so return the basepoint
        | Int
degp Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = a -> Maybe a
forall a. a -> Maybe a
Just a
l
        -- if the poly is a non-zero constant, no root is present
        | Int
degp Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = Maybe a
forall a. Maybe a
Nothing
        -- if the polynomial has degree 1, can calculate the root exactly
        | Int
degp Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = a -> Maybe a
forall a. a -> Maybe a
Just (-([a] -> a
forall a. HasCallStack => [a] -> a
head [a]
ps a -> a -> a
forall a. Fractional a => a -> a -> a
/ [a] -> a
forall a. HasCallStack => [a] -> a
last [a]
ps)) -- p0 + p1x = 0 => x = -p0/p1
        | a
eps a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 = String -> Maybe a
forall a. HasCallStack => String -> a
error String
"Invalid precision value"
        | Bool
otherwise = a -> (a, a) -> (a, a) -> Poly a -> Maybe a
(Fractional a, Eq a, Num a, Ord a) =>
a -> (a, a) -> (a, a) -> Poly a -> Maybe a
bisect a
eps (a
l, a
u) (Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval Poly a
p a
l, Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval Poly a
p a
u) Poly a
p
      where
        ps :: [a]
ps = Poly a -> [a]
forall a. Poly a -> [a]
toCoefficients Poly a
p
        degp :: Int
degp = Poly a -> Int
forall a. (Eq a, Num a) => Poly a -> Int
degree Poly a
p
        {- We bisect the interval exploiting the Intermediate Value Theorem:
        if a polynomial has different signs at the ends of an interval, it must be zero somewhere in the interval.
        If there is no change of sign, use countRoots to find which side of the interval the root is on.
        -}
        bisect
            :: (Fractional a, Eq a, Num a, Ord a) => a -> (a, a) -> (a, a) -> Poly a -> Maybe a
        bisect :: (Fractional a, Eq a, Num a, Ord a) =>
a -> (a, a) -> (a, a) -> Poly a -> Maybe a
bisect a
e (a
x, a
y) (a
px, a
py) Poly a
p'
            -- if we already have a root, choose it
            | a
px a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = a -> Maybe a
forall a. a -> Maybe a
Just a
x
            | a
py a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = a -> Maybe a
forall a. a -> Maybe a
Just a
y
            | a
pmid a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = a -> Maybe a
forall a. a -> Maybe a
Just a
mid
            -- when the interval is small enough, stop:
            -- the root is in this interval, so take the mid point
            | a
width a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
e = a -> Maybe a
forall a. a -> Maybe a
Just a
mid
            -- choose the lower half, if the polynomial has different signs at the ends
            | a -> a
forall a. Num a => a -> a
signum a
px a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a -> a
forall a. Num a => a -> a
signum a
pmid = a -> (a, a) -> (a, a) -> Poly a -> Maybe a
(Fractional a, Eq a, Num a, Ord a) =>
a -> (a, a) -> (a, a) -> Poly a -> Maybe a
bisect a
e (a
x, a
mid) (a
px, a
pmid) Poly a
p'
            -- choose the upper half, if the polynomial has different signs at the ends
            | a -> a
forall a. Num a => a -> a
signum a
py a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a -> a
forall a. Num a => a -> a
signum a
pmid = a -> (a, a) -> (a, a) -> Poly a -> Maybe a
(Fractional a, Eq a, Num a, Ord a) =>
a -> (a, a) -> (a, a) -> Poly a -> Maybe a
bisect a
e (a
mid, a
y) (a
pmid, a
py) Poly a
p'
            -- no sign change found, so we resort to counting roots
            | (a, a, Poly a) -> Int
forall a. (Fractional a, Ord a) => (a, a, Poly a) -> Int
countRoots (a
x, a
mid, Poly a
p') Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 = a -> (a, a) -> (a, a) -> Poly a -> Maybe a
(Fractional a, Eq a, Num a, Ord a) =>
a -> (a, a) -> (a, a) -> Poly a -> Maybe a
bisect a
e (a
x, a
mid) (a
px, a
pmid) Poly a
p'
            | (a, a, Poly a) -> Int
forall a. (Fractional a, Ord a) => (a, a, Poly a) -> Int
countRoots (a
mid, a
y, Poly a
p') Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 = a -> (a, a) -> (a, a) -> Poly a -> Maybe a
(Fractional a, Eq a, Num a, Ord a) =>
a -> (a, a) -> (a, a) -> Poly a -> Maybe a
bisect a
e (a
mid, a
y) (a
pmid, a
py) Poly a
p'
            | Bool
otherwise = Maybe a
forall a. Maybe a
Nothing
          where
            width :: a
width = a
y a -> a -> a
forall a. Num a => a -> a -> a
- a
x
            mid :: a
mid = a
x a -> a -> a
forall a. Num a => a -> a -> a
+ a
width a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
2
            pmid :: a
pmid = Poly a -> a -> a
forall a. Num a => Poly a -> a -> a
eval Poly a
p' a
mid

{-| We are seeking the point at which a polynomial has a specific value.:
subtract the value we are looking for so that we seek a zero crossing
-}
root -- TODO: this should probably be called something else such as 'findCrossing'
    :: (Ord a, Num a, Eq a, Fractional a)
    => a
    -> a
    -> (a, a)
    -> Poly a
    -> Maybe a
root :: forall a.
(Ord a, Num a, Eq a, Fractional a) =>
a -> a -> (a, a) -> Poly a -> Maybe a
root a
e a
x (a
l, a
u) Poly a
p = a -> (a, a) -> Poly a -> Maybe a
forall a.
(Fractional a, Eq a, Num a, Ord a) =>
a -> (a, a) -> Poly a -> Maybe a
findRoot a
e (a
l, a
u) (Poly a
p Poly a -> Poly a -> Poly a
forall a. Num a => a -> a -> a
- a -> Poly a
forall a. a -> Poly a
constant a
x)

-- | Greatest monic common divisor of two polynomials.
gcdPoly
    :: forall a. (Fractional a, Eq a, Num a, Ord a) => Poly a -> Poly a -> Poly a
gcdPoly :: forall a.
(Fractional a, Eq a, Num a, Ord a) =>
Poly a -> Poly a -> Poly a
gcdPoly Poly a
a Poly a
b = if Poly a
b Poly a -> Poly a -> Bool
forall a. Eq a => a -> a -> Bool
== Poly a
forall a. Num a => Poly a
zero then Poly a
a else Poly a -> Poly a
makeMonic (Poly a -> Poly a -> Poly a
forall a.
(Fractional a, Eq a, Num a, Ord a) =>
Poly a -> Poly a -> Poly a
gcdPoly Poly a
b (Poly a -> Poly a -> Poly a
polyRemainder Poly a
a Poly a
b))
  where
    makeMonic :: Poly a -> Poly a
    makeMonic :: Poly a -> Poly a
makeMonic (Poly [a]
as) = a -> Poly a -> Poly a
forall a. Num a => a -> Poly a -> Poly a
scale (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ [a] -> a
forall a. HasCallStack => [a] -> a
last [a]
as) ([a] -> Poly a
forall a. [a] -> Poly a
Poly [a]
as)
    polyRemainder :: Poly a -> Poly a -> Poly a
    polyRemainder :: Poly a -> Poly a -> Poly a
polyRemainder Poly a
x Poly a
y = (Poly a, Poly a) -> Poly a
forall a b. (a, b) -> b
snd (Poly a -> Poly a -> (Poly a, Poly a)
forall a.
(Fractional a, Eq a, Ord a) =>
Poly a -> Poly a -> (Poly a, Poly a)
euclidianDivision Poly a
x Poly a
y)

{-|
We compute the square-free factorisation of a polynomial using Yun's algorithm.
Yun, David Y.Y. (1976). "On square-free decomposition algorithms".
SYMSAC '76 Proceedings of the third ACM Symposium on Symbolic and Algebraic Computation.
Association for Computing Machinery. pp. 26–35. doi:10.1145/800205.806320. ISBN 978-1-4503-7790-4. S2CID 12861227.
https://dl.acm.org/doi/10.1145/800205.806320
G <- gcd (P, P')
C1 <- P / G
D1 <- P' / G - C1'
until Ci = 1 do
    Pi <- gcd (Ci, Di)
    Ci+1 <- Ci/Pi
    Di+1 <- Di / Ai - Ci+1'
-}
squareFreeFactorisation
    :: (Fractional a, Eq a, Num a, Ord a) => Poly a -> [Poly a]
squareFreeFactorisation :: forall a. (Fractional a, Eq a, Num a, Ord a) => Poly a -> [Poly a]
squareFreeFactorisation Poly a
p = 
    -- if p has degree <= 1 it can have no factors but itself
    if Poly a -> Int
forall a. (Eq a, Num a) => Poly a -> Int
degree Poly a
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1 then [Poly a
p] else Poly a -> Poly a -> [Poly a]
forall {a}. (Fractional a, Ord a) => Poly a -> Poly a -> [Poly a]
go Poly a
c1 Poly a
d1  
  where
    diffP :: Poly a
diffP = Poly a -> Poly a
forall a. Num a => Poly a -> Poly a
differentiate Poly a
p
    g0 :: Poly a
g0 = Poly a -> Poly a -> Poly a
forall a.
(Fractional a, Eq a, Num a, Ord a) =>
Poly a -> Poly a -> Poly a
gcdPoly Poly a
p Poly a
diffP
    c1 :: Poly a
c1 = Poly a
p Poly a -> Poly a -> Poly a
forall {a}. (Fractional a, Ord a) => Poly a -> Poly a -> Poly a
`divide` Poly a
g0
    d1 :: Poly a
d1 = (Poly a
diffP Poly a -> Poly a -> Poly a
forall {a}. (Fractional a, Ord a) => Poly a -> Poly a -> Poly a
`divide` Poly a
g0) Poly a -> Poly a -> Poly a
forall a. Num a => a -> a -> a
- Poly a -> Poly a
forall a. Num a => Poly a -> Poly a
differentiate Poly a
c1
    divide :: Poly a -> Poly a -> Poly a
divide Poly a
x Poly a
y = (Poly a, Poly a) -> Poly a
forall a b. (a, b) -> a
fst (Poly a -> Poly a -> (Poly a, Poly a)
forall a.
(Fractional a, Eq a, Ord a) =>
Poly a -> Poly a -> (Poly a, Poly a)
euclidianDivision Poly a
x Poly a
y)
    go :: Poly a -> Poly a -> [Poly a]
go Poly a
c Poly a
d
        | Poly a
c Poly a -> Poly a -> Bool
forall a. Eq a => a -> a -> Bool
== a -> Poly a
forall a. a -> Poly a
constant a
1 = [] -- terminate the recursion
        | Poly a
a' Poly a -> Poly a -> Bool
forall a. Eq a => a -> a -> Bool
== a -> Poly a
forall a. a -> Poly a
constant a
1 = Poly a -> Poly a -> [Poly a]
go Poly a
c' Poly a
d' -- skip over the constant polynomial
        | Bool
otherwise = Poly a
a' Poly a -> [Poly a] -> [Poly a]
forall a. a -> [a] -> [a]
: Poly a -> Poly a -> [Poly a]
go Poly a
c' Poly a
d'
      where
        a' :: Poly a
a' = Poly a -> Poly a -> Poly a
forall a.
(Fractional a, Eq a, Num a, Ord a) =>
Poly a -> Poly a -> Poly a
gcdPoly Poly a
c Poly a
d
        c' :: Poly a
c' = Poly a
c Poly a -> Poly a -> Poly a
forall {a}. (Fractional a, Ord a) => Poly a -> Poly a -> Poly a
`divide` Poly a
a'
        d' :: Poly a
d' = (Poly a
d Poly a -> Poly a -> Poly a
forall {a}. (Fractional a, Ord a) => Poly a -> Poly a -> Poly a
`divide` Poly a
a') Poly a -> Poly a -> Poly a
forall a. Num a => a -> a -> a
- Poly a -> Poly a
forall a. Num a => Poly a -> Poly a
differentiate Poly a
c'