presburger 0.3 → 0.4
raw patch · 11 files changed
+1593/−671 lines, 11 files
Files
- presburger.cabal +22/−6
- src/Data/Integer/OldPresburger.hs +673/−0
- src/Data/Integer/Presburger.hs +3/−665
- src/Data/Integer/Presburger/Form.hs +213/−0
- src/Data/Integer/Presburger/HOAS.hs +125/−0
- src/Data/Integer/Presburger/ModArith.hs +30/−0
- src/Data/Integer/Presburger/Notation.hs +47/−0
- src/Data/Integer/Presburger/Prop.hs +193/−0
- src/Data/Integer/Presburger/SolveDiv.hs +100/−0
- src/Data/Integer/Presburger/Term.hs +142/−0
- src/Data/Integer/Presburger/Utils.hs +45/−0
presburger.cabal view
@@ -1,5 +1,5 @@ Name: presburger-Version: 0.3+Version: 0.4 License: BSD3 License-file: LICENSE Author: Iavor S. Diatchki@@ -8,11 +8,27 @@ Category: Algorithms Synopsis: Cooper's decision procedure for Presburger arithmetic. Description: Cooper's decision procedure for Presburger arithmetic.-hs-source-dirs: src-Build-Depends: base, containers, pretty Build-type: Simple-Exposed-modules: Data.Integer.Presburger+Cabal-version: >= 1.6 -Extensions:-GHC-options: -O2 -Wall+library+ 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++ GHC-options: -O2 -Wall++source-repository head+ type: git+ location: git://github.com/yav/presburger.git
+ src/Data/Integer/OldPresburger.hs view
@@ -0,0 +1,673 @@+{-| 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 view
@@ -6,668 +6,6 @@ * title: "Theorem Proving in Arithmetic without Multiplication" * year: 1972 -}-module Data.Integer.Presburger- ( 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--+module Data.Integer.Presburger (module X) where+ +import Data.Integer.Presburger.HOAS as X
+ src/Data/Integer/Presburger/Form.hs view
@@ -0,0 +1,213 @@+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 view
@@ -0,0 +1,125 @@+{-# 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 view
@@ -0,0 +1,30 @@+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 view
@@ -0,0 +1,47 @@+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 view
@@ -0,0 +1,193 @@+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 view
@@ -0,0 +1,100 @@+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 view
@@ -0,0 +1,142 @@+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 view
@@ -0,0 +1,45 @@+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+