diff --git a/src/Theory/LinearArithmatic/Constraint.hs b/src/Theory/LinearArithmatic/Constraint.hs index 6198e9f..72b21b5 100644 --- a/src/Theory/LinearArithmatic/Constraint.hs +++ b/src/Theory/LinearArithmatic/Constraint.hs @@ -47,6 +47,16 @@ data AffineExpr = AffineExpr { getConstant :: Rational , getCoeffMap :: LinearExpr } +instance Num AffineExpr where + AffineExpr constA coeffsA + AffineExpr constB coeffsB = + AffineExpr (constA+constB) (M.unionWith (+) coeffsA coeffsB) + (*) = error "AffineExpr not closed under multiplication" + abs = error "TODO" + signum = error "TODO" + fromInteger n = AffineExpr (fromInteger n) M.empty + negate (AffineExpr constant coeffs) = + AffineExpr (negate constant) (M.map negate coeffs) + instance Show AffineExpr where show (AffineExpr constant coeff_map) = let diff --git a/src/Theory/NonLinearRealArithmatic/CAD.hs b/src/Theory/NonLinearRealArithmatic/CAD.hs index f46f169..1314a5b 100644 --- a/src/Theory/NonLinearRealArithmatic/CAD.hs +++ b/src/Theory/NonLinearRealArithmatic/CAD.hs @@ -116,8 +116,8 @@ split = \case , Closed end end ] -isolateRoots :: forall a. (Fractional a, Ord a, Show a) => UnivariatePolynomial a -> [Interval a] -isolateRoots polynomial = concatMap go $ maybeToList $ cauchyBound polynomial +isolateRootsIn :: forall a. (Fractional a, Ord a, Show a) => Interval a -> UnivariatePolynomial a -> [Interval a] +isolateRootsIn initial_interval polynomial = go initial_interval where sturm_seq = sturmSequence polynomial @@ -129,3 +129,9 @@ isolateRoots polynomial = concatMap go $ maybeToList $ cauchyBound polynomial | otherwise = undefined where root_count = countRootsIn sturm_seq interval + +isolateRoots :: forall a. (Fractional a, Ord a, Show a) => UnivariatePolynomial a -> [Interval a] +isolateRoots polynomial = + case cauchyBound polynomial of + Nothing -> [] + Just interval -> isolateRootsIn interval polynomial diff --git a/src/Theory/NonLinearRealArithmatic/Constraint.hs b/src/Theory/NonLinearRealArithmatic/Constraint.hs index 2ba8b9e..e35af97 100644 --- a/src/Theory/NonLinearRealArithmatic/Constraint.hs +++ b/src/Theory/NonLinearRealArithmatic/Constraint.hs @@ -1,4 +1,3 @@ - module Theory.NonLinearRealArithmatic.Constraint ( Constraint , ConstraintRelation(..) diff --git a/src/Theory/NonLinearRealArithmatic/Polynomial.hs b/src/Theory/NonLinearRealArithmatic/Polynomial.hs index c74cc86..cb849dc 100644 --- a/src/Theory/NonLinearRealArithmatic/Polynomial.hs +++ b/src/Theory/NonLinearRealArithmatic/Polynomial.hs @@ -33,6 +33,7 @@ import Data.Maybe ( maybeToList, fromMaybe ) import Prelude hiding (map) import Control.Arrow ((&&&)) import Control.Exception (assert) +import qualified Theory.LinearArithmatic.Constraint as LIA -- | Map of variables to integer exponents. type Monomial = IntMap Int @@ -256,3 +257,33 @@ class Num a => Assignable a where eval :: Assignment a -> Polynomial a -> a eval assignment polynomial = sum (evalTerm assignment <$> getTerms polynomial) +type UnivariatePolynomial a = [ (Int, a) ] + +{-| + Extracts list of exponent/coefficient pairs, sorted by exponent. + Returns Nothing, if the polynomial is multivariate. +-} +toUnivariate :: Polynomial a -> Maybe (UnivariatePolynomial a) +toUnivariate polynomial = + case varsIn polynomial of + [ var ] -> Just + $ List.sortOn fst + $ fmap (exponentOf var . getMonomial &&& getCoeff) + $ getTerms polynomial + zero_or_two_or_more_vars -> Nothing + +-- | If possible, converts polynomial to affine expression. Otherwise returns `Nothing`. +toAffineExpr :: Polynomial Rational -> Maybe LIA.AffineExpr +toAffineExpr (Polynomial terms) = + let + from_term :: Term Rational -> Maybe LIA.AffineExpr + from_term (Term coeff monomial) = + case M.toList monomial of + -- constant term: + [] -> Just $ LIA.AffineExpr coeff M.empty + -- single variable with exponent one: + [ (1, var) ] -> Just $ LIA.AffineExpr 0 (M.singleton var coeff) + -- more than one variable and/or exponents greater one: + _non_affine_term -> Nothing + in + sum <$> traverse from_term terms diff --git a/src/Theory/NonLinearRealArithmatic/Subtropical.hs b/src/Theory/NonLinearRealArithmatic/Subtropical.hs index bcb3f81..263ae74 100644 --- a/src/Theory/NonLinearRealArithmatic/Subtropical.hs +++ b/src/Theory/NonLinearRealArithmatic/Subtropical.hs @@ -50,7 +50,7 @@ frame polynomial = undefined -- TODO $ L.filter ((/= 0) . getCoeff) (getTerms polynomial) findDominatingDirection :: (Num a, Ord a) => Polynomial a -> Maybe (Assignment Int) -findDominatingDirection terms = undefined +findDominatingDirection terms = error "TODO: compute using Simplex" where pos_terms = filter ((> 0) . getCoeff) (getTerms terms) @@ -132,7 +132,7 @@ intermediateRoot polynomial neg_sol pos_sol = Just [ (0, c), (1, b), (2, a) ] -> [ (-b + sqrt (b^2 - 4*a*c)) / (2*a) , (-b - sqrt (b^2 - 4*a*c)) / (2*a) ] - -- TODO: higher degree polynomial ==> use bisection? + -- TODO: higher degree polynomial ==> use CAD Just higher_degree_polynomial -> error "TODO: subtropical does not support higher degree polynomials yet" t_solution :: Assignment a