-- SPDX-License-Identifier: MPL-2.0

{- | Univariate polynomials with ascending coefficients: @[a0,a1,...]@ means
@a0+a1*x+...@. Trailing exact zeros are removed; zero is @[]@ with degree
'Nothing'. Exact algorithms use Rational. Floating root solvers use Complex
Double with explicit tolerance and iteration budgets, returning failure when
residuals or corrections do not converge.
-}
module HMath.Polynomial (
    Polynomial,
    polynomial,
    coefficients,
    polynomialDegree,
    evaluatePolynomial,
    addPolynomial,
    subtractPolynomial,
    scalePolynomial,
    multiplyPolynomial,
    composePolynomial,
    dividePolynomial,
    monicPolynomial,
    gcdPolynomial,
    differentiatePolynomial,
    integratePolynomial,
    squareFreeFactors,
    rationalRoots,
    factorRationalPolynomial,
    sturmSequence,
    countRealRoots,
    isolateRealRoots,
    durandKerner,
    aberthEhrlich,
    companionMatrix,
    lagrangeInterpolation,
    newtonInterpolation,
    chebyshevPolynomial,
    legendrePolynomial,
    hermitePolynomial,
    laguerrePolynomial,
    bernsteinPolynomial,
    polynomialRegression,
    chebyshevApproximation,
    PowerSeries,
    powerSeries,
    seriesCoefficients,
    addSeries,
    multiplySeries,
    reciprocalSeries,
    divideSeries,
    composeSeries,
    differentiateSeries,
    integrateSeries,
    expSeries,
    logSeries,
    RationalFunction,
    rationalFunction,
    rationalNumerator,
    rationalDenominator,
    evaluateRationalFunction,
    addRationalFunction,
    multiplyRationalFunction,
    partialFractions,
) where

import Control.Monad (foldM)
import Data.Complex (Complex ((:+)), magnitude, mkPolar)
import Data.List (nub, sort)
import Data.Ratio (denominator, numerator, (%))
import HMath.Error (MathError (..))
import qualified HMath.Matrix as M
import HMath.Matrix.QR (leastSquares)
import HMath.NumberTheory (divisors)

-- | Canonical, finite ascending coefficient sequence.
newtype Polynomial a = Polynomial [a] deriving (Polynomial a -> Polynomial a -> Bool
(Polynomial a -> Polynomial a -> Bool)
-> (Polynomial a -> Polynomial a -> Bool) -> Eq (Polynomial a)
forall a. Eq a => Polynomial a -> Polynomial a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => Polynomial a -> Polynomial a -> Bool
== :: Polynomial a -> Polynomial a -> Bool
$c/= :: forall a. Eq a => Polynomial a -> Polynomial a -> Bool
/= :: Polynomial a -> Polynomial a -> Bool
Eq, Int -> Polynomial a -> ShowS
[Polynomial a] -> ShowS
Polynomial a -> String
(Int -> Polynomial a -> ShowS)
-> (Polynomial a -> String)
-> ([Polynomial a] -> ShowS)
-> Show (Polynomial a)
forall a. Show a => Int -> Polynomial a -> ShowS
forall a. Show a => [Polynomial a] -> ShowS
forall a. Show a => Polynomial a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Polynomial a -> ShowS
showsPrec :: Int -> Polynomial a -> ShowS
$cshow :: forall a. Show a => Polynomial a -> String
show :: Polynomial a -> String
$cshowList :: forall a. Show a => [Polynomial a] -> ShowS
showList :: [Polynomial a] -> ShowS
Show)

-- | Construct a polynomial, deleting trailing exact zeros (no floating tolerance).
polynomial :: (Eq a, Num a) => [a] -> Polynomial a
polynomial :: forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial = [a] -> Polynomial a
forall a. [a] -> Polynomial a
Polynomial ([a] -> Polynomial a) -> ([a] -> [a]) -> [a] -> Polynomial a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [a]
forall a. [a] -> [a]
reverse ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0) ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [a]
forall a. [a] -> [a]
reverse

-- | Ascending coefficients, with no trailing zero; zero gives an empty list.
coefficients :: Polynomial a -> [a]
coefficients :: forall a. Polynomial a -> [a]
coefficients (Polynomial [a]
xs) = [a]
xs

-- | Degree of a nonzero polynomial; zero has no degree.
polynomialDegree :: Polynomial a -> Maybe Int
polynomialDegree :: forall a. Polynomial a -> Maybe Int
polynomialDegree (Polynomial []) = Maybe Int
forall a. Maybe a
Nothing
polynomialDegree (Polynomial [a]
xs) = Int -> Maybe Int
forall a. a -> Maybe a
Just ([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)

-- | Horner evaluation, with the constant coefficient first in storage.
evaluatePolynomial :: (Num a) => Polynomial a -> a -> a
evaluatePolynomial :: forall a. Num a => Polynomial a -> a -> a
evaluatePolynomial (Polynomial [a]
xs) 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
a a
b -> a
a a -> a -> a
forall a. Num a => a -> a -> a
+ a
x a -> a -> a
forall a. Num a => a -> a -> a
* a
b) a
0 [a]
xs

zipPad :: (Num a) => (a -> a -> a) -> [a] -> [a] -> [a]
zipPad :: forall a. Num a => (a -> a -> a) -> [a] -> [a] -> [a]
zipPad a -> a -> a
f [] [a]
ys = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
f a
0) [a]
ys
zipPad a -> a -> a
f [a]
xs [] = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
`f` a
0) [a]
xs
zipPad a -> a -> a
f (a
x : [a]
xs) (a
y : [a]
ys) = a -> a -> a
f a
x a
y a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a -> a) -> [a] -> [a] -> [a]
forall a. Num a => (a -> a -> a) -> [a] -> [a] -> [a]
zipPad a -> a -> a
f [a]
xs [a]
ys

-- | Coefficientwise addition with canonical exact-zero normalization.
addPolynomial :: (Eq a, Num a) => Polynomial a -> Polynomial a -> Polynomial a
addPolynomial :: forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
addPolynomial (Polynomial [a]
a) (Polynomial [a]
b) = [a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial ((a -> a -> a) -> [a] -> [a] -> [a]
forall a. Num a => (a -> a -> a) -> [a] -> [a] -> [a]
zipPad a -> a -> a
forall a. Num a => a -> a -> a
(+) [a]
a [a]
b)

-- | Coefficientwise subtraction.
subtractPolynomial :: (Eq a, Num a) => Polynomial a -> Polynomial a -> Polynomial a
subtractPolynomial :: forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
subtractPolynomial (Polynomial [a]
a) (Polynomial [a]
b) = [a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial ((a -> a -> a) -> [a] -> [a] -> [a]
forall a. Num a => (a -> a -> a) -> [a] -> [a] -> [a]
zipPad (-) [a]
a [a]
b)

-- | Multiply every coefficient by a scalar.
scalePolynomial :: (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial :: forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial a
a (Polynomial [a]
xs) = [a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
a a -> a -> a
forall a. Num a => a -> a -> a
*) [a]
xs)

-- | Exact schoolbook convolution, including zero polynomials.
multiplyPolynomial :: (Eq a, Num a) => Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial :: forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial (Polynomial [a]
a) (Polynomial [a]
b) = [a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial ([a] -> [a]
go [a]
a)
  where
    go :: [a] -> [a]
go [] = []
    go (a
x : [a]
xs) = (a -> a -> a) -> [a] -> [a] -> [a]
forall a. Num a => (a -> a -> a) -> [a] -> [a] -> [a]
zipPad a -> a -> a
forall a. Num a => a -> a -> a
(+) ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
x a -> a -> a
forall a. Num a => a -> a -> a
*) [a]
b) (a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
go [a]
xs)

-- | Composition @p(q(x))@ by polynomial Horner evaluation.
composePolynomial :: (Eq a, Num a) => Polynomial a -> Polynomial a -> Polynomial a
composePolynomial :: forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
composePolynomial (Polynomial [a]
xs) Polynomial a
q =
    (a -> Polynomial a -> Polynomial a)
-> Polynomial a -> [a] -> Polynomial 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
a Polynomial a
r ->
            Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
addPolynomial
                ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
a])
                (Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial Polynomial a
q Polynomial a
r)
        )
        ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [])
        [a]
xs

{- | Long division over a field. Exact coefficient equality controls termination;
numerical cancellation is truncated at the eliminated leading coefficient.
Zero denominator is rejected.
-}
dividePolynomial ::
    (Eq a, Fractional a) =>
    Polynomial a ->
    Polynomial a ->
    Either MathError (Polynomial a, Polynomial a)
dividePolynomial :: forall a.
(Eq a, Fractional a) =>
Polynomial a
-> Polynomial a -> Either MathError (Polynomial a, Polynomial a)
dividePolynomial Polynomial a
_ (Polynomial []) = MathError -> Either MathError (Polynomial a, Polynomial a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"zero polynomial divisor")
dividePolynomial (Polynomial [a]
xs) (Polynomial [a]
ys) = (Polynomial a, Polynomial a)
-> Either MathError (Polynomial a, Polynomial a)
forall a b. b -> Either a b
Right ([a] -> Polynomial a -> (Polynomial a, Polynomial a)
go [a]
xs ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial []))
  where
    n :: Int
n = [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
ys
    lead :: a
lead = (a -> a -> a) -> a -> [a] -> a
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\a
_ a
x -> a
x) a
0 [a]
ys
    go :: [a] -> Polynomial a -> (Polynomial a, Polynomial a)
go [a]
r Polynomial a
q
        | [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
r Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n = (Polynomial a
q, [a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a]
r)
        | Bool
otherwise =
            let k :: Int
k = [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
r Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
n
                a :: a
a = (a -> a -> a) -> a -> [a] -> a
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\a
_ a
x -> a
x) a
0 [a]
r a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
lead
                term :: Polynomial a
term = [a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial (Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
k a
0 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
a])
                eliminated :: [a]
eliminated = (a -> a -> a) -> [a] -> [a] -> [a]
forall a. Num a => (a -> a -> a) -> [a] -> [a] -> [a]
zipPad (-) [a]
r (Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
k a
0 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
a a -> a -> a
forall a. Num a => a -> a -> a
*) [a]
ys)
                lower :: Polynomial a
