ghc-typelits-presburger 0.1.1.1 → 0.2.0.0
raw patch · 2 files changed
+1030/−49 lines, 2 filesdep +containersdep +prettydep −presburgerdep ~basedep ~equational-reasoningdep ~ghcPVP ok
version bump matches the API change (PVP)
Dependencies added: containers, pretty
Dependencies removed: presburger
Dependency ranges changed: base, equational-reasoning, ghc, ghc-typelits-presburger, reflection, singletons
API changes (from Hackage documentation)
Files
- ghc-typelits-presburger.cabal +49/−49
- presburger/src/Data/Integer/SAT.hs +981/−0
ghc-typelits-presburger.cabal view
@@ -1,56 +1,56 @@-name: ghc-typelits-presburger-version: 0.1.1.1-cabal-version: >=1.10-build-type: Simple-license: BSD3-license-file: LICENSE-copyright: 2015 (c) Hiromi ISHII-maintainer: konn.jinro _at_ gmail.com-homepage: https://github.com/konn/ghc-typelits-presburger#readme-synopsis: Presburger Arithmetic Solver for GHC Type-level natural numbers.-description:- @ghc-typelits-presburger@ augments GHC type-system with Presburger Arithmetic Solver for Type-level natural numbers.- You can use by adding this package to @build-depends@ and add the following pragma to the head of .hs files:- .- > OPTIONS_GHC -fplugin GHC.TypeLits.Presburger-category: Math, Type System-author: Hiromi ISHII-tested-with: GHC ==8.0.2 GHC ==8.2.2--source-repository head- type: git- location: https://github.com/konn/ghc-typelits-presburger+name: ghc-typelits-presburger+version: 0.2.0.0+synopsis: Presburger Arithmetic Solver for GHC Type-level natural numbers.+description: + @ghc-typelits-presburger@ augments GHC type-system with Presburger Arithmetic Solver for Type-level natural numbers.+ You can use by adding this package to @build-depends@ and add the following pragma to the head of .hs files:+ .+ > OPTIONS_GHC -fplugin GHC.TypeLits.Presburger+ +homepage: https://github.com/konn/ghc-typelits-presburger#readme+license: BSD3+license-file: LICENSE+author: Hiromi ISHII+maintainer: konn.jinro _at_ gmail.com+copyright: 2015 (c) Hiromi ISHII+category: Math, Type System+build-type: Simple+-- extra-source-files:+cabal-version: >=1.10+tested-with: GHC == 8.0.2, GHC == 8.2.2 flag examples- default: False+ Default: False library- exposed-modules:- GHC.TypeLits.Presburger- build-depends:- base >=4.7 && <5,- ghc >=7.10 && <8.4,- ghc-tcplugins-extra ==0.2.*,- presburger ==1.3.*,- equational-reasoning >=0.4.0.0 && <0.6,- reflection >=2.1.2 && <2.2- default-language: Haskell2010- hs-source-dirs: src- other-modules:- GHC.Compat- ghc-options: -Wall+ ghc-options: -Wall+ hs-source-dirs: src, presburger/src+ exposed-modules: GHC.TypeLits.Presburger+ other-modules: GHC.Compat+ Data.Integer.SAT+ build-depends: base >= 4.7 && < 5+ , ghc >= 7.10 && < 8.5+ , ghc-tcplugins-extra >= 0.2 && < 0.3+ , equational-reasoning >= 0.4.0.0+ , containers+ , pretty+ , reflection+ default-language: Haskell2010 executable simple-arith- - if !flag(examples)- buildable: False- main-is: simple-arith.hs- build-depends:- base >=4.10.1.0 && <4.11,- ghc-typelits-presburger >=0.1.1.1 && <0.2,- equational-reasoning >=0.5.0.0 && <0.6,- singletons >=2.3.1 && <2.4- default-language: Haskell2010- hs-source-dirs: examples- ghc-options: -Wall+ if flag(examples)+ buildable: True+ else+ buildable: False+ hs-source-dirs: examples+ Main-is: simple-arith.hs+ build-depends: base+ , ghc-typelits-presburger+ , equational-reasoning+ , singletons+ default-language: Haskell2010+ ghc-options: -Wall +source-repository head+ type: git+ location: https://github.com/konn/ghc-typelits-presburger
+ presburger/src/Data/Integer/SAT.hs view
@@ -0,0 +1,981 @@+{-# LANGUAGE BangPatterns, PatternGuards, Trustworthy #-}+{-|+This module implements a decision procedure for quantifier-free linear+arithmetic. The algorithm is based on the following paper:++ An Online Proof-Producing Decision Procedure for+ Mixed-Integer Linear Arithmetic+ by+ Sergey Berezin, Vijay Ganesh, and David L. Dill+-}+module Data.Integer.SAT+ ( PropSet+ , noProps+ , checkSat+ , assert+ , Prop(..)+ , Expr(..)+ , BoundType(..)+ , getExprBound+ , getExprRange+ , Name+ , toName+ , fromName+ -- * Iterators+ , allSolutions+ , slnCurrent+ , slnNextVal+ , slnNextVar+ , slnEnumerate+++ -- * Debug+ , dotPropSet+ , sizePropSet+ , allInerts+ , ppInerts++ -- * For QuickCheck+ , iPickBounded+ , Bound(..)+ , tConst+ ) where++import Debug.Trace++import Control.Applicative (Alternative (..), Applicative (..), (<$>))+import Control.Monad (MonadPlus (..), ap, guard, liftM)+import Data.List (partition)+import Data.Map (Map)+import qualified Data.Map as Map+import Data.Maybe (fromMaybe, mapMaybe, maybeToList)+import Prelude hiding ((<>))+import Text.PrettyPrint++infixr 2 :||+infixr 3 :&&+infix 4 :==, :/=, :<, :<=, :>, :>=+infixl 6 :+, :-+infixl 7 :*++--------------------------------------------------------------------------------+-- Solver interface++-- | A collection of propositions.+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++-- | Add a new proposition to an existing collection.+assert :: Prop -> PropSet -> PropSet+assert p (State rws) = State $ fmap snd $ m =<< rws+ where S m = prop p++-- | Extract a model from a consistent set of propositions.+-- Returns 'Nothing' if the assertions have no model.+-- If a variable does not appear in the assignment, then it is 0 (?).+checkSat :: PropSet -> Maybe [(Int,Integer)]+checkSat (State m) = go m+ where+ go None = mzero+ 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+getExprBound bt e (State s) =+ do let S m = expr e+ check (t,s1) = iTermBound bt t (inerts s1)+ bs <- mapM check $ toList $ s >>= m+ case bs of+ [] -> Nothing+ _ -> Just (maximum bs)++-- | Compute the range of possible values for an expression.+-- Returns `Nothing` if the bound is not known.+getExprRange :: Expr -> PropSet -> Maybe [Integer]+getExprRange e (State s) =+ do let S m = expr e+ check (t,s1) = do l <- iTermBound Lower t (inerts s1)+ u <- iTermBound Upper t (inerts s1)+ return (l,u)+ bs <- mapM check $ toList $ s >>= m+ case bs of+ [] -> Nothing+ _ -> let (ls,us) = unzip bs+ in Just [ x | x <- [ minimum ls .. maximum us ] ]++++-- | The type of proposition.+data Prop = PTrue+ | PFalse+ | Prop :|| Prop+ | Prop :&& Prop+ | Not Prop+ | Expr :== Expr+ | Expr :/= Expr+ | Expr :< Expr+ | Expr :> Expr+ | Expr :<= Expr+ | Expr :>= Expr+ deriving (Read,Show)++-- | The type of integer expressions.+-- Variable names must be non-negative.+data Expr = Expr :+ Expr -- ^ Addition+ | Expr :- Expr -- ^ Subtraction+ | Integer :* Expr -- ^ Multiplication by a constant+ | Negate Expr -- ^ Negation+ | Var Name -- ^ Variable+ | K Integer -- ^ Constant+ | If Prop Expr Expr -- ^ A conditional expression+ | Div Expr Integer -- ^ Division, rounds down+ | Mod Expr Integer -- ^ Non-negative remainder+ deriving (Read,Show)++prop :: Prop -> S ()+prop PTrue = return ()+prop PFalse = mzero+prop (p1 :|| p2) = prop p1 `mplus` prop p2+prop (p1 :&& p2) = prop p1 >> prop p2+prop (Not p) = prop (neg p)+ where+ neg PTrue = PFalse+ neg PFalse = PTrue+ neg (p1 :&& p2) = neg p1 :|| neg p2+ neg (p1 :|| p2) = neg p1 :&& neg p2+ neg (Not q) = q+ neg (e1 :== e2) = e1 :/= e2+ neg (e1 :/= e2) = e1 :== e2+ neg (e1 :< e2) = e1 :>= e2+ neg (e1 :<= e2) = e1 :> e2+ neg (e1 :> e2) = e1 :<= e2+ neg (e1 :>= e2) = e1 :< e2++prop (e1 :== e2) = do t1 <- expr e1+ t2 <- expr e2+ solveIs0 (t1 |-| t2)++prop (e1 :/= e2) = do t1 <- expr e1+ t2 <- expr e2+ let t = t1 |-| t2+ solveIsNeg t `orElse` solveIsNeg (tNeg t)++prop (e1 :< e2) = do t1 <- expr e1+ t2 <- expr e2+ solveIsNeg (t1 |-| t2)++prop (e1 :<= e2) = do t1 <- expr e1+ t2 <- expr e2+ let t = t1 |-| t2 |-| tConst 1+ solveIsNeg t++prop (e1 :> e2) = prop (e2 :< e1)+prop (e1 :>= e2) = prop (e2 :<= e1)+++expr :: Expr -> S Term+expr (e1 :+ e2) = (|+|) <$> expr e1 <*> expr e2+expr (e1 :- e2) = (|-|) <$> expr e1 <*> expr e2+expr (k :* e2) = (k |*|) <$> expr e2+expr (Negate e) = tNeg <$> expr e+expr (Var x) = pure (tVar x)+expr (K x) = pure (tConst x)+expr (If p e1 e2) = do x <- newVar+ prop (p :&& Var x :== e1 :|| Not p :&& Var x :== e2)+ return (tVar x)+expr (Div e k) = fmap fst $ exprDivMod e k+expr (Mod e k) = fmap snd $ exprDivMod e k++exprDivMod :: Expr -> Integer -> S (Term,Term)+exprDivMod e k =+ do guard (k /= 0) -- Always unsat+ q <- newVar+ r <- newVar+ let er = Var r+ prop (k :* Var q :+ er :== e :&& er :< K k :&& K 0 :<= er)+ return (tVar q, tVar r)++++++--------------------------------------------------------------------------------++data RW = RW { nameSource :: !Int+ , inerts :: Inerts+ } deriving Show++initRW :: RW+initRW = RW { nameSource = 0, inerts = iNone }++--------------------------------------------------------------------------------+-- Constraints and Bound on Variables++ctLt :: Term -> Term -> Term+ctLt t1 t2 = t1 |-| t2++ctEq :: Term -> Term -> Term+ctEq t1 t2 = t1 |-| t2++data Bound = Bound Integer Term -- ^ The integer is strictly positive+ deriving Show++data BoundType = Lower | Upper+ deriving Show++toCt :: BoundType -> Name -> Bound -> Term+toCt Lower x (Bound c t) = ctLt t (c |*| tVar x)+toCt Upper x (Bound c t) = ctLt (c |*| tVar x) t++++--------------------------------------------------------------------------------+-- Inert set++-- | The inert contains the solver state on one possible path.+data Inerts = Inerts+ { bounds :: NameMap ([Bound],[Bound])+ -- ^ Known lower and upper bounds for variables.+ -- Each bound @(c,t)@ in the first list asserts that @t < c * x@+ -- Each bound @(c,t)@ in the second list asserts that @c * x < t@++ , solved :: NameMap Term+ -- ^ Definitions for resolved variables.+ -- 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+ , solved = Map.empty+ }++-- | Rewrite a term using the definitions from an inert set.+iApSubst :: Inerts -> Term -> Term+iApSubst i t = foldr apS t $ Map.toList $ solved i+ where apS (x,t1) t2 = tLet x t1 t2++-- | Add a definition. Upper and lower bound constraints that mention+-- the variable are "kicked-out" so that they can be reinserted in the+-- context of the new knowledge.+--+-- * Assumes substitution has already been applied.+--+-- * The kicked-out constraints are NOT rewritten, this happens+-- when they get inserted in the work queue.++iSolved :: Name -> Term -> Inerts -> ([Term], Inerts)+iSolved x t i =+ ( kickedOut+ , Inerts { bounds = otherBounds+ , solved = Map.insert x t $ Map.map (tLet x t) $ solved i+ }+ )+ where+ (kickedOut, otherBounds) =++ -- First, we eliminate all entries for `x`+ let (mb, mp1) = Map.updateLookupWithKey (\_ _ -> Nothing) x (bounds i)++ -- Next, we elminate all constraints that mentiond `x` in bounds+ mp2 = Map.mapWithKey extractBounds mp1++ in ( [ ct | (lbs,ubs) <- maybeToList mb+ , ct <- map (toCt Lower x) lbs ++ map (toCt Upper x) ubs ]+ +++ [ ct | (_,cts) <- Map.elems mp2, ct <- cts ]++ , fmap fst mp2+ )++ extractBounds y (lbs,ubs) =+ let (lbsStay, lbsKick) = partition stay lbs+ (ubsStay, ubsKick) = partition stay ubs+ in ( (lbsStay,ubsStay)+ , map (toCt Lower y) lbsKick +++ map (toCt Upper y) ubsKick+ )++ 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+iPickBounded _ [] = Nothing+iPickBounded bt bs =+ do xs <- mapM (normBound bt) bs+ return $ case bt of+ 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 (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-1) c)+++-- | The largest (resp. least) upper (resp. lower) bound on a term+-- that will satisfy the model+iTermBound :: BoundType -> Term -> Inerts -> Maybe Integer+iTermBound bt (T k xs) is = do ks <- mapM summand (Map.toList xs)+ return $ sum $ k : ks+ where+ summand (x,c) = fmap (c *) (iVarBound (newBt c) x is)+ newBt c = if c > 0 then bt else case bt of+ Lower -> Upper+ Upper -> Lower+++-- | The largest (resp. least) upper (resp. lower) bound on a variable+-- that will satisfy the model.+iVarBound :: BoundType -> Name -> Inerts -> Maybe Integer+iVarBound bt x is+ | Just t <- Map.lookup x (solved is) = iTermBound bt t is++iVarBound bt x is =+ do both <- Map.lookup x (bounds is)+ case mapMaybe fromBound (chooseBounds both) of+ [] -> Nothing+ bs -> return (combineBounds bs)+ where+ fromBound (Bound c t) = fmap (scaleBound c) (iTermBound bt t is)++ combineBounds = case bt of+ Upper -> minimum+ Lower -> maximum++ chooseBounds = case bt of+ Upper -> snd+ Lower -> fst++ scaleBound c b = case bt of+ Upper -> div (b-1) c+ Lower -> div b c + 1+++++iModel :: Inerts -> [(Name,Integer)]+iModel i = goBounds [] (bounds i)+ where+ goBounds su mp =+ case Map.maxViewWithKey mp of+ Nothing -> goEqs su $ Map.toList $ solved i+ Just ((x,(lbs0,ubs0)), mp1) ->+ let lbs = [ Bound c (tLetNums su t) | Bound c t <- lbs0 ]+ ubs = [ Bound c (tLetNums su t) | Bound c t <- ubs0 ]+ sln = fromMaybe 0+ $ mplus (iPickBounded Lower lbs) (iPickBounded Upper ubs)+ in goBounds ((x,sln) : su) mp1++ goEqs su [] = su+ goEqs su ((x,t) : more) =+ let t1 = tLetNums su t+ vs = tVarList t1+ su1 = [ (v,0) | v <- vs ] ++ (x,tConstPart t1) : su+ in goEqs su1 more+++--------------------------------------------------------------------------------+-- 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++ -- A == 0+ | Just a <- isConst t = guard (a == 0)++ -- A + B * x = 0+ | Just (a,b,x) <- tIsOneVar t =+ case divMod (-a) b of+ (q,0) -> addDef x (tConst q)+ _ -> mzero++ -- 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 = solveIs0 s++ -- See Section 3.1 of paper for details.+ -- We obtain an equivalent formulation but with smaller coefficients.+ | Just (ak,xk,s) <- tLeastAbsCoeff t =+ do let m = abs ak + 1+ v <- newVar+ let sgn = signum ak+ soln = (negate sgn * m) |*| tVar v+ |+| tMapCoeff (\c -> sgn * modulus c m) s+ addDef xk soln++ let upd i = div (2*i + m) (2*m) + modulus i m+ solveIs0 (negate (abs ak) |*| tVar v |+| tMapCoeff upd s)++ | otherwise = error "solveIs0: unreachable"++modulus :: Integer -> Integer -> Integer+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++ -- A < 0+ | Just a <- isConst t = guard (a < 0)++ -- A * S < 0+ | Just (_,s) <- tFactor t = solveIsNeg s++ -- See Section 5.1 of the paper+ | Just (xc,x,s) <- tLeastVar t =++ do ctrs <- if xc < 0+ -- -XC*x + S < 0+ -- S < XC*x+ then do ubs <- getBounds Upper x+ let b = negate xc+ beta = s+ addBound Lower x (Bound b beta)+ return [ (a,alpha,b,beta) | Bound a alpha <- ubs ]+ -- XC*x + S < 0+ -- XC*x < -S+ else do lbs <- getBounds Lower x+ let a = xc+ alpha = tNeg s+ addBound Upper x (Bound a alpha)+ return [ (a,alpha,b,beta) | Bound b beta <- lbs ]++ -- See Note [Shadows]+ mapM_ (\(a,alpha,b,beta) ->+ do let real = ctLt (a |*| beta) (b |*| alpha)+ dark = ctLt (tConst (a * b)) (b |*| alpha |-| a |*| beta)+ gray = [ ctEq (b |*| tVar x) (tConst i |+| beta)+ | i <- [ 1 .. b - 1 ] ]+ 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]++ P: beta < b * x+ Q: a * x < alpha++real: a * beta < b * alpha++ beta < b * x -- from P+ a * beta < a * b * x -- (a *)+ a * beta < b * alpha -- comm. and Q+++dark: b * alpha - a * beta > a * b+++gray: b * x = beta + 1 \/+ b * x = beta + 2 \/+ ...+ b * x = beta + (b-1)++We stop at @b - 1@ because if:++> b * x >= beta + b+> a * b * x >= a * (beta + b) -- (a *)+> a * b * x >= a * beta + a * b -- distrib.+> b * alpha > a * beta + a * b -- comm. and Q+> b * alpha - a * beta > a * b -- subtract (a * beta)++which is covered by the dark shadow.+-}+++--------------------------------------------------------------------------------+-- Monads++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+ go (Choice xs ys) zs = go xs (go ys zs)+ go (One x) xs = x : xs+ go None xs = xs+++instance Monad Answer where+ return a = One a+ fail _ = None+ None >>= _ = None+ One a >>= k = k a+ Choice m1 m2 >>= k = mplus (m1 >>= k) (m2 >>= k)++instance Alternative Answer where+ empty = mzero+ (<|>) = mplus++instance MonadPlus Answer where+ mzero = None+ mplus None x = x+ -- mplus (Choice x y) z = mplus x (mplus y z)+ mplus x y = Choice x y++instance Functor Answer where+ fmap _ None = None+ fmap f (One x) = One (f x)+ fmap f (Choice x1 x2) = Choice (fmap f x1) (fmap f x2)++instance Applicative Answer where+ pure = return+ (<*>) = ap+++newtype S a = S (RW -> Answer (a,RW))++instance Monad S where+ return a = S $ \s -> return (a,s)+ S m >>= k = S $ \s -> do (a,s1) <- m s+ let S m1 = k a+ m1 s1++instance Alternative S where+ empty = mzero+ (<|>) = mplus++instance MonadPlus S where+ mzero = S $ \_ -> mzero+ mplus (S m1) (S m2) = S $ \s -> mplus (m1 s) (m2 s)++instance Functor S where+ fmap = liftM++instance Applicative S where+ pure = return+ (<*>) = ap++updS :: (RW -> (a,RW)) -> S a+updS f = S $ \s -> return (f s)++updS_ :: (RW -> RW) -> S ()+updS_ f = updS $ \rw -> ((),f rw)++get :: (RW -> a) -> S a+get f = updS $ \rw -> (f rw, rw)++newVar :: S Name+newVar = updS $ \rw -> ( SysName (nameSource rw)+ , rw { nameSource = nameSource rw + 1 }+ )++-- | 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+ Nothing -> []+ Just bs -> case f of+ Lower -> fst bs+ Upper -> snd bs++addBound :: BoundType -> Name -> Bound -> S ()+addBound bt x b = updS_ $ \rw ->+ let i = inerts rw+ entry = case bt of+ Lower -> ([b],[])+ Upper -> ([],[b])+ jn (newL,newU) (oldL,oldU) = (newL++oldL, newU++oldU)+ in rw { inerts = i { bounds = Map.insertWith jn x entry (bounds i) }}++-- | Add a new definition.+-- Assumes substitution has already been applied+addDef :: Name -> Term -> S ()+addDef x t =+ do newWork <- updS $ \rw -> let (newWork,newInerts) = iSolved x t (inerts rw)+ in (newWork, rw { inerts = newInerts })+ mapM_ solveIsNeg newWork++apSubst :: Term -> S Term+apSubst t =+ do i <- get inerts+ return (iApSubst i t)+++++--------------------------------------------------------------------------------+++data Name = UserName !Int | SysName !Int+ deriving (Read,Show,Eq,Ord)++ppName :: Name -> Doc+ppName (UserName x) = text "u" <> int x+ppName (SysName x) = text "s" <> int x++toName :: Int -> Name+toName = UserName++fromName :: Name -> Maybe Int+fromName (UserName x) = Just x+fromName (SysName _) = Nothing+++++type NameMap = Map Name++-- | The type of terms. The integer is the constant part of the term,+-- and the `Map` maps variables (represented by @Int@ to their coefficients).+-- The term is a sum of its parts.+-- INVARIANT: the `Map` does not map anything to 0.+data Term = T !Integer (NameMap Integer)+ deriving (Eq,Ord)++infixl 6 |+|, |-|+infixr 7 |*|++-- | A constant term.+tConst :: Integer -> Term+tConst k = T k Map.empty++-- | Construct a term with a single variable.+tVar :: Name -> Term+tVar x = T 0 (Map.singleton x 1)++(|+|) :: Term -> Term -> Term+T n1 m1 |+| T n2 m2 = T (n1 + n2)+ $ if Map.null m1 then m2 else+ if Map.null m2 then m1 else+ Map.filter (/= 0) $ Map.unionWith (+) m1 m2++(|*|) :: Integer -> Term -> Term+0 |*| _ = tConst 0+1 |*| t = t+k |*| T n m = T (k * n) (fmap (k *) m)++tNeg :: Term -> Term+tNeg t = (-1) |*| t++(|-|) :: Term -> Term -> Term+t1 |-| t2 = t1 |+| tNeg t2+++-- | Replace a variable with a term.+tLet :: Name -> Term -> Term -> Term+tLet x t1 t2 = let (a,t) = tSplitVar x t2+ in a |*| t1 |+| t++-- | Replace a variable with a constant.+tLetNum :: Name -> Integer -> Term -> Term+tLetNum x k t = let (c,T n m) = tSplitVar x t+ in T (c * k + n) m++-- | Replace the given variables with constants.+tLetNums :: [(Name,Integer)] -> Term -> Term+tLetNums xs t = foldr (\(x,i) t1 -> tLetNum x i t1) t xs+++++instance Show Term where+ showsPrec c t = showsPrec c (show (ppTerm t))++ppTerm :: Term -> Doc+ppTerm (T k m) =+ case Map.toList m of+ [] -> integer k+ xs | k /= 0 -> hsep (integer k : map ppProd xs)+ x : xs -> hsep (ppFst x : map ppProd xs)++ where+ ppFst (x,1) = ppName x+ ppFst (x,-1) = text "-" <> ppName x+ ppFst (x,n) = ppMul n x++ ppProd (x,1) = text "+" <+> ppName x+ ppProd (x,-1) = text "-" <+> ppName x+ ppProd (x,n) | n > 0 = text "+" <+> ppMul n x+ | otherwise = text "-" <+> ppMul (abs n) x++ ppMul n x = integer n <+> text "*" <+> ppName x++-- | Remove a variable from the term, and return its coefficient.+-- If the variable is not present in the term, the coefficient is 0.+tSplitVar :: Name -> Term -> (Integer, Term)+tSplitVar x t@(T n m) =+ case Map.updateLookupWithKey (\_ _ -> Nothing) x m of+ (Nothing,_) -> (0,t)+ (Just k,m1) -> (k, T n m1)++-- | Does the term contain this varibale?+tHasVar :: Name -> Term -> Bool+tHasVar x (T _ m) = Map.member x m++-- | Is this terms just an integer.+isConst :: Term -> Maybe Integer+isConst (T n m)+ | Map.null m = Just n+ | otherwise = Nothing++tConstPart :: Term -> Integer+tConstPart (T n _) = n++-- | Returns: @Just (a, b, x)@ if the term is the form: @a + b * x@+tIsOneVar :: Term -> Maybe (Integer, Integer, Name)+tIsOneVar (T a m) = case Map.toList m of+ [ (x,b) ] -> Just (a, b, x)+ _ -> Nothing++-- | Spots terms that contain variables with unit coefficients+-- (i.e., of the form @x + t@ or @t - x@).+-- Returns (coeff, var, rest of term)+tGetSimpleCoeff :: Term -> Maybe (Integer, Name, Term)+tGetSimpleCoeff (T a m) =+ do let (m1,m2) = Map.partition (\x -> x == 1 || x == -1) m+ ((x,xc), m3) <- Map.minViewWithKey m1+ return (xc, x, T a (Map.union m3 m2))++tVarList :: Term -> [Name]+tVarList (T _ m) = Map.keys m+++-- | Try to factor-out a common consant (> 1) from a term.+-- For example, @2 + 4x@ becomes @2 * (1 + 2x)@.+tFactor :: Term -> Maybe (Integer, Term)+tFactor (T c m) =+ do d <- common (c : Map.elems m)+ return (d, T (div c d) (fmap (`div` d) m))+ where+ common :: [Integer] -> Maybe Integer+ common [] = Nothing+ common [x] = Just x+ common (x : y : zs) =+ case gcd x y of+ 1 -> Nothing+ n -> common (n : zs)++-- | Extract a variable with a coefficient whose absolute value is minimal.+tLeastAbsCoeff :: Term -> Maybe (Integer, Name, Term)+tLeastAbsCoeff (T c m) = do (xc,x,m1) <- Map.foldWithKey step Nothing m+ return (xc, x, T c m1)+ where+ step x xc Nothing = Just (xc, x, Map.delete x m)+ step x xc (Just (yc,_,_))+ | abs xc < abs yc = Just (xc, x, Map.delete x m)+ step _ _ it = it++-- | Extract the least variable from a term+tLeastVar :: Term -> Maybe (Integer, Name, Term)+tLeastVar (T c m) =+ do ((x,xc), m1) <- Map.minViewWithKey m+ return (xc, x, T c m1)++-- | Apply a function to all coefficients, including the constnat+tMapCoeff :: (Integer -> Integer) -> Term -> Term+tMapCoeff f (T c m) = T (f c) (fmap f m)+++++++