-- SPDX-License-Identifier: MPL-2.0

{- | Bounded symbolic real algebra over the common 'Expr' language.

This is a reference CAS with explicit unsupported results, not a decision
procedure for arbitrary elementary expressions. Variables are real; assumptions
are caller assertions, checked for contradictions and retained as proof premises.
Cancellation and log/power rewrites require domain proofs. Generic expansion
and ring collection preserve partial-function domains by refusing to erase an
unproved undefined subexpression. All public transformations validate input and
output against caller limits. Polynomial operations additionally bound degrees
and rational coefficient growth before entering existing exact algorithms.

Integration handles rational polynomials, affine elementary primitives,
polynomial-times-exp/sin/cos by parts, and rational partial fractions with linear
or simple quadratic denominators. Limits and Laurent jets are exact for rational
functions; regular elementary Taylor jets use symbolic derivatives. General
transcendental limits, non-elementary integrals, higher-dimensional eigensystems,
and branch-dependent identities return 'Unsupported'. There are no hidden
numerical guesses or silently discarded branch conditions.
-}
module HMath.Symbolic (
    module HMath.Expression,
    Property (..),
    Assumptions,
    assumptions,
    noAssumptions,
    proves,
    SymbolicOptions (..),
    defaultSymbolicOptions,
    normalize,
    simplify,
    expand,
    collectTerms,
    factorPolynomial,
    PolynomialFactorization (..),
    rationalSimplify,
    trigonometricSimplify,
    logarithmicSimplify,
    substitute,
    Pattern (..),
    matchPattern,
    RewriteRule (..),
    applyRewrite,
    differentiate,
    gradient,
    jacobian,
    hessian,
    integrateBasic,
    integrateByParts,
    integrateSubstitution,
    integrateRational,
    Direction (..),
    LimitResult (..),
    symbolicLimit,
    Series (..),
    taylorSeries,
    maclaurinSeries,
    laurentSeries,
    formalPowerSeries,
    multiplyFormalSeries,
    reciprocalFormalSeries,
    toPolynomial,
    fromPolynomial,
    exactPolynomialRoots,
    numericalPolynomialRoots,
    RelationOp (..),
    Relation (..),
    RelationResult (..),
    simplifyRelation,
    solveLinearInequality,
    LinearSolution (..),
    solveLinearEquations,
    SymbolicMatrix,
    symbolicMatrix,
    matrixExpressions,
    matrixAdd,
    matrixMultiply,
    symbolicDeterminant,
    Guarded (..),
    symbolicEigenvalues,
    numericalSubstitution,
    toLaTeX,
    toMathML,
) where

import Control.Monad (foldM)
import Data.Complex (Complex ((:+)))
import qualified Data.Foldable as Foldable (foldl')
import Data.List (intercalate, nub, sort)
import qualified Data.Map.Strict as Map
import Data.Ratio (denominator, numerator)
import qualified Data.Set as Set
import HMath.Expression
import qualified HMath.Polynomial as P

failure :: ErrorKind -> Expr -> String -> Either ExprError a
failure :: forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
k Expr
e = ExprError -> Either ExprError a
forall a b. a -> Either a b
Left (ExprError -> Either ExprError a)
-> (String -> ExprError) -> String -> Either ExprError a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ErrorKind -> SourcePos -> String -> ExprError
ExprError ErrorKind
k (Expr -> SourcePos
expressionPosition Expr
e)
unsupported :: Expr -> String -> Either ExprError a
unsupported :: forall a. Expr -> String -> Either ExprError a
unsupported = ErrorKind -> Expr -> String -> Either ExprError a
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
Unsupported

{- | Symbolic domain fact such as real,positive,nonzero,or sign restrictions, used to guard
domain-sensitive identities.
-}
data Property = Real | NonZero | Positive | NonNegative | Integral | Negative | NonPositive
    deriving (Property -> Property -> Bool
(Property -> Property -> Bool)
-> (Property -> Property -> Bool) -> Eq Property
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Property -> Property -> Bool
== :: Property -> Property -> Bool
$c/= :: Property -> Property -> Bool
/= :: Property -> Property -> Bool
Eq, Eq Property
Eq Property =>
(Property -> Property -> Ordering)
-> (Property -> Property -> Bool)
-> (Property -> Property -> Bool)
-> (Property -> Property -> Bool)
-> (Property -> Property -> Bool)
-> (Property -> Property -> Property)
-> (Property -> Property -> Property)
-> Ord Property
Property -> Property -> Bool
Property -> Property -> Ordering
Property -> Property -> Property
forall a.
Eq a =>
(a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: Property -> Property -> Ordering
compare :: Property -> Property -> Ordering
$c< :: Property -> Property -> Bool
< :: Property -> Property -> Bool
$c<= :: Property -> Property -> Bool
<= :: Property -> Property -> Bool
$c> :: Property -> Property -> Bool
> :: Property -> Property -> Bool
$c>= :: Property -> Property -> Bool
>= :: Property -> Property -> Bool
$cmax :: Property -> Property -> Property
max :: Property -> Property -> Property
$cmin :: Property -> Property -> Property
min :: Property -> Property -> Property
Ord, Int -> Property -> ShowS
[Property] -> ShowS
Property -> String
(Int -> Property -> ShowS)
-> (Property -> String) -> ([Property] -> ShowS) -> Show Property
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Property -> ShowS
showsPrec :: Int -> Property -> ShowS
$cshow :: Property -> String
show :: Property -> String
$cshowList :: [Property] -> ShowS
showList :: [Property] -> ShowS
Show)

{- | Validated variable-property environment. Build through the assumptions API;
noAssumptions proves no additional variable facts.
-}
newtype Assumptions = Assumptions (Map.Map String [Property]) deriving (Assumptions -> Assumptions -> Bool
(Assumptions -> Assumptions -> Bool)
-> (Assumptions -> Assumptions -> Bool) -> Eq Assumptions
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Assumptions -> Assumptions -> Bool
== :: Assumptions -> Assumptions -> Bool
$c/= :: Assumptions -> Assumptions -> Bool
/= :: Assumptions -> Assumptions -> Bool
Eq, Int -> Assumptions -> ShowS
[Assumptions] -> ShowS
Assumptions -> String
(Int -> Assumptions -> ShowS)
-> (Assumptions -> String)
-> ([Assumptions] -> ShowS)
-> Show Assumptions
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Assumptions -> ShowS
showsPrec :: Int -> Assumptions -> ShowS
$cshow :: Assumptions -> String
show :: Assumptions -> String
$cshowList :: [Assumptions] -> ShowS
showList :: [Assumptions] -> ShowS
Show)

{- | Construct an empty variable-property environment. Simplifications requiring
positivity/nonzero facts remain guarded or unchanged without proof.
-}
noAssumptions :: Assumptions
noAssumptions :: Assumptions
noAssumptions = Map String [Property] -> Assumptions
Assumptions Map String [Property]
forall k a. Map k a
Map.empty

{- | Only variable assertions are accepted. Positive/negative/zero constraints
are checked for contradictions, including repeated entries for a variable.
-}
assumptions :: [(String, [Property])] -> Either ExprError Assumptions
assumptions :: [(String, [Property])] -> Either ExprError Assumptions
assumptions [(String, [Property])]
entries = do
    let merged :: Map String [Property]
merged = ([Property] -> [Property] -> [Property])
-> [(String, [Property])] -> Map String [Property]
forall k a. Ord k => (a -> a -> a) -> [(k, a)] -> Map k a
Map.fromListWith [Property] -> [Property] -> [Property]
forall a. [a] -> [a] -> [a]
(++) [(String, [Property])]
entries
    ((String, [Property]) -> Either ExprError ())
-> [(String, [Property])] -> Either ExprError ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (String, [Property]) -> Either ExprError ()
forall {t :: * -> *}.
Foldable t =>
(String, t Property) -> Either ExprError ()
check (Map String [Property] -> [(String, [Property])]
forall k a. Map k a -> [(k, a)]
Map.toList Map String [Property]
merged)
    Assumptions -> Either ExprError Assumptions
forall a b. b -> Either a b
Right (Map String [Property] -> Assumptions
Assumptions (([Property] -> [Property])
-> Map String [Property] -> Map String [Property]
forall a b k. (a -> b) -> Map k a -> Map k b
Map.map ([Property] -> [Property]
forall a. Eq a => [a] -> [a]
nub ([Property] -> [Property])
-> ([Property] -> [Property]) -> [Property] -> [Property]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Property] -> [Property]
forall a. Ord a => [a] -> [a]
sort) Map String [Property]
merged))
  where
    check :: (String, t Property) -> Either ExprError ()
check (String
name, t Property
ps) = do
        Limits -> Expr -> Either ExprError ()
validateExpression Limits
defaultLimits (String -> Expr
variable String
name)
        if (Property
Positive Property -> t Property -> Bool
forall a. Eq a => a -> t a -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` t Property
ps Bool -> Bool -> Bool
&& (Property
Negative Property -> t Property -> Bool
forall a. Eq a => a -> t a -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` t Property
ps Bool -> Bool -> Bool
|| Property
NonPositive Property -> t Property -> Bool
forall a. Eq a => a -> t a -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` t Property
ps))
            Bool -> Bool -> Bool
|| (Property
Negative Property -> t Property -> Bool
forall a. Eq a => a -> t a -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` t Property
ps Bool -> Bool -> Bool
&& Property
NonNegative Property -> t Property -> Bool
forall a. Eq a => a -> t a -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` t Property
ps)
            Bool -> Bool -> Bool
|| (Property
NonZero Property -> t Property -> Bool
forall a. Eq a => a -> t a -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` t Property
ps Bool -> Bool -> Bool
&& Property
NonNegative Property -> t Property -> Bool
forall a. Eq a => a -> t a -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` t Property
ps Bool -> Bool -> Bool
&& Property
NonPositive Property -> t Property -> Bool
forall a. Eq a => a -> t a -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` t Property
ps)
            then ErrorKind -> Expr -> String -> Either ExprError ()
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
DomainFailure (String -> Expr
variable String
name) String
"contradictory assumptions"
            else () -> Either ExprError ()
forall a b. b -> Either a b
Right ()

{- | Expression budgets together with permitted variable assumptions for checked symbolic
operations.
-}
data SymbolicOptions = SymbolicOptions
    { SymbolicOptions -> Limits
symbolicLimits :: Limits
    -- ^ Expression resource/numeric budgets used by symbolic operations.
    , SymbolicOptions -> Assumptions
symbolicAssumptions :: Assumptions
    -- ^ Variable facts permitted when proving domain-sensitive rewrites.
    }
    deriving (SymbolicOptions -> SymbolicOptions -> Bool
(SymbolicOptions -> SymbolicOptions -> Bool)
-> (SymbolicOptions -> SymbolicOptions -> Bool)
-> Eq SymbolicOptions
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: SymbolicOptions -> SymbolicOptions -> Bool
== :: SymbolicOptions -> SymbolicOptions -> Bool
$c/= :: SymbolicOptions -> SymbolicOptions -> Bool
/= :: SymbolicOptions -> SymbolicOptions -> Bool
Eq, Int -> SymbolicOptions -> ShowS
[SymbolicOptions] -> ShowS
SymbolicOptions -> String
(Int -> SymbolicOptions -> ShowS)
-> (SymbolicOptions -> String)
-> ([SymbolicOptions] -> ShowS)
-> Show SymbolicOptions
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> SymbolicOptions -> ShowS
showsPrec :: Int -> SymbolicOptions -> ShowS
$cshow :: SymbolicOptions -> String
show :: SymbolicOptions -> String
$cshowList :: [SymbolicOptions] -> ShowS
showList :: [SymbolicOptions] -> ShowS
Show)

{- | Default expression budgets with no assumptions; domain-sensitive rules
require their preconditions to be proved before rewriting.
-}
defaultSymbolicOptions :: SymbolicOptions
defaultSymbolicOptions :: SymbolicOptions
defaultSymbolicOptions = Limits -> Assumptions -> SymbolicOptions
SymbolicOptions Limits
defaultLimits Assumptions
noAssumptions

{- | Given assumptions,requested Property,then Expr, conservatively test whether the fact is
established. False means not proved, not a proof of its negation.
-}
proves :: Assumptions -> Property -> Expr -> Bool
proves :: Assumptions -> Property -> Expr -> Bool
proves a :: Assumptions
a@(Assumptions Map String [Property]
names) Property
prop Expr
e = case Expr -> ExprNode
expressionNode Expr
e of
    Literal Rational
q -> case Property
prop of
        Property
Real -> Bool
True
        Property
NonZero -> Rational
q Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0
        Property
Positive -> Rational
q Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
0
        Property
NonNegative -> Rational
q Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
0
        Property
Integral -> Rational -> Integer
forall a. Ratio a -> a
denominator Rational
q Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1
        Property
Negative -> Rational
q Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Rational
0
        Property
NonPositive -> Rational
q Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
0
    Variable String
n -> Property
prop Property -> Property -> Bool
forall a. Eq a => a -> a -> Bool
== Property
Real Bool -> Bool -> Bool
|| (Property -> Bool) -> [Property] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Property -> Property -> Bool
implies Property
prop) ([Property] -> String -> Map String [Property] -> [Property]
forall k a. Ord k => a -> k -> Map k a -> a
Map.findWithDefault [] String
n Map String [Property]
names)
    Constant NamedConstant
_ -> Property
prop Property -> [Property] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Property
Real, Property
NonZero, Property
Positive, Property
NonNegative]
    Unary UnaryOp
Negate Expr
x -> case Property
prop of
        Property
Positive -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Negative Expr
x
        Property
Negative -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x
        Property
NonNegative -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonPositive Expr
x
        Property
NonPositive -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonNegative Expr
x
        Property
_ -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
prop Expr
x
    Binary BinaryOp
Multiply Expr
x Expr
y -> case Property
prop of
        Property
NonZero -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
x Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
y
        Property
Positive ->
            (Expr
x Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Expr
y Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
x)
                Bool -> Bool -> Bool
|| (Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
y)
                Bool -> Bool -> Bool
|| (Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Negative Expr
x Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Negative Expr
y)
        Property
NonNegative ->
            (Expr
x Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Expr
y Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
x)
                Bool -> Bool -> Bool
|| (Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonNegative Expr
x Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonNegative Expr
y)
                Bool -> Bool -> Bool
|| (Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonPositive Expr
x Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonPositive Expr
y)
        Property
Integral -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Integral Expr
x Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Integral Expr
y
        Property
Real -> Assumptions -> Expr -> Bool
defined Assumptions
a Expr
e
        Property
_ -> Bool
False
    Binary BinaryOp
Add Expr
x Expr
y -> case Property
prop of
        Property
NonZero -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
e Bool -> Bool -> Bool
|| Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Negative Expr
e
        Property
Positive ->
            (Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonNegative Expr
y) Bool -> Bool -> Bool
|| (Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonNegative Expr
x Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
y)
        Property
NonNegative -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonNegative Expr
x Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonNegative Expr
y
        Property
Integral -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Integral Expr
x Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Integral Expr
y
        Property
Real -> Assumptions -> Expr -> Bool
defined Assumptions
a Expr
e
        Property
_ -> Bool
False
    Binary BinaryOp
Divide Expr
x Expr
y ->
        Property
prop Property -> Property -> Bool
forall a. Eq a => a -> a -> Bool
== Property
NonZero Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
e Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
x Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
y
            Bool -> Bool -> Bool
|| Property
prop Property -> Property -> Bool
forall a. Eq a => a -> a -> Bool
== Property
Real Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
e
    Binary BinaryOp
Power Expr
x (Expr SourcePos
_ (Literal Rational
q)) ->
        (Property
prop Property -> Property -> Bool
forall a. Eq a => a -> a -> Bool
== Property
Real Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
e)
            Bool -> Bool -> Bool
|| (Property
prop Property -> Property -> Bool
forall a. Eq a => a -> a -> Bool
== Property
NonZero Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
x Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
e)
            Bool -> Bool -> Bool
|| ( Property
prop Property -> Property -> Bool
forall a. Eq a => a -> a -> Bool
== Property
Positive
                    Bool -> Bool -> Bool
&& ((Rational -> Integer
forall a. Ratio a -> a
denominator Rational
q Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1 Bool -> Bool -> Bool
&& Integer -> Bool
forall a. Integral a => a -> Bool
even (Rational -> Integer
forall a. Ratio a -> a
numerator Rational
q) Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
x) Bool -> Bool -> Bool
|| Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x)
               )
            Bool -> Bool -> Bool
|| (Property
prop Property -> Property -> Bool
forall a. Eq a => a -> a -> Bool
== Property
NonNegative Bool -> Bool -> Bool
&& Rational -> Integer
forall a. Ratio a -> a
denominator Rational
q Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1 Bool -> Bool -> Bool
&& Integer -> Bool
forall a. Integral a => a -> Bool
even (Rational -> Integer
forall a. Ratio a -> a
numerator Rational
q) Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
e)
    Call String
"exp" [Expr
x] -> Property
prop Property -> [Property] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Property
Real, Property
Positive, Property
NonZero, Property
NonNegative] Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
x
    Call String
"abs" [Expr
x] ->
        (Property
prop Property -> [Property] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Property
Real, Property
NonNegative] Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
x)
            Bool -> Bool -> Bool
|| (Property
prop Property -> [Property] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Property
Positive, Property
NonZero] Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
x)
    Call String
"sqrt" [Expr
x] ->
        (Property
prop Property -> [Property] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Property
Real, Property
NonNegative] Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonNegative Expr
x)
            Bool -> Bool -> Bool
|| (Property
prop Property -> [Property] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Property
Positive, Property
NonZero] Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x)
    ExprNode
_ -> Property
prop Property -> Property -> Bool
forall a. Eq a => a -> a -> Bool
== Property
Real Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
e
  where
    implies :: Property -> Property -> Bool
implies Property
p Property
q =
        Property
p Property -> Property -> Bool
forall a. Eq a => a -> a -> Bool
== Property
q Bool -> Bool -> Bool
|| case Property
q of
            Property
Positive -> Property
p Property -> [Property] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Property
Real, Property
NonZero, Property
NonNegative]
            Property
Negative -> Property
p Property -> [Property] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Property
Real, Property
NonZero, Property
NonPositive]
            Property
Integral -> Property
p Property -> Property -> Bool
forall a. Eq a => a -> a -> Bool
== Property
Real
            Property
_ -> Bool
False

defined :: Assumptions -> Expr -> Bool
defined :: Assumptions -> Expr -> Bool
defined Assumptions
a (Expr SourcePos
_ ExprNode
n) = case ExprNode
n of
    Literal Rational
_ -> Bool
True
    Variable String
_ -> Bool
True
    Constant NamedConstant
_ -> Bool
True
    Unary UnaryOp
_ Expr
x -> Assumptions -> Expr -> Bool
defined Assumptions
a Expr
x
    Binary BinaryOp
Divide Expr
x Expr
y -> Assumptions -> Expr -> Bool
defined Assumptions
a Expr
x Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
y Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
y
    Binary BinaryOp
Power Expr
x (Expr SourcePos
_ (Literal Rational
q)) -> Assumptions -> Expr -> Bool
defined Assumptions
a Expr
x Bool -> Bool -> Bool
&& (if Rational -> Integer
forall a. Ratio a -> a
denominator Rational
q Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1 then Rational
q Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
0 Bool -> Bool -> Bool
|| Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
x else Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x)
    Binary BinaryOp
Power Expr
x Expr
y -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
y
    Binary BinaryOp
_ Expr
x Expr
y -> Assumptions -> Expr -> Bool
defined Assumptions
a Expr
x Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
y
    Call String
f [Expr]
xs
        | String
f
            String -> [String] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [ String
"sin"
                   , String
"cos"
                   , String
"exp"
                   , String
"abs"
                   , String
"signum"
                   , String
"floor"
                   , String
"ceil"
                   , String
"atan"
                   , String
"sinh"
                   , String
"cosh"
                   , String
"tanh"
                   , String
"min"
                   , String
"max"
                   ] ->
            (Expr -> Bool) -> [Expr] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Assumptions -> Expr -> Bool
