diff --git a/presburger.cabal b/presburger.cabal
--- a/presburger.cabal
+++ b/presburger.cabal
@@ -1,5 +1,5 @@
 Name:           presburger
-Version:        1.1
+Version:        1.2
 License:        BSD3
 License-file:   LICENSE
 Author:         Iavor S. Diatchki
@@ -10,7 +10,7 @@
 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
+Cabal-version:  >= 1.8
 
 library
   Build-Depends:  base < 10, containers, pretty
@@ -23,4 +23,10 @@
 source-repository head
   type: git
   location: git://github.com/yav/presburger.git
+
+Test-Suite pressburger-qc-tests
+  type: exitcode-stdio-1.0
+  hs-source-dirs: tests
+  main-is: qc.hs
+  build-depends: base, presburger == 1.2, QuickCheck
 
diff --git a/src/Data/Integer/SAT.hs b/src/Data/Integer/SAT.hs
--- a/src/Data/Integer/SAT.hs
+++ b/src/Data/Integer/SAT.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE Safe, PatternGuards #-}
+{-# LANGUAGE Trustworthy, PatternGuards, BangPatterns #-}
 {-|
 This module implements a decision procedure for quantifier-free linear
 arithmetic.  The algorithm is based on the following paper:
@@ -21,15 +21,35 @@
   , Name
   , toName
   , fromName
+  -- * Iterators
+  , allSolutions
+  , slnCurrent
+  , slnNextVal
+  , slnNextVar
+  , slnEnumerate
+
+
+  -- * Debug
+  , dotPropSet
+  , sizePropSet
+  , allInerts
+  , ppInerts
+
+  -- * For QuickCheck
+  , iPickBounded
+  , Bound(..)
+  , tConst
   ) where
 
+import Debug.Trace
+
 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)
+import           Control.Monad(liftM,ap,MonadPlus(..),guard)
+import           Text.PrettyPrint
 
 infixr 2 :||
 infixr 3 :&&
@@ -44,6 +64,12 @@
 newtype PropSet = State (Answer RW)
                   deriving Show
 
+dotPropSet :: PropSet -> Doc
+dotPropSet (State a) = dotAnswer (ppInerts . inerts) a
+
+sizePropSet :: PropSet -> (Integer,Integer,Integer)
+sizePropSet (State a) = answerSize a
+
 -- | An empty collection of propositions.
 noProps :: PropSet
 noProps = State $ return initRW
@@ -63,6 +89,13 @@
   go (One rw)        = return [ (x,v) | (UserName x, v) <- iModel (inerts rw) ]
   go (Choice m1 m2)  = mplus (go m1) (go m2)
 
+allInerts :: PropSet -> [Inerts]
+allInerts (State m) = map inerts (toList m)
+
+allSolutions :: PropSet -> [Solutions]
+allSolutions = map startIter . allInerts
+
+
 -- | Computes bounds on the expression that are compatible with the model.
 -- Returns `Nothing` if the bound is not known.
 getExprBound :: BoundType -> Expr -> PropSet -> Maybe Integer
@@ -138,21 +171,21 @@
 
 prop (e1 :== e2) = do t1 <- expr e1
                       t2 <- expr e2
-                      enqAndGo qZeroTerms (t1 |-| t2)
+                      solveIs0 (t1 |-| t2)
 
 prop (e1 :/= e2)  = do t1 <- expr e1
                        t2 <- expr e2
                        let t = t1 |-| t2
-                       enqAndGo qNegTerms t `mplus` enqAndGo qNegTerms (tNeg t)
+                       solveIsNeg t `orElse` solveIsNeg (tNeg t)
 
 prop (e1 :< e2)   = do t1 <- expr e1
                        t2 <- expr e2
-                       enqAndGo qNegTerms (t1 |-| t2)
+                       solveIsNeg (t1 |-| t2)
 
 prop (e1 :<= e2)  = do t1 <- expr e1
                        t2 <- expr e2
-                       let t = t1 |-| t2
-                       enqAndGo qZeroTerms t `mplus` enqAndGo qNegTerms t
+                       let t = t1 |-| t2 |-| tConst 1
+                       solveIsNeg t
 
 prop (e1 :> e2)   = prop (e2 :<  e1)
 prop (e1 :>= e2)  = prop (e2 :<= e1)
@@ -187,47 +220,11 @@
 --------------------------------------------------------------------------------
 
 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 })
+initRW = RW { nameSource = 0, inerts = iNone }
 
 --------------------------------------------------------------------------------
 -- Constraints and Bound on Variables
