-- SPDX-License-Identifier: MPL-2.0

{- | Checked Double root solvers. Absolute residual and mixed absolute/relative
step tolerances are separate. A small step alone never certifies a root.
All function and derivative calls count against the evaluation budget.
Sign-changing isolation finds simple roots; exact polynomial Sturm isolation
additionally finds even/repeated roots with multiplicity. Callback functions
must be deterministic and total on points visited by the chosen method.
-}
module HMath.Root (
    RootOptions (..),
    defaultRootOptions,
    RootResult (..),
    Bracket,
    bracket,
    BracketMethod (..),
    solveBracket,
    bracketScan,
    isolateIntervals,
    DerivativeMethod (..),
    solveDerivative,
    secant,
    steffensen,
    FixedPointMethod (..),
    fixedPoint,
    anderson,
    SystemMethod (..),
    SystemRootResult (..),
    solveSystem,
    polynomialRoots,
    realPolynomialRoots,
) where

import Control.Monad (foldM)
import Data.Complex (Complex)
import Data.List (sortOn)
import HMath.Error
import qualified HMath.Matrix as M
import qualified HMath.Matrix.LU as LU
import qualified HMath.Polynomial as P

{- | Root-solver iteration/evaluation budgets and residual/step tolerances, initialized with
defaultRootOptions.
-}
data RootOptions = RootOptions
    { RootOptions -> Int
rootMaxIterations :: Int
    -- ^ Maximum scalar-root iteration count.
    , RootOptions -> Int
rootMaxEvaluations :: Int
    -- ^ Maximum scalar function/derivative evaluation budget.
    , RootOptions -> Double
rootTolerance :: Double
    -- ^ Requested finite positive residual tolerance.
    , RootOptions -> Double
rootStepTolerance :: Double
    -- ^ Requested finite positive parameter-step tolerance.
    }
    deriving (RootOptions -> RootOptions -> Bool
(RootOptions -> RootOptions -> Bool)
-> (RootOptions -> RootOptions -> Bool) -> Eq RootOptions
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RootOptions -> RootOptions -> Bool
== :: RootOptions -> RootOptions -> Bool
$c/= :: RootOptions -> RootOptions -> Bool
/= :: RootOptions -> RootOptions -> Bool
Eq, Int -> RootOptions -> ShowS
[RootOptions] -> ShowS
RootOptions -> String
(Int -> RootOptions -> ShowS)
-> (RootOptions -> String)
-> ([RootOptions] -> ShowS)
-> Show RootOptions
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> RootOptions -> ShowS
showsPrec :: Int -> RootOptions -> ShowS
$cshow :: RootOptions -> String
show :: RootOptions -> String
$cshowList :: [RootOptions] -> ShowS
showList :: [RootOptions] -> ShowS
Show)

{- | Default root-solving limits: 200 iterations, 4000 evaluations, residual tolerance 1e-10
and step tolerance 1e-12. Solvers validate options and report their own convergence
diagnostics.
-}
defaultRootOptions :: RootOptions
defaultRootOptions :: RootOptions
defaultRootOptions = Int -> Int -> Double -> Double -> RootOptions
RootOptions Int
200 Int
4000 Double
1e-10 Double
1e-12

{- | Scalar root estimate with residual,optional certified bracket error bound,and work
counts.
-}
data RootResult = RootResult
    { RootResult -> Double
rootValue :: Double
    -- ^ Returned scalar root estimate.
    , RootResult -> Double
rootResidual :: Double
    -- ^ Absolute function residual at rootValue.
    , RootResult -> Maybe Double
rootErrorBound :: Maybe Double
    -- ^ Optional bracket-based absolute root error bound; Nothing when the method cannot certify one.
    , RootResult -> Int
rootIterations :: Int
    -- ^ Completed scalar-root iterations.
    , RootResult -> Int
rootEvaluations :: Int
    -- ^ Recorded function/derivative evaluation count.
    }
    deriving (RootResult -> RootResult -> Bool
(RootResult -> RootResult -> Bool)
-> (RootResult -> RootResult -> Bool) -> Eq RootResult
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RootResult -> RootResult -> Bool
== :: RootResult -> RootResult -> Bool
$c/= :: RootResult -> RootResult -> Bool
/= :: RootResult -> RootResult -> Bool
Eq, Int -> RootResult -> ShowS
[RootResult] -> ShowS
RootResult -> String
(Int -> RootResult -> ShowS)
-> (RootResult -> String)
-> ([RootResult] -> ShowS)
-> Show RootResult
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> RootResult -> ShowS
showsPrec :: Int -> RootResult -> ShowS
$cshow :: RootResult -> String
show :: RootResult -> String
$cshowList :: [RootResult] -> ShowS
showList :: [RootResult] -> ShowS
Show)

{- | Validated scalar bracketing state used by checked bracketed root solvers; construct
through the bracketing API.
-}
data Bracket = Bracket Double Double Double Double deriving (Bracket -> Bracket -> Bool
(Bracket -> Bracket -> Bool)
-> (Bracket -> Bracket -> Bool) -> Eq Bracket
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Bracket -> Bracket -> Bool
== :: Bracket -> Bracket -> Bool
$c/= :: Bracket -> Bracket -> Bool
/= :: Bracket -> Bracket -> Bool
Eq, Int -> Bracket -> ShowS
[Bracket] -> ShowS
Bracket -> String
(Int -> Bracket -> ShowS)
-> (Bracket -> String) -> ([Bracket] -> ShowS) -> Show Bracket
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Bracket -> ShowS
showsPrec :: Int -> Bracket -> ShowS
$cshow :: Bracket -> String
show :: Bracket -> String
$cshowList :: [Bracket] -> ShowS
showList :: [Bracket] -> ShowS
Show)

