presburger 0.4 → 1.0
raw patch · 12 files changed
+827/−1592 lines, 12 filesdep ~basedep ~containers
Dependency ranges changed: base, containers
Files
- presburger.cabal +5/−13
- src/Data/Integer/OldPresburger.hs +0/−673
- src/Data/Integer/Presburger.hs +0/−11
- src/Data/Integer/Presburger/Form.hs +0/−213
- src/Data/Integer/Presburger/HOAS.hs +0/−125
- src/Data/Integer/Presburger/ModArith.hs +0/−30
- src/Data/Integer/Presburger/Notation.hs +0/−47
- src/Data/Integer/Presburger/Prop.hs +0/−193
- src/Data/Integer/Presburger/SolveDiv.hs +0/−100
- src/Data/Integer/Presburger/Term.hs +0/−142
- src/Data/Integer/Presburger/Utils.hs +0/−45
- src/Data/Integer/SAT.hs +822/−0
presburger.cabal view
@@ -1,13 +1,14 @@ Name: presburger-Version: 0.4+Version: 1.0 License: BSD3 License-file: LICENSE Author: Iavor S. Diatchki Homepage: http://github.com/yav/presburger Maintainer: diatchki@galois.com Category: Algorithms-Synopsis: Cooper's decision procedure for Presburger arithmetic.-Description: Cooper's decision procedure for Presburger arithmetic.+Synopsis: A decision procedure for quantifier-free linear arithmetic.+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 @@ -15,16 +16,7 @@ Build-Depends: base < 10, containers, pretty hs-source-dirs: src Exposed-modules:- Data.Integer.Presburger- Data.Integer.OldPresburger- Data.Integer.Presburger.Term- Data.Integer.Presburger.Prop- Data.Integer.Presburger.Form- Data.Integer.Presburger.SolveDiv- Data.Integer.Presburger.Notation- Data.Integer.Presburger.HOAS- Data.Integer.Presburger.ModArith- Data.Integer.Presburger.Utils+ Data.Integer.SAT GHC-options: -O2 -Wall
− src/Data/Integer/OldPresburger.hs
@@ -1,673 +0,0 @@-{-| This module implements Cooper's algorithm for deciding- first order formulas over integers with addition.--Based on the paper:- * author: D.C.Cooper- * title: "Theorem Proving in Arithmetic without Multiplication"- * year: 1972--}-module Data.Integer.OldPresburger- ( check, simplify, Formula(..), Term, (.*), is_constant- , PP(..)- ) where---import qualified Data.IntMap as Map-import Data.Maybe(fromMaybe)-import Data.List(nub,foldl')-import Control.Monad(mplus,guard)-import Prelude hiding (LT,EQ)--import Text.PrettyPrint.HughesPJ----- | Check if a formula is true.-check :: Formula -> Bool-check f = eval_form (pre (True,0) f)--simplify :: Formula -> Formula-simplify f = invert (pre (True,0) f)---- Sugar --------------------------------------------------------------------------infixl 3 :/\:-infixl 2 :\/:-infixr 1 :=>:--infix 4 :<:, :<=:, :>:, :>=:, :=:, :/=:, :|----- Forst-oreder formulas for Presburger arithmetic.-data Formula = Formula :/\: Formula- | Formula :\/: Formula- | Formula :=>: Formula- | Not Formula- | Exists (Term -> Formula)- | Forall (Term -> Formula)- | TRUE- | FALSE- | Term :<: Term- | Term :>: Term- | Term :<=: Term- | Term :>=: Term- | Term :=: Term- | Term :/=: Term- | Integer :| Term--pre :: (Bool,Int) -> Formula -> Form-pre n form = case form of- f1 :/\: f2 -> and' (pre n f1) (pre n f2)- f1 :\/: f2 -> or' (pre n f1) (pre n f2)- f1 :=>: f2 -> pre n (Not f1 :\/: f2)- Exists f -> pre_ex (top,x + 1) [x] (f (var x))- where (top,x) = n- Forall f -> pre n (Not (Exists (Not . f)))- TRUE -> tt'- FALSE -> ff'- t1 :<: t2 -> lt' t1 t2- t1 :>: t2 -> lt' t2 t1- t1 :<=: t2 -> leq' t1 t2- t1 :>=: t2 -> leq' t2 t1- t1 :=: t2 -> eq' t1 t2- t1 :/=: t2 -> neq' t1 t2- k :| t -> divs' k t- Not form1 -> case form1 of- Not f -> pre n f- Forall f -> pre n (Exists (Not . f))- _ -> not' (pre n form1)--pre_ex :: (Bool,Int) -> [Name] -> Formula -> Form-pre_ex (top,n) xs form = case form of- Exists f -> pre_ex (top,n+1) (n:xs) (f (var n))- f1 :\/: f2 -> or' (pre_ex (top,n) xs f1) (pre_ex (top,n) xs f2)- Not form1 ->- case form1 of- Not form2 -> pre_ex (top,n) xs form2- Forall f -> pre_ex (top,n) xs (Exists (Not . f))- p :/\: q -> pre_ex (top,n) xs (Not p :\/: Not q)- _ -> exists_many top xs (pre (False,n) form)- _ -> exists_many top xs (pre (False,n) form)--invert :: Form -> Formula-invert form = case form of- Conn And f1 f2 -> invert f1 :/\: invert f2- Conn Or f1 f2 -> invert f1 :\/: invert f2- Prop prop -> case prop of- Pred FF True :> [] -> FALSE- Pred FF False :> [] -> TRUE- Pred LT True :> [t1,t2] -> t1 :<: t2- Pred LT False :> [t1,t2] -> t1 :>=: t2- Pred LEQ True :> [t1,t2] -> t1 :<=: t2- Pred LEQ False :> [t1,t2] -> t1 :>: t2- Pred EQ True :> [t1,t2] -> t1 :=: t2- Pred EQ False :> [t1,t2] -> t1 :/=: t2- Pred (Divs n) True :> [t] -> n :| t- Pred (Divs n) False :> [t] -> Not (n :| t)- _ -> error "(bug) Type error in 'invert'"----- Terms -------------------------------------------------------------------------- | Terms of Presburger arithmetic.--- Term are created by using the 'Num' class.--- WARNING: Presburger arithmetic only supports multiplication--- by a constant, trying to create invalid terms will result--- in a run-time error. A more type-safe alternative is to--- use the '(.*)' operator.-data Term = Term (Map.IntMap Integer) Integer---type Name = Int---- | @split_term x (n * x + t1) = (n,t1)@--- @x@ does not occur in @t1@-split_term :: Name -> Term -> (Integer,Term)-split_term x (Term m n) = (fromMaybe 0 c, Term m1 n)- where (c,m1) = Map.updateLookupWithKey (\_ _ -> Nothing) x m--var :: Name -> Term-var x = Term (Map.singleton x 1) 0--num :: Integer -> Term-num n = Term Map.empty n-------------------------------------------------------------------------------------instance Eq Term where- t1 == t2 = is_constant (t1 - t2) == Just 0--instance Num Term where- fromInteger n = Term Map.empty n-- Term m1 n1 + Term m2 n2 = Term (Map.unionWith (+) m1 m2) (n1 + n2)-- negate (Term m n) = Term (Map.map negate m) (negate n)-- t1 * t2 = case fmap (.* t2) (is_constant t1) `mplus`- fmap (.* t1) (is_constant t2) of- Just t -> t- Nothing -> error $ unlines [ "[(*) @ Term] Non-linear product:"- , " *** " ++ show t1- , " *** " ++ show t2- ]- signum t = case is_constant t of- Just n -> num (signum n)- Nothing -> error $ unlines [ "[signum @ Term]: Non-constant:"- , " *** " ++ show t- ]-- abs t = case is_constant t of- Just n -> num (abs n)- Nothing -> error $ unlines [ "[abs @ Term]: Non-constant:"- , " *** " ++ show t- ]----- | Check if a term is a constant (i.e., contains no variables).--- If so, then we return the constant, otherwise we return 'Nothing'.-is_constant :: Term -> Maybe Integer-is_constant (Term m n) = guard (all (0 ==) (Map.elems m)) >> return n--(.*) :: Integer -> Term -> Term-0 .* _ = 0-1 .* t = t-k .* Term m n = Term (Map.map (k *) m) (k * n)----- Formulas ----------------------------------------------------------------------data PredSym = FF | LT | LEQ | EQ | Divs Integer {- +ve -}-data Pred = Pred PredSym Bool -- Bool: positive (i.e. non-negated)?-data Prop = Pred :> [Term]-data Conn = And | Or deriving Eq-data Form = Conn Conn Form Form | Prop Prop--abs_form :: Form -> ([Prop],[Prop] -> Form)-abs_form fo = let (ps,skel) = loop [] fo- in (reverse ps, fst . skel)- where loop ps (Conn c p q) =- let (ps1,f1) = loop ps p- (ps2,f2) = loop ps1 q- in (ps2, \fs -> let (p1,fs1) = f1 fs- (p2,fs2) = f2 fs1- in (Conn c p1 p2, fs2))- loop ps (Prop p) = (p:ps, \(f:fs) -> (Prop f,fs))---not' :: Form -> Form-not' (Conn c t1 t2) = Conn (not_conn c) (not' t1) (not' t2)-not' (Prop p) = Prop (not_prop p)--ff' :: Form-ff' = Prop $ Pred FF True :>[]--tt' :: Form-tt' = Prop $ Pred FF False :>[]--lt' :: Term -> Term -> Form-lt' t1 t2 = Prop $ Pred LT True :> [t1,t2]--leq' :: Term -> Term -> Form-leq' t1 t2 = Prop $ Pred LEQ True :> [t1,t2]--eq' :: Term -> Term -> Form-eq' t1 t2 = Prop $ Pred EQ True :> [t1,t2]--neq' :: Term -> Term -> Form-neq' t1 t2 = Prop $ Pred EQ False :> [t1,t2]--and' :: Form -> Form -> Form-and' p q = Conn And p q--or' :: Form -> Form -> Form-or' p q = Conn Or p q--divs' :: Integer -> Term -> Form-divs' n t = Prop $ Pred (Divs n) True :> [t]--ors' :: [Form] -> Form-ors' [] = ff'-ors' xs = foldr1 or' xs--not_conn :: Conn -> Conn-not_conn And = Or-not_conn Or = And--not_prop :: Prop -> Prop-not_prop (f :> ts) = not_pred f :> ts--not_pred :: Pred -> Pred-not_pred (Pred p pos) = Pred p (not pos)------ Eliminating existential quantifiers -------------------------------------------data NormProp = Ind Prop- | L Pred Term--norm2 :: Name -> Integer -> Pred -> Term -> Term -> (Integer,NormProp)-norm2 x final_k p t1 t2- | k1 == k2 = (1, Ind (p :> [t1',t2']))- | k1 > k2 = (abs k, L p t)- | otherwise = (abs k, L p' t)-- where (k1,t1') = split_term x t1- (k2,t2') = split_term x t2-- k = k1 - k2- t = (final_k `div` k) .* (t2' - t1') -- only used when k /= 0-- p' = case p of- Pred LT b -> Pred LEQ (not b)- Pred LEQ b -> Pred LT (not b)- _ -> p--norm1 :: Name -> Integer -> Pred -> Term -> (Integer,NormProp)-norm1 x final_k p@(Pred (Divs d) b) t- | k == 0 = (1, Ind (p :> [t]))- | otherwise = (abs k, L ps (l .* t'))-- where (k,t') = split_term x t- l = final_k `div` k- ps = Pred (Divs (d * abs l)) b--norm1 _ _ _ _ = error "(bug) norm1 applied to a non-unary operator"---norm_prop :: Name -> Integer -> Prop -> (Integer,NormProp)-norm_prop _ _ p@(_ :> []) = (1,Ind p)-norm_prop x final_k (p :> [t]) = norm1 x final_k p t-norm_prop x final_k (p :> [t1,t2]) = norm2 x final_k p t1 t2-norm_prop _ _ _ = error "(bug) norm_prop on arity > 2"---- The integer is "length as - length bs"-a_b_sets :: (Integer,[Term],[Term]) -> NormProp -> (Integer,[Term],[Term])-a_b_sets (o,as,bs) p = case p of- Ind _ -> (o,as,bs)-- L (Pred op True) t ->- case op of- LT -> (1 + o , t : as, bs)- LEQ -> (1 + o , (t+1) : as, bs)- EQ -> (o , (t+1) : as, (t-1) : bs)- _ -> (o , as, bs)-- L (Pred op False) t ->- case op of- LT -> (o - 1 , as, (t-1) : bs)- LEQ -> (o - 1 , as, t : bs)- EQ -> (o , t : as, t : bs)- _ -> (o , as, bs)---analyze_props :: Name -> [Prop] -> ( [NormProp]- , Integer -- scale- , Integer -- bound- , Either [Term] [Term] -- A set or B set- )-analyze_props x ps = (ps1, final_k, bnd, if o < 0 then Left as else Right bs)- where (ks,ps1) = unzip $ map (norm_prop x final_k) ps- final_k = lcms ks- (o,as,bs) = foldl' a_b_sets (0,[],[]) ps1- bnd = lcms (final_k : [ d | L (Pred (Divs d) _) _ <- ps1 ])--from_bool :: Bool -> Prop-from_bool True = Pred FF False :> []-from_bool False = Pred FF True :> []--neg_inf :: NormProp -> Term -> Prop-neg_inf prop t = case prop of- Ind p -> p- L ps@(Pred op pos) t1 -> case op of- LT -> from_bool pos- LEQ -> from_bool pos- EQ -> from_bool (not pos)- Divs {} -> ps :> [t + t1]- FF -> error "(bug) FF in NormPred"--pos_inf :: NormProp -> Term -> Prop-pos_inf prop t = case prop of- Ind p -> p- L ps@(Pred op pos) t1 -> case op of- LT -> from_bool (not pos)- LEQ -> from_bool (not pos)- EQ -> from_bool (not pos)- Divs {} -> ps :> [t + t1]- FF -> error "(bug) FF in NormPred"--normal :: NormProp -> Term -> Prop-normal prop t = case prop of- Ind p -> p- L ps@(Pred (Divs {}) _) t1 -> ps :> [t + t1]- L ps t1 -> ps :> [t,t1]---data Ex = Ex [(Name,Integer)]- [Constraint]- [Prop]--exists_many :: Bool -> [Name] -> Form -> Form-exists_many top xs f = ors'- $ map exp_f- $ foldr (concatMap . ex_step) [Ex [] [] ps] (nub xs)- where (ps,skel) = abs_form f- exp_f = if top then expand_top skel else expand skel---ex_step :: Name -> Ex -> [Ex]-ex_step x (Ex xs ds ps) = case as_or_bs of- Left as ->- ( let arg = negate (var x)- in Ex ((x,d) : xs) (constr arg) (map (`pos_inf` arg) ps1)- ) : [ let arg = a - var x- in Ex ((x,d) : xs) (constr arg) (map (`normal` arg) ps1) | a <- as ]-- Right bs ->- ( let arg = var x- in Ex ((x,d) : xs) (constr arg) (map (`neg_inf` arg) ps1)- ) : [ let arg = b + var x- in Ex ((x,d) : xs) (constr arg) (map (`normal` arg) ps1) | b <- bs ]-- where (ps1,k,d',as_or_bs) = analyze_props x ps- d = lcms (d' : map fst ds)- constr t = if k == 1 then ds else (k,t) : ds---expand_top :: ([Prop] -> Form) -> Ex -> Form-expand_top skel (Ex xs ds ps) =- ors' [ skel (map (subst_prop env) ps) | env <- elim xs ds ]--expand :: ([Prop] -> Form) -> Ex -> Form-expand skel (Ex xs ds ps) =- ors' [ foldr and' (skel (map (subst_prop env) ps)) (map (`ctr` env) ds)- | env <- envs xs ]-- where envs [] = [ Map.empty ]- envs ((x,bnd):qs) = [ Map.insert x v env- | env <- envs qs, v <- [ 1 .. bnd ] ]-- ctr (k,t) env = Prop (Pred (Divs k) True :> [ subst_term env t ])----type Env = Map.IntMap Integer--subst_prop :: Env -> Prop -> Prop-subst_prop env (p :> ts) = p :> map (subst_term env) ts--subst_term :: Env -> Term -> Term-subst_term env (Term m n) =- let (xs,vs) = unzip $ Map.toList $ Map.intersectionWith (*) env m- in Term (foldl' (flip Map.delete) m xs) (foldl' (+) n vs)------- Evaluation ---------------------------------------------------------------------- The meanings of formulas.-eval_form :: Form -> Bool-eval_form (Conn c p q) = eval_conn c (eval_form p) (eval_form q)-eval_form (Prop p) = eval_prop p---- The meanings of connectives.-eval_conn :: Conn -> Bool -> Bool -> Bool-eval_conn And = (&&)-eval_conn Or = (||)---- The meanings of atomic propositions.-eval_prop :: Prop -> Bool-eval_prop (Pred p pos :> ts) = if pos then res else not res- where res = eval_pred p (map eval_term ts)---- The meanings of predicate symbols.-eval_pred :: PredSym -> [Integer] -> Bool-eval_pred p ts = case (p,ts) of- (FF, []) -> False- (Divs d, [k]) -> divides d k- (LT, [x,y]) -> x < y- (LEQ, [x,y]) -> x <= y- (EQ, [x,y]) -> x == y- _ -> error "Type error"---- We define: "d | a" as "exists y. d * y = a"-divides :: Integral a => a -> a -> Bool-0 `divides` 0 = True-0 `divides` _ = False-x `divides` y = mod y x == 0---- The meaning of a term with no free variables.--- NOTE: We do not check that there are no free variables.-eval_term :: Term -> Integer-eval_term (Term _ k) = k---- The meaning of a term with free variables-eval_term_env :: Term -> Env -> Integer-eval_term_env (Term m k) env = sum (k : map eval_var (Map.toList m))- where eval_var (x,c) = case Map.lookup x env of- Nothing -> error "free var"- Just v -> c * v-------------------------------------------------------------------------------------- Solving divides constraints ---------------------------------------------------- See the paper's appendix.----- | let (p,q,r) = extended_gcd x y--- in (x * p + y * q = r) && (gcd x y = r)-extended_gcd :: Integral a => a -> a -> (a,a,a)-extended_gcd arg1 arg2 = loop arg1 arg2 0 1 1 0- where loop a b x lastx y lasty- | b /= 0 = let (q,b') = divMod a b- x' = lastx - q * x- y' = lasty - q * y- in x' `seq` y' `seq` loop b b' x' x y' y- | otherwise = (lastx,lasty,a)---type Constraint = (Integer,Term)-type VarConstraint = (Integer,Integer,Term)---- m | (x * a1 + b1) /\ (n | x * a2 + b2)-theorem1 :: VarConstraint -> VarConstraint -> (VarConstraint, Constraint)-theorem1 (m,a1,b1) (n,a2,b2) = (new_x, new_other)- where new_x = (m * n, d, (p*n) .* b1 + (q * m) .* b2)- new_other = (d, a2 .* b1 - a1 .* b2)-- (p,q,d) = extended_gcd (a1 * n) (a2 * m)---- solutions for x in [1 .. bnd] of: m | x * a + b-theorem2 :: Integer -> (Integer,Integer,Integer) -> [Integer]-theorem2 bnd (m,a,b)- | r == 0 = [ t * k - c | t <- [ lower .. upper ] ]- | otherwise = []- where k = div m d- c = p * qu- (p,_,d) = extended_gcd a m- (qu,r) = divMod b d-- (lower',r1) = divMod (1 + c) k- lower = if r1 == 0 then lower' else lower' + 1 -- hmm- upper = div (bnd + c) k-- -- lower and upper:- -- t * k - c = 1 --> t = (1 + c) / k- -- t * k - c = bnd --> t = (bnd + c) / k-----elim :: [(Name,Integer)] -> [Constraint] -> [ Env ]-elim [] ts = if all chk ts then [ Map.empty ] else []- where chk (x,t) = divides x (eval_term t)-elim ((x,bnd):xs) cs = do env <- elim xs cs1- v <- case mb of- Nothing -> [ 1 .. bnd ]- Just (a,b,t) ->- theorem2 bnd (a,b,eval_term_env t env)- return (Map.insert x v env)-- where (mb,cs1) = elim_var x cs-----elim_var :: Name -> [Constraint] -> (Maybe VarConstraint, [Constraint])-elim_var x cs = case foldl' part ([],[]) cs of- ([], have_not) -> (Nothing, have_not)- (h : hs, have_not) -> let (c,hn) = step h hs have_not- in (Just c,hn)- where part s@(have,have_not) c@(m,t)- | m == 1 = s- | a == 0 = (have , c:have_not)- | otherwise = ((m,a,b):have, have_not)- where (a,b) = split_term x t-- step :: VarConstraint -> [VarConstraint] -> [Constraint]- -> (VarConstraint,[Constraint])- step h [] ns = (h,ns)- step h (h1:hs) ns = step h2 hs (n : ns)- where (h2,n) = theorem1 h h1---- Misc -------------------------------------------------------------------------lcms :: Integral a => [a] -> a-lcms xs = foldr lcm 1 xs----- Pretty Printing ---------------------------------------------------------------class PP a where- pp :: a -> Doc---var_name :: Name -> String-var_name x = let (a,b) = divMod x 26- rest = if a == 0 then "" else show a- in toEnum (97 + b) : rest--instance Show Term where show x = show (pp x)-instance PP Term where- pp (Term m k) | isEmpty vars = text (show k)- | k == 0 = vars- | k > 0 = vars <+> char '+' <+> text (show k)- | otherwise = vars <+> char '-' <+> text (show $ abs k)- where ppvar (x,n) = sign <+> co <+> text (var_name x)- where (sign,co)- | n == -1 = (char '-', empty)- | n < 0 = (char '-', text (show (abs n)) <+> char '*')- | n == 1 = (char '+', empty)- | otherwise = (char '+', text (show n) <+> char '*')- first_var (x,1) = text (var_name x)- first_var (x,-1) = char '-' <> text (var_name x)- first_var (x,n) = text (show n) <+> char '*' <+> text (var_name x)-- vars = case filter ((/= 0) . snd) (Map.toList m) of- [] -> empty- v : vs -> first_var v <+> hsep (map ppvar vs)----- 4: wrap term, not--- 3: wrap and--- 2: wrap or--- 1: wrap implies, quantifiers-instance PP Formula where- pp = pp1 0 -- ' 0 0- where- pp1 :: Int -> Formula -> Doc- pp1 p form = case form of- _ :/\: _ -> hang (text "/\\") 2 (loop form)- where loop (f1 :/\: f2) = loop f1 $$ loop f2- loop f = pp f-- _ :\/: _ -> hang (text "\\/") 2 (loop form)- where loop (f1 :\/: f2) = loop f1 $$ loop f2- loop f = pp f-- _ -> pp' 0 p form---- pp' :: Int -> Name -> Formula -> Doc- pp' n p form = case form of- f1 :/\: f2 | n < 3 -> pp' 2 p f1 <+> text "/\\" <+> pp' 2 p f2- f1 :\/: f2 | n < 2 -> pp' 1 p f1 <+> text "\\/" <+> pp' 1 p f2- f1 :=>: f2 | n < 1 -> pp' 1 p f1 <+> text "=>" <+> pp' 0 p f2- Not f | n < 4 -> text "Not" <+> pp' 4 p f- Exists {} | n < 1 -> pp_ex (text "exists") p form- where pp_ex d q (Exists g) = pp_ex (d <+> text (var_name q))- (q+1) (g (var q))- pp_ex d q g = d <> text "." <+> pp' 0 q g-- Forall {} | n < 1 -> pp_ex (text "forall") p form- where pp_ex d q (Forall g) = pp_ex (d <+> text (var_name q))- (q+1) (g (var q))- pp_ex d q g = d <> text "." <+> pp' 0 q g- TRUE -> text "true"- FALSE -> text "false"- t1 :<: t2 | n < 4 -> pp t1 <+> text "<" <+> pp t2- t1 :>: t2 | n < 4 -> pp t1 <+> text ">" <+> pp t2- t1 :<=: t2 | n < 4 -> pp t1 <+> text "<=" <+> pp t2- t1 :>=: t2 | n < 4 -> pp t1 <+> text ">=" <+> pp t2- t1 :=: t2 | n < 4 -> pp t1 <+> text "=" <+> pp t2- t1 :/=: t2 | n < 4 -> pp t1 <+> text "/=" <+> pp t2- k :| t1 | n < 4 -> text (show k) <+> text "|" <+> pp t1- _ -> parens (pp' 0 p form)--instance Show Formula where show = show . pp----instance PP PredSym where- pp p = case p of- FF -> text "false"- LT -> text "<"- LEQ -> text "<="- EQ -> text "==="- Divs n -> text (show n) <+> text "|"--instance PP Pred where- pp (Pred p True) = pp p- pp (Pred p False) = case p of- FF -> text "true"- LT -> text ">="- LEQ -> text ">"- EQ -> text "=/="- Divs n -> text (show n) <+> text "/|"--instance Show Prop where show = show . pp-instance PP Prop where- pp (p :> [t1,t2]) = pp t1 <+> pp p <+> pp t2- pp (p :> ts) = pp p <+> hsep (map pp ts)---instance PP Conn where- pp And = text "/\\"- pp Or = text "\\/"--instance PP Form where- pp me@(Conn c _ _) = hang (pp c) 2 (vcat $ map pp $ jn me [])- where jn (Conn c1 p1 q1) fs | c == c1 = jn p1 (jn q1 fs)- jn f fs = f : fs- pp (Prop p) = pp p--instance PP NormProp where- pp (Ind p) = pp p- pp (L p@(Pred (Divs {}) _) t) = pp p <+> text "_ +" <+> pp t- pp (L p t) = text "_" <+> pp p <+> pp t--instance Show NormProp where show = show . pp--instance PP Ex where- pp (Ex xs ps ss) = hang (text "OR" <+> hsep (map quant xs)) 2- ( text "!" <+> hsep (map (parens . divs) ps)- $$ vcat (map pp ss)- )- where quant (x,n) = parens $ text (var_name x) <> colon <> text (show n)- divs (x,t) = text (show x) <+> text "|" <+> pp t--
− src/Data/Integer/Presburger.hs
@@ -1,11 +0,0 @@-{-| This module implements Cooper's algorithm for deciding- first order formulas over integers with addition.--Based on the paper:- * author: D.C.Cooper- * title: "Theorem Proving in Arithmetic without Multiplication"- * year: 1972--}-module Data.Integer.Presburger (module X) where- -import Data.Integer.Presburger.HOAS as X
− src/Data/Integer/Presburger/Form.hs
@@ -1,213 +0,0 @@-module Data.Integer.Presburger.Form- ( module Data.Integer.Presburger.Form- , module Data.Integer.Presburger.Prop- ) where--import Data.Integer.Presburger.Prop-import Data.Integer.Presburger.SolveDiv--check :: Form (Prop PosP) -> Bool-check f = eval_form f env_empty---data Conn = And | Or deriving Eq-data Form p = Node !Conn (Form p) (Form p)- | Leaf !p-- -- A special form of disjunction. Bool = negated?- | Ex Bool (Name,Integer) (Form p)--instance Functor Form where- fmap f (Node c f1 f2) = Node c (fmap f f1) (fmap f f2)- fmap f (Ex b xs g) = Ex b xs (fmap f g)- fmap f (Leaf p) = Leaf (f p)--form_lcm :: Form (NormProp CVarP) -> Integer-form_lcm (Node _ f1 f2) = lcm (form_lcm f1) (form_lcm f2)-form_lcm (Leaf p) = case p of- Ind {} -> 1- Norm p1 -> coeff (prop p1)-form_lcm (Ex _ _ f) = form_lcm f----form_scale :: Name -> Form (Prop PosP) -> Form (NormProp VarP)-form_scale x form- | k /= 1 = Node And (Leaf $ Norm $ Prop False $ NDivides k 0) sf- | otherwise = sf- where- nf = fmap (norm x) form- k = form_lcm nf- sf = fmap leaf nf-- leaf p = case p of- Ind p1 -> Ind p1- Norm p1 -> Norm (scale k p1)----- The integer is "length as - length bs"-a_b_sets :: (Integer,[Term],[Term]) -> NormProp VarP -> (Integer,[Term],[Term])-a_b_sets (o,as,bs) p = case p of- Ind _ -> (o,as,bs)- Norm (Prop _ (NDivides {})) -> (o,as,bs)-- -- positive- Norm (Prop False (NBin op t)) ->- case op of- LessThan -> (1 + o , t : as, bs)- LessThanEqual -> (1 + o , (t+1) : as, bs)- Equal -> (o , (t+1) : as, (t-1) : bs)-- -- negative- Norm (Prop True (NBin op t)) ->- case op of- LessThan -> (o - 1 , as, (t-1) : bs)- LessThanEqual -> (o - 1 , as, t : bs)- Equal -> (o , t : as, t : bs)---form_pos_inf :: Term -> Form (NormProp VarP) -> Form (Prop PosP)-form_pos_inf t form = fmap leaf form- where leaf p = case p of- Ind p1 -> p1- Norm p1 -> pos_inf t p1--form_neg_inf :: Term -> Form (NormProp VarP) -> Form (Prop PosP)-form_neg_inf t form = fmap leaf form- where leaf p = case p of- Ind p1 -> p1- Norm p1 -> neg_inf t p1--form_no_inf :: Term -> Form (NormProp VarP) -> Form (Prop PosP)-form_no_inf t form = fmap leaf form- where leaf p = case p of- Ind p1 -> p1- Norm p1 -> normal t p1---neg :: Form (Prop PosP) -> Form (Prop PosP)-neg (Node And f1 f2) = Node Or (neg f1) (neg f2)-neg (Node Or f1 f2) = Node And (neg f1) (neg f2)-neg (Ex b x f) = Ex (not b) x f-neg (Leaf (Prop b p)) = Leaf (Prop (not b) p)---simplify :: Form (Prop PosP) -> Form (Prop PosP)-simplify (Node c f1 f2) =- case simplify f1 of- r@(Leaf (Prop n FF)) | n && c == Or- || not n && c == And -> r- | otherwise -> simplify f2- r1 -> case simplify f2 of- r@(Leaf (Prop n FF)) | n && c == Or- || not n && c == And -> r- | otherwise -> r1- r2 -> Node c r1 r2----simplify (Ex False (x,1) f) = simplify (subst_form x 1 f)-simplify (Ex True (x,1) f) = simplify (neg (subst_form x 1 f))--simplify (Ex b x f) = case simplify f of- Leaf (Prop n FF) -> Leaf (Prop (not (b == n)) FF)- f1 -> Ex b x f1- -simplify (Leaf l) = Leaf (simplify_prop l)----ex_step :: Name -> Form (Prop PosP) -> Form (Prop PosP)-ex_step x (Node Or f1 f2) = Node Or (ex_step x f1) (ex_step x f2)-ex_step x f- | as_minus_bs >= 0 = thm_as as x norm_f- | otherwise = thm_bs bs x norm_f- - where - norm_f :: Form (NormProp VarP)- norm_f = form_scale x f-- (as_minus_bs, as, bs) = loop (0,[],[]) norm_f-- loop s (Node _ f1 f2) = loop (loop s f1) f2- loop s (Ex _ _ f1) = loop s f1- loop s (Leaf p) = a_b_sets s p----thm_as :: [Term] -> Name -> Form (NormProp VarP) -> Form (Prop PosP)-thm_as as x f = simplify $- foldr1 (Node Or) $ map (Ex False (x, form_bound f))- $ form_pos_inf (negate (var x)) f- : [ form_no_inf (a - var x) f | a <- as ]--thm_bs :: [Term] -> Name -> Form (NormProp VarP) -> Form (Prop PosP)-thm_bs bs x f = simplify $- foldr1 (Node Or) $ map (Ex False (x, form_bound f))- $ form_neg_inf (var x) f- : [ form_no_inf (b + var x) f | b <- bs ]---form_bound :: Form (NormProp VarP) -> Integer-form_bound (Node _ f1 f2) = lcm (form_bound f1) (form_bound f2)-form_bound (Leaf p) = case p of- Norm (Prop _ (NDivides n _)) -> n- _ -> 1-form_bound (Ex _ _ f) = form_bound f----- Evaluation ---------------------------------------------------------------------- The meanings of formulas.-eval_form :: Form (Prop PosP) -> Env -> Bool-eval_form (Node c p q) env = eval_conn c (eval_form p env) (eval_form q env)-eval_form (Leaf p) env = eval_prop p env-eval_form (Ex b x f) env0 =- case splt f [x] of- (xs,cs,f1) -> let v = any (eval_form f1) (elim env0 xs cs)- in if b then not v else v- where splt (Ex False y f1) ys = splt f1 (y:ys)- splt f1 ys = case split_divs f1 of- (ds,f2) -> (ys,ds,f2)- --split_ands :: Form p -> [Form p]-split_ands form = loop form []- where loop (Node And f1 f2) fs = loop f1 (loop f2 fs)- loop f fs = f : fs--split_divs :: Form (Prop PosP) -> ([DivCtr], Form (Prop PosP))-split_divs form = loop (split_ands form) ([], Leaf (Prop True FF))- where- loop (Leaf (Prop False (Divides n t)) : fs) (cs, f)- = loop fs (Divs n t : cs, f)- loop (f:fs) (cs, f1) = loop fs (cs, Node And f f1)- loop [] cs = cs----- The meanings of connectives.-eval_conn :: Conn -> Bool -> Bool -> Bool-eval_conn And = (&&)-eval_conn Or = (||)--subst_form :: Name -> Integer -> Form (Prop PosP) -> Form (Prop PosP)-subst_form x n f = fmap (subst_prop x n) f-----------------------------------------------------------------------------------instance PP Conn where- pp And = text "/\\"- pp Or = text "\\/"--instance PP p => PP (Form p) where- pp me@(Node c _ _) = hang (pp c) 2 (vcat $ map pp $ jn me [])- where jn (Node c1 p1 q1) fs | c == c1 = jn p1 (jn q1 fs)- jn f fs = f : fs- pp (Leaf p) = pp p-- pp (Ex n q f) = hang (how <+> quant q <> text ".") 2 (pp f)- where quant (x,b) = text (var_name x) <+> text "<=" <+> text (show b)- how = (if n then text "Not" else empty) <+> text "Ex"----
− src/Data/Integer/Presburger/HOAS.hs
@@ -1,125 +0,0 @@-{-# LANGUAGE FlexibleInstances #-} - -module Data.Integer.Presburger.HOAS - ( Formula(..), check, translate - , Quant, exists, forall - , Term, (.*), is_constant - , PP(..) - ) where - -import Data.Integer.Presburger.Form hiding (check) -import qualified Data.Integer.Presburger.Form as F - -check :: Formula -> Bool -check f = F.check (translate f) - - -infixl 3 :/\: -infixl 2 :\/: -infixr 1 :=>: -infix 0 :<=>: - -infix 4 :<:, :<=:, :>:, :>=:, :=:, :/=:, :| - --- Forst-oreder formulas for Presburger arithmetic. -data Formula = Formula :/\: Formula - | Formula :\/: Formula - | Formula :=>: Formula - | Formula :<=>: Formula - | Not Formula - | Exists (Term -> Formula) - | Forall (Term -> Formula) - | TRUE - | FALSE - | Term :<: Term - | Term :>: Term - | Term :<=: Term - | Term :>=: Term - | Term :=: Term - | Term :/=: Term - | Integer :| Term - -translate :: Formula -> Form (Prop PosP) -translate = loop 0 - where loop n form = case form of - Exists f -> ex_step n (loop (n+1) (f (var n))) - Forall f -> loop n (Not (Exists (Not . f))) - Not f -> neg (loop n f) - f1 :=>: f2 -> loop n (f2 :\/: Not f1) - f1 :<=>: f2 -> loop n (f1 :/\: f2 :\/: Not f1 :/\: Not f2) - f1 :/\: f2 -> Node And (loop n f1) (loop n f2) - f1 :\/: f2 -> Node Or (loop n f1) (loop n f2) - - FALSE -> Leaf (Prop False FF) - t1 :=: t2 -> Leaf (Prop False (Bin Equal t1 t2)) - t1 :<: t2 -> Leaf (Prop False (Bin LessThan t1 t2)) - t1 :<=: t2 -> Leaf (Prop False (Bin LessThanEqual t1 t2)) - - TRUE -> Leaf (Prop True FF) - t1 :/=: t2 -> Leaf (Prop True (Bin Equal t1 t2)) - t1 :>=: t2 -> Leaf (Prop True (Bin LessThan t1 t2)) - t1 :>: t2 -> Leaf (Prop True (Bin LessThanEqual t1 t2)) - - d :| t -> Leaf (Prop False (Divides d t)) - -class Quant t where - quant :: ((Term -> Formula) -> Formula) -> t -> Formula - -instance Quant Formula where - quant _ p = p - -instance Quant t => Quant (Term -> t) where - quant q p = q (\x -> quant q (p x)) - -exists, forall :: Quant t => t -> Formula -exists p = quant Exists p -forall p = quant Forall p - --- 4: wrap term, not --- 3: wrap and --- 2: wrap or --- 1: wrap implies, quantifiers -instance PP Formula where - pp = pp1 0 -- ' 0 0 - where - pp1 :: Int -> Formula -> Doc - pp1 p form = case form of - _ :/\: _ -> hang (text "/\\") 2 (loop form) - where loop (f1 :/\: f2) = loop f1 $$ loop f2 - loop f = pp f - - _ :\/: _ -> hang (text "\\/") 2 (loop form) - where loop (f1 :\/: f2) = loop f1 $$ loop f2 - loop f = pp f - - _ -> pp' 0 p form - - - - pp' :: Int -> Name -> Formula -> Doc - pp' n p form = case form of - f1 :/\: f2 | n < 3 -> pp' 2 p f1 <+> text "/\\" <+> pp' 2 p f2 - f1 :\/: f2 | n < 2 -> pp' 1 p f1 <+> text "\\/" <+> pp' 1 p f2 - f1 :=>: f2 | n < 1 -> pp' 1 p f1 <+> text "=>" <+> pp' 0 p f2 - f1 :<=>: f2 | n < 1 -> pp' 1 p f1 <+> text "=>" <+> pp' 0 p f2 - Not f | n < 4 -> text "Not" <+> pp' 4 p f - Exists {} | n < 1 -> pp_ex (text "exists") p form - where pp_ex d q (Exists g) = pp_ex (d <+> text (var_name q)) - (q+1) (g (var q)) - pp_ex d q g = d <> text "." <+> pp' 0 q g - - Forall {} | n < 1 -> pp_ex (text "forall") p form - where pp_ex d q (Forall g) = pp_ex (d <+> text (var_name q)) - (q+1) (g (var q)) - pp_ex d q g = d <> text "." <+> pp' 0 q g - TRUE -> text "true" - FALSE -> text "false" - t1 :<: t2 | n < 4 -> pp t1 <+> text "<" <+> pp t2 - t1 :>: t2 | n < 4 -> pp t1 <+> text ">" <+> pp t2 - t1 :<=: t2 | n < 4 -> pp t1 <+> text "<=" <+> pp t2 - t1 :>=: t2 | n < 4 -> pp t1 <+> text ">=" <+> pp t2 - t1 :=: t2 | n < 4 -> pp t1 <+> text "=" <+> pp t2 - t1 :/=: t2 | n < 4 -> pp t1 <+> text "/=" <+> pp t2 - k :| t1 | n < 4 -> text (show k) <+> text "|" <+> pp t1 - _ -> parens (pp' 0 p form) -
− src/Data/Integer/Presburger/ModArith.hs
@@ -1,30 +0,0 @@-module Data.Integer.Presburger.ModArith where - -import Data.Integer.Presburger - -is_nat :: Term -> Formula -is_nat t = 0 :<=: t - -is_reminder :: Integer -> Term -> Formula -is_reminder d r = is_nat r :/\: r :<: fromIntegral d - --- | divMod t d == (q,r) -div_mod_is :: Term -> Integer -> Term -> Term -> Formula -div_mod_is t d q r = is_reminder d r :/\: d .* q + r :=: t - --- | mod t d == r -mod_is :: Term -> Integer -> Term -> Formula -mod_is t d r = is_reminder d r :/\: d :| (t - r) - -bin_op_mod :: Integer -> (Term -> Term -> Term) - -> Term -> Term -> Term -> Formula -bin_op_mod d f t1 t2 t3 = mod_is (f t1 t2) d t3 - -add_mod, sub_mod, mul_mod :: Integer -> Term -> Term -> Term -> Formula -add_mod d = bin_op_mod d (+) -sub_mod d = bin_op_mod d (-) -mul_mod d = bin_op_mod d (*) - - - -
− src/Data/Integer/Presburger/Notation.hs
@@ -1,47 +0,0 @@-module Data.Integer.Presburger.Notation- ( check- , module Data.Integer.Presburger.Notation- ) where--import Data.Integer.Presburger.Form-import Prelude hiding ((<),(<=),(==),(/=),(>),(>=), not, (&&), (||))-import qualified Prelude as P--type Formula = Form (Prop PosP)--infixr 2 ||-infixr 3 &&-infix 4 <, <=, ==, >, >=, /=----(&&), (||) :: Formula -> Formula -> Formula-f1 && f2 = Node And f1 f2-f1 || f2 = Node Or f1 f2--(<) :: Term -> Term -> Formula-t1 < t2 = Leaf $ Prop False $ Bin LessThan t1 t2--(<=) :: Term -> Term -> Formula-t1 <= t2 = Leaf $ Prop False $ Bin LessThanEqual t1 t2--(==) :: Term -> Term -> Formula-t1 == t2 = Leaf $ Prop False $ Bin Equal t1 t2--exists :: Name -> Formula -> Formula-exists x f = ex_step x f--not :: Formula -> Formula-not = neg--(>) :: Term -> Term -> Formula-t1 > t2 = not (t1 <= t2)--(>=) :: Term -> Term -> Formula-t1 >= t2 = not (t1 < t2)--(/=) :: Term -> Term -> Formula-t1 /= t2 = not (t1 == t2)--forall :: Name -> Formula -> Formula-forall x f = not (exists x (not f))
− src/Data/Integer/Presburger/Prop.hs
@@ -1,193 +0,0 @@-module Data.Integer.Presburger.Prop- ( module Data.Integer.Presburger.Prop- , module X- ) where--import Data.Integer.Presburger.Term as X---- | Possibly negated propositions.--- For example, we would express "t1 not equal to t2" like this:--- @Prop { negated = True, prop = Bin Equal t1 t2 }@-data Prop p = Prop { negated :: !Bool, prop :: !p }---- | A proposition normalized with respect to a particular variable.-data NormProp p = Ind (Prop PosP) -- ^ Independent of variable.- | Norm (Prop p) -- ^ Depends on variable.---- | Basic binary relations.-data RelOp = Equal | LessThan | LessThanEqual deriving Eq---- | Basic propositions.-data PosP = Bin !RelOp Term Term | Divides !Integer Term | FF---- | Propositions specialized to say something about a particular variable.-data VarP = NBin !RelOp Term -- ^ x `op` t- | NDivides !Integer Term -- ^ n | x + t---- | Propositions specialized for a variable with a coefficient.--- For example: 4 * x = t--- @CVarP { coeff = 4, pprop = NBin Equal t }@-data CVarP = CVarP { coeff :: !Integer, pprop :: !VarP }----- | Rewrite a propositions as a predicate about a specific variable.-norm :: Name -> Prop PosP -> NormProp CVarP-norm x p = case prop p of-- Bin op t1 t2- | k1 == k2 -> Ind p { prop = Bin op t1' t2' }- | k1 > k2 -> Norm p { prop = CVarP (k1 - k2) (NBin op (t2' - t1')) }- | otherwise -> Norm Prop { prop = CVarP (k2 - k1) (NBin op' (t1' - t2'))- , negated = neg'- }- - where (k1,t1') = split_term x t1 -- t1 = k1 * x + t1'- (k2,t2') = split_term x t2 -- t2 = k2 * x + t2'-- (neg',op') = case op of- Equal -> (negated p, Equal)- LessThan -> (not (negated p), LessThanEqual)- LessThanEqual -> (not (negated p), LessThan)- - -- a < t --> same- -- Not (a < t) --> same- -- t < a --> Not (a =< t)- -- Not (t < a) --> a =< t--- Divides n t1- | k1 == 0 -> Ind p- | k1 > 0 -> Norm p { prop = CVarP k1 (NDivides n t1') }- | otherwise -> Norm p { prop = CVarP (negate k1) (NDivides n (negate t1'))}- where(k1,t1') = split_term x t1 -- t1 = k1 * x + t1'-- FF -> Ind p----- | Eliminate variable coefficients by scaling the properties appropriately.-scale :: Integer -> Prop CVarP -> Prop VarP-scale k p =- let np = prop p- sc = k `div` coeff np- in p { prop = case pprop np of- NBin op t -> NBin op (sc .* t)- NDivides n t -> NDivides (sc * n) (sc .* t)- }----- | Evaluate a proposition for a sufficiently small value of--- the variable, if possible. Otherwise, substitute the given term.-neg_inf :: Term -> Prop VarP -> Prop PosP-neg_inf t p = case prop p of- NBin Equal _ -> Prop { negated = negated p, prop = FF }- NBin _ _ -> Prop { negated = not (negated p), prop = FF }- NDivides n t1 -> p { prop = Divides n (t + t1) }---- | Evaluate a proposition for a sufficiently large value of--- the variable, if possible. Otherwise, substitute the given term.-pos_inf :: Term -> Prop VarP -> Prop PosP-pos_inf t p = case prop p of- NDivides n t1 -> p { prop = Divides n (t + t1) }- _ -> Prop { negated = negated p, prop = FF }----- | Evaluate a proposition with a given term for the variable.-normal :: Term -> Prop VarP -> Prop PosP-normal t p = case prop p of- NBin op t1 -> p { prop = Bin op t t1 }- NDivides n t1 -> p { prop = Divides n (t + t1) }--------------------------------------------------------------------------------------- | The meanings of atomic propositions-eval_prop :: Prop PosP -> Env -> Bool-eval_prop (Prop neg p) env = if neg then not res else res- where res = case p of- FF -> False- Divides n t -> divides n (eval_term t env)- Bin op t1 t2 -> bin_op op (eval_term t1 env) (eval_term t2 env)- --bin_op :: RelOp -> Integer -> Integer -> Bool-bin_op op x y = case op of- Equal -> x == y- LessThan -> x < y- LessThanEqual -> x <= y---- | Replace a variable with a constant, in a property.-subst_prop :: Name -> Integer -> Prop PosP -> Prop PosP-subst_prop x n (Prop b p) =- case p of- Bin op t1 t2 -> Prop b (Bin op (subst_term x n t1) (subst_term x n t2))- Divides k t -> Prop b (Divides k (subst_term x n t))- FF -> Prop b FF--simplify_prop :: Prop PosP -> Prop PosP-simplify_prop it@(Prop b p) =- case p of- Divides n t -> case is_constant t of- Just v -> Prop (b /= divides n v) FF- Nothing -> it- Bin Equal t1 t2 | not b && t1 == t2 -> Prop True FF- Bin op t1 t2 -> case (is_constant t1, is_constant t2) of- (Just v1, Just v2) -> Prop (b /= bin_op op v1 v2) FF- _ -> it- FF -> it------------------------------------------------------------------------------------class SignPP t where- pp_neg :: Bool -> t -> Doc---instance SignPP RelOp where-- pp_neg False r = case r of- Equal -> text "=="- LessThan -> text "<"- LessThanEqual -> text "<="-- pp_neg True r = case r of- Equal -> text "/="- LessThan -> text ">="- LessThanEqual -> text ">"---pp_neg_div :: Bool -> Doc-pp_neg_div False = text "|"-pp_neg_div True = text "/|"---instance SignPP PosP where- pp_neg n (Bin op t1 t2) = pp t1 <+> pp_neg n op <+> pp t2- pp_neg n (Divides d t) = text (show d) <+> pp_neg_div n <+> pp t- pp_neg n FF = text (if n then "True" else "False")---instance SignPP VarP where- pp_neg n (NBin op t) = text "_" <+> pp_neg n op <+> pp t- pp_neg n (NDivides d t) = text (show d) <+> pp_neg_div n- <+> text "_ +" <+> pp t---instance SignPP CVarP where- pp_neg n p = case pprop p of- NBin op t -> it <+> pp_neg n op <+> pp t- NDivides d t -> text (show d) <+> pp_neg_div n- <+> it <+> text "+" <+> pp t- where it | c == 1 = text "_"- | c == (-1) = text "- _"- | otherwise = text (show c) <+> text "* _"-- c = coeff p - --instance SignPP p => PP (Prop p) where- pp p = pp_neg (negated p) (prop p)---instance SignPP p => PP (NormProp p) where- pp (Ind p) = pp p- pp (Norm p) = pp p-
− src/Data/Integer/Presburger/SolveDiv.hs
@@ -1,100 +0,0 @@-module Data.Integer.Presburger.SolveDiv- ( DivCtr(..), Env, elim- ) where--import Data.Integer.Presburger.Term-import Data.List(foldl')----- | A general "divisible by" constraint.-data DivCtr = Divs !Integer !Term----- | Given some variables with bounds on them, and a set of--- "divisible by" constraints, we produce all possible assignments--- to the variables that are in bounds, and satisfy the constraints.-elim :: Env -> [(Name,Integer)] -> [DivCtr] -> [ Env ]-elim env0 [] ts = if all chk ts then [ env0 ] else []- where chk (Divs x t) = divides x (eval_term t env0)-elim env0 ((x,bnd):xs) cs = do let (mb,cs1) = elim_var x cs- env <- elim env0 xs cs1- v <- case mb of- Nothing -> [ 1 .. bnd ]- Just (NDivides a b t) ->- theorem2 bnd (a,b,eval_term t env)- return (env_extend x v env)------ | "divisible by" constraint on a variable with a coefficient.-data VarDivCtr = NDivides { divisor :: !Integer- , coeff :: !Integer- , rest :: !Term- }----- | This theorem combines two "divisible by" contratints on a single--- variable, into a single constraint on the variable, and a generic--- "divisible by" constraint that does not mention the variable.-theorem1 :: VarDivCtr -> VarDivCtr -> (VarDivCtr, DivCtr)-theorem1 NDivides { divisor = m, coeff = a1, rest = b1 }- NDivides { divisor = n, coeff = a2, rest = b2 }- = (new_x, new_other)-- where (p,q,d) = extended_gcd (a1 * n) (a2 * m)-- new_x = NDivides { divisor = m * n- , coeff = d- , rest = (p * n) .* b1 + (q * m) .* b2- }-- new_other = Divs d (a2 .* b1 - a1 .* b2)----- | Repeatedly apply theorem 1 to a set of constraints,--- to split them into a single constraint on the variable,--- and additional constraints that do not mention the varibale.-elim_var :: Name -> [DivCtr] -> (Maybe VarDivCtr, [DivCtr])-elim_var x cs = case foldl' part ([],[]) cs of- ([], have_not) -> (Nothing, have_not)- (h : hs, have_not) -> let (c,hn) = step h hs have_not- in (Just c,hn)-- where part s@(have,have_not) c@(Divs m t)- | m == 1 = s -- ignore "divisible by 1" constraints.- | a == 0 = (have , c : have_not)- | otherwise = (NDivides m a b : have, have_not)- where (a,b) = split_term x t -- t = a * x + b-- step :: VarDivCtr -> [VarDivCtr] -> [DivCtr] -> (VarDivCtr,[DivCtr])- step h [] ns = (h,ns)- step h (h1:hs) ns = step h2 hs (n : ns)- where (h2,n) = theorem1 h h1------ | This theorem produces the solutions for a "divisible by" constraint--- on a variable, where the "rest" term is a constant.--- We peoduce only the solutions that are in the range [1 .. bnd]------ solutions for x in [1 .. bnd] of: m | x * a + b-theorem2 :: Integer -> (Integer,Integer,Integer) -> [Integer]-theorem2 bnd (m,a,b)- | r == 0 = [ t * k - c | t <- [ lower .. upper ] ]- | otherwise = []- where k = div m d- c = p * qu- (p,_,d) = extended_gcd a m- (qu,r) = divMod b d-- (lower',r1) = divMod (1 + c) k- lower = if r1 == 0 then lower' else lower' + 1 -- hmm- upper = div (bnd + c) k-- -- lower and upper:- -- t * k - c = 1 --> t = (1 + c) / k- -- t * k - c = bnd --> t = (bnd + c) / k----
− src/Data/Integer/Presburger/Term.hs
@@ -1,142 +0,0 @@-module Data.Integer.Presburger.Term- ( Term, Name, split_term, is_constant, (.*), var, num- , Env, env_empty, env_extend- , eval_term, subst_term- , var_name- , module U- ) where--import Data.Integer.Presburger.Utils as U--import qualified Data.IntMap as Map-import Data.Maybe(fromMaybe)-import Control.Monad(mplus,guard)----- | We represent the names of variables in terms as integers.-type Name = Int---- | Terms of Presburger arithmetic.--- Term are created by using the 'Num' class.--- WARNING: Presburger arithmetic only supports multiplication--- by a constant, trying to create invalid terms will result--- in a run-time error. A more type-safe alternative is to--- use the '(.*)' operator.-data Term = Term (Map.IntMap Integer) Integer----- | @split_term x (n * x + t1) = (n,t1)@--- @x@ does not occur in @t1@-split_term :: Name -> Term -> (Integer,Term)-split_term x (Term m n) = (fromMaybe 0 c, Term m1 n)- where (c,m1) = Map.updateLookupWithKey (\_ _ -> Nothing) x m--var :: Name -> Term-var x = Term (Map.singleton x 1) 0--num :: Integer -> Term-num n = Term Map.empty n----- Evaluation -------------------------------------------------------------------newtype Env = Env (Map.IntMap Integer)--env_empty :: Env-env_empty = Env (Map.empty)--env_extend :: Name -> Integer -> Env -> Env-env_extend x v (Env m) = Env (Map.insert x v m)---- The meaning of a term with free variables--- If the term contains free variables that are not defined, then--- we assume that these variables are 0.-eval_term :: Term -> Env -> Integer-eval_term (Term m k) (Env env) = sum (k : map eval_var (Map.toList m))- where eval_var (x,c) = case Map.lookup x env of- Nothing -> 0- Just v -> c * v--subst_term :: Name -> Integer -> Term -> Term-subst_term x n t = case split_term x t of- (c, Term m k) -> Term m (k + c * n)------------------------------------------------------------------------------------instance Eq Term where- t1 == t2 = is_constant (t1 - t2) == Just 0--instance Num Term where- fromInteger n = Term Map.empty n-- Term m1 n1 + Term m2 n2 = Term (Map.unionWith (+) m1 m2) (n1 + n2)-- negate (Term m n) = Term (Map.map negate m) (negate n)-- t1 * t2 = case fmap (.* t2) (is_constant t1) `mplus`- fmap (.* t1) (is_constant t2) of- Just t -> t- Nothing -> error $ unlines [ "[(*) @ Term] Non-linear product:"- , " *** " ++ show t1- , " *** " ++ show t2- ]- signum t = case is_constant t of- Just n -> num (signum n)- Nothing -> error $ unlines [ "[signum @ Term]: Non-constant:"- , " *** " ++ show t- ]-- abs t = case is_constant t of- Just n -> num (abs n)- Nothing -> error $ unlines [ "[abs @ Term]: Non-constant:"- , " *** " ++ show t- ]----- | Check if a term is a constant (i.e., contains no variables).--- If so, then we return the constant, otherwise we return 'Nothing'.-is_constant :: Term -> Maybe Integer-is_constant (Term m n) = guard (all (0 ==) (Map.elems m)) >> return n--(.*) :: Integer -> Term -> Term-0 .* _ = 0-1 .* t = t-k .* Term m n = Term (Map.map (k *) m) (k * n)---var_name :: Name -> String-var_name x = let (a,b) = divMod x 26- rest = if a == 0 then "" else show a- in toEnum (97 + b) : rest--instance Show Term where show x = show (pp x)-instance PP Term where- pp (Term m k) | isEmpty vars = text (show k)- | k == 0 = vars- | k > 0 = vars <+> char '+' <+> text (show k)- | otherwise = vars <+> char '-' <+> text (show $ abs k)- where ppvar (x,n) = sign <+> co <+> text (var_name x)- where (sign,co)- | n == -1 = (char '-', empty)- | n < 0 = (char '-', text (show (abs n)) <+> char '*')- | n == 1 = (char '+', empty)- | otherwise = (char '+', text (show n) <+> char '*')- first_var (x,1) = text (var_name x)- first_var (x,-1) = char '-' <> text (var_name x)- first_var (x,n) = text (show n) <+> char '*' <+> text (var_name x)-- vars = case filter ((/= 0) . snd) (Map.toList m) of- [] -> empty- v : vs -> first_var v <+> hsep (map ppvar vs)---instance PP Env where- pp (Env e) = vcat (map sh (Map.toList e))- where sh (x,y) = text (var_name x) <+> text "=" <+> text (show y)--------
− src/Data/Integer/Presburger/Utils.hs
@@ -1,45 +0,0 @@-module Data.Integer.Presburger.Utils- ( module Data.Integer.Presburger.Utils- , module PP- ) where--import Text.PrettyPrint.HughesPJ as PP-----lcms :: Integral a => [a] -> a-lcms xs = foldr lcm 1 xs---groupEither :: [Either a b] -> ([a],[b])-groupEither xs = foldr cons ([],[]) xs- where cons (Left a) (as,bs) = (a:as,bs)- cons (Right b) (as,bs) = (as,b:bs)--mapEither :: (a -> Either x y) -> [a] -> ([x],[y])-mapEither f xs = groupEither (map f xs)----- | let (p,q,r) = extended_gcd x y--- in (x * p + y * q = r) && (gcd x y = r)-extended_gcd :: Integral a => a -> a -> (a,a,a)-extended_gcd arg1 arg2 = loop arg1 arg2 0 1 1 0- where loop a b x lastx y lasty- | b /= 0 = let (q,b') = divMod a b- x' = lastx - q * x- y' = lasty - q * y- in x' `seq` y' `seq` loop b b' x' x y' y- | otherwise = (lastx,lasty,a)----- We define: "d | a" as "exists y. d * y = a"-divides :: Integral a => a -> a -> Bool-0 `divides` 0 = True-0 `divides` _ = False-x `divides` y = mod y x == 0---class PP a where- pp :: a -> Doc-
+ src/Data/Integer/SAT.hs view
@@ -0,0 +1,822 @@+{-# LANGUAGE Safe, PatternGuards #-}+{-|+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+ ) where++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)++infixr 2 :||+infixr 3 :&&+infix 4 :==, :/=, :<, :<=, :>, :>=+infixl 6 :+, :-+infixl 7 :*++--------------------------------------------------------------------------------+-- Solver interface++-- | A collection of propositions.+newtype PropSet = State (Answer RW)+ deriving Show++-- | 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)++-- | 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+ enqAndGo qZeroTerms (t1 |-| t2)++prop (e1 :/= e2) = do t1 <- expr e1+ t2 <- expr e2+ let t = t1 |-| t2+ enqAndGo qNegTerms t `mplus` enqAndGo qNegTerms (tNeg t)++prop (e1 :< e2) = do t1 <- expr e1+ t2 <- expr e2+ enqAndGo qNegTerms (t1 |-| t2)++prop (e1 :<= e2) = do t1 <- expr e1+ t2 <- expr e2+ let t = t1 |-| t2+ enqAndGo qZeroTerms t `mplus` enqAndGo qNegTerms 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+ , 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 })++--------------------------------------------------------------------------------+-- 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 variabless.+ -- These form an idempotent substitution.+ } deriving Show+++-- | 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 kciked-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 a list of lower (resp. upper) bounds, compute the least (resp. largest)+-- value that satisfies them all.+iPickBounded :: BoundType -> [Bound] -> Maybe Integer+iPickBounded bt bs = go bs Nothing+ where+ go [] mb = mb+ go (Bound c t : more) mb =+ do k <- isConst t+ let t1 = maybe k (combine k) mb+ go more $ Just $ compute t1 c++ combine = case bt of+ Lower -> max+ Upper -> min++ compute v c = case bt of+ Lower -> div v c + 1+ Upper -> let (q,r) = divMod v c+ in if r == 0 then q - 1 else q+++-- | 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++-- | 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+ | 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++ -- 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+ addWork qZeroTerms (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)+++-- | 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 = addWork qNegTerms 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 ] ]+ addWork qNegTerms real+ msum (addWork qNegTerms dark : map (addWork qZeroTerms) gray)+ ) ctrs++ | otherwise = error "solveIsNeg: unreachable"+++{- 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++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 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 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 }+ )++-- | 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+ 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 = 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 }+ }++enqAndGo :: Field Term -> Term -> S ()+enqAndGo q t =+ do i <- get inerts+ addWork q $ iApSubst i t+ solveAll+++++--------------------------------------------------------------------------------+++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)+++++++