@@ -265,7 +262,20 @@
     -- These form an idempotent substitution.
   } deriving Show
 
+ppInerts :: Inerts -> Doc
+ppInerts is = vcat $ [ ppLower x b | (x,(ls,_)) <- bnds, b <- ls ] ++
+                     [ ppUpper x b | (x,(_,us)) <- bnds, b <- us ] ++
+                     [ ppEq e      | e <- Map.toList (solved is) ]
+  where
+  bnds = Map.toList (bounds is)
 
+  ppT c x                = ppTerm (c |*| tVar x)
+  ppLower x (Bound c t)  = ppTerm t <+> text "<" <+> ppT c x
+  ppUpper x (Bound c t)  = ppT c x  <+> text "<" <+> ppTerm t
+  ppEq (x,t)             = ppName x <+> text "=" <+> ppTerm t
+
+
+
 -- | An empty inert set.
 iNone :: Inerts
 iNone = Inerts { bounds = Map.empty
@@ -321,6 +331,114 @@
   stay (Bound _ bnd) = not (tHasVar x bnd)
 
 
+-- | Given some lower and upper bounds, find the interval the satisfies them.
+-- Note the upper and lower bounds are strict (i.e., < and >)
+boundInterval :: [Bound] -> [Bound] -> Maybe (Maybe Integer, Maybe Integer)
+boundInterval lbs ubs =
+  do ls <- mapM (normBound Lower) lbs
+     us <- mapM (normBound Upper) ubs
+     let lb = case ls of
+                [] -> Nothing
+                _  -> Just (maximum ls + 1)
+         ub = case us of
+                [] -> Nothing
+                _  -> Just (minimum us - 1)
+     case (lb,ub) of
+       (Just l, Just u) -> guard (l <= u)
+       _                -> return ()
+     return (lb,ub)
+  where
+  normBound Lower (Bound c t) = do k <- isConst t
+                                   return (div (k + c - 1) c)
+  normBound Upper (Bound c t) = do k <- isConst t
+                                   return (div k c)
+
+data Solutions = Done
+               | TopVar Name Integer (Maybe Integer) (Maybe Integer) Inerts
+               | FixedVar Name Integer Solutions
+                  deriving Show
+
+slnCurrent :: Solutions -> [(Int,Integer)]
+slnCurrent s = [ (x,v) | (UserName x, v) <- go s ]
+  where
+  go Done                = []
+  go (TopVar x v _ _ is) = (x, v) : iModel (iLet x v is)
+  go (FixedVar x v i)    = (x, v) : go i
+
+-- | Replace occurances of a variable with an integer.
+-- WARNING: The integer should be a valid value for the variable.
+iLet :: Name -> Integer -> Inerts -> Inerts
+iLet x v is = Inerts { bounds = fmap updBs (bounds is)
+                     , solved = fmap (tLetNum x v) (solved is) }
+  where
+  updB (Bound c t) = Bound c (tLetNum x v t)
+  updBs (ls,us)    = (map updB ls, map updB us)
+
+
+startIter :: Inerts -> Solutions
+startIter is =
+  case Map.maxViewWithKey (bounds is) of
+    Nothing ->
+      case Map.maxViewWithKey (solved is) of
+        Nothing -> Done
+        Just ((x,t), mp1) ->
+          case [ y | y <- tVarList t ] of
+            y : _ -> TopVar y 0 Nothing Nothing is
+            [] -> let v = tConstPart t
+                  in TopVar x v (Just v) (Just v) $ is { solved = mp1 }
+    Just ((x,(lbs,ubs)), mp1) ->
+      case [ y | Bound _ t <- lbs ++ ubs, y <- tVarList t ] of
+        y : _ -> TopVar y 0 Nothing Nothing is
+        [] -> case boundInterval lbs ubs of
+                Nothing -> error "bug: cannot compute interval?"
+                Just (lb,ub) ->
+                  let v = fromMaybe 0 (mplus lb ub)
+                  in TopVar x v lb ub $ is { bounds = mp1 }
+
+slnEnumerate :: Solutions -> [ Solutions ]
+slnEnumerate s0 = go s0 []
+  where
+  go s k  = case slnNextVar s of
+              Nothing -> hor s k
+              Just s1 -> go s1 $ case slnNextVal s of
+                                   Nothing -> k
+                                   Just s2 -> go s2 k
+
+  hor s k = s
+          : case slnNextVal s of
+              Nothing -> k
+              Just s1 -> hor s1 k
+
+slnNextVal :: Solutions -> Maybe Solutions
+slnNextVal Done = Nothing
+slnNextVal (FixedVar x v i) = FixedVar x v `fmap` slnNextVal i
+slnNextVal it@(TopVar _ _ lb _ _) =
+  case lb of
+    Just _  -> slnNextValWith (+1) it
+    Nothing -> slnNextValWith (subtract 1) it
+
+
+slnNextValWith :: (Integer -> Integer) -> Solutions -> Maybe Solutions
+slnNextValWith _ Done = Nothing
+slnNextValWith f (FixedVar x v i) = FixedVar x v `fmap` slnNextValWith f i
+slnNextValWith f (TopVar x v lb ub is) =
+  do let v1 = f v
+     case lb of
+       Just l  -> guard (l <= v1)
+       Nothing -> return ()
+     case ub of
+       Just u  -> guard (v1 <= u)
+       Nothing -> return ()
+     return $ TopVar x v1 lb ub is
+
+slnNextVar :: Solutions -> Maybe Solutions
+slnNextVar Done = Nothing
+slnNextVar (TopVar x v _ _ is) = Just $ FixedVar x v $ startIter $ iLet x v is
+slnNextVar (FixedVar x v i)    = FixedVar x v `fmap` slnNextVar i
+
+
+
+
 -- Given a list of lower (resp. upper) bounds, compute the least (resp. largest)
 -- value that satisfies them all.
 iPickBounded :: BoundType -> [Bound] -> Maybe Integer
@@ -328,13 +446,23 @@
 iPickBounded bt bs =
   do xs <- mapM (normBound bt) bs
      return $ case bt of
-                Lower -> maximum xs + 1
-                Upper -> minimum xs - 1
+                Lower -> maximum xs
+                Upper -> minimum xs
   where
+  -- t < c*x
+  -- <=> t+1 <= c*x
+  -- <=> (t+1)/c <= x
+  -- <=> ceil((t+1)/c) <= x
+  -- <=> t `div` c + 1 <= x
   normBound Lower (Bound c t) = do k <- isConst t
-                                   return (div (k + c - 1) c)
+                                   return (k `div` c + 1)
+  -- c*x < t
+  -- <=> c*x <= t-1
+  -- <=> x   <= (t-1)/c
+  -- <=> x   <= floor((t-1)/c)
+  -- <=> x   <= (t-1) `div` c
   normBound Upper (Bound c t) = do k <- isConst t
-                                   return (div k c)
+                                   return (div (k-1) c)
 
 
 -- | The largest (resp. least) upper (resp. lower) bound on a term
@@ -377,6 +505,7 @@
 
 
 
+
 iModel :: Inerts -> [(Name,Integer)]
 iModel i = goBounds [] (bounds i)
   where
@@ -401,10 +530,13 @@
 --------------------------------------------------------------------------------
 -- Solving constraints
 
+solveIs0 :: Term -> S ()
+solveIs0 t = solveIs0' =<< apSubst t
+
 -- | Solve a constraint if the form @t = 0@.
 -- Assumes substitution has already been applied.
-solveIs0 :: Term -> S ()
-solveIs0 t
+solveIs0' :: Term -> S ()
+solveIs0' t
 
   -- A == 0
   | Just a <- isConst t = guard (a == 0)
@@ -415,12 +547,13 @@
       (q,0) -> addDef x (tConst q)
       _     -> mzero
 
-  -- x + S = 0
+  --  x + S = 0
+  -- -x + S = 0
   | Just (xc,x,s) <- tGetSimpleCoeff t =
     addDef x (if xc > 0 then tNeg s else s)
 
   -- A * S = 0
-  | Just (_, s) <- tFactor t  = addWork qZeroTerms s
+  | Just (_, s) <- tFactor t  = solveIs0 s
 
   -- See Section 3.1 of paper for details.
   -- We obtain an equivalent formulation but with smaller coefficients.
@@ -433,7 +566,7 @@
          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)