lower = [a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
r Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) [a]
eliminated)
             in [a] -> Polynomial a -> (Polynomial a, Polynomial a)
go (Polynomial a -> [a]
forall a. Polynomial a -> [a]
coefficients Polynomial a
lower) (Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
addPolynomial Polynomial a
q Polynomial a
term)

-- | Normalize a nonzero polynomial to leading coefficient one; zero remains zero.
monicPolynomial :: (Eq a, Fractional a) => Polynomial a -> Polynomial a
monicPolynomial :: forall a. (Eq a, Fractional a) => Polynomial a -> Polynomial a
monicPolynomial p :: Polynomial a
p@(Polynomial []) = Polynomial a
p
monicPolynomial p :: Polynomial a
p@(Polynomial [a]
xs) = a -> Polynomial a -> Polynomial a
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial (a -> a
forall a. Fractional a => a -> a
recip ((a -> a -> a) -> a -> [a] -> a
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\a
_ a
a -> a
a) a
0 [a]
xs)) Polynomial a
p

-- | Monic Euclidean GCD over an exact field; @gcd(0,0)=0@.
gcdPolynomial :: (Eq a, Fractional a) => Polynomial a -> Polynomial a -> Polynomial a
gcdPolynomial :: forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
gcdPolynomial Polynomial a
a (Polynomial []) = Polynomial a -> Polynomial a
forall a. (Eq a, Fractional a) => Polynomial a -> Polynomial a
monicPolynomial Polynomial a
a
gcdPolynomial Polynomial a
a Polynomial a
b = case Polynomial a
-> Polynomial a -> Either MathError (Polynomial a, Polynomial a)
forall a.
(Eq a, Fractional a) =>
Polynomial a
-> Polynomial a -> Either MathError (Polynomial a, Polynomial a)
dividePolynomial Polynomial a
a Polynomial a
b of
    Right (Polynomial a
_, Polynomial a
r) -> Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
gcdPolynomial Polynomial a
b Polynomial a
r
    Left MathError
_ -> Polynomial a -> Polynomial a
forall a. (Eq a, Fractional a) => Polynomial a -> Polynomial a
monicPolynomial Polynomial a
a

-- | Formal derivative.
differentiatePolynomial :: (Eq a, Num a) => Polynomial a -> Polynomial a
differentiatePolynomial :: forall a. (Eq a, Num a) => Polynomial a -> Polynomial a
differentiatePolynomial (Polynomial [a]
xs) = [a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial ((Integer -> a -> a) -> [Integer] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Integer
i a
a -> Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
i a -> a -> a
forall a. Num a => a -> a -> a
* a
a) [Integer
1 ..] (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
1 [a]
xs))

-- | Formal integral with the supplied integration constant.
integratePolynomial :: (Eq a, Fractional a) => a -> Polynomial a -> Polynomial a
integratePolynomial :: forall a. (Eq a, Fractional a) => a -> Polynomial a -> Polynomial a
integratePolynomial a
c (Polynomial [a]
xs) = [a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial (a
c a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (Integer -> a -> a) -> [Integer] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Integer
i a
a -> a
a a -> a -> a
forall a. Fractional a => a -> a -> a
/ Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
i) [Integer
1 ..] [a]
xs)

{- | Monic square-free factors with positive multiplicities over Rational.
The leading coefficient is omitted; constant polynomials have no factors.
-}
squareFreeFactors :: Polynomial Rational -> [(Polynomial Rational, Int)]
squareFreeFactors :: Polynomial Rational -> [(Polynomial Rational, Int)]
squareFreeFactors Polynomial Rational
p
    | Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 = []
    | Bool
otherwise = Int
-> Polynomial Rational
-> Polynomial Rational
-> [(Polynomial Rational, Int)]
forall {a} {t}.
(Eq a, Fractional a, Num t) =>
t -> Polynomial a -> Polynomial a -> [(Polynomial a, t)]
go Int
1 Polynomial Rational
w Polynomial Rational
c
  where
    f :: Polynomial Rational
f = Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Fractional a) => Polynomial a -> Polynomial a
monicPolynomial Polynomial Rational
p
    c :: Polynomial Rational
c = Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
gcdPolynomial Polynomial Rational
f (Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Num a) => Polynomial a -> Polynomial a
differentiatePolynomial Polynomial Rational
f)
    w :: Polynomial Rational
w = Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
quotient Polynomial Rational
f Polynomial Rational
c
    go :: t -> Polynomial a -> Polynomial a -> [(Polynomial a, t)]
go t
i Polynomial a
a Polynomial a
b
        | Polynomial a -> Int
forall a. Polynomial a -> Int
deg Polynomial a
a Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 = []
        | Bool
otherwise =
            let y :: Polynomial a
y = Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
gcdPolynomial Polynomial a
a Polynomial a
b; z :: Polynomial a
z = Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
quotient Polynomial a
a Polynomial a
y
             in (if Polynomial a -> Int
forall a. Polynomial a -> Int
deg Polynomial a
z Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 then [(Polynomial a
z, t
i)] else []) [(Polynomial a, t)] -> [(Polynomial a, t)] -> [(Polynomial a, t)]
forall a. [a] -> [a] -> [a]
++ t -> Polynomial a -> Polynomial a -> [(Polynomial a, t)]
go (t
i t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) Polynomial a
y (Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
quotient Polynomial a
b Polynomial a
y)

deg :: Polynomial a -> Int
deg :: forall a. Polynomial a -> Int
deg Polynomial a
p = Int -> (Int -> Int) -> Maybe Int -> Int
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (-Int
1) Int -> Int
forall a. a -> a
id (Polynomial a -> Maybe Int
forall a. Polynomial a -> Maybe Int
polynomialDegree Polynomial a
p)

quotient :: (Eq a, Fractional a) => Polynomial a -> Polynomial a -> Polynomial a
quotient :: forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
quotient Polynomial a
a Polynomial a
b = (MathError -> Polynomial a)
-> ((Polynomial a, Polynomial a) -> Polynomial a)
-> Either MathError (Polynomial a, Polynomial a)
-> Polynomial a
forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either (Polynomial a -> MathError -> Polynomial a
forall a b. a -> b -> a
const ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [])) (Polynomial a, Polynomial a) -> Polynomial a
forall a b. (a, b) -> a
fst (Polynomial a
-> Polynomial a -> Either MathError (Polynomial a, Polynomial a)
forall a.
(Eq a, Fractional a) =>
Polynomial a
-> Polynomial a -> Either MathError (Polynomial a, Polynomial a)
dividePolynomial Polynomial a
a Polynomial a
b)

{- | Rational roots with exact multiplicities using the rational-root theorem.
Budget applies to integer factorization of numerator and leading coefficient.
-}
rationalRoots :: Int -> Polynomial Rational -> Either MathError [(Rational, Int)]
rationalRoots :: Int -> Polynomial Rational -> Either MathError [(Rational, Int)]
rationalRoots Int
budget Polynomial Rational
p
    | Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError [(Rational, Int)]
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"negative rational-root budget")
    | Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError [(Rational, Int)]
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"zero polynomial has infinitely many roots")
    | Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = [(Rational, Int)] -> Either MathError [(Rational, Int)]
forall a b. b -> Either a b
Right []
    | Bool
otherwise = do
        let xs :: [Integer]
xs = Polynomial Rational -> [Integer]
integerCoefficients Polynomial Rational
p
            zeroCount :: Int
zeroCount = [Integer] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0) [Integer]
xs)
            rest :: [Integer]
rest = Int -> [Integer] -> [Integer]
forall a. Int -> [a] -> [a]
drop Int
zeroCount [Integer]
xs
            c :: Integer
c = case [Integer]
rest of Integer
a : [Integer]
_ -> Integer -> Integer
forall a. Num a => a -> a
abs Integer
a; [] -> Integer
1
            l :: Integer
l = Integer -> Integer
forall a. Num a => a -> a
abs ((Integer -> Integer -> Integer) -> Integer -> [Integer] -> Integer
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\Integer
_ Integer
a -> Integer
a) Integer
1 [Integer]
rest)
        ns <- Int -> Integer -> Either MathError [Integer]
divisors Int
budget Integer
c
        ds <- divisors budget l
        let candidates = [Rational] -> [Rational]
forall a. Eq a => [a] -> [a]
nub [Integer
sign Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
d | Integer
sign <- [-Integer
1, Integer
1], Integer
n <- [Integer]
ns, Integer
d <- [Integer]
ds]
            multiplicity a
r Polynomial a
q t
k
                | Polynomial a -> Int
forall a. Polynomial a -> Int
deg Polynomial a
q Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Polynomial a -> a -> a
forall a. Num a => Polynomial a -> a -> a
evaluatePolynomial Polynomial a
q a
r a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 =
                    a -> Polynomial a -> t -> t
multiplicity a
r (Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
quotient Polynomial a
q ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [-a
r, a
1])) (t
k t -> t -> t
forall a. Num a => a -> a -> a
+ t
1)
                | Bool
otherwise = t
k
        pure
            ( (if zeroCount > 0 then [(0, zeroCount)] else [])
                ++ [(r, k) | r <- sort candidates, let k = Rational -> Polynomial Rational -> Int -> Int
forall {a} {t}.
(Eq a, Fractional a, Num t) =>
a -> Polynomial a -> t -> t
multiplicity Rational
r Polynomial Rational
p Int
0, k > 0]
            )

integerCoefficients :: Polynomial Rational -> [Integer]
integerCoefficients :: Polynomial Rational -> [Integer]
integerCoefficients Polynomial Rational
p =
    let xs :: [Rational]
xs = Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
coefficients Polynomial Rational
p
        den :: Integer
den = (Integer -> Integer -> Integer) -> Integer -> [Integer] -> Integer
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
lcm Integer
1 ((Rational -> Integer) -> [Rational] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map Rational -> Integer
forall a. Ratio a -> a
denominator [Rational]
xs)
        ints :: [Integer]
ints = (Rational -> Integer) -> [Rational] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (\Rational
a -> Rational -> Integer
forall a. Ratio a -> a
numerator (Rational
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Integer -> Rational
forall a. Num a => Integer -> a
fromInteger Integer
den)) [Rational]
xs
        g :: Integer