{- | Bracket-preserving root strategy such as bisection,Brent,Ridder,or regula falsi
variants.
-}
data BracketMethod = Bisection | RegulaFalsi | Illinois | Brent | Ridder deriving (BracketMethod -> BracketMethod -> Bool
(BracketMethod -> BracketMethod -> Bool)
-> (BracketMethod -> BracketMethod -> Bool) -> Eq BracketMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: BracketMethod -> BracketMethod -> Bool
== :: BracketMethod -> BracketMethod -> Bool
$c/= :: BracketMethod -> BracketMethod -> Bool
/= :: BracketMethod -> BracketMethod -> Bool
Eq, Int -> BracketMethod -> ShowS
[BracketMethod] -> ShowS
BracketMethod -> String
(Int -> BracketMethod -> ShowS)
-> (BracketMethod -> String)
-> ([BracketMethod] -> ShowS)
-> Show BracketMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> BracketMethod -> ShowS
showsPrec :: Int -> BracketMethod -> ShowS
$cshow :: BracketMethod -> String
show :: BracketMethod -> String
$cshowList :: [BracketMethod] -> ShowS
showList :: [BracketMethod] -> ShowS
Show)

{- | Derivative-based scalar-root strategy, including Newton,Halley,and Householder methods
with their required derivatives.
-}
data DerivativeMethod = Newton | Halley | Householder3 deriving (DerivativeMethod -> DerivativeMethod -> Bool
(DerivativeMethod -> DerivativeMethod -> Bool)
-> (DerivativeMethod -> DerivativeMethod -> Bool)
-> Eq DerivativeMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: DerivativeMethod -> DerivativeMethod -> Bool
== :: DerivativeMethod -> DerivativeMethod -> Bool
$c/= :: DerivativeMethod -> DerivativeMethod -> Bool
/= :: DerivativeMethod -> DerivativeMethod -> Bool
Eq, Int -> DerivativeMethod -> ShowS
[DerivativeMethod] -> ShowS
DerivativeMethod -> String
(Int -> DerivativeMethod -> ShowS)
-> (DerivativeMethod -> String)
-> ([DerivativeMethod] -> ShowS)
-> Show DerivativeMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> DerivativeMethod -> ShowS
showsPrec :: Int -> DerivativeMethod -> ShowS
$cshow :: DerivativeMethod -> String
show :: DerivativeMethod -> String
$cshowList :: [DerivativeMethod] -> ShowS
showList :: [DerivativeMethod] -> ShowS
Show)

-- | Fixed-point acceleration strategy, selecting Picard,Aitken,or Steffensen iteration.
data FixedPointMethod = Picard | Aitken | SteffensenFixed deriving (FixedPointMethod -> FixedPointMethod -> Bool
(FixedPointMethod -> FixedPointMethod -> Bool)
-> (FixedPointMethod -> FixedPointMethod -> Bool)
-> Eq FixedPointMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: FixedPointMethod -> FixedPointMethod -> Bool
== :: FixedPointMethod -> FixedPointMethod -> Bool
$c/= :: FixedPointMethod -> FixedPointMethod -> Bool
/= :: FixedPointMethod -> FixedPointMethod -> Bool
Eq, Int -> FixedPointMethod -> ShowS
[FixedPointMethod] -> ShowS
FixedPointMethod -> String
(Int -> FixedPointMethod -> ShowS)
-> (FixedPointMethod -> String)
-> ([FixedPointMethod] -> ShowS)
-> Show FixedPointMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> FixedPointMethod -> ShowS
showsPrec :: Int -> FixedPointMethod -> ShowS
$cshow :: FixedPointMethod -> String
show :: FixedPointMethod -> String
$cshowList :: [FixedPointMethod] -> ShowS
showList :: [FixedPointMethod] -> ShowS
Show)

-- | Multivariate root strategy, selecting checked Newton or Broyden updates.
data SystemMethod = MultivariateNewton | NewtonTrustRegion Double | Broyden deriving (SystemMethod -> SystemMethod -> Bool
(SystemMethod -> SystemMethod -> Bool)
-> (SystemMethod -> SystemMethod -> Bool) -> Eq SystemMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: SystemMethod -> SystemMethod -> Bool
== :: SystemMethod -> SystemMethod -> Bool
$c/= :: SystemMethod -> SystemMethod -> Bool
/= :: SystemMethod -> SystemMethod -> Bool
Eq, Int -> SystemMethod -> ShowS
[SystemMethod] -> ShowS
SystemMethod -> String
(Int -> SystemMethod -> ShowS)
-> (SystemMethod -> String)
-> ([SystemMethod] -> ShowS)
-> Show SystemMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> SystemMethod -> ShowS
showsPrec :: Int -> SystemMethod -> ShowS
$cshow :: SystemMethod -> String
show :: SystemMethod -> String
$cshowList :: [SystemMethod] -> ShowS
showList :: [SystemMethod] -> ShowS
Show)