+         solveIs0 (negate (abs ak) |*| tVar v |+| tMapCoeff upd s)
 
   | otherwise = error "solveIs0: unreachable"
 
@@ -441,16 +574,20 @@
 modulus a m = a - m * div (2 * a + m) (2 * m)
 
 
+solveIsNeg :: Term -> S ()
+solveIsNeg t = solveIsNeg' =<< apSubst t
+
+
 -- | Solve a constraint of the form @t < 0@.
 -- Assumes that substitution has been applied
-solveIsNeg :: Term -> S ()
-solveIsNeg t
+solveIsNeg' :: Term -> S ()
+solveIsNeg' t
 
   -- A < 0
   | Just a <- isConst t = guard (a < 0)
 
   -- A * S < 0
-  |Just (_,s) <- tFactor t = addWork qNegTerms s
+  | Just (_,s) <- tFactor t = solveIsNeg s
 
   -- See Section 5.1 of the paper
   | Just (xc,x,s) <- tLeastVar t =
@@ -477,12 +614,14 @@
                  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)
+             solveIsNeg real
+             foldl orElse (solveIsNeg dark) (map solveIs0 gray)
              ) ctrs
 
   | otherwise = error "solveIsNeg: unreachable"
 
+orElse :: S () -> S () -> S ()
+orElse x y = mplus x y
 
 {- Note [Shadows]
 
@@ -522,6 +661,40 @@
 data Answer a = None | One a | Choice (Answer a) (Answer a)
                 deriving Show
 
+
+answerSize :: Answer a -> (Integer,Integer,Integer)
+answerSize = go 0 0 0
+  where
+  go !n !o !c ans =
+    case ans of
+      None  -> (n+1, o, c)
+      One _ -> (n, o + 1, c)
+      Choice x y ->
+        case go n o (c+1) x of
+          (n',o',c') -> go n' o' c' y
+
+
+dotAnswer :: (a -> Doc) -> Answer a -> Doc
+dotAnswer pp g0 = vcat [text "digraph {", nest 2 (fst $ go 0 g0), text "}"]
+  where
+  node x d            = integer x <+> brackets (text "label=" <> text (show d))
+                                                              <> semi
+  edge x y            = integer x <+> text "->" <+> integer y
+
+  go x None           = let x' = x + 1
+                        in seq x' ( node x "", x' )
+  go x (One a)        = let x' = x + 1
+                        in seq x' ( node x (show (pp a)), x' )
+  go x (Choice c1 c2) = let x'       = x + 1
+                            (ls1,x1) = go x' c1
+                            (ls2,x2) = go x1    c2
+                        in seq x'
+                           ( vcat [ node x "|"
+                                  , edge x x'
+                                  , edge x x1
+                                  , ls1
+                                  , ls2
+                                  ], x2 )
 toList :: Answer a -> [a]
 toList a = go a []
   where