g = (Integer -> Integer -> Integer) -> Integer -> [Integer] -> Integer
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
gcd Integer
0 [Integer]
ints
     in if Integer
g Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 then [] else (Integer -> Integer) -> [Integer] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
g) [Integer]
ints

{- | Complete rational factorization by bounded Kronecker interpolation. Returns
the original leading scalar and monic irreducible factors with multiplicity.
Every candidate is divided exactly. Exhausting the candidate budget reports
failure, never an unsupported claim of irreducibility. Intended for small
exact polynomials; coefficient divisor enumeration also uses this budget.
-}
factorRationalPolynomial ::
    Int ->
    Polynomial Rational ->
    Either MathError (Rational, [(Polynomial Rational, Int)])
factorRationalPolynomial :: Int
-> Polynomial Rational
-> Either MathError (Rational, [(Polynomial Rational, Int)])
factorRationalPolynomial Int
budget Polynomial Rational
p
    | Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError
-> Either MathError (Rational, [(Polynomial Rational, Int)])
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"negative factorization budget")
    | Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError
-> Either MathError (Rational, [(Polynomial Rational, Int)])
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"zero polynomial cannot be factored")
    | Bool
otherwise = do
        let leading :: Rational
leading = (Rational -> Rational -> Rational)
-> Rational -> [Rational] -> Rational
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\Rational
_ Rational
x -> Rational
x) Rational
0 (Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
coefficients Polynomial Rational
p)
        (fs, _) <- Int
-> Polynomial Rational
-> Either MathError ([Polynomial Rational], Int)
forall {t}.
(Eq t, Num t) =>
t
-> Polynomial Rational
-> Either MathError ([Polynomial Rational], t)
factor Int
budget (Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Fractional a) => Polynomial a -> Polynomial a
monicPolynomial Polynomial Rational
p)
        pure (leading, combineFactors fs)
  where
    combineFactors :: [a] -> [(a, Int)]
combineFactors [] = []
    combineFactors (a
f : [a]
fs) = (a
f, Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ [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. Eq a => a -> a -> Bool
== a
f) [a]
fs)) (a, Int) -> [(a, Int)] -> [(a, Int)]
forall a. a -> [a] -> [a]
: [a] -> [(a, Int)]
combineFactors ((a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
f) [a]
fs)
    factor :: t
-> Polynomial Rational
-> Either MathError ([Polynomial Rational], t)
factor t
left Polynomial Rational
f
        | Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
f Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 = ([Polynomial Rational], t)
-> Either MathError ([Polynomial Rational], t)
forall a b. b -> Either a b
Right ([], t
left)
        | Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
f Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = ([Polynomial Rational], t)
-> Either MathError ([Polynomial Rational], t)
forall a b. b -> Either a b
Right ([Polynomial Rational
f], t
left)
        | Bool
otherwise = do
            found <- t
-> Polynomial Rational
-> Int
-> Either MathError (Maybe (Polynomial Rational), t)
forall {b}.
(Eq b, Num b) =>
b
-> Polynomial Rational
-> Int
-> Either MathError (Maybe (Polynomial Rational), b)
findFactor t
left Polynomial Rational
f Int
1
            case found of
                (Maybe (Polynomial Rational)
Nothing, t
remain) -> ([Polynomial Rational], t)
-> Either MathError ([Polynomial Rational], t)
forall a b. b -> Either a b
Right ([Polynomial Rational
f], t
remain)
                (Just Polynomial Rational
g, t
remain) -> do
                    (a, l) <- t
-> Polynomial Rational
-> Either MathError ([Polynomial Rational], t)
factor t
remain (Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Fractional a) => Polynomial a -> Polynomial a
monicPolynomial Polynomial Rational
g)
                    (b, r) <- factor l (monicPolynomial (quotient f g))
                    pure (a ++ b, r)
    findFactor :: b
-> Polynomial Rational
-> Int
-> Either MathError (Maybe (Polynomial Rational), b)
findFactor b
left Polynomial Rational
f Int
k
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
f Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2 = (Maybe (Polynomial Rational), b)
-> Either MathError (Maybe (Polynomial Rational), b)
forall a b. b -> Either a b
Right (Maybe (Polynomial Rational)
forall a. Maybe a
Nothing, b
left)
        | Bool
otherwise = do
            let ints :: Polynomial Rational
ints = [Rational] -> Polynomial Rational
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial ((Integer -> Rational) -> [Integer] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Polynomial Rational -> [Integer]
integerCoefficients Polynomial Rational
f))
                samplePoints :: [Rational]
samplePoints = Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
take (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Rational
x | Rational
x <- Rational
0 Rational -> [Rational] -> [Rational]
forall a. a -> [a] -> [a]
: [[Rational]] -> [Rational]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Rational
j, -Rational
j] | Rational
j <- [Rational
1 ..]], Polynomial Rational -> Rational -> Rational
forall a. Num a => Polynomial a -> a -> a
evaluatePolynomial Polynomial Rational
ints Rational
x Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0]
            values <- (Rational -> Either MathError [Integer])
-> [Rational] -> Either MathError [[Integer]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM (\Rational
x -> Int -> Integer -> Either MathError [Integer]
divisors Int
budget (Integer -> Integer
forall a. Num a => a -> a
abs (Rational -> Integer
forall a. Ratio a -> a
numerator (Polynomial Rational -> Rational -> Rational
forall a. Num a => Polynomial a -> a -> a
evaluatePolynomial Polynomial Rational
ints Rational
x)))) [Rational]
samplePoints
            let choices = ([Integer] -> [Rational]) -> [[Integer]] -> [[Rational]]
forall a b. (a -> b) -> [a] -> [b]
map (\[Integer]
ds -> (Integer -> [Rational]) -> [Integer] -> [Rational]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (\Integer
d -> [Integer -> Rational
forall a. Num a => Integer -> a
fromInteger Integer
d, Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (-Integer
d)]) [Integer]
ds) [[Integer]]
values
            attempt left (sequence choices)
      where
        points :: [Rational]
points =
            Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
take
                (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                [ Rational
x
                | Rational
x <- Rational
0 Rational -> [Rational] -> [Rational]
forall a. a -> [a] -> [a]
: [[Rational]] -> [Rational]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Rational
j, -Rational
j] | Rational
j <- [Rational
1 ..]]
                , Polynomial Rational -> Rational -> Rational
forall a. Num a => Polynomial a -> a -> a
evaluatePolynomial ([Rational] -> Polynomial Rational
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial ((Integer -> Rational) -> [Integer] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Polynomial Rational -> [Integer]
integerCoefficients Polynomial Rational
f))) Rational
x Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0
                ]
        attempt :: b
-> [[Rational]]
-> Either MathError (Maybe (Polynomial Rational), b)
attempt b
remain [] = b
-> Polynomial Rational
-> Int
-> Either MathError (Maybe (Polynomial Rational), b)
findFactor b
remain Polynomial Rational
f (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
        attempt b
0 [[Rational]]
_ = MathError -> Either MathError (Maybe (Polynomial Rational), b)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"Kronecker factorization" Int
budget)
        attempt b
remain ([Rational]
ys : [[Rational]]
rest) = do
            g <- [(Rational, Rational)] -> Either MathError (Polynomial Rational)
forall a.
(Eq a, Fractional a) =>
[(a, a)] -> Either MathError (Polynomial a)
lagrangeInterpolation ([Rational] -> [Rational] -> [(Rational, Rational)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Rational]
points [Rational]
ys)
            if deg g > 0 && deg g <= k && all ((== 1) . denominator) (coefficients g)
                then do
                    (_, r) <- dividePolynomial f g
                    if deg r < 0 then Right (Just g, remain - 1) else attempt (remain - 1) rest
                else attempt (remain - 1) rest

{- | Sturm chain over exact rationals, after removing the derivative GCD.
Zero polynomial is invalid. Square-free normalization ensures a repeated root
cannot make every chain member vanish and destroy one-sided variation counts.
-}
sturmSequence :: Polynomial Rational -> Either MathError [Polynomial Rational]
sturmSequence :: Polynomial Rational -> Either MathError [Polynomial Rational]
sturmSequence Polynomial Rational
p
    | Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError [Polynomial Rational]
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"Sturm chain of zero polynomial")
    | Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = [Polynomial Rational] -> Either MathError [Polynomial Rational]
forall a b. b -> Either a b
Right [Polynomial Rational
p]
    | Bool
otherwise = [Polynomial Rational] -> Either MathError [Polynomial Rational]
forall a b. b -> Either a b
Right (Polynomial Rational
squareFree Polynomial Rational
-> [Polynomial Rational] -> [Polynomial Rational]
forall a. a -> [a] -> [a]
: Polynomial Rational -> Polynomial Rational -> [Polynomial Rational]
forall {a}.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> [Polynomial a]
go Polynomial Rational
squareFree (Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Num a) => Polynomial a -> Polynomial a
differentiatePolynomial Polynomial Rational
squareFree))
  where
    squareFree :: Polynomial Rational
squareFree = Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
quotient Polynomial Rational
p (Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
gcdPolynomial Polynomial Rational
p (Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Num a) => Polynomial a -> Polynomial a
differentiatePolynomial Polynomial Rational
p))
    go :: Polynomial a -> Polynomial a -> [Polynomial a]
go Polynomial a
a Polynomial a
b
        | Polynomial a -> Int
forall a. Polynomial a -> Int
deg Polynomial a
b Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = []
        | Bool
otherwise =
            Polynomial a
b Polynomial a -> [Polynomial a] -> [Polynomial a]
forall a. a -> [a] -> [a]
: case Polynomial a
-> Polynomial a -> Either MathError (Polynomial a, Polynomial a)
forall a.
(Eq a, Fractional a) =>
Polynomial a
-> Polynomial a -> Either MathError (Polynomial a, Polynomial a)
dividePolynomial Polynomial a
a Polynomial a
b of
                Right (Polynomial a
_, Polynomial a
r) -> Polynomial a -> Polynomial a -> [Polynomial a]
go Polynomial a
b (a -> Polynomial a -> Polynomial a
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial (-a
1) Polynomial a
r)
                Left MathError
_ -> []