-- | Returned multivariate root vector with residual norm and iteration/evaluation counts.
data SystemRootResult = SystemRootResult
    { SystemRootResult -> [Double]
systemRoot :: [Double]
    -- ^ Returned multivariate root parameter vector.
    , SystemRootResult -> Double
systemResidual :: Double
    -- ^ Norm of the residual function at systemRoot.
    , SystemRootResult -> Int
systemIterations :: Int
    -- ^ Completed nonlinear-system iterations.
    , SystemRootResult -> Int
systemEvaluations :: Int
    -- ^ Recorded system function/evaluation count.
    }
    deriving (SystemRootResult -> SystemRootResult -> Bool
(SystemRootResult -> SystemRootResult -> Bool)
-> (SystemRootResult -> SystemRootResult -> Bool)
-> Eq SystemRootResult
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: SystemRootResult -> SystemRootResult -> Bool
== :: SystemRootResult -> SystemRootResult -> Bool
$c/= :: SystemRootResult -> SystemRootResult -> Bool
/= :: SystemRootResult -> SystemRootResult -> Bool
Eq, Int -> SystemRootResult -> ShowS
[SystemRootResult] -> ShowS
SystemRootResult -> String
(Int -> SystemRootResult -> ShowS)
-> (SystemRootResult -> String)
-> ([SystemRootResult] -> ShowS)
-> Show SystemRootResult
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> SystemRootResult -> ShowS
showsPrec :: Int -> SystemRootResult -> ShowS
$cshow :: SystemRootResult -> String
show :: SystemRootResult -> String
$cshowList :: [SystemRootResult] -> ShowS
showList :: [SystemRootResult] -> ShowS
Show)

finite :: Double -> Bool
finite :: Double -> Bool
finite 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)
valid :: RootOptions -> Either MathError ()
valid :: RootOptions -> Either MathError ()
valid RootOptions
o
    | RootOptions -> Int
rootMaxIterations RootOptions
o Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 Bool -> Bool -> Bool
|| RootOptions -> Int
rootMaxEvaluations RootOptions
o Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 =
        MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"root budgets must be positive")
    | Bool -> Bool
not ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [RootOptions -> Double
rootTolerance RootOptions
o, RootOptions -> Double
rootStepTolerance RootOptions
o])
        Bool -> Bool -> Bool
|| RootOptions -> Double
rootTolerance RootOptions
o Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0
        Bool -> Bool -> Bool
|| RootOptions -> Double
rootStepTolerance RootOptions
o Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 =
        MathError -> Either MathError ()
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | Bool
otherwise = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
call :: RootOptions -> (Double -> Double) -> Int -> Double -> Either MathError (Double, Int)
call :: RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
f Int
n Double
x
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= RootOptions -> Int
rootMaxEvaluations RootOptions
o = MathError -> Either MathError (Double, Int)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"root evaluation budget" Int
n)
    | Bool -> Bool
not (Double -> Bool
finite Double
x) = MathError -> Either MathError (Double, Int)
forall a b. a -> Either a b
Left MathError
NonFiniteResult
    | Bool
otherwise = let y :: Double
y = Double -> Double
f Double
x in if Double -> Bool
finite Double
y then (Double, Int) -> Either MathError (Double, Int)
forall a b. b -> Either a b
Right (Double
y, Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) else MathError -> Either MathError (Double, Int)
forall a b. a -> Either a b
Left MathError
NonFiniteResult

-- | Validate finite a<b and endpoint roots or opposite signs; construction never assumes a discontinuity is a root.
bracket :: (Double -> Double) -> Double -> Double -> Either MathError Bracket
bracket :: (Double -> Double) -> Double -> Double -> Either MathError Bracket
bracket Double -> Double
f Double
a Double
b
    | Bool -> Bool
not ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double
a, Double
b]) Bool -> Bool -> Bool
|| Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
b = MathError -> Either MathError Bracket
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"root bracket requires finite a < b")
    | Bool -> Bool
not ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double
fa, Double
fb]) = MathError -> Either MathError Bracket
forall a b. a -> Either a b
Left MathError
NonFiniteResult
    | Double
fa Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0 Bool -> Bool -> Bool
&& Double
fb Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0 Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
signum Double
fa Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double -> Double
forall a. Num a => a -> a
signum Double
fb =
        MathError -> Either MathError Bracket
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"root endpoints do not bracket a sign change")
    | Bool
otherwise = Bracket -> Either MathError Bracket
forall a b. b -> Either a b
Right (Double -> Double -> Double -> Double -> Bracket
Bracket Double
a Double
b Double
fa Double
fb)
  where
    fa :: Double
fa = Double -> Double
f Double
a; fb :: Double
fb = Double -> Double
f Double
b

{- | Brackets are re-evaluated with the supplied function, preventing stale or
mismatched endpoint caches from certifying a different function's root.
-}
solveBracket ::
    BracketMethod -> RootOptions -> (Double -> Double) -> Bracket -> Either MathError RootResult
solveBracket :: BracketMethod
-> RootOptions
-> (Double -> Double)
-> Bracket
-> Either MathError RootResult
solveBracket BracketMethod
method RootOptions
o Double -> Double
f (Bracket Double
a0 Double
b0 Double
_ Double
_) = do
    RootOptions -> Either MathError ()
valid RootOptions
o
    (fa, n1) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
f Int
0 Double
a0
    (fb, n2) <- call o f n1 b0
    _ <- bracket (\Double
x -> if Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
a0 then Double
fa else Double
fb) a0 b0
    loop 0 n2 a0 b0 fa fb fa fb Nothing (0 :: Int)
  where
    loop :: Int
-> Int
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Maybe (Double, Double)
-> t
-> Either MathError RootResult
loop Int
k Int
n Double
a Double
b Double
fa Double
fb Double
wa Double
wb Maybe (Double, Double)
old t
side
        | Double -> Double
forall a. Num a => a -> a
abs Double
fa Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= RootOptions -> Double
rootTolerance RootOptions
o = RootResult -> Either MathError RootResult
forall a b. b -> Either a b
Right (Double -> Double -> Maybe Double -> Int -> Int -> RootResult
RootResult Double
a (Double -> Double
forall a. Num a => a -> a
abs Double
fa) (Double -> Maybe Double
forall a. a -> Maybe a
Just (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a)) Int
k Int
n)
        | Double -> Double