@@ -586,20 +759,6 @@
                        , 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
@@ -620,19 +779,15 @@
 -- | 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 }
-        }
+addDef x t =
+  do newWork <- updS $ \rw -> let (newWork,newInerts) = iSolved x t (inerts rw)
+                              in (newWork, rw { inerts = newInerts })
+     mapM_ solveIsNeg newWork
 
-enqAndGo :: Field Term -> Term -> S ()
-enqAndGo q t =
+apSubst :: Term -> S Term
+apSubst t =
   do i <- get inerts
-     addWork q $ iApSubst i t
-     solveAll
+     return (iApSubst i t)
 
 
 
diff --git a/tests/qc.hs b/tests/qc.hs
new file mode 100644
--- /dev/null
+++ b/tests/qc.hs
@@ -0,0 +1,36 @@
+{-# LANGUAGE TemplateHaskell #-}
+import Data.Integer.SAT
+
+import Test.QuickCheck
+import System.Exit
+
+instance Arbitrary BoundType where
+  arbitrary = elements [Lower, Upper]
+
+withBounds :: Testable prop =>
+  BoundType -> [(Positive Integer, Integer)] -> (Integer -> prop) -> Property
+withBounds kind bs prop =
+  counterexample (show (map toBound bs)) $
+  case iPickBounded kind (map toBound bs) of
+    Nothing -> property Discard
+    Just n -> counterexample (show n) (property (prop n))
+  where
+    toBound (Positive c, t) = Bound c (tConst t)
+
+prop_lower, prop_upper :: [(Positive Integer, Integer)] -> Property
+prop_lower bs =
+  withBounds Lower bs $ \n ->
+    and [t <  c * n     | (Positive c, t) <- bs] &&
+    or  [t >= c * (n-1) | (Positive c, t) <- bs]
+prop_upper bs =
+  withBounds Upper bs $ \n ->
+    and [c * n     < t  | (Positive c, t) <- bs] &&
+    or  [c * (n+1) >= t | (Positive c, t) <- bs]
+
+-- This is so that the Template Haskell below can see the above properties.
+$(return [])
+
+main :: IO ()
+main = do ok <- $(quickCheckAll)
+          if ok then exitSuccess else exitFailure
+