variations :: [Rational] -> Int
variations :: [Rational] -> Int
variations [Rational]
xs = [()] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [() | (Rational
a, Rational
b) <- [Rational] -> [Rational] -> [(Rational, Rational)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Rational]
nz (Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
drop Int
1 [Rational]
nz), Rational -> Rational
forall a. Num a => a -> a
signum Rational
a Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational -> Rational
forall a. Num a => a -> a
signum Rational
b]
  where
    nz :: [Rational]
nz = (Rational -> Bool) -> [Rational] -> [Rational]
forall a. (a -> Bool) -> [a] -> [a]
filter (Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0) [Rational]
xs

{- | Number of distinct real roots in @(a,b]@. Endpoints that are roots of the
polynomial are rejected to make the interval convention unambiguous.
-}
countRealRoots :: Polynomial Rational -> Rational -> Rational -> Either MathError Int
countRealRoots :: Polynomial Rational -> Rational -> Rational -> Either MathError Int
countRealRoots Polynomial Rational
p Rational
a Rational
b
    | Rational
a Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
b = MathError -> Either MathError Int
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"root-count interval must increase")
    | Polynomial Rational -> Rational -> Rational
forall a. Num a => Polynomial a -> a -> a
evaluatePolynomial Polynomial Rational
p Rational
a Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0 Bool -> Bool -> Bool
|| Polynomial Rational -> Rational -> Rational
forall a. Num a => Polynomial a -> a -> a
evaluatePolynomial Polynomial Rational
p Rational
b Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0 =
        MathError -> Either MathError Int
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"root-count endpoint is a root")
    | Bool
otherwise = do
        chain <- Polynomial Rational -> Either MathError [Polynomial Rational]
sturmSequence Polynomial Rational
p
        pure
            ( variations (map (`evaluatePolynomial` a) chain)
                - variations (map (`evaluatePolynomial` b) chain)
            )

{- | Isolate every distinct real root in rational intervals with disjoint interiors and width
at most tolerance. Intervals have the root-count convention (a,b]; shared endpoints
belong to the interval on their left.
Bisection budget bounds visited intervals; all arithmetic is exact.
-}
isolateRealRoots ::
    Int ->
    Rational ->
    Polynomial Rational ->
    Either MathError [(Rational, Rational)]
isolateRealRoots :: Int
-> Rational
-> Polynomial Rational
-> Either MathError [(Rational, Rational)]
isolateRealRoots Int
budget Rational
tolerance Polynomial Rational
p
    | Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError [(Rational, Rational)]
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"negative isolation budget")
    | Rational
tolerance Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
0 = MathError -> Either MathError [(Rational, Rational)]
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError [(Rational, Rational)]
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"zero polynomial has infinitely many roots")
    | Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = [(Rational, Rational)] -> Either MathError [(Rational, Rational)]
forall a b. b -> Either a b
Right []
    | Bool
otherwise = do
        let f :: Polynomial Rational
f = Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Fractional a) => Polynomial a -> Polynomial a
monicPolynomial Polynomial Rational
p
            bound :: Rational
bound = Rational
1 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ [Rational] -> Rational
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Rational
0 Rational -> [Rational] -> [Rational]
forall a. a -> [a] -> [a]
: (Rational -> Rational) -> [Rational] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map Rational -> Rational
forall a. Num a => a -> a
abs (Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
take (Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
f) (Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
coefficients Polynomial Rational
f)))
        chain <- Polynomial Rational -> Either MathError [Polynomial Rational]
sturmSequence Polynomial Rational
f
        fst <$> split chain budget (-bound) bound
  where
    var :: [Polynomial Rational] -> Rational -> Int
var [Polynomial Rational]
chain Rational
x = [Rational] -> Int
variations ((Polynomial Rational -> Rational)
-> [Polynomial Rational] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Polynomial Rational -> Rational -> Rational
forall a. Num a => Polynomial a -> a -> a
`evaluatePolynomial` Rational
x) [Polynomial Rational]
chain)
    split :: [Polynomial Rational]
-> t
-> Rational
-> Rational
-> Either MathError ([(Rational, Rational)], t)
split [Polynomial Rational]
_ t
0 Rational
_ Rational
_ = MathError -> Either MathError ([(Rational, Rational)], t)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"Sturm isolation" Int
budget)
    split [Polynomial Rational]
chain t
left Rational
a Rational
b
        | Int
count Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = ([(Rational, Rational)], t)
-> Either MathError ([(Rational, Rational)], t)
forall a b. b -> Either a b
Right ([], t
left t -> t -> t
forall a. Num a => a -> a -> a
- t
1)
        | Int
count Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 Bool -> Bool -> Bool
&& Rational
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
a Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
tolerance = ([(Rational, Rational)], t)
-> Either MathError ([(Rational, Rational)], t)
forall a b. b -> Either a b
Right ([(Rational
a, Rational
b)], t
left t -> t -> t
forall a. Num a => a -> a -> a
- t
1)
        | Bool
otherwise = do
            let mid :: Rational
mid = (Rational
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
b) Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
2
            (lo, l) <- [Polynomial Rational]
-> t
-> Rational
-> Rational
-> Either MathError ([(Rational, Rational)], t)
split [Polynomial Rational]
chain (t
left t -> t -> t
forall a. Num a => a -> a -> a
- t
1) Rational
a Rational
mid
            (hi, r) <- split chain l mid b
            pure (lo ++ hi, r)
      where
        count :: Int
count = [Polynomial Rational] -> Rational -> Int
var [Polynomial Rational]
chain Rational
a Int -> Int -> Int
forall a. Num a => a -> a -> a
- [Polynomial Rational] -> Rational -> Int
var [Polynomial Rational]
chain Rational
b

{- | Durand--Kerner simultaneous complex root iteration. Distinct initial roots
are placed on a Cauchy-bound circle; repeated roots may converge slowly.
-}
durandKerner ::
    Int ->
    Double ->
    Polynomial (Complex Double) ->
    Either MathError [Complex Double]
durandKerner :: Int
-> Double
-> Polynomial (Complex Double)
-> Either MathError [Complex Double]
durandKerner = Bool
-> Int
-> Double
-> Polynomial (Complex Double)
-> Either MathError [Complex Double]
complexRoots Bool
False

{- | Aberth--Ehrlich simultaneous iteration, combining Newton corrections with
mutual-root repulsion. All roots are verified against a scaled residual.
-}
aberthEhrlich ::
    Int ->
    Double ->
    Polynomial (Complex Double) ->
    Either MathError [Complex Double]
aberthEhrlich :: Int
-> Double
-> Polynomial (Complex Double)
-> Either MathError [Complex Double]
aberthEhrlich = Bool
-> Int
-> Double
-> Polynomial (Complex Double)
-> Either MathError [Complex Double]
complexRoots Bool
True

complexRoots ::
    Bool ->
    Int ->
    Double ->
    Polynomial (Complex Double) ->
    Either MathError [Complex Double]
complexRoots :: Bool
-> Int
-> Double
-> Polynomial (Complex Double)
-> Either MathError [Complex Double]
complexRoots Bool
useAberth Int
budget Double
tolerance Polynomial (Complex Double)
p
    | Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError [Complex Double]
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"negative root iteration budget")
    | Double
tolerance Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
tolerance Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
tolerance = MathError -> Either MathError [Complex Double]
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | (Complex Double -> Bool) -> [Complex Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Bool -> Bool
not (Bool -> Bool)
-> (Complex Double -> Bool) -> Complex Double -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Complex Double -> Bool
finiteComplex) (Polynomial (Complex Double) -> [Complex Double]
forall a. Polynomial a -> [a]
coefficients Polynomial (Complex Double)
p) = MathError -> Either MathError [Complex Double]
forall a b. a -> Either a b
Left MathError
NonFiniteResult
    | Polynomial (Complex Double) -> Int
forall a. Polynomial a -> Int
deg Polynomial (Complex Double)
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError [Complex Double]
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"zero polynomial has infinitely many roots")
    | Polynomial (Complex Double) -> Int
forall a. Polynomial a -> Int
deg Polynomial (Complex Double)
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = [Complex Double] -> Either MathError [Complex Double]
forall a b. b -> Either a b
Right []
    | (Complex Double -> Bool) -> [Complex Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Bool -> Bool
not (Bool -> Bool)
-> (Complex Double -> Bool) -> Complex Double -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Complex Double -> Bool
finiteComplex) (Polynomial (Complex Double) -> [Complex Double]
forall a. Polynomial a -> [a]
coefficients Polynomial (Complex Double)
f) = MathError -> Either MathError [Complex Double]
forall a b. a -> Either a b
Left MathError
NonFiniteResult
    | Bool
otherwise = Int -> [Complex Double] -> Either MathError [Complex Double]
forall {t}.
(Eq t, Num t) =>
t -> [Complex Double] -> Either MathError [Complex Double]
go Int
budget [Complex Double]
initial
  where
    f :: Polynomial (Complex Double)
f = Polynomial (Complex Double) -> Polynomial (Complex Double)
forall a. (Eq a, Fractional a) => Polynomial a -> Polynomial a
monicPolynomial Polynomial (Complex Double)
p
    n :: Int
n = Polynomial (Complex Double) -> Int
forall a. Polynomial a -> Int
deg Polynomial (Complex Double)
f
    df :: Polynomial (Complex Double)
df = Polynomial (Complex Double) -> Polynomial (Complex Double)
forall a. (Eq a, Num a) => Polynomial a -> Polynomial a
differentiatePolynomial Polynomial (Complex Double)
f
    radius :: Double
radius = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Complex Double -> Double) -> [Complex Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude (Int -> [Complex Double] -> [Complex Double]
forall a. Int -> [a] -> [a]
take Int
n (Polynomial (Complex Double) -> [Complex Double]
forall a. Polynomial a -> [a]
coefficients Polynomial (Complex Double)
f)))
    initial :: [Complex Double]
initial = [Double -> Double -> Complex Double
forall a. Floating a => a -> a -> Complex a
mkPolar Double
radius (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.271828) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    residual :: Complex Double -> Bool
residual Complex Double
z =
        let coefficientScale :: Double
coefficientScale = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
1 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Complex Double -> Double) -> [Complex Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude (Polynomial (Complex Double) -> [Complex Double]
forall a. Polynomial a -> [a]
coefficients Polynomial (Complex Double)
f))
            scaled :: [Complex Double]
