diff --git a/presburger.cabal b/presburger.cabal --- a/presburger.cabal +++ b/presburger.cabal @@ -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 diff --git a/src/Data/Integer/OldPresburger.hs b/src/Data/Integer/OldPresburger.hs new file mode 100644 --- /dev/null +++ b/src/Data/Integer/OldPresburger.hs @@ -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 + + diff --git a/src/Data/Integer/Presburger.hs b/src/Data/Integer/Presburger.hs --- a/src/Data/Integer/Presburger.hs +++ b/src/Data/Integer/Presburger.hs @@ -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 diff --git a/src/Data/Integer/Presburger/Form.hs b/src/Data/Integer/Presburger/Form.hs new file mode 100644 --- /dev/null +++ b/src/Data/Integer/Presburger/Form.hs @@ -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" + + + + diff --git a/src/Data/Integer/Presburger/HOAS.hs b/src/Data/Integer/Presburger/HOAS.hs new file mode 100644 --- /dev/null +++ b/src/Data/Integer/Presburger/HOAS.hs @@ -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) + diff --git a/src/Data/Integer/Presburger/ModArith.hs b/src/Data/Integer/Presburger/ModArith.hs new file mode 100644 --- /dev/null +++ b/src/Data/Integer/Presburger/ModArith.hs @@ -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 (*) + + + + diff --git a/src/Data/Integer/Presburger/Notation.hs b/src/Data/Integer/Presburger/Notation.hs new file mode 100644 --- /dev/null +++ b/src/Data/Integer/Presburger/Notation.hs @@ -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)) diff --git a/src/Data/Integer/Presburger/Prop.hs b/src/Data/Integer/Presburger/Prop.hs new file mode 100644 --- /dev/null +++ b/src/Data/Integer/Presburger/Prop.hs @@ -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 + diff --git a/src/Data/Integer/Presburger/SolveDiv.hs b/src/Data/Integer/Presburger/SolveDiv.hs new file mode 100644 --- /dev/null +++ b/src/Data/Integer/Presburger/SolveDiv.hs @@ -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 + + + + diff --git a/src/Data/Integer/Presburger/Term.hs b/src/Data/Integer/Presburger/Term.hs new file mode 100644 --- /dev/null +++ b/src/Data/Integer/Presburger/Term.hs @@ -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) + + + + + + + + diff --git a/src/Data/Integer/Presburger/Utils.hs b/src/Data/Integer/Presburger/Utils.hs new file mode 100644 --- /dev/null +++ b/src/Data/Integer/Presburger/Utils.hs @@ -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 +