From 865e82a0639ac285b4ceb02fd6fc9299d9c66fbe Mon Sep 17 00:00:00 2001 From: Niklas Gruhn Date: Sat, 30 Dec 2023 01:07:33 +0100 Subject: [PATCH] tinker: Cylindrical Algebraic Decomposition --- SMT.cabal | 1 + src/Theory/NonLinearRealArithmatic/CAD.hs | 126 +++++++++++++++++ .../NonLinearRealArithmatic/Polynomial.hs | 44 +----- .../NonLinearRealArithmatic/Subtropical.hs | 13 +- .../UnivariatePolynomial.hs | 127 ++++++++++++++++++ src/Utils.hs | 12 ++ 6 files changed, 280 insertions(+), 43 deletions(-) create mode 100644 src/Theory/NonLinearRealArithmatic/UnivariatePolynomial.hs diff --git a/SMT.cabal b/SMT.cabal index e39dcf9..be3827d 100644 --- a/SMT.cabal +++ b/SMT.cabal @@ -24,6 +24,7 @@ library , Theory.NonLinearRealArithmatic.Interval , Theory.NonLinearRealArithmatic.IntervalUnion , Theory.NonLinearRealArithmatic.Polynomial + , Theory.NonLinearRealArithmatic.UnivariatePolynomial , Theory.NonLinearRealArithmatic.Constraint , Theory.NonLinearRealArithmatic.BoundedFloating , Theory.NonLinearRealArithmatic.IntervalConstraintPropagation diff --git a/src/Theory/NonLinearRealArithmatic/CAD.hs b/src/Theory/NonLinearRealArithmatic/CAD.hs index 53bdc57..f46f169 100644 --- a/src/Theory/NonLinearRealArithmatic/CAD.hs +++ b/src/Theory/NonLinearRealArithmatic/CAD.hs @@ -3,3 +3,129 @@ -- module Theory.NonLinearRealArithmatic.CAD where +import Theory.NonLinearRealArithmatic.UnivariatePolynomial (UnivariatePolynomial, viewLeadTerm) +import qualified Theory.NonLinearRealArithmatic.UnivariatePolynomial as Uni +import Utils (adjacentPairs, count) +import Data.Maybe (maybeToList) +import Debug.Trace + +data Interval a = Open a a | Closed a a + +instance Show a => Show (Interval a) where + show (Open start end) = "(" ++ show start ++ "," ++ show end ++ ")" + show (Closed start end) = "[" ++ show start ++ "," ++ show end ++ "]" + +getLowerBound :: Interval a -> a +getLowerBound = \case + Open lb _ -> lb + Closed lb _ -> lb + +getUpperBound :: Interval a -> a +getUpperBound = \case + Open _ ub -> ub + Closed _ ub -> ub + +isOpen :: Interval a -> Bool +isOpen (Open _ _) = True +isOpen (Closed _ _) = False + +isClosed :: Interval a -> Bool +isClosed = not . isOpen + +{-| + The Cauchy bound isolates the roots of a univariate polynomial down + to a closed interval. For example, a Cauchy bound of 15, means that + all roots are somewhere in the interval: + + [-15, 15] + + Returns `Nothing` if the polynomial only has a constant term. + Notice, if the polynomial is just a non-zero constant, it has no roots. + If it is zero, it's nothing but roots. +-} +cauchyBound :: (Fractional a, Ord a) => UnivariatePolynomial a -> Maybe (Interval a) +cauchyBound polynomial = do + ((_, max_degree_coeff), rest_polynomial) <- Uni.viewLeadTerm polynomial + let other_coeffs = map snd $ Uni.toList rest_polynomial + let bound = 1 + maximum [ abs (coeff / max_degree_coeff) | coeff <- other_coeffs ] + return $ Closed (-bound) bound + +-- | With the Sturm sequence, we can count the number of roots in an interval. +sturmSequence :: forall a. (Fractional a, Ord a) => UnivariatePolynomial a -> [UnivariatePolynomial a] +sturmSequence polynomial = takeWhile (/= 0) $ go polynomial (Uni.derivative polynomial) + where + go :: UnivariatePolynomial a -> UnivariatePolynomial a -> [UnivariatePolynomial a] + go p1 p2 = p1 : go p2 (negate $ snd $ p1 `Uni.divide` p2) + +countSignChanges :: (Eq a, Num a) => [a] -> Int +countSignChanges as = count (uncurry (/=)) $ adjacentPairs signs + where signs = filter (/= 0) $ map signum as + +-- | Count the roots of a univariate polynomial using Sturm sequence in a given interval. +countRootsIn :: (Eq a, Num a) => [UnivariatePolynomial a] -> Interval a -> Int +countRootsIn sturm_seq interval = + let + polynomial = head sturm_seq + + start = getLowerBound interval + end = getUpperBound interval + + sign_changes_start = countSignChanges $ map (Uni.eval start) sturm_seq + sign_changes_end = countSignChanges $ map (Uni.eval end) sturm_seq + + root_count = sign_changes_start - sign_changes_end + in + -- With the Sturm sequence we technically count the roots in the interval + -- + -- (start, end] + -- + -- That is, the interval is *open* on the side of the lower bound and *closed* + -- on the side of the upper bound. + case interval of + -- So to get the root count for the closed interval [start,end] we have to add + -- one, if `start` itself is a root. + Closed _ _ -> + if Uni.eval start polynomial == 0 then + root_count + 1 + else + root_count + + -- For the root count of the open interval (start,end) we have to subtract one + -- if `end` is itself a root. + Open _ _ -> + if Uni.eval end polynomial == 0 then + root_count - 1 + else + root_count + +split :: Fractional a => Interval a -> [Interval a] +split = \case + Open start end -> + let mid = (start + end) / 2 in + [ Open start mid + , Closed mid mid + , Open mid end + ] + + Closed start end -> + let mid = (start + end) / 2 in + [ Closed start start + , Open start mid + , Closed mid mid + , Open mid end + , Closed end end + ] + +isolateRoots :: forall a. (Fractional a, Ord a, Show a) => UnivariatePolynomial a -> [Interval a] +isolateRoots polynomial = concatMap go $ maybeToList $ cauchyBound polynomial + where + sturm_seq = sturmSequence polynomial + + go :: Interval a -> [Interval a] + go interval + | root_count == 0 = [] + | root_count == 1 = [ interval ] + | root_count >= 2 = concatMap go $ split interval + | otherwise = undefined + where + root_count = countRootsIn sturm_seq interval diff --git a/src/Theory/NonLinearRealArithmatic/Polynomial.hs b/src/Theory/NonLinearRealArithmatic/Polynomial.hs index e70f508..c74cc86 100644 --- a/src/Theory/NonLinearRealArithmatic/Polynomial.hs +++ b/src/Theory/NonLinearRealArithmatic/Polynomial.hs @@ -7,6 +7,7 @@ module Theory.NonLinearRealArithmatic.Polynomial , Monomial , mkPolynomial , exponentOf + , coefficientOf , Term(..) , modifyCoeff , modifyMonomial @@ -16,16 +17,11 @@ module Theory.NonLinearRealArithmatic.Polynomial , fromExpr , toExpr , map - , degree - , termDegree , Assignable(..) , Assignment - , UnivariatePolynomial - , toUnivariate ) where import Theory.NonLinearRealArithmatic.Expr (Expr(..), Var, BinaryOp(..), UnaryOp(..)) -import qualified Data.List as L import qualified Data.IntMap as M import qualified Data.IntSet as S import qualified Data.List as List @@ -33,9 +29,10 @@ import Data.IntMap ( IntMap ) import Data.Function ( on ) import Data.Containers.ListUtils ( nubOrd ) import Control.Monad ( guard ) -import Data.Maybe ( maybeToList ) +import Data.Maybe ( maybeToList, fromMaybe ) import Prelude hiding (map) import Control.Arrow ((&&&)) +import Control.Exception (assert) -- | Map of variables to integer exponents. type Monomial = IntMap Int @@ -63,6 +60,10 @@ modifyMonomial f (Term coeff monomial) = Term coeff (f monomial) newtype Polynomial a = Polynomial { getTerms :: [Term a] } +coefficientOf :: Num a => Monomial -> Polynomial a -> a +coefficientOf monomial = + maybe 0 getCoeff . List.find ((== monomial) . getMonomial) . getTerms + {-| In a normalized polynomial: @@ -223,23 +224,6 @@ isLinear monomial = is_zero || is_unit is_zero = null non_zero_exponents is_unit = [1] == non_zero_exponents -termDegree :: Term a -> Int -termDegree = sum . getMonomial - -degree :: (Num a, Ord a) => Polynomial a -> Int -degree = maximum . fmap termDegree . getTerms - --- TODO: does that make sense for multivariate polynomials? -cauchyBound :: (Fractional a, Ord a) => Polynomial a -> a -cauchyBound polynomial = 1 + maximum [ abs (coeff / top_coeff) | coeff <- coeffs ] - where - -- Extract highest degree term with non-zero coefficient. - (top_coeff : coeffs) = - filter (/= 0) - $ fmap getCoeff - $ List.sortOn (negate . termDegree) - $ getTerms polynomial - varsIn :: Polynomial a -> [Var] varsIn (Polynomial terms) = nubOrd (terms >>= M.keys . getMonomial) @@ -272,17 +256,3 @@ 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 - $ L.sortOn fst - $ fmap (exponentOf var . getMonomial &&& getCoeff) - $ getTerms polynomial - zero_or_two_or_more_vars -> Nothing diff --git a/src/Theory/NonLinearRealArithmatic/Subtropical.hs b/src/Theory/NonLinearRealArithmatic/Subtropical.hs index 1dd85e0..51f4425 100644 --- a/src/Theory/NonLinearRealArithmatic/Subtropical.hs +++ b/src/Theory/NonLinearRealArithmatic/Subtropical.hs @@ -8,8 +8,9 @@ import qualified Data.IntMap as M import qualified Data.IntMap.Merge.Lazy as MM import Theory.NonLinearRealArithmatic.Constraint (Constraint, ConstraintRelation (..)) import Theory.NonLinearRealArithmatic.Expr (Expr (..), Var, BinaryOp (..), substitute) -import Theory.NonLinearRealArithmatic.Polynomial +import Theory.NonLinearRealArithmatic.Polynomial (Monomial, Polynomial (..), Assignment, Assignable (eval), Term (..), fromExpr, toExpr, varsIn) import Control.Arrow ((&&&)) +import qualified Theory.NonLinearRealArithmatic.UnivariatePolynomial as Univariate {-| The frame of a polynomial is a set of points, obtained from the @@ -114,19 +115,19 @@ intermediateRoot polynomial neg_sol pos_sol = t_roots :: [a] t_roots = - case toUnivariate t_polynomial of + case Univariate.toList <$> Univariate.fromPolynomial t_polynomial of Nothing -> error "Polynomial should be univariate by construction" -- linear polynomial ==> solve directly for t - Just [ (0, c), (1, b) ] -> [ - b/c ] + Just [ (0, c), (1, b) ] -> [- b/c ] -- quadratic polynomial ==> apply quadratic formula - Just [ (0, c), (1, b), (2, a) ] -> + 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? Just higher_degree_polynomial -> error "TODO: subtropical does not support higher degree polynomials yet" t_solution :: Assignment a - t_solution = + t_solution = case L.find (\t' -> 0 <= t' && t' <= 1) t_roots of Just value -> M.singleton t value Nothing -> error "No solution for t between 0 and 1" @@ -141,7 +142,7 @@ subtropical (Equals, polynomial) = let -- assignment that maps all variables to one one :: Assignment a - one = M.fromList $ zip (varsIn polynomial) (repeat 1) + one = M.fromList $ map (, 1) (varsIn polynomial) go :: Assignment a -> Polynomial a -> Maybe (Assignment a) go neg_sol polynomial = intermediateRoot polynomial neg_sol <$> positiveSolution polynomial diff --git a/src/Theory/NonLinearRealArithmatic/UnivariatePolynomial.hs b/src/Theory/NonLinearRealArithmatic/UnivariatePolynomial.hs new file mode 100644 index 0000000..5aea713 --- /dev/null +++ b/src/Theory/NonLinearRealArithmatic/UnivariatePolynomial.hs @@ -0,0 +1,127 @@ +module Theory.NonLinearRealArithmatic.UnivariatePolynomial + ( term + , constant + -- WARN: don't expose type constructor to make sure invariants are enforced + , UnivariatePolynomial + , viewLeadTerm + , lookupLeadTerm + , toList + , fromList + , fromPolynomial + , toPolynomial + , derivative + , divide + , eval + ) where + +import Theory.NonLinearRealArithmatic.Polynomial (Polynomial) +import qualified Theory.NonLinearRealArithmatic.Polynomial as P +import Data.IntMap (IntMap) +import qualified Data.IntMap as M +import Control.Exception (assert) +import qualified Data.List as List +import Utils (assertM, adjacentPairs, count) + +-- TODO: property test: never contains zero coefficients +newtype UnivariatePolynomial a = Univariate { getTerms :: IntMap a } + deriving Eq + +instance Show a => Show (UnivariatePolynomial a) where + show (Univariate terms) = + List.intercalate " + " [ show coeff ++ "x^" ++ show exp | (exp, coeff) <- M.toDescList terms ] + +term :: (Eq a, Num a) => Int -> a -> UnivariatePolynomial a +term exp coeff + | coeff == 0 = Univariate M.empty + | otherwise = Univariate $ M.singleton exp coeff + +constant :: (Eq a, Num a) => a -> UnivariatePolynomial a +constant = term 0 + +-- | Extract term with largest exponent. +viewLeadTerm :: (Eq a, Num a) => UnivariatePolynomial a -> Maybe ((Int, a), UnivariatePolynomial a) +viewLeadTerm polynomial = do + ((exp, coeff), rest_terms) <- M.maxViewWithKey $ getTerms polynomial + -- assume that polynomial is normalized and that all terms with zero coefficient are removed: + assertM (coeff /= 0) + return ((exp, coeff), Univariate rest_terms) + +lookupLeadTerm :: (Eq a, Num a) => UnivariatePolynomial a -> Maybe (Int, a) +lookupLeadTerm = fmap fst . viewLeadTerm + +-- | Largest exponent of term with non-zero coefficient in the polynomial. +degree :: (Eq a, Num a, Show a) => UnivariatePolynomial a -> Int +degree = maybe 0 fst . lookupLeadTerm + +instance (Num a, Eq a) => Num (UnivariatePolynomial a) where + Univariate p + Univariate q = Univariate $ M.filter (/= 0) $ M.unionWith (+) p q + + Univariate p * Univariate q = Univariate $ M.fromList $ do + (exp_p, coeff_p) <- M.toList p + (exp_q, coeff_q) <- M.toList q + assertM $ coeff_p /= 0 && coeff_q /= 0 + return (exp_p + exp_q, coeff_p * coeff_q) + + abs = error "TODO: not implemented" + + signum = error "TODO: not implemented" + + fromInteger = constant . fromInteger + + negate (Univariate p) = Univariate $ M.map negate p + +-- | From list of exponent/coefficient pairs. +fromList :: [(Int, a)] -> UnivariatePolynomial a +fromList = Univariate . M.fromList + +-- | So list of exponent/coefficient pairs, ascendingly ordered by exponents. +toList :: UnivariatePolynomial a -> [(Int, a)] +toList = M.toAscList . getTerms + +-- | Extracts univariate polynomial. Returns `Nothing` if the polynomial is multivariate. +fromPolynomial :: forall a. (Num a, Eq a) => Polynomial a -> Maybe (UnivariatePolynomial a) +fromPolynomial polynomial = + case P.varsIn polynomial of + -- polynomial is just a constant: + [] -> Just $ constant $ P.coefficientOf M.empty polynomial + + -- polynomial has exactly one variable: + [ var ] -> Just $ Univariate $ M.fromList $ do + P.Term coeff monomial <- P.getTerms polynomial + let exp = P.exponentOf var monomial + return (exp, coeff) + + -- polynomial has two or more variables: + _two_or_more_variables -> Nothing + +toPolynomial :: (Num a, Ord a) => UnivariatePolynomial a -> Polynomial a +toPolynomial polynomial = P.mkPolynomial $ do + (exp, coeff) <- M.toList $ getTerms polynomial + return $ P.Term coeff (M.singleton 0 exp) + +eval :: Num a => a -> UnivariatePolynomial a -> a +eval x (Univariate terms) = sum [ coeff * x^exp | (exp, coeff) <- M.toList terms ] + +derivative :: forall a. Num a => UnivariatePolynomial a -> UnivariatePolynomial a +derivative (Univariate polynomial) = Univariate $ M.fromList + [ (exp-1, coeff * fromIntegral exp) | (exp, coeff) <- M.toList polynomial, exp > 0 ] + +-- | Perform polynomial division and return both quotient and remainder. +-- +-- TODO: property test division/multiplication. +divide :: forall a. (Fractional a, Eq a) => UnivariatePolynomial a -> UnivariatePolynomial a -> (UnivariatePolynomial a, UnivariatePolynomial a) +divide dividend divisor = go 0 dividend + where + go :: UnivariatePolynomial a -> UnivariatePolynomial a -> (UnivariatePolynomial a, UnivariatePolynomial a) + go quotient remainder = + case (lookupLeadTerm remainder, lookupLeadTerm divisor) of + (_, Nothing) -> error "division by zero in polynomial division" + + (Nothing, _) -> (quotient, 0) + + (Just (rem_exp, rem_coeff), Just (div_exp, div_coeff)) -> + if rem_exp < div_exp then + (quotient, remainder) + else + let quot_term = term (rem_exp - div_exp) (rem_coeff / div_coeff) + in go (quotient + quot_term) (remainder - quot_term * divisor) diff --git a/src/Utils.hs b/src/Utils.hs index 43c1206..4c70ad9 100644 --- a/src/Utils.hs +++ b/src/Utils.hs @@ -20,3 +20,15 @@ takeWhileJust = catMaybes . takeWhile isJust combinations :: [a] -> [(a,a)] combinations [] = [] combinations (a:as) = map (a,) as ++ combinations as + +assertM :: Monad m => Bool -> m () +assertM condition + | condition = return () + | otherwise = error "assertion failure" + +count :: (a -> Bool) -> [a] -> Int +count p = length . filter p + +adjacentPairs :: [a] -> [(a, a)] +adjacentPairs as = zip as (tail as) +