scaled = (Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ (Double
coefficientScale Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0)) (Polynomial (Complex Double) -> [Complex Double]
forall a. Polynomial a -> [a]
coefficients Polynomial (Complex Double)
f)
            reverseEvaluation :: Bool
reverseEvaluation = Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude Complex Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1
            point :: Complex Double
point = if Bool
reverseEvaluation then Complex Double
1 Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ Complex Double
z else Complex Double
z
            cs :: [Complex Double]
cs = if Bool
reverseEvaluation then [Complex Double] -> [Complex Double]
forall a. [a] -> [a]
reverse [Complex Double]
scaled else [Complex Double]
scaled
            value :: Double
value = Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude (Polynomial (Complex Double) -> Complex Double -> Complex Double
forall a. Num a => Polynomial a -> a -> a
evaluatePolynomial ([Complex Double] -> Polynomial (Complex Double)
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [Complex Double]
cs) Complex Double
point)
            bound :: Double
bound = (Complex Double -> Double -> Double)
-> Double -> [Complex Double] -> Double
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\Complex Double
a Double
b -> Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude Complex Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude Complex Double
point Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) Double
0 [Complex Double]
cs
         in (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
x -> Bool -> Bool
not (Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x)) [Double
value, Double
bound] Bool -> Bool -> Bool
&& Double
value Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
tolerance Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
bound
    go :: t -> [Complex Double] -> Either MathError [Complex Double]
go t
_ [Complex Double]
zs | (Complex Double -> Bool) -> [Complex Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Complex Double -> Bool
residual [Complex Double]
zs = [Complex Double] -> Either MathError [Complex Double]
forall a b. b -> Either a b
Right [Complex Double]
zs
    go t
0 [Complex Double]
_ = MathError -> Either MathError [Complex Double]
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence (if Bool
useAberth then String
"Aberth-Ehrlich" else String
"Durand-Kerner") Int
budget)
    go t
left [Complex Double]
zs = do
        next <- ((Int, Complex Double) -> Either MathError (Complex Double))
-> [(Int, Complex Double)] -> Either MathError [Complex Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM ([Complex Double]
-> (Int, Complex Double) -> Either MathError (Complex Double)
update [Complex Double]
zs) ([Int] -> [Complex Double] -> [(Int, Complex Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Complex Double]
zs)
        if all finiteComplex next then go (left - 1) next else Left NonFiniteResult
    update :: [Complex Double]
-> (Int, Complex Double) -> Either MathError (Complex Double)
update [Complex Double]
zs (Int
i, Complex Double
z) =
        let others :: [Complex Double]
others = [Complex Double
w | (Int
j, Complex Double
w) <- [Int] -> [Complex Double] -> [(Int, Complex Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip ([Int
0 ..] :: [Int]) [Complex Double]
zs, Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i]
            value :: Complex Double
value = Polynomial (Complex Double) -> Complex Double -> Complex Double
forall a. Num a => Polynomial a -> a -> a
evaluatePolynomial Polynomial (Complex Double)
f Complex Double
z
            denom :: Complex Double
denom =
                if Bool
useAberth
                    then
                        Polynomial (Complex Double) -> Complex Double -> Complex Double
forall a. Num a => Polynomial a -> a -> a
evaluatePolynomial Polynomial (Complex Double)
df Complex Double
z
                            Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
- Complex Double
value Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* [Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Complex Double
1 Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ (Complex Double
z Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
- Complex Double
w) | Complex Double
w <- [Complex Double]
others]
                    else [Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Complex Double
z Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
- Complex Double
w | Complex Double
w <- [Complex Double]
others]
         in if Complex Double
denom Complex Double -> Complex Double -> Bool
forall a. Eq a => a -> a -> Bool
== Complex Double
0
                then MathError -> Either MathError (Complex Double)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"colliding polynomial roots" Int
budget)
                else Complex Double -> Either MathError (Complex Double)
forall a b. b -> Either a b
Right (Complex Double
z Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
- Complex Double
value Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ Complex Double
denom)

finiteComplex :: Complex Double -> Bool
finiteComplex :: Complex Double -> Bool
finiteComplex (Double
a :+ Double
b) = (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
x -> Bool -> Bool
not (Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x)) [Double
a, Double
b]

{- | Frobenius companion matrix of a nonconstant polynomial, with ones on the
subdiagonal and negated monic coefficients in the last column.
-}
companionMatrix :: (Eq a, Fractional a) => Polynomial a -> Either MathError (M.Matrix a)
companionMatrix :: forall a.
(Eq a, Fractional a) =>
Polynomial a -> Either MathError (Matrix a)
companionMatrix Polynomial a
p
    | Polynomial a -> Int
forall a. Polynomial a -> Int
deg Polynomial a
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"companion matrix needs positive degree")
    | Bool
otherwise =
        Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked
            Int
n
            Int
n
            [Int -> Int -> a
entry Int
i Int
j | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
  where
    n :: Int
n = Polynomial a -> Int
forall a. Polynomial a -> Int
deg Polynomial a
p
    cs :: [a]
cs = Polynomial a -> [a]
forall a. Polynomial a -> [a]
coefficients (Polynomial a -> Polynomial a
forall a. (Eq a, Fractional a) => Polynomial a -> Polynomial a
monicPolynomial Polynomial a
p)
    entry :: Int -> Int -> a
entry Int
i Int
j
        | Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1 = -([a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
1 (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
i [a]
cs)))
        | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 = a
1
        | Bool
otherwise = a
0

-- | Exact Lagrange interpolating polynomial. Duplicate nodes are rejected.
lagrangeInterpolation :: (Eq a, Fractional a) => [(a, a)] -> Either MathError (Polynomial a)
lagrangeInterpolation :: forall a.
(Eq a, Fractional a) =>
[(a, a)] -> Either MathError (Polynomial a)
lagrangeInterpolation [(a, a)]
points
    | [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([a] -> [a]
forall a. Eq a => [a] -> [a]
nub (((a, a) -> a) -> [(a, a)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a, a) -> a
forall a b. (a, b) -> a
fst [(a, a)]
points)) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [(a, a)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(a, a)]
points =
        MathError -> Either MathError (Polynomial a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"duplicate interpolation nodes")
    | Bool
otherwise =
        Polynomial a -> Either MathError (Polynomial a)
forall a b. b -> Either a b
Right
            ( (Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a -> [Polynomial a] -> Polynomial a
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
                Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
addPolynomial
                ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [])
                [ a -> Polynomial a -> Polynomial a
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial
                    (a
y a -> a -> a
forall a. Fractional a => a -> a -> a
/ [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
z | (Int
j, (a
z, a
_)) <- [(Int, (a, a))]
indexed, Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i])
                    ( (Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a -> [Polynomial a] -> Polynomial a
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
                        Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial
                        ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
1])
                        [[a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [-a
z, a
1] | (Int
j, (a
z, a
_)) <- [(Int, (a, a))]
indexed, Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i]
                    )
                | (Int
i, (a
x, a
y)) <- [(Int, (a, a))]
indexed
                ]
            )
  where
    indexed :: [(Int, (a, a))]
indexed = [Int] -> [(a, a)] -> [(Int, (a, a))]
forall a b. [a] -> [b] -> [(a, b)]
zip ([Int
0 ..] :: [Int]) [(a, a)]
points

-- | Newton divided differences, returned in the standard power basis.
newtonInterpolation :: (Eq a, Fractional a) => [(a, a)] -> Either MathError (Polynomial a)
newtonInterpolation :: forall a.
(Eq a, Fractional a) =>
[(a, a)] -> Either MathError (Polynomial a)
newtonInterpolation [(a, a)]
points
    | [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([a] -> [a]
forall a. Eq a => [a] -> [a]
nub [a]
xs) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs = MathError -> Either MathError (Polynomial a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"duplicate interpolation nodes")
    | Bool
otherwise = Polynomial a -> Either MathError (Polynomial a)
forall a b. b -> Either a b
Right (Int -> [a] -> Polynomial a -> Polynomial a -> Polynomial a
go Int
0 (((a, a) -> a) -> [(a, a)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a, a) -> a
forall a b. (a, b) -> b
snd [(a, a)]
points) ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
1]) ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial []))
  where
    xs :: [a]
xs = ((a, a) -> a) -> [(a, a)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a, a) -> a
forall a b. (a, b) -> a
fst [(a, a)]
points
    go :: Int -> [a] -> Polynomial a -> Polynomial a -> Polynomial a
go Int
_ [] Polynomial a
_ Polynomial a
acc = Polynomial a
acc
    go Int
k [a]
ys Polynomial a
basis Polynomial a
acc =
        let lead :: a
lead = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
1 [a]
ys)
            result :: Polynomial a
result = Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
addPolynomial Polynomial a
acc (a -> Polynomial a -> Polynomial a
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial a
lead Polynomial a
basis)
            next :: [a]
next = (a -> a -> (a, a) -> a) -> [a] -> [a] -> [(a, a)] -> [a]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\a
a a
b (a
x, a
z) -> (a
b a -> a -> a
forall a. Num a => a -> a -> a
- a
a) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
z a -> a -> a
forall a. Num a => a -> a -> a
- a
x)) [a]
ys (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
1 [a]
ys) ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
xs (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
xs))
            node :: a
node = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
1 (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
k [a]
xs))
         in Int -> [a] -> Polynomial a -> Polynomial a -> Polynomial a
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
next (Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial Polynomial a
basis ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [-a
node, a
1])) Polynomial a
result

orthogonal ::
    (Fractional a) =>
    Int ->
    (Int -> Polynomial a -> Polynomial a -> Polynomial a) ->
    Polynomial a ->
    Polynomial a ->
    Either MathError (Polynomial a)