defined Assumptions
a) [Expr]
xs Bool -> Bool -> Bool
&& [Expr] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Expr]
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== (if String
f String -> [String] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [String
"min", String
"max"] then Int
2 else Int
1)
    Call String
"log" [Expr
x] -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x
    Call String
"sqrt" [Expr
x] -> Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonNegative Expr
x
    Call String
_ [Expr]
_ -> Bool
False
    Let String
_ Expr
_ Expr
_ -> Bool
False

checked :: SymbolicOptions -> Expr -> Either ExprError Expr
checked :: SymbolicOptions -> Expr -> Either ExprError Expr
checked SymbolicOptions
o Expr
e = Limits -> Expr -> Either ExprError ()
validateExpression (SymbolicOptions -> Limits
symbolicLimits SymbolicOptions
o) Expr
e Either ExprError ()
-> Either ExprError Expr -> Either ExprError Expr
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right Expr
e
inputs :: SymbolicOptions -> [Expr] -> Either ExprError ()
inputs :: SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [Expr]
es = Limits -> Either ExprError ()
validateLimits (SymbolicOptions -> Limits
symbolicLimits SymbolicOptions
o) Either ExprError () -> Either ExprError () -> Either ExprError ()
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> (Expr -> Either ExprError ()) -> [Expr] -> Either ExprError ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (Limits -> Expr -> Either ExprError ()
validateExpression (SymbolicOptions -> Limits
symbolicLimits SymbolicOptions
o)) [Expr]
es
qvalue :: Expr -> Maybe Rational
qvalue :: Expr -> Maybe Rational
qvalue (Expr SourcePos
_ (Literal Rational
q)) = Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
q
qvalue Expr
_ = Maybe Rational
forall a. Maybe a
Nothing
plus, minus, times, divide, power :: Expr -> Expr -> Expr
plus :: Expr -> Expr -> Expr
plus = BinaryOp -> Expr -> Expr -> Expr
binary BinaryOp
Add
minus :: Expr -> Expr -> Expr
minus = BinaryOp -> Expr -> Expr -> Expr
binary BinaryOp
Subtract
times :: Expr -> Expr -> Expr
times = BinaryOp -> Expr -> Expr -> Expr
binary BinaryOp
Multiply
divide :: Expr -> Expr -> Expr
divide = BinaryOp -> Expr -> Expr -> Expr
binary BinaryOp
Divide
power :: Expr -> Expr -> Expr
power = BinaryOp -> Expr -> Expr -> Expr
binary BinaryOp
Power
neg :: Expr -> Expr
neg :: Expr -> Expr
neg Expr
e = case Expr -> Maybe Rational
qvalue Expr
e of Just Rational
q -> Rational -> Expr
literal (-Rational
q); Maybe Rational
Nothing -> UnaryOp -> Expr -> Expr
unary UnaryOp
Negate Expr
e

{- | Bottom-up folding and conservative canonical ring collection. It never
rewrites x/x or removes an unproved partial function from a zero product.
-}
normalize :: SymbolicOptions -> Expr -> Either ExprError Expr
normalize :: SymbolicOptions -> Expr -> Either ExprError Expr
normalize = SymbolicOptions -> Expr -> Either ExprError Expr
simplify

{- | Given SymbolicOptions then Expr, simplify within its resource limits and supplied
assumptions, preserving domain-sensitive conditions. Return ExprError for
invalid/budget-exceeding input/output.
-}
simplify :: SymbolicOptions -> Expr -> Either ExprError Expr
simplify :: SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o Expr
e = SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [Expr
e] Either ExprError ()
-> Either ExprError Expr -> Either ExprError Expr
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Expr -> Either ExprError Expr
go Expr
e Either ExprError Expr
-> (Expr -> Either ExprError Expr) -> Either ExprError Expr
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= SymbolicOptions -> Expr -> Either ExprError Expr
checked SymbolicOptions
o
  where
    a :: Assumptions
a = SymbolicOptions -> Assumptions
symbolicAssumptions SymbolicOptions
o
    go :: Expr -> Either ExprError Expr
go (Expr SourcePos
p ExprNode
n) = do
        result <- case ExprNode
n of
            Unary UnaryOp
Negate Expr
x -> do
                u <- Expr -> Either ExprError Expr
go Expr
x
                Right (case expressionNode u of Literal Rational
q -> Rational -> Expr
literal (-Rational
q); Unary UnaryOp
Negate Expr
y -> Expr
y; ExprNode
_ -> Expr -> Expr
neg Expr
u)
            Binary BinaryOp
op Expr
x Expr
y -> do
                u <- Expr -> Either ExprError Expr
go Expr
x
                v <- go y
                foldNode p op u v
            Call String
f [Expr]
xs -> do ys <- (Expr -> Either ExprError Expr)
-> [Expr] -> Either ExprError [Expr]
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 Expr -> Either ExprError Expr
go [Expr]
xs; Right (functionFold f ys)
            Let String
name Expr
x Expr
body -> do
                u <- Expr -> Either ExprError Expr
go Expr
x
                v <- go body
                if defined a u
                    then substitute (symbolicLimits o) [(name, u)] v >>= go
                    else Right (Expr p (Let name u v))
            ExprNode
_ -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (SourcePos -> ExprNode -> Expr
Expr SourcePos
p ExprNode
n)
        checked o result
    foldNode :: SourcePos -> BinaryOp -> Expr -> Expr -> Either a Expr
foldNode SourcePos
p BinaryOp
op Expr
u Expr
v = case (Expr -> Maybe Rational
qvalue Expr
u, Expr -> Maybe Rational
qvalue Expr
v) of
        (Just Rational
_, Just Rational
_) -> case Limits
-> Precision
-> Registry
-> Bindings
-> Expr
-> Either ExprError Value
evaluateExpression
            (SymbolicOptions -> Limits
symbolicLimits SymbolicOptions
o)
            Precision
ExactRational
            Registry
defaultRegistry
            Bindings
forall k a. Map k a
Map.empty
            (SourcePos -> ExprNode -> Expr
Expr SourcePos
p (BinaryOp -> Expr -> Expr -> ExprNode
Binary BinaryOp
op Expr
u Expr
v)) of
            Right (ExactValue Rational
q) -> Expr -> Either a Expr
forall a b. b -> Either a b
Right (Rational -> Expr
literal Rational
q)
            Either ExprError Value
_ -> BinaryOp -> Expr -> Expr -> Either a Expr
forall {a}. BinaryOp -> Expr -> Expr -> Either a Expr
fallback BinaryOp
op Expr
u Expr
v
        (Maybe Rational, Maybe Rational)
_ -> BinaryOp -> Expr -> Expr -> Either a Expr
forall {a}. BinaryOp -> Expr -> Expr -> Either a Expr
fallback BinaryOp
op Expr
u Expr
v
    fallback :: BinaryOp -> Expr -> Expr -> Either a Expr
fallback BinaryOp
op Expr
u Expr
v = Expr -> Either a Expr
forall a b. b -> Either a b
Right (Expr -> Either a Expr) -> Expr -> Either a Expr
forall a b. (a -> b) -> a -> b
$ case BinaryOp
op of
        BinaryOp
Add
            | Expr -> Maybe Rational
qvalue Expr
u Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
0 -> Expr
v
            | Expr -> Maybe Rational
qvalue Expr
v Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
0 -> Expr
u
            | Assumptions -> Expr -> Bool
defined Assumptions
a Expr
u Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
v -> Expr -> Expr
collect (Expr -> Expr -> Expr
plus Expr
u Expr
v)
        BinaryOp
Subtract
            | Expr -> Maybe Rational
qvalue Expr
v Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
0 -> Expr
u
            | Expr
u Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Expr
v Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
u -> Rational -> Expr
literal Rational
0
            | Assumptions -> Expr -> Bool
defined Assumptions
a Expr
u Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
v -> Expr -> Expr
collect (Expr -> Expr -> Expr
minus Expr
u Expr
v)
        BinaryOp
Multiply
            | Expr -> Maybe Rational
qvalue Expr
u Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
1 -> Expr
v
            | Expr -> Maybe Rational
qvalue Expr
v Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
1 -> Expr
u
            | Expr -> Maybe Rational
qvalue Expr
u Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
0 Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
v -> Rational -> Expr
literal Rational
0
            | Expr -> Maybe Rational
qvalue Expr
v Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
0 Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
u -> Rational -> Expr
literal Rational
0
            | Expr -> Maybe Rational
qvalue Expr
u Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just (-Rational
1) -> Expr -> Expr
neg Expr
v
            | Expr -> Maybe Rational
qvalue Expr
v Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just (-Rational
1) -> Expr -> Expr
neg Expr
u
            | Assumptions -> Expr -> Bool
defined Assumptions
a Expr
u Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined Assumptions
a Expr
v -> Expr -> Expr -> Expr
productCanonical Expr
u Expr
v
        BinaryOp
Divide
            | Expr -> Maybe Rational
qvalue Expr
v Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
1 -> Expr
u
            | Expr
u Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Expr
v Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
u -> Rational -> Expr
literal Rational
1
            | Expr -> Maybe Rational
qvalue Expr
u Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
0 Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
v -> Rational -> Expr
literal Rational
0
        BinaryOp
Power
            | Expr -> Maybe Rational
qvalue Expr
v Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
1 -> Expr
u
            | Expr -> Maybe Rational
qvalue Expr
v Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
0 Bool -> Bool -> Bool
&& Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
u -> Rational -> Expr
literal Rational
1
        BinaryOp
_ -> BinaryOp -> Expr -> Expr -> Expr
binary BinaryOp
op Expr
u Expr
v
    functionFold :: String -> [Expr] -> Expr
functionFold String
f [Expr]
xs = case (String
f, [Expr]
xs) of
        (String
"sin", [Expr
x]) | Expr -> Maybe Rational
qvalue Expr
x Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
0 -> Rational -> Expr
literal Rational
0
        (String
"cos", [Expr
x]) | Expr -> Maybe Rational
qvalue Expr
x Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
0 -> Rational -> Expr
literal Rational
1
        (String
"exp", [Expr
x]) | Expr -> Maybe Rational
qvalue Expr
x Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
0 -> Rational -> Expr
literal Rational
1
        (String
"log", [Expr
x]) | Expr -> Maybe Rational
qvalue Expr
x Maybe Rational -> Maybe Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Maybe Rational
forall a. a -> Maybe a
Just Rational
1 -> Rational -> Expr
literal Rational
0
        (String
"log", [Expr SourcePos
_ (Call String
"exp" [Expr
x])]) | Assumptions -> Expr -> Bool
defined Assumptions
a Expr
x -> Expr
x
        (String
"exp", [Expr SourcePos
_ (Call String
"log" [Expr
x])]) | Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x -> Expr
x
        (String
"abs", [Expr
x]) | Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonNegative Expr
x -> Expr
x
        (String
"abs", [Expr
x]) | Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonPositive Expr
x -> Expr -> Expr
neg Expr
x
        (String
"sqrt", [Expr
x]) -> case Limits
-> Precision
-> Registry
-> Bindings
-> Expr
-> Either ExprError Value
evaluateExpression (SymbolicOptions -> Limits
symbolicLimits SymbolicOptions
o) Precision
ExactRational Registry
defaultRegistry Bindings
forall k a. Map k a
Map.empty (String -> [Expr] -> Expr
call String
f [Expr]
xs) of
            Right (ExactValue Rational
q) -> Rational -> Expr
literal Rational
q
            Either ExprError Value
_ -> String -> [Expr] -> Expr
call String
"sqrt" [Expr
x]
        (String, [Expr])
_ -> String -> [Expr] -> Expr
call String
f [Expr]
xs
    collect :: Expr -> Expr
collect Expr
expr = [(Expr, Rational)] -> Expr
assemble (Map Expr Rational -> [(Expr, Rational)]
forall k a. Map k a -> [(k, a)]
Map.toAscList (Expr -> Map Expr Rational
terms Expr
expr))
    terms :: Expr -> Map Expr Rational
terms (Expr SourcePos
_ (Binary BinaryOp
Add Expr
x Expr
y)) = (Rational -> Rational -> Rational)
-> Map Expr Rational -> Map Expr Rational -> Map Expr Rational
forall k a. Ord k => (a -> a -> a) -> Map k a -> Map k a -> Map k a
Map.unionWith Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
(+) (Expr -> Map Expr Rational
terms Expr
x) (Expr -> Map Expr Rational
terms Expr
y)
    terms (Expr SourcePos
_ (Binary BinaryOp
Subtract Expr
x Expr
y)) = (Rational -> Rational -> Rational)
-> Map Expr Rational -> Map Expr Rational -> Map Expr Rational
forall k a. Ord k => (a -> a -> a) -> Map k a -> Map k a -> Map k a
Map.unionWith Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
(+) (Expr -> Map Expr Rational
terms Expr
x) ((Rational -> Rational) -> Map Expr Rational -> Map Expr Rational
forall a b k. (a -> b) -> Map k a -> Map k b
Map.map Rational -> Rational
forall a. Num a => a -> a
negate (Expr -> Map Expr Rational
terms Expr
y))
    terms (Expr SourcePos
_ (Unary UnaryOp
Negate Expr
x)) = (Rational -> Rational) -> Map Expr Rational -> Map Expr Rational
forall a b k. (a -> b) -> Map k a -> Map k b
Map.map Rational -> Rational
forall a. Num a => a -> a
negate (Expr -> Map Expr Rational
terms Expr
x)
    terms Expr
x = let (Rational
c, Expr
body) = Expr -> (Rational, Expr)
coefficient Expr
x in Expr -> Rational -> Map Expr Rational
forall k a. k -> a -> Map k a
Map.singleton Expr
body Rational
c
    coefficient :: Expr -> (Rational, Expr)
coefficient Expr
x = let (Rational
q, [Expr]
fs) = Expr -> (Rational, [Expr])
productParts Expr
x in (Rational
q, [Expr] -> Expr
productBody [Expr]
fs)
    productParts :: Expr -> (Rational, [Expr])
productParts Expr
x = case Expr -> ExprNode
expressionNode Expr
x of
        Literal Rational
q -> (Rational
q, [])
        Unary UnaryOp
Negate Expr
y -> let (Rational
q, [Expr]
fs) = Expr -> (Rational, [Expr])
productParts Expr
y in (-Rational
q, [Expr]
fs)
        Binary BinaryOp
Multiply Expr
u Expr
v -> let (Rational
q, [Expr]
fs) = Expr -> (Rational, [Expr])
productParts Expr
u; (Rational
r, [Expr]
gs) = Expr -> (Rational, [Expr])
productParts Expr
v in (Rational
q Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
r, [Expr]
fs [Expr] -> [Expr] -> [Expr]
forall a. [a] -> [a] -> [a]
++ [Expr]
gs)
        ExprNode
_ -> (Rational
1, [Expr
x])
    productBody :: [Expr] -> Expr
productBody [Expr]
fs = case [Expr] -> [Expr]
forall a. Ord a => [a] -> [a]
sort [Expr]
fs of [] -> Rational -> Expr
literal Rational
1; Expr
y : [Expr]
ys -> (Expr -> Expr -> Expr) -> Expr -> [Expr] -> Expr
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl' Expr -> Expr -> Expr
times Expr
y [Expr]
ys
    productCanonical :: Expr -> Expr -> Expr
productCanonical Expr
u Expr
v =
        let (Rational
q, [Expr]
fs) = Expr -> (Rational, [Expr])
productParts (Expr -> Expr -> Expr
times Expr
u Expr
v); body :: Expr
body = [Expr] -> Expr
productBody [Expr]
fs
         in if [Expr] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Expr]
fs
                then Rational -> Expr
literal Rational
q
                else if Rational
q Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
1 then Expr
body else if Rational
q Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== (-Rational
1) then Expr -> Expr
neg Expr
body else Expr -> Expr -> Expr
times (Rational -> Expr
literal Rational
q) Expr
body
    assemble :: [(Expr, Rational)] -> Expr
assemble [(Expr, Rational)]
entries = case ((Expr, Rational) -> Expr) -> [(Expr, Rational)] -> [Expr]
forall a b. (a -> b) -> [a] -> [b]
map (Expr, Rational) -> Expr
term (((Expr, Rational) -> Bool)
-> [(Expr, Rational)] -> [(Expr, Rational)]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0) (Rational -> Bool)
-> ((Expr, Rational) -> Rational) -> (Expr, Rational) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Expr, Rational) -> Rational
forall a b. (a, b) -> b
snd) [(Expr, Rational)]
entries) of
        [] -> Rational -> Expr
literal Rational
0
        Expr
x : [Expr]
xs -> (Expr -> Expr -> Expr) -> Expr -> [Expr] -> Expr
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl' Expr -> Expr -> Expr
plus Expr
x [Expr]
xs
    term :: (Expr, Rational) -> Expr
term (Expr
x, Rational
q)
        | Expr
x Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Expr
literal Rational
1 = Rational -> Expr
literal Rational
q
        | Rational
q Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
1 = Expr
x
        | Rational
q Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== (-Rational
1) = Expr -> Expr
neg Expr
x
        | Bool
otherwise = Expr -> Expr -> Expr
times (Rational -> Expr
literal Rational
q) Expr
x

{- | Alias of simplify, combining supported like terms under the same assumption and budget
checks.
-}
collectTerms :: SymbolicOptions -> Expr -> Either ExprError Expr
collectTerms :: SymbolicOptions -> Expr -> Either ExprError Expr
collectTerms = SymbolicOptions -> Expr -> Either ExprError Expr
simplify

{- | Distribute products over sums and nonnegative integer powers <=32, with a
checked output size at every convolution. Other powers remain symbolic.
-}
expand :: SymbolicOptions -> Expr -> Either ExprError Expr
expand :: SymbolicOptions -> Expr -> Either ExprError Expr
expand SymbolicOptions
o Expr
e = SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [Expr
e] Either ExprError ()
-> Either ExprError Expr -> Either ExprError Expr
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Expr -> Either ExprError Expr
go Expr
e Either ExprError Expr
-> (Expr -> Either ExprError Expr) -> Either ExprError Expr
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o
  where
    go :: Expr -> Either ExprError Expr
go x :: Expr
x@(Expr SourcePos
_ ExprNode
n) = case ExprNode
n of
        Binary BinaryOp
Add Expr
a Expr
b -> BinaryOp -> Expr -> Expr -> Either ExprError Expr
join BinaryOp
Add Expr
a Expr
b
        Binary BinaryOp
Subtract Expr
a Expr
b -> BinaryOp -> Expr -> Expr -> Either ExprError Expr
join BinaryOp
Subtract Expr
a Expr
b
        Binary BinaryOp
Multiply Expr
a Expr
b -> do u <- Expr -> Either ExprError Expr
go Expr
a; v <- go b; distribute u v
        Binary BinaryOp
Power Expr
a (Expr SourcePos
_ (Literal Rational
k)) | Rational -> Integer
forall a. Ratio a -> a
denominator Rational
k Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1 Bool -> Bool -> Bool
&& Rational
k Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
1 Bool -> Bool -> Bool
&& Rational
k Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
32 -> do
            u <- Expr -> Either ExprError Expr
go Expr
a
            foldM (\Expr
acc Integer
_ -> Expr -> Expr -> Either ExprError Expr
distribute Expr
acc Expr
u) (literal 1) [1 .. numerator k]
        Binary BinaryOp
op Expr
a Expr
b -> BinaryOp -> Expr -> Expr -> Either ExprError Expr
join BinaryOp
op Expr
a Expr
b
        Unary UnaryOp
op Expr
a -> Expr -> Either ExprError Expr
go Expr
a Either ExprError Expr
-> (Expr -> Either ExprError Expr) -> Either ExprError Expr
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= SymbolicOptions -> Expr -> Either ExprError Expr
checked SymbolicOptions
o (Expr -> Either ExprError Expr)
-> (Expr -> Expr) -> Expr -> Either ExprError Expr
forall b c a. (b -> c) -> (a -> b) -> a -> c
. UnaryOp -> Expr -> Expr
unary UnaryOp
op
        Call String