forall a. Num a => a -> a
abs Double
fb Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= RootOptions -> Double
rootTolerance RootOptions
o = RootResult -> Either MathError RootResult
forall a b. b -> Either a b
Right (Double -> Double -> Maybe Double -> Int -> Int -> RootResult
RootResult Double
b (Double -> Double
forall a. Num a => a -> a
abs Double
fb) (Double -> Maybe Double
forall a. a -> Maybe a
Just (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a)) Int
k Int
n)
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= RootOptions -> Int
rootMaxIterations RootOptions
o = MathError -> Either MathError RootResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence (BracketMethod -> String
forall a. Show a => a -> String
show BracketMethod
method) Int
k)
        | Bool
otherwise = do
            let mid :: Double
mid = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                interp :: Double
interp = if Double
wb Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
wa then Double
mid else Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
wa Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
wb Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
wa)
                guarded :: Double -> Double
guarded Double
z = if Double -> Bool
finite Double
z Bool -> Bool -> Bool
&& Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.01 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Bool -> Bool -> Bool
&& Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.01 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) then Double
z else Double
mid
            (x, fx, n') <- case BracketMethod
method of
                BracketMethod
Bisection -> do (y, m) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
f Int
n Double
mid; Right (mid, y, m)
                BracketMethod
RegulaFalsi -> do
                    let z :: Double
z = if Double
interp Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
a Bool -> Bool -> Bool
&& Double
interp Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
b then Double
interp else Double
mid
                    (y, m) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
f Int
n Double
z
                    Right (z, y, m)
                BracketMethod
Illinois -> do
                    let z :: Double
z = if Double
interp Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
a Bool -> Bool -> Bool
&& Double
interp Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
b then Double
interp else Double
mid
                    (y, m) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
f Int
n Double
z
                    Right (z, y, m)
                BracketMethod
Brent -> do
                    -- Inverse quadratic interpolation with the previous trial;
                    -- a central safeguard guarantees interval contraction.
                    let z :: Double
z = case Maybe (Double, Double)
old of
                            Just (Double
c, Double
fc)
                                | Double
fa Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
fb Bool -> Bool -> Bool
&& Double
fa Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
fc Bool -> Bool -> Bool
&& Double
fb Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
fc ->
                                    Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fb Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Double
fa Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fb) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
fa Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fc))
                                        Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fa Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Double
fb Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fa) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
fb Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fc))
                                        Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fa Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fb Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Double
fc Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fa) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
fc Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fb))
                            Maybe (Double, Double)
_ -> Double
interp
                        trial :: Double
trial = Double -> Double
guarded Double
z
                    (y, m) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
f Int
n Double
trial
                    Right (trial, y, m)
                BracketMethod
Ridder -> do
                    (fm, m) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
f Int
n Double
mid
                    let scale = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double -> Double
forall a. Num a => a -> a
abs Double
fa, Double -> Double
forall a. Num a => a -> a
abs Double
fb, Double -> Double
forall a. Num a => a -> a
abs Double
fm]
                        disc = (Double
fm Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
fa Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
fb Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale)
                        trial =
                            if Double
scale Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
|| Double
disc Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0
                                then Double
mid
                                else Double -> Double
guarded (Double
mid Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
mid Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Num a => a -> a
signum (Double
fa Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fb) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
fm Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
disc)
                    if trial == mid then Right (mid, fm, m) else do (y, q) <- call o f m trial; Right (trial, y, q)
            if abs fx <= rootTolerance o
                then Right (RootResult x (abs fx) (Just (b - a)) (k + 1) n')
                else
                    if x == a || x == b || b - a <= rootStepTolerance o * (1 + abs x)
                        then
                            Left
                                (NonConvergence "bracket collapsed without residual convergence (possible discontinuity)" (k + 1))
                        else
                            if signum fa /= signum fx
                                then
                                    loop
                                        (k + 1)
                                        n'
                                        a
                                        x
                                        fa
                                        fx
                                        (if method == Illinois && side == 1 then wa / 2 else fa)
                                        fx
                                        (Just (b, fb))
                                        1
                                else
                                    loop
                                        (k + 1)
                                        n'
                                        x
                                        b
                                        fx
                                        fb
                                        fx
                                        (if method == Illinois && side == (-1) then wb / 2 else fb)
                                        (Just (a, fa))
                                        (-1)

-- | Evaluate a uniform finite grid under the evaluation budget and return every endpoint/sign-change bracket.
bracketScan ::
    RootOptions -> (Double -> Double) -> Double -> Double -> Int -> Either MathError [Bracket]
bracketScan :: RootOptions
-> (Double -> Double)
-> Double
-> Double
-> Int
-> Either MathError [Bracket]
bracketScan RootOptions
o Double -> Double
f Double
a Double
b Int
cells = do
    RootOptions -> Either MathError ()
valid RootOptions
o
    if Int
cells Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1
        Bool -> Bool -> Bool
|| Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
cells Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Integer
forall a. Integral a => a -> Integer
toInteger (RootOptions -> Int
rootMaxEvaluations RootOptions
o)
        Bool -> Bool -> Bool
|| Bool -> Bool
not ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double
a, Double
b])
        Bool -> Bool -> Bool
|| Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
b
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"invalid root scan grid or budget")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    samples <-
        (Int -> Either MathError (Double, Double))
-> [Int] -> Either MathError [(Double, 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
i ->
                let x :: Double
x = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
cells)
                 in do (y, _) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
f Int
i Double
x; Right (x, y)
            )
            [Int
0 .. Int
cells]
    Right
        [ Bracket x y fx fy
        | ((x, fx), (y, fy)) <- zip samples (drop 1 samples)
        , fx == 0 || fy == 0 || signum fx /= signum fy
        ]

-- | Refine sign-change scan brackets and deduplicate endpoint roots. Even roots require polynomial isolation.
isolateIntervals ::
    RootOptions -> (Double -> Double) -> Double -> Double -> Int -> Either MathError [RootResult]