orthogonal :: forall a.
Fractional a =>
Int
-> (Int -> Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a
-> Polynomial a
-> Either MathError (Polynomial a)
orthogonal Int
n Int -> Polynomial a -> Polynomial a -> Polynomial a
step Polynomial a
p0 Polynomial a
p1
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError (Polynomial a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"negative polynomial degree")
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = Polynomial a -> Either MathError (Polynomial a)
forall a b. b -> Either a b
Right Polynomial a
p0
    | Bool
otherwise = Polynomial a -> Either MathError (Polynomial a)
forall a b. b -> Either a b
Right (Int -> Polynomial a -> Polynomial a -> Polynomial a
go Int
1 Polynomial a
p0 Polynomial a
p1)
  where
    go :: Int -> Polynomial a -> Polynomial a -> Polynomial a
go Int
k Polynomial a
a Polynomial a
b
        | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n = Polynomial a
b
        | Bool
otherwise = Int -> Polynomial a -> Polynomial a -> Polynomial a
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Polynomial a
b (Int -> Polynomial a -> Polynomial a -> Polynomial a
step Int
k Polynomial a
a Polynomial a
b)

-- | Chebyshev T_n, normalized by T_0=1, T_1=x.
chebyshevPolynomial :: (Eq a, Fractional a) => Int -> Either MathError (Polynomial a)
chebyshevPolynomial :: forall a.
(Eq a, Fractional a) =>
Int -> Either MathError (Polynomial a)
chebyshevPolynomial Int
n =
    Int
-> (Int -> Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a
-> Polynomial a
-> Either MathError (Polynomial a)
forall a.
Fractional a =>
Int
-> (Int -> Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a
-> Polynomial a
-> Either MathError (Polynomial a)
orthogonal
        Int
n
        ( \Int
_ Polynomial a
a Polynomial a
b ->
            Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
subtractPolynomial
                (Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
0, a
2]) Polynomial a
b)
                Polynomial a
a
        )
        ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
1])
        ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
0, a
1])

-- | Legendre P_n, normalized by P_n(1)=1.
legendrePolynomial :: (Eq a, Fractional a) => Int -> Either MathError (Polynomial a)
legendrePolynomial :: forall a.
(Eq a, Fractional a) =>
Int -> Either MathError (Polynomial a)
legendrePolynomial Int
n =
    Int
-> (Int -> Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a
-> Polynomial a
-> Either MathError (Polynomial a)
forall a.
Fractional a =>
Int
-> (Int -> Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a
-> Polynomial a
-> Either MathError (Polynomial a)
orthogonal
        Int
n
        ( \Int
k Polynomial a
a Polynomial a
b ->
            a -> Polynomial a -> Polynomial a
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial
                (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))
                ( Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
subtractPolynomial
                    (Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
0, Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)]) Polynomial a
b)
                    (a -> Polynomial a -> Polynomial a
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial (Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k) Polynomial a
a)
                )
        )
        ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
1])
        ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
0, a
1])

-- | Physicists' Hermite H_n: H_0=1, H_1=2x.
hermitePolynomial :: (Eq a, Fractional a) => Int -> Either MathError (Polynomial a)
hermitePolynomial :: forall a.
(Eq a, Fractional a) =>
Int -> Either MathError (Polynomial a)
hermitePolynomial Int
n =
    Int
-> (Int -> Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a
-> Polynomial a
-> Either MathError (Polynomial a)
forall a.
Fractional a =>
Int
-> (Int -> Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a
-> Polynomial a
-> Either MathError (Polynomial a)
orthogonal
        Int
n
        ( \Int
k Polynomial a
a Polynomial a
b ->
            Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
subtractPolynomial
                (Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
0, a
2]) Polynomial a
b)
                (a -> Polynomial a -> Polynomial a
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial (Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k)) Polynomial a
a)
        )
        ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
1])
        ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
0, a
2])

-- | Ordinary Laguerre L_n with alpha=0 and L_n(0)=1.
laguerrePolynomial :: (Eq a, Fractional a) => Int -> Either MathError (Polynomial a)
laguerrePolynomial :: forall a.
(Eq a, Fractional a) =>
Int -> Either MathError (Polynomial a)
laguerrePolynomial Int
n =
    Int
-> (Int -> Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a
-> Polynomial a
-> Either MathError (Polynomial a)
forall a.
Fractional a =>
Int
-> (Int -> Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a
-> Polynomial a
-> Either MathError (Polynomial a)
orthogonal
        Int
n
        ( \Int
k Polynomial a
a Polynomial a
b ->
            a -> Polynomial a -> Polynomial a
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial
                (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))
                ( Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
subtractPolynomial
                    (Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1), -a
1]) Polynomial a
b)
                    (a -> Polynomial a -> Polynomial a
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial (Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k) Polynomial a
a)
                )
        )
        ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
1])
        ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
1, -a
1])

-- | Bernstein basis @binomial(n,k)*x^k*(1-x)^(n-k)@ for 0<=k<=n.
bernsteinPolynomial :: (Eq a, Fractional a) => Int -> Int -> Either MathError (Polynomial a)
bernsteinPolynomial :: forall a.
(Eq a, Fractional a) =>
Int -> Int -> Either MathError (Polynomial a)
bernsteinPolynomial Int
n Int
k
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n = MathError -> Either MathError (Polynomial a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"invalid Bernstein indices")
    | Bool
otherwise =
        Polynomial a -> Either MathError (Polynomial a)
forall a b. b -> Either a b
Right
            ( a -> Polynomial a -> Polynomial a
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial
                (Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
choose)
                ( Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial
                    ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial (Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
k a
0 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
1]))
                    (Polynomial a -> Int -> Polynomial a
forall {a}. (Eq a, Num a) => Polynomial a -> Int -> Polynomial a
power ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
1, -a
1]) (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k))
                )
            )
  where
    choose :: Integer
choose = [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
n 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 -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n] Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
1 .. Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
k]
    power :: Polynomial a -> Int -> Polynomial a
power Polynomial a
p Int
j = (Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a -> [Polynomial a] -> Polynomial a
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
1]) (Int -> Polynomial a -> [Polynomial a]
forall a. Int -> a -> [a]
replicate Int
j Polynomial a
p)

{- | Least-squares power-basis regression using column-pivoted Householder QR.
Degree must be nonnegative; undersampling, nonfinite data and rank deficiency
are reported by the checked design-matrix solver.
-}
polynomialRegression :: (RealFloat a) => Int -> [(a, a)] -> Either MathError (Polynomial a)
polynomialRegression :: forall a.
RealFloat a =>
Int -> [(a, a)] -> Either MathError (Polynomial a)
polynomialRegression Int
degree [(a, a)]
points
    | Int
degree Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError (Polynomial a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"negative regression degree")
    | Int
degree Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
forall a. Bounded a => a
maxBound = MathError -> Either MathError (Polynomial a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"regression degree too large")
    | Bool
otherwise = do
        matrix <-
            Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked
                ([(a, a)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(a, a)]
points)
                (Int
degree Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                (((a, a) -> [a]) -> [(a, a)] -> [a]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
degree Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) ([a] -> [a]) -> ((a, a) -> [a]) -> (a, a) -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> [a]
forall {a}. Num a => a -> [a]
iterateMultiply (a -> [a]) -> ((a, a) -> a) -> (a, a) -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a, a) -> a
forall a b. (a, b) -> a
fst) [(a, a)]
points)
        polynomial <$> leastSquares matrix (map snd points)
  where
    iterateMultiply :: a -> [a]
iterateMultiply a
x = (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. Num a => a -> a -> a
* a
x) a
1

{- | Truncated Chebyshev interpolation at first-kind nodes on [-1,1], returned
in power basis. Uses degree+1 evaluations; nonfinite values are rejected.
-}
chebyshevApproximation :: Int -> (Double -> Double) -> Either MathError (Polynomial Double)
chebyshevApproximation :: Int -> (Double -> Double) -> Either MathError (Polynomial Double)
chebyshevApproximation Int
degree Double -> Double
f
    | Int
degree Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
degree Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
forall a. Bounded a => a
maxBound = MathError -> Either MathError (Polynomial Double)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"invalid approximation degree")
    | (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\Double
x -> Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x) [Double]
ys = MathError -> Either MathError (Polynomial Double)
forall a b. a -> Either a b
Left MathError
NonFiniteResult
    | Bool
otherwise = do
        basis <- (Int -> Either MathError (Polynomial Double))
-> [Int] -> Either MathError [Polynomial Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM Int -> Either MathError (Polynomial Double)
forall a.
(Eq a, Fractional a) =>
Int -> Either MathError (Polynomial a)
chebyshevPolynomial [Int
0 .. Int
degree]
        pure (foldl addPolynomial (polynomial []) (zipWith scalePolynomial weights basis))
  where
    n :: Int
n = Int
degree Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
    angles :: [Double]
angles = [Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.5) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n | Int
j <- [Int
0 .. Int
degree]]
    ys :: [Double]
ys = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double
f (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
forall a. Floating a => a -> a
cos) [Double]
angles
    weights :: [Double]
weights =
        [ (if Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Double
1 else Double
2)
            Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
            Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
y Double
theta -> Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
theta)) [Double]
ys [Double]
angles)
        | Int
k <- [Int
0 .. Int
degree]
        ]

{- | A finite Taylor jet, coefficients ascending and zeros retained. Its length
is the truncation order; operations preserve the shorter operand's order.
-}
newtype PowerSeries a = PowerSeries [a] deriving (PowerSeries a -> PowerSeries a -> Bool
(PowerSeries a -> PowerSeries a -> Bool)
-> (PowerSeries a -> PowerSeries a -> Bool) -> Eq (PowerSeries a)
forall a. Eq a => PowerSeries a -> PowerSeries a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => PowerSeries a -> PowerSeries a -> Bool
== :: PowerSeries a -> PowerSeries a -> Bool
$c/= :: forall a. Eq a => PowerSeries a -> PowerSeries a -> Bool
/= :: PowerSeries a -> PowerSeries a -> Bool
Eq, Int -> PowerSeries a -> ShowS
[PowerSeries a] -> ShowS
PowerSeries a -> String
(Int -> PowerSeries a -> ShowS)
-> (PowerSeries a -> String)
-> ([PowerSeries a] -> ShowS)
-> Show (PowerSeries a)
forall a. Show a => Int -> PowerSeries a -> ShowS
forall a. Show a => [PowerSeries a] -> ShowS
forall a. Show a => PowerSeries a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> PowerSeries a -> ShowS
showsPrec :: Int -> PowerSeries a -> ShowS
$cshow :: forall a. Show a => PowerSeries a -> String
show :: PowerSeries a -> String
$cshowList :: forall a. Show a => [PowerSeries a] -> ShowS
showList :: [PowerSeries a] -> ShowS
Show)