name [Expr]
xs -> (Expr -> Either ExprError Expr)
-> [Expr] -> Either ExprError [Expr]
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 Expr -> Either ExprError Expr
go [Expr]
xs Either ExprError [Expr]
-> ([Expr] -> Either ExprError Expr) -> Either ExprError Expr
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= SymbolicOptions -> Expr -> Either ExprError Expr
checked SymbolicOptions
o (Expr -> Either ExprError Expr)
-> ([Expr] -> Expr) -> [Expr] -> Either ExprError Expr
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> [Expr] -> Expr
call String
name
        Let String
binding Expr
a Expr
b -> do u <- Expr -> Either ExprError Expr
go Expr
a; v <- go b; checked o (letIn binding u v)
        ExprNode
_ -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right Expr
x
    join :: BinaryOp -> Expr -> Expr -> Either ExprError Expr
join BinaryOp
op Expr
a Expr
b = do u <- Expr -> Either ExprError Expr
go Expr
a; v <- go b; checked o (binary op u v)
    distribute :: Expr -> Expr -> Either ExprError Expr
distribute (Expr SourcePos
_ (Binary BinaryOp
Add Expr
a Expr
b)) Expr
c = do x <- Expr -> Expr -> Either ExprError Expr
distribute Expr
a Expr
c; y <- distribute b c; checked o (plus x y)
    distribute (Expr SourcePos
_ (Binary BinaryOp
Subtract Expr
a Expr
b)) Expr
c = do x <- Expr -> Expr -> Either ExprError Expr
distribute Expr
a Expr
c; y <- distribute b c; checked o (minus x y)
    distribute Expr
a (Expr SourcePos
_ (Binary BinaryOp
Add Expr
b Expr
c)) = do x <- Expr -> Expr -> Either ExprError Expr
distribute Expr
a Expr
b; y <- distribute a c; checked o (plus x y)
    distribute Expr
a (Expr SourcePos
_ (Binary BinaryOp
Subtract Expr
b Expr
c)) = do x <- Expr -> Expr -> Either ExprError Expr
distribute Expr
a Expr
b; y <- distribute a c; checked o (minus x y)
    distribute Expr
a Expr
b = SymbolicOptions -> Expr -> Either ExprError Expr
checked SymbolicOptions
o (Expr -> Expr -> Expr
times Expr
a Expr
b)

{- | Alpha-renaming makes simultaneous substitution capture-safe. Replacement
expressions are inserted once; substitution does not recurse into them.
-}
substitute :: Limits -> [(String, Expr)] -> Expr -> Either ExprError Expr
substitute :: Limits -> [(String, Expr)] -> Expr -> Either ExprError Expr
substitute Limits
limits [(String, Expr)]
replacements Expr
expression = do
    Limits -> Expr -> Either ExprError ()
validateExpression Limits
limits Expr
expression
    if [(String, Expr)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [(String, Expr)] -> [(String, Expr)]
forall a. Int -> [a] -> [a]
take (Limits -> Int
maximumNodes Limits
limits Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [(String, Expr)]
replacements) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Limits -> Int
maximumNodes Limits
limits
        then ErrorKind -> Expr -> String -> Either ExprError ()
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
MemoryLimit Expr
expression String
"substitution table limit"
        else () -> Either ExprError ()
forall a b. b -> Either a b
Right ()
    ((String, Expr) -> Either ExprError ())
-> [(String, Expr)] -> Either ExprError ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_
        (\(String
n, Expr
e) -> Limits -> Expr -> Either ExprError ()
validateExpression Limits
limits (String -> Expr
variable String
n) Either ExprError () -> Either ExprError () -> Either ExprError ()
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Limits -> Expr -> Either ExprError ()
validateExpression Limits
limits Expr
e)
        [(String, Expr)]
replacements
    Map String Expr -> Expr -> Either ExprError Expr
go ([(String, Expr)] -> Map String Expr
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList [(String, Expr)]
replacements) Expr
expression Either ExprError Expr
-> (Expr -> Either ExprError Expr) -> Either ExprError Expr
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= \Expr
e -> Limits -> Expr -> Either ExprError ()
validateExpression Limits
limits Expr
e Either ExprError ()
-> Either ExprError Expr -> Either ExprError Expr
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right Expr
e
  where
    go :: Map String Expr -> Expr -> Either ExprError Expr
go Map String Expr
env e :: Expr
e@(Expr SourcePos
p ExprNode
node) = case ExprNode
node of
        Variable String
n -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Expr -> String -> Map String Expr -> Expr
forall k a. Ord k => a -> k -> Map k a -> a
Map.findWithDefault Expr
e String
n Map String Expr
env)
        Unary UnaryOp
op Expr
a -> SourcePos -> ExprNode -> Expr
Expr SourcePos
p (ExprNode -> Expr) -> (Expr -> ExprNode) -> Expr -> Expr
forall b c a. (b -> c) -> (a -> b) -> a -> c
. UnaryOp -> Expr -> ExprNode
Unary UnaryOp
op (Expr -> Expr) -> Either ExprError Expr -> Either ExprError Expr
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Map String Expr -> Expr -> Either ExprError Expr
go Map String Expr
env Expr
a
        Binary BinaryOp
op Expr
a Expr
b -> do x <- Map String Expr -> Expr -> Either ExprError Expr
go Map String Expr
env Expr
a; y <- go env b; bounded (Expr p (Binary op x y))
        Call String
name [Expr]
xs -> (Expr -> Either ExprError Expr)
-> [Expr] -> Either ExprError [Expr]
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 (Map String Expr -> Expr -> Either ExprError Expr
go Map String Expr
env) [Expr]
xs Either ExprError [Expr]
-> ([Expr] -> Either ExprError Expr) -> Either ExprError Expr
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Expr -> Either ExprError Expr
bounded (Expr -> Either ExprError Expr)
-> ([Expr] -> Expr) -> [Expr] -> Either ExprError Expr
forall b c a. (b -> c) -> (a -> b) -> a -> c
. SourcePos -> ExprNode -> Expr
Expr SourcePos
p (ExprNode -> Expr) -> ([Expr] -> ExprNode) -> [Expr] -> Expr
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> [Expr] -> ExprNode
Call String
name
        Let String
n Expr
a Expr
body -> do
            x <- Map String Expr -> Expr -> Either ExprError Expr
go Map String Expr
env Expr
a
            let active = String -> Map String Expr -> Map String Expr
forall k a. Ord k => k -> Map k a -> Map k a
Map.delete String
n Map String Expr
env
                inserted = [Set String] -> Set String
forall (f :: * -> *) a. (Foldable f, Ord a) => f (Set a) -> Set a
Set.unions ((Expr -> Set String) -> [Expr] -> [Set String]
forall a b. (a -> b) -> [a] -> [b]
map ([String] -> Set String
forall a. Ord a => [a] -> Set a
Set.fromList ([String] -> Set String)
-> (Expr -> [String]) -> Expr -> Set String
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Expr -> [String]
freeVariables) (Map String Expr -> [Expr]
forall k a. Map k a -> [a]
Map.elems Map String Expr
active))
            if Set.member n inserted
                then do
                    let occupied = [Set String] -> Set String
forall (f :: * -> *) a. (Foldable f, Ord a) => f (Set a) -> Set a
Set.unions [Set String
inserted, [String] -> Set String
forall a. Ord a => [a] -> Set a
Set.fromList (Expr -> [String]
allNames Expr
body), Map String Expr -> Set String
forall k a. Map k a -> Set k
Map.keysSet Map String Expr
active]
                        fresh = Set String -> String -> Int -> String
forall {t}. (Num t, Show t) => Set String -> String -> t -> String
freshName Set String
occupied String
n (Int
0 :: Int)
                    renamed <- go (Map.singleton n (variable fresh)) body
                    y <- go active renamed
                    bounded (Expr p (Let fresh x y))
                else do y <- go active body; bounded (Expr p (Let n x y))
        ExprNode
_ -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right Expr
e
    bounded :: Expr -> Either ExprError Expr
bounded Expr
e = Limits -> Expr -> Either ExprError ()
validateExpression Limits
limits Expr
e Either ExprError ()
-> Either ExprError Expr -> Either ExprError Expr
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right Expr
e
    freshName :: Set String -> String -> t -> String
freshName Set String
occupied String
base t
i =
        let n :: String
n = String
base String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"_" String -> ShowS
forall a. [a] -> [a] -> [a]
++ t -> String
forall a. Show a => a -> String
show t
i
         in if String -> Set String -> Bool
forall a. Ord a => a -> Set a -> Bool
Set.member String
n Set String
occupied then Set String -> String -> t -> String
freshName Set String
occupied String
base (t
i t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) else String
n
    allNames :: Expr -> [String]
allNames (Expr SourcePos
_ ExprNode
node) = case ExprNode
node of
        Variable String
n -> [String
n]
        Unary UnaryOp
_ Expr
a -> Expr -> [String]
allNames Expr
a
        Binary BinaryOp
_ Expr
a Expr
b -> Expr -> [String]
allNames Expr
a [String] -> [String] -> [String]
forall a. [a] -> [a] -> [a]
++ Expr -> [String]
allNames Expr
b
        Call String
n [Expr]
xs -> String
n String -> [String] -> [String]
forall a. a -> [a] -> [a]
: (Expr -> [String]) -> [Expr] -> [String]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap Expr -> [String]
allNames [Expr]
xs
        Let String
n Expr
a Expr
b -> String
n String -> [String] -> [String]
forall a. a -> [a] -> [a]
: Expr -> [String]
allNames Expr
a [String] -> [String] -> [String]
forall a. [a] -> [a] -> [a]
++ Expr -> [String]
allNames Expr
b
        ExprNode
_ -> []

{- | Symbolic rewrite pattern with literal/variable/operator/call structure and binding
placeholders.
-}
data Pattern
    = Any String
    | Exactly Expr
    | UnaryPattern UnaryOp Pattern
    | BinaryPattern BinaryOp Pattern Pattern
    | CallPattern String [Pattern]
    deriving (Pattern -> Pattern -> Bool
(Pattern -> Pattern -> Bool)
-> (Pattern -> Pattern -> Bool) -> Eq Pattern
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Pattern -> Pattern -> Bool
== :: Pattern -> Pattern -> Bool
$c/= :: Pattern -> Pattern -> Bool
/= :: Pattern -> Pattern -> Bool
Eq, Int -> Pattern -> ShowS
[Pattern] -> ShowS
Pattern -> String
(Int -> Pattern -> ShowS)
-> (Pattern -> String) -> ([Pattern] -> ShowS) -> Show Pattern
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Pattern -> ShowS
showsPrec :: Int -> Pattern -> ShowS
$cshow :: Pattern -> String
show :: Pattern -> String
$cshowList :: [Pattern] -> ShowS
showList :: [Pattern] -> ShowS
Show)

-- | Repeated wildcard names must match structurally identical subtrees.
matchPattern :: Limits -> Pattern -> Expr -> Either ExprError (Maybe [(String, Expr)])
matchPattern :: Limits
-> Pattern -> Expr -> Either ExprError (Maybe [(String, Expr)])
matchPattern Limits
limits Pattern
pat Expr
expression = do
    Limits -> Expr -> Either ExprError ()
validateExpression Limits
limits Expr
expression
    (_, result) <- Int
-> Int
-> Map String Expr
-> Pattern
-> Expr
-> Either ExprError (Int, Maybe (Map String Expr))
forall {t}.
(Ord t, Num t) =>
Int
-> t
-> Map String Expr
-> Pattern
-> Expr
-> Either ExprError (t, Maybe (Map String Expr))
go Int
1 (Limits -> Int
maximumSteps Limits
limits) Map String Expr
forall k a. Map k a
Map.empty Pattern
pat Expr
expression
    Right (Map.toAscList <$> result)
  where
    go :: Int
-> t
-> Map String Expr
-> Pattern
-> Expr
-> Either ExprError (t, Maybe (Map String Expr))
go Int
d t
fuel Map String Expr
env Pattern
patternNode Expr
expr
        | t
fuel t -> t -> Bool
forall a. Ord a => a -> a -> Bool
<= t
0 = ErrorKind
-> Expr -> String -> Either ExprError (t, Maybe (Map String Expr))
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
StepLimit Expr
expr String
"pattern matching budget"
        | Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Limits -> Int
maximumDepth Limits
limits = ErrorKind
-> Expr -> String -> Either ExprError (t, Maybe (Map String Expr))
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
DepthLimit Expr
expr String
"pattern matching depth"
        | Bool
otherwise = case Pattern
patternNode of
            Any String
name -> do
                Limits -> Expr -> Either ExprError ()
validateExpression Limits
limits (String -> Expr
variable String
name)
                (t, Maybe (Map String Expr))
-> Either ExprError (t, Maybe (Map String Expr))
forall a b. b -> Either a b
Right
                    ( t
fuel t -> t -> t
forall a. Num a => a -> a -> a
- t
1
                    , case String -> Map String Expr -> Maybe Expr
forall k a. Ord k => k -> Map k a -> Maybe a
Map.lookup String
name Map String Expr
env of
                        Maybe Expr
Nothing -> Map String Expr -> Maybe (Map String Expr)
forall a. a -> Maybe a
Just (String -> Expr -> Map String Expr -> Map String Expr
forall k a. Ord k => k -> a -> Map k a -> Map k a
Map.insert String
name Expr
expr Map String Expr
env)
                        Just Expr
x -> if Expr
x Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Expr
expr then Map String Expr -> Maybe (Map String Expr)
forall a. a -> Maybe a
Just Map String Expr
env else Maybe (Map String Expr)
forall a. Maybe a
Nothing
                    )
            Exactly Expr
expected ->
                Limits -> Expr -> Either ExprError ()
validateExpression Limits
limits Expr
expected
                    Either ExprError ()
-> Either ExprError (t, Maybe (Map String Expr))
-> Either ExprError (t, Maybe (Map String Expr))
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> (t, Maybe (Map String Expr))
-> Either ExprError (t, Maybe (Map String Expr))
forall a b. b -> Either a b
Right (t
fuel t -> t -> t
forall a. Num a => a -> a -> a
- t
1, if Expr
expected Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Expr
expr then Map String Expr -> Maybe (Map String Expr)
forall a. a -> Maybe a
Just Map String Expr
env else Maybe (Map String Expr)
forall a. Maybe a
Nothing)
            UnaryPattern UnaryOp
op Pattern
a -> case Expr -> ExprNode
expressionNode Expr
expr of
                Unary UnaryOp
op1 Expr
x | UnaryOp
op UnaryOp -> UnaryOp -> Bool
forall a. Eq a => a -> a -> Bool
== UnaryOp
op1 -> Int
-> t
-> Map String Expr
-> Pattern
-> Expr
-> Either ExprError (t, Maybe (Map String Expr))
go (Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (t
fuel t -> t -> t
forall a. Num a => a -> a -> a
- t
1) Map String Expr
env Pattern
a Expr
x
                ExprNode
_ -> (t, Maybe (Map String Expr))
-> Either ExprError (t, Maybe (Map String Expr))
forall a b. b -> Either a b
Right (t
fuel t -> t -> t
forall a. Num a => a -> a -> a
- t
1, Maybe (Map String Expr)
forall a. Maybe a
Nothing)
            BinaryPattern BinaryOp
op Pattern
a Pattern
b -> case Expr -> ExprNode
expressionNode Expr
expr of
                Binary BinaryOp
op1 Expr
x Expr
y | BinaryOp
op BinaryOp -> BinaryOp -> Bool
forall a. Eq a => a -> a -> Bool
== BinaryOp
op1 -> do
                    (remaining, first) <- Int
-> t
-> Map String Expr
-> Pattern
-> Expr
-> Either ExprError (t, Maybe (Map String Expr))
go (Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (t
fuel t -> t -> t
forall a. Num a => a -> a -> a
- t
1) Map String Expr
env Pattern
a Expr
x
                    case first of { Maybe (Map String Expr)
Nothing -> (t, Maybe (Map String Expr))
-> Either ExprError (t, Maybe (Map String Expr))
forall a b. b -> Either a b
Right (t
remaining, Maybe (Map String Expr)
forall a. Maybe a
Nothing); Just Map String Expr
next -> Int
-> t
-> Map String Expr
-> Pattern
-> Expr
-> Either ExprError (t, Maybe (Map String Expr))
go (Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) t
remaining Map String Expr
next Pattern
b Expr
y }
                ExprNode
_ -> (t, Maybe (Map String Expr))
-> Either ExprError (t, Maybe (Map String Expr))
forall a b. b -> Either a b
Right (t
fuel t -> t -> t
forall a. Num a => a -> a -> a
- t
1, Maybe (Map String Expr)
forall a. Maybe a
Nothing)
            CallPattern String
name [Pattern]
ps -> case Expr -> ExprNode
expressionNode Expr
expr of
                Call String
name1 [Expr]
xs
                    | String
name String -> String -> Bool
forall a. Eq a => a -> a -> Bool
== String
name1 Bool -> Bool -> Bool
&& [Pattern] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [Pattern] -> [Pattern]
forall a. Int -> [a] -> [a]
take (Limits -> Int
maximumNodes Limits
limits Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Pattern]
ps) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Expr] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Expr]
xs ->
                        ((t, Maybe (Map String Expr))
 -> (Pattern, Expr)
 -> Either ExprError (t, Maybe (Map String Expr)))
-> (t, Maybe (Map String Expr))
-> [(Pattern, Expr)]
-> Either ExprError (t, Maybe (Map String Expr))
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM (Int
-> (t, Maybe (Map String Expr))
-> (Pattern, Expr)
-> Either ExprError (t, Maybe (Map String Expr))
step (Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)) (t
fuel t -> t -> t
forall a. Num a => a -> a -> a
- t
1, Map String Expr -> Maybe (Map String Expr)
forall a. a -> Maybe a
Just Map String Expr
env) ([Pattern] -> [Expr] -> [(Pattern, Expr)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Pattern]
ps [Expr]
xs)
                ExprNode
_ -> (t, Maybe (Map String Expr))
-> Either ExprError (t, Maybe (Map String Expr))
forall a b. b -> Either a b
Right (t
fuel t -> t -> t
forall a. Num a => a -> a -> a
- t
1, Maybe (Map String Expr)
forall a. Maybe a
Nothing)
    step :: Int
-> (t, Maybe (Map String Expr))
-> (Pattern, Expr)
-> Either ExprError (t, Maybe (Map String Expr))
step Int
_ (t
fuel, Maybe (Map String Expr)
Nothing) (Pattern, Expr)
_ = (t, Maybe (Map String Expr))
-> Either ExprError (t, Maybe (Map String Expr))
forall a b. b -> Either a b
Right (t
fuel, Maybe (Map String Expr)
forall a. Maybe a
Nothing)
    step Int
d (t
fuel, Just Map String Expr
env) (Pattern
p, Expr
e) = Int
-> t
-> Map String Expr
-> Pattern
-> Expr
-> Either ExprError (t, Maybe (Map String Expr))
go Int
d t
fuel Map String Expr
env Pattern
p Expr
e

{- | Pattern,replacement expression,and property requirements needed for a valid bounded
symbolic rewrite.
-}
data RewriteRule = RewriteRule
    { RewriteRule -> Pattern
rewritePattern :: Pattern
    -- ^ Pattern to match before applying this rewrite rule.
    , RewriteRule -> Expr
rewriteReplacement :: Expr
    -- ^ Replacement expression, with matched bindings substituted.
    , RewriteRule -> [(String, Property)]
rewriteRequirements :: [(String, Property)]
    -- ^ Variable properties required to justify the rewrite.
    }
    deriving (RewriteRule -> RewriteRule -> Bool
(RewriteRule -> RewriteRule -> Bool)
-> (RewriteRule -> RewriteRule -> Bool) -> Eq RewriteRule
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RewriteRule -> RewriteRule -> Bool
== :: RewriteRule -> RewriteRule -> Bool
$c/= :: RewriteRule -> RewriteRule -> Bool
/= :: RewriteRule -> RewriteRule -> Bool
Eq, Int -> RewriteRule -> ShowS
[RewriteRule] -> ShowS
RewriteRule -> String
(Int -> RewriteRule -> ShowS)
-> (RewriteRule -> String)
-> ([RewriteRule] -> ShowS)
-> Show RewriteRule
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> RewriteRule -> ShowS
showsPrec :: Int -> RewriteRule -> ShowS
$cshow :: RewriteRule -> String
show :: RewriteRule -> String
$cshowList :: [RewriteRule] -> ShowS
showList :: [RewriteRule] -> ShowS
Show)

