diff --git a/presburger.cabal b/presburger.cabal
--- a/presburger.cabal
+++ b/presburger.cabal
@@ -1,13 +1,14 @@
 Name:           presburger
-Version:        0.4
+Version:        1.0
 License:        BSD3
 License-file:   LICENSE
 Author:         Iavor S. Diatchki
 Homepage:       http://github.com/yav/presburger
 Maintainer:     diatchki@galois.com
 Category:       Algorithms
-Synopsis:       Cooper's decision procedure for Presburger arithmetic.
-Description:    Cooper's decision procedure for Presburger arithmetic.
+Synopsis:       A decision procedure for quantifier-free linear arithmetic.
+Description:    The decision procedure is based on the algorithm used in
+                CVC4, which is itself based on the Omega test.
 Build-type:     Simple
 Cabal-version:  >= 1.6
 
@@ -15,16 +16,7 @@
   Build-Depends:  base < 10, containers, pretty
   hs-source-dirs: src
   Exposed-modules:
-    Data.Integer.Presburger
-    Data.Integer.OldPresburger
-    Data.Integer.Presburger.Term
-    Data.Integer.Presburger.Prop
-    Data.Integer.Presburger.Form
-    Data.Integer.Presburger.SolveDiv
-    Data.Integer.Presburger.Notation
-    Data.Integer.Presburger.HOAS
-    Data.Integer.Presburger.ModArith
-    Data.Integer.Presburger.Utils
+    Data.Integer.SAT
 
   GHC-options:    -O2 -Wall
 
