tinker: Cylindrical Algebraic Decomposition

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

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

src/Theory/NonLinearRealArithmatic/Polynomial.hs

@@ -7,6 +7,7 @@ module Theory.NonLinearRealArithmatic.Polynomial
   , Monomial
   , mkPolynomial
   , exponentOf
+  , coefficientOf
   , Term(..)
   , modifyCoeff
   , modifyMonomial
@@ -16,26 +17,22 @@ 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
 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

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

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)