{- | Apply one caller-supplied rule at the root only. Rules are trusted algebraic
assertions; requirements are checked against the matched expressions.
-}
applyRewrite :: SymbolicOptions -> RewriteRule -> Expr -> Either ExprError Expr
applyRewrite :: SymbolicOptions -> RewriteRule -> Expr -> Either ExprError Expr
applyRewrite SymbolicOptions
o RewriteRule
rule Expr
e = do
    result <- Limits
-> Pattern -> Expr -> Either ExprError (Maybe [(String, Expr)])
matchPattern (SymbolicOptions -> Limits
symbolicLimits SymbolicOptions
o) (RewriteRule -> Pattern
rewritePattern RewriteRule
rule) Expr
e
    case result of
        Maybe [(String, Expr)]
Nothing -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right Expr
e
        Just [(String, Expr)]
bindings ->
            if ((String, Property) -> Bool) -> [(String, Property)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all
                (\(String
n, Property
p) -> Bool -> (Expr -> Bool) -> Maybe Expr -> Bool
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Bool
False (Assumptions -> Property -> Expr -> Bool
proves (SymbolicOptions -> Assumptions
symbolicAssumptions SymbolicOptions
o) Property
p) (String -> [(String, Expr)] -> Maybe Expr
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup String
n [(String, Expr)]
bindings))
                (RewriteRule -> [(String, Property)]
rewriteRequirements RewriteRule
rule)
                then Limits -> [(String, Expr)] -> Expr -> Either ExprError Expr
substitute (SymbolicOptions -> Limits
symbolicLimits SymbolicOptions
o) [(String, Expr)]
bindings (RewriteRule -> Expr
rewriteReplacement RewriteRule
rule)
                else ErrorKind -> Expr -> String -> Either ExprError Expr
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
DomainFailure Expr
e String
"rewrite assumptions are not established"

{- | Differentiate supported elementary real functions. Noninteger powers and
logarithms require a positive argument; abs requires a nonzero argument.
-}
differentiate :: SymbolicOptions -> String -> Expr -> Either ExprError Expr
differentiate :: SymbolicOptions -> String -> Expr -> Either ExprError Expr
differentiate SymbolicOptions
o String
name Expr
expression = do
    SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [String -> Expr
variable String
name, Expr
expression]
    Expr -> Either ExprError Expr
go Expr
expression Either ExprError Expr
-> (Expr -> Either ExprError Expr) -> Either ExprError Expr
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o
  where
    a :: Assumptions
a = SymbolicOptions -> Assumptions
symbolicAssumptions SymbolicOptions
o
    go :: Expr -> Either ExprError Expr
go e :: Expr
e@(Expr SourcePos
_ ExprNode
node) = do
        result <- case ExprNode
node of
            Literal Rational
_ -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Rational -> Expr
literal Rational
0)
            Constant NamedConstant
_ -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Rational -> Expr
literal Rational
0)
            Variable String
n -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Rational -> Expr
literal (if String
n String -> String -> Bool
forall a. Eq a => a -> a -> Bool
== String
name then Rational
1 else Rational
0))
            Unary UnaryOp
Negate Expr
x -> Expr -> Expr
neg (Expr -> Expr) -> Either ExprError Expr -> Either ExprError Expr
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Expr -> Either ExprError Expr
go Expr
x
            Binary BinaryOp
Add Expr
x Expr
y -> Expr -> Expr -> Expr
plus (Expr -> Expr -> Expr)
-> Either ExprError Expr -> Either ExprError (Expr -> Expr)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Expr -> Either ExprError Expr
go Expr
x Either ExprError (Expr -> Expr)
-> Either ExprError Expr -> Either ExprError Expr
forall a b.
Either ExprError (a -> b)
-> Either ExprError a -> Either ExprError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Expr -> Either ExprError Expr
go Expr
y
            Binary BinaryOp
Subtract Expr
x Expr
y -> Expr -> Expr -> Expr
minus (Expr -> Expr -> Expr)
-> Either ExprError Expr -> Either ExprError (Expr -> Expr)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Expr -> Either ExprError Expr
go Expr
x Either ExprError (Expr -> Expr)
-> Either ExprError Expr -> Either ExprError Expr
forall a b.
Either ExprError (a -> b)
-> Either ExprError a -> Either ExprError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Expr -> Either ExprError Expr
go Expr
y
            Binary BinaryOp
Multiply Expr
x Expr
y -> do dx <- Expr -> Either ExprError Expr
go Expr
x; dy <- go y; Right (plus (times dx y) (times x dy))
            Binary BinaryOp
Divide Expr
x Expr
y -> do dx <- Expr -> Either ExprError Expr
go Expr
x; dy <- go y; Right (divide (minus (times dx y) (times x dy)) (power y (literal 2)))
            Binary BinaryOp
Power Expr
x (Expr SourcePos
_ (Literal Rational
n)) | Rational -> Integer
forall a. Ratio a -> a
denominator Rational
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1 Bool -> Bool -> Bool
|| Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x -> do
                dx <- Expr -> Either ExprError Expr
go Expr
x
                if n == 0
                    then
                        if proves a NonZero x then Right (literal 0) else unsupported e "zero power requires nonzero base"
                    else if n == 1 then Right dx else Right (times (times (literal n) (power x (literal (n - 1)))) dx)
            Binary BinaryOp
Power Expr
x Expr
y | Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x -> do
                dx <- Expr -> Either ExprError Expr
go Expr
x
                dy <- go y
                Right (times e (plus (times dy (call "log" [x])) (divide (times y dx) x)))
            Call String
f [Expr
x] -> do
                dx <- Expr -> Either ExprError Expr
go Expr
x
                outer <- case f of
                    String
"sin" -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (String -> [Expr] -> Expr
call String
"cos" [Expr
x])
                    String
"cos" -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Expr -> Expr
neg (String -> [Expr] -> Expr
call String
"sin" [Expr
x]))
                    String
"tan" -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Expr -> Expr -> Expr
divide (Rational -> Expr
literal Rational
1) (Expr -> Expr -> Expr
power (String -> [Expr] -> Expr
call String
"cos" [Expr
x]) (Rational -> Expr
literal Rational
2)))
                    String
"exp" -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right Expr
e
                    String
"log" | Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Expr -> Expr -> Expr
divide (Rational -> Expr
literal Rational
1) Expr
x)
                    String
"sqrt" | Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
Positive Expr
x -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Expr -> Expr -> Expr
divide (Rational -> Expr
literal Rational
1) (Expr -> Expr -> Expr
times (Rational -> Expr
literal Rational
2) Expr
e))
                    String
"abs" | Assumptions -> Property -> Expr -> Bool
proves Assumptions
a Property
NonZero Expr
x -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (String -> [Expr] -> Expr
call String
"signum" [Expr
x])
                    String
"atan" -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Expr -> Expr -> Expr
divide (Rational -> Expr
literal Rational
1) (Expr -> Expr -> Expr
plus (Rational -> Expr
literal Rational
1) (Expr -> Expr -> Expr
power Expr
x (Rational -> Expr
literal Rational
2))))
                    String
"sinh" -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (String -> [Expr] -> Expr
call String
"cosh" [Expr
x])
                    String
"cosh" -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (String -> [Expr] -> Expr
call String
"sinh" [Expr
x])
                    String
"tanh" -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Expr -> Expr -> Expr
minus (Rational -> Expr
literal Rational
1) (Expr -> Expr -> Expr
power Expr
e (Rational -> Expr
literal Rational
2)))
                    String
_ -> Expr -> String -> Either ExprError Expr
forall a. Expr -> String -> Either ExprError a
unsupported Expr
e String
"unsupported derivative or missing domain assumption"
                Right (times outer dx)
            Let String
n Expr
x Expr
body | Assumptions -> Expr -> Bool
defined Assumptions
a Expr
x -> Limits -> [(String, Expr)] -> Expr -> Either ExprError Expr
substitute (SymbolicOptions -> Limits
symbolicLimits SymbolicOptions
o) [(String
n, Expr
x)] Expr
body Either ExprError Expr
-> (Expr -> Either ExprError Expr) -> Either ExprError Expr
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Expr -> Either ExprError Expr
go
            Let String
_ Expr
_ Expr
_ -> Expr -> String -> Either ExprError Expr
forall a. Expr -> String -> Either ExprError a
unsupported Expr
e String
"derivative of strict let requires a defined initializer"
            ExprNode
_ -> Expr -> String -> Either ExprError Expr
forall a. Expr -> String -> Either ExprError a
unsupported Expr
e String
"unsupported symbolic derivative"
        checked o result

{- | Given symbolic options,ordered variable names,then scalar expression, return symbolic
partial derivatives in that order. At most 32 variables are supported;
expression/domain/resource failures return ExprError.
-}
gradient :: SymbolicOptions -> [String] -> Expr -> Either ExprError [Expr]
gradient :: SymbolicOptions -> [String] -> Expr -> Either ExprError [Expr]
gradient SymbolicOptions
o [String]
names Expr
e = SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [Expr
e] Either ExprError () -> Either ExprError () -> Either ExprError ()
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> SymbolicOptions -> [String] -> Either ExprError ()
boundedVariables SymbolicOptions
o [String]
names Either ExprError ()
-> Either ExprError [Expr] -> Either ExprError [Expr]
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> (String -> Either ExprError Expr)
-> [String] -> Either ExprError [Expr]
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 (\String
n -> SymbolicOptions -> String -> Expr -> Either ExprError Expr
differentiate SymbolicOptions
o String
n Expr
e) [String]
names

{- | Given options,ordered variable names,then output-expression list, return one gradient
row per output. Validate variable/expression budgets; failures return ExprError.
-}
jacobian :: SymbolicOptions -> [String] -> [Expr] -> Either ExprError [[Expr]]
jacobian :: SymbolicOptions -> [String] -> [Expr] -> Either ExprError [[Expr]]
jacobian SymbolicOptions
o [String]
names [Expr]
es = SymbolicOptions -> [String] -> Either ExprError ()
boundedVariables SymbolicOptions
o [String]
names Either ExprError () -> Either ExprError () -> Either ExprError ()
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> SymbolicOptions -> [Expr] -> Either ExprError ()
boundedExpressions SymbolicOptions
o [Expr]
es Either ExprError ()
-> Either ExprError [[Expr]] -> Either ExprError [[Expr]]
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> (Expr -> Either ExprError [Expr])
-> [Expr] -> Either ExprError [[Expr]]
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 (SymbolicOptions -> [String] -> Expr -> Either ExprError [Expr]
gradient SymbolicOptions
o [String]
names) [Expr]
es

{- | Given options,ordered variables,then scalar Expr, return derivatives of each gradient
component in the same row/column variable order. Validation/differentiation/resource
failures propagate.
-}
hessian :: SymbolicOptions -> [String] -> Expr -> Either ExprError [[Expr]]
hessian :: SymbolicOptions -> [String] -> Expr -> Either ExprError [[Expr]]
hessian SymbolicOptions
o [String]
names Expr
e = SymbolicOptions -> [String] -> Expr -> Either ExprError [Expr]
gradient SymbolicOptions
o [String]
names Expr
e Either ExprError [Expr]
-> ([Expr] -> Either ExprError [[Expr]])
-> Either ExprError [[Expr]]
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= (Expr -> Either ExprError [Expr])
-> [Expr] -> Either ExprError [[Expr]]
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 (SymbolicOptions -> [String] -> Expr -> Either ExprError [Expr]
gradient SymbolicOptions
o [String]
names)

boundedVariables :: SymbolicOptions -> [String] -> Either ExprError ()
boundedVariables :: SymbolicOptions -> [String] -> Either ExprError ()
boundedVariables SymbolicOptions
o [String]
names
    | [String] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [String] -> [String]
forall a. Int -> [a] -> [a]
take Int
33 [String]
names) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
32 = ErrorKind -> Expr -> String -> Either ExprError ()
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
MemoryLimit (Rational -> Expr
literal Rational
0) String
"at most 32 symbolic variables"
    | [String] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([String] -> [String]
forall a. Eq a => [a] -> [a]
nub [String]
names) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [String] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [String]
names =
        ErrorKind -> Expr -> String -> Either ExprError ()
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
DomainFailure (Rational -> Expr
literal Rational
0) String
"duplicate symbolic variables"
    | Bool
otherwise = SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o ((String -> Expr) -> [String] -> [Expr]
forall a b. (a -> b) -> [a] -> [b]
map String -> Expr
variable [String]
names)
boundedExpressions :: SymbolicOptions -> [Expr] -> Either ExprError ()
boundedExpressions :: SymbolicOptions -> [Expr] -> Either ExprError ()
boundedExpressions SymbolicOptions
o [Expr]
es
    | [Expr] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [Expr] -> [Expr]
forall a. Int -> [a] -> [a]
take Int
65 [Expr]
es) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
64 = ErrorKind -> Expr -> String -> Either ExprError ()
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
MemoryLimit (Rational -> Expr
literal Rational
0) String
"at most 64 symbolic outputs"
    | Bool
otherwise = SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [Expr]
es

-- Polynomial admission checks size/degree/coefficient bits at every operation.

{- | Given options,variable name,then Expr, convert a supported univariate exact Rational
polynomial. Other variables/functions, invalid exponents or resource-limit violations
return ExprError.
-}
toPolynomial :: SymbolicOptions -> String -> Expr -> Either ExprError (P.Polynomial Rational)
toPolynomial :: SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name Expr
expression = SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [String -> Expr
variable String
name, Expr
expression] Either ExprError ()
-> Either ExprError (Polynomial Rational)
-> Either ExprError (Polynomial Rational)
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Expr -> Either ExprError (Polynomial Rational)
go Expr
expression
  where
    go :: Expr -> Either ExprError (Polynomial Rational)
go e :: Expr
e@(Expr SourcePos
_ ExprNode
node) = case ExprNode
node of
        Literal Rational
q -> Polynomial Rational -> Either ExprError (Polynomial Rational)
check ([Rational] -> Polynomial Rational
forall a. (Eq a, Num a) => [a] -> Polynomial a
P.polynomial [Rational
q])
        Variable String
n | String
n String -> String -> Bool
forall a. Eq a => a -> a -> Bool
== String
name -> Polynomial Rational -> Either ExprError (Polynomial Rational)
check ([Rational] -> Polynomial Rational
forall a. (Eq a, Num a) => [a] -> Polynomial a
P.polynomial [Rational
0, Rational
1])
        Unary UnaryOp
Negate Expr
x -> Expr -> Either ExprError (Polynomial Rational)
go Expr
x Either ExprError (Polynomial Rational)
-> (Polynomial Rational -> Either ExprError (Polynomial Rational))
-> Either ExprError (Polynomial Rational)
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Polynomial Rational -> Either ExprError (Polynomial Rational)
check (Polynomial Rational -> Either ExprError (Polynomial Rational))
-> (Polynomial Rational -> Polynomial Rational)
-> Polynomial Rational
-> Either ExprError (Polynomial Rational)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
P.scalePolynomial (-Rational
1)
        Binary BinaryOp
Add Expr
x Expr
y -> (Polynomial Rational -> Polynomial Rational -> Polynomial Rational)
-> Expr -> Expr -> Either ExprError (Polynomial Rational)
combine Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
P.addPolynomial Expr
x Expr
y
        Binary BinaryOp
Subtract Expr
x Expr
y -> (Polynomial Rational -> Polynomial Rational -> Polynomial Rational)
-> Expr -> Expr -> Either ExprError (Polynomial Rational)
combine Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
P.subtractPolynomial Expr
x Expr
y
        Binary BinaryOp
Multiply Expr
x Expr
y -> (Polynomial Rational -> Polynomial Rational -> Polynomial Rational)
-> Expr -> Expr -> Either ExprError (Polynomial Rational)
combine Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
P.multiplyPolynomial Expr
x Expr
y
        Binary BinaryOp
Divide Expr
x (Expr SourcePos
_ (Literal Rational
q)) | Rational
q Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0 -> Expr -> Either ExprError (Polynomial Rational)
go Expr
x Either ExprError (Polynomial Rational)
-> (Polynomial Rational -> Either ExprError (Polynomial Rational))
-> Either ExprError (Polynomial Rational)
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Polynomial Rational -> Either ExprError (Polynomial Rational)
check (Polynomial Rational -> Either ExprError (Polynomial Rational))
-> (Polynomial Rational -> Polynomial Rational)
-> Polynomial Rational
-> Either ExprError (Polynomial Rational)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Num a) => a -> Polynomial a -> Polynomial a
P.scalePolynomial (Rational -> Rational
forall a. Fractional a => a -> a
recip Rational
q)
        Binary BinaryOp
Power Expr
x (Expr SourcePos
_ (Literal Rational
q)) | Rational -> Integer
forall a. Ratio a -> a
denominator Rational
q Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1 Bool -> Bool -> Bool
&& Rational
q Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
1 Bool -> Bool -> Bool
&& Rational
q Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
64 -> do
            p <- Expr -> Either ExprError (Polynomial Rational)
go Expr
x
            foldM (\Polynomial Rational
acc Integer
_ -> Polynomial Rational -> Either ExprError (Polynomial Rational)
check (Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
P.multiplyPolynomial Polynomial Rational
acc Polynomial Rational
p)) (P.polynomial [1]) [1 .. numerator q]
        ExprNode
_ -> Expr -> String -> Either ExprError (Polynomial Rational)
forall a. Expr -> String -> Either ExprError a
unsupported Expr
e String
"expression is not an admitted rational polynomial"
    combine :: (Polynomial Rational -> Polynomial Rational -> Polynomial Rational)
-> Expr -> Expr -> Either ExprError (Polynomial Rational)
combine Polynomial Rational -> Polynomial Rational -> Polynomial Rational
f Expr
x Expr
y = do
        a <- Expr -> Either ExprError (Polynomial Rational)
go Expr
x
        b <- go y
        -- Multiplication cannot allocate more than 129 coefficients.
        check (f a b)
    check :: Polynomial Rational -> Either ExprError (Polynomial Rational)
check Polynomial Rational
p
        | Int -> (Int -> Int) -> Maybe Int -> Int
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Int
0 Int -> Int
forall a. a -> a
id (Polynomial Rational -> Maybe Int
forall a. Polynomial a -> Maybe Int
P.polynomialDegree Polynomial Rational
p) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
64 =
            ErrorKind
-> Expr -> String -> Either ExprError (Polynomial Rational)
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
MemoryLimit Expr
expression String
"polynomial degree limit 64"
        | Bool
otherwise =
            (Rational -> Either ExprError ())
-> [Rational] -> Either ExprError ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Rational
q -> Limits -> Expr -> Either ExprError ()
validateExpression (SymbolicOptions -> Limits
symbolicLimits SymbolicOptions
o) (Rational -> Expr
literal Rational
q)) (Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
P.coefficients Polynomial Rational
p) Either ExprError ()
-> Either ExprError (Polynomial Rational)
-> Either ExprError (Polynomial Rational)
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Polynomial Rational -> Either ExprError (Polynomial Rational)
forall a b. b -> Either a b
Right Polynomial Rational
p

{- | Given variable name then exact Rational Polynomial, build an Expr in Horner form. This
pure builder does not itself validate SymbolicOptions limits.
-}
fromPolynomial :: String -> P.Polynomial Rational -> Expr
fromPolynomial :: String -> Polynomial Rational -> Expr
fromPolynomial String
name Polynomial Rational
p = (Rational -> Expr -> Expr) -> Expr -> [Rational] -> Expr
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\Rational
a Expr
acc -> Expr -> Expr -> Expr
plus (Rational -> Expr
literal Rational
a) (Expr -> Expr -> Expr
times (String -> Expr
variable String
name) Expr
acc)) (Rational -> Expr
literal Rational
0) (Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
P.coefficients Polynomial Rational
p)

