presburger 1.1 → 1.2
raw patch · 3 files changed
+283/−86 lines, 3 filesdep +QuickCheckdep +presburgerdep ~base
Dependencies added: QuickCheck, presburger
Dependency ranges changed: base
Files
- presburger.cabal +8/−2
- src/Data/Integer/SAT.hs +239/−84
- tests/qc.hs +36/−0
presburger.cabal view
@@ -1,5 +1,5 @@ Name: presburger-Version: 1.1+Version: 1.2 License: BSD3 License-file: LICENSE Author: Iavor S. Diatchki@@ -10,7 +10,7 @@ Description: The decision procedure is based on the algorithm used in CVC4, which is itself based on the Omega test. Build-type: Simple-Cabal-version: >= 1.6+Cabal-version: >= 1.8 library Build-Depends: base < 10, containers, pretty@@ -23,4 +23,10 @@ source-repository head type: git location: git://github.com/yav/presburger.git++Test-Suite pressburger-qc-tests+ type: exitcode-stdio-1.0+ hs-source-dirs: tests+ main-is: qc.hs+ build-depends: base, presburger == 1.2, QuickCheck
src/Data/Integer/SAT.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE Safe, PatternGuards #-}+{-# LANGUAGE Trustworthy, PatternGuards, BangPatterns #-} {-| This module implements a decision procedure for quantifier-free linear arithmetic. The algorithm is based on the following paper:@@ -21,15 +21,35 @@ , Name , toName , fromName+ -- * Iterators+ , allSolutions+ , slnCurrent+ , slnNextVal+ , slnNextVar+ , slnEnumerate+++ -- * Debug+ , dotPropSet+ , sizePropSet+ , allInerts+ , ppInerts++ -- * For QuickCheck+ , iPickBounded+ , Bound(..)+ , tConst ) where +import Debug.Trace+ import Data.Map (Map) import qualified Data.Map as Map import Data.List(partition) import Data.Maybe(maybeToList,fromMaybe,mapMaybe) import Control.Applicative(Applicative(..), (<$>))-import Control.Monad(liftM,ap,MonadPlus(..),msum,guard)-import Text.PrettyPrint(Doc,(<+>), (<>), integer, int, hsep, text)+import Control.Monad(liftM,ap,MonadPlus(..),guard)+import Text.PrettyPrint infixr 2 :|| infixr 3 :&&@@ -44,6 +64,12 @@ newtype PropSet = State (Answer RW) deriving Show +dotPropSet :: PropSet -> Doc+dotPropSet (State a) = dotAnswer (ppInerts . inerts) a++sizePropSet :: PropSet -> (Integer,Integer,Integer)+sizePropSet (State a) = answerSize a+ -- | An empty collection of propositions. noProps :: PropSet noProps = State $ return initRW@@ -63,6 +89,13 @@ go (One rw) = return [ (x,v) | (UserName x, v) <- iModel (inerts rw) ] go (Choice m1 m2) = mplus (go m1) (go m2) +allInerts :: PropSet -> [Inerts]+allInerts (State m) = map inerts (toList m)++allSolutions :: PropSet -> [Solutions]+allSolutions = map startIter . allInerts++ -- | Computes bounds on the expression that are compatible with the model. -- Returns `Nothing` if the bound is not known. getExprBound :: BoundType -> Expr -> PropSet -> Maybe Integer@@ -138,21 +171,21 @@ prop (e1 :== e2) = do t1 <- expr e1 t2 <- expr e2- enqAndGo qZeroTerms (t1 |-| t2)+ solveIs0 (t1 |-| t2) prop (e1 :/= e2) = do t1 <- expr e1 t2 <- expr e2 let t = t1 |-| t2- enqAndGo qNegTerms t `mplus` enqAndGo qNegTerms (tNeg t)+ solveIsNeg t `orElse` solveIsNeg (tNeg t) prop (e1 :< e2) = do t1 <- expr e1 t2 <- expr e2- enqAndGo qNegTerms (t1 |-| t2)+ solveIsNeg (t1 |-| t2) prop (e1 :<= e2) = do t1 <- expr e1 t2 <- expr e2- let t = t1 |-| t2- enqAndGo qZeroTerms t `mplus` enqAndGo qNegTerms t+ let t = t1 |-| t2 |-| tConst 1+ solveIsNeg t prop (e1 :> e2) = prop (e2 :< e1) prop (e1 :>= e2) = prop (e2 :<= e1)@@ -187,47 +220,11 @@ -------------------------------------------------------------------------------- data RW = RW { nameSource :: !Int- , todo :: WorkQ , inerts :: Inerts } deriving Show initRW :: RW-initRW = RW { nameSource = 0, todo = qEmpty, inerts = iNone }--solveAll :: S ()-solveAll =- do mbEq <- getWork qZeroTerms- case mbEq of- Just p -> solveIs0 p >> solveAll- Nothing ->- do mbLt <- getWork qNegTerms- case mbLt of- Just p -> solveIsNeg p >> solveAll- Nothing -> return ()-------------------------------------------------------------------------------------- The work queue--data WorkQ = WorkQ { zeroTerms :: [Term] -- ^ t == 0- , negTerms :: [Term] -- ^ t < 0- } deriving Show--qEmpty :: WorkQ-qEmpty = WorkQ { zeroTerms = [], negTerms = [] }--qLet :: Name -> Term -> WorkQ -> WorkQ-qLet x t q = WorkQ { zeroTerms = map (tLet x t) (zeroTerms q)- , negTerms = map (tLet x t) (negTerms q)- }--type Field t = (WorkQ -> [t], [t] -> WorkQ -> WorkQ)--qZeroTerms :: Field Term-qZeroTerms = (zeroTerms, \a q -> q { zeroTerms = a })--qNegTerms :: Field Term-qNegTerms = (negTerms, \a q -> q { negTerms = a })+initRW = RW { nameSource = 0, inerts = iNone } -------------------------------------------------------------------------------- -- Constraints and Bound on Variables@@ -265,7 +262,20 @@ -- These form an idempotent substitution. } deriving Show +ppInerts :: Inerts -> Doc+ppInerts is = vcat $ [ ppLower x b | (x,(ls,_)) <- bnds, b <- ls ] +++ [ ppUpper x b | (x,(_,us)) <- bnds, b <- us ] +++ [ ppEq e | e <- Map.toList (solved is) ]+ where+ bnds = Map.toList (bounds is) + ppT c x = ppTerm (c |*| tVar x)+ ppLower x (Bound c t) = ppTerm t <+> text "<" <+> ppT c x+ ppUpper x (Bound c t) = ppT c x <+> text "<" <+> ppTerm t+ ppEq (x,t) = ppName x <+> text "=" <+> ppTerm t+++ -- | An empty inert set. iNone :: Inerts iNone = Inerts { bounds = Map.empty@@ -321,6 +331,114 @@ stay (Bound _ bnd) = not (tHasVar x bnd) +-- | Given some lower and upper bounds, find the interval the satisfies them.+-- Note the upper and lower bounds are strict (i.e., < and >)+boundInterval :: [Bound] -> [Bound] -> Maybe (Maybe Integer, Maybe Integer)+boundInterval lbs ubs =+ do ls <- mapM (normBound Lower) lbs+ us <- mapM (normBound Upper) ubs+ let lb = case ls of+ [] -> Nothing+ _ -> Just (maximum ls + 1)+ ub = case us of+ [] -> Nothing+ _ -> Just (minimum us - 1)+ case (lb,ub) of+ (Just l, Just u) -> guard (l <= u)+ _ -> return ()+ return (lb,ub)+ where+ normBound Lower (Bound c t) = do k <- isConst t+ return (div (k + c - 1) c)+ normBound Upper (Bound c t) = do k <- isConst t+ return (div k c)++data Solutions = Done+ | TopVar Name Integer (Maybe Integer) (Maybe Integer) Inerts+ | FixedVar Name Integer Solutions+ deriving Show++slnCurrent :: Solutions -> [(Int,Integer)]+slnCurrent s = [ (x,v) | (UserName x, v) <- go s ]+ where+ go Done = []+ go (TopVar x v _ _ is) = (x, v) : iModel (iLet x v is)+ go (FixedVar x v i) = (x, v) : go i++-- | Replace occurances of a variable with an integer.+-- WARNING: The integer should be a valid value for the variable.+iLet :: Name -> Integer -> Inerts -> Inerts+iLet x v is = Inerts { bounds = fmap updBs (bounds is)+ , solved = fmap (tLetNum x v) (solved is) }+ where+ updB (Bound c t) = Bound c (tLetNum x v t)+ updBs (ls,us) = (map updB ls, map updB us)+++startIter :: Inerts -> Solutions+startIter is =+ case Map.maxViewWithKey (bounds is) of+ Nothing ->+ case Map.maxViewWithKey (solved is) of+ Nothing -> Done+ Just ((x,t), mp1) ->+ case [ y | y <- tVarList t ] of+ y : _ -> TopVar y 0 Nothing Nothing is+ [] -> let v = tConstPart t+ in TopVar x v (Just v) (Just v) $ is { solved = mp1 }+ Just ((x,(lbs,ubs)), mp1) ->+ case [ y | Bound _ t <- lbs ++ ubs, y <- tVarList t ] of+ y : _ -> TopVar y 0 Nothing Nothing is+ [] -> case boundInterval lbs ubs of+ Nothing -> error "bug: cannot compute interval?"+ Just (lb,ub) ->+ let v = fromMaybe 0 (mplus lb ub)+ in TopVar x v lb ub $ is { bounds = mp1 }++slnEnumerate :: Solutions -> [ Solutions ]+slnEnumerate s0 = go s0 []+ where+ go s k = case slnNextVar s of+ Nothing -> hor s k+ Just s1 -> go s1 $ case slnNextVal s of+ Nothing -> k+ Just s2 -> go s2 k++ hor s k = s+ : case slnNextVal s of+ Nothing -> k+ Just s1 -> hor s1 k++slnNextVal :: Solutions -> Maybe Solutions+slnNextVal Done = Nothing+slnNextVal (FixedVar x v i) = FixedVar x v `fmap` slnNextVal i+slnNextVal it@(TopVar _ _ lb _ _) =+ case lb of+ Just _ -> slnNextValWith (+1) it+ Nothing -> slnNextValWith (subtract 1) it+++slnNextValWith :: (Integer -> Integer) -> Solutions -> Maybe Solutions+slnNextValWith _ Done = Nothing+slnNextValWith f (FixedVar x v i) = FixedVar x v `fmap` slnNextValWith f i+slnNextValWith f (TopVar x v lb ub is) =+ do let v1 = f v+ case lb of+ Just l -> guard (l <= v1)+ Nothing -> return ()+ case ub of+ Just u -> guard (v1 <= u)+ Nothing -> return ()+ return $ TopVar x v1 lb ub is++slnNextVar :: Solutions -> Maybe Solutions+slnNextVar Done = Nothing+slnNextVar (TopVar x v _ _ is) = Just $ FixedVar x v $ startIter $ iLet x v is+slnNextVar (FixedVar x v i) = FixedVar x v `fmap` slnNextVar i++++ -- Given a list of lower (resp. upper) bounds, compute the least (resp. largest) -- value that satisfies them all. iPickBounded :: BoundType -> [Bound] -> Maybe Integer@@ -328,13 +446,23 @@ iPickBounded bt bs = do xs <- mapM (normBound bt) bs return $ case bt of- Lower -> maximum xs + 1- Upper -> minimum xs - 1+ Lower -> maximum xs+ Upper -> minimum xs where+ -- t < c*x+ -- <=> t+1 <= c*x+ -- <=> (t+1)/c <= x+ -- <=> ceil((t+1)/c) <= x+ -- <=> t `div` c + 1 <= x normBound Lower (Bound c t) = do k <- isConst t- return (div (k + c - 1) c)+ return (k `div` c + 1)+ -- c*x < t+ -- <=> c*x <= t-1+ -- <=> x <= (t-1)/c+ -- <=> x <= floor((t-1)/c)+ -- <=> x <= (t-1) `div` c normBound Upper (Bound c t) = do k <- isConst t- return (div k c)+ return (div (k-1) c) -- | The largest (resp. least) upper (resp. lower) bound on a term@@ -377,6 +505,7 @@ + iModel :: Inerts -> [(Name,Integer)] iModel i = goBounds [] (bounds i) where@@ -401,10 +530,13 @@ -------------------------------------------------------------------------------- -- Solving constraints +solveIs0 :: Term -> S ()+solveIs0 t = solveIs0' =<< apSubst t+ -- | Solve a constraint if the form @t = 0@. -- Assumes substitution has already been applied.-solveIs0 :: Term -> S ()-solveIs0 t+solveIs0' :: Term -> S ()+solveIs0' t -- A == 0 | Just a <- isConst t = guard (a == 0)@@ -415,12 +547,13 @@ (q,0) -> addDef x (tConst q) _ -> mzero - -- x + S = 0+ -- x + S = 0+ -- -x + S = 0 | Just (xc,x,s) <- tGetSimpleCoeff t = addDef x (if xc > 0 then tNeg s else s) -- A * S = 0- | Just (_, s) <- tFactor t = addWork qZeroTerms s+ | Just (_, s) <- tFactor t = solveIs0 s -- See Section 3.1 of paper for details. -- We obtain an equivalent formulation but with smaller coefficients.@@ -433,7 +566,7 @@ addDef xk soln let upd i = div (2*i + m) (2*m) + modulus i m- addWork qZeroTerms (negate (abs ak) |*| tVar v |+| tMapCoeff upd s)+ solveIs0 (negate (abs ak) |*| tVar v |+| tMapCoeff upd s) | otherwise = error "solveIs0: unreachable" @@ -441,16 +574,20 @@ modulus a m = a - m * div (2 * a + m) (2 * m) +solveIsNeg :: Term -> S ()+solveIsNeg t = solveIsNeg' =<< apSubst t++ -- | Solve a constraint of the form @t < 0@. -- Assumes that substitution has been applied-solveIsNeg :: Term -> S ()-solveIsNeg t+solveIsNeg' :: Term -> S ()+solveIsNeg' t -- A < 0 | Just a <- isConst t = guard (a < 0) -- A * S < 0- |Just (_,s) <- tFactor t = addWork qNegTerms s+ | Just (_,s) <- tFactor t = solveIsNeg s -- See Section 5.1 of the paper | Just (xc,x,s) <- tLeastVar t =@@ -477,12 +614,14 @@ dark = ctLt (tConst (a * b)) (b |*| alpha |-| a |*| beta) gray = [ ctEq (b |*| tVar x) (tConst i |+| beta) | i <- [ 1 .. b - 1 ] ]- addWork qNegTerms real- msum (addWork qNegTerms dark : map (addWork qZeroTerms) gray)+ solveIsNeg real+ foldl orElse (solveIsNeg dark) (map solveIs0 gray) ) ctrs | otherwise = error "solveIsNeg: unreachable" +orElse :: S () -> S () -> S ()+orElse x y = mplus x y {- Note [Shadows] @@ -522,6 +661,40 @@ data Answer a = None | One a | Choice (Answer a) (Answer a) deriving Show ++answerSize :: Answer a -> (Integer,Integer,Integer)+answerSize = go 0 0 0+ where+ go !n !o !c ans =+ case ans of+ None -> (n+1, o, c)+ One _ -> (n, o + 1, c)+ Choice x y ->+ case go n o (c+1) x of+ (n',o',c') -> go n' o' c' y+++dotAnswer :: (a -> Doc) -> Answer a -> Doc+dotAnswer pp g0 = vcat [text "digraph {", nest 2 (fst $ go 0 g0), text "}"]+ where+ node x d = integer x <+> brackets (text "label=" <> text (show d))+ <> semi+ edge x y = integer x <+> text "->" <+> integer y++ go x None = let x' = x + 1+ in seq x' ( node x "", x' )+ go x (One a) = let x' = x + 1+ in seq x' ( node x (show (pp a)), x' )+ go x (Choice c1 c2) = let x' = x + 1+ (ls1,x1) = go x' c1+ (ls2,x2) = go x1 c2+ in seq x'+ ( vcat [ node x "|"+ , edge x x'+ , edge x x1+ , ls1+ , ls2+ ], x2 ) toList :: Answer a -> [a] toList a = go a [] where@@ -586,20 +759,6 @@ , rw { nameSource = nameSource rw + 1 } ) --- | Try to get a new item from the work queue.-getWork :: Field t -> S (Maybe t)-getWork (getF,setF) = updS $ \rw ->- let work = todo rw- in case getF work of- [] -> (Nothing, rw)- t : ts -> (Just t, rw { todo = setF ts work })---- | Add a new item to the work queue.-addWork :: Field t -> t -> S ()-addWork (getF,setF) t = updS_ $ \rw ->- let work = todo rw- in rw { todo = setF (t : getF work) work }- -- | Get lower ('fst'), or upper ('snd') bounds for a variable. getBounds :: BoundType -> Name -> S [Bound] getBounds f x = get $ \rw -> case Map.lookup x $ bounds $ inerts rw of@@ -620,19 +779,15 @@ -- | Add a new definition. -- Assumes substitution has already been applied addDef :: Name -> Term -> S ()-addDef x t = updS_ $ \rw ->- let (newWork,newInerts) = iSolved x t (inerts rw)- in rw { inerts = newInerts- , todo = qLet x t $- let work = todo rw- in work { negTerms = newWork ++ negTerms work }- }+addDef x t =+ do newWork <- updS $ \rw -> let (newWork,newInerts) = iSolved x t (inerts rw)+ in (newWork, rw { inerts = newInerts })+ mapM_ solveIsNeg newWork -enqAndGo :: Field Term -> Term -> S ()-enqAndGo q t =+apSubst :: Term -> S Term+apSubst t = do i <- get inerts- addWork q $ iApSubst i t- solveAll+ return (iApSubst i t)
+ tests/qc.hs view
@@ -0,0 +1,36 @@+{-# LANGUAGE TemplateHaskell #-}+import Data.Integer.SAT++import Test.QuickCheck+import System.Exit++instance Arbitrary BoundType where+ arbitrary = elements [Lower, Upper]++withBounds :: Testable prop =>+ BoundType -> [(Positive Integer, Integer)] -> (Integer -> prop) -> Property+withBounds kind bs prop =+ counterexample (show (map toBound bs)) $+ case iPickBounded kind (map toBound bs) of+ Nothing -> property Discard+ Just n -> counterexample (show n) (property (prop n))+ where+ toBound (Positive c, t) = Bound c (tConst t)++prop_lower, prop_upper :: [(Positive Integer, Integer)] -> Property+prop_lower bs =+ withBounds Lower bs $ \n ->+ and [t < c * n | (Positive c, t) <- bs] &&+ or [t >= c * (n-1) | (Positive c, t) <- bs]+prop_upper bs =+ withBounds Upper bs $ \n ->+ and [c * n < t | (Positive c, t) <- bs] &&+ or [c * (n+1) >= t | (Positive c, t) <- bs]++-- This is so that the Template Haskell below can see the above properties.+$(return [])++main :: IO ()+main = do ok <- $(quickCheckAll)+ if ok then exitSuccess else exitFailure+