-- | Construct a jet with explicit coefficients, including trailing zeros.
powerSeries :: [a] -> PowerSeries a
powerSeries :: forall a. [a] -> PowerSeries a
powerSeries = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PowerSeries

-- | Stored Taylor coefficients; coefficient n multiplies x^n, not x^n/n!.
seriesCoefficients :: PowerSeries a -> [a]
seriesCoefficients :: forall a. PowerSeries a -> [a]
seriesCoefficients (PowerSeries [a]
xs) = [a]
xs

-- | Add jets to their common truncation order.
addSeries :: (Num a) => PowerSeries a -> PowerSeries a -> PowerSeries a
addSeries :: forall a. Num a => PowerSeries a -> PowerSeries a -> PowerSeries a
addSeries (PowerSeries [a]
a) (PowerSeries [a]
b) = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PowerSeries ((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]
a [a]
b)

-- | Multiply jets by truncated convolution.
multiplySeries :: (Num a) => PowerSeries a -> PowerSeries a -> PowerSeries a
multiplySeries :: forall a. Num a => PowerSeries a -> PowerSeries a -> PowerSeries a
multiplySeries (PowerSeries [a]
a) (PowerSeries [a]
b) =
    [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PowerSeries
        [ [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((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
(*) (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
a) ([a] -> [a]
forall a. [a] -> [a]
reverse (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
b)))
        | Int
k <- [Int
0 .. Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
a) ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
b) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
        ]

-- | Reciprocal jet from coefficient recurrence; constant term must be nonzero.
reciprocalSeries :: (Eq a, Fractional a) => PowerSeries a -> Either MathError (PowerSeries a)
reciprocalSeries :: forall a.
(Eq a, Fractional a) =>
PowerSeries a -> Either MathError (PowerSeries a)
reciprocalSeries (PowerSeries []) = MathError -> Either MathError (PowerSeries a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"empty reciprocal series")
reciprocalSeries (PowerSeries (a
a : [a]
as))
    | a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = MathError -> Either MathError (PowerSeries a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"zero constant reciprocal series")
    | Bool
otherwise = PowerSeries a -> Either MathError (PowerSeries a)
forall a b. b -> Either a b
Right ([a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PowerSeries [a]
result)
  where
    result :: [a]
result = a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
a a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [-[a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((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
(*) (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
k [a]
as) ([a] -> [a]
forall a. [a] -> [a]
reverse (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
k [a]
result))) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
a | Int
k <- [Int
1 .. [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
as]]

-- | Divide jets to their common order; denominator constant must be nonzero.
divideSeries ::
    (Eq a, Fractional a) => PowerSeries a -> PowerSeries a -> Either MathError (PowerSeries a)
divideSeries :: forall a.
(Eq a, Fractional a) =>
PowerSeries a -> PowerSeries a -> Either MathError (PowerSeries a)
divideSeries PowerSeries a
a PowerSeries a
b = PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Num a => PowerSeries a -> PowerSeries a -> PowerSeries a
multiplySeries PowerSeries a
a (PowerSeries a -> PowerSeries a)
-> Either MathError (PowerSeries a)
-> Either MathError (PowerSeries a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> PowerSeries a -> Either MathError (PowerSeries a)
forall a.
(Eq a, Fractional a) =>
PowerSeries a -> Either MathError (PowerSeries a)
reciprocalSeries PowerSeries a
b

{- | Formal composition to the inner jet's order. Inner constant must be zero,
so unknown higher outer coefficients cannot affect the returned jet.
-}
composeSeries :: (Eq a, Num a) => PowerSeries a -> PowerSeries a -> Either MathError (PowerSeries a)
composeSeries :: forall a.
(Eq a, Num a) =>
PowerSeries a -> PowerSeries a -> Either MathError (PowerSeries a)
composeSeries (PowerSeries [a]
outer) inner :: PowerSeries a
inner@(PowerSeries [a]
xs)
    | (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0) (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
1 [a]
xs) = MathError -> Either MathError (PowerSeries a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"formal composition needs zero inner constant")
    | [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
outer Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs = MathError -> Either MathError (PowerSeries a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"outer jet order is too short")
    | Bool
otherwise = PowerSeries a -> Either MathError (PowerSeries a)
forall a b. b -> Either a b
Right ((a -> PowerSeries a -> PowerSeries a)
-> PowerSeries a -> [a] -> PowerSeries 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 -> PowerSeries a -> PowerSeries a
step ([a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PowerSeries (Int -> a -> [a]
forall a. Int -> a -> [a]
replicate ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs) a
0)) (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs) [a]
outer))
  where
    step :: a -> PowerSeries a -> PowerSeries a
step a
a PowerSeries a
acc = PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Num a => PowerSeries a -> PowerSeries a -> PowerSeries a
addSeries ([a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PowerSeries (a
a a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> a -> [a]
forall a. Int -> a -> [a]
replicate (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 ([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)) a
0)) (PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Num a => PowerSeries a -> PowerSeries a -> PowerSeries a
multiplySeries PowerSeries a
inner PowerSeries a
acc)

-- | Derivative jet has one fewer coefficient.
differentiateSeries :: (Num a) => PowerSeries a -> PowerSeries a
differentiateSeries :: forall a. Num a => PowerSeries a -> PowerSeries a
differentiateSeries (PowerSeries [a]
xs) = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PowerSeries ((Integer -> a -> a) -> [Integer] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Integer
k a
a -> Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
k a -> a -> a
forall a. Num a => a -> a -> a
* a
a) [Integer
1 ..] (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
1 [a]
xs))

-- | Integral jet with supplied constant, one coefficient longer.
integrateSeries :: (Fractional a) => a -> PowerSeries a -> PowerSeries a
integrateSeries :: forall a. Fractional a => a -> PowerSeries a -> PowerSeries a
integrateSeries a
c (PowerSeries [a]
xs) = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PowerSeries (a
c a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (Integer -> a -> a) -> [Integer] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Integer
k a
a -> a
a a -> a -> a
forall a. Fractional a => a -> a -> a
/ Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
k) [Integer
1 ..] [a]
xs)

{- | Exponential of a finite jet from b'=a'b; constant uses the coefficient
type's exp. The output has the same truncation order.
-}
expSeries :: (Floating a) => PowerSeries a -> PowerSeries a
expSeries :: forall a. Floating a => PowerSeries a -> PowerSeries a
expSeries (PowerSeries []) = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PowerSeries []
expSeries (PowerSeries (a
a : [a]
as)) = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PowerSeries [a]
bs
  where
    bs :: [a]
bs =
        a -> a
forall a. Floating a => a -> a
exp a
a
            a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [ [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                    [ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k a -> a -> a
forall a. Num a => a -> a -> a
* a
ak a -> a -> a
forall a. Num a => a -> a -> a
* a
bj
                    | (Int
k, a
ak) <- [Int] -> [a] -> [(Int, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1 .. Int
n] (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
n [a]
as)
                    , let bj :: a
bj = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
1 (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k) [a]
bs))
                    ]
                    a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
              | Int
n <- [Int
1 .. [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
as]
              ]

{- | Formal log through derivative/division/integration. Requires positive real
constant term, with exact truncation matching the input jet.
-}
logSeries :: (Ord a, Floating a) => PowerSeries a -> Either MathError (PowerSeries a)
logSeries :: forall a.
(Ord a, Floating a) =>
PowerSeries a -> Either MathError (PowerSeries a)
logSeries (PowerSeries []) = MathError -> Either MathError (PowerSeries a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"empty logarithm series")
logSeries a :: PowerSeries a
a@(PowerSeries (a
c : [a]
_))
    | a
c a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 = MathError -> Either MathError (PowerSeries a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"log series needs positive constant")
    | Bool
otherwise = do
        inv <- PowerSeries a -> Either MathError (PowerSeries a)
forall a.
(Eq a, Fractional a) =>
PowerSeries a -> Either MathError (PowerSeries a)
reciprocalSeries PowerSeries a
a
        let deriv = PowerSeries a -> PowerSeries a
forall a. Num a => PowerSeries a -> PowerSeries a
differentiateSeries PowerSeries a
a
            ratio = PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Num a => PowerSeries a -> PowerSeries a -> PowerSeries a
multiplySeries PowerSeries a
deriv PowerSeries a
inv
        pure (integrateSeries (log c) ratio)

-- | Reduced exact rational function with monic denominator.
data RationalFunction = RationalFunction (Polynomial Rational) (Polynomial Rational)
    deriving (RationalFunction -> RationalFunction -> Bool
(RationalFunction -> RationalFunction -> Bool)
-> (RationalFunction -> RationalFunction -> Bool)
-> Eq RationalFunction
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RationalFunction -> RationalFunction -> Bool
== :: RationalFunction -> RationalFunction -> Bool
$c/= :: RationalFunction -> RationalFunction -> Bool
/= :: RationalFunction -> RationalFunction -> Bool
Eq, Int -> RationalFunction -> ShowS
[RationalFunction] -> ShowS
RationalFunction -> String
(Int -> RationalFunction -> ShowS)
-> (RationalFunction -> String)
-> ([RationalFunction] -> ShowS)
-> Show RationalFunction
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> RationalFunction -> ShowS
showsPrec :: Int -> RationalFunction -> ShowS
$cshow :: RationalFunction -> String
show :: RationalFunction -> String
$cshowList :: [RationalFunction] -> ShowS
showList :: [RationalFunction] -> ShowS
Show)

-- | Normalize by polynomial GCD and reject a zero denominator.
rationalFunction :: Polynomial Rational -> Polynomial Rational -> Either MathError RationalFunction
rationalFunction :: Polynomial Rational
-> Polynomial Rational -> Either MathError RationalFunction
rationalFunction Polynomial Rational
p Polynomial Rational
q
    | Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
q Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError RationalFunction
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"zero rational-function denominator")
    | Bool
otherwise =
        let g :: Polynomial Rational
g = Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
gcdPolynomial Polynomial Rational
p Polynomial Rational
q
            a :: Polynomial Rational
a = Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
quotient Polynomial Rational
p Polynomial Rational
g
            b :: Polynomial Rational
b = Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
quotient Polynomial Rational
q Polynomial Rational
g
            lead :: Rational
lead = (Rational -> Rational -> Rational)
-> Rational -> [Rational] -> Rational
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\Rational
_ Rational
x -> Rational
x) Rational
1 (Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
coefficients Polynomial Rational
b)
         in RationalFunction -> Either MathError RationalFunction
forall a b. b -> Either a b
Right (Polynomial Rational -> Polynomial Rational -> RationalFunction
RationalFunction (Rational -> Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
scalePolynomial (Rational
1 Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
lead) Polynomial Rational
a) (Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Fractional a) => Polynomial a -> Polynomial a
monicPolynomial Polynomial Rational
b))

-- | Reduced numerator.
rationalNumerator :: RationalFunction -> Polynomial Rational
rationalNumerator :: RationalFunction -> Polynomial Rational
rationalNumerator (RationalFunction Polynomial Rational
p Polynomial Rational
_) = Polynomial Rational
p

-- | Monic reduced denominator.
rationalDenominator :: RationalFunction -> Polynomial Rational
rationalDenominator :: RationalFunction -> Polynomial Rational
rationalDenominator (RationalFunction Polynomial Rational
_ Polynomial Rational
q) = Polynomial Rational
q

-- | Exact evaluation, rejecting poles.
evaluateRationalFunction :: RationalFunction -> Rational -> Either MathError Rational
evaluateRationalFunction :: RationalFunction -> Rational -> Either MathError Rational
evaluateRationalFunction (RationalFunction Polynomial Rational
p Polynomial Rational
q) Rational
x
    | Rational
den Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0 = MathError -> Either MathError Rational
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"rational-function pole")
    | Bool
otherwise = Rational -> Either MathError Rational
forall a b. b -> Either a b
Right (Polynomial Rational -> Rational -> Rational
forall a. Num a => Polynomial a -> a -> a
evaluatePolynomial Polynomial Rational
p Rational
x Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
den)
  where
    den :: Rational
den = Polynomial Rational -> Rational -> Rational
forall a. Num a => Polynomial a -> a -> a
evaluatePolynomial Polynomial Rational
q Rational
x

-- | Exact rational-function addition with cancellation.
addRationalFunction :: RationalFunction -> RationalFunction -> Either MathError RationalFunction
addRationalFunction :: RationalFunction
-> RationalFunction -> Either MathError RationalFunction
addRationalFunction (RationalFunction Polynomial Rational
a Polynomial Rational
b) (RationalFunction Polynomial Rational
c Polynomial Rational
d) =
    Polynomial Rational
-> Polynomial Rational -> Either MathError RationalFunction
rationalFunction
        (Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
addPolynomial (Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial Polynomial Rational
a Polynomial Rational
d) (Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial Polynomial Rational
c Polynomial Rational
b))
        (Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial Polynomial Rational
b Polynomial Rational
d)

-- | Exact rational-function multiplication with cancellation.
multiplyRationalFunction ::
    RationalFunction -> RationalFunction -> Either MathError RationalFunction
multiplyRationalFunction :: RationalFunction
-> RationalFunction -> Either MathError RationalFunction
multiplyRationalFunction (RationalFunction Polynomial Rational
a Polynomial Rational
b) (RationalFunction Polynomial Rational
c Polynomial Rational
d) =
    Polynomial Rational
-> Polynomial Rational -> Either MathError RationalFunction
rationalFunction (Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial Polynomial Rational
a Polynomial Rational
c) (Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial Polynomial Rational
b Polynomial Rational
d)

{- | Exact partial fractions over Q. Each term is @(numerator,monic irreducible
denominator,power)@ and numerator degree is smaller than denominator degree.
Supports repeated and nonlinear factors. Kronecker factorization and linear
elimination are bounded by the supplied candidate budget.
-}
partialFractions ::
    Int ->
    RationalFunction ->
    Either MathError (Polynomial Rational, [(Polynomial Rational, Polynomial Rational, Int)])
partialFractions :: Int
-> RationalFunction
-> Either
     MathError
     (Polynomial Rational,
      [(Polynomial Rational, Polynomial Rational, Int)])
partialFractions Int
budget (RationalFunction Polynomial Rational
p Polynomial Rational
q) = do
    (whole, remainder) <- Polynomial Rational
-> Polynomial Rational
-> Either MathError (Polynomial Rational, Polynomial Rational)
forall a.
(Eq a, Fractional a) =>
Polynomial a
-> Polynomial a -> Either MathError (Polynomial a, Polynomial a)
dividePolynomial Polynomial Rational
p Polynomial Rational
q
    (_, factors) <- factorRationalPolynomial budget q
    let basis = [(Polynomial Rational
f, Int
k, Int
j, Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Fractional a) =>
Polynomial a -> Polynomial a -> Polynomial a
quotient Polynomial Rational
q (Polynomial Rational -> Int -> Polynomial Rational
forall {a}. (Eq a, Num a) => Polynomial a -> Int -> Polynomial a
power Polynomial Rational
f Int
k)) | (Polynomial Rational
f, Int
m) <- [(Polynomial Rational, Int)]
factors, Int
k <- [Int
1 .. Int
m], Int
j <- [Int
0 .. Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
f Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        columns = [Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
coefficients (Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial ([Rational] -> Polynomial Rational
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial (Int -> Rational -> [Rational]
forall a. Int -> a -> [a]
replicate Int
j Rational
0 [Rational] -> [Rational] -> [Rational]
forall a. [a] -> [a] -> [a]
++ [Rational
1])) Polynomial Rational
b) | (Polynomial Rational
_, Int
_, Int
j, Polynomial Rational
b) <- [(Polynomial Rational, Int, Int, Polynomial Rational)]
basis]
        n = Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
q
        rows =
            [ [[Rational] -> Rational
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
take Int
1 (Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
drop Int
i [Rational]
col)) | [Rational]
col <- [[Rational]]
columns] [Rational] -> [Rational] -> [Rational]
forall a. [a] -> [a] -> [a]
++ [[Rational] -> Rational
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
take Int
1 (Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
drop Int
i (Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
coefficients Polynomial Rational
remainder)))]
            | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            ]
    values <- gaussian rows
    let result = [([Rational] -> Polynomial Rational
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [Polynomial Rational -> Int -> Int -> Rational
valueFor Polynomial Rational
f Int
k Int
j | Int
j <- [Int
0 .. Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
f Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]], Polynomial Rational
f, Int
k) | (Polynomial Rational
f, Int
m) <- [(Polynomial Rational, Int)]
factors, Int
k <- [Int
1 .. Int
m]]
        valueFor Polynomial Rational
f Int
k Int
j = [Rational] -> Rational
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Rational
v | ((Polynomial Rational
g, Int
l, Int
h, Polynomial Rational
_), Rational
v) <- [(Polynomial Rational, Int, Int, Polynomial Rational)]
-> [Rational]
-> [((Polynomial Rational, Int, Int, Polynomial Rational),
     Rational)]
forall a b. [a] -> [b] -> [(a, b)]
zip [(Polynomial Rational, Int, Int, Polynomial Rational)]
basis [Rational]
values, Polynomial Rational
g Polynomial Rational -> Polynomial Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Polynomial Rational
f, Int
l Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k, Int
h Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j]
    pure (whole, filter (\(Polynomial Rational
a, Polynomial Rational
_, Int
_) -> Polynomial Rational -> Int
forall a. Polynomial a -> Int
deg Polynomial Rational
a Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0) result)
  where
    power :: Polynomial a -> Int -> Polynomial a
power Polynomial a
f Int
k = (Polynomial a -> Polynomial a -> Polynomial a)
-> Polynomial a -> [Polynomial a] -> Polynomial a
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl Polynomial a -> Polynomial a -> Polynomial a
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
multiplyPolynomial ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
polynomial [a
1]) (Int -> Polynomial a -> [Polynomial a]
forall a. Int -> a -> [a]
replicate Int
k Polynomial a
f)

gaussian :: [[Rational]] -> Either MathError [Rational]
gaussian :: [[Rational]] -> Either MathError [Rational]
gaussian [] = [Rational] -> Either MathError [Rational]
forall a b. b -> Either a b
Right []
gaussian [[Rational]]
rows = do
    reduced <- ([[Rational]] -> Int -> Either MathError [[Rational]])
-> [[Rational]] -> [Int] -> Either MathError [[Rational]]
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM [[Rational]] -> Int -> Either MathError [[Rational]]
forall {a}.
(Eq a, Fractional a) =>
[[a]] -> Int -> Either MathError [[a]]
eliminate [[Rational]]
rows [Int
0 .. [[Rational]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Rational]]
rows Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
    pure [sum (take 1 (drop (length rows) row)) | row <- reduced]
  where
    eliminate :: [[a]] -> Int -> Either MathError [[a]]
eliminate [[a]]
rs Int
k = case [(Int
i, [a]
row) | (Int
i, [a]
row) <- [Int] -> [[a]] -> [(Int, [a])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[a]]
rs, Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
k, [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
1 (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
k [a]
row)) a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0] of
        [] -> MathError -> Either MathError [[a]]
forall a b. a -> Either a b
Left (Int -> MathError
SingularMatrix Int
k)
        (Int
i, [a]
pivot) : [(Int, [a])]
_ ->
            let chosen :: [[a]]
chosen = [if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k then [a]
pivot else if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i then Int -> [[a]] -> [a]
forall {a}. Int -> [[a]] -> [a]
atRow Int
k [[a]]
rs else [a]
row | (Int
j, [a]
row) <- [Int] -> [[a]] -> [(Int, [a])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[a]]
rs]
                val :: a
val = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
1 (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
k [a]
pivot))
                normalized :: [a]
normalized = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
val) [a]
pivot
             in [[a]] -> Either MathError [[a]]
forall a b. b -> Either a b
Right
                    [ if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k
                        then [a]
normalized
                        else
                            (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith
                                (-)
                                [a]
row
                                ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Num a => a -> a -> a
* ([a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
1 (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
k [a]
row)))) [a]
normalized)
                    | (Int
j, [a]
row) <- [Int] -> [[a]] -> [(Int, [a])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[a]]
chosen
                    ]
    atRow :: Int -> [[a]] -> [a]
atRow Int
k [[a]]
rs = [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat (Int -> [[a]] -> [[a]]
forall a. Int -> [a] -> [a]
take Int
1 (Int -> [[a]] -> [[a]]
forall a. Int -> [a] -> [a]
drop Int
k [[a]]
rs))