{- | Exact leading coefficient and symbolic factor/multiplicity pairs; unfactored remainder
behavior follows the factorization function.
-}
data PolynomialFactorization = PolynomialFactorization
    { PolynomialFactorization -> Rational
factorLeading :: Rational
    -- ^ Exact scalar leading factor outside the polynomial factors.
    , PolynomialFactorization -> [(Expr, Int)]
factorTerms :: [(Expr, Int)]
    -- ^ Symbolic factors paired with positive multiplicities.
    }
    deriving (PolynomialFactorization -> PolynomialFactorization -> Bool
(PolynomialFactorization -> PolynomialFactorization -> Bool)
-> (PolynomialFactorization -> PolynomialFactorization -> Bool)
-> Eq PolynomialFactorization
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: PolynomialFactorization -> PolynomialFactorization -> Bool
== :: PolynomialFactorization -> PolynomialFactorization -> Bool
$c/= :: PolynomialFactorization -> PolynomialFactorization -> Bool
/= :: PolynomialFactorization -> PolynomialFactorization -> Bool
Eq, Int -> PolynomialFactorization -> ShowS
[PolynomialFactorization] -> ShowS
PolynomialFactorization -> String
(Int -> PolynomialFactorization -> ShowS)
-> (PolynomialFactorization -> String)
-> ([PolynomialFactorization] -> ShowS)
-> Show PolynomialFactorization
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> PolynomialFactorization -> ShowS
showsPrec :: Int -> PolynomialFactorization -> ShowS
$cshow :: PolynomialFactorization -> String
show :: PolynomialFactorization -> String
$cshowList :: [PolynomialFactorization] -> ShowS
showList :: [PolynomialFactorization] -> ShowS
Show)

-- | Exact factorization over Q via the existing bounded Kronecker algorithm.
factorPolynomial :: SymbolicOptions -> String -> Expr -> Either ExprError PolynomialFactorization
factorPolynomial :: SymbolicOptions
-> String -> Expr -> Either ExprError PolynomialFactorization
factorPolynomial SymbolicOptions
o String
name Expr
e = do
    p <- SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name Expr
e
    (c, fs) <-
        mathResult e (P.factorRationalPolynomial (min 100000 (maximumSteps (symbolicLimits o))) p)
    terms <- mapM (\(Polynomial Rational
q, Int
k) -> do x <- SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o (String -> Polynomial Rational -> Expr
fromPolynomial String
name Polynomial Rational
q); Right (x, k)) fs
    Right (PolynomialFactorization c terms)

mathResult :: (Show e) => Expr -> Either e a -> Either ExprError a
mathResult :: forall e a. Show e => Expr -> Either e a -> Either ExprError a
mathResult Expr
e = (e -> Either ExprError a)
-> (a -> Either ExprError a) -> Either e a -> Either ExprError a
forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either (Expr -> String -> Either ExprError a
forall a. Expr -> String -> Either ExprError a
unsupported Expr
e (String -> Either ExprError a)
-> (e -> String) -> e -> Either ExprError a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. e -> String
forall a. Show a => a -> String
show) a -> Either ExprError a
forall a b. b -> Either a b
Right

{- | Cancel a rational quotient only when its original denominator is proved
nonzero; this preserves holes that ordinary field cancellation would erase.
-}
rationalSimplify :: SymbolicOptions -> String -> Expr -> Either ExprError Expr
rationalSimplify :: SymbolicOptions -> String -> Expr -> Either ExprError Expr
rationalSimplify SymbolicOptions
o String
name Expr
e = case Expr -> ExprNode
expressionNode Expr
e of
    Binary BinaryOp
Divide Expr
n Expr
d | Assumptions -> Property -> Expr -> Bool
proves (SymbolicOptions -> Assumptions
symbolicAssumptions SymbolicOptions
o) Property
NonZero Expr
d -> do
        p <- SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name Expr
n
        q <- toPolynomial o name d
        r <- mathResult e (P.rationalFunction p q)
        simplify
            o
            (divide (fromPolynomial name (P.rationalNumerator r)) (fromPolynomial name (P.rationalDenominator r)))
    ExprNode
_ -> SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o Expr
e

{- | Given options then Expr, apply supported trigonometric identities such as matching
sin-squared plus cos-squared terms, after input validation. This is a bounded rule set,
not a general identity prover.
-}
trigonometricSimplify :: SymbolicOptions -> Expr -> Either ExprError Expr
trigonometricSimplify :: SymbolicOptions -> Expr -> Either ExprError Expr
trigonometricSimplify SymbolicOptions
o Expr
e = do
    SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [Expr
e]
    let squared :: String -> Expr -> Maybe Expr
squared String
name Expr
x = case Expr -> ExprNode
expressionNode Expr
x of
            Binary BinaryOp
Power (Expr SourcePos
_ (Call String
n [Expr
a])) (Expr SourcePos
_ (Literal Rational
2)) | String
name String -> String -> Bool
forall a. Eq a => a -> a -> Bool
== String
n -> Expr -> Maybe Expr
forall a. a -> Maybe a
Just Expr
a
            ExprNode
_ -> Maybe Expr
forall a. Maybe a
Nothing
        identity :: Expr -> Expr -> Bool
identity Expr
x Expr
y = case (String -> Expr -> Maybe Expr
squared String
"sin" Expr
x, String -> Expr -> Maybe Expr
squared String
"cos" Expr
y) of
            (Just Expr
a, Just Expr
b) -> Expr
a Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Expr
b Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined (SymbolicOptions -> Assumptions
symbolicAssumptions SymbolicOptions
o) Expr
a
            (Maybe Expr, Maybe Expr)
_ -> Bool
False
    case Expr -> ExprNode
expressionNode Expr
e of
        Binary BinaryOp
Add Expr
a Expr
b | Expr -> Expr -> Bool
identity Expr
a Expr
b Bool -> Bool -> Bool
|| Expr -> Expr -> Bool
identity Expr
b Expr
a -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Rational -> Expr
literal Rational
1)
        ExprNode
_ -> SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o Expr
e

{- | Given options then Expr, apply log product/quotient/power rules only when the supplied
assumptions establish required real-domain facts. Unproved transformations are left
unchanged; budget errors propagate.
-}
logarithmicSimplify :: SymbolicOptions -> Expr -> Either ExprError Expr
logarithmicSimplify :: SymbolicOptions -> Expr -> Either ExprError Expr
logarithmicSimplify SymbolicOptions
o Expr
e = do
    SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [Expr
e]
    case Expr -> ExprNode
expressionNode Expr
e of
        Call String
"log" [Expr SourcePos
_ (Binary BinaryOp
Multiply Expr
a Expr
b)]
            | (Expr -> Bool) -> [Expr] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Assumptions -> Property -> Expr -> Bool
proves (SymbolicOptions -> Assumptions
symbolicAssumptions SymbolicOptions
o) Property
Positive) [Expr
a, Expr
b] ->
                SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o (Expr -> Expr -> Expr
plus (String -> [Expr] -> Expr
call String
"log" [Expr
a]) (String -> [Expr] -> Expr
call String
"log" [Expr
b]))
        ExprNode
_ -> SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o Expr
e

{- | Polynomial primitives and affine exp/sin/cos and reciprocal powers. The
logarithmic primitive log(abs x) has the original nonzero domain.
-}
integrateBasic :: SymbolicOptions -> String -> Expr -> Either ExprError Expr
integrateBasic :: SymbolicOptions -> String -> Expr -> Either ExprError Expr
integrateBasic SymbolicOptions
o String
name Expr
expression = do
    SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [String -> Expr
variable String
name, Expr
expression]
    Expr -> Either ExprError Expr
primitive Expr
expression Either ExprError Expr
-> (Expr -> Either ExprError Expr) -> Either ExprError Expr
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o
  where
    x :: Expr
x = String -> Expr
variable String
name
    primitive :: Expr -> Either ExprError Expr
primitive Expr
e = case SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name Expr
e of
        Right Polynomial Rational
p -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (String -> Polynomial Rational -> Expr
fromPolynomial String
name (Rational -> Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Fractional a) => a -> Polynomial a -> Polynomial a
P.integratePolynomial Rational
0 Polynomial Rational
p))
        Left ExprError
_ -> case Expr -> ExprNode
expressionNode Expr
e of
            Binary BinaryOp
Add Expr
a Expr
b -> Expr -> Expr -> Expr
plus (Expr -> Expr -> Expr)
-> Either ExprError Expr -> Either ExprError (Expr -> Expr)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Expr -> Either ExprError Expr
primitive Expr
a Either ExprError (Expr -> Expr)
-> Either ExprError Expr -> Either ExprError Expr
forall a b.
Either ExprError (a -> b)
-> Either ExprError a -> Either ExprError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Expr -> Either ExprError Expr
primitive Expr
b
            Binary BinaryOp
Subtract Expr
a Expr
b -> Expr -> Expr -> Expr
minus (Expr -> Expr -> Expr)
-> Either ExprError Expr -> Either ExprError (Expr -> Expr)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Expr -> Either ExprError Expr
primitive Expr
a Either ExprError (Expr -> Expr)
-> Either ExprError Expr -> Either ExprError Expr
forall a b.
Either ExprError (a -> b)
-> Either ExprError a -> Either ExprError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Expr -> Either ExprError Expr
primitive Expr
b
            Unary UnaryOp
Negate Expr
a -> Expr -> Expr
neg (Expr -> Expr) -> Either ExprError Expr -> Either ExprError Expr
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Expr -> Either ExprError Expr
primitive Expr
a
            Binary BinaryOp
Multiply Expr
a Expr
b | String
name String -> [String] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` Expr -> [String]
freeVariables Expr
a -> Expr -> Expr -> Expr
times Expr
a (Expr -> Expr) -> Either ExprError Expr -> Either ExprError Expr
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Expr -> Either ExprError Expr
primitive Expr
b
            Binary BinaryOp
Multiply Expr
a Expr
b | String
name String -> [String] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` Expr -> [String]
freeVariables Expr
b -> (Expr -> Expr -> Expr
`times` Expr
b) (Expr -> Expr) -> Either ExprError Expr -> Either ExprError Expr
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Expr -> Either ExprError Expr
primitive Expr
a
            Binary BinaryOp
Divide (Expr SourcePos
_ (Literal Rational
c)) Expr
a -> Rational -> Expr -> Either ExprError Expr
reciprocalPrimitive Rational
c Expr
a
            Binary BinaryOp
Power Expr
a (Expr SourcePos
_ (Literal Rational
q)) -> do
                slope <- Expr -> Either ExprError Rational
affineSlope Expr
a
                if q == (-1)
                    then Right (divide (call "log" [call "abs" [a]]) (literal slope))
                    else
                        if denominator q == 1 || proves (symbolicAssumptions o) Positive a
                            then Right (divide (power a (literal (q + 1))) (literal (slope * (q + 1))))
                            else unsupported e "fractional power primitive needs positive base"
            Call String
f [Expr
a] | String
f String -> [String] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [String
"exp", String
"sin", String
"cos"] -> do
                slope <- Expr -> Either ExprError Rational
affineSlope Expr
a
                Right
                    (divide (case f of String
"exp" -> Expr
e; String
"sin" -> Expr -> Expr
neg (String -> [Expr] -> Expr
call String
"cos" [Expr
a]); String
_ -> String -> [Expr] -> Expr
call String
"sin" [Expr
a]) (literal slope))
            ExprNode
_
                | String
name String -> [String] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` Expr -> [String]
freeVariables Expr
e Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined (SymbolicOptions -> Assumptions
symbolicAssumptions SymbolicOptions
o) Expr
e -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Expr -> Expr -> Expr
times Expr
e Expr
x)
                | Bool
otherwise -> Expr -> String -> Either ExprError Expr
forall a. Expr -> String -> Either ExprError a
unsupported Expr
e String
"no supported elementary primitive"
    reciprocalPrimitive :: Rational -> Expr -> Either ExprError Expr
reciprocalPrimitive Rational
c Expr
a = do slope <- Expr -> Either ExprError Rational
affineSlope Expr
a; Right (times (literal (c / slope)) (call "log" [call "abs" [a]]))
    affineSlope :: Expr -> Either ExprError Rational
affineSlope Expr
a = do
        p <- SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name Expr
a
        case P.coefficients p of
            [Rational
_, Rational
m] | Rational
m Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0 -> Rational -> Either ExprError Rational
forall a b. b -> Either a b
Right Rational
m
            [Rational]
_ -> Expr -> String -> Either ExprError Rational
forall a. Expr -> String -> Either ExprError a
unsupported Expr
a String
"primitive requires nonconstant affine argument"

{- | Integrate a supplied u*dv by recursively reducing polynomial degree when
v is an affine exponential/sine/cosine primitive. Degree is capped at 16.
-}
integrateByParts :: SymbolicOptions -> String -> Expr -> Expr -> Either ExprError Expr
integrateByParts :: SymbolicOptions -> String -> Expr -> Expr -> Either ExprError Expr
integrateByParts SymbolicOptions
o String
name Expr
u Expr
dv = do
    p <- SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name Expr
u
    if maybe 0 id (P.polynomialDegree p) > 16
        then failure MemoryLimit u "integration-by-parts degree limit 16"
        else parts (0 :: Int) u dv >>= simplify o
  where
    parts :: t -> Expr -> Expr -> Either ExprError Expr
parts t
n Expr
a Expr
b
        | t
n t -> t -> Bool
forall a. Ord a => a -> a -> Bool
> t
20 = ErrorKind -> Expr -> String -> Either ExprError Expr
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
StepLimit Expr
a String
"integration-by-parts recursion limit"
        | Bool
otherwise = do
            v <- SymbolicOptions -> String -> Expr -> Either ExprError Expr
integrateBasic SymbolicOptions
o String
name Expr
b
            du <- differentiate o name a
            if qvalue du == Just 0
                then checked o (times a v)
                else do
                    residual <- case integrateBasic o name (times du v) of
                        Right Expr
r -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right Expr
r
                        Left ExprError
_ -> case Expr -> ExprNode
expressionNode Expr
v of
                            Binary BinaryOp
Divide Expr
g (Expr SourcePos
_ (Literal Rational
c)) -> Expr -> Expr -> Expr
divide (Expr -> Expr -> Expr)
-> Either ExprError Expr -> Either ExprError (Expr -> Expr)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> t -> Expr -> Expr -> Either ExprError Expr
parts (t
n t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) Expr
du Expr
g Either ExprError (Expr -> Expr)
-> Either ExprError Expr -> Either ExprError Expr
forall a b.
Either ExprError (a -> b)
-> Either ExprError a -> Either ExprError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Expr -> Either ExprError Expr
forall a. a -> Either ExprError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Rational -> Expr
literal Rational
c)
                            Unary UnaryOp
Negate Expr
g -> Expr -> Expr
neg (Expr -> Expr) -> Either ExprError Expr -> Either ExprError Expr
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> t -> Expr -> Expr -> Either ExprError Expr
parts (t
n t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) Expr
du Expr
g
                            ExprNode
_ -> t -> Expr -> Expr -> Either ExprError Expr
parts (t
n t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) Expr
du Expr
v
                    checked o (minus (times a v) residual)

{- | Verified substitution: callers supply an outer expression in a fresh
variable t; the integrand must simplify to outer(inner)*inner'.
-}
integrateSubstitution ::
    SymbolicOptions -> String -> String -> Expr -> Expr -> Expr -> Either ExprError Expr
integrateSubstitution :: SymbolicOptions
-> String
-> String
-> Expr
-> Expr
-> Expr
-> Either ExprError Expr
integrateSubstitution SymbolicOptions
o String
name String
temporary Expr
inner Expr
outer Expr
integrand = do
    if String
name String -> String -> Bool
forall a. Eq a => a -> a -> Bool
== String
temporary Bool -> Bool -> Bool
|| String
temporary String -> [String] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` Expr -> [String]
freeVariables Expr
inner
        then ErrorKind -> Expr -> String -> Either ExprError ()
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
DomainFailure Expr
inner String
"substitution variable must be fresh"
        else () -> Either ExprError ()
forall a b. b -> Either a b
Right ()
    d <- SymbolicOptions -> String -> Expr -> Either ExprError Expr
differentiate SymbolicOptions
o String
name Expr
inner
    composed <- substitute (symbolicLimits o) [(temporary, inner)] outer
    difference <- expand o (minus integrand (times composed d)) >>= simplify o
    if qvalue difference /= Just 0
        then unsupported integrand "integrand does not match the supplied substitution"
        else do
            primitive <- integrateBasic o temporary outer
            substitute (symbolicLimits o) [(temporary, inner)] primitive >>= simplify o

{- | Exact partial fractions; supports all repeated linear poles and simple
irreducible quadratic factors. Unsupported higher quadratic powers are explicit.
-}
integrateRational :: SymbolicOptions -> String -> Expr -> Expr -> Either ExprError Expr
integrateRational :: SymbolicOptions -> String -> Expr -> Expr -> Either ExprError Expr
integrateRational SymbolicOptions
o String
name Expr
numeratorExpr Expr
denominatorExpr = do
    p <- SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name Expr
numeratorExpr
    q <- toPolynomial o name denominatorExpr
    rf <- mathResult denominatorExpr (P.rationalFunction p q)
    (whole, fs) <-
        mathResult denominatorExpr (P.partialFractions (min 100000 (maximumSteps (symbolicLimits o))) rf)
    terms <- mapM term fs
    simplify o (Foldable.foldl' plus (fromPolynomial name (P.integratePolynomial 0 whole)) terms)
  where
    x :: Expr
x = String -> Expr
variable String
name
    term :: (Polynomial Rational, Polynomial Rational, a)
-> Either ExprError Expr
term (Polynomial Rational
n, Polynomial Rational
d, a
k) = case (Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
P.coefficients Polynomial Rational
n, Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
P.coefficients Polynomial Rational
d) of
        ([Rational
c], [Rational
b, Rational
1]) ->
            let linear :: Expr
linear = Expr -> Expr -> Expr
plus Expr
x (Rational -> Expr
literal Rational
b)
             in Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right
                    ( if a
k a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
1
                        then Expr -> Expr -> Expr
times (Rational -> Expr
literal Rational
c) (String -> [Expr] -> Expr
call String
"log" [String -> [Expr] -> Expr
call String
"abs" [Expr
linear]])
                        else Expr -> Expr -> Expr
times (Rational -> Expr
literal (Rational
c Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ a -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a
k))) (Expr -> Expr -> Expr
power Expr
linear (Rational -> Expr
literal (a -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a
k))))
                    )
        ([Rational]
coeff, [Rational
c, Rational
b, Rational
1]) | a
k a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
1 Bool -> Bool -> Bool
&& Rational
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
b Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Rational
4 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
c -> do
            let a :: Rational
a = case [Rational]
coeff of Rational
_ : Rational
z : [Rational]
_ -> Rational
z; [Rational]
_ -> Rational
0
                b0 :: Rational
b0 = case [Rational]
coeff of Rational
z : [Rational]
_ -> Rational
z; [Rational]
_ -> Rational
0
                disc :: Expr
disc = Rational -> Expr
literal (Rational
4 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
c Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
b)
                root :: Expr
root = String -> [Expr] -> Expr
call String
"sqrt" [Expr
disc]
                quadratic :: Expr
quadratic = String -> Polynomial Rational -> Expr
fromPolynomial String
name Polynomial Rational
d
                logTerm :: Expr
logTerm = Expr -> Expr -> Expr
times (Rational -> Expr
literal (Rational
a Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
2)) (String -> [Expr] -> Expr
call String
"log" [Expr
quadratic])
                atanTerm :: Expr
atanTerm =
                    Expr -> Expr -> Expr
times
                        (Expr -> Expr -> Expr
divide (Rational -> Expr
literal (Rational
2 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
b0 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
b)) Expr
root)
                        (String -> [Expr] -> Expr
call String
"atan" [Expr -> Expr -> Expr
divide (Expr -> Expr -> Expr
plus (Expr -> Expr -> Expr
times (Rational -> Expr
literal Rational
2) Expr
x) (Rational -> Expr
literal Rational
b)) Expr
root])
            SymbolicOptions -> Expr -> Either ExprError Expr