isolateIntervals :: RootOptions
-> (Double -> Double)
-> Double
-> Double
-> Int
-> Either MathError [RootResult]
isolateIntervals RootOptions
o Double -> Double
f Double
a Double
b Int
cells = do
    bs <- RootOptions
-> (Double -> Double)
-> Double
-> Double
-> Int
-> Either MathError [Bracket]
bracketScan RootOptions
o Double -> Double
f Double
a Double
b Int
cells
    rs <- mapM (solveBracket Brent o f) bs
    Right (dedup (sortOn rootValue rs))
  where
    dedup :: [RootResult] -> [RootResult]
dedup [] = []
    dedup (RootResult
x : [RootResult]
xs) =
        RootResult
x
            RootResult -> [RootResult] -> [RootResult]
forall a. a -> [a] -> [a]
: [RootResult] -> [RootResult]
dedup
                ( (RootResult -> Bool) -> [RootResult] -> [RootResult]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile
                    (\RootResult
y -> Double -> Double
forall a. Num a => a -> a
abs (RootResult -> Double
rootValue RootResult
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- RootResult -> Double
rootValue RootResult
x) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= RootOptions -> Double
rootStepTolerance RootOptions
o Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Num a => a -> a
abs (RootResult -> Double
rootValue RootResult
x)))
                    [RootResult]
xs
                )

{- | Derivatives [f', f'', f'''] are required through order 1, 2, or 3.
Householder3 is the order-four rational Householder update.
-}
solveDerivative ::
    DerivativeMethod ->
    RootOptions ->
    (Double -> Double) ->
    [Double -> Double] ->
    Double ->
    Either MathError RootResult
solveDerivative :: DerivativeMethod
-> RootOptions
-> (Double -> Double)
-> [Double -> Double]
-> Double
-> Either MathError RootResult
solveDerivative DerivativeMethod
method RootOptions
o Double -> Double
f [Double -> Double]
derivatives Double
x0 = do
    RootOptions -> Either MathError ()
valid RootOptions
o
    if [Double -> Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double -> Double]
derivatives Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
required
        then MathError -> Either MathError RootResult
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"insufficient analytic derivatives")
        else Int -> Int -> Double -> Either MathError RootResult
loop Int
0 Int
0 Double
x0
  where
    required :: Int
required = case DerivativeMethod
method of DerivativeMethod
Newton -> Int
1; DerivativeMethod
Halley -> Int
2; DerivativeMethod
Householder3 -> Int
3
    loop :: Int -> Int -> Double -> Either MathError RootResult
loop Int
k Int
n Double
x = do
        (y, n1) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
f Int
n Double
x
        if abs y <= rootTolerance o
            then Right (RootResult x (abs y) Nothing k n1)
            else
                if k >= rootMaxIterations o
                    then Left (NonConvergence (show method) k)
                    else do
                        (ds, n2) <-
                            foldM
                                (\([Double]
vs, Int
m) Double -> Double
d -> do (v, q) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
d Int
m Double
x; Right (vs ++ [v], q))
                                ([], n1)
                                (take required derivatives)
                        let d = [Double]
ds [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0
                        if d == 0 then Left (DomainError "zero root derivative") else Right ()
                        let u = Double
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
d
                            delta = case DerivativeMethod
method of
                                DerivativeMethod
Newton -> Double
u
                                DerivativeMethod
Halley -> Double
u Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
ds [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d))
                                DerivativeMethod