diff --git a/src/Data/Integer/OldPresburger.hs b/src/Data/Integer/OldPresburger.hs
deleted file mode 100644
--- a/src/Data/Integer/OldPresburger.hs
+++ /dev/null
@@ -1,673 +0,0 @@
-{-| This module implements Cooper's algorithm for deciding
-    first order formulas over integers with addition.
-
-Based on the paper:
- * author: D.C.Cooper
- * title:  "Theorem Proving in Arithmetic without Multiplication"
- * year:   1972
--}
-module Data.Integer.OldPresburger
-  ( check, simplify, Formula(..), Term, (.*), is_constant
-  , PP(..)
-  ) where
-
-
-import qualified Data.IntMap as Map
-import Data.Maybe(fromMaybe)
-import Data.List(nub,foldl')
-import Control.Monad(mplus,guard)
-import Prelude hiding (LT,EQ)
-
-import Text.PrettyPrint.HughesPJ
-
-
--- | Check if a formula is true.
-check :: Formula -> Bool
-check f = eval_form (pre (True,0) f)
-
-simplify :: Formula -> Formula
-simplify f = invert (pre (True,0) f)
-
--- Sugar -----------------------------------------------------------------------
-
-
-infixl 3 :/\:
-infixl 2 :\/:
-infixr 1 :=>:
-
-infix 4 :<:, :<=:, :>:, :>=:, :=:, :/=:, :|
-
-
--- Forst-oreder formulas for Presburger arithmetic.
-data Formula  = Formula :/\: Formula
-              | Formula :\/: Formula
-              | Formula :=>: Formula
-              | Not Formula
-              | Exists (Term -> Formula)
-              | Forall (Term -> Formula)
-              | TRUE
-              | FALSE
-              | Term :<:   Term
-              | Term :>:   Term
-              | Term :<=:  Term
-              | Term :>=:  Term
-              | Term :=:   Term
-              | Term :/=:  Term
-              | Integer :| Term
-
-pre :: (Bool,Int) -> Formula -> Form
-pre n form = case form of
-  f1 :/\: f2        -> and' (pre n f1) (pre n f2)
-  f1 :\/: f2        -> or'  (pre n f1) (pre n f2)
-  f1 :=>: f2        -> pre n (Not f1 :\/: f2)
-  Exists f          -> pre_ex (top,x + 1) [x] (f (var x))
-    where (top,x) = n
-  Forall f          -> pre n (Not (Exists (Not . f)))
-  TRUE              -> tt'
-  FALSE             -> ff'
-  t1 :<: t2         -> lt' t1 t2
-  t1 :>: t2         -> lt' t2 t1
-  t1 :<=: t2        -> leq' t1 t2
-  t1 :>=: t2        -> leq' t2 t1
-  t1 :=: t2         -> eq' t1 t2
-  t1 :/=: t2        -> neq' t1 t2
-  k :| t            -> divs' k t
-  Not form1 -> case form1 of
-    Not f           -> pre n f
-    Forall f        -> pre n (Exists (Not . f))
-    _               -> not' (pre n form1)
-
-pre_ex :: (Bool,Int) -> [Name] -> Formula -> Form
-pre_ex (top,n) xs form = case form of
-  Exists f          -> pre_ex (top,n+1) (n:xs) (f (var n))
-  f1 :\/: f2        -> or' (pre_ex (top,n) xs f1) (pre_ex (top,n) xs f2)
-  Not form1 ->
-    case form1 of
-      Not form2     -> pre_ex (top,n) xs form2
-      Forall f      -> pre_ex (top,n) xs (Exists (Not . f))
-      p :/\: q      -> pre_ex (top,n) xs (Not p :\/: Not q)
-      _             -> exists_many top xs (pre (False,n) form)
-  _                 -> exists_many top xs (pre (False,n) form)
-
-invert :: Form -> Formula
-invert form = case form of
-  Conn And f1 f2 -> invert f1 :/\: invert f2
-  Conn Or  f1 f2 -> invert f1 :\/: invert f2
-  Prop prop -> case prop of
-    Pred FF   True      :> []      -> FALSE
-    Pred FF   False     :> []      -> TRUE
-    Pred LT   True      :> [t1,t2] -> t1 :<: t2
-    Pred LT   False     :> [t1,t2] -> t1 :>=: t2
-    Pred LEQ  True      :> [t1,t2] -> t1 :<=: t2
-    Pred LEQ  False     :> [t1,t2] -> t1 :>: t2
-    Pred EQ   True      :> [t1,t2] -> t1 :=: t2
-    Pred EQ   False     :> [t1,t2] -> t1 :/=: t2
-    Pred (Divs n) True  :> [t]     -> n :| t
-    Pred (Divs n) False :> [t]     -> Not (n :| t)
-    _ -> error "(bug) Type error in 'invert'"
-
-
--- Terms ----------------------------------------------------------------------
-
--- | Terms of Presburger arithmetic.
--- Term are created by using the 'Num' class.
--- WARNING: Presburger arithmetic only supports multiplication
--- by a constant, trying to create invalid terms will result
--- in a run-time error.  A more type-safe alternative is to
--- use the '(.*)' operator.
-data Term           = Term (Map.IntMap Integer) Integer
-
-
-type Name           = Int
-
--- | @split_term x (n * x + t1) = (n,t1)@
--- @x@ does not occur in @t1@
-split_term         :: Name -> Term -> (Integer,Term)
-split_term x (Term m n) = (fromMaybe 0 c, Term m1 n)
-  where (c,m1) = Map.updateLookupWithKey (\_ _ -> Nothing) x m
-
-var                :: Name -> Term
-var x               = Term (Map.singleton x 1) 0
-
-num                :: Integer -> Term
-num n               = Term Map.empty n
-
-
---------------------------------------------------------------------------------
-
-instance Eq Term where
-  t1 == t2  = is_constant (t1 - t2) == Just 0
-
-instance Num Term where
-  fromInteger n             = Term Map.empty n
-
-  Term m1 n1 + Term m2 n2   = Term (Map.unionWith (+) m1 m2) (n1 + n2)
-
-  negate (Term m n)         = Term (Map.map negate m) (negate n)
-
-  t1 * t2  = case fmap (.* t2) (is_constant t1) `mplus`
-                  fmap (.* t1) (is_constant t2) of
-               Just t  -> t
-               Nothing -> error $ unlines [ "[(*) @ Term] Non-linear product:"
-                                          , "  *** " ++ show t1
-                                          , "  *** " ++ show t2
-                                          ]
-  signum t  = case is_constant t of
-                Just n  -> num (signum n)
-                Nothing -> error $ unlines [ "[signum @ Term]: Non-constant:"
-                                           , " *** " ++ show t
-                                           ]
-
-  abs t     = case is_constant t of
-                Just n  -> num (abs n)
-                Nothing -> error $ unlines [ "[abs @ Term]: Non-constant:"
-                                           , " *** " ++ show t
-                                           ]
-
-
--- | Check if a term is a constant (i.e., contains no variables).
--- If so, then we return the constant, otherwise we return 'Nothing'.
-is_constant :: Term -> Maybe Integer
-is_constant (Term m n) = guard (all (0 ==) (Map.elems m)) >> return n
-
-(.*) :: Integer -> Term -> Term
-0 .* _        = 0
-1 .* t        = t
-k .* Term m n = Term (Map.map (k *) m) (k * n)
-
-
--- Formulas --------------------------------------------------------------------
-
-data PredSym    = FF | LT | LEQ | EQ | Divs Integer {- +ve -}
-data Pred       = Pred PredSym Bool -- Bool: positive (i.e. non-negated)?
-data Prop       = Pred :> [Term]
-data Conn       = And | Or deriving Eq
-data Form       = Conn Conn Form Form | Prop Prop
-
-abs_form       :: Form -> ([Prop],[Prop] -> Form)
-abs_form fo     = let (ps,skel) = loop [] fo
-                  in (reverse ps, fst . skel)
-  where loop ps (Conn c p q) =
-          let (ps1,f1) = loop ps p
-              (ps2,f2) = loop ps1 q
-          in (ps2, \fs -> let (p1,fs1) = f1 fs
-                              (p2,fs2) = f2 fs1
-                          in (Conn c p1 p2, fs2))
-        loop ps (Prop p) = (p:ps, \(f:fs) -> (Prop f,fs))
-
-
-not' :: Form -> Form
-not' (Conn c t1 t2) = Conn (not_conn c) (not' t1) (not' t2)
-not' (Prop p)       = Prop (not_prop p)
-
-ff' :: Form
-ff' = Prop $ Pred FF True :>[]
-
-tt' :: Form
-tt' = Prop $ Pred FF False :>[]
-
-lt' :: Term -> Term -> Form
-lt' t1 t2 = Prop $ Pred LT True :> [t1,t2]
-
-leq' :: Term -> Term -> Form
-leq' t1 t2 = Prop $ Pred LEQ True :> [t1,t2]
-
-eq' :: Term -> Term -> Form
-eq' t1 t2 = Prop $ Pred EQ True :> [t1,t2]
-
-neq' :: Term -> Term -> Form
-neq' t1 t2 = Prop $ Pred EQ False :> [t1,t2]
-
-and' :: Form -> Form -> Form
-and' p q = Conn And p q
-
-or' :: Form -> Form -> Form
-or' p q = Conn Or p q
-
-divs' :: Integer -> Term -> Form
-divs' n t = Prop $ Pred (Divs n) True :> [t]
-
-ors' :: [Form] -> Form
-ors' [] = ff'
-ors' xs = foldr1 or' xs
-
-not_conn :: Conn -> Conn
-not_conn And = Or
-not_conn Or  = And
-
-not_prop :: Prop -> Prop
-not_prop (f :> ts) = not_pred f :> ts
-
-not_pred :: Pred -> Pred
-not_pred (Pred p pos) = Pred p (not pos)
-
-
-
--- Eliminating existential quantifiers -----------------------------------------
-
-data NormProp = Ind Prop
-              | L Pred Term
-
-norm2 :: Name -> Integer -> Pred -> Term -> Term -> (Integer,NormProp)
-norm2 x final_k p t1 t2
-  | k1 == k2   = (1, Ind (p :> [t1',t2']))
-  | k1 > k2    = (abs k, L p t)
-  | otherwise  = (abs k, L p' t)
-
-  where (k1,t1') = split_term x t1
-        (k2,t2') = split_term x t2
-
-        k   = k1 - k2
-        t   = (final_k `div` k) .* (t2' - t1')   -- only used when k /= 0
-
-        p'  = case p of
-                Pred LT b  -> Pred LEQ (not b)
-                Pred LEQ b -> Pred LT (not b)
-                _          -> p
-
-norm1 :: Name -> Integer -> Pred -> Term -> (Integer,NormProp)
-norm1 x final_k p@(Pred (Divs d) b) t
-  | k == 0    = (1, Ind (p :> [t]))
-  | otherwise = (abs k, L ps (l .* t'))
-
-  where (k,t')  = split_term x t
-        l       = final_k `div` k
-        ps      = Pred (Divs (d * abs l)) b
-
-norm1 _ _ _ _ = error "(bug) norm1 applied to a non-unary operator"
-
-
-norm_prop :: Name -> Integer -> Prop -> (Integer,NormProp)
-norm_prop _ _ p@(_ :> [])           = (1,Ind p)
-norm_prop x final_k (p :> [t])      = norm1 x final_k p t
-norm_prop x final_k (p :> [t1,t2])  = norm2 x final_k p t1 t2
-norm_prop _ _ _                     = error "(bug) norm_prop on arity > 2"
-
--- The integer is "length as - length bs"
-a_b_sets :: (Integer,[Term],[Term]) -> NormProp -> (Integer,[Term],[Term])
-a_b_sets (o,as,bs) p = case p of
-  Ind _ -> (o,as,bs)
-
-  L (Pred op True) t ->
-    case op of
-      LT  -> (1 + o , t     : as,         bs)
-      LEQ -> (1 + o , (t+1) : as,         bs)
-      EQ  -> (o     , (t+1) : as, (t-1) : bs)
-      _   -> (o     ,         as,         bs)
-
-  L (Pred op False) t ->
-    case op of
-      LT  -> (o - 1 ,         as, (t-1) : bs)
-      LEQ -> (o - 1 ,         as, t     : bs)
-      EQ  -> (o     , t     : as, t     : bs)
-      _   -> (o     ,         as,         bs)
-
-
-analyze_props :: Name -> [Prop] -> ( [NormProp]
-                                   , Integer    -- scale
-                                   , Integer    -- bound
-                                   , Either [Term] [Term]  -- A set or B set
-                                   )
-analyze_props x ps = (ps1, final_k, bnd, if o < 0 then Left as else Right bs)
-  where (ks,ps1)  = unzip $ map (norm_prop x final_k) ps
-        final_k   = lcms ks
-        (o,as,bs) = foldl' a_b_sets (0,[],[]) ps1
-        bnd       = lcms (final_k : [ d | L (Pred (Divs d) _) _ <- ps1 ])
-
-from_bool :: Bool -> Prop
-from_bool True  = Pred FF False :> []
-from_bool False = Pred FF True :> []
-
-neg_inf :: NormProp -> Term -> Prop
-neg_inf prop t = case prop of
-  Ind p -> p
-  L ps@(Pred op pos) t1 -> case op of
-    LT      -> from_bool pos
-    LEQ     -> from_bool pos
-    EQ      -> from_bool (not pos)
-    Divs {} -> ps :> [t + t1]
-    FF      -> error "(bug) FF in NormPred"
-
-pos_inf :: NormProp -> Term -> Prop
-pos_inf prop t = case prop of
-  Ind p -> p
-  L ps@(Pred op pos) t1 -> case op of
-    LT      -> from_bool (not pos)
-    LEQ     -> from_bool (not pos)
-    EQ      -> from_bool (not pos)
-    Divs {} -> ps :> [t + t1]
-    FF      -> error "(bug) FF in NormPred"
-
-normal :: NormProp -> Term -> Prop
-normal prop t = case prop of
-  Ind p -> p
-  L ps@(Pred (Divs {}) _) t1  -> ps :> [t + t1]
-  L ps t1                     -> ps :> [t,t1]
-
-
-data Ex = Ex [(Name,Integer)]
-             [Constraint]
-             [Prop]
-
-exists_many :: Bool -> [Name] -> Form -> Form
-exists_many top xs f  = ors'
-                  $ map exp_f
-                  $ foldr (concatMap . ex_step) [Ex [] [] ps] (nub xs)
-  where (ps,skel) = abs_form f
-        exp_f = if top then expand_top skel else expand skel
-
-
-ex_step :: Name -> Ex -> [Ex]
-ex_step x (Ex xs ds ps) = case as_or_bs of
-  Left as ->
-    ( let arg = negate (var x)
-      in Ex ((x,d) : xs) (constr arg) (map (`pos_inf` arg) ps1)
-    ) : [ let arg = a - var x
-          in Ex ((x,d) : xs) (constr arg) (map (`normal` arg) ps1) | a <- as ]
-
-  Right bs ->
-    ( let arg = var x
-      in Ex ((x,d) : xs) (constr arg) (map (`neg_inf` arg) ps1)
-    ) : [ let arg = b + var x
-          in Ex ((x,d) : xs) (constr arg) (map (`normal` arg) ps1) | b <- bs ]
-
-  where (ps1,k,d',as_or_bs) = analyze_props x ps
-        d = lcms (d' : map fst ds)
-        constr t = if k == 1 then ds else (k,t) : ds
-
-
-expand_top :: ([Prop] -> Form) -> Ex -> Form
-expand_top skel (Ex xs ds ps) =
-  ors' [ skel (map (subst_prop env) ps) | env <- elim xs ds ]
-
-expand :: ([Prop] -> Form) -> Ex -> Form
-expand skel (Ex xs ds ps) =
-  ors' [ foldr and' (skel (map (subst_prop env) ps)) (map (`ctr` env) ds)
-            | env <- envs xs ]
-
-  where envs []         = [ Map.empty ]
-        envs ((x,bnd):qs) = [ Map.insert x v env
-                                      | env <- envs qs, v <- [ 1 .. bnd ] ]
-
-        ctr (k,t) env = Prop (Pred (Divs k) True :> [ subst_term env t ])
-
-
-
-type Env = Map.IntMap Integer
-
-subst_prop :: Env -> Prop -> Prop
-subst_prop env (p :> ts) = p :> map (subst_term env) ts
-
-subst_term :: Env -> Term -> Term
-subst_term env (Term m n) =
-  let (xs,vs) = unzip $ Map.toList $ Map.intersectionWith (*) env m
-  in Term (foldl' (flip Map.delete) m xs) (foldl' (+) n vs)
-
-
-
-
--- Evaluation ------------------------------------------------------------------
-
--- The meanings of formulas.
-eval_form :: Form -> Bool
-eval_form (Conn c p q) = eval_conn c (eval_form p) (eval_form q)
-eval_form (Prop p)     = eval_prop p
-
--- The meanings of connectives.
-eval_conn :: Conn -> Bool -> Bool -> Bool
-eval_conn And = (&&)
-eval_conn Or  = (||)
-
--- The meanings of atomic propositions.
-eval_prop :: Prop -> Bool
-eval_prop (Pred p pos :> ts) = if pos then res else not res
-  where res = eval_pred p (map eval_term ts)
-
--- The meanings of predicate symbols.
-eval_pred :: PredSym -> [Integer] -> Bool
-eval_pred p ts = case (p,ts) of
-  (FF,     [])    -> False
-  (Divs d, [k])   -> divides d k
-  (LT,     [x,y]) -> x < y
-  (LEQ,    [x,y]) -> x <= y
-  (EQ,     [x,y]) -> x == y
-  _               -> error "Type error"
-
--- We define: "d | a" as "exists y. d * y = a"
-divides :: Integral a => a -> a -> Bool
-0 `divides` 0 = True
-0 `divides` _ = False
-x `divides` y = mod y x == 0
-
--- The meaning of a term with no free variables.
--- NOTE: We do not check that there are no free variables.
-eval_term :: Term -> Integer
-eval_term (Term _ k) = k
-
--- The meaning of a term with free variables
-eval_term_env :: Term -> Env -> Integer
-eval_term_env (Term m k) env = sum (k : map eval_var (Map.toList m))
-  where eval_var (x,c) = case Map.lookup x env of
-                           Nothing -> error "free var"
-                           Just v  -> c * v
---------------------------------------------------------------------------------
-
-
--- Solving divides constraints -------------------------------------------------
--- See the paper's appendix.
-
-
--- | let (p,q,r) = extended_gcd x y
---   in (x * p + y * q = r)  &&  (gcd x y = r)
-extended_gcd :: Integral a => a -> a -> (a,a,a)
-extended_gcd arg1 arg2 = loop arg1 arg2 0 1 1 0
-  where loop a b x lastx y lasty
-          | b /= 0    = let (q,b') = divMod a b
-                            x'     = lastx - q * x
-                            y'     = lasty - q * y
-                        in x' `seq` y' `seq` loop b b' x' x y' y
-          | otherwise = (lastx,lasty,a)
-
-
-type Constraint     = (Integer,Term)
-type VarConstraint  = (Integer,Integer,Term)
-
--- m | (x * a1 + b1) /\ (n | x * a2 + b2)
-theorem1 :: VarConstraint -> VarConstraint -> (VarConstraint, Constraint)
-theorem1 (m,a1,b1) (n,a2,b2) = (new_x, new_other)
-  where new_x     = (m * n, d, (p*n) .* b1 + (q * m) .* b2)
-        new_other = (d, a2 .* b1 - a1 .* b2)
-
-        (p,q,d)   = extended_gcd (a1 * n) (a2 * m)
-
--- solutions for x in [1 .. bnd] of: m | x * a + b
-theorem2 :: Integer -> (Integer,Integer,Integer) -> [Integer]
-theorem2 bnd (m,a,b)
-  | r == 0      = [ t * k - c | t <- [ lower .. upper ] ]
-  | otherwise   = []
-  where k           = div m d
-        c           = p * qu
-        (p,_,d)     = extended_gcd a m
-        (qu,r)      = divMod b d
-
-        (lower',r1) = divMod (1 + c) k
-        lower       = if r1 == 0 then lower' else lower' + 1  -- hmm
-        upper       = div (bnd + c) k
-
-  -- lower and upper:
-  -- t * k - c = 1   --> t = (1 + c) / k
-  -- t * k - c = bnd --> t = (bnd + c) / k
-
-
-
-
-elim :: [(Name,Integer)] -> [Constraint] -> [ Env ]
-elim [] ts = if all chk ts then [ Map.empty ] else []
-  where chk (x,t) = divides x (eval_term t)
-elim ((x,bnd):xs) cs = do env <- elim xs cs1
-                          v <- case mb of
-                                 Nothing      -> [ 1 .. bnd ]
-                                 Just (a,b,t) ->
-                                   theorem2 bnd (a,b,eval_term_env t env)
-                          return (Map.insert x v env)
-
-  where (mb,cs1) = elim_var x cs
-
-
-
-
-elim_var :: Name -> [Constraint] -> (Maybe VarConstraint, [Constraint])
-elim_var x cs = case foldl' part ([],[]) cs of
-                  ([], have_not)     -> (Nothing, have_not)
-                  (h : hs, have_not) -> let (c,hn) = step h hs have_not
-                                        in (Just c,hn)
-  where part s@(have,have_not) c@(m,t)
-          | m == 1      = s
-          | a == 0      = (have        , c:have_not)
-          | otherwise   = ((m,a,b):have,   have_not)
-            where (a,b) = split_term x t
-
-        step :: VarConstraint -> [VarConstraint] -> [Constraint]
-             -> (VarConstraint,[Constraint])
-        step h [] ns      = (h,ns)
-        step h (h1:hs) ns = step h2 hs (n : ns)
-          where (h2,n) = theorem1 h h1
-
--- Misc -----------------------------------------------------------------------
-
-lcms :: Integral a => [a] -> a
-lcms xs = foldr lcm 1 xs
-
-
--- Pretty Printing -------------------------------------------------------------
-
-class PP a where
-  pp :: a -> Doc
-
-
-var_name           :: Name -> String
-var_name x          = let (a,b) = divMod x 26
-                          rest = if a == 0 then "" else show a
-                      in toEnum (97 + b) : rest
-
-instance Show Term where show x = show (pp x)
-instance PP Term where
-  pp (Term m k) | isEmpty vars  = text (show k)
-                | k == 0        = vars
-                | k > 0         = vars <+> char '+' <+> text (show k)
-                | otherwise     = vars <+> char '-' <+> text (show $ abs k)
-    where ppvar (x,n) = sign <+> co <+> text (var_name x)
-            where (sign,co)
-                     | n == -1    = (char '-', empty)
-                     | n < 0      = (char '-', text (show (abs n)) <+> char '*')
-                     | n == 1     = (char '+', empty)
-                     | otherwise  = (char '+', text (show n) <+> char '*')
-          first_var (x,1)  = text (var_name x)
-          first_var (x,-1) = char '-' <> text (var_name x)
-          first_var (x,n)  = text (show n) <+> char '*' <+> text (var_name x)
-
-          vars = case filter ((/= 0) . snd) (Map.toList m) of
-                   []     -> empty
-                   v : vs -> first_var v <+> hsep (map ppvar vs)
-
-
--- 4: wrap term, not
--- 3: wrap and
--- 2: wrap or
--- 1: wrap implies, quantifiers
-instance PP Formula where
-  pp = pp1 0 -- ' 0 0
-    where
-    pp1 :: Int -> Formula -> Doc
-    pp1 p form = case form of
-      _ :/\: _ -> hang (text "/\\") 2 (loop form)
-        where loop (f1 :/\: f2) = loop f1 $$ loop f2
-              loop f            = pp f
-
-      _ :\/: _ -> hang (text "\\/") 2 (loop form)
-        where loop (f1 :\/: f2) = loop f1 $$ loop f2
-              loop f            = pp f
-
-      _ -> pp' 0 p form
-
-
-
-    pp' :: Int -> Name -> Formula -> Doc
-    pp' n p form = case form of
-      f1 :/\: f2 | n < 3  -> pp' 2 p f1 <+> text "/\\" <+> pp' 2 p f2
-      f1 :\/: f2 | n < 2  -> pp' 1 p f1 <+> text "\\/" <+> pp' 1 p f2
-      f1 :=>: f2 | n < 1  -> pp' 1 p f1 <+> text "=>" <+> pp' 0 p f2
-      Not f      | n < 4  -> text "Not" <+> pp' 4 p f
-      Exists {}  | n < 1  -> pp_ex (text "exists") p form
-        where pp_ex d q (Exists g) = pp_ex (d <+> text (var_name q))
-                                                          (q+1) (g (var q))
-              pp_ex d q g          = d <> text "." <+> pp' 0 q g
-
-      Forall {} | n < 1 -> pp_ex (text "forall") p form
-        where pp_ex d q (Forall g) = pp_ex (d <+> text (var_name q))
-                                                          (q+1) (g (var q))
-              pp_ex d q g          = d <> text "." <+> pp' 0 q g
-      TRUE        -> text "true"
-      FALSE       -> text "false"
-      t1 :<:  t2 | n < 4  -> pp t1 <+> text "<"  <+> pp t2
-      t1 :>:  t2 | n < 4  -> pp t1 <+> text ">"  <+> pp t2
-      t1 :<=: t2 | n < 4  -> pp t1 <+> text "<=" <+> pp t2
-      t1 :>=: t2 | n < 4  -> pp t1 <+> text ">=" <+> pp t2
-      t1 :=:  t2 | n < 4  -> pp t1 <+> text "="  <+> pp t2
-      t1 :/=: t2 | n < 4  -> pp t1 <+> text "/=" <+> pp t2
-      k :| t1    | n < 4  -> text (show k) <+> text "|" <+> pp t1
-      _ -> parens (pp' 0 p form)
-
-instance Show Formula where show = show . pp
-
-
-
-instance PP PredSym where
-  pp p = case p of
-    FF      -> text "false"
-    LT      -> text "<"
-    LEQ     -> text "<="
-    EQ      -> text "==="
-    Divs n  -> text (show n) <+> text "|"
-
-instance PP Pred where
-  pp (Pred p True) = pp p
-  pp (Pred p False) = case p of
-    FF      -> text "true"
-    LT      -> text ">="
-    LEQ     -> text ">"
-    EQ      -> text "=/="
-    Divs n  -> text (show n) <+> text "/|"
-
-instance Show Prop where show = show . pp
-instance PP Prop where
-  pp (p :> [t1,t2]) = pp t1 <+> pp p <+> pp t2
-  pp (p :> ts)      = pp p <+> hsep (map pp ts)
-
-
-instance PP Conn where
-  pp And  = text "/\\"
-  pp Or   = text "\\/"
-
-instance PP Form where
-  pp me@(Conn c _ _) = hang (pp c) 2 (vcat $ map pp $ jn me [])
-    where jn (Conn c1 p1 q1) fs | c == c1 = jn p1 (jn q1 fs)
-          jn f fs = f : fs
-  pp (Prop p)     = pp p
-
-instance PP NormProp where
-  pp (Ind p)  = pp p
-  pp (L p@(Pred (Divs {}) _) t) = pp p <+> text "_ +" <+> pp t
-  pp (L p t)                    = text "_" <+> pp p <+> pp t
-
-instance Show NormProp where show = show . pp
-
-instance PP Ex where
-  pp (Ex xs ps ss) = hang (text "OR" <+> hsep (map quant xs)) 2
-             ( text "!" <+> hsep (map (parens . divs) ps)
-            $$ vcat (map pp ss)
-             )
-    where quant (x,n) = parens $ text (var_name x) <> colon <> text (show n)
-          divs (x,t)  = text (show x) <+> text "|" <+> pp t
-
-
diff --git a/src/Data/Integer/Presburger.hs b/src/Data/Integer/Presburger.hs
deleted file mode 100644
--- a/src/Data/Integer/Presburger.hs
+++ /dev/null
@@ -1,11 +0,0 @@
-{-| This module implements Cooper's algorithm for deciding
-    first order formulas over integers with addition.
-
-Based on the paper:
- * author: D.C.Cooper
- * title:  "Theorem Proving in Arithmetic without Multiplication"
- * year:   1972
--}
-module Data.Integer.Presburger (module X) where
-  
-import Data.Integer.Presburger.HOAS as X
diff --git a/src/Data/Integer/Presburger/Form.hs b/src/Data/Integer/Presburger/Form.hs
deleted file mode 100644
--- a/src/Data/Integer/Presburger/Form.hs
+++ /dev/null
@@ -1,213 +0,0 @@
-module Data.Integer.Presburger.Form
-  ( module Data.Integer.Presburger.Form
-  , module Data.Integer.Presburger.Prop
-  ) where
-
-import Data.Integer.Presburger.Prop
-import Data.Integer.Presburger.SolveDiv
-
-check :: Form (Prop PosP) -> Bool
-check f = eval_form f env_empty
-
-
-data Conn       = And | Or deriving Eq
-data Form p     = Node !Conn (Form p) (Form p)
-                | Leaf !p
-
-                -- A special form of disjunction. Bool = negated?
-                | Ex Bool (Name,Integer) (Form p)
-
-instance Functor Form where
-  fmap f (Node c f1 f2)    = Node c (fmap f f1) (fmap f f2)
-  fmap f (Ex b xs g)       = Ex b xs (fmap f g)
-  fmap f (Leaf p)          = Leaf (f p)
-
-form_lcm                  :: Form (NormProp CVarP) -> Integer
-form_lcm (Node _ f1 f2)    = lcm (form_lcm f1) (form_lcm f2)
-form_lcm (Leaf p)          = case p of
-                               Ind {}  -> 1
-                               Norm p1 -> coeff (prop p1)
-form_lcm (Ex _ _ f)        = form_lcm f
-
-
-
-form_scale  :: Name -> Form (Prop PosP) -> Form (NormProp VarP)
-form_scale x form
-  | k /= 1    = Node And (Leaf $ Norm $ Prop False $ NDivides k 0) sf
-  | otherwise = sf
-  where
-  nf  = fmap (norm x) form
-  k   = form_lcm nf
-  sf  = fmap leaf nf
-
-  leaf p = case p of
-             Ind p1  -> Ind p1
-             Norm p1 -> Norm (scale k p1)
-
-
--- The integer is "length as - length bs"
-a_b_sets :: (Integer,[Term],[Term]) -> NormProp VarP -> (Integer,[Term],[Term])
-a_b_sets (o,as,bs) p = case p of
-  Ind _                       -> (o,as,bs)
-  Norm (Prop _ (NDivides {})) -> (o,as,bs)
-
-  -- positive
-  Norm (Prop False (NBin op t)) ->
-    case op of
-      LessThan      -> (1 + o , t     : as,         bs)
-      LessThanEqual -> (1 + o , (t+1) : as,         bs)
-      Equal         -> (o     , (t+1) : as, (t-1) : bs)
-
-  -- negative
-  Norm (Prop True (NBin op t)) ->
-    case op of
-      LessThan      -> (o - 1 ,         as, (t-1) : bs)
-      LessThanEqual -> (o - 1 ,         as, t     : bs)
-      Equal         -> (o     , t     : as, t     : bs)
-
-
-form_pos_inf :: Term -> Form (NormProp VarP) -> Form (Prop PosP)
-form_pos_inf t form = fmap leaf form
-  where leaf p = case p of
-                   Ind p1  -> p1
-                   Norm p1 -> pos_inf t p1
-
-form_neg_inf :: Term -> Form (NormProp VarP) -> Form (Prop PosP)
-form_neg_inf t form = fmap leaf form
-  where leaf p  = case p of
-                    Ind p1  -> p1
-                    Norm p1 -> neg_inf t p1
-
-form_no_inf :: Term -> Form (NormProp VarP) -> Form (Prop PosP)
-form_no_inf t form  = fmap leaf form
-  where leaf p  = case p of
-                    Ind p1  -> p1
-                    Norm p1 -> normal t p1
-
-
-neg :: Form (Prop PosP) -> Form (Prop PosP)
-neg (Node And f1 f2)  = Node Or (neg f1) (neg f2)
-neg (Node Or f1 f2)   = Node And (neg f1) (neg f2)
-neg (Ex b x f)        = Ex (not b) x f
-neg (Leaf (Prop b p)) = Leaf (Prop (not b) p)
-
-
-simplify :: Form (Prop PosP) -> Form (Prop PosP)
-simplify (Node c f1 f2) =
-  case simplify f1 of
-    r@(Leaf (Prop n FF)) | n && c == Or
-                        || not n && c == And -> r
-                         | otherwise -> simplify f2
-    r1 -> case simplify f2 of
-            r@(Leaf (Prop n FF)) | n && c == Or
-                                || not n && c == And -> r
-                                 | otherwise -> r1
-            r2 -> Node c r1 r2
-
-
-
-simplify (Ex False (x,1) f) = simplify (subst_form x 1 f)
-simplify (Ex True (x,1) f)  = simplify (neg (subst_form x 1 f))
-
-simplify (Ex b x f) = case simplify f of
-                        Leaf (Prop n FF) -> Leaf (Prop (not (b == n)) FF)
-                        f1               -> Ex b x f1
-                              
-simplify (Leaf l) = Leaf (simplify_prop l)
-
-
-
-ex_step :: Name -> Form (Prop PosP) -> Form (Prop PosP)
-ex_step x (Node Or f1 f2) = Node Or (ex_step x f1) (ex_step x f2)
-ex_step x f
-  | as_minus_bs >= 0    = thm_as as x norm_f
-  | otherwise           = thm_bs bs x norm_f
-  
-  where 
-  norm_f               :: Form (NormProp VarP)
-  norm_f                = form_scale x f
-
-  (as_minus_bs, as, bs) = loop (0,[],[]) norm_f
-
-  loop s (Node _ f1 f2) = loop (loop s f1) f2
-  loop s (Ex _ _ f1)    = loop s f1
-  loop s (Leaf p)       = a_b_sets s p
-
-
-
-thm_as :: [Term] -> Name -> Form (NormProp VarP) -> Form (Prop PosP)
-thm_as as x f = simplify $
-  foldr1 (Node Or) $ map (Ex False (x, form_bound f))
-                   $ form_pos_inf (negate (var x)) f
-                   : [ form_no_inf (a - var x) f | a <- as ]
-
-thm_bs :: [Term] -> Name -> Form (NormProp VarP) -> Form (Prop PosP)
-thm_bs bs x f = simplify $
-  foldr1 (Node Or) $ map (Ex False (x, form_bound f))
-                   $ form_neg_inf (var x) f
-                   : [ form_no_inf (b + var x) f | b <- bs ]
-
-
-form_bound                :: Form (NormProp VarP) -> Integer
-form_bound (Node _ f1 f2)  = lcm (form_bound f1) (form_bound f2)
-form_bound (Leaf p)        = case p of
-                               Norm (Prop _ (NDivides n _)) -> n
-                               _ -> 1
-form_bound (Ex _ _ f)      = form_bound f
-
-
--- Evaluation ------------------------------------------------------------------
-
--- The meanings of formulas.
-eval_form :: Form (Prop PosP) -> Env -> Bool
-eval_form (Node c p q) env  = eval_conn c (eval_form p env) (eval_form q env)
-eval_form (Leaf p) env      = eval_prop p env
-eval_form (Ex b x f) env0 =
-  case splt f [x] of
-    (xs,cs,f1) -> let v = any (eval_form f1) (elim env0 xs cs)
-                  in if b then not v else v
-  where splt (Ex False y f1) ys = splt f1 (y:ys)
-        splt f1 ys          = case split_divs f1 of
-                                 (ds,f2) -> (ys,ds,f2)
-        
-
-split_ands :: Form p -> [Form p]
-split_ands form = loop form []
-  where loop (Node And f1 f2) fs  = loop f1 (loop f2 fs)
-        loop f fs                 = f : fs
-
-split_divs :: Form (Prop PosP) -> ([DivCtr], Form (Prop PosP))
-split_divs form = loop (split_ands form) ([], Leaf (Prop True FF))
-  where
-  loop (Leaf (Prop False (Divides n t)) : fs) (cs, f)
-                              = loop fs (Divs n t : cs, f)
-  loop (f:fs) (cs, f1)        = loop fs (cs, Node And f f1)
-  loop [] cs                  = cs
-
-
--- The meanings of connectives.
-eval_conn :: Conn -> Bool -> Bool -> Bool
-eval_conn And = (&&)
-eval_conn Or  = (||)
-
-subst_form :: Name -> Integer -> Form (Prop PosP) -> Form (Prop PosP)
-subst_form x n f = fmap (subst_prop x n) f
---------------------------------------------------------------------------------
-
-instance PP Conn where
-  pp And  = text "/\\"
-  pp Or   = text "\\/"
-
-instance PP p => PP (Form p) where
-  pp me@(Node c _ _) = hang (pp c) 2 (vcat $ map pp $ jn me [])
-    where jn (Node c1 p1 q1) fs | c == c1 = jn p1 (jn q1 fs)
-          jn f fs = f : fs
-  pp (Leaf p)     = pp p
-
-  pp (Ex n q f) = hang (how <+> quant q <> text ".") 2 (pp f)
-    where quant (x,b) = text (var_name x) <+> text "<=" <+> text (show b)
-          how = (if n then text "Not" else empty) <+> text "Ex"
-
-
-
-
diff --git a/src/Data/Integer/Presburger/HOAS.hs b/src/Data/Integer/Presburger/HOAS.hs
deleted file mode 100644
--- a/src/Data/Integer/Presburger/HOAS.hs
+++ /dev/null
@@ -1,125 +0,0 @@
-{-# LANGUAGE FlexibleInstances #-}
-
-module Data.Integer.Presburger.HOAS
-  ( Formula(..), check, translate
-  , Quant, exists, forall
-  , Term, (.*), is_constant
-  , PP(..)
-  ) where
-
-import Data.Integer.Presburger.Form hiding (check)
-import qualified Data.Integer.Presburger.Form as F
-
-check :: Formula -> Bool
-check f = F.check (translate f)
-
-
-infixl 3 :/\:
-infixl 2 :\/:
-infixr 1 :=>:
-infix  0 :<=>:
-
-infix 4 :<:, :<=:, :>:, :>=:, :=:, :/=:, :|
-
--- Forst-oreder formulas for Presburger arithmetic.
-data Formula  = Formula :/\: Formula
-              | Formula :\/: Formula
-              | Formula :=>: Formula
-              | Formula :<=>: Formula
-              | Not Formula
-              | Exists (Term -> Formula)
-              | Forall (Term -> Formula)
-              | TRUE
-              | FALSE
-              | Term :<:   Term
-              | Term :>:   Term
-              | Term :<=:  Term
-              | Term :>=:  Term
-              | Term :=:   Term
-              | Term :/=:  Term
-              | Integer :| Term
-
-translate :: Formula -> Form (Prop PosP)
-translate = loop 0
-  where loop n form = case form of
-          Exists f    -> ex_step n (loop (n+1) (f (var n)))
-          Forall f    -> loop n (Not (Exists (Not . f)))
-          Not f       -> neg (loop n f)
-          f1 :=>: f2  -> loop n (f2 :\/: Not f1)
-          f1 :<=>: f2 -> loop n (f1 :/\: f2 :\/: Not f1 :/\: Not f2)
-          f1 :/\: f2  -> Node And (loop n f1) (loop n f2)
-          f1 :\/: f2  -> Node Or  (loop n f1) (loop n f2)
-          
-          FALSE       -> Leaf (Prop False FF)
-          t1 :=: t2   -> Leaf (Prop False (Bin Equal t1 t2))
-          t1 :<: t2   -> Leaf (Prop False (Bin LessThan t1 t2))
-          t1 :<=: t2  -> Leaf (Prop False (Bin LessThanEqual t1 t2))
-
-          TRUE        -> Leaf (Prop True FF)
-          t1 :/=: t2  -> Leaf (Prop True (Bin Equal t1 t2))
-          t1 :>=: t2  -> Leaf (Prop True (Bin LessThan t1 t2))
-          t1 :>: t2   -> Leaf (Prop True (Bin LessThanEqual t1 t2))
-            
-          d :| t      -> Leaf (Prop False (Divides d t))
-
-class Quant t where
-  quant :: ((Term -> Formula) -> Formula) -> t -> Formula
-
-instance Quant Formula where
-  quant _ p = p
-
-instance Quant t => Quant (Term -> t) where
-  quant q p = q (\x -> quant q (p x))
-
-exists, forall :: Quant t => t -> Formula
-exists p  = quant Exists p
-forall p  = quant Forall p
-
--- 4: wrap term, not
--- 3: wrap and
--- 2: wrap or
--- 1: wrap implies, quantifiers
-instance PP Formula where
-  pp = pp1 0 -- ' 0 0
-    where
-    pp1 :: Int -> Formula -> Doc
-    pp1 p form = case form of
-      _ :/\: _ -> hang (text "/\\") 2 (loop form)
-        where loop (f1 :/\: f2) = loop f1 $$ loop f2
-              loop f            = pp f
-
-      _ :\/: _ -> hang (text "\\/") 2 (loop form)
-        where loop (f1 :\/: f2) = loop f1 $$ loop f2
-              loop f            = pp f
-
-      _ -> pp' 0 p form
-
-
-
-    pp' :: Int -> Name -> Formula -> Doc
-    pp' n p form = case form of
-      f1 :/\: f2 | n < 3  -> pp' 2 p f1 <+> text "/\\" <+> pp' 2 p f2
-      f1 :\/: f2 | n < 2  -> pp' 1 p f1 <+> text "\\/" <+> pp' 1 p f2
-      f1 :=>: f2 | n < 1  -> pp' 1 p f1 <+> text "=>" <+> pp' 0 p f2
-      f1 :<=>: f2 | n < 1  -> pp' 1 p f1 <+> text "=>" <+> pp' 0 p f2
-      Not f      | n < 4  -> text "Not" <+> pp' 4 p f
-      Exists {}  | n < 1  -> pp_ex (text "exists") p form
-        where pp_ex d q (Exists g) = pp_ex (d <+> text (var_name q))
-                                                          (q+1) (g (var q))
-              pp_ex d q g          = d <> text "." <+> pp' 0 q g
-
-      Forall {} | n < 1 -> pp_ex (text "forall") p form
-        where pp_ex d q (Forall g) = pp_ex (d <+> text (var_name q))
-                                                          (q+1) (g (var q))
-              pp_ex d q g          = d <> text "." <+> pp' 0 q g
-      TRUE        -> text "true"
-      FALSE       -> text "false"
-      t1 :<:  t2 | n < 4  -> pp t1 <+> text "<"  <+> pp t2
-      t1 :>:  t2 | n < 4  -> pp t1 <+> text ">"  <+> pp t2
-      t1 :<=: t2 | n < 4  -> pp t1 <+> text "<=" <+> pp t2
-      t1 :>=: t2 | n < 4  -> pp t1 <+> text ">=" <+> pp t2
-      t1 :=:  t2 | n < 4  -> pp t1 <+> text "="  <+> pp t2
-      t1 :/=: t2 | n < 4  -> pp t1 <+> text "/=" <+> pp t2
-      k :| t1    | n < 4  -> text (show k) <+> text "|" <+> pp t1
-      _ -> parens (pp' 0 p form)
-
diff --git a/src/Data/Integer/Presburger/ModArith.hs b/src/Data/Integer/Presburger/ModArith.hs
deleted file mode 100644
--- a/src/Data/Integer/Presburger/ModArith.hs
+++ /dev/null
@@ -1,30 +0,0 @@
-module Data.Integer.Presburger.ModArith where
-
-import Data.Integer.Presburger
-
-is_nat         :: Term -> Formula
-is_nat t        = 0 :<=: t
-
-is_reminder    :: Integer -> Term -> Formula
-is_reminder d r = is_nat r :/\: r :<: fromIntegral d
-
--- | divMod t d == (q,r)
-div_mod_is     :: Term -> Integer -> Term -> Term -> Formula
-div_mod_is t d q r = is_reminder d r :/\: d .* q + r :=: t
-
--- | mod t d == r
-mod_is         :: Term -> Integer -> Term -> Formula
-mod_is t d r    = is_reminder d r :/\: d :| (t - r)
-
-bin_op_mod :: Integer -> (Term -> Term -> Term)
-           -> Term -> Term -> Term -> Formula
-bin_op_mod d f t1 t2 t3 = mod_is (f t1 t2) d t3
-
-add_mod, sub_mod, mul_mod :: Integer -> Term -> Term -> Term -> Formula
-add_mod d = bin_op_mod d (+)
-sub_mod d = bin_op_mod d (-)
-mul_mod d = bin_op_mod d (*)
-
-
-
-
diff --git a/src/Data/Integer/Presburger/Notation.hs b/src/Data/Integer/Presburger/Notation.hs
deleted file mode 100644
--- a/src/Data/Integer/Presburger/Notation.hs
+++ /dev/null
@@ -1,47 +0,0 @@
-module Data.Integer.Presburger.Notation
-  ( check
-  , module Data.Integer.Presburger.Notation
-  ) where
-
-import Data.Integer.Presburger.Form
-import Prelude hiding ((<),(<=),(==),(/=),(>),(>=), not, (&&), (||))
-import qualified Prelude as P
-
-type Formula = Form (Prop PosP)
-
-infixr 2 ||
-infixr 3 &&
-infix 4 <, <=, ==, >, >=, /=
-
-
-
-(&&), (||) :: Formula -> Formula -> Formula
-f1 && f2 = Node And f1 f2
-f1 || f2 = Node Or f1 f2
-
-(<) :: Term -> Term -> Formula
-t1 < t2 = Leaf $ Prop False $ Bin LessThan t1 t2
-
-(<=) :: Term -> Term -> Formula
-t1 <= t2 = Leaf $ Prop False $ Bin LessThanEqual t1 t2
-
-(==) :: Term -> Term -> Formula
-t1 == t2 = Leaf $ Prop False $ Bin Equal t1 t2
-
-exists :: Name -> Formula -> Formula
-exists x f = ex_step x f
-
-not :: Formula -> Formula
-not = neg
-
-(>) :: Term -> Term -> Formula
-t1 > t2 = not (t1 <= t2)
-
-(>=) :: Term -> Term -> Formula
-t1 >= t2 = not (t1 < t2)
-
-(/=) :: Term -> Term -> Formula
-t1 /= t2  = not (t1 == t2)
-
-forall :: Name -> Formula -> Formula
-forall x f = not (exists x (not f))
diff --git a/src/Data/Integer/Presburger/Prop.hs b/src/Data/Integer/Presburger/Prop.hs
deleted file mode 100644
--- a/src/Data/Integer/Presburger/Prop.hs
+++ /dev/null
@@ -1,193 +0,0 @@
-module Data.Integer.Presburger.Prop
-  ( module Data.Integer.Presburger.Prop
-  , module X
-  ) where
-
-import Data.Integer.Presburger.Term as X
-
--- | Possibly negated propositions.
--- For example, we would express "t1 not equal to t2" like this:
--- @Prop { negated = True, prop = Bin Equal t1 t2 }@
-data Prop p   = Prop { negated :: !Bool, prop :: !p }
-
--- | A proposition normalized with respect to a particular variable.
-data NormProp p = Ind (Prop PosP)   -- ^ Independent of variable.
-                | Norm (Prop p)     -- ^ Depends on variable.
-
--- | Basic binary relations.
-data RelOp    = Equal | LessThan | LessThanEqual deriving Eq
-
--- | Basic propositions.
-data PosP     = Bin !RelOp Term Term | Divides !Integer Term | FF
-
--- | Propositions specialized to say something about a particular variable.
-data VarP     = NBin !RelOp Term        -- ^ x `op` t
-              | NDivides !Integer Term  -- ^ n | x + t
-
--- | Propositions specialized for a variable with a coefficient.
--- For example: 4 * x = t
--- @CVarP { coeff = 4, pprop = NBin Equal t }@
-data CVarP    = CVarP { coeff :: !Integer, pprop :: !VarP }
-
-
--- | Rewrite a propositions as a predicate about a specific variable.
-norm :: Name -> Prop PosP -> NormProp CVarP
-norm x p = case prop p of
-
-  Bin op t1 t2
-    | k1 == k2    -> Ind  p    { prop = Bin op t1' t2' }
-    | k1 > k2     -> Norm p    { prop = CVarP (k1 - k2) (NBin op (t2' - t1')) }
-    | otherwise   -> Norm Prop { prop = CVarP (k2 - k1) (NBin op' (t1' - t2'))
-                               , negated = neg'
-                               }
-                          
-    where (k1,t1')  = split_term x t1   -- t1 = k1 * x + t1'
-          (k2,t2')  = split_term x t2   -- t2 = k2 * x + t2'
-
-          (neg',op') = case op of
-                         Equal         -> (negated p, Equal)
-                         LessThan      -> (not (negated p), LessThanEqual)
-                         LessThanEqual -> (not (negated p), LessThan)
- 
-    -- a < t        --> same
-    -- Not (a < t)  --> same
-    -- t < a        --> Not (a =< t)
-    -- Not (t < a)  --> a =< t
-
-
-  Divides n t1
-    | k1 == 0    -> Ind  p
-    | k1 > 0     -> Norm p { prop = CVarP k1 (NDivides n t1') }
-    | otherwise  -> Norm p { prop = CVarP (negate k1) (NDivides n (negate t1'))}
-    where(k1,t1') = split_term x t1     -- t1 = k1 * x + t1'
-
-  FF -> Ind p
-
-
--- | Eliminate variable coefficients by scaling the properties appropriately.
-scale :: Integer -> Prop CVarP -> Prop VarP
-scale k p =
-  let np = prop p
-      sc = k `div` coeff np
-  in p { prop = case pprop np of
-                  NBin op t    -> NBin op (sc .* t)
-                  NDivides n t -> NDivides (sc * n) (sc .* t)
-       }
-
-
--- | Evaluate a proposition for a sufficiently small value of
--- the variable, if possible.  Otherwise, substitute the given term.
-neg_inf :: Term -> Prop VarP -> Prop PosP
-neg_inf t p = case prop p of
-  NBin Equal _  -> Prop { negated = negated p, prop = FF }
-  NBin _ _      -> Prop { negated = not (negated p), prop = FF }
-  NDivides n t1 -> p    { prop = Divides n (t + t1) }
-
--- | Evaluate a proposition for a sufficiently large value of
--- the variable, if possible.  Otherwise, substitute the given term.
-pos_inf :: Term -> Prop VarP -> Prop PosP
-pos_inf t p = case prop p of
-  NDivides n t1 -> p    { prop = Divides n (t + t1) }
-  _             -> Prop { negated = negated p, prop = FF }
-
-
--- | Evaluate a proposition with a given term for the variable.
-normal :: Term -> Prop VarP -> Prop PosP
-normal t p = case prop p of
-  NBin op t1    -> p { prop = Bin op t t1 }
-  NDivides n t1 -> p { prop = Divides n (t + t1) }
-
-
---------------------------------------------------------------------------------
-
--- | The meanings of atomic propositions
-eval_prop :: Prop PosP -> Env -> Bool
-eval_prop (Prop neg p) env = if neg then not res else res
-  where res = case p of
-                FF -> False
-                Divides n t  -> divides n (eval_term t env)
-                Bin op t1 t2 -> bin_op op (eval_term t1 env) (eval_term t2 env)
-                  
-
-bin_op :: RelOp -> Integer -> Integer -> Bool
-bin_op op x y = case op of
-                  Equal         -> x == y
-                  LessThan      -> x < y
-                  LessThanEqual -> x <= y
-
--- | Replace a variable with a constant, in a property.
-subst_prop :: Name -> Integer -> Prop PosP -> Prop PosP
-subst_prop x n (Prop b p) =
-  case p of
-    Bin op t1 t2 -> Prop b (Bin op (subst_term x n t1) (subst_term x n t2))
-    Divides k t  -> Prop b (Divides k (subst_term x n t))
-    FF           -> Prop b FF
-
-simplify_prop :: Prop PosP -> Prop PosP
-simplify_prop it@(Prop b p) =
-  case p of
-    Divides n t -> case is_constant t of
-                      Just v -> Prop (b /= divides n v) FF
-                      Nothing -> it
-    Bin Equal t1 t2 | not b && t1 == t2 -> Prop True FF
-    Bin op t1 t2 -> case (is_constant t1, is_constant t2) of
-                      (Just v1, Just v2) -> Prop (b /= bin_op op v1 v2) FF
-                      _ -> it
-    FF -> it
-
---------------------------------------------------------------------------------
-
-class SignPP t where
-  pp_neg :: Bool -> t -> Doc
-
-
-instance SignPP RelOp where
-
-  pp_neg False r = case r of
-    Equal         -> text "=="
-    LessThan      -> text "<"
-    LessThanEqual -> text "<="
-
-  pp_neg True r = case r of
-    Equal         -> text "/="
-    LessThan      -> text ">="
-    LessThanEqual -> text ">"
-
-
-pp_neg_div :: Bool -> Doc
-pp_neg_div False  = text "|"
-pp_neg_div True   = text "/|"
-
-
-instance SignPP PosP where
-  pp_neg n (Bin op t1 t2) = pp t1         <+> pp_neg n op  <+> pp t2
-  pp_neg n (Divides d t)  = text (show d) <+> pp_neg_div n <+> pp t
-  pp_neg n FF             = text (if n then "True" else "False")
-
-
-instance SignPP VarP where
-  pp_neg n (NBin op t)    = text "_" <+> pp_neg n op  <+> pp t
-  pp_neg n (NDivides d t) = text (show d) <+> pp_neg_div n
-                                          <+> text "_ +" <+> pp t
-
-
-instance SignPP CVarP where
-  pp_neg n p = case pprop p of
-    NBin op t     -> it <+> pp_neg n op  <+> pp t
-    NDivides d t  -> text (show d) <+> pp_neg_div n
-                                   <+> it <+> text "+" <+> pp t
-    where it  | c == 1    = text "_"
-              | c == (-1) = text "- _"
-              | otherwise = text (show c) <+> text "* _"
-
-          c = coeff p 
-               
-
-instance SignPP p => PP (Prop p) where
-  pp p  = pp_neg (negated p) (prop p)
-
-
-instance SignPP p => PP (NormProp p) where
-  pp (Ind p)  = pp p
-  pp (Norm p) = pp p
-
diff --git a/src/Data/Integer/Presburger/SolveDiv.hs b/src/Data/Integer/Presburger/SolveDiv.hs
deleted file mode 100644
--- a/src/Data/Integer/Presburger/SolveDiv.hs
+++ /dev/null
@@ -1,100 +0,0 @@
-module Data.Integer.Presburger.SolveDiv
-  ( DivCtr(..), Env, elim
-  ) where
-
-import Data.Integer.Presburger.Term
-import Data.List(foldl')
-
-
--- | A general "divisible by" constraint.
-data DivCtr     = Divs !Integer !Term
-
-
--- | Given some variables with bounds on them, and a set of
--- "divisible by" constraints, we produce all possible assignments
--- to the variables that are in bounds, and satisfy the constraints.
-elim :: Env -> [(Name,Integer)] -> [DivCtr] -> [ Env ]
-elim env0 [] ts = if all chk ts then [ env0 ] else []
-  where chk (Divs x t) = divides x (eval_term t env0)
-elim env0 ((x,bnd):xs) cs = do let (mb,cs1) = elim_var x cs
-                               env <- elim env0 xs cs1
-                               v <- case mb of
-                                      Nothing -> [ 1 .. bnd ]
-                                      Just (NDivides a b t) ->
-                                        theorem2 bnd (a,b,eval_term t env)
-                               return (env_extend x v env)
-
-
-
--- | "divisible by" constraint on a variable with a coefficient.
-data VarDivCtr  = NDivides { divisor  :: !Integer
-                           , coeff    :: !Integer
-                           , rest     :: !Term
-                           }
-
-
--- | This theorem combines two "divisible by" contratints on a single
--- variable, into a single constraint on the variable, and a generic
--- "divisible by" constraint that does not mention the variable.
-theorem1 :: VarDivCtr -> VarDivCtr -> (VarDivCtr, DivCtr)
-theorem1 NDivides { divisor = m, coeff = a1, rest = b1 }
-         NDivides { divisor = n, coeff = a2, rest = b2 }
-  = (new_x, new_other)
-
-  where (p,q,d)   = extended_gcd (a1 * n) (a2 * m)
-
-        new_x     = NDivides { divisor = m * n
-                             , coeff   = d
-                             , rest    = (p * n) .* b1 + (q * m) .* b2
-                             }
-
-        new_other = Divs d (a2 .* b1 - a1 .* b2)
-
-
--- | Repeatedly apply theorem 1 to a set of constraints,
--- to split them into a single constraint on the variable,
--- and additional constraints that do not mention the varibale.
-elim_var :: Name -> [DivCtr] -> (Maybe VarDivCtr, [DivCtr])
-elim_var x cs = case foldl' part ([],[]) cs of
-                  ([], have_not)     -> (Nothing, have_not)
-                  (h : hs, have_not) -> let (c,hn) = step h hs have_not
-                                        in (Just c,hn)
-
-  where part s@(have,have_not) c@(Divs m t)
-          | m == 1      = s -- ignore "divisible by 1" constraints.
-          | a == 0      = (have                 , c : have_not)
-          | otherwise   = (NDivides m a b : have,     have_not)
-            where (a,b) = split_term x t  -- t = a * x + b
-
-        step :: VarDivCtr -> [VarDivCtr] -> [DivCtr] -> (VarDivCtr,[DivCtr])
-        step h [] ns      = (h,ns)
-        step h (h1:hs) ns = step h2 hs (n : ns)
-          where (h2,n) = theorem1 h h1
-
-
-
--- | This theorem produces the solutions for a "divisible by" constraint
--- on a variable, where the "rest" term is a constant.
--- We peoduce only the solutions that are in the range [1 .. bnd]
---
--- solutions for x in [1 .. bnd] of: m | x * a + b
-theorem2 :: Integer -> (Integer,Integer,Integer) -> [Integer]
-theorem2 bnd (m,a,b)
-  | r == 0      = [ t * k - c | t <- [ lower .. upper ] ]
-  | otherwise   = []
-  where k           = div m d
-        c           = p * qu
-        (p,_,d)     = extended_gcd a m
-        (qu,r)      = divMod b d
-
-        (lower',r1) = divMod (1 + c) k
-        lower       = if r1 == 0 then lower' else lower' + 1  -- hmm
-        upper       = div (bnd + c) k
-
-  -- lower and upper:
-  -- t * k - c = 1   --> t = (1 + c) / k
-  -- t * k - c = bnd --> t = (bnd + c) / k
-
-
-
-
diff --git a/src/Data/Integer/Presburger/Term.hs b/src/Data/Integer/Presburger/Term.hs
deleted file mode 100644
--- a/src/Data/Integer/Presburger/Term.hs
+++ /dev/null
@@ -1,142 +0,0 @@
-module Data.Integer.Presburger.Term
-  ( Term, Name, split_term, is_constant, (.*), var, num
-  , Env, env_empty, env_extend
-  , eval_term, subst_term
-  , var_name
-  , module U
-  ) where
-
-import Data.Integer.Presburger.Utils as U
-
-import qualified Data.IntMap as Map
-import Data.Maybe(fromMaybe)
-import Control.Monad(mplus,guard)
-
-
--- | We represent the names of variables in terms as integers.
-type Name           = Int
-
--- | Terms of Presburger arithmetic.
--- Term are created by using the 'Num' class.
--- WARNING: Presburger arithmetic only supports multiplication
--- by a constant, trying to create invalid terms will result
--- in a run-time error.  A more type-safe alternative is to
--- use the '(.*)' operator.
-data Term           = Term (Map.IntMap Integer) Integer
-
-
--- | @split_term x (n * x + t1) = (n,t1)@
--- @x@ does not occur in @t1@
-split_term         :: Name -> Term -> (Integer,Term)
-split_term x (Term m n) = (fromMaybe 0 c, Term m1 n)
-  where (c,m1) = Map.updateLookupWithKey (\_ _ -> Nothing) x m
-
-var                :: Name -> Term
-var x               = Term (Map.singleton x 1) 0
-
-num                :: Integer -> Term
-num n               = Term Map.empty n
-
-
--- Evaluation ------------------------------------------------------------------
-newtype Env = Env (Map.IntMap Integer)
-
-env_empty :: Env
-env_empty = Env (Map.empty)
-
-env_extend :: Name -> Integer -> Env -> Env
-env_extend x v (Env m) = Env (Map.insert x v m)
-
--- The meaning of a term with free variables
--- If the term contains free variables that are not defined, then
--- we assume that these variables are 0.
-eval_term :: Term -> Env -> Integer
-eval_term (Term m k) (Env env) = sum (k : map eval_var (Map.toList m))
-  where eval_var (x,c) = case Map.lookup x env of
-                           Nothing -> 0
-                           Just v  -> c * v
-
-subst_term :: Name -> Integer -> Term -> Term
-subst_term x n t = case split_term x t of
-                     (c, Term m k) -> Term m (k + c * n)
-
---------------------------------------------------------------------------------
-
-instance Eq Term where
-  t1 == t2  = is_constant (t1 - t2) == Just 0
-
-instance Num Term where
-  fromInteger n             = Term Map.empty n
-
-  Term m1 n1 + Term m2 n2   = Term (Map.unionWith (+) m1 m2) (n1 + n2)
-
-  negate (Term m n)         = Term (Map.map negate m) (negate n)
-
-  t1 * t2  = case fmap (.* t2) (is_constant t1) `mplus`
-                  fmap (.* t1) (is_constant t2) of
-               Just t  -> t
-               Nothing -> error $ unlines [ "[(*) @ Term] Non-linear product:"
-                                          , "  *** " ++ show t1
-                                          , "  *** " ++ show t2
-                                          ]
-  signum t  = case is_constant t of
-                Just n  -> num (signum n)
-                Nothing -> error $ unlines [ "[signum @ Term]: Non-constant:"
-                                           , " *** " ++ show t
-                                           ]
-
-  abs t     = case is_constant t of
-                Just n  -> num (abs n)
-                Nothing -> error $ unlines [ "[abs @ Term]: Non-constant:"
-                                           , " *** " ++ show t
-                                           ]
-
-
--- | Check if a term is a constant (i.e., contains no variables).
--- If so, then we return the constant, otherwise we return 'Nothing'.
-is_constant :: Term -> Maybe Integer
-is_constant (Term m n) = guard (all (0 ==) (Map.elems m)) >> return n
-
-(.*) :: Integer -> Term -> Term
-0 .* _        = 0
-1 .* t        = t
-k .* Term m n = Term (Map.map (k *) m) (k * n)
-
-
-var_name           :: Name -> String
-var_name x          = let (a,b) = divMod x 26
-                          rest = if a == 0 then "" else show a
-                      in toEnum (97 + b) : rest
-
-instance Show Term where show x = show (pp x)
-instance PP Term where
-  pp (Term m k) | isEmpty vars  = text (show k)
-                | k == 0        = vars
-                | k > 0         = vars <+> char '+' <+> text (show k)
-                | otherwise     = vars <+> char '-' <+> text (show $ abs k)
-    where ppvar (x,n) = sign <+> co <+> text (var_name x)
-            where (sign,co)
-                     | n == -1    = (char '-', empty)
-                     | n < 0      = (char '-', text (show (abs n)) <+> char '*')
-                     | n == 1     = (char '+', empty)
-                     | otherwise  = (char '+', text (show n) <+> char '*')
-          first_var (x,1)  = text (var_name x)
-          first_var (x,-1) = char '-' <> text (var_name x)
-          first_var (x,n)  = text (show n) <+> char '*' <+> text (var_name x)
-
-          vars = case filter ((/= 0) . snd) (Map.toList m) of
-                   []     -> empty
-                   v : vs -> first_var v <+> hsep (map ppvar vs)
-
-
-instance PP Env where
-  pp (Env e)  = vcat (map sh (Map.toList e))
-    where sh (x,y)  = text (var_name x) <+> text "=" <+> text (show y)
-
-
-
-
-
-
-
-
diff --git a/src/Data/Integer/Presburger/Utils.hs b/src/Data/Integer/Presburger/Utils.hs
deleted file mode 100644
--- a/src/Data/Integer/Presburger/Utils.hs
+++ /dev/null
@@ -1,45 +0,0 @@
-module Data.Integer.Presburger.Utils
-  ( module Data.Integer.Presburger.Utils
-  , module PP
-  ) where
-
-import Text.PrettyPrint.HughesPJ as PP
-
-
-
-
-lcms :: Integral a => [a] -> a
-lcms xs = foldr lcm 1 xs
-
-
-groupEither :: [Either a b] -> ([a],[b])
-groupEither xs = foldr cons ([],[]) xs
-  where cons (Left a)  (as,bs) = (a:as,bs)
-        cons (Right b) (as,bs) = (as,b:bs)
-
-mapEither :: (a -> Either x y) -> [a] -> ([x],[y])
-mapEither f xs = groupEither (map f xs)
-
-
--- | let (p,q,r) = extended_gcd x y
---   in (x * p + y * q = r)  &&  (gcd x y = r)
-extended_gcd :: Integral a => a -> a -> (a,a,a)
-extended_gcd arg1 arg2 = loop arg1 arg2 0 1 1 0
-  where loop a b x lastx y lasty
-          | b /= 0    = let (q,b') = divMod a b
-                            x'     = lastx - q * x
-                            y'     = lasty - q * y
-                        in x' `seq` y' `seq` loop b b' x' x y' y
-          | otherwise = (lastx,lasty,a)
-
-
--- We define: "d | a" as "exists y. d * y = a"
-divides :: Integral a => a -> a -> Bool
-0 `divides` 0 = True
-0 `divides` _ = False
-x `divides` y = mod y x == 0
-
-
-class PP a where
-  pp :: a -> Doc
-
diff --git a/src/Data/Integer/SAT.hs b/src/Data/Integer/SAT.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Integer/SAT.hs
@@ -0,0 +1,822 @@
+{-# LANGUAGE Safe, PatternGuards #-}
+{-|
+This module implements a decision procedure for quantifier-free linear
+arithmetic.  The algorithm is based on the following paper:
+
+  An Online Proof-Producing Decision Procedure for
+  Mixed-Integer Linear Arithmetic
+  by
+  Sergey Berezin, Vijay Ganesh, and David L. Dill
+-}
+module Data.Integer.SAT
+  ( PropSet
+  , noProps
+  , checkSat
+  , assert
+  , Prop(..)
+  , Expr(..)
+  , BoundType(..)
+  , getExprBound
+  , getExprRange
+  , Name
+  , toName
+  , fromName
+  ) where
+
+import           Data.Map (Map)
+import qualified Data.Map as Map
+import           Data.List(partition)
+import           Data.Maybe(maybeToList,fromMaybe,mapMaybe)
+import           Control.Applicative(Applicative(..), (<$>))
+import           Control.Monad(liftM,ap,MonadPlus(..),msum,guard)
+import           Text.PrettyPrint(Doc,(<+>), (<>), integer, int, hsep, text)
+
+infixr 2 :||
+infixr 3 :&&
+infix  4 :==, :/=, :<, :<=, :>, :>=
+infixl 6 :+, :-
+infixl 7 :*
+
+--------------------------------------------------------------------------------
+-- Solver interface
+
+-- | A collection of propositions.
+newtype PropSet = State (Answer RW)
+                  deriving Show
+
+-- | An empty collection of propositions.
+noProps :: PropSet
+noProps = State $ return initRW
+
+-- | Add a new proposition to an existing collection.
+assert :: Prop -> PropSet -> PropSet
+assert p (State rws) = State $ fmap snd $ m =<< rws
+  where S m = prop p
+
+-- | Extract a model from a consistent set of propositions.
+-- Returns 'Nothing' if the assertions have no model.
+-- If a variable does not appear in the assignment, then it is 0 (?).
+checkSat :: PropSet -> Maybe [(Int,Integer)]
+checkSat (State m) = go m
+  where
+  go None            = mzero
+  go (One rw)        = return [ (x,v) | (UserName x, v) <- iModel (inerts rw) ]
+  go (Choice m1 m2)  = mplus (go m1) (go m2)
+
+-- | Computes bounds on the expression that are compatible with the model.
+-- Returns `Nothing` if the bound is not known.
+getExprBound :: BoundType -> Expr -> PropSet -> Maybe Integer
+getExprBound bt e (State s) =
+  do let S m          = expr e
+         check (t,s1) = iTermBound bt t (inerts s1)
+     bs <- mapM check $ toList $ s >>= m
+     case bs of
+       [] -> Nothing
+       _  -> Just (maximum bs)
+
+-- | Compute the range of possible values for an expression.
+-- Returns `Nothing` if the bound is not known.
+getExprRange :: Expr -> PropSet -> Maybe [Integer]
+getExprRange e (State s) =
+  do let S m          = expr e
+         check (t,s1) = do l <- iTermBound Lower t (inerts s1)
+                           u <- iTermBound Upper t (inerts s1)
+                           return (l,u)
+     bs <- mapM check $ toList $ s >>= m
+     case bs of
+       [] -> Nothing
+       _  -> let (ls,us) = unzip bs
+             in Just [ x | x <- [ minimum ls .. maximum us ] ]
+
+
+
+-- | The type of proposition.
+data Prop = PTrue
+          | PFalse
+          | Prop :|| Prop
+          | Prop :&& Prop
+          | Not Prop
+          | Expr :== Expr
+          | Expr :/= Expr
+          | Expr :<  Expr
+          | Expr :>  Expr
+          | Expr :<= Expr
+          | Expr :>= Expr
+            deriving (Read,Show)
+
+-- | The type of integer expressions.
+-- Variable names must be non-negative.
+data Expr = Expr :+ Expr          -- ^ Addition
+          | Expr :- Expr          -- ^ Subtraction
+          | Integer :* Expr       -- ^ Multiplication by a constant
+          | Negate Expr           -- ^ Negation
+          | Var Name              -- ^ Variable
+          | K Integer             -- ^ Constant
+          | If Prop Expr Expr     -- ^ A conditional expression
+          | Div Expr Integer      -- ^ Division, rounds down
+          | Mod Expr Integer      -- ^ Non-negative remainder
+            deriving (Read,Show)
+
+prop :: Prop -> S ()
+prop PTrue       = return ()
+prop PFalse      = mzero
+prop (p1 :|| p2) = prop p1 `mplus` prop p2
+prop (p1 :&& p2) = prop p1 >> prop p2
+prop (Not p)     = prop (neg p)
+  where
+  neg PTrue       = PFalse
+  neg PFalse      = PTrue
+  neg (p1 :&& p2) = neg p1 :|| neg p2
+  neg (p1 :|| p2) = neg p1 :&& neg p2
+  neg (Not q)     = q
+  neg (e1 :== e2) = e1 :/= e2
+  neg (e1 :/= e2) = e1 :== e2
+  neg (e1 :<  e2) = e1 :>= e2
+  neg (e1 :<= e2) = e1 :>  e2
+  neg (e1 :>  e2) = e1 :<= e2
+  neg (e1 :>= e2) = e1 :<  e2
+
+prop (e1 :== e2) = do t1 <- expr e1
+                      t2 <- expr e2
+                      enqAndGo qZeroTerms (t1 |-| t2)
+
+prop (e1 :/= e2)  = do t1 <- expr e1
+                       t2 <- expr e2
+                       let t = t1 |-| t2
+                       enqAndGo qNegTerms t `mplus` enqAndGo qNegTerms (tNeg t)
+
+prop (e1 :< e2)   = do t1 <- expr e1
+                       t2 <- expr e2
+                       enqAndGo qNegTerms (t1 |-| t2)
+
+prop (e1 :<= e2)  = do t1 <- expr e1
+                       t2 <- expr e2
+                       let t = t1 |-| t2
+                       enqAndGo qZeroTerms t `mplus` enqAndGo qNegTerms t
+
+prop (e1 :> e2)   = prop (e2 :<  e1)
+prop (e1 :>= e2)  = prop (e2 :<= e1)
+
+
+expr :: Expr -> S Term
+expr (e1 :+ e2)   = (|+|)   <$> expr e1 <*> expr e2
+expr (e1 :- e2)   = (|-|)   <$> expr e1 <*> expr e2
+expr (k  :* e2)   = (k |*|) <$> expr e2
+expr (Negate e)   = tNeg    <$> expr e
+expr (Var x)      = pure (tVar x)
+expr (K x)        = pure (tConst x)
+expr (If p e1 e2) = do x <- newVar
+                       prop (p :&& Var x :== e1 :|| Not p :&& Var x :== e2)
+                       return (tVar x)
+expr (Div e k)    = fmap fst $ exprDivMod e k
+expr (Mod e k)    = fmap snd $ exprDivMod e k
+
+exprDivMod :: Expr -> Integer -> S (Term,Term)
+exprDivMod e k =
+  do guard (k /= 0) -- Always unsat
+     q <- newVar
+     r <- newVar
+     let er = Var r
+     prop (k :* Var q :+ er :== e :&& er :< K k :&& K 0 :<= er)
+     return (tVar q, tVar r)
+
+
+
+
+
+--------------------------------------------------------------------------------
+
+data RW = RW { nameSource :: !Int
+             , todo       :: WorkQ
+             , inerts     :: Inerts
+             } deriving Show
+
+initRW :: RW
+initRW = RW { nameSource = 0, todo = qEmpty, inerts = iNone }
+
+solveAll :: S ()
+solveAll =
+  do mbEq <- getWork qZeroTerms
+     case mbEq of
+       Just p  -> solveIs0 p >> solveAll
+       Nothing ->
+         do mbLt <- getWork qNegTerms
+            case mbLt of
+              Just p  -> solveIsNeg p >> solveAll
+              Nothing -> return ()
+
+
+--------------------------------------------------------------------------------
+-- The work queue
+
+data WorkQ = WorkQ { zeroTerms     :: [Term]    -- ^ t == 0
+                   , negTerms      :: [Term]    -- ^ t <  0
+                   } deriving Show
+
+qEmpty :: WorkQ
+qEmpty = WorkQ { zeroTerms = [], negTerms = [] }
+
+qLet :: Name -> Term -> WorkQ -> WorkQ
+qLet x t q = WorkQ { zeroTerms      = map (tLet x t) (zeroTerms q)
+                   , negTerms       = map (tLet x t) (negTerms  q)
+                   }
+
+type Field t = (WorkQ -> [t], [t] -> WorkQ -> WorkQ)
+
+qZeroTerms :: Field Term
+qZeroTerms = (zeroTerms, \a q -> q { zeroTerms = a })
+
+qNegTerms :: Field Term
+qNegTerms = (negTerms, \a q -> q { negTerms = a })
+
+--------------------------------------------------------------------------------
+-- Constraints and Bound on Variables
+
+ctLt :: Term -> Term -> Term
+ctLt t1 t2 = t1 |-| t2
+
+ctEq :: Term -> Term -> Term
+ctEq t1 t2 = t1 |-| t2
+
+data Bound      = Bound Integer Term  -- ^ The integer is strictly positive
+                  deriving Show
+
+data BoundType  = Lower | Upper
+                  deriving Show
+
+toCt :: BoundType -> Name -> Bound -> Term
+toCt Lower x (Bound c t) = ctLt t              (c |*| tVar x)
+toCt Upper x (Bound c t) = ctLt (c |*| tVar x) t
+
+
+
+--------------------------------------------------------------------------------
+-- Inert set
+
+-- | The inert contains the solver state on one possible path.
+data Inerts = Inerts
+  { bounds :: NameMap ([Bound],[Bound])
+    -- ^ Known lower and upper bounds for variables.
+    -- Each bound @(c,t)@ in the first list asserts that  @t < c * x@
+    -- Each bound @(c,t)@ in the second list asserts that @c * x < t@
+
+  , solved :: NameMap Term
+    -- ^ Definitions for resolved variabless.
+    -- These form an idempotent substitution.
+  } deriving Show
+
+
+-- | An empty inert set.
+iNone :: Inerts
+iNone = Inerts { bounds = Map.empty
+               , solved = Map.empty
+               }
+
+-- | Rewrite a term using the definitions from an inert set.
+iApSubst :: Inerts -> Term -> Term
+iApSubst i t = foldr apS t $ Map.toList $ solved i
+  where apS (x,t1) t2 = tLet x t1 t2
+
+-- | Add a definition.  Upper and lower bound constraints that mention
+-- the variable are "kicked-out" so that they can be reinserted in the
+-- context of the new knowledge.
+--
+--    * Assumes substitution has already been applied.
+--
+--    * The kciked-out constraints are NOT rewritten, this happens
+--      when they get inserted in the work queue.
+
+iSolved :: Name -> Term -> Inerts -> ([Term], Inerts)
+iSolved x t i =
+  ( kickedOut
+  , Inerts { bounds = otherBounds
+           , solved = Map.insert x t $ Map.map (tLet x t) $ solved i
+           }
+  )
+  where
+  (kickedOut, otherBounds) =
+
+        -- First, we eliminate all entries for `x`
+    let (mb, mp1) = Map.updateLookupWithKey (\_ _ -> Nothing) x (bounds i)
+
+        -- Next, we elminate all constraints that mentiond `x` in bounds
+        mp2 = Map.mapWithKey extractBounds mp1
+
+    in ( [ ct | (lbs,ubs) <- maybeToList mb
+              ,  ct <- map (toCt Lower x) lbs ++ map (toCt Upper x) ubs ]
+         ++
+         [ ct | (_,cts) <- Map.elems mp2, ct <- cts ]
+
+       , fmap fst mp2
+       )
+
+  extractBounds y (lbs,ubs) =
+    let (lbsStay, lbsKick) = partition stay lbs
+        (ubsStay, ubsKick) = partition stay ubs
+    in ( (lbsStay,ubsStay)
+       , map (toCt Lower y) lbsKick ++
+         map (toCt Upper y) ubsKick
+       )
+
+  stay (Bound _ bnd) = not (tHasVar x bnd)
+
+
+-- Given a list of lower (resp. upper) bounds, compute the least (resp. largest)
+-- value that satisfies them all.
+iPickBounded :: BoundType -> [Bound] -> Maybe Integer
+iPickBounded bt bs = go bs Nothing
+  where
+  go [] mb = mb
+  go (Bound c t : more) mb =
+    do k <- isConst t
+       let t1 = maybe k (combine k) mb
+       go more $ Just $ compute t1 c
+
+  combine = case bt of
+              Lower -> max
+              Upper -> min
+
+  compute v c = case bt of
+                  Lower -> div v c + 1
+                  Upper -> let (q,r) = divMod v c
+                           in if r == 0 then q - 1 else q
+
+
+-- | The largest (resp. least) upper (resp. lower) bound on a term
+-- that will satisfy the model
+iTermBound :: BoundType -> Term -> Inerts -> Maybe Integer
+iTermBound bt (T k xs) is = do ks <- mapM summand (Map.toList xs)
+                               return $ sum $ k : ks
+  where
+  summand (x,c) = fmap (c *) (iVarBound (newBt c) x is)
+  newBt c = if c > 0 then bt else case bt of
+                                    Lower -> Upper
+                                    Upper -> Lower
+
+
+-- | The largest (resp. least) upper (resp. lower) bound on a variable
+-- that will satisfy the model.
+iVarBound :: BoundType -> Name -> Inerts -> Maybe Integer
+iVarBound bt x is
+  | Just t <- Map.lookup x (solved is) = iTermBound bt t is
+
+iVarBound bt x is =
+  do both <- Map.lookup x (bounds is)
+     case mapMaybe fromBound (chooseBounds both) of
+       [] -> Nothing
+       bs -> return (combineBounds bs)
+  where
+  fromBound (Bound c t) = fmap (scaleBound c) (iTermBound bt t is)
+
+  combineBounds = case bt of
+                    Upper -> minimum
+                    Lower -> maximum
+
+  chooseBounds = case bt of
+                   Upper -> snd
+                   Lower -> fst
+
+  scaleBound c b = case bt of
+                     Upper -> div (b-1) c
+                     Lower -> div b c + 1
+
+
+
+iModel :: Inerts -> [(Name,Integer)]
+iModel i = goBounds [] (bounds i)
+  where
+  goBounds su mp =
+    case Map.maxViewWithKey mp of
+      Nothing -> goEqs su $ Map.toList $ solved i
+      Just ((x,(lbs0,ubs0)), mp1) ->
+        let lbs = [ Bound c (tLetNums su t) | Bound c t <- lbs0 ]
+            ubs = [ Bound c (tLetNums su t) | Bound c t <- ubs0 ]
+            sln = fromMaybe 0
+                $ mplus (iPickBounded Lower lbs) (iPickBounded Upper ubs)
+        in goBounds ((x,sln) : su) mp1
+
+  goEqs su [] = su
+  goEqs su ((x,t) : more) =
+    let t1  = tLetNums su t
+        vs  = tVarList t1
+        su1 = [ (v,0) | v <- vs ] ++ (x,tConstPart t1) : su
+    in goEqs su1 more
+
+
+--------------------------------------------------------------------------------
+-- Solving constraints
+
+-- | Solve a constraint if the form @t = 0@.
+-- Assumes substitution has already been applied.
+solveIs0 :: Term -> S ()
+solveIs0 t
+
+  -- A == 0
+  | Just a <- isConst t = guard (a == 0)
+
+  -- A + B * x = 0
+  | Just (a,b,x) <- tIsOneVar t =
+    case divMod (-a) b of
+      (q,0) -> addDef x (tConst q)
+      _     -> mzero
+
+  -- x + S = 0
+  | Just (xc,x,s) <- tGetSimpleCoeff t =
+    addDef x (if xc > 0 then tNeg s else s)
+
+  -- A * S = 0
+  | Just (_, s) <- tFactor t  = addWork qZeroTerms s
+
+  -- See Section 3.1 of paper for details.
+  -- We obtain an equivalent formulation but with smaller coefficients.
+  | Just (ak,xk,s) <- tLeastAbsCoeff t =
+      do let m = abs ak + 1
+         v <- newVar
+         let sgn  = signum ak
+             soln =     (negate sgn * m) |*| tVar v
+                    |+| tMapCoeff (\c -> sgn * modulus c m) s
+         addDef xk soln
+
+         let upd i = div (2*i + m) (2*m) + modulus i m
+         addWork qZeroTerms (negate (abs ak) |*| tVar v |+| tMapCoeff upd s)
+
+  | otherwise = error "solveIs0: unreachable"
+
+modulus :: Integer -> Integer -> Integer
+modulus a m = a - m * div (2 * a + m) (2 * m)
+
+
+-- | Solve a constraint of the form @t < 0@.
+-- Assumes that substitution has been applied
+solveIsNeg :: Term -> S ()
+solveIsNeg t
+
+  -- A < 0
+  | Just a <- isConst t = guard (a < 0)
+
+  -- A * S < 0
+  |Just (_,s) <- tFactor t = addWork qNegTerms s
+
+  -- See Section 5.1 of the paper
+  | Just (xc,x,s) <- tLeastVar t =
+
+    do ctrs <- if xc < 0
+               -- -XC*x + S < 0
+               -- S < XC*x
+               then do ubs <- getBounds Upper x
+                       let b    = negate xc
+                           beta = s
+                       addBound Lower x (Bound b beta)
+                       return [ (a,alpha,b,beta) | Bound a alpha <- ubs ]
+               -- XC*x + S < 0
+               -- XC*x < -S
+               else do lbs <- getBounds Lower x
+                       let a     = xc
+                           alpha = tNeg s
+                       addBound Upper x (Bound a alpha)
+                       return [ (a,alpha,b,beta) | Bound b beta <- lbs ]
+
+      -- See Note [Shadows]
+       mapM_ (\(a,alpha,b,beta) ->
+          do let real = ctLt (a |*| beta) (b |*| alpha)
+                 dark = ctLt (tConst (a * b)) (b |*| alpha |-| a |*| beta)
+                 gray = [ ctEq (b |*| tVar x) (tConst i |+| beta)
+                                                      | i <- [ 1 .. b - 1 ] ]
+             addWork qNegTerms real
+             msum (addWork qNegTerms dark : map (addWork qZeroTerms) gray)
+             ) ctrs
+
+  | otherwise = error "solveIsNeg: unreachable"
+
+
+{- Note [Shadows]
+
+  P: beta < b * x
+  Q: a * x < alpha
+
+real: a * beta < b * alpha
+
+  beta     < b * x      -- from P
+  a * beta < a * b * x  -- (a *)
+  a * beta < b * alpha  -- comm. and Q
+
+
+dark: b * alpha - a * beta > a * b
+
+
+gray: b * x = beta + 1 \/
+      b * x = beta + 2 \/
+      ...
+      b * x = beta + (b-1)
+
+We stop at @b - 1@ because if:
+
+> b * x                >= beta + b
+> a * b * x            >= a * (beta + b)     -- (a *)
+> a * b * x            >= a * beta + a * b   -- distrib.
+> b * alpha            >  a * beta + a * b   -- comm. and Q
+> b * alpha - a * beta > a * b               -- subtract (a * beta)
+
+which is covered by the dark shadow.
+-}
+
+
+--------------------------------------------------------------------------------
+-- Monads
+
+data Answer a = None | One a | Choice (Answer a) (Answer a)
+                deriving Show
+
+toList :: Answer a -> [a]
+toList a = go a []
+  where
+  go (Choice xs ys) zs = go xs (go ys zs)
+  go (One x) xs        = x : xs
+  go None xs           = xs
+
+
+instance Monad Answer where
+  return a           = One a
+  fail _             = None
+  None >>= _         = None
+  One a >>= k        = k a
+  Choice m1 m2 >>= k = mplus (m1 >>= k) (m2 >>= k)
+
+instance MonadPlus Answer where
+  mzero                = None
+  mplus None x         = x
+  -- mplus (Choice x y) z = mplus x (mplus y z)
+  mplus x y            = Choice x y
+
+instance Functor Answer where
+  fmap _ None           = None
+  fmap f (One x)        = One (f x)
+  fmap f (Choice x1 x2) = Choice (fmap f x1) (fmap f x2)
+
+instance Applicative Answer where
+  pure  = return
+  (<*>) = ap
+
+
+newtype S a = S (RW -> Answer (a,RW))
+
+instance Monad S where
+  return a      = S $ \s -> return (a,s)
+  S m >>= k     = S $ \s -> do (a,s1) <- m s
+                               let S m1 = k a
+                               m1 s1
+
+instance MonadPlus S where
+  mzero               = S $ \_ -> mzero
+  mplus (S m1) (S m2) = S $ \s -> mplus (m1 s) (m2 s)
+
+instance Functor S where
+  fmap = liftM
+
+instance Applicative S where
+  pure  = return
+  (<*>) = ap
+
+updS :: (RW -> (a,RW)) -> S a
+updS f = S $ \s -> return (f s)
+
+updS_ :: (RW -> RW) -> S ()
+updS_ f = updS $ \rw -> ((),f rw)
+
+get :: (RW -> a) -> S a
+get f = updS $ \rw -> (f rw, rw)
+
+newVar :: S Name
+newVar = updS $ \rw -> ( SysName (nameSource rw)
+                       , rw { nameSource = nameSource rw + 1 }
+                       )
+
+-- | Try to get a new item from the work queue.
+getWork :: Field t -> S (Maybe t)
+getWork (getF,setF) = updS $ \rw ->
+  let work = todo rw
+  in case getF work of
+       []     -> (Nothing, rw)
+       t : ts -> (Just t,  rw { todo = setF ts work })
+
+-- | Add a new item to the work queue.
+addWork :: Field t -> t -> S ()
+addWork (getF,setF) t = updS_ $ \rw ->
+  let work = todo rw
+  in rw { todo = setF (t : getF work) work }
+
+-- | Get lower ('fst'), or upper ('snd') bounds for a variable.
+getBounds :: BoundType -> Name -> S [Bound]
+getBounds f x = get $ \rw -> case Map.lookup x $ bounds $ inerts rw of
+                               Nothing -> []
+                               Just bs -> case f of
+                                            Lower -> fst bs
+                                            Upper -> snd bs
+
+addBound :: BoundType -> Name -> Bound -> S ()
+addBound bt x b = updS_ $ \rw ->
+  let i = inerts rw
+      entry = case bt of
+                Lower -> ([b],[])
+                Upper -> ([],[b])
+      jn (newL,newU) (oldL,oldU) = (newL++oldL, newU++oldU)
+  in rw { inerts = i { bounds = Map.insertWith jn x entry (bounds i) }}
+
+-- | Add a new definition.
+-- Assumes substitution has already been applied
+addDef :: Name -> Term -> S ()
+addDef x t = updS_ $ \rw ->
+  let (newWork,newInerts) = iSolved x t (inerts rw)
+  in rw { inerts = newInerts
+        , todo   = qLet x t $
+                     let work = todo rw
+                     in work { negTerms = newWork ++ negTerms work }
+        }
+
+enqAndGo :: Field Term -> Term -> S ()
+enqAndGo q t =
+  do i <- get inerts
+     addWork q $ iApSubst i t
+     solveAll
+
+
+
+
+--------------------------------------------------------------------------------
+
+
+data Name = UserName !Int | SysName !Int
+            deriving (Read,Show,Eq,Ord)
+
+ppName :: Name -> Doc
+ppName (UserName x) = text "u" <> int x
+ppName (SysName x)  = text "s" <> int x
+
+toName :: Int -> Name
+toName = UserName
+
+fromName :: Name -> Maybe Int
+fromName (UserName x) = Just x
+fromName (SysName _)  = Nothing
+
+
+
+
+type NameMap = Map Name
+
+-- | The type of terms.  The integer is the constant part of the term,
+-- and the `Map` maps variables (represented by @Int@ to their coefficients).
+-- The term is a sum of its parts.
+-- INVARIANT: the `Map` does not map anything to 0.
+data Term = T !Integer (NameMap Integer)
+              deriving (Eq,Ord)
+
+infixl 6 |+|, |-|
+infixr 7 |*|
+
+-- | A constant term.
+tConst :: Integer -> Term
+tConst k = T k Map.empty
+
+-- | Construct a term with a single variable.
+tVar :: Name -> Term
+tVar x = T 0 (Map.singleton x 1)
+
+(|+|) :: Term -> Term -> Term
+T n1 m1 |+| T n2 m2 = T (n1 + n2)
+                    $ if Map.null m1 then m2 else
+                      if Map.null m2 then m1 else
+                      Map.filter (/= 0) $ Map.unionWith (+) m1 m2
+
+(|*|) :: Integer -> Term -> Term
+0 |*| _     = tConst 0
+1 |*| t     = t
+k |*| T n m = T (k * n) (fmap (k *) m)
+
+tNeg :: Term -> Term
+tNeg t = (-1) |*| t
+
+(|-|) :: Term -> Term -> Term
+t1 |-| t2 = t1 |+| tNeg t2
+
+
+-- | Replace a variable with a term.
+tLet :: Name -> Term -> Term -> Term
+tLet x t1 t2 = let (a,t) = tSplitVar x t2
+               in a |*| t1 |+| t
+
+-- | Replace a variable with a constant.
+tLetNum :: Name -> Integer -> Term -> Term
+tLetNum x k t = let (c,T n m) = tSplitVar x t
+                in T (c * k + n) m
+
+-- | Replace the given variables with constants.
+tLetNums :: [(Name,Integer)] -> Term -> Term
+tLetNums xs t = foldr (\(x,i) t1 -> tLetNum x i t1) t xs
+
+
+
+
+instance Show Term where
+  showsPrec c t = showsPrec c (show (ppTerm t))
+
+ppTerm :: Term -> Doc
+ppTerm (T k m) =
+  case Map.toList m of
+    [] -> integer k
+    xs | k /= 0 -> hsep (integer k : map ppProd xs)
+    x : xs      -> hsep (ppFst x   : map ppProd xs)
+
+  where
+  ppFst (x,1)   = ppName x
+  ppFst (x,-1)  = text "-" <> ppName x
+  ppFst (x,n)   = ppMul n x
+
+  ppProd (x,1)  = text "+" <+> ppName x
+  ppProd (x,-1) = text "-" <+> ppName x
+  ppProd (x,n) | n > 0      = text "+" <+> ppMul n x
+               | otherwise  = text "-" <+> ppMul (abs n) x
+
+  ppMul n x = integer n <+> text "*" <+> ppName x
+
+-- | Remove a variable from the term, and return its coefficient.
+-- If the variable is not present in the term, the coefficient is 0.
+tSplitVar :: Name -> Term -> (Integer, Term)
+tSplitVar x t@(T n m) =
+  case Map.updateLookupWithKey (\_ _ -> Nothing) x m of
+    (Nothing,_) -> (0,t)
+    (Just k,m1) -> (k, T n m1)
+
+-- | Does the term contain this varibale?
+tHasVar :: Name -> Term -> Bool
+tHasVar x (T _ m) = Map.member x m
+
+-- | Is this terms just an integer.
+isConst :: Term -> Maybe Integer
+isConst (T n m)
+  | Map.null m  = Just n
+  | otherwise   = Nothing
+
+tConstPart :: Term -> Integer
+tConstPart (T n _) = n
+
+-- | Returns: @Just (a, b, x)@ if the term is the form: @a + b * x@
+tIsOneVar :: Term -> Maybe (Integer, Integer, Name)
+tIsOneVar (T a m) = case Map.toList m of
+                      [ (x,b) ] -> Just (a, b, x)
+                      _         -> Nothing
+
+-- | Spots terms that contain variables with unit coefficients
+-- (i.e., of the form @x + t@ or @t - x@).
+-- Returns (coeff, var, rest of term)
+tGetSimpleCoeff :: Term -> Maybe (Integer, Name, Term)
+tGetSimpleCoeff (T a m) =
+  do let (m1,m2) = Map.partition (\x -> x == 1 || x == -1) m
+     ((x,xc), m3) <- Map.minViewWithKey m1
+     return (xc, x, T a (Map.union m3 m2))
+
+tVarList :: Term -> [Name]
+tVarList (T _ m) = Map.keys m
+
+
+-- | Try to factor-out a common consant (> 1) from a term.
+-- For example, @2 + 4x@ becomes @2 * (1 + 2x)@.
+tFactor :: Term -> Maybe (Integer, Term)
+tFactor (T c m) =
+  do d <- common (c : Map.elems m)
+     return (d, T (div c d) (fmap (`div` d) m))
+  where
+  common :: [Integer] -> Maybe Integer
+  common []  = Nothing
+  common [x] = Just x
+  common (x : y : zs) =
+    case gcd x y of
+      1 -> Nothing
+      n -> common (n : zs)
+
+-- | Extract a variable with a coefficient whose absolute value is minimal.
+tLeastAbsCoeff :: Term -> Maybe (Integer, Name, Term)
+tLeastAbsCoeff (T c m) = do (xc,x,m1) <- Map.foldWithKey step Nothing m
+                            return (xc, x, T c m1)
+  where
+  step x xc Nothing   = Just (xc, x, Map.delete x m)
+  step x xc (Just (yc,_,_))
+    | abs xc < abs yc = Just (xc, x, Map.delete x m)
+  step _ _ it         = it
+
+-- | Extract the least variable from a term
+tLeastVar :: Term -> Maybe (Integer, Name, Term)
+tLeastVar (T c m) =
+  do ((x,xc), m1) <- Map.minViewWithKey m
+     return (xc, x, T c m1)
+
+-- | Apply a function to all coefficients, including the constnat
+tMapCoeff :: (Integer -> Integer) -> Term -> Term
+tMapCoeff f (T c m) = T (f c) (fmap f m)
+
+
+
+
+
+
+