checked SymbolicOptions
o (Expr -> Expr -> Expr
plus Expr
logTerm Expr
atanTerm)
        ([Rational], [Rational])
_ ->
            Expr -> String -> Either ExprError Expr
forall a. Expr -> String -> Either ExprError a
unsupported
                Expr
denominatorExpr
                String
"partial fraction primitive requires linear or simple positive quadratic denominator"

-- | Limit approach direction: from left,from right,or two-sided.
data Direction = FromLeft | FromRight | TwoSided deriving (Direction -> Direction -> Bool
(Direction -> Direction -> Bool)
-> (Direction -> Direction -> Bool) -> Eq Direction
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Direction -> Direction -> Bool
== :: Direction -> Direction -> Bool
$c/= :: Direction -> Direction -> Bool
/= :: Direction -> Direction -> Bool
Eq, Int -> Direction -> ShowS
[Direction] -> ShowS
Direction -> String
(Int -> Direction -> ShowS)
-> (Direction -> String)
-> ([Direction] -> ShowS)
-> Show Direction
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Direction -> ShowS
showsPrec :: Int -> Direction -> ShowS
$cshow :: Direction -> String
show :: Direction -> String
$cshowList :: [Direction] -> ShowS
showList :: [Direction] -> ShowS
Show)

{- | Symbolic finite/infinite limit classification,including an explicit does-not-exist or
unresolved result.
-}
data LimitResult = FiniteLimit Rational | PositiveInfinity | NegativeInfinity | DoesNotExist
    deriving (LimitResult -> LimitResult -> Bool
(LimitResult -> LimitResult -> Bool)
-> (LimitResult -> LimitResult -> Bool) -> Eq LimitResult
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: LimitResult -> LimitResult -> Bool
== :: LimitResult -> LimitResult -> Bool
$c/= :: LimitResult -> LimitResult -> Bool
/= :: LimitResult -> LimitResult -> Bool
Eq, Int -> LimitResult -> ShowS
[LimitResult] -> ShowS
LimitResult -> String
(Int -> LimitResult -> ShowS)
-> (LimitResult -> String)
-> ([LimitResult] -> ShowS)
-> Show LimitResult
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> LimitResult -> ShowS
showsPrec :: Int -> LimitResult -> ShowS
$cshow :: LimitResult -> String
show :: LimitResult -> String
$cshowList :: [LimitResult] -> ShowS
showList :: [LimitResult] -> ShowS
Show)

{- | Exact finite-point rational limits by root multiplicity. Non-rational
transcendental limits are explicitly unsupported.
-}
symbolicLimit ::
    SymbolicOptions -> String -> Rational -> Direction -> Expr -> Either ExprError LimitResult
symbolicLimit :: SymbolicOptions
-> String
-> Rational
-> Direction
-> Expr
-> Either ExprError LimitResult
symbolicLimit SymbolicOptions
o String
name Rational
center Direction
direction Expr
e = do
    SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [Expr
e, Rational -> Expr
literal Rational
center]
    (n, d) <- SymbolicOptions
-> String
-> Expr
-> Either ExprError (Polynomial Rational, Polynomial Rational)
rationalPolynomials SymbolicOptions
o String
name Expr
e
    let (kn, an) = leadingAt center n; (kd, ad) = leadingAt center d
    if ad == 0
        then failure DomainFailure e "identically zero denominator"
        else
            if an == 0 || kn > kd
                then finite 0
                else
                    if kn == kd
                        then finite (an / ad)
                        else
                            let pole = Int
kd Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
kn
                                signRight = Rational -> Rational
forall a. Num a => a -> a
signum (Rational
an Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
ad)
                                signed a
s = if a
s a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 then LimitResult
PositiveInfinity else LimitResult
NegativeInfinity
                             in Right
                                    ( case direction of
                                        Direction
FromRight -> Rational -> LimitResult
forall {a}. (Ord a, Num a) => a -> LimitResult
signed Rational
signRight
                                        Direction
FromLeft -> Rational -> LimitResult
forall {a}. (Ord a, Num a) => a -> LimitResult
signed (if Int -> Bool
forall a. Integral a => a -> Bool
odd Int
pole then -Rational
signRight else Rational
signRight)
                                        Direction
TwoSided -> if Int -> Bool
forall a. Integral a => a -> Bool
odd Int
pole then LimitResult
DoesNotExist else Rational -> LimitResult
forall {a}. (Ord a, Num a) => a -> LimitResult
signed Rational
signRight
                                    )
  where
    finite :: Rational -> Either ExprError LimitResult
finite Rational
q = SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [Rational -> Expr
literal Rational
q] Either ExprError ()
-> Either ExprError LimitResult -> Either ExprError LimitResult
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> LimitResult -> Either ExprError LimitResult
forall a b. b -> Either a b
Right (Rational -> LimitResult
FiniteLimit Rational
q)

rationalPolynomials ::
    SymbolicOptions -> String -> Expr -> Either ExprError (P.Polynomial Rational, P.Polynomial Rational)
rationalPolynomials :: SymbolicOptions
-> String
-> Expr
-> Either ExprError (Polynomial Rational, Polynomial Rational)
rationalPolynomials SymbolicOptions
o String
name Expr
e = case Expr -> ExprNode
expressionNode Expr
e of
    Binary BinaryOp
Divide Expr
n Expr
d -> (,) (Polynomial Rational
 -> Polynomial Rational
 -> (Polynomial Rational, Polynomial Rational))
-> Either ExprError (Polynomial Rational)
-> Either
     ExprError
     (Polynomial Rational -> (Polynomial Rational, Polynomial Rational))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name Expr
n Either
  ExprError
  (Polynomial Rational -> (Polynomial Rational, Polynomial Rational))
-> Either ExprError (Polynomial Rational)
-> Either ExprError (Polynomial Rational, Polynomial Rational)
forall a b.
Either ExprError (a -> b)
-> Either ExprError a -> Either ExprError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name Expr
d
    ExprNode
_ -> do p <- SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name Expr
e; Right (p, P.polynomial [1])
leadingAt :: Rational -> P.Polynomial Rational -> (Int, Rational)
leadingAt :: Rational -> Polynomial Rational -> (Int, Rational)
leadingAt Rational
center Polynomial Rational
p = Int -> Polynomial Rational -> (Int, Rational)
forall {t}. Integral t => t -> Polynomial Rational -> (t, Rational)
go Int
0 Polynomial Rational
p
  where
    go :: t -> Polynomial Rational -> (t, Rational)
go t
k Polynomial Rational
q
        | [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
q) = (t
k, Rational
0)
        | Rational
value Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0 = (t
k, Rational
value)
        | Bool
otherwise = t -> Polynomial Rational -> (t, Rational)
go (t
k t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) (Polynomial Rational -> Polynomial Rational
forall a. (Eq a, Num a) => Polynomial a -> Polynomial a
P.differentiatePolynomial Polynomial Rational
q)
      where
        value :: Rational
value = Polynomial Rational -> Rational -> Rational
forall a. Num a => Polynomial a -> a -> a
P.evaluatePolynomial Polynomial Rational
q Rational
center Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Integer -> Rational
forall a. Num a => Integer -> a
fromInteger ([Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
1 .. t -> Integer
forall a. Integral a => a -> Integer
toInteger t
k])

{- | Finite symbolic Taylor/Laurent coefficient list with exact center and starting power;
no analytic convergence is implied by storage alone.
-}
data Series = Series
    { Series -> Rational
seriesCenter :: Rational
    -- ^ Exact Rational expansion center.
    , Series -> Int
seriesMinimumPower :: Int
    -- ^ Power corresponding to the first coefficient, including negative powers for Laurent series.
    , Series -> [Expr]
seriesTerms :: [Expr]
    -- ^ Coefficient expressions ordered from seriesMinimumPower upward.
    }
    deriving (Series -> Series -> Bool
(Series -> Series -> Bool)
-> (Series -> Series -> Bool) -> Eq Series
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Series -> Series -> Bool
== :: Series -> Series -> Bool
$c/= :: Series -> Series -> Bool
/= :: Series -> Series -> Bool
Eq, Int -> Series -> ShowS
[Series] -> ShowS
Series -> String
(Int -> Series -> ShowS)
-> (Series -> String) -> ([Series] -> ShowS) -> Show Series
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Series -> ShowS
showsPrec :: Int -> Series -> ShowS
$cshow :: Series -> String
show :: Series -> String
$cshowList :: [Series] -> ShowS
showList :: [Series] -> ShowS
Show)

{- | Truncated Taylor coefficients f^(k)(c)/k!, order 0..20. Missing regularity
assumptions or a pole produce an error rather than a fictitious Taylor jet.
-}
taylorSeries :: SymbolicOptions -> String -> Rational -> Int -> Expr -> Either ExprError Series
taylorSeries :: SymbolicOptions
-> String -> Rational -> Int -> Expr -> Either ExprError Series
taylorSeries SymbolicOptions
o String
name Rational
center Int
order Expr
e = do
    Expr -> Int -> Either ExprError ()
seriesOrder Expr
e Int
order
    SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [Expr
e, Rational -> Expr
literal Rational
center]
    coefficients <- Int -> Integer -> Expr -> Either ExprError [Expr]
go Int
0 Integer
1 Expr
e
    Right (Series center 0 coefficients)
  where
    go :: Int -> Integer -> Expr -> Either ExprError [Expr]
go Int
k Integer
factorial Expr
current = do
        value <- Limits -> [(String, Expr)] -> Expr -> Either ExprError Expr
substitute (SymbolicOptions -> Limits
symbolicLimits SymbolicOptions
o) [(String
name, Rational -> Expr
literal Rational
center)] Expr
current Either ExprError Expr
-> (Expr -> Either ExprError Expr) -> Either ExprError Expr
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o
        if not (defined noAssumptions value)
            then unsupported value "Taylor coefficient is singular or domain unproved"
            else Right ()
        c <- simplify o (divide value (literal (fromInteger factorial)))
        if k == order
            then Right [c]
            else do
                next <- differentiate o name current
                more <- go (k + 1) (factorial * toInteger (k + 1)) next
                Right (c : more)

{- | Given options,variable name,order,then Expr, compute the Taylor series centered at zero
via taylorSeries. Unsupported symbolic derivatives or exceeded order/resource bounds
return ExprError.
-}
maclaurinSeries :: SymbolicOptions -> String -> Int -> Expr -> Either ExprError Series
maclaurinSeries :: SymbolicOptions -> String -> Int -> Expr -> Either ExprError Series
maclaurinSeries SymbolicOptions
o String
name = SymbolicOptions
-> String -> Rational -> Int -> Expr -> Either ExprError Series
taylorSeries SymbolicOptions
o String
name Rational
0

{- | Rational Laurent coefficients at a rational center, from the lowest power
through that power+order. Exact finite jets reuse HMath.Polynomial division.
-}
laurentSeries :: SymbolicOptions -> String -> Rational -> Int -> Expr -> Either ExprError Series
laurentSeries :: SymbolicOptions
-> String -> Rational -> Int -> Expr -> Either ExprError Series
laurentSeries SymbolicOptions
o String
name Rational
center Int
order Expr
e = do
    Expr -> Int -> Either ExprError ()
seriesOrder Expr
e Int
order
    SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [Expr
e, Rational -> Expr
literal Rational
center]
    (p, q) <- SymbolicOptions
-> String
-> Expr
-> Either ExprError (Polynomial Rational, Polynomial Rational)
rationalPolynomials SymbolicOptions
o String
name Expr
e
    let shift = [Rational] -> Polynomial Rational
forall a. (Eq a, Num a) => [a] -> Polynomial a
P.polynomial [Rational
center, Rational
1]
        as = Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
P.coefficients (Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
P.composePolynomial Polynomial Rational
p Polynomial Rational
shift)
        bs = Polynomial Rational -> [Rational]
forall a. Polynomial a -> [a]
P.coefficients (Polynomial Rational -> Polynomial Rational -> Polynomial Rational
forall a.
(Eq a, Num a) =>
Polynomial a -> Polynomial a -> Polynomial a
P.composePolynomial Polynomial Rational
q Polynomial Rational
shift)
        na = [Rational] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Rational -> Bool) -> [Rational] -> [Rational]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0) [Rational]
as)
        nb = [Rational] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Rational -> Bool) -> [Rational] -> [Rational]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0) [Rational]
bs)
        pad [a]
xs = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
order Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) ([a]
xs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
0)
    if null bs
        then failure DomainFailure e "zero Laurent denominator"
        else
            if null as
                then Right (Series center 0 (replicate (order + 1) (literal 0)))
                else do
                    jet <-
                        mathResult e (P.divideSeries (P.powerSeries (pad (drop na as))) (P.powerSeries (pad (drop nb bs))))
                    let coefficients = Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
take (Int
order Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (PowerSeries Rational -> [Rational]
forall a. PowerSeries a -> [a]
P.seriesCoefficients PowerSeries Rational
jet)
                    inputs o (map literal coefficients)
                    Right (Series center (na - nb) (map literal coefficients))

seriesOrder :: Expr -> Int -> Either ExprError ()
seriesOrder :: Expr -> Int -> Either ExprError ()
seriesOrder Expr
e Int
n
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
20 = ErrorKind -> Expr -> String -> Either ExprError ()
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
MemoryLimit Expr
e String
"series order must be 0..20"
    | Bool
otherwise = () -> Either ExprError ()
forall a b. b -> Either a b
Right ()

{- | Given options then Rational coefficient list, construct a finite formal jet with at
most 64 coefficients, subject to expression numeric budgets. No analytic convergence is
claimed.
-}
formalPowerSeries :: SymbolicOptions -> [Rational] -> Either ExprError (P.PowerSeries Rational)
formalPowerSeries :: SymbolicOptions
-> [Rational] -> Either ExprError (PowerSeries Rational)
formalPowerSeries SymbolicOptions
o [Rational]
coefficients
    | [Rational] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
take Int
65 [Rational]
coefficients) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
64 =
        ErrorKind
-> Expr -> String -> Either ExprError (PowerSeries Rational)
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
MemoryLimit (Rational -> Expr
literal Rational
0) String
"formal jet supports at most 64 coefficients"
    | Bool
otherwise = SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o ((Rational -> Expr) -> [Rational] -> [Expr]
forall a b. (a -> b) -> [a] -> [b]
map Rational -> Expr
literal [Rational]
coefficients) Either ExprError ()
-> Either ExprError (PowerSeries Rational)
-> Either ExprError (PowerSeries Rational)
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> PowerSeries Rational -> Either ExprError (PowerSeries Rational)
forall a b. b -> Either a b
Right ([Rational] -> PowerSeries Rational
forall a. [a] -> PowerSeries a
P.powerSeries [Rational]
coefficients)

{- | Given options then two formal Rational jets, validate and convolve coefficients under
the 64-coefficient cap. Oversized input/output or numeric-limit failures return
ExprError.
-}
multiplyFormalSeries ::
    SymbolicOptions ->
    P.PowerSeries Rational ->
    P.PowerSeries Rational ->
    Either ExprError (P.PowerSeries Rational)
multiplyFormalSeries :: SymbolicOptions
-> PowerSeries Rational
-> PowerSeries Rational
-> Either ExprError (PowerSeries Rational)
multiplyFormalSeries SymbolicOptions
o PowerSeries Rational
a PowerSeries Rational
b = do
    _ <- SymbolicOptions
-> [Rational] -> Either ExprError (PowerSeries Rational)
formalPowerSeries SymbolicOptions
o (PowerSeries Rational -> [Rational]
forall a. PowerSeries a -> [a]
P.seriesCoefficients PowerSeries Rational
a)
    _ <- formalPowerSeries o (P.seriesCoefficients b)
    formalPowerSeries o (P.seriesCoefficients (P.multiplySeries a b))

{- | Given options then formal Rational series with nonzero constant term, compute its
reciprocal jet. Invalid constant/size or numeric budgets return ExprError.
-}
reciprocalFormalSeries ::
    SymbolicOptions -> P.PowerSeries Rational -> Either ExprError (P.PowerSeries Rational)
reciprocalFormalSeries :: SymbolicOptions
-> PowerSeries Rational -> Either ExprError (PowerSeries Rational)
reciprocalFormalSeries SymbolicOptions
o PowerSeries Rational
a = do
    _ <- SymbolicOptions
-> [Rational] -> Either ExprError (PowerSeries Rational)
formalPowerSeries SymbolicOptions
o (PowerSeries Rational -> [Rational]
forall a. PowerSeries a -> [a]
P.seriesCoefficients PowerSeries Rational
a)
    result <- mathResult (literal 0) (P.reciprocalSeries a)
    formalPowerSeries o (P.seriesCoefficients result)

{- | Given options,variable name,then polynomial Expr, return rational roots with
multiplicities plus the remaining deflated Expr. Irrational/non-rational roots remain
in the remainder; unsupported conversion or work-budget exhaustion returns ExprError.
-}
exactPolynomialRoots ::
    SymbolicOptions -> String -> Expr -> Either ExprError ([(Rational, Int)], Expr)
exactPolynomialRoots :: SymbolicOptions
-> String -> Expr -> Either ExprError ([(Rational, Int)], Expr)
exactPolynomialRoots SymbolicOptions
o String
name Expr
e = do
    p <- SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name Expr
e
    roots <- mathResult e (P.rationalRoots (min 100000 (maximumSteps (symbolicLimits o))) p)
    let deflate Polynomial a
q (a
r, a
k) =
            (Polynomial a -> a -> Polynomial a)
-> Polynomial a -> [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
Foldable.foldl'
                ( \Polynomial a
current 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)
P.dividePolynomial Polynomial a
current ([a] -> Polynomial a
forall a. (Eq a, Num a) => [a] -> Polynomial a
P.polynomial [-a
r, a
1]) of
                    Right (Polynomial a
result, Polynomial a
_) -> Polynomial a
result
                    Left MathError
_ -> Polynomial a
current
                )
                Polynomial a