Householder3 -> Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
ds [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
ds [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
ds [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
6 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d))
                            x' = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
delta
                        if not (finite x')
                            then Left NonFiniteResult
                            else
                                if x' == x
                                    then Left (NonConvergence "derivative iteration stagnated" k)
                                    else loop (k + 1) n2 x'

-- | Unbracketed two-point secant iteration; zero slopes, stagnation, and exhausted budgets fail explicitly.
secant :: RootOptions -> (Double -> Double) -> Double -> Double -> Either MathError RootResult
secant :: RootOptions
-> (Double -> Double)
-> Double
-> Double
-> Either MathError RootResult
secant RootOptions
o Double -> Double
f Double
a Double
b = do
    RootOptions -> Either MathError ()
valid RootOptions
o
    (fa, n) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
f Int
0 Double
a
    loop 0 n a fa b
  where
    loop :: Int
-> Int -> Double -> Double -> Double -> Either MathError RootResult
loop Int
k Int
n Double
x Double
fx Double
y = do
        (fy, n') <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
f Int
n Double
y
        if abs fy <= rootTolerance o
            then Right (RootResult y (abs fy) Nothing k n')
            else
                if k >= rootMaxIterations o
                    then Left (NonConvergence "secant" k)
                    else
                        if fy == fx
                            then Left (DomainError "secant zero slope")
                            else
                                let z = Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fy Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
fy Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fx)
                                 in if z == y then Left (NonConvergence "secant stagnation" k) else loop (k + 1) n' y fy z

-- | Derivative-free Steffensen root update using f(x+f(x)); termination checks the original residual.
steffensen :: RootOptions -> (Double -> Double) -> Double -> Either MathError RootResult
steffensen :: RootOptions
-> (Double -> Double) -> Double -> Either MathError RootResult
steffensen RootOptions
o Double -> Double
f Double
x0 = do RootOptions -> Either MathError ()
valid RootOptions
o; Int -> Int -> Double -> Either MathError RootResult
loop Int
0 Int
0 Double
x0
  where
    loop :: Int -> Int -> Double -> Either MathError RootResult
loop Int
k Int
n Double
x = do
        (fx, m) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
f Int
n Double
x
        if abs fx <= rootTolerance o
            then Right (RootResult x (abs fx) Nothing k m)
            else
                if k >= rootMaxIterations o
                    then Left (NonConvergence "Steffensen" k)
                    else do
                        (fy, q) <- call o f m (x + fx)
                        let den = Double
fy Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fx
                        if den == 0
                            then Left (DomainError "Steffensen zero denominator")
                            else loop (k + 1) q (x - fx * fx / den)

-- | Picard iteration, non-restarted Aitken sequence acceleration, or restarted Steffensen acceleration.
fixedPoint ::
    FixedPointMethod -> RootOptions -> (Double -> Double) -> Double -> Either MathError RootResult
fixedPoint :: FixedPointMethod
-> RootOptions
-> (Double -> Double)
-> Double
-> Either MathError RootResult
fixedPoint FixedPointMethod
method RootOptions
o Double -> Double
g Double
x0 = do
    RootOptions -> Either MathError ()
valid RootOptions
o
    if FixedPointMethod
method FixedPointMethod -> FixedPointMethod -> Bool
forall a. Eq a => a -> a -> Bool
== FixedPointMethod
Aitken then Int -> Int -> Double -> Either MathError RootResult
aitkenLoop Int
0 Int
0 Double
x0 else Int -> Int -> Double -> Either MathError RootResult
loop Int
0 Int
0 Double
x0
  where
    loop :: Int -> Int -> Double -> Either MathError RootResult
loop Int
k Int
n Double
x = do
        (y, m) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
g Int
n Double
x
        let r = Double -> Double
forall a. Num a => a -> a
abs (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x)
        if r <= rootTolerance o
            then Right (RootResult x r Nothing k m)
            else
                if k >= rootMaxIterations o
                    then Left (NonConvergence (show method) k)
                    else case method of
                        FixedPointMethod
Picard -> Int -> Int -> Double -> Either MathError RootResult
loop (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int
m Double
y
                        FixedPointMethod
Aitken -> Int -> Int -> Double -> Double -> Either MathError RootResult
accelerate Int
k Int
m Double
x Double
y
                        FixedPointMethod
SteffensenFixed -> Int -> Int -> Double -> Double -> Either MathError RootResult
accelerate Int
k Int
m Double
x Double
y
    -- Aitken accelerates the raw sequence without restarting; Steffensen
    -- restarts at each Δ² estimate. Each tests the original fixed-point residual.
    aitkenLoop :: Int -> Int -> Double -> Either MathError RootResult
aitkenLoop Int
k Int
n Double
x = do
        (y, m) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
g Int
n Double
x
        (z, q) <- call o g m y
        let denominator = Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x
            estimate = if Double
denominator Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
y else Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
denominator
        (value, count) <- call o g q estimate
        let residual = Double -> Double
forall a. Num a => a -> a
abs (Double
value Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
estimate)
        if residual <= rootTolerance o
            then Right (RootResult estimate residual Nothing k count)
            else
                if k >= rootMaxIterations o
                    then Left (NonConvergence "Aitken raw sequence" k)
                    else aitkenLoop (k + 1) count y
    accelerate :: Int -> Int -> Double -> Double -> Either MathError RootResult
accelerate Int
k Int
n Double
x Double
y = do
        (z, m) <- RootOptions
-> (Double -> Double)
-> Int
-> Double
-> Either MathError (Double, Int)
call RootOptions
o Double -> Double
g Int
n Double
y
        let d = Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x; next = if Double
d Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
y else Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
d
        loop (k + 1) m next

norm :: [Double] -> Double
norm :: [Double] -> Double
norm [Double]
xs = [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]
: (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
xs)
vectorCall ::
    RootOptions -> ([Double] -> [Double]) -> Int -> [Double] -> Either MathError ([Double], Int)
vectorCall :: RootOptions
-> ([Double] -> [Double])
-> Int
-> [Double]
-> Either MathError ([Double], Int)
vectorCall RootOptions
o [Double] -> [Double]
f Int
n [Double]
x
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= RootOptions -> Int
rootMaxEvaluations RootOptions
o = MathError -> Either MathError ([Double], Int)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"system evaluation budget" Int
n)
    | Bool -> Bool
not ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
x) = MathError -> Either MathError ([Double], Int)
forall a b. a -> Either a b
Left MathError
NonFiniteResult
    | Bool
otherwise =
        let y :: [Double]
y = [Double] -> [Double]
f [Double]
x
         in if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x
                then MathError -> Either MathError ([Double], Int)
forall a b. a -> Either a b
Left (Int -> Int -> MathError
VectorLengthMismatch ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x) ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y))
                else if (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
y then ([Double], Int) -> Either MathError ([Double], Int)
forall a b. b -> Either a b
Right ([Double]
y, Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) else MathError -> Either MathError ([Double], Int)
forall a b. a -> Either a b
Left MathError
NonFiniteResult
linearSolve :: [[Double]] -> [Double] -> Either MathError [Double]
linearSolve :: [[Double]] -> [Double] -> Either MathError [Double]
linearSolve [[Double]]
a [Double]
b = do m <- Int -> Int -> [Double] -> Either MathError (Matrix Double)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
a) ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
a) ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
a); LU.solve m b

{- | Anderson mixing with depth m and Tikhonov regularization. Residuals are
g(x)-x. Solving the constrained residual Gram system is a bounded reference
implementation; it does not pretend to use a low-storage QR update.
-}
anderson ::
    Int ->
    Double ->
    RootOptions ->
    ([Double] -> [Double]) ->
    [Double] ->
    Either MathError SystemRootResult
anderson :: Int
-> Double
-> RootOptions
-> ([Double] -> [Double])
-> [Double]
-> Either MathError SystemRootResult
anderson Int
depth Double
regularization RootOptions
o [Double] -> [Double]
g [Double]
x0 = do
    RootOptions -> Either MathError ()
valid RootOptions
o
    if Int
depth Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 Bool -> Bool -> Bool
|| Double
regularization Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| Bool -> Bool
not (Double -> Bool
finite Double
regularization) Bool -> Bool -> Bool
|| [Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
x0
        then MathError -> Either MathError SystemRootResult
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"Anderson depth, regularization, or dimension")
        else Int
-> Int
-> [([Double], [Double])]
-> [Double]
-> Either MathError SystemRootResult
loop Int
0 Int
0 [] [Double]
x0
  where
    loop :: Int
-> Int
-> [([Double], [Double])]
-> [Double]
-> Either MathError SystemRootResult
loop Int
k Int
n [([Double], [Double])]
hist [Double]
x = do
        (gx, m) <- RootOptions
-> ([Double] -> [Double])
-> Int
-> [Double]
-> Either MathError ([Double], Int)
vectorCall RootOptions
o [Double] -> [Double]
g Int
n [Double]
x
        let r = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
gx [Double]
x
        if norm r <= rootTolerance o
            then Right (SystemRootResult x (norm r) k m)
            else
                if k >= rootMaxIterations o
                    then Left (NonConvergence "Anderson" k)
                    else do
                        let hs = Int -> [([Double], [Double])] -> [([Double], [Double])]
forall a. Int -> [a] -> [a]
take Int
depth (([Double]
gx, [Double]
r) ([Double], [Double])
-> [([Double], [Double])] -> [([Double], [Double])]
forall a. a -> [a] -> [a]
: [([Double], [Double])]
hist)
                            p = [([Double], [Double])] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [([Double], [Double])]
hs
                            dot [a]
a [a]
b = [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
(*) [a]
a [a]
b)
                            gram =
                                [ [[Double] -> [Double] -> Double
forall {a}. Num a => [a] -> [a] -> a
dot [Double]
a [Double]
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Double
regularization else Double
0 | (Int
j, ([Double]
_, [Double]
b)) <- [Int] -> [([Double], [Double])] -> [(Int, ([Double], [Double]))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [([Double], [Double])]
hs]
                                | (Int
i, ([Double]
_, [Double]
a)) <- [Int] -> [([Double], [Double])] -> [(Int, ([Double], [Double]))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [([Double], [Double])]
hs
                                ]
                            aug = ([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
1]) [[Double]]
gram [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ [Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
p Double
1 [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
0]]
                        coeffs <- linearSolve aug (replicate p 0 ++ [1])
                        let next = [[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
a [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) | (Double
c, ([Double]
a, [Double]
_)) <- [Double]
-> [([Double], [Double])] -> [(Double, ([Double], [Double]))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
coeffs [([Double], [Double])]
hs] | Int
j <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                        loop (k + 1) m hs next

-- | Solve a square analytic-Jacobian system using Newton, trust-radius Newton, or rank-one Broyden updates.
solveSystem ::
    SystemMethod ->
    RootOptions ->
    ([Double] -> [Double]) ->
    ([Double] -> [[Double]]) ->
    [Double] ->
    Either MathError SystemRootResult
solveSystem :: SystemMethod
-> RootOptions
-> ([Double] -> [Double])
-> ([Double] -> [[Double]])
-> [Double]
-> Either MathError SystemRootResult
solveSystem SystemMethod
method RootOptions
o [Double] -> [Double]
f [Double] -> [[Double]]
jac [Double]
x0 = do
    RootOptions -> Either MathError ()
valid RootOptions
o
    if [Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
x0 then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"empty nonlinear system") else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    case SystemMethod
method of
        NewtonTrustRegion Double
radius | Double
radius Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| Bool -> Bool
not (Double -> Bool
finite Double
radius) -> MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"positive trust radius required")
        SystemMethod
_ -> () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    let radius :: Double
radius = case SystemMethod
method of NewtonTrustRegion Double
r -> Double
r; SystemMethod
_ -> Double
1
    (fx, n) <- RootOptions
-> ([Double] -> [Double])
-> Int
-> [Double]
-> Either MathError ([Double], Int)
vectorCall RootOptions
o [Double] -> [Double]
f Int
0 [Double]
x0
    if norm fx <= rootTolerance o
        then Right (SystemRootResult x0 (norm fx) 0 n)
        else do
            if n >= rootMaxEvaluations o then Left (NonConvergence "Jacobian evaluation budget" n) else Right ()
            j <- checkedJac x0
            loop 0 (n + 1) radius x0 fx j
  where
    checkedJac :: [Double] -> Either MathError [[Double]]
checkedJac [Double]
x =
        let j :: [[Double]]
j = [Double] -> [[Double]]
jac [Double]
x; n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x
         in if [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n Bool -> Bool -> Bool
|| ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
j
                then MathError -> Either MathError [[Double]]
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"Jacobian shape")
                else if (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
j) then [[Double]] -> Either MathError [[Double]]
forall a b. b -> Either a b
Right [[Double]]
j else MathError -> Either MathError [[Double]]
forall a b. a -> Either a b
Left MathError
NonFiniteResult
    loop :: Int
-> Int
-> Double
-> [Double]
-> [Double]
-> [[Double]]
-> Either MathError SystemRootResult
loop Int
k Int
n Double
radius [Double]
x [Double]
fx [[Double]]
j
        | [Double] -> Double
norm [Double]
fx Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= RootOptions -> Double
rootTolerance RootOptions
o = SystemRootResult -> Either MathError SystemRootResult
forall a b. b -> Either a b
Right ([Double] -> Double -> Int -> Int -> SystemRootResult
SystemRootResult [Double]
x ([Double] -> Double
norm [Double]
fx) Int
k Int
n)
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= RootOptions -> Int
rootMaxIterations RootOptions
o Bool -> Bool -> Bool
|| Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= RootOptions -> Int
rootMaxEvaluations RootOptions
o = MathError -> Either MathError SystemRootResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"nonlinear system" Int
k)
        | Bool
otherwise = do
            raw <- [[Double]] -> [Double] -> Either MathError [Double]
linearSolve [[Double]]
j ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
negate [Double]
fx)
            let magnitude = [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]
: (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
raw)
                unitNorm = if Double
magnitude Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [(Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
magnitude) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) | Double
z <- [Double]
raw])
                step = case SystemMethod
method of
                    NewtonTrustRegion Double
_
                        | Double
magnitude Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
magnitude Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
radius Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
unitNorm ->
                            (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Double
z -> (Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
magnitude) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
radius Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
unitNorm)) [Double]
raw
                    SystemMethod
_ -> [Double]
raw
                candidate = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) [Double]
x [Double]
step
            (fy, n1) <- vectorCall o f n candidate
            case method of
                SystemMethod
_ | [Double] -> Double
norm [Double]
fy Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= RootOptions -> Double
rootTolerance RootOptions
o -> SystemRootResult -> Either MathError SystemRootResult
forall a b. b -> Either a b
Right ([Double] -> Double -> Int -> Int -> SystemRootResult
SystemRootResult [Double]
candidate ([Double] -> Double
norm [Double]
fy) (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int
n1)
                NewtonTrustRegion Double
_ | [Double] -> Double
norm [Double]
fy Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= [Double] -> Double
norm [Double]
fx -> Int
-> Int
-> Double
-> [Double]
-> [Double]
-> [[Double]]
-> Either MathError SystemRootResult
loop (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int
n1 (Double
radius Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) [Double]
x [Double]
fx [[Double]]
j
                SystemMethod
_ -> do
                    (j', n2) <- case SystemMethod
method of
                        SystemMethod
Broyden -> do
                            let y :: [Double]
y = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
fy [Double]
fx
                                js :: [Double]
js = ([Double] -> Double) -> [[Double]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([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 b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
step) [[Double]]
j
                                stepMagnitude :: Double
stepMagnitude = [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]
: (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
step)
                                unit :: [Double]
unit = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
stepMagnitude) [Double]
step
                                denominator :: Double
denominator = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Double
z -> Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z) [Double]
unit)
                            if Double
stepMagnitude Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
                                then MathError -> Either MathError ([[Double]], Int)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"Broyden zero step" Int
k)
                                else do
                                    let updated :: [[Double]]
updated =
                                            ([Double] -> Double -> [Double])
-> [[Double]] -> [Double] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith
                                                (\[Double]
row Double
delta -> (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) [Double]
row ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((Double
delta Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
stepMagnitude) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
denominator)) [Double]
unit))
                                                [[Double]]
j
                                                ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
y [Double]
js)
                                    if (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
updated) then ([[Double]], Int) -> Either MathError ([[Double]], Int)
forall a b. b -> Either a b
Right ([[Double]]
updated, Int
n1) else MathError -> Either MathError ([[Double]], Int)
forall a b. a -> Either a b
Left MathError
NonFiniteResult
                        SystemMethod
_ ->
                            if [Double] -> Double
norm [Double]
fy Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= RootOptions -> Double
rootTolerance RootOptions
o
                                then ([[Double]], Int) -> Either MathError ([[Double]], Int)
forall a b. b -> Either a b
Right ([[Double]]
j, Int
n1)
                                else do
                                    if Int
n1 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= RootOptions -> Int
rootMaxEvaluations RootOptions
o
                                        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"Jacobian evaluation budget" Int
n1)
                                        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
                                    a <- [Double] -> Either MathError [[Double]]
checkedJac [Double]
candidate
                                    Right (a, n1 + 1)
                    loop (k + 1) n2 (min (radius * 1.5) 1e6) candidate fy j'

-- | Reuse the checked Aberth-Ehrlich simultaneous complex solver for an ascending-coefficient polynomial.
polynomialRoots :: RootOptions -> P.Polynomial (Complex Double) -> Either MathError [Complex Double]
polynomialRoots :: RootOptions
-> Polynomial (Complex Double) -> Either MathError [Complex Double]
polynomialRoots RootOptions
o Polynomial (Complex Double)
p = do RootOptions -> Either MathError ()
valid RootOptions
o; Int
-> Double
-> Polynomial (Complex Double)
-> Either MathError [Complex Double]
P.aberthEhrlich (RootOptions -> Int
rootMaxIterations RootOptions
o) (RootOptions -> Double
rootTolerance RootOptions
o) Polynomial (Complex Double)
p

-- | Use exact square-free factorization and Sturm isolation; return rational intervals and integer multiplicities.
realPolynomialRoots ::
    Int -> Rational -> P.Polynomial Rational -> Either MathError [(Rational, Rational, Int)]
realPolynomialRoots :: Int
-> Rational
-> Polynomial Rational
-> Either MathError [(Rational, Rational, Int)]
realPolynomialRoots Int
budget Rational
tolerance Polynomial Rational
p = do
    if [Rational] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null (Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
P.coefficients Polynomial Rational
p)
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"zero polynomial has infinitely many roots")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    pieces <-
        ((Polynomial Rational, Int)
 -> Either MathError [(Rational, Rational, Int)])
-> [(Polynomial Rational, Int)]
-> Either MathError [[(Rational, Rational, Int)]]
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
            ( \(Polynomial Rational
factor, Int
multiplicity) -> do
                intervals <- Int
-> Rational
-> Polynomial Rational
-> Either MathError [(Rational, Rational)]
P.isolateRealRoots Int
budget Rational
tolerance Polynomial Rational
factor
                Right [(a, b, multiplicity) | (a, b) <- intervals]
            )
            (Polynomial Rational -> [(Polynomial Rational, Int)]
P.squareFreeFactors Polynomial Rational
p)
    Right (sortOn (\(Rational
a, Rational
_, Int
_) -> Rational
a) (concat pieces))