q
                [a
1 .. a
k]
    residual <- simplify o (fromPolynomial name (Foldable.foldl' deflate p roots))
    Right (roots, residual)

-- | Bounded Durand-Kerner roots, degree <=32; nonconvergence is explicit.
numericalPolynomialRoots :: SymbolicOptions -> String -> Expr -> Either ExprError [Complex Double]
numericalPolynomialRoots :: SymbolicOptions
-> String -> Expr -> Either ExprError [Complex Double]
numericalPolynomialRoots SymbolicOptions
o String
name Expr
e = do
    p <- SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name Expr
e
    if maybe 0 id (P.polynomialDegree p) > 32
        then failure MemoryLimit e "numerical root degree limit 32"
        else
            mathResult
                e
                ( P.durandKerner
                    (min 2000 (maximumSteps (symbolicLimits o)))
                    1e-10
                    (P.polynomial [fromRational q :+ 0 | q <- P.coefficients p])
                )

-- | Comparison relation selector: equality,inequality,or ordered less/greater comparisons.
data RelationOp = Equal | NotEqual | Less | LessEqual | Greater | GreaterEqual deriving (RelationOp -> RelationOp -> Bool
(RelationOp -> RelationOp -> Bool)
-> (RelationOp -> RelationOp -> Bool) -> Eq RelationOp
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RelationOp -> RelationOp -> Bool
== :: RelationOp -> RelationOp -> Bool
$c/= :: RelationOp -> RelationOp -> Bool
/= :: RelationOp -> RelationOp -> Bool
Eq, Int -> RelationOp -> ShowS
[RelationOp] -> ShowS
RelationOp -> String
(Int -> RelationOp -> ShowS)
-> (RelationOp -> String)
-> ([RelationOp] -> ShowS)
-> Show RelationOp
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> RelationOp -> ShowS
showsPrec :: Int -> RelationOp -> ShowS
$cshow :: RelationOp -> String
show :: RelationOp -> String
$cshowList :: [RelationOp] -> ShowS
showList :: [RelationOp] -> ShowS
Show)

{- | Symbolic comparison of two expressions; use simplifyRelation to decide or preserve
unresolved conditions.
-}
data Relation = Relation RelationOp Expr Expr deriving (Relation -> Relation -> Bool
(Relation -> Relation -> Bool)
-> (Relation -> Relation -> Bool) -> Eq Relation
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Relation -> Relation -> Bool
== :: Relation -> Relation -> Bool
$c/= :: Relation -> Relation -> Bool
/= :: Relation -> Relation -> Bool
Eq, Int -> Relation -> ShowS
[Relation] -> ShowS
Relation -> String
(Int -> Relation -> ShowS)
-> (Relation -> String) -> ([Relation] -> ShowS) -> Show Relation
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Relation -> ShowS
showsPrec :: Int -> Relation -> ShowS
$cshow :: Relation -> String
show :: Relation -> String
$cshowList :: [Relation] -> ShowS
showList :: [Relation] -> ShowS
Show)

-- | Decided truth value or unresolved symbolic relation under the provided assumptions.
data RelationResult = Truth Bool | Unresolved Relation deriving (RelationResult -> RelationResult -> Bool
(RelationResult -> RelationResult -> Bool)
-> (RelationResult -> RelationResult -> Bool) -> Eq RelationResult
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RelationResult -> RelationResult -> Bool
== :: RelationResult -> RelationResult -> Bool
$c/= :: RelationResult -> RelationResult -> Bool
/= :: RelationResult -> RelationResult -> Bool
Eq, Int -> RelationResult -> ShowS
[RelationResult] -> ShowS
RelationResult -> String
(Int -> RelationResult -> ShowS)
-> (RelationResult -> String)
-> ([RelationResult] -> ShowS)
-> Show RelationResult
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> RelationResult -> ShowS
showsPrec :: Int -> RelationResult -> ShowS
$cshow :: RelationResult -> String
show :: RelationResult -> String
$cshowList :: [RelationResult] -> ShowS
showList :: [RelationResult] -> ShowS
Show)

{- | Given options then Relation, simplify both sides and return a Boolean Truth when
decidable, otherwise a symbolic unresolved relation/conditions. Invalid expressions or
resource failures return ExprError.
-}
simplifyRelation :: SymbolicOptions -> Relation -> Either ExprError RelationResult
simplifyRelation :: SymbolicOptions -> Relation -> Either ExprError RelationResult
simplifyRelation SymbolicOptions
o (Relation RelationOp
op Expr
a Expr
b) = do
    x <- SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o Expr
a
    y <- simplify o b
    case (qvalue x, qvalue y) of
        (Just Rational
u, Just Rational
v) -> RelationResult -> Either ExprError RelationResult
forall a b. b -> Either a b
Right (Bool -> RelationResult
Truth (RelationOp -> Rational -> Rational -> Bool
forall a. Ord a => RelationOp -> a -> a -> Bool
compareRelation RelationOp
op Rational
u Rational
v))
        (Maybe Rational, Maybe Rational)
_
            | Expr
x Expr -> Expr -> Bool
forall a. Eq a => a -> a -> Bool
== Expr
y Bool -> Bool -> Bool
&& Assumptions -> Expr -> Bool
defined (SymbolicOptions -> Assumptions
symbolicAssumptions SymbolicOptions
o) Expr
x ->
                RelationResult -> Either ExprError RelationResult
forall a b. b -> Either a b
Right (Bool -> RelationResult
Truth (RelationOp
op RelationOp -> [RelationOp] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [RelationOp
Equal, RelationOp
LessEqual, RelationOp
GreaterEqual]))
            | Bool
otherwise -> RelationResult -> Either ExprError RelationResult
forall a b. b -> Either a b
Right (Relation -> RelationResult
Unresolved (RelationOp -> Expr -> Expr -> Relation
Relation RelationOp
op Expr
x Expr
y))

compareRelation :: (Ord a) => RelationOp -> a -> a -> Bool
compareRelation :: forall a. Ord a => RelationOp -> a -> a -> Bool
compareRelation RelationOp
Equal = a -> a -> Bool
forall a. Eq a => a -> a -> Bool
(==)
compareRelation RelationOp
NotEqual = a -> a -> Bool
forall a. Eq a => a -> a -> Bool
(/=)
compareRelation RelationOp
Less = a -> a -> Bool
forall a. Ord a => a -> a -> Bool
(<)
compareRelation RelationOp
LessEqual = a -> a -> Bool
forall a. Ord a => a -> a -> Bool
(<=)
compareRelation RelationOp
Greater = a -> a -> Bool
forall a. Ord a => a -> a -> Bool
(>)
compareRelation RelationOp
GreaterEqual = a -> a -> Bool
forall a. Ord a => a -> a -> Bool
(>=)

{- | Isolate a single affine variable; multiplying by a negative coefficient
flips the inequality. Equality and zero slopes are handled explicitly.
-}
solveLinearInequality :: SymbolicOptions -> String -> Relation -> Either ExprError RelationResult
solveLinearInequality :: SymbolicOptions
-> String -> Relation -> Either ExprError RelationResult
solveLinearInequality SymbolicOptions
o String
name (Relation RelationOp
op Expr
a Expr
b) = do
    p <- SymbolicOptions
-> String -> Expr -> Either ExprError (Polynomial Rational)
toPolynomial SymbolicOptions
o String
name (Expr -> Expr -> Expr
minus Expr
a Expr
b)
    case P.coefficients p of
        [] -> RelationResult -> Either ExprError RelationResult
forall a b. b -> Either a b
Right (Bool -> RelationResult
Truth (RelationOp -> Rational -> Rational -> Bool
forall a. Ord a => RelationOp -> a -> a -> Bool
compareRelation RelationOp
op (Rational
0 :: Rational) Rational
0))
        [Rational
c] -> RelationResult -> Either ExprError RelationResult
forall a b. b -> Either a b
Right (Bool -> RelationResult
Truth (RelationOp -> Rational -> Rational -> Bool
forall a. Ord a => RelationOp -> a -> a -> Bool
compareRelation RelationOp
op Rational
c Rational
0))
        [Rational
c, Rational
m] -> do
            bound <- SymbolicOptions -> Expr -> Either ExprError Expr
checked SymbolicOptions
o (Rational -> Expr
literal (-Rational
c Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
m))
            Right (Unresolved (Relation (if m < 0 then flipOp op else op) (variable name) bound))
        [Rational]
_ -> Expr -> String -> Either ExprError RelationResult
forall a. Expr -> String -> Either ExprError a
unsupported Expr
a String
"inequality must be affine in one variable"
  where
    flipOp :: RelationOp -> RelationOp
flipOp RelationOp
Less = RelationOp
Greater
    flipOp RelationOp
LessEqual = RelationOp
GreaterEqual
    flipOp RelationOp
Greater = RelationOp
Less
    flipOp RelationOp
GreaterEqual = RelationOp
LessEqual
    flipOp RelationOp
x = RelationOp
x

{- | Symbolic linear-system outcome distinguishing a determined/parametric solution from
inconsistency or unresolved computation.
-}
data LinearSolution = Inconsistent | ParametricSolution [(String, Expr)] [String]
    deriving (LinearSolution -> LinearSolution -> Bool
(LinearSolution -> LinearSolution -> Bool)
-> (LinearSolution -> LinearSolution -> Bool) -> Eq LinearSolution
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: LinearSolution -> LinearSolution -> Bool
== :: LinearSolution -> LinearSolution -> Bool
$c/= :: LinearSolution -> LinearSolution -> Bool
/= :: LinearSolution -> LinearSolution -> Bool
Eq, Int -> LinearSolution -> ShowS
[LinearSolution] -> ShowS
LinearSolution -> String
(Int -> LinearSolution -> ShowS)
-> (LinearSolution -> String)
-> ([LinearSolution] -> ShowS)
-> Show LinearSolution
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> LinearSolution -> ShowS
showsPrec :: Int -> LinearSolution -> ShowS
$cshow :: LinearSolution -> String
show :: LinearSolution -> String
$cshowList :: [LinearSolution] -> ShowS
showList :: [LinearSolution] -> ShowS
Show)

{- | Exact rational RREF for at most 16 variables/equations. Free variables keep
their original names; a zero row with nonzero RHS certifies inconsistency.
-}
solveLinearEquations ::
    SymbolicOptions -> [String] -> [(Expr, Expr)] -> Either ExprError LinearSolution
solveLinearEquations :: SymbolicOptions
-> [String] -> [(Expr, Expr)] -> Either ExprError LinearSolution
solveLinearEquations SymbolicOptions
o [String]
names [(Expr, Expr)]
equations = do
    SymbolicOptions -> [String] -> Either ExprError ()
boundedVariables SymbolicOptions
o [String]
names
    if [String] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [String]
names Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
16 Bool -> Bool -> Bool
|| [(Expr, Expr)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [(Expr, Expr)] -> [(Expr, Expr)]
forall a. Int -> [a] -> [a]
take Int
17 [(Expr, Expr)]
equations) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
16
        then ErrorKind -> Expr -> String -> Either ExprError ()
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
MemoryLimit (Rational -> Expr
literal Rational
0) String
"linear system limit 16"
        else () -> Either ExprError ()
forall a b. b -> Either a b
Right ()
    rows <- ((Expr, Expr) -> Either ExprError [Rational])
-> [(Expr, Expr)] -> Either ExprError [[Rational]]
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 (Expr, Expr) -> Either ExprError [Rational]
row [(Expr, Expr)]
equations
    (reduced, pivots) <- rref 0 0 rows []
    if any (\[Rational]
r -> (Rational -> Bool) -> [Rational] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0) (Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
take ([String] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [String]
names) [Rational]
r) Bool -> Bool -> Bool
&& [Rational] -> Rational
lastValue [Rational]
r Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0) reduced
        then Right Inconsistent
        else do
            let free = [Int
i | Int
i <- [Int
0 .. [String] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [String]
names Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
i Int -> [Int] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` ((Int, Int) -> Int) -> [(Int, Int)] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Int) -> Int
forall a b. (a, b) -> b
snd [(Int, Int)]
pivots]
            solved <-
                mapM
                    ( \(Int
r, Int
c) -> do
                        let values :: [Rational]
values = [[Rational]]
reduced [[Rational]] -> Int -> [Rational]
forall a. HasCallStack => [a] -> Int -> a
!! Int
r
                            expr :: Expr
expr =
                                (Expr -> Expr -> Expr) -> Expr -> [Expr] -> Expr
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl'
                                    Expr -> Expr -> Expr
minus
                                    (Rational -> Expr
literal ([Rational] -> Rational
lastValue [Rational]
values))
                                    [Expr -> Expr -> Expr
times (Rational -> Expr
literal ([Rational]
values [Rational] -> Int -> Rational
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)) (String -> Expr
variable ([String]
names [String] -> Int -> String
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)) | Int
j <- [Int]
free, [Rational]
values [Rational] -> Int -> Rational
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0]
                        value <- SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o Expr
expr
                        Right (names !! c, value)
                    )
                    pivots
            Right (ParametricSolution solved [names !! i | i <- free])
  where
    lastValue :: [Rational] -> Rational
lastValue = (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
Foldable.foldl' (\Rational
_ Rational
x -> Rational
x) Rational
0
    row :: (Expr, Expr) -> Either ExprError [Rational]
row (Expr
a, Expr
b) = do
        expanded <- SymbolicOptions -> Expr -> Either ExprError Expr
expand SymbolicOptions
o (Expr -> Expr -> Expr
minus Expr
a Expr
b)
        (c, coeff) <- affine expanded
        let r = [Rational -> String -> Map String Rational -> Rational
forall k a. Ord k => a -> k -> Map k a -> a
Map.findWithDefault Rational
0 String
n Map String Rational
coeff | String
n <- [String]
names] [Rational] -> [Rational] -> [Rational]
forall a. [a] -> [a] -> [a]
++ [-Rational
c]
        inputs o (map literal r)
        Right r
    affine :: Expr -> Either ExprError (Rational, Map String Rational)
affine Expr
e = case Expr -> ExprNode
expressionNode Expr
e of
        Literal Rational
q -> (Rational, Map String Rational)
-> Either ExprError (Rational, Map String Rational)
forall a b. b -> Either a b
Right (Rational
q, Map String Rational
forall k a. Map k a
Map.empty)
        Variable String
n | String
n String -> [String] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [String]
names -> (Rational, Map String Rational)
-> Either ExprError (Rational, Map String Rational)
forall a b. b -> Either a b
Right (Rational
0, String -> Rational -> Map String Rational
forall k a. k -> a -> Map k a
Map.singleton String
n Rational
1)
        Unary UnaryOp
Negate Expr
x -> do (c, m) <- Expr -> Either ExprError (Rational, Map String Rational)
affine Expr
x; Right (-c, Map.map negate m)
        Binary BinaryOp
Add Expr
x Expr
y -> Rational
-> Expr -> Expr -> Either ExprError (Rational, Map String Rational)
combine Rational
1 Expr
x Expr
y
        Binary BinaryOp
Subtract Expr
x Expr
y -> Rational
-> Expr -> Expr -> Either ExprError (Rational, Map String Rational)
combine (-Rational
1) Expr
x Expr
y
        Binary BinaryOp
Multiply (Expr SourcePos
_ (Literal Rational
q)) Expr
x -> do (c, m) <- Expr -> Either ExprError (Rational, Map String Rational)
affine Expr
x; Right (q * c, Map.map (q *) m)
        Binary BinaryOp
Multiply Expr
x (Expr SourcePos
_ (Literal Rational
q)) -> Expr -> Either ExprError (Rational, Map String Rational)
affine (Expr -> Expr -> Expr
times (Rational -> Expr
literal Rational
q) Expr
x)
        Binary BinaryOp
Divide Expr
x (Expr SourcePos
_ (Literal Rational
q)) | Rational
q Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0 -> Expr -> Either ExprError (Rational, Map String Rational)
affine (Expr -> Expr -> Expr
times (Rational -> Expr
literal (Rational -> Rational
forall a. Fractional a => a -> a
recip Rational
q)) Expr
x)
        ExprNode
_ -> Expr -> String -> Either ExprError (Rational, Map String Rational)
forall a. Expr -> String -> Either ExprError a
unsupported Expr
e String
"linear equations require rational affine expressions"
    combine :: Rational
-> Expr -> Expr -> Either ExprError (Rational, Map String Rational)
combine Rational
sign Expr
x Expr
y = do
        (c, m) <- Expr -> Either ExprError (Rational, Map String Rational)
affine Expr
x
        (d, n) <- affine y
        Right (c + sign * d, Map.unionWith (+) m (Map.map (sign *) n))
    rref :: Int
-> Int
-> [[Rational]]
-> [(Int, Int)]
-> Either ExprError ([[Rational]], [(Int, Int)])
rref Int
r Int
c [[Rational]]
rows [(Int, Int)]
pivots
        | Int
r Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= [[Rational]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Rational]]
rows Bool -> Bool -> Bool
|| Int
c Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= [String] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [String]
names = ([[Rational]], [(Int, Int)])
-> Either ExprError ([[Rational]], [(Int, Int)])
forall a b. b -> Either a b
Right ([[Rational]]
rows, [(Int, Int)]
pivots)
        | Bool
otherwise = case [Int
i | Int
i <- [Int
r .. [[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], [[Rational]]
rows [[Rational]] -> Int -> [Rational]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Rational] -> Int -> Rational
forall a. HasCallStack => [a] -> Int -> a
!! Int
c Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0] of
            [] -> Int
-> Int
-> [[Rational]]
-> [(Int, Int)]
-> Either ExprError ([[Rational]], [(Int, Int)])
rref Int
r (Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [[Rational]]
rows [(Int, Int)]
pivots
            Int
i : [Int]
_ -> do
                let pivot :: [Rational]
pivot = (Rational -> Rational) -> [Rational] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ ([[Rational]]
rows [[Rational]] -> Int -> [Rational]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Rational] -> Int -> Rational
forall a. HasCallStack => [a] -> Int -> a
!! Int
c)) ([[Rational]]
rows [[Rational]] -> Int -> [Rational]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
                    swapped :: [[Rational]]
swapped =
                        [if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
r then [[Rational]]
rows [[Rational]] -> Int -> [Rational]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i else if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i then [[Rational]]
rows [[Rational]] -> Int -> [Rational]
forall a. HasCallStack => [a] -> Int -> a
!! Int
r else [Rational]
value | (Int
j, [Rational]
value) <- [Int] -> [[Rational]] -> [(Int, [Rational])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[Rational]]
rows]
                    next :: [[Rational]]
next =
                        [ if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
r then [Rational]
pivot else (Rational -> Rational -> Rational)
-> [Rational] -> [Rational] -> [Rational]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Rational]
value ((Rational -> Rational) -> [Rational] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map (([Rational]
value [Rational] -> Int -> Rational
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
*) [Rational]
pivot)
                        | (Int
j, [Rational]
value) <- [Int] -> [[Rational]] -> [(Int, [Rational])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[Rational]]
swapped
                        ]
                SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o ((Rational -> Expr) -> [Rational] -> [Expr]
forall a b. (a -> b) -> [a] -> [b]
map Rational -> Expr
literal ([[Rational]] -> [Rational]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Rational]]
next))
                Int
-> Int
-> [[Rational]]
-> [(Int, Int)]
-> Either ExprError ([[Rational]], [(Int, Int)])
rref (Int
r Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [[Rational]]
next ([(Int, Int)]
pivots [(Int, Int)] -> [(Int, Int)] -> [(Int, Int)]
forall a. [a] -> [a] -> [a]
++ [(Int
r, Int
c)])

-- | Checked bounded rectangular expression matrix, constructed with symbolicMatrix.
newtype SymbolicMatrix = SymbolicMatrix [[Expr]] deriving (SymbolicMatrix -> SymbolicMatrix -> Bool
(SymbolicMatrix -> SymbolicMatrix -> Bool)
-> (SymbolicMatrix -> SymbolicMatrix -> Bool) -> Eq SymbolicMatrix
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: SymbolicMatrix -> SymbolicMatrix -> Bool
== :: SymbolicMatrix -> SymbolicMatrix -> Bool
$c/= :: SymbolicMatrix -> SymbolicMatrix -> Bool
/= :: SymbolicMatrix -> SymbolicMatrix -> Bool
Eq, Int -> SymbolicMatrix -> ShowS
[SymbolicMatrix] -> ShowS
SymbolicMatrix -> String
(Int -> SymbolicMatrix -> ShowS)
-> (SymbolicMatrix -> String)
-> ([SymbolicMatrix] -> ShowS)
-> Show SymbolicMatrix
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> SymbolicMatrix -> ShowS
showsPrec :: Int -> SymbolicMatrix -> ShowS
$cshow :: SymbolicMatrix -> String
show :: SymbolicMatrix -> String
$cshowList :: [SymbolicMatrix] -> ShowS
showList :: [SymbolicMatrix] -> ShowS
Show)

{- | Given options then nonempty rectangular Expr rows, construct a checked symbolic matrix
bounded to 8-by-8. Empty/ragged/oversized or invalid expressions fail.
-}
symbolicMatrix :: SymbolicOptions -> [[Expr]] -> Either ExprError SymbolicMatrix
symbolicMatrix :: SymbolicOptions -> [[Expr]] -> Either ExprError SymbolicMatrix
symbolicMatrix SymbolicOptions
o [[Expr]]
rows
    | [[Expr]] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [[Expr]]
rows
        Bool -> Bool -> Bool
|| [[Expr]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [[Expr]] -> [[Expr]]
forall a. Int -> [a] -> [a]
take Int
9 [[Expr]]
rows) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
8
        Bool -> Bool -> Bool
|| ([Expr] -> Bool) -> [[Expr]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
width) (Int -> Bool) -> ([Expr] -> Int) -> [Expr] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Expr] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Expr] -> Int) -> ([Expr] -> [Expr]) -> [Expr] -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [Expr] -> [Expr]
forall a. Int -> [a] -> [a]
take Int
9) [[Expr]]
rows
        Bool -> Bool -> Bool
|| Int
width Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
        Bool -> Bool -> Bool
|| Int
width Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
8 =
        ErrorKind -> Expr -> String -> Either ExprError SymbolicMatrix
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
DomainFailure (Rational -> Expr
literal Rational
0) String
"symbolic matrix requires rectangular 1..8 dimensions"
    | Bool
otherwise = SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o ([[Expr]] -> [Expr]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Expr]]
rows) Either ExprError ()
-> Either ExprError SymbolicMatrix
-> Either ExprError SymbolicMatrix
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> SymbolicMatrix -> Either ExprError SymbolicMatrix
forall a b. b -> Either a b
Right ([[Expr]] -> SymbolicMatrix
SymbolicMatrix [[Expr]]
rows)
  where
    width :: Int
width = case [[Expr]]
rows of [Expr]
r : [[Expr]]
_ -> [Expr] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [Expr] -> [Expr]
forall a. Int -> [a] -> [a]
take Int
9 [Expr]
r); [] -> Int
0

-- | Extract row-major symbolic expressions without evaluation or simplification.
matrixExpressions :: SymbolicMatrix -> [[Expr]]
matrixExpressions :: SymbolicMatrix -> [[Expr]]
matrixExpressions (SymbolicMatrix [[Expr]]
rows) = [[Expr]]
rows

{- | Given options then equally shaped symbolic matrices, add and simplify corresponding
expressions. Shape mismatch or symbolic resource/domain failure returns ExprError.
-}
matrixAdd :: SymbolicOptions -> SymbolicMatrix -> SymbolicMatrix -> Either ExprError SymbolicMatrix
matrixAdd :: SymbolicOptions
-> SymbolicMatrix
-> SymbolicMatrix
-> Either ExprError SymbolicMatrix
matrixAdd SymbolicOptions
o (SymbolicMatrix [[Expr]]
a) (SymbolicMatrix [[Expr]]
b)
    | [[Expr]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Expr]]
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [[Expr]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Expr]]
b Bool -> Bool -> Bool
|| ([Expr] -> Int) -> [[Expr]] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map [Expr] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Expr]]
a [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
/= ([Expr] -> Int) -> [[Expr]] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map [Expr] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Expr]]
b =
        ErrorKind -> Expr -> String -> Either ExprError SymbolicMatrix
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
DomainFailure (Rational -> Expr
literal Rational
0) String
"matrix addition shape mismatch"
    | Bool
otherwise = (([Expr], [Expr]) -> Either ExprError [Expr])
-> [([Expr], [Expr])] -> Either ExprError [[Expr]]
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 ((Expr -> Either ExprError Expr)
-> [Expr] -> Either ExprError [Expr]
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 (SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o) ([Expr] -> Either ExprError [Expr])
-> (([Expr], [Expr]) -> [Expr])
-> ([Expr], [Expr])
-> Either ExprError [Expr]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Expr] -> [Expr] -> [Expr]) -> ([Expr], [Expr]) -> [Expr]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry ((Expr -> Expr -> Expr) -> [Expr] -> [Expr] -> [Expr]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Expr -> Expr -> Expr
plus)) ([[Expr]] -> [[Expr]] -> [([Expr], [Expr])]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Expr]]
a [[Expr]]
b) Either ExprError [[Expr]]
-> ([[Expr]] -> Either ExprError SymbolicMatrix)
-> Either ExprError SymbolicMatrix
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= SymbolicOptions -> [[Expr]] -> Either ExprError SymbolicMatrix
symbolicMatrix SymbolicOptions
o

{- | Given options then compatible symbolic matrices, form simplified dot-product entries.
Inner-dimension mismatch or expression/output limits return ExprError.
-}
matrixMultiply ::
    SymbolicOptions -> SymbolicMatrix -> SymbolicMatrix -> Either ExprError SymbolicMatrix
matrixMultiply :: SymbolicOptions
-> SymbolicMatrix
-> SymbolicMatrix
-> Either ExprError SymbolicMatrix
matrixMultiply SymbolicOptions
o (SymbolicMatrix [[Expr]]
a) (SymbolicMatrix [[Expr]]
b)
    | [[Expr]] -> Int
forall {t :: * -> *} {a}. Foldable t => [t a] -> Int
columns [[Expr]]
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [[Expr]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Expr]]
b = ErrorKind -> Expr -> String -> Either ExprError SymbolicMatrix
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
DomainFailure (Rational -> Expr
literal Rational
0) String
"matrix multiplication shape mismatch"
    | Bool
otherwise =
        ([Expr] -> Either ExprError [Expr])
-> [[Expr]] -> Either ExprError [[Expr]]
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
            ( \[Expr]
row ->
                (Int -> Either ExprError Expr) -> [Int] -> Either ExprError [Expr]
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
j -> SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o ((Expr -> Expr -> Expr) -> Expr -> [Expr] -> Expr
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl' Expr -> Expr -> Expr
plus (Rational -> Expr
literal Rational
0) ((Expr -> Expr -> Expr) -> [Expr] -> [Expr] -> [Expr]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Expr -> Expr -> Expr
times [Expr]
row [[Expr]
r [Expr] -> Int -> Expr
forall a. HasCallStack => [a] -> Int -> a
!! Int
j | [Expr]
r <- [[Expr]]
b])))
                    [Int
0 .. [[Expr]] -> Int
forall {t :: * -> *} {a}. Foldable t => [t a] -> Int
columns [[Expr]]
b Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            )
            [[Expr]]
a
            Either ExprError [[Expr]]
-> ([[Expr]] -> Either ExprError SymbolicMatrix)
-> Either ExprError SymbolicMatrix
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= SymbolicOptions -> [[Expr]] -> Either ExprError SymbolicMatrix
symbolicMatrix SymbolicOptions
o
  where
    columns :: [t a] -> Int
columns [t a]
rows = case [t a]
rows of t a
r : [t a]
_ -> t a -> Int
forall a. t a -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length t a
r; [] -> Int
0

{- | Cofactor determinant, at most 6x6. This bounded reference algorithm permits
arbitrary symbolic entries without assuming nonzero division pivots.
-}
symbolicDeterminant :: SymbolicOptions -> SymbolicMatrix -> Either ExprError Expr
symbolicDeterminant :: SymbolicOptions -> SymbolicMatrix -> Either ExprError Expr
symbolicDeterminant SymbolicOptions
o (SymbolicMatrix [[Expr]]
rows)
    | [[Expr]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Expr]]
rows Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
6 Bool -> Bool -> Bool
|| ([Expr] -> Bool) -> [[Expr]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [[Expr]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Expr]]
rows) (Int -> Bool) -> ([Expr] -> Int) -> [Expr] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Expr] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Expr]]
rows =
        ErrorKind -> Expr -> String -> Either ExprError Expr
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
DomainFailure (Rational -> Expr
literal Rational
0) String
"determinant requires square matrix of order <=6"
    | Bool
otherwise = [[Expr]] -> Either ExprError Expr
det [[Expr]]
rows Either ExprError Expr
-> (Expr -> Either ExprError Expr) -> Either ExprError Expr
forall a b.
Either ExprError a
-> (a -> Either ExprError b) -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o
  where
    det :: [[Expr]] -> Either ExprError Expr
det [] = Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right (Rational -> Expr
literal Rational
1)
    det [[Expr]
r] = case [Expr]
r of [Expr
x] -> Expr -> Either ExprError Expr
forall a b. b -> Either a b
Right Expr
x; [Expr]
_ -> ErrorKind -> Expr -> String -> Either ExprError Expr
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
DomainFailure (Rational -> Expr
literal Rational
0) String
"invalid determinant row"
    det ([Expr]
r : [[Expr]]
rs) = do
        terms <-
            ((Int, Expr) -> Either ExprError Expr)
-> [(Int, Expr)] -> Either ExprError [Expr]
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
j, Expr
x) -> do
                    sub <- [[Expr]] -> Either ExprError Expr
det [Int -> [Expr] -> [Expr]
forall a. Int -> [a] -> [a]
take Int
j [Expr]
row [Expr] -> [Expr] -> [Expr]
forall a. [a] -> [a] -> [a]
++ Int -> [Expr] -> [Expr]
forall a. Int -> [a] -> [a]
drop (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Expr]
row | [Expr]
row <- [[Expr]]
rs]
                    checked o (times (literal (if even j then 1 else -1)) (times x sub))
                )
                ([Int] -> [Expr] -> [(Int, Expr)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Expr]
r)
        checked o (Foldable.foldl' plus (literal 0) terms)

-- | Value paired with mathematical property conditions required for its unconditional use.
data Guarded a = Guarded
    { forall a. Guarded a -> a
guardedValue :: a
    -- ^ Result whose validity depends on requiredConditions.
    , forall a. Guarded a -> [(Expr, Property)]
requiredConditions :: [(Expr, Property)]
    -- ^ Expression/property facts that must hold before using guardedValue unconditionally.
    }
    deriving (Guarded a -> Guarded a -> Bool
(Guarded a -> Guarded a -> Bool)
-> (Guarded a -> Guarded a -> Bool) -> Eq (Guarded a)
forall a. Eq a => Guarded a -> Guarded a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => Guarded a -> Guarded a -> Bool
== :: Guarded a -> Guarded a -> Bool
$c/= :: forall a. Eq a => Guarded a -> Guarded a -> Bool
/= :: Guarded a -> Guarded a -> Bool
Eq, Int -> Guarded a -> ShowS
[Guarded a] -> ShowS
Guarded a -> String
(Int -> Guarded a -> ShowS)
-> (Guarded a -> String)
-> ([Guarded a] -> ShowS)
-> Show (Guarded a)
forall a. Show a => Int -> Guarded a -> ShowS
forall a. Show a => [Guarded a] -> ShowS
forall a. Show a => Guarded a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Guarded a -> ShowS
showsPrec :: Int -> Guarded a -> ShowS
$cshow :: forall a. Show a => Guarded a -> String
show :: Guarded a -> String
$cshowList :: forall a. Show a => [Guarded a] -> ShowS
showList :: [Guarded a] -> ShowS
Show)

{- | Exact 1x1/2x2 real eigenvalue formulas; an unproved discriminant is returned
as an explicit NonNegative condition. Negative numeric discriminants fail.
-}
symbolicEigenvalues :: SymbolicOptions -> SymbolicMatrix -> Either ExprError (Guarded [Expr])
symbolicEigenvalues :: SymbolicOptions
-> SymbolicMatrix -> Either ExprError (Guarded [Expr])
symbolicEigenvalues SymbolicOptions
o (SymbolicMatrix [[Expr]]
rows) = case [[Expr]]
rows of
    [[Expr
x]] -> SymbolicOptions -> [Expr] -> Either ExprError ()
inputs SymbolicOptions
o [Expr
x] Either ExprError ()
-> Either ExprError (Guarded [Expr])
-> Either ExprError (Guarded [Expr])
forall a b.
Either ExprError a -> Either ExprError b -> Either ExprError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Guarded [Expr] -> Either ExprError (Guarded [Expr])
forall a b. b -> Either a b
Right ([Expr] -> [(Expr, Property)] -> Guarded [Expr]
forall a. a -> [(Expr, Property)] -> Guarded a
Guarded [Expr
x] [])
    [[Expr
a, Expr
b], [Expr
c, Expr
d]] -> do
        trace <- SymbolicOptions -> Expr -> Either ExprError Expr
simplify SymbolicOptions
o (Expr -> Expr -> Expr
plus Expr
a Expr
d)
        disc <- expand o (plus (power (minus a d) (literal 2)) (times (literal 4) (times b c)))
        if maybe False (< 0) (qvalue disc)
            then unsupported disc "complex eigenvalues require a complex symbolic domain"
            else do
                let root = String -> [Expr] -> Expr
call String
"sqrt" [Expr
disc]
                lo <- simplify o (divide (minus trace root) (literal 2))
                hi <- simplify o (divide (plus trace root) (literal 2))
                Right
                    ( Guarded
                        [lo, hi]
                        (if proves (symbolicAssumptions o) NonNegative disc then [] else [(disc, NonNegative)])
                    )
    [[Expr]]
_ -> Expr -> String -> Either ExprError (Guarded [Expr])
forall a. Expr -> String -> Either ExprError a
unsupported (Rational -> Expr
literal Rational
0) String
"symbolic eigenvalues support order 1 or 2"

{- | Given options,precision,bindings list,then Expr, evaluate using the default builtin
registry and symbolic limits. Later duplicate binding names replace earlier ones.
Evaluation/name/domain/resource errors remain ExprError.
-}
numericalSubstitution ::
    SymbolicOptions -> Precision -> [(String, Value)] -> Expr -> Either ExprError Value
numericalSubstitution :: SymbolicOptions
-> Precision -> [(String, Value)] -> Expr -> Either ExprError Value
numericalSubstitution SymbolicOptions
o Precision
mode [(String, Value)]
bindings = Limits
-> Precision
-> Registry
-> Bindings
-> Expr
-> Either ExprError Value
evaluateExpression (SymbolicOptions -> Limits
symbolicLimits SymbolicOptions
o) Precision
mode Registry
defaultRegistry ([(String, Value)] -> Bindings
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList [(String, Value)]
bindings)

-- | TeX with escaped identifiers; expressions are never interpreted as raw TeX.
toLaTeX :: Limits -> Expr -> Either ExprError String
toLaTeX :: Limits -> Expr -> Either ExprError String
toLaTeX Limits
limits Expr
e = Limits -> (Expr -> String) -> Expr -> Either ExprError String
export Limits
limits Expr -> String
tex Expr
e
  where
    tex :: Expr -> String
tex (Expr SourcePos
_ ExprNode
node) = case ExprNode
node of
        Literal Rational
q
            | Rational -> Integer
forall a. Ratio a -> a
denominator Rational
q Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1 -> Integer -> String
forall a. Show a => a -> String
show (Rational -> Integer
forall a. Ratio a -> a
numerator Rational
q)
            | Bool
otherwise -> String
"\\frac{" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show (Rational -> Integer
forall a. Ratio a -> a
numerator Rational
q) String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"}{" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show (Rational -> Integer
forall a. Ratio a -> a
denominator Rational
q) String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"}"
        Variable String
n -> String
"\\mathrm{" String -> ShowS
forall a. [a] -> [a] -> [a]
++ (Char -> String) -> ShowS
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap Char -> String
escape String
n String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"}"
        Constant NamedConstant
Pi -> String
"\\pi"
        Constant NamedConstant
Euler -> String
"e"
        Constant NamedConstant
Tau -> String
"2\\pi"
        Unary UnaryOp
Negate Expr
a -> String
"-\\left(" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
tex Expr
a String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"\\right)"
        Binary BinaryOp
Divide Expr
a Expr
b -> String
"\\frac{" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
tex Expr
a String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"}{" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
tex Expr
b String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"}"
        Binary BinaryOp
Power Expr
a Expr
b -> String
"{\\left(" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
tex Expr
a String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"\\right)}^{" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
tex Expr
b String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"}"
        Binary BinaryOp
op Expr
a Expr
b ->
            String
"\\left(" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
tex Expr
a String -> ShowS
forall a. [a] -> [a] -> [a]
++ (if BinaryOp
op BinaryOp -> BinaryOp -> Bool
forall a. Eq a => a -> a -> Bool
== BinaryOp
Multiply then String
"\\cdot " else BinaryOp -> String
operatorText BinaryOp
op) String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
tex Expr
b String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"\\right)"
        Call String
n [Expr]
xs ->
            String
"\\operatorname{" String -> ShowS
forall a. [a] -> [a] -> [a]
++ (Char -> String) -> ShowS
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap Char -> String
escape String
n String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"}\\left(" String -> ShowS
forall a. [a] -> [a] -> [a]
++ String -> [String] -> String
forall a. [a] -> [[a]] -> [a]
intercalate String
"," ((Expr -> String) -> [Expr] -> [String]
forall a b. (a -> b) -> [a] -> [b]
map Expr -> String
tex [Expr]
xs) String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"\\right)"
        Let String
n Expr
a Expr
b ->
            String
"\\operatorname{let}\\;\\mathrm{"
                String -> ShowS
forall a. [a] -> [a] -> [a]
++ (Char -> String) -> ShowS
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap Char -> String
escape String
n
                String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"}="
                String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
tex Expr
a
                String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"\\;\\operatorname{in}\\;"
                String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
tex Expr
b
    escape :: Char -> String
escape Char
'_' = String
"\\_"
    escape Char
c = [Char
c]

operatorText :: BinaryOp -> String
operatorText :: BinaryOp -> String
operatorText BinaryOp
Add = String
"+"
operatorText BinaryOp
Subtract = String
"-"
operatorText BinaryOp
Multiply = String
"*"
operatorText BinaryOp
Divide = String
"/"
operatorText BinaryOp
Power = String
"^"

-- | Content MathML, including an explicit lambda application for lexical let.
toMathML :: Limits -> Expr -> Either ExprError String
toMathML :: Limits -> Expr -> Either ExprError String
toMathML Limits
limits Expr
e =
    Limits -> (Expr -> String) -> Expr -> Either ExprError String
export Limits
limits (\Expr
x -> String
"<math xmlns=\"http://www.w3.org/1998/Math/MathML\">" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
xml Expr
x String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"</math>") Expr
e
  where
    xml :: Expr -> String
xml (Expr SourcePos
_ ExprNode
node) = case ExprNode
node of
        Literal Rational
q
            | Rational -> Integer
forall a. Ratio a -> a
denominator Rational
q Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1 -> String
"<cn>" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show (Rational -> Integer
forall a. Ratio a -> a
numerator Rational
q) String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"</cn>"
            | Bool
otherwise ->
                String
"<cn type=\"rational\">" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show (Rational -> Integer
forall a. Ratio a -> a
numerator Rational
q) String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"<sep/>" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show (Rational -> Integer
forall a. Ratio a -> a
denominator Rational
q) String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"</cn>"
        Variable String
n -> String
"<ci>" String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
n String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"</ci>"
        Constant NamedConstant
Pi -> String
"<pi/>"
        Constant NamedConstant
Euler -> String
"<exponentiale/>"
        Constant NamedConstant
Tau -> String
"<apply><times/><cn>2</cn><pi/></apply>"
        Unary UnaryOp
Negate Expr
a -> String
"<apply><minus/>" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
xml Expr
a String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"</apply>"
        Binary BinaryOp
op Expr
a Expr
b -> String
"<apply><" String -> ShowS
forall a. [a] -> [a] -> [a]
++ BinaryOp -> String
tag BinaryOp
op String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"/>" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
xml Expr
a String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
xml Expr
b String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"</apply>"
        Call String
n [Expr]
xs -> String
"<apply><ci>" String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
n String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"</ci>" String -> ShowS
forall a. [a] -> [a] -> [a]
++ (Expr -> String) -> [Expr] -> String
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap Expr -> String
xml [Expr]
xs String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"</apply>"
        Let String
n Expr
a Expr
b ->
            String
"<apply><lambda><bvar><ci>" String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
n String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"</ci></bvar>" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
xml Expr
b String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"</lambda>" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Expr -> String
xml Expr
a String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"</apply>"
    tag :: BinaryOp -> String
tag BinaryOp
Add = String
"plus"
    tag BinaryOp
Subtract = String
"minus"
    tag BinaryOp
Multiply = String
"times"
    tag BinaryOp
Divide = String
"divide"
    tag BinaryOp
Power = String
"power"

export :: Limits -> (Expr -> String) -> Expr -> Either ExprError String
export :: Limits -> (Expr -> String) -> Expr -> Either ExprError String
export Limits
limits Expr -> String
render Expr
e = do
    Limits -> Expr -> Either ExprError ()
validateExpression Limits
limits Expr
e
    let text :: String
text = Expr -> String
render Expr
e
    if String -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> ShowS
forall a. Int -> [a] -> [a]
take (Limits -> Int
maximumInputChars Limits
limits Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) String
text) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Limits -> Int
maximumInputChars Limits
limits
        then ErrorKind -> Expr -> String -> Either ExprError String
forall a. ErrorKind -> Expr -> String -> Either ExprError a
failure ErrorKind
MemoryLimit Expr
e String
"export character budget"
        else String -> Either ExprError String
forall a b. b -> Either a b
Right String
text