diff --git a/Data/Presburger/Omega/Expr.hs b/Data/Presburger/Omega/Expr.hs
new file mode 100644
--- /dev/null
+++ b/Data/Presburger/Omega/Expr.hs
@@ -0,0 +1,958 @@
+
+-- | Expressions are the high-level interface for creating Presburger
+-- formulae.  As in Presburger arithmetic, expressions can
+-- represent addition, subtraction, quantification, inequalities, and boolean
+-- operators.
+--
+-- Expressions allow formulas to be input in a freeform manner.  When
+-- converted to a formula with 'expToFormula', they will be simplified to a
+-- form that the underlying library can use.
+-- Multplication is unrestricted; however, if an
+-- expression involves the product of two non-constant terms, it cannot be
+-- converted to a formula.
+-- 
+-- This module handles expressions and converts them to formulas.
+-- Sets and relations are managed by the "Data.Presburger.Omega.Set"
+-- and "Data.Presburger.Omega.Rel" modules.
+
+{-# OPTIONS_GHC -XBangPatterns
+                -XTypeFamilies
+                -XEmptyDataDecls
+                -XFlexibleInstances
+                -XFlexibleContexts
+                -XUndecidableInstances #-}
+module Data.Presburger.Omega.Expr
+    (-- * Expressions
+     Exp, IntExp, BoolExp,
+     Var,
+
+     -- ** Construction
+     nthVariable, takeFreeVariables, takeFreeVariables',
+     varE, nthVarE, intE, boolE, trueE, falseE, negateE,
+     sumE, prodE, notE, conjE, disjE,
+     (|&&|),
+     sumOfProductsE,
+     (|+|), (|-|), (|*|), (*|),
+     isZeroE, isNonnegativeE,
+     (|==|), (|/=|), (|>|), (|>=|), (|<|), (|<=|),
+     forallE, existsE,
+
+     -- ** Internal data structures
+     --
+     -- | These are exported to allow other modules to build the low-level
+     -- representation of expressions, and avoid the cost of simplifying
+     -- expressions.  Normally, the 'Exp' functions are sufficient.
+     Expr, IntExpr, BoolExpr,
+     PredOp(..),
+     wrapExpr, wrapSimplifiedExpr,
+     varExpr, sumOfProductsExpr, conjExpr, disjExpr, testExpr, existsExpr,
+
+     -- ** Operations on expressions
+     expEqual,
+     expToFormula,
+
+     -- ** Manipulating variables
+     rename,
+     adjustBindings,
+     variablesWithinRange,
+    )
+where
+
+import Control.Monad
+import Data.IORef
+import Data.List
+import Data.Maybe
+import qualified Data.IntMap as IntMap
+import Data.IntMap(IntMap)
+import qualified Data.Set as Set
+import Data.Set(Set)
+import Data.Unique
+import Debug.Trace
+import System.IO.Unsafe
+
+import Data.Presburger.Omega.LowLevel
+
+infixl 7 |*|
+infixl 7 *|
+infixl 6 |+|, |-|
+infix 4 |>|, |>=|, |<|, |<=|, |==|, |/=|
+infixr 3 |&&|
+
+-- | Integer and boolean-valued expressions.
+
+-- Expressions can convert themselves to a normalized form, under the hood.
+-- Normalization is only done when necessary.  The IORef allows an expression
+-- to be updated with its normalized form.
+newtype Exp t = Exp (IORef (ExprBox t))
+
+type IntExp = Exp Int
+type BoolExp = Exp Bool
+
+instance Show (Exp Int) where
+    showsPrec n e =
+        showsIntExprPrec emptyShowsEnv n (getSimplifiedExpr e)
+
+instance Show (Exp Bool) where
+    showsPrec n e =
+        showsBoolExprPrec emptyShowsEnv n (getSimplifiedExpr e)
+
+-- | We keep track of whether an expression is simplified.
+data ExprBox t =
+    ExprBox
+    { isSimplified :: {-# UNPACK #-} !Bool
+    , expression   :: !(Expr t)
+    }
+
+-- | Get an expression, without trying to simplify it.
+getExpr :: Exp t -> Expr t
+getExpr (Exp ref) = expression $ unsafePerformIO $ readIORef ref
+
+-- | Get the simplified form of an expression.
+getSimplifiedExpr :: Exp t -> Expr t
+getSimplifiedExpr (Exp ref) =
+    unsafePerformIO $ readIORef ref >>= simplifyAndUpdate
+    where
+      simplifyAndUpdate (ExprBox True  e) = return e
+      simplifyAndUpdate (ExprBox False e) =
+          let e'     = simplify e
+              newBox = ExprBox True e'
+          in do writeIORef ref $! newBox
+                return e'
+
+-- | Wrap an expression.
+wrapExpr :: Expr t -> Exp t
+wrapExpr e = Exp $ unsafePerformIO $ newIORef (ExprBox False e)
+
+-- | Wrap an expression that is known to be in simplified form.
+wrapSimplifiedExpr :: Expr t -> Exp t
+wrapSimplifiedExpr e = Exp $ unsafePerformIO $ newIORef (ExprBox True e)
+
+-- 'wrap' and 'get' are inverses
+{-# RULES
+"wrap/getExpr"		 forall x. getExpr (wrapExpr x) = x
+"wrapSimplified/getExpr" forall x. getExpr (wrapSimplifiedExpr x) = x
+"wrap/getSimplifiedExpr" forall x. getSimplifiedExpr (wrapSimplifiedExpr x) = x
+ #-}
+
+-- | Variables.  Variables are represented internally by de Bruijn indices.
+
+-- Variables are represented by a de Bruijn index.  The "innermost" variable
+-- is zero, and outer variables have higher indices.
+
+-- The 'Quantified' constructor is used temporarily when building a quantified
+-- expression.  It is only seen by 'rename' and 'adjustBindings'.
+data Var = Bound {-# UNPACK #-} !Int
+         | Quantified !Unique
+           deriving(Eq, Ord)
+
+-- | Produce the Nth bound variable.  Zero is the innermost variable index.
+nthVariable :: Int -> Var
+nthVariable = Bound
+
+-- | Construct a new quantified variable.
+newQuantified :: IO Var
+newQuantified = do u <- newUnique
+                   return (Quantified u)
+
+freeVariables :: [Var]
+freeVariables = map Bound [0..]
+
+-- | Produce a set of variables to use as free variables in an expression.
+-- This produces the list @[nthVariable 0, nthVariable 1, ...]@
+takeFreeVariables :: Int -> [Var]
+takeFreeVariables n = take n freeVariables
+
+-- | Like 'takeFreeVariables', but produce the expression corresponding to
+-- each variable.
+takeFreeVariables' :: Int -> [IntExp]
+takeFreeVariables' n = map varE $ take n freeVariables
+
+-------------------------------------------------------------------------------
+-- Building expressions
+
+varE :: Var -> IntExp
+varE v = wrapExpr $ VarE v
+
+nthVarE :: Int -> IntExp
+nthVarE n = varE (nthVariable n)
+
+intE :: Int -> IntExp
+intE n = wrapExpr $ LitE n
+
+boolE :: Bool -> BoolExp
+boolE b = wrapExpr $ LitE b
+
+trueE, falseE :: BoolExp
+trueE = boolE True
+falseE = boolE False
+
+-- | Multiplication by -1
+negateE :: IntExp -> IntExp
+negateE e = wrapExpr $ CAUE Prod (-1) [getExpr e]
+
+-- | Summation
+sumE :: [IntExp] -> IntExp
+sumE es = wrapExpr $ CAUE Sum 0 $ map getExpr es
+
+-- | Multiplication
+prodE :: [IntExp] -> IntExp
+prodE es = wrapExpr $ CAUE Prod 1 $ map getExpr es
+
+-- | Logical negation
+notE :: BoolExp -> BoolExp
+notE e = wrapExpr $ NotE $ getExpr e
+
+-- | Conjunction
+conjE :: [BoolExp] -> BoolExp
+conjE es = wrapExpr $ CAUE Conj True $ map getExpr es
+
+-- | Disjunction
+disjE :: [BoolExp] -> BoolExp
+disjE es = wrapExpr $ CAUE Disj False $ map getExpr es
+
+-- | Conjunction
+(|&&|) :: BoolExp -> BoolExp -> BoolExp
+e |&&| f = wrapExpr $ CAUE Conj True [getExpr e, getExpr f]
+
+-- | Add
+(|+|) :: IntExp -> IntExp -> IntExp
+e |+| f = sumE [e, f]
+
+-- | Subtract
+(|-|) :: IntExp -> IntExp -> IntExp
+e |-| f = sumE [e, negateE f]
+
+-- | Multiply
+(|*|) :: IntExp -> IntExp -> IntExp
+e |*| f = prodE [e, f]
+
+-- | Multiply by an integer
+(*|) :: Int -> IntExp -> IntExp
+n *| f = wrapExpr $ CAUE Prod n [getExpr f]
+
+-- | Test whether an integer expression is zero
+isZeroE :: IntExp -> BoolExp
+isZeroE e = wrapExpr $ PredE IsZero $ getExpr e
+
+-- | Test whether an integer expression is nonnegative
+isNonnegativeE :: IntExp -> BoolExp
+isNonnegativeE e = wrapExpr $ PredE IsGEZ $ getExpr e
+
+-- | Equality test
+(|==|) :: IntExp -> IntExp -> BoolExp
+e |==| f = isZeroE (e |-| f)
+
+-- | Inequality test
+(|/=|) :: IntExp -> IntExp -> BoolExp
+e |/=| f = disjE [e |>| f, e |<| f]
+
+-- | Greater than
+(|>|) :: IntExp -> IntExp -> BoolExp
+e |>| f = isNonnegativeE (wrapExpr $ CAUE Sum (-1) [getExpr e, getExpr $ negateE f])
+
+-- | Less than
+(|<|) :: IntExp -> IntExp -> BoolExp
+e |<| f = f |>| e
+
+-- | Greater than or equal
+(|>=|) :: IntExp -> IntExp -> BoolExp
+e |>=| f = isNonnegativeE (e |-| f)
+
+-- | Less than or equal
+(|<=|) :: IntExp -> IntExp -> BoolExp
+e |<=| f = f |>=| e
+
+-- | Build a universally quantified formula.
+forallE :: (Var -> Exp t) -> Exp t
+forallE f = wrapExpr $ QuantE Forall $ getExpr $ withFreshVariable f
+
+-- | Build an existentially quantified formula.
+existsE :: (Var -> Exp t) -> Exp t
+existsE f = wrapExpr $ QuantE Exists $ getExpr $ withFreshVariable f
+
+-- | Use a fresh variable in an expression.  After the expression is
+-- constructed, rename/adjust variable indices so that the fresh variable
+-- has index 0 and all other free variables' indices are incremented
+-- by 1.
+withFreshVariable :: (Var -> Exp t) -> Exp t
+withFreshVariable f =unsafePerformIO $ do
+  v <- newQuantified
+  return $ rename v (Bound 0) $ adjustBindings 0 1 $ f v
+
+-------------------------------------------------------------------------------
+
+-- | The internal representation of expressions.
+data Expr t where
+    -- Application of a commutative and associative operator
+    CAUE :: !(CAUOp t)          -- operator
+         -> !t                  -- literal operand
+         -> [Expr t]             -- other operands
+         -> Expr t
+
+    -- A predicate on an integer expression
+    PredE :: !PredOp            -- operator
+          -> Expr Int            -- integer operand
+          -> Expr Bool
+
+    -- Boolean negation
+    NotE :: Expr Bool -> Expr Bool
+
+    -- A literal
+    LitE :: !t -> Expr t
+
+    -- A variable.  Only integer-valued variables are permitted.
+    VarE :: !Var -> Expr Int
+
+    -- An expression quantified over an integer variable
+    QuantE :: !Quantifier -> Expr t -> Expr t
+
+type IntExpr = Expr Int
+type BoolExpr = Expr Bool
+
+-- | A commutative and associative operator with a unit.
+-- The type parameter 't' gives the operator's parameter and return type.
+data CAUOp t where
+    Sum  :: CAUOp Int
+    Prod :: CAUOp Int
+    Conj :: CAUOp Bool 
+    Disj :: CAUOp Bool
+
+instance Eq (CAUOp t) where
+    Sum  == Sum  = True
+    Prod == Prod = True
+    Conj == Conj = True
+    Disj == Disj = True
+    _    == _    = False
+
+instance Show (CAUOp t) where
+    show Sum  = "Sum"
+    show Prod = "Prod"
+    show Conj = "Conj"
+    show Disj = "Disj"
+
+-- | A predicate on an integer expresion.
+data PredOp = IsZero | IsGEZ
+              deriving(Eq, Show)
+
+-- Quantifiers.
+data Quantifier = Forall | Exists
+                  deriving(Eq, Show)
+
+varExpr :: Var -> IntExpr
+varExpr = VarE
+
+-- | Create a sum of products expression
+sumOfProductsE :: Int           -- ^ constant part of sum
+               -> [(Int, [Var])] -- ^ product terms
+               -> IntExp
+sumOfProductsE n prods = wrapSimplifiedExpr $ CAUE Sum n $ map prod prods
+    where
+      prod (n, vars) = CAUE Prod n $ map VarE vars
+
+sumOfProductsExpr :: Int           -- ^ constant part of sum
+                  -> [(Int, [Var])] -- ^ product terms
+                  -> IntExpr
+sumOfProductsExpr n prods = CAUE Sum n $ map prod prods
+    where
+      prod (n, vars) = CAUE Prod n $ map VarE vars
+
+testExpr :: PredOp -> IntExpr -> BoolExpr
+testExpr p e = PredE p e
+
+conjExpr :: [BoolExpr] -> BoolExpr
+conjExpr = CAUE Conj True
+
+disjExpr :: [BoolExpr] -> BoolExpr
+disjExpr = CAUE Disj False
+
+existsExpr :: BoolExpr -> BoolExpr
+existsExpr e = QuantE Exists e
+
+-------------------------------------------------------------------------------
+
+isLitE :: Expr t -> Bool
+isLitE (LitE _) = True
+isLitE _        = False
+
+deconstructProduct :: IntExpr -> Term Int
+deconstructProduct (CAUE Prod n xs) = (n, xs)
+deconstructProduct e                = (unit Prod, [e])
+
+rebuildProduct :: Term Int -> Expr Int
+rebuildProduct (1, [e]) = e
+rebuildProduct (n, es)  = CAUE Prod n es
+
+deconstructSum :: Expr Int -> Term Int
+deconstructSum (CAUE Sum n xs) = (n, xs)
+deconstructSum e               = (unit Sum, [e])
+
+rebuildSum :: Term Int -> Expr Int
+rebuildSum (1, [e]) = e
+rebuildSum (n, es)  = CAUE Sum n es
+
+-- Get the 'equality' operator for type t.
+cauEq :: CAUOp t -> t -> t -> Bool
+cauEq Sum  = (==)
+cauEq Prod = (==)
+cauEq Conj = (==)
+cauEq Disj = (==)
+
+-- Get the 'shows' operator for type t.
+cauShows :: CAUOp t -> t -> ShowS
+cauShows Sum  = shows
+cauShows Prod = shows
+cauShows Conj = shows
+cauShows Disj = shows
+
+-- Get the zero for a CAU op (if one exists)
+zero :: CAUOp t -> Maybe t
+zero Sum  = Nothing
+zero Prod = Just 0
+zero Conj = Just False
+zero Disj = Just True
+
+-- Get the unit for a CAU op
+unit :: CAUOp t -> t
+unit Sum  = 0
+unit Prod = 1
+unit Conj = True
+unit Disj = False
+
+-- | True if the literal is the operator's zero.
+isZeroOf :: t -> CAUOp t -> Bool
+l `isZeroOf` op = case zero op
+                  of Nothing -> False
+                     Just z  -> cauEq op l z
+
+-- | True if the literal is the operator's unit.
+isUnitOf :: t -> CAUOp t -> Bool
+l `isUnitOf` op = cauEq op (unit op) l
+
+-- Evaluate an operator on a list of literals
+evalCAUOp :: CAUOp t -> [t] -> t
+evalCAUOp Sum  = sum
+evalCAUOp Prod = product
+evalCAUOp Conj = and
+evalCAUOp Disj = or
+
+-- Evaluate a predicate
+evalPred :: PredOp -> Int -> Bool
+evalPred IsZero = (0 ==)
+evalPred IsGEZ  = (0 <=)
+
+-------------------------------------------------------------------------------
+-- Showing expressions
+
+appPrec = 10
+mulPrec = 7
+addPrec = 6
+relPrec = 4
+lamPrec = 0
+
+-- An environment for showing expressions.
+--
+-- Quantified variables are shown as lambda-bound variables.  This structure
+-- keeps track of lambda-bound variable names and how to show them.
+
+data ShowsEnv =
+    ShowsEnv
+    { -- How to show the n_th bound variable, given a precedence context
+      showNthVar :: [Int -> ShowS]
+      -- Number of bound variables we know about.
+      --   numBound e == length (showNthVar e)
+    , numBound   :: !Int
+      -- Names for new bound variables
+    , varNames   :: [ShowS]
+    }
+
+emptyShowsEnv =
+    ShowsEnv
+    { showNthVar = []
+    , numBound = 0
+    , varNames = map showChar $
+                 ['x', 'y', 'z'] ++
+                 ['a' .. 'w'] ++
+                 [error "out of variable names"]
+    }
+
+-- Add a variable binding to the environment
+bindVariable :: ShowsEnv -> (ShowS, ShowsEnv)
+bindVariable env =
+    case varNames env
+    of nm : nms ->
+           let env' = ShowsEnv
+                      { showNthVar = showVar nm : showNthVar env
+                      , numBound   = 1 + numBound env
+                      , varNames   = nms
+                      }
+           in (nm, env')
+    where
+      -- Showing a variable produces "varE varName"
+      showVar nm n = showParen (n >= appPrec) $ showString "varE " . nm
+
+showsVarPrec :: ShowsEnv -> Int -> Var -> ShowS
+showsVarPrec env prec (Bound i) =
+    if i < numBound env
+    then (showNthVar env !! i) prec
+    else shift (numBound env)
+    where
+      -- The variable is not bound locally, so show its constructor.
+      -- We have to subtract an offset to account for the local variable
+      -- bindings, basically undoing the shift that 'withFreshVariable'
+      -- applies.
+      shift n = showParen (prec >= appPrec) $
+                    showString "nthVarE " . shows (i-n)
+
+-- Unique is not showable, but users shouldn't see quantified variables anyway
+showsVarPrec _ _ (Quantified u) = showString "(Quantified _)"
+
+showsInt :: Int -> ShowS
+showsInt n | n >= 0    = showString "intE " . shows n
+           | otherwise = showString "intE " . showParen True (shows n)
+
+
+showsIntExprPrec :: ShowsEnv -> Int -> IntExpr -> ShowS
+showsIntExprPrec env n expression =
+    case expression
+    of CAUE Sum lit es  -> showParen (n >= addPrec) $ showSum env lit es
+       CAUE Prod lit es -> showParen (n >= mulPrec) $ showProd env lit es
+       LitE l           -> showParen (n >= appPrec) $
+                           showsInt l
+       VarE v           -> showsVarPrec env n v
+       QuantE q e       -> showParen (n >= appPrec) $
+                           showQuantifier showsIntExprPrec env q e
+
+showsBoolExprPrec :: ShowsEnv -> Int -> BoolExpr -> ShowS
+showsBoolExprPrec env n expression =
+    case expression
+    of CAUE Conj lit es
+           | lit        -> let texts = map (showsBoolExprPrec env 0) es
+                           in texts `showSepBy` showString " |&&| "
+           | otherwise  -> showString "falseE"
+       CAUE Disj lit es
+           | lit        -> showString "trueE"
+           | otherwise  -> let texts = map (showsBoolExprPrec env 0) es
+                           in showParen (n >= appPrec) $
+                              showString "disjE " . showsList texts
+       PredE p e        -> let operator =
+                                   case p
+                                   of IsZero -> showString "isZeroE "
+                                      IsGEZ  -> showString "isNonnegativeE "
+                           in showParen (n >= appPrec) $
+                              operator . showsIntExprPrec env appPrec e
+       NotE e           -> showString "notE " . showsBoolExprPrec env appPrec e
+       LitE True        -> showString "trueE"
+       LitE False       -> showString "falseE"
+       QuantE q e       -> showParen (n >= appPrec) $
+                           showQuantifier showsBoolExprPrec env q e
+
+-- Show a sum term
+showSum env lit es =
+    -- The first element of the summation gets shown a little differently.
+    -- There are a couple of cases, depending on what is the first element.
+    if lit == 0
+    then case es
+         of e : es' -> showsIntExprPrec env addPrec e . showSumTail es'
+            []      -> showsInt 0
+    else showsInt lit . showSumTail es
+    where
+      -- Show the tail of a sum term.  Each expression is preceded by
+      -- the |+| or |-| operator.
+      showSumTail es = foldr (.) id $ map showSumTailElement es
+
+      showSumTailElement e =
+          case deconstructProduct e
+          of (1, es)             -> add . showProd env 1 es
+             (-1, es)            -> sub . showProd env 1 es
+             (n, es) | n >= 0    -> add . showProd env n es
+                     | otherwise -> sub . showProd env (negate n) es
+
+      add = showString " |+| "
+      sub = showString " |-| "
+
+-- Show a product term
+showProd env lit es =
+    let text = map (showsIntExprPrec env mulPrec) es
+        textLit = if lit == 1
+                  then id
+                  else showsPrec mulPrec lit . showString " *| "
+    in textLit . (text `showSepBy` showString " |*| ")
+        where
+      showMulOperator = showString " |*| "
+
+-- Show a list in [,,] syntax
+showsList :: [ShowS] -> ShowS
+showsList ss z =
+    showChar '[' $
+    foldr ($) (showChar ']' $ z) (intersperse (showString ", ") ss)
+
+-- Show a list with a separator interspersed
+showSepBy :: [ShowS] -> ShowS -> ShowS
+xs `showSepBy` sep = foldr (.) id (intersperse sep xs)
+
+-- Show a quantified expression, e.g. (forallE. (x + 1))
+showQuantifier :: (ShowsEnv -> Int -> Expr t -> ShowS)
+               -> ShowsEnv -> Quantifier -> Expr t -> ShowS
+showQuantifier showExpr env q e =
+    let quantifier = case q
+                     of Forall -> showString "forallE $ \\"
+                        Exists -> showString "existsE $ \\"
+
+        -- Take a new variable name
+        (varName, env') = bindVariable env
+
+    in quantifier . varName . showString " -> " . showExpr env' lamPrec e
+
+-------------------------------------------------------------------------------
+-- Syntactic equality on expressions
+
+-- | Decide whether two expressions are syntactically equal, modulo
+-- commutativity, associativity, and alpha-renaming.
+expEqual :: Eq t => Expr t -> Expr t -> Bool
+expEqual expr1 expr2 =
+    case (expr1, expr2)
+    of (CAUE op1 l1 es1, CAUE op2 l2 es2) ->
+          op1 == op2 && l1 == l2 && expListsEqual es1 es2
+
+       (PredE op1 e1, PredE op2 e2) ->
+          op1 == op2 && expEqual e1 e2
+
+       (NotE e1, NotE e2) -> expEqual e1 e2
+
+       (LitE l1, LitE l2) -> l1 == l2
+
+       (VarE v1, VarE v2) -> v1 == v2
+
+       (QuantE q1 e1, QuantE q2 e2) ->
+          q1 == q2 && expEqual e1 e2
+
+       (_, _) -> False          -- Different constructors
+
+-- Decide whether two unordered expression lists are equal.
+-- For each element of the first list, find
+-- a matching element of the second list and repeat.
+expListsEqual :: Eq t => [Expr t] -> [Expr t] -> Bool
+expListsEqual (e:es1) es2 =
+    case findEqualExpr e es2
+    of Just (_, es2') -> expListsEqual es1 es2'
+       Nothing        -> False
+
+expListsEqual [] [] = True      -- All elements matched
+expListsEqual [] _  = False     -- Some leftover elements in es2
+
+-- Find an equal expression in the list.
+findEqualExpr :: Eq t => Expr t -> [Expr t] -> Maybe (Expr t, [Expr t])
+findEqualExpr searchE es = go es id
+    where
+      go (e:es) h | expEqual searchE e = Just (e, h es)
+                  | otherwise          = go es (h . (e:))
+      go []     _                      = Nothing
+
+-------------------------------------------------------------------------------
+-- Simplification rules
+
+-- This is the main rule for simplifying an expression.
+--
+-- First, subexpressions are simplified (simplifyRec).
+-- Then "basic" simplifications are performed.  These restructure the
+-- current term, but no other terms.
+-- Then complex simplifications are performed that restructure the current
+-- term and subtems.
+
+-- | Normalize an expression.
+simplify :: Expr t -> Expr t
+simplify e =
+    complexSimplifications $ basicSimplifications $ simplifyRec e
+
+simplifyRec :: Expr t -> Expr t
+simplifyRec expr =
+    case expr
+    of CAUE op lit es -> CAUE op lit $ map simplify es
+       PredE op e1 -> PredE op $ simplify e1
+       NotE e -> NotE $ simplify e
+       LitE _ -> expr
+       VarE v -> expr
+       QuantE q e -> QuantE q $ simplify e 
+
+basicSimplifications :: Expr t -> Expr t
+basicSimplifications = zus . peval . flatten
+
+-- Some complex simplifications require steps of simplification to be re-run.
+complexSimplifications :: Expr t -> Expr t
+complexSimplifications e =
+    case e
+    of CAUE Sum _ _  -> basicSimplifications $ collect e
+       CAUE Prod _ _ -> posToSop e
+       _             -> e
+
+-- Convert a product of sums to a sum of products.  If conversion happens,
+-- simplification is re-run.
+
+posToSop :: Expr Int -> Expr Int
+posToSop expr@(CAUE Prod n es)
+    | all (isSingletonList . snd) terms =
+        -- If no terms are sums, then the expression is unchanged
+        expr
+
+    | otherwise =
+          let -- Make a list of lists.
+              -- The expression corresponds to
+              --   product (map sum terms')
+              terms' = [LitE n] : map mkTermList terms
+
+              -- The cartesian product converts this to a sum of products.
+              sop    = sequence terms'
+              expr'  = CAUE Sum 0 (map (CAUE Prod 1) sop)
+          in simplify expr'
+    where
+      terms = map deconstructSum es
+      mkTermList (n, es) = LitE n : es
+      isSingletonList [_] = True
+      isSingletonList _   = False
+
+posToSop expr = expr            -- Terms other than products are not modified
+
+-- Flatten nested CA expressions
+flatten :: forall t. Expr t -> Expr t
+flatten (CAUE op lit es) = CAUE op lit (flat es)
+    where
+      -- Wherever a nested CA expression with the same operator appears,
+      -- include its terms in the list
+      flat :: [Expr t] -> [Expr t]
+      flat (e:es) = case e
+                    of CAUE op2 lit2 es2
+                           | op == op2 -> LitE lit2 : es2 ++ flat es
+                       _ -> e:flat es
+      flat []     = []
+flatten e = e
+
+-- Partially evaluate an expression
+peval :: Expr t -> Expr t
+peval exp@(CAUE op l es) =
+    case partition isLitE es
+    of ([], _)         -> exp
+       (lits, notLits) -> let literals = l : map fromLitE lits
+                          in CAUE op (evalCAUOp op literals) notLits
+    where
+      fromLitE (LitE l) = l
+      fromLitE _        = error "peval: unexpected expression"
+
+peval exp@(PredE op e) =
+    case e
+    of LitE l -> LitE $ evalPred op l
+       _      -> exp
+
+peval exp@(NotE e) =
+    case e
+    of LitE l -> LitE $ not l
+       _      -> exp
+
+peval e = e
+
+-- Zero, unit, singleton rules.  May eliminate an
+-- expression here.
+zus :: Expr t -> Expr t
+zus exp@(CAUE op l es) =
+    case es
+    of [] -> LitE l
+       [e] | l `isZeroOf` op -> LitE l -- zero * x = zero
+           | l `isUnitOf` op -> e      -- unit * x = x
+           | otherwise       -> exp    -- no simplificaiton
+       _ | l `isZeroOf` op   -> LitE l -- zero * x = zero
+         | otherwise         -> exp    -- no simplification
+
+zus e = e
+
+-- Given a sum of products, collect terms that differ only in their
+-- constant multiplier.
+--
+-- For example:
+--
+--  collect (2xy + 3x - 3xy)
+--  becomes (-1)xy + 3x
+
+type Term t = (t, [Expr t])
+
+collect :: Expr Int -> Expr Int
+collect (CAUE Sum literal es) =
+    let es' = map simplify $
+              map rebuildProduct $
+              collectTerms $
+              map deconstructProduct es
+    in CAUE Sum literal es'
+
+    where
+      collectTerms :: [Term Int] -> [Term Int]
+      collectTerms (t:ts) =
+          case collectTerm t ts of (t', ts') -> t':collectTerms ts'
+      collectTerms [] = []
+
+      -- Collect together all terms from the list that differ from
+      -- the first term only in their multiplier.  The collected terms'
+      -- multipliers are summed.  The result is the collected term
+      -- and the unused terms from the list.
+      collectTerm :: Term Int -> [Term Int] -> (Term Int, [Term Int])
+      collectTerm (factor, t) terms =
+          let (equalTerms, terms') = partition (sameTerms t) terms
+              factor'              = factor + sum (map fst equalTerms)
+          in ((factor', t), terms')
+
+      -- Decide whether the expression lists are equal.
+      sameTerms t (_, t') = expListsEqual t t'
+
+collect e = e                   -- Terms other than sums do not change
+
+-------------------------------------------------------------------------------
+-- Converting expressions to formulas
+
+-- | Look up a variable in a list.  The variable's position is its
+-- de Bruijn index.
+
+lookupVar :: Int -> [VarHandle] -> VarHandle
+lookupVar n (v : vars) | n > 0  = lookupVar (n - 1) vars
+                       | n == 0 = v
+                       | otherwise = error "lookupVar: negative index"
+
+lookupVar _ [] = error "lookupVar: variable index out of range"
+
+-- | Convert a boolean expression to a formula.
+--
+-- The expression must be a Presburger formula.  In particular, if an
+-- expression involves the product of two non-constant terms, it cannot be
+-- converted to a formula.  The library
+-- internally simplifies expressions to sum-of-products form, so complex
+-- expressions are valid as long as each simplified product has at most
+-- one variable.
+-- The library currently cannot create a set or relation if any
+-- integer expressions contain quantifiers, but this restriction could be
+-- lifted in the future.
+
+expToFormula :: [VarHandle]     -- ^ Free variables
+             -> BoolExp         -- ^ Expression to convert
+             -> Formula
+expToFormula freeVars e = exprToFormula freeVars (getSimplifiedExpr e)
+
+exprToFormula :: [VarHandle]     -- ^ Free variables
+              -> BoolExpr        -- ^ Expression to convert
+              -> Formula
+exprToFormula freeVars expr =
+    case expr
+    of CAUE op lit es
+           | lit `isUnitOf` op ->
+               case op
+               of Conj -> conjunction $ map (exprToFormula freeVars) es
+                  Disj -> disjunction $ map (exprToFormula freeVars) es
+                  _    -> expToFormulaError "unhandled operator"
+           | otherwise ->
+               -- This boolean literal overrides all other terms
+               if lit then true else false
+
+       PredE op e ->
+           case sumToConstraint freeVars e
+           of (terms, constant) ->
+                  case op
+                  of IsZero -> equality terms constant
+                     IsGEZ  -> inequality terms constant
+
+       NotE e -> negation $ exprToFormula freeVars e
+
+       LitE True  -> true
+       LitE False -> false
+
+       QuantE q e -> let body v = exprToFormula (v:freeVars) e
+                     in case q
+                        of Forall -> qForall body
+                           Exists -> qExists body
+
+-- | Convert an integer term to a coefficients for an equality or
+-- inequality constraint.
+sumToConstraint :: [VarHandle]  -- ^ free variables
+                -> IntExpr      -- ^ expression to convert
+                -> ([Coefficient], Int)
+sumToConstraint freeVars expr =
+    case deconstructSum expr
+    of (constant, terms) -> (map deconstructTerm terms, constant)
+    where
+      deconstructTerm :: IntExpr -> Coefficient
+      deconstructTerm expr =
+          case deconstructProduct expr
+          of (n, [VarE (Bound i)]) -> Coefficient (lookupVar i freeVars) n
+             _ -> expToFormulaError "expression is non-affine"
+
+expToFormulaError :: String -> a
+expToFormulaError s = error $ "expToFormula: " ++ s
+
+-- | Substitute a single variable in an expression.
+rename  :: Var               -- ^ variable to replace
+        -> Var               -- ^ its replacement
+        -> Exp t             -- ^ expression to rename
+        -> Exp t             -- ^ renamed expression
+rename v1 v2 e = wrapExpr $ renameExpr v1 v2 $ getExpr e
+
+renameExpr :: Var               -- ^ variable to replace
+           -> Var               -- ^ its replacement
+           -> Expr t            -- ^ expression to rename
+           -> Expr t            -- ^ renamed expression
+renameExpr !v1 v2 expr = rn expr
+    where
+      rn :: forall t. Expr t -> Expr t
+      rn (CAUE op lit es) = CAUE op lit $ map rn es
+      rn (PredE op e)     = PredE op $ rn e
+      rn (NotE e)         = NotE $ rn e
+      rn expr@(LitE _)    = expr
+      rn expr@(VarE v)    | v == v1   = VarE v2
+                          | otherwise = expr
+      rn (QuantE q e)     = QuantE q $ renameExpr (bump v1) (bump v2) e
+
+      -- Increment a de Bruijn index
+      bump (Bound n)        = Bound (n+1)
+      bump v@(Quantified _) = v
+
+-- | Adjust bound variable bindings by adding an offset to all bound variable
+-- indices beyond a given level.
+adjustBindings :: Int           -- ^ first variable to change
+               -> Int           -- ^ Amount to shift by
+               -> Exp t         -- ^ Input expression
+               -> Exp t         -- ^ Adjusted expression
+adjustBindings firstBound shift e =
+    wrapExpr $ adjustBindingsExpr firstBound shift $ getExpr e
+
+adjustBindingsExpr :: Int       -- ^ first variable to change
+                   -> Int       -- ^ Amount to shift by
+                   -> Expr t    -- ^ Input expression
+                   -> Expr t    -- ^ Adjusted expression
+adjustBindingsExpr !firstBound !shift e = adj e
+    where
+      adj :: Expr t -> Expr t
+      adj (CAUE op lit es) = CAUE op lit (map adj es)
+      adj (PredE op e)     = PredE op (adj e)
+      adj (NotE e)         = NotE (adj e)
+      adj expr@(LitE _)    = expr
+      adj expr@(VarE v)    = case v
+                             of Bound n
+                                    | n >= firstBound ->
+                                        VarE $ Bound (n + shift)
+                                    | otherwise ->
+                                        expr
+                                Quantified _ -> expr
+      adj (QuantE q e)     = QuantE q $
+                             adjustBindingsExpr (firstBound + 1) shift e
+
+-- | True if the expression has no more than the specified number
+-- of free variables.
+variablesWithinRange :: Int -> Exp t -> Bool
+variablesWithinRange n e = check n $ getExpr e
+    where
+      check :: Int -> Expr t -> Bool
+      check n e = check' e
+          where
+            check' :: Expr t -> Bool
+            check' (CAUE _ _ es)         = all check' es
+            check' (PredE _ e)           = check' e
+            check' (NotE e)              = check' e
+            check' (LitE _)              = True
+            check' (VarE (Bound i))      = i < n
+            check' (VarE (Quantified _)) = quantifiedVar
+            check' (QuantE _ e)          = check (n+1) e
+
+      quantifiedVar = error "Unexpected quantified variable"
diff --git a/Data/Presburger/Omega/LowLevel.hsc b/Data/Presburger/Omega/LowLevel.hsc
new file mode 100644
--- /dev/null
+++ b/Data/Presburger/Omega/LowLevel.hsc
@@ -0,0 +1,990 @@
+
+-- | This module provides a low-level interface for creating,
+-- manipulating, and querying Presburger arithmetic formulae.
+-- The real work is done by the C++ Omega library
+-- (<http://github.com/davewathaverford/the-omega-project>).
+--
+-- The main data types are 'OmegaSet' and 'OmegaRel', which use a formula
+-- to define a set or relation, respectively, on integer-valued points in
+-- Cartesian space.
+-- A typical use involves creating a Presburger arithmetic 'Formula', using
+-- it to create a set or relation, and then querying the set or relation.
+ 
+{-# OPTIONS_GHC -XForeignFunctionInterface -fwarn-incomplete-patterns
+                -XEmptyDataDecls -XRankNTypes -XMultiParamTypeClasses
+                -XFunctionalDependencies -XTypeSynonymInstances
+                -XFlexibleInstances -XFlexibleContexts #-}
+
+module Data.Presburger.Omega.LowLevel
+    (-- * Sets and relations
+     Presburger,
+     OmegaSet, newOmegaSet,
+     OmegaRel, newOmegaRel,
+
+     -- * Inspecting sets and relations directly
+     queryDNFSet, queryDNFRelation,
+
+     -- * Queries on sets and relations
+     lowerBoundSatisfiable, upperBoundSatisfiable,
+     obviousTautology, definiteTautology,
+     exact, inexact, unknown,
+
+     -- * Creating new sets and relations from old ones
+
+     -- ** Bounds
+     upperBound, lowerBound,
+
+     -- ** Binary operations
+     equal, union, intersection, composition,
+     restrictDomain, restrictRange,
+     difference, crossProduct, 
+     Effort(..),
+     gist,
+
+     -- ** Unary operations
+     transitiveClosure,
+     domain, range, inverse, complement,
+     deltas, approximate,
+
+     -- * Constructing formulas
+     Formula,
+     true, false,
+     conjunction, disjunction, negation,
+     VarHandle,
+     qForall, qExists,
+     Coefficient(..),
+     inequality, equality
+     )
+where
+
+#include "C_omega.h"
+
+#let alignof x = "%d", __alignof__(x)
+
+import Control.Monad
+import Data.Int
+import Data.List(findIndex)
+import Data.Word
+import Foreign.C
+import Foreign.ForeignPtr
+import qualified Foreign.Marshal.Alloc as ForeignAlloc
+import Foreign.Marshal.Array
+import Foreign.Ptr
+import Foreign.Storable
+import System.IO.Unsafe(unsafePerformIO)
+
+-------------------------------------------------------------------------------
+-- Data types, classes, and functions imported from C++
+
+-- External data types, these have the same name as in C.
+data Relation                   -- A set or relation
+data Form                       -- A logic formula (Formula)
+data F_And                      -- A conjunction
+data F_Declaration              -- A forall or exists formula
+data Var_Decl                   -- A handle to a variable
+data DNF_Iterator               -- Iterator over a DNF clause
+data Conjunct                   -- One conjunct within a DNF clause
+data Tuple_Iterator a           -- Iterator over a Tuple (in Omega library)
+data EQ_Iterator                -- Iterator over a set of EQ constraints
+data EQ_Handle                  -- Handle to an EQ constraint
+data GEQ_Iterator               -- Iterator over a set of GEQ constraints
+data GEQ_Handle                 -- Handle to a GEQ constraint
+data Constr_Vars_Iter           -- Iterate over coefficients in a constraint
+
+class Constraint a              -- The 'Constraint' base class
+
+instance Constraint EQ_Handle
+instance Constraint GEQ_Handle
+
+-- Pointers to external data types
+type C_Relation         = Ptr Relation
+type C_Form             = Ptr Form
+type C_And              = Ptr F_And
+type C_Quantifier       = Ptr F_Declaration
+type C_Var              = Ptr Var_Decl
+type C_DNF_Iterator     = Ptr DNF_Iterator
+type C_Conjunct         = Ptr Conjunct
+type C_Tuple_Iterator a = Ptr (Tuple_Iterator a)
+type C_EQ_Iterator      = Ptr EQ_Iterator
+type C_EQ_Handle        = Ptr EQ_Handle
+type C_GEQ_Iterator     = Ptr GEQ_Iterator
+type C_GEQ_Handle       = Ptr GEQ_Handle
+type C_Constr_Vars_Iter = Ptr Constr_Vars_Iter
+
+-- | The 'gist' routine takes a parameter specifying how much effort to
+-- put into generating a good result.  Higher effort takes more time.
+-- It's unspecified what a given effort level does.
+data Effort = Light | Moderate | Strenuous
+              deriving (Eq, Show, Enum)
+
+-- Everything containing a formula is an instance of class Logical
+class Logical f where
+    add_and    :: f -> IO C_Form
+    add_or     :: f -> IO C_Form
+    add_not    :: f -> IO C_Form
+    add_forall :: f -> IO C_Quantifier
+    add_exists :: f -> IO C_Quantifier
+    convert_to_and :: f -> IO C_And
+    finalize   :: f -> IO ()
+
+instance Logical C_Relation where
+    add_and    = hsw_relation_add_and
+    add_or     = hsw_relation_add_or
+    add_not    = hsw_relation_add_not
+    add_forall = hsw_relation_add_forall
+    add_exists = hsw_relation_add_exists
+    -- We take advantage of the fact that C_And is a subclass of C_Form
+    -- here, and simply cast the pointer.
+    convert_to_and r = liftM castPtr $ hsw_relation_add_and r
+    finalize   = hsw_relation_finalize
+
+instance Logical C_Form where
+    add_and    = hsw_formula_add_and
+    add_or     = hsw_formula_add_or
+    add_not    = hsw_formula_add_not
+    add_forall = hsw_formula_add_forall
+    add_exists = hsw_formula_add_exists
+    convert_to_and = hsw_formula_to_and
+    finalize   = hsw_formula_finalize
+
+-- C_And is a subclass of C_Form and implements all its methods.
+-- Consequently, we simply cast to C_Form
+instance Logical C_And where
+    add_and    = hsw_formula_add_and . castPtr
+    add_or     = hsw_formula_add_or . castPtr
+    add_not    = hsw_formula_add_not . castPtr
+    add_forall = hsw_formula_add_forall . castPtr
+    add_exists = hsw_formula_add_exists . castPtr
+    convert_to_and = return
+    finalize   = hsw_formula_finalize . castPtr
+
+-- C_Quantifier is a subclass of C_Form and implements all its methods.
+-- Consequently, we simply cast to C_Form
+instance Logical C_Quantifier where
+    add_and    = hsw_formula_add_and . castPtr
+    add_or     = hsw_formula_add_or . castPtr
+    add_not    = hsw_formula_add_not . castPtr
+    add_forall = hsw_formula_add_forall . castPtr
+    add_exists = hsw_formula_add_exists . castPtr
+    convert_to_and = hsw_formula_to_and . castPtr
+    finalize   = hsw_formula_finalize . castPtr
+
+-- Used for freeing data that was allocated in C
+foreign import ccall safe free :: Ptr a -> IO ()
+
+foreign import ccall safe hsw_new_relation
+    :: CInt -> CInt -> IO C_Relation
+foreign import ccall safe hsw_new_set
+    :: CInt -> IO C_Relation
+foreign import ccall safe hsw_free_relation
+    :: C_Relation -> IO ()
+foreign import ccall "&hsw_free_relation" ptr_to_free_relation
+    :: FunPtr (C_Relation -> IO ())
+foreign import ccall safe hsw_relation_show
+    :: C_Relation -> IO CString
+foreign import ccall safe hsw_num_input_vars
+    :: C_Relation -> IO CInt
+foreign import ccall safe hsw_num_output_vars
+    :: C_Relation -> IO CInt
+foreign import ccall safe hsw_num_set_vars
+    :: C_Relation -> IO CInt
+foreign import ccall safe hsw_input_var
+    :: C_Relation -> CInt -> IO C_Var
+foreign import ccall safe hsw_output_var
+    :: C_Relation -> CInt -> IO C_Var
+foreign import ccall safe hsw_set_var
+    :: C_Relation -> CInt -> IO C_Var
+foreign import ccall safe hsw_is_lower_bound_satisfiable
+    :: C_Relation -> IO Bool
+foreign import ccall safe hsw_is_upper_bound_satisfiable
+    :: C_Relation -> IO Bool
+foreign import ccall safe hsw_is_obvious_tautology
+    :: C_Relation -> IO Bool
+foreign import ccall safe hsw_is_definite_tautology
+    :: C_Relation -> IO Bool
+foreign import ccall safe hsw_is_exact
+    :: C_Relation -> IO Bool
+foreign import ccall safe hsw_is_inexact
+    :: C_Relation -> IO Bool
+foreign import ccall safe hsw_is_unknown
+    :: C_Relation -> IO Bool
+foreign import ccall safe hsw_upper_bound
+    :: C_Relation -> IO C_Relation
+foreign import ccall safe hsw_lower_bound
+    :: C_Relation -> IO C_Relation
+foreign import ccall safe hsw_equal
+    :: C_Relation -> C_Relation -> IO CInt
+foreign import ccall safe hsw_union
+    :: C_Relation -> C_Relation -> IO C_Relation
+foreign import ccall safe hsw_intersection
+    :: C_Relation -> C_Relation -> IO C_Relation
+foreign import ccall safe hsw_composition
+    :: C_Relation -> C_Relation -> IO C_Relation
+foreign import ccall safe hsw_restrict_domain
+    :: C_Relation -> C_Relation -> IO C_Relation
+foreign import ccall safe hsw_restrict_range
+    :: C_Relation -> C_Relation -> IO C_Relation
+foreign import ccall safe hsw_difference
+    :: C_Relation -> C_Relation -> IO C_Relation
+foreign import ccall safe hsw_cross_product
+    :: C_Relation -> C_Relation -> IO C_Relation
+foreign import ccall safe hsw_gist
+    :: C_Relation -> C_Relation -> CInt -> IO C_Relation
+foreign import ccall safe hsw_transitive_closure
+    :: C_Relation -> IO C_Relation
+foreign import ccall safe hsw_domain
+    :: C_Relation -> IO C_Relation
+foreign import ccall safe hsw_range
+    :: C_Relation -> IO C_Relation
+foreign import ccall safe hsw_inverse
+    :: C_Relation -> IO C_Relation
+foreign import ccall safe hsw_complement
+    :: C_Relation -> IO C_Relation
+foreign import ccall safe hsw_deltas
+    :: C_Relation -> IO C_Relation
+foreign import ccall safe hsw_approximate
+    :: C_Relation -> IO C_Relation
+
+foreign import ccall safe hsw_relation_add_and
+    :: C_Relation -> IO C_Form
+foreign import ccall safe hsw_relation_add_or
+    :: C_Relation -> IO C_Form
+foreign import ccall safe hsw_relation_add_not
+    :: C_Relation -> IO C_Form
+foreign import ccall safe hsw_relation_add_forall
+    :: C_Relation -> IO C_Quantifier
+foreign import ccall safe hsw_relation_add_exists
+    :: C_Relation -> IO C_Quantifier
+foreign import ccall safe hsw_relation_finalize
+    :: C_Relation -> IO ()
+
+-- These functions take formula pointer arguments
+foreign import ccall safe hsw_formula_add_and
+    :: C_Form -> IO C_Form
+foreign import ccall safe hsw_formula_add_or
+    :: C_Form -> IO C_Form
+foreign import ccall safe hsw_formula_add_not
+    :: C_Form -> IO C_Form
+foreign import ccall safe hsw_formula_add_forall
+    :: C_Form -> IO C_Quantifier
+foreign import ccall safe hsw_formula_add_exists
+    :: C_Form -> IO C_Quantifier
+foreign import ccall safe hsw_formula_finalize
+    :: C_Form -> IO ()
+
+foreign import ccall safe hsw_declaration_declare
+    :: C_Quantifier -> IO C_Var
+
+-- If the argument is a C_And, the argument is returned;
+-- otherwise, add_and is called
+foreign import ccall safe hsw_formula_to_and
+    :: C_Form -> IO C_And
+
+foreign import ccall safe hsw_add_constraint
+    :: C_And -> Bool -> CInt -> Ptr CInt -> Ptr C_Var -> CInt -> IO ()
+
+foreign import ccall safe separate_relation_dimensions
+    :: Ptr C_Relation -> C_Relation -> IO ()
+
+-- Look at the internal representation of a set
+foreign import ccall safe hsw_query_dnf
+    :: C_Relation -> IO C_DNF_Iterator
+foreign import ccall safe hsw_dnf_iterator_next
+    :: C_DNF_Iterator -> IO C_Conjunct
+foreign import ccall safe hsw_dnf_iterator_free
+    :: C_DNF_Iterator -> IO ()
+
+-- For inspecting Omega data structures
+foreign import ccall safe hsw_get_conjunct_variables
+    :: C_Conjunct -> IO (C_Tuple_Iterator C_Var)
+foreign import ccall safe hsw_tuple_iterator_next
+    :: (C_Tuple_Iterator (Ptr a)) -> IO (Ptr a)
+foreign import ccall safe hsw_tuple_iterator_free
+    :: (C_Tuple_Iterator a) -> IO ()
+
+foreign import ccall safe hsw_get_eqs
+    :: C_Conjunct -> IO C_EQ_Iterator
+foreign import ccall safe hsw_eqs_next
+    :: C_EQ_Iterator -> IO C_EQ_Handle
+foreign import ccall safe hsw_eqs_free
+    :: C_EQ_Iterator -> IO ()
+foreign import ccall safe hsw_eq_handle_free
+    :: C_EQ_Handle -> IO ()
+
+foreign import ccall safe hsw_get_geqs
+    :: C_Conjunct -> IO C_GEQ_Iterator
+foreign import ccall safe hsw_geqs_next
+    :: C_GEQ_Iterator -> IO C_GEQ_Handle
+foreign import ccall safe hsw_geqs_free
+    :: C_GEQ_Iterator -> IO ()
+foreign import ccall safe hsw_geq_handle_free
+    :: C_GEQ_Handle -> IO ()
+
+foreign import ccall safe hsw_constraint_get_const
+    :: Ptr a -> IO #{type coefficient_t}
+foreign import ccall safe hsw_constraint_get_coefficients
+    :: Ptr a -> IO C_Constr_Vars_Iter
+foreign import ccall safe hsw_constr_vars_next
+    :: Ptr Coefficient -> C_Constr_Vars_Iter -> IO Bool 
+foreign import ccall safe hsw_constr_vars_free
+    :: C_Constr_Vars_Iter -> IO ()
+
+
+-- For debugging
+foreign import ccall safe hsw_debug_print_eq
+    :: C_EQ_Handle -> IO ()
+foreign import ccall safe hsw_debug_print_geq
+    :: C_GEQ_Handle -> IO ()
+
+-------------------------------------------------------------------------------
+-- Marshalling from Omega Library to Haskell
+
+-- A C++ iterator
+class Iterator i a | i -> a where next :: i -> IO a
+
+instance Iterator C_DNF_Iterator C_Conjunct where
+    next = hsw_dnf_iterator_next
+
+instance Iterator (C_Tuple_Iterator (Ptr a)) (Ptr a) where
+    next = hsw_tuple_iterator_next
+
+instance Iterator C_EQ_Iterator C_EQ_Handle where
+    next = hsw_eqs_next
+
+instance Iterator C_GEQ_Iterator C_GEQ_Handle where
+    next = hsw_geqs_next
+
+-- Imperatively accumulate over the contents of the iterator
+foreach :: Iterator i (Ptr b) => (a -> Ptr b -> IO a) -> a -> i -> IO a
+foreach f x iter = visit x
+    where
+      visit x = do y <- next iter
+                   if y == nullPtr
+                     then return x
+                     else visit =<< f x y
+
+-- Iterate through each conjunct in a DNF clause
+iterateDNF :: (a -> C_Conjunct -> IO a) -> a -> C_Relation -> IO a
+iterateDNF f x rel = do
+  iter <- hsw_query_dnf rel
+  x' <- foreach f x iter
+  hsw_dnf_iterator_free iter
+  return x'
+
+-- Iterate through the variables in a conjunct
+iterateConjVars :: (a -> C_Var -> IO a) -> a -> C_Conjunct -> IO a
+iterateConjVars f x conj = do
+  iter <- hsw_get_conjunct_variables conj
+  x' <- foreach f x iter
+  hsw_tuple_iterator_free iter
+  return x'
+
+-- Iterate through the equality constraints in a conjunct
+iterateEqs :: (a -> C_EQ_Handle -> IO a) -> a -> C_Conjunct -> IO a
+iterateEqs f x conj = do
+  iter <- hsw_get_eqs conj
+  x' <- foreach wrapped_f x iter
+  hsw_eqs_free iter
+  return x'
+    where
+      -- This wrapper just makes sure the handle is freed after use
+      wrapped_f x eqHdl = do 
+        x' <- f x eqHdl
+        hsw_eq_handle_free eqHdl
+        return x'
+
+-- Iterate through the inequality constraints in a conjunct
+iterateGeqs :: (a -> C_GEQ_Handle -> IO a) -> a -> C_Conjunct -> IO a
+iterateGeqs f x conj = do
+  iter <- hsw_get_geqs conj
+  x' <- foreach wrapped_f x iter
+  hsw_geqs_free iter
+  return x'
+    where
+      -- This wrapper just makes sure the handle is freed after use
+      wrapped_f x geqHdl = do 
+        x' <- f x geqHdl
+        hsw_geq_handle_free geqHdl
+        return x'
+
+-- Read the coefficients from a Constraint
+peekConstraintVars :: Constraint a => Ptr a -> IO [Coefficient]
+peekConstraintVars cst = do
+  iter <- hsw_constraint_get_coefficients cst
+
+  -- Allocate temporary storage on the C side for some data
+  c_var_info <- ForeignAlloc.malloc
+
+  -- Traverse and pull out values
+  coeffs <- getCoefficients iter c_var_info []
+
+  -- Free the temporary storage and the iterator
+  ForeignAlloc.free c_var_info
+  hsw_constr_vars_free iter
+
+  return coeffs
+    where
+      getCoefficients iter c_var_info coeffs = do
+
+        -- Read one coefficient
+        ok <- hsw_constr_vars_next c_var_info iter
+
+        -- If it returned false, we're done
+        if not ok then return coeffs else do
+        coeff <- peek c_var_info
+        getCoefficients iter c_var_info (coeff:coeffs)
+
+-- Read the list of input variables [N, N-1 ... 1].
+-- This will probably crash if the number of variables is not specified
+-- correctly.
+peekInputVars, peekOutputVars, peekSetVars
+    :: CInt -> C_Relation -> IO [VarHandle]
+peekInputVars n rel =
+    mapM (liftM VarHandle . hsw_input_var rel) [n, n - 1 .. 1]
+
+peekOutputVars n rel =
+    mapM (liftM VarHandle . hsw_output_var rel) [n, n - 1 .. 1]
+
+peekSetVars n rel =
+    mapM (liftM VarHandle . hsw_set_var rel) [n, n - 1 .. 1]
+
+-- Helper function to read a constraint.
+queryConstraint :: Constraint c =>
+                   ([Coefficient] -> Int -> a -> a) -- Accumulating function
+                -> a            -- Initial value
+                -> Ptr c        -- Constraint to query
+                -> IO a
+queryConstraint f acc eq = do
+  coefficients <- peekConstraintVars eq
+  constant <- hsw_constraint_get_const eq
+  return $ f coefficients (fromIntegral constant) acc
+
+-------------------------------------------------------------------------------
+-- Exported interface
+
+-- | Data types containing Presburger formulae.
+class Presburger a where
+    -- | Extract the pointer from a formula
+    pPtr :: a -> ForeignPtr Relation
+
+    -- | Test whether two sets or relations have the same arity
+    sameArity :: a -> a -> Bool
+
+    -- | Convert a raw pointer to an OmegaSet or OmegaRel
+    fromPtr :: Ptr Relation -> IO a
+
+-- Use a wrapped relation or set
+withPresburger :: Presburger a => a -> (C_Relation -> IO b) -> IO b
+withPresburger p = withForeignPtr (pPtr p)
+
+-- Use two wrapped relations or sets
+withPresburger2 :: (Presburger a, Presburger b) =>
+                   a -> b -> (C_Relation -> C_Relation -> IO c) -> IO c
+withPresburger2 p q f = withForeignPtr (pPtr p) $ \ptr ->
+                        withForeignPtr (pPtr q) $ \ptr2 ->
+                        f ptr ptr2
+
+-- | A set of points in Z^n.
+-- This is a wrapper around the Omega library's Relation type.  
+data OmegaSet = OmegaSet { sPtr :: {-# UNPACK #-} !(ForeignPtr Relation)
+                         , sDom :: [VarHandle]
+                         }
+
+instance Show OmegaSet where
+    show rel = unsafePerformIO $ withPresburger rel $ \ptr -> do
+        -- Call hsw_relation_show to get a C string, then convert to String
+        cStr <- hsw_relation_show ptr
+        str  <- peekCString cStr
+        free cStr
+        return str
+
+instance Presburger OmegaSet where
+    pPtr = sPtr
+
+    sameArity s1 s2 =
+      length (sDom s1) == length (sDom s2)
+
+    fromPtr ptr = do
+      numVars <- hsw_num_set_vars ptr
+      varIDs <- peekSetVars numVars ptr
+      wrapOmegaSet ptr varIDs
+
+-- Convert a raw set pointer to an OmegaSet
+wrapOmegaSet :: C_Relation -> [VarHandle] -> IO OmegaSet
+wrapOmegaSet ptr dom = do
+  foreignptr <- newForeignPtr ptr_to_free_relation ptr
+  return $! OmegaSet { sPtr = foreignptr
+                     , sDom = dom
+                     }
+
+-- | Create an Omega set.  The first parameter is the number of dimensions
+-- the set inhabits.  The second parameter builds a formula
+-- defining the set's members. The set's members are those points
+-- for which the formula evaluates to True.
+newOmegaSet :: Int              -- ^ Dimensionality of the space that the set
+                                -- inhabits
+            -> ([VarHandle] -> Formula) -- ^ Set members
+            -> IO OmegaSet
+newOmegaSet numVars init = do
+  rel <- hsw_new_set (fromIntegral numVars)
+
+  -- Look up the ID for each variable in the tuple.  Variables are ordered
+  -- from last to first because the last variable is "innermost," has
+  -- de Bruijn index 1, and belongs at position 1 in the list.
+  freeVarIDs <- peekSetVars (fromIntegral numVars) rel
+
+  runFD (init freeVarIDs) rel
+  wrapOmegaSet rel freeVarIDs
+
+-- | Inspect a set's low-level representation directly.  This function
+-- takes care of data structure traversal and relies on other routines to
+-- interpret the data.
+--
+-- All three accumulating functions take the set variables as their
+-- first parameter, and any existentially quantified variables as
+-- their second parameter.  The set variables are returned along with
+-- a result value.
+queryDNFSet :: ([VarHandle] -> [VarHandle] -> [Coefficient] -> Int -> a -> a)
+               -- ^ Accumulating function for equality constraints
+            -> a                -- ^ Initial value for equality constraints
+            -> ([VarHandle] -> [VarHandle] -> [Coefficient] -> Int -> b -> b)
+               -- ^ Accumulating function for inequality constraints
+            -> b                -- ^ Initial value for inequality constraints
+            -> ([VarHandle] -> [VarHandle] -> a -> b -> c -> c)
+               -- ^ Accumulating function for conjuncts
+            -> c                -- ^ Initial value for conjuncts
+            -> OmegaSet         -- ^ Set to query
+            -> IO ([VarHandle], c)
+queryDNFSet readEq unitEq readGeq unitGeq readConj unitConj s = do
+    conjuncts <- withPresburger s $ iterateDNF doConjunct unitConj
+    return (sDom s, conjuncts)
+    where
+      doConjunct acc conjunct = do
+        -- Find existentially bound variables in this conjunct, which
+        -- Omega calls "wildcard variables"
+        wildcardVars <- iterateConjVars findWildcards [] conjunct
+        let wc = map VarHandle wildcardVars
+
+        -- For each EQ relation, get the relation
+        eqs <- iterateEqs (queryConstraint $ readEq (sDom s) wc)
+               unitEq conjunct
+
+        -- For each GE relation, get the relation
+        geqs <- iterateGeqs (queryConstraint $ readGeq (sDom s) wc)
+                unitGeq conjunct
+
+        return $ readConj (sDom s) wc eqs geqs acc
+
+      findWildcards acc var =
+          -- Is this an input variable?
+          case findIndex (var ==) (map unVarHandle $ sDom s)
+          of Just n  -> return acc
+             Nothing -> -- Otherwise, assume it's a wildcard
+                        -- FIXME: call into C to check the variable's kind
+                        return $ var : acc
+
+-- | A relation from points in a /domain/ Z^m
+-- to points in a /range/ Z^n.
+-- This is a wrapper around the Omega library's Relation type.
+--
+-- A relation can be considered as just a set of points in Z^(m+n).
+-- However, many routines treat the domain and range differently.
+
+data OmegaRel = OmegaRel { rPtr :: {-# UNPACK #-} !(ForeignPtr Relation)
+                         , rDom :: [VarHandle]
+                         , rRng :: [VarHandle]
+                         }
+
+instance Show OmegaRel where
+    show rel = unsafePerformIO $ withPresburger rel $ \ptr -> do
+        -- Call hsw_relation_show to get a C string, then convert to String
+        cStr <- hsw_relation_show ptr
+        str  <- peekCString cStr
+        free cStr
+        return str
+
+instance Presburger OmegaRel where
+    pPtr = rPtr
+
+    sameArity r1 r2 =
+      length (rDom r1) == length (rDom r2) &&
+      length (rRng r1) == length (rRng r2)
+
+    fromPtr ptr = do
+      numOutputs <- hsw_num_output_vars ptr
+      outputVarIDs <- peekOutputVars numOutputs ptr
+
+      numInputs <- hsw_num_input_vars ptr
+      inputVarIDs <- peekInputVars numInputs ptr
+
+      wrapOmegaRel ptr inputVarIDs outputVarIDs
+
+-- Convert a raw relation pointer to an OmegaSet
+wrapOmegaRel :: C_Relation -> [VarHandle] -> [VarHandle] -> IO OmegaRel
+wrapOmegaRel ptr dom rng = do
+  foreignptr <- newForeignPtr ptr_to_free_relation ptr
+  return $! OmegaRel { rPtr = foreignptr
+                     , rDom = dom
+                     , rRng = rng }
+
+-- | Create an Omega relation.  The first two parameters are the number
+-- of dimensions of the domain and range, respectively.  The third parameter
+-- builds a formula defining the relation.  Two points are related iff the
+-- formula evaluates to True on those points.
+newOmegaRel :: Int              -- ^ Dimensionality of the domain
+            -> Int              -- ^ Dimensionality of the range
+            -> ([VarHandle] -> [VarHandle] -> Formula)
+                                -- ^ The relation
+            -> IO OmegaRel
+newOmegaRel numInputs numOutputs init = do
+  rel <- hsw_new_relation (fromIntegral numInputs) (fromIntegral numOutputs)
+
+  -- Look up the IDs for the input and output variables.
+  outputVarIds <- peekOutputVars (fromIntegral numOutputs) rel
+  inputVarIds <- peekInputVars (fromIntegral numInputs) rel
+
+  runFD (init inputVarIds outputVarIds) rel
+  wrapOmegaRel rel inputVarIds outputVarIds
+
+-- | Inspect a relation's low-level representation directly.  This function
+-- takes care of data structure traversal and relies on other routines to
+-- interpret the data.
+--
+-- All three accumulating functions take the relation's input and
+-- output variables as their first and second parameters, respectively,
+-- and any existentially quantified variables as
+-- their second parameter.  The relation's input and output variables are
+-- returned along with a result value.
+queryDNFRelation :: ([VarHandle] -> [VarHandle] -> [VarHandle] -> [Coefficient] -> Int -> a -> a)
+                    -- ^ Accumulating function for equality constraints
+                 -> a           -- ^ Initial value for equality constraints
+                 -> ([VarHandle] -> [VarHandle] -> [VarHandle] -> [Coefficient] -> Int -> b -> b)
+                    -- ^ Accumulating function for inequality constraints
+                 -> b           -- ^ Initial value for inequality constraints
+                 -> ([VarHandle] -> [VarHandle] -> [VarHandle] -> a -> b -> c -> c)
+                    -- ^ Accumulating function for conjuncts
+                 -> c           -- ^ Initial value for conjuncts
+                 -> OmegaRel    -- ^ Relation to query
+                 -> IO ([VarHandle], [VarHandle], c) -- ^ Input variables, output variables, and result
+queryDNFRelation readEq unitEq readGeq unitGeq readConj unitConj r = do
+    conjuncts <- withPresburger r $ iterateDNF doConjunct unitConj
+    return (rDom r, rRng r, conjuncts)
+    where
+      doConjunct acc conjunct = do
+        -- Find existentially bound variables in this conjunct, which
+        -- Omega calls "wildcard variables"
+        wildcardVars <- iterateConjVars findWildcards [] conjunct
+        let wc = map VarHandle wildcardVars
+
+        -- For each EQ relation, get the relation
+        eqs <- iterateEqs (queryConstraint $ readEq (rDom r) (rRng r) wc)
+               unitEq conjunct
+
+        -- For each GE relation, get the relation
+        geqs <- iterateGeqs (queryConstraint $ readGeq (rDom r) (rRng r) wc)
+                unitGeq conjunct
+
+        return $ readConj (rDom r) (rRng r) wc eqs geqs acc
+
+      findWildcards acc var =
+          -- Is this an input variable?
+          case findIndex (var ==) (map unVarHandle $ rDom r ++ rRng r)
+          of Just n  -> return acc
+             Nothing -> -- Otherwise, assume it's a wildcard
+                        -- FIXME: call into C to check the variable's kind
+                        return $ var : acc
+
+-- | Get a list of relations, one per output variable, with the same
+-- input and output dimensions as the original, but whose constraints
+-- mention only one output variable and no existential constraints.
+--
+-- This function is needed to create a high-level Rel from a low-level
+-- OmegaRel.
+separateRelationDimensions :: OmegaRel -> IO [OmegaRel]
+separateRelationDimensions r = do
+  -- Allocate an array to store outputs
+  allocaArray numOutputs $ \outputArray -> do
+    -- Call into C
+    withPresburger r $ \ptr -> separate_relation_dimensions outputArray ptr
+
+    -- Wrap each output as a relation
+    mapM readRelation =<< peekArray numOutputs outputArray
+
+    where
+      numInputs  = length $ rDom r
+      numOutputs = length $ rRng r
+
+      readRelation rel = do
+        -- Look up the IDs for the input and output variables.
+        outputVarIds <- peekOutputVars (fromIntegral numOutputs) rel
+        inputVarIds <- peekInputVars (fromIntegral numInputs) rel
+
+        wrapOmegaRel rel inputVarIds outputVarIds
+
+-------------------------------------------------------------------------------
+-- Queries
+
+-- | Determine a lower bound on whether the formula is satisfiable.
+-- The lower bound is based on treating all UNKNOWN constraints as false.
+lowerBoundSatisfiable :: Presburger a => a -> IO Bool
+
+-- | Determine an upper bound on whether the formula is satisfiable.
+-- The lower bound is based on treating all UNKNOWN constraints as false.
+upperBoundSatisfiable :: Presburger a => a -> IO Bool
+
+-- | Use simple, fast tests to decide whether the formula is a tautology.
+obviousTautology      :: Presburger a => a -> IO Bool
+
+-- | True if the formula is a tautology.
+definiteTautology     :: Presburger a => a -> IO Bool
+
+-- | True if the formula has no UNKNOWN constraints.
+exact                 :: Presburger a => a -> IO Bool
+
+-- | True if the formula has UNKNOWN constraints.
+inexact               :: Presburger a => a -> IO Bool
+
+-- | True if the formula is UNKNOWN.
+unknown               :: Presburger a => a -> IO Bool
+
+lowerBoundSatisfiable rel = withPresburger rel hsw_is_lower_bound_satisfiable
+upperBoundSatisfiable rel = withPresburger rel hsw_is_upper_bound_satisfiable
+obviousTautology rel      = withPresburger rel hsw_is_obvious_tautology
+definiteTautology rel     = withPresburger rel hsw_is_definite_tautology
+exact rel                 = withPresburger rel hsw_is_exact
+inexact rel               = withPresburger rel hsw_is_inexact
+unknown rel               = withPresburger rel hsw_is_unknown
+
+-------------------------------------------------------------------------------
+-- Creating new sets and relations from old ones
+
+-- | Compute the upper bound of a set or relation by setting all UNKNOWN
+-- constraints to true.
+upperBound :: Presburger a => a -> IO a
+upperBound rel = fromPtr =<< withPresburger rel hsw_upper_bound
+
+-- | Compute the lower bound of a set or relation by setting all UNKNOWN
+-- constraints to false.
+lowerBound :: Presburger a => a -> IO a
+lowerBound rel = fromPtr =<< withPresburger rel hsw_lower_bound
+
+-- | Test whether two sets or relations are equal.
+-- The sets or relations must have the same arity.
+--
+-- The answer is precise if both arguments are 'exact'.
+-- If either argument is inexact, this function returns @False@.
+equal :: Presburger a => a -> a -> IO Bool
+equal rel1 rel2
+    | sameArity rel1 rel2 = do
+          eq <- withPresburger2 rel1 rel2 hsw_equal
+          return $! eq /= 0
+    | otherwise = error "equal: arguments have different arities"
+
+-- | Compute the union of two sets or relations.  The sets or relations
+-- must have the same arity.
+union :: Presburger a => a -> a -> IO a
+union rel1 rel2
+    | sameArity rel1 rel2 =
+          fromPtr =<< withPresburger2 rel1 rel2 hsw_union 
+    | otherwise = error "union: arguments have different arities"
+
+-- | Compute the intersection of two sets or relations.  The sets or relations
+-- must have the same arity.
+intersection :: Presburger a => a -> a -> IO a
+intersection rel1 rel2
+    | sameArity rel1 rel2 =
+          fromPtr =<< withPresburger2 rel1 rel2 hsw_intersection
+    | otherwise = error "intersection: arguments have different arities"
+
+-- | Compute the composition of two sets or relations.  The
+-- first relation's domain must be the same dimension as the second's range.
+composition :: OmegaRel -> OmegaRel -> IO OmegaRel
+composition rel1 rel2
+    | length (rDom rel1) == length (rRng rel2) =
+          fromPtr =<< withPresburger2 rel1 rel2 hsw_composition
+    | otherwise = error "composition: argument arities do not agree"
+
+restrictDomain :: OmegaRel -> OmegaSet -> IO OmegaRel
+restrictDomain rel1 set
+    | length (rDom rel1) == length (sDom set) =
+          fromPtr =<< withPresburger2 rel1 set hsw_restrict_domain
+    | otherwise = error "restrictDomain: argument arities do not agree"
+
+restrictRange :: OmegaRel -> OmegaSet -> IO OmegaRel
+restrictRange rel1 set
+    | length (rDom rel1) == length (sDom set) =
+          fromPtr =<< withPresburger2 rel1 set hsw_restrict_range
+    | otherwise = error "restrictRange: argument arities do not agree"
+
+difference :: Presburger a => a -> a -> IO a
+difference rel1 rel2
+    | sameArity rel1 rel2 =
+        fromPtr =<< withPresburger2 rel1 rel2 hsw_difference
+    | otherwise = error "difference: arguments have different arities"
+
+crossProduct :: OmegaSet -> OmegaSet -> IO OmegaRel
+crossProduct set1 set2 =
+    fromPtr =<< withPresburger2 set1 set2 hsw_cross_product
+
+-- | Get the gist of a set or relation, given some background truth.  The
+-- gist operator uses heuristics to make a set or relation simpler, while
+-- still retaining sufficient information to regenerate the original by
+-- re-introducing the background truth.  The sets or relations
+-- must have the same arity.
+--
+-- Given @x@ computed by
+--
+-- > x <- intersection given =<< gist effort r given
+--
+-- we have @x == r@.
+gist :: Presburger a => Effort -> a -> a -> IO a
+gist effort rel given
+    | sameArity rel given =
+        withPresburger2 rel given $ \ptr ptrGiven ->
+          fromPtr =<< hsw_gist ptr ptrGiven (fromIntegral $ fromEnum effort)
+    | otherwise = error "gist: arguments have different arities"
+
+-- | Get the transitive closure of a relation.  In some cases, the transitive
+-- closure cannot be computed exactly, in which case a lower bound is
+-- returned.
+transitiveClosure :: OmegaRel -> IO OmegaRel
+transitiveClosure rel = fromPtr =<< withPresburger rel hsw_transitive_closure
+
+-- | Get the domain of a relation.
+domain :: OmegaRel -> IO OmegaSet
+domain rel = fromPtr =<< withPresburger rel hsw_domain
+
+-- | Get the range of a relation.
+range :: OmegaRel -> IO OmegaSet
+range rel = fromPtr =<< withPresburger rel hsw_range
+
+-- | Get the inverse of a relation.
+inverse :: OmegaRel -> IO OmegaRel
+inverse rel = fromPtr =<< withPresburger rel hsw_inverse
+
+-- | Get the complement of a set or relation.
+complement :: Presburger a => a -> IO a
+complement rel = fromPtr =<< withPresburger rel hsw_complement
+
+-- | Get the deltas of a relation.
+-- The relation's input dimensionality must be the same as its output
+-- dimensionality.
+deltas :: OmegaRel -> IO OmegaSet
+deltas rel
+    | length (rDom rel) == length (rRng rel) =
+        fromPtr =<< withPresburger rel hsw_deltas
+    | otherwise =
+        error "deltas: relation has different input and output dimensionality"
+
+-- | Approximate a set or relation by allowing all existentially quantified
+-- variables to take on rational values.  This allows these variables to be
+-- eliminated from the formula.
+approximate :: Presburger a => a -> IO a
+approximate rel = fromPtr =<< withPresburger rel hsw_approximate
+
+-------------------------------------------------------------------------------
+-- Formulae
+
+-- | A boolean-valued Presburger formula.
+
+-- This is actually a function that builds a Presburger formula.
+newtype Formula = FD {runFD :: forall a. Logical a => a -> IO ()}
+
+-- | Logical conjunction (and).
+conjunction :: [Formula] -> Formula
+conjunction formulaDefs = FD $ \f -> do
+  newF <- add_and f
+  mapM_ (\func -> runFD func newF) formulaDefs
+  finalize newF
+
+-- | Logical disjunction (or).
+disjunction :: [Formula] -> Formula
+disjunction formulaDefs = FD $ \f -> do
+  newF <- add_or f
+  mapM_ (\func -> runFD func newF) formulaDefs
+  finalize newF
+
+-- | Logical negation.
+negation :: Formula -> Formula
+negation formulaDef = FD $ \f -> do
+  newF <- add_not f
+  runFD formulaDef newF
+  finalize newF
+
+-- | Universal quantification.  The 'VarHandle' parameter is the variable
+-- bound by the quantifier.
+qForall :: (VarHandle -> Formula) -> Formula
+qForall makeBody = FD $ \f -> do
+  newFormula <- add_forall f
+  localVar <- hsw_declaration_declare newFormula
+  runFD (makeBody (VarHandle localVar)) newFormula
+  finalize newFormula
+
+-- | Existential quantification.  The 'VarHandle' parameter is the variable
+-- bound by the quantifier.
+qExists :: (VarHandle -> Formula) -> Formula
+qExists makeBody = FD $ \f -> do
+  newFormula <- add_exists f
+  localVar <- hsw_declaration_declare newFormula
+  runFD (makeBody (VarHandle localVar)) newFormula
+  finalize newFormula
+
+-- Add an equality or inequality constraint to a conjunction.
+addConstraint :: Bool -> [Coefficient] -> Int -> C_And -> IO ()
+addConstraint kind terms constant formula = do
+  let numTerms     = length terms
+      numTermsCInt = fromIntegral numTerms
+      constantCInt = fromIntegral constant
+      coefficients = map (fromIntegral . coeffValue) terms
+      variables    = map ((\(VarHandle h) -> h) . coeffVar) terms
+
+  -- Marshal the coefficients and variables to C as arrays
+  withArray coefficients $ \coeffPtr ->
+      withArray variables $ \varPtr ->
+          -- then, call code to set the constraint
+          hsw_add_constraint formula kind numTermsCInt coeffPtr varPtr constantCInt
+
+-- | Construct an inequation of the form @a*x + b*y + ... + d >= 0@.
+inequality :: [Coefficient] -> Int -> Formula
+inequality terms constant = FD $ \formula ->
+    addConstraint False terms constant =<< convert_to_and formula
+
+-- | Construct an equation of the form @a*x + b*y + ... + d = 0@.
+equality :: [Coefficient] -> Int -> Formula
+equality terms constant = FD $ \formula ->
+    addConstraint True terms constant =<< convert_to_and formula
+
+-- | Truth.
+true :: Formula
+true = equality [] 0
+
+-- | Falsity.
+false :: Formula
+false = equality [] 1
+
+-- | A variable in a formula.
+
+-- These data structures are owned by OmegaSet or OmegaRel instances,
+-- which take care of allocation and deallocation.
+newtype VarHandle = VarHandle { unVarHandle :: C_Var } deriving(Eq)
+
+-- | An integer-valued term @n*v@ in a formula.
+data Coefficient =
+    Coefficient { coeffVar :: {-# UNPACK #-} !VarHandle
+                , coeffValue :: {-# UNPACK #-} !Int
+                }
+
+instance Show Coefficient where
+    show (Coefficient v n) =
+        "(" ++ show n ++ " * " ++ show (unVarHandle v) ++ ")"
+
+instance Storable Coefficient where
+    sizeOf _ = #{size Variable_Info_struct}
+    alignment _ = #{alignof Variable_Info_struct}
+    peek p = do
+      var  <- #{peek Variable_Info_struct, var} p :: IO C_Var
+      coef <- #{peek Variable_Info_struct, coef} p :: IO #{type coefficient_t}
+      return $ Coefficient { coeffVar = VarHandle var
+                           , coeffValue = fromIntegral coef
+                           }
+
diff --git a/Data/Presburger/Omega/Rel.hs b/Data/Presburger/Omega/Rel.hs
new file mode 100644
--- /dev/null
+++ b/Data/Presburger/Omega/Rel.hs
@@ -0,0 +1,351 @@
+
+-- | Relations whose members are represented compactly using a
+-- Presburger arithmetic formula.  This is a high-level interface to
+-- 'OmegaRel'.
+--
+-- This module is intended to be imported qualified, e.g.
+--
+-- > import qualified Data.Presburger.Omega.Rel as WRel
+
+module Data.Presburger.Omega.Rel
+    (Rel,
+     -- * Building relations
+     rel, functionalRel, fromOmegaRel,
+
+     -- * Operations on relations
+     toOmegaRel,
+
+     -- ** Inspecting
+     inputDimension, outputDimension,
+     predicate,
+     lowerBoundSatisfiable,
+     upperBoundSatisfiable,
+     obviousTautology,
+     definiteTautology,
+     exact,
+     inexact,
+     unknown,
+     equal,
+
+     -- ** Bounds
+     upperBound, lowerBound,
+
+     -- ** Binary operations
+     union, intersection, composition, join,
+     restrictDomain, restrictRange,
+     difference, crossProduct,
+     Effort(..),
+     gist,
+
+     -- ** Unary operations
+     transitiveClosure,
+     domain, range,
+     inverse,
+     complement,
+     deltas,
+     approximate
+    )
+where
+
+import System.IO.Unsafe
+
+import Data.Presburger.Omega.Expr
+import qualified Data.Presburger.Omega.LowLevel as L
+import Data.Presburger.Omega.LowLevel(OmegaRel, Effort(..))
+import Data.Presburger.Omega.SetRel
+import qualified Data.Presburger.Omega.Set as Set
+import Data.Presburger.Omega.Set(Set)
+
+-- | A relation from points in a /domain/ Z^m to points in a /range/ Z^n.
+--
+-- A relation can be considered just a set of points in Z^(m+n).  However,
+-- many functions that operate on relations treat the domain and range
+-- differently.
+
+-- Variables are referenced by de Bruijn index.  The order is:
+-- [dom_1, dom_2 ... dom_n, rng_1, rng_2 ... rng_m]
+-- where rng_1 has the lowest index and dom_m the highest.
+data Rel = Rel
+    { relInpDim :: !Int         -- ^ number of variables in the input
+    , relOutDim :: !Int         -- ^ the function from input to output
+    , relFun    :: BoolExp      -- ^ function defining the relation
+    , relOmegaRel :: OmegaRel   -- ^ low-level representation of this relation
+    }
+
+instance Show Rel where
+    -- Generate a call to 'rel'
+    showsPrec n r = showParen (n >= 10) $
+                    showString "rel " .
+                    shows (relInpDim r) .
+                    showChar ' ' .
+                    shows (relOutDim r) .
+                    showChar ' ' .
+                    showsPrec 10 (relFun r)
+        where
+          showChar c = (c:)
+
+-- | Create a relation whose members are defined by a predicate.
+--
+-- The expression should have @m+n@ free variables, where @m@ and @n@ are
+-- the first two parameters.  The first @m@
+-- variables refer to the domain, and the remaining variables refer to
+-- the range.
+
+rel :: Int                      -- ^ Dimensionality of the domain
+    -> Int                      -- ^ Dimensionality of the range
+    -> BoolExp                  -- ^ Predicate defining the relation
+    -> Rel
+rel inDim outDim expr
+    | variablesWithinRange (inDim + outDim) expr =
+        Rel
+        { relInpDim   = inDim
+        , relOutDim   = outDim
+        , relFun      = expr
+        , relOmegaRel = unsafePerformIO $ mkOmegaRel inDim outDim expr
+        }
+    | otherwise = error "rel: Variables out of range"
+
+mkOmegaRel inDim outDim expr =
+    L.newOmegaRel inDim outDim $ \dom rng -> expToFormula (dom ++ rng) expr
+
+-- | Create a relation where each output is a function of the inputs.
+--
+-- Each expression should have @m@ free variables, where @m@
+-- is the first parameter.
+--
+-- For example, the relation @{(x, y) -> (y, x) | x > 0 && y > 0}@ is
+--
+-- > let [x, y] = takeFreeVariables' 2
+-- > in functionalRel 2 [y, x] (conjE [y |>| intE 0, x |>| intE 0])
+
+functionalRel :: Int            -- ^ Dimensionality of the domain
+              -> [IntExp]       -- ^ Function relating domain to range
+              -> BoolExp        -- ^ Predicate restricting the domain
+              -> Rel
+functionalRel dim range domain
+    | all (variablesWithinRange dim) range &&
+      variablesWithinRange dim domain =
+        Rel
+        { relInpDim   = dim
+        , relOutDim   = length range
+        , relFun      = relationPredicate
+        , relOmegaRel = unsafePerformIO $
+                        mkFunctionalOmegaRel dim range domain
+        }
+    | otherwise = error "functionalRel: Variables out of range"
+    where
+      -- construct the expression domain && rangeVar1 == rangeExp1 && ...
+      relationPredicate =
+          conjE (domain : zipWith outputPredicate [dim..] range)
+
+      outputPredicate index expr =
+          varE (nthVariable index) |==| expr
+
+-- To make an omega relation, we combine the range variables and the domain
+-- into one big happy formula, with the conjunction
+-- @domain /\ rangeVar1 == rangeExp1 /\ ... /\ rangeVarN == rangeExpN@.
+
+mkFunctionalOmegaRel :: Int -> [IntExp] -> BoolExp -> IO OmegaRel
+mkFunctionalOmegaRel dim range domain =
+    L.newOmegaRel dim (length range) $ \dom rng ->
+        L.conjunction (domainConstraint dom : rangeConstraints dom rng)
+    where
+      domainConstraint dom = expToFormula dom domain
+
+      rangeConstraints dom rng = zipWith (rangeConstraint dom) range rng
+
+      -- To make a range constraint, we first add the range variable
+      -- as the outermost bound variable, then convert this expression to an
+      -- equality constraint (rangeVar == ...), then convert 
+      rangeConstraint dom expr rngVar =
+          let -- Add the range variable as the outermost bound variable
+              vars = dom ++ [rngVar]
+
+              -- Turn the range formula into an equality constraint
+              -- (rngVar == ...)
+              expr' = expr |==| varE (nthVariable dim)
+
+          in expToFormula vars expr'
+
+-- | Convert an 'OmegaRel' to a 'Rel'.
+fromOmegaRel :: OmegaRel -> IO Rel
+fromOmegaRel orel = do
+  (dim, range, expr) <- relToExpression orel
+  return $ Rel
+             { relInpDim   = dim
+             , relOutDim   = range
+             , relFun      = expr
+             , relOmegaRel = orel
+             }
+
+-- | Internal function to convert an 'OmegaRel' to a 'Rel', when we know
+-- the relation's dimensions.
+omegaRelToRel :: Int -> Int -> OmegaRel -> IO Rel
+omegaRelToRel inpDim outDim orel = return $
+    Rel
+    { relInpDim   = inpDim
+    , relOutDim   = outDim
+    , relFun      = unsafePerformIO $ do (_, _, expr) <- relToExpression orel
+                                         return $ expr
+    , relOmegaRel = orel
+    }
+
+-------------------------------------------------------------------------------
+-- Operations on relations
+
+-- Some helper functions
+useRel :: (OmegaRel -> IO a) -> Rel -> a
+useRel f r = unsafePerformIO $ f $ relOmegaRel r
+
+useRelRel :: (OmegaRel -> IO OmegaRel) -> Int -> Int -> Rel -> Rel
+useRelRel f inpDim outDim r = unsafePerformIO $ do
+  omegaRelToRel inpDim outDim =<< f (relOmegaRel r)
+
+useRel2 :: (OmegaRel -> OmegaRel -> IO a) -> Rel -> Rel -> a
+useRel2 f r1 r2 = unsafePerformIO $ f (relOmegaRel r1) (relOmegaRel r2)
+
+useRel2Rel :: (OmegaRel -> OmegaRel -> IO OmegaRel)
+           -> Int -> Int -> Rel -> Rel -> Rel
+useRel2Rel f inpDim outDim r1 r2 = unsafePerformIO $ do
+  omegaRelToRel inpDim outDim =<< f (relOmegaRel r1) (relOmegaRel r2)
+
+-- | Get the dimensionality of a relation's domain
+inputDimension :: Rel -> Int
+inputDimension = relInpDim
+
+-- | Get the dimensionality of a relation's range
+outputDimension :: Rel -> Int
+outputDimension = relOutDim
+
+-- | Convert a 'Rel' to an 'OmegaRel'.
+toOmegaRel :: Rel -> OmegaRel
+toOmegaRel = relOmegaRel
+
+-- | Get the predicate defining a relation.
+predicate :: Rel -> BoolExp
+predicate = relFun
+
+domain :: Rel -> Set
+domain r = useRel (\ptr -> Set.fromOmegaSet =<< L.domain ptr) r
+
+range :: Rel -> Set
+range r = useRel (\ptr -> Set.fromOmegaSet =<< L.range ptr) r
+
+lowerBoundSatisfiable :: Rel -> Bool
+lowerBoundSatisfiable = useRel L.lowerBoundSatisfiable
+
+upperBoundSatisfiable :: Rel -> Bool
+upperBoundSatisfiable = useRel L.upperBoundSatisfiable
+
+obviousTautology :: Rel -> Bool
+obviousTautology = useRel L.obviousTautology
+
+definiteTautology :: Rel -> Bool
+definiteTautology = useRel L.definiteTautology
+
+exact :: Rel -> Bool
+exact = useRel L.exact
+
+inexact :: Rel -> Bool
+inexact = useRel L.inexact
+
+unknown :: Rel -> Bool
+unknown = useRel L.unknown
+
+upperBound :: Rel -> Rel
+upperBound r = useRelRel L.upperBound (relInpDim r) (relOutDim r) r
+
+lowerBound :: Rel -> Rel
+lowerBound r = useRelRel L.lowerBound (relInpDim r) (relOutDim r) r
+
+-- | Test whether two relations are equal.
+-- The relations must have the same dimension
+-- (@inputDimension r1 == inputDimension r2 && outputDimension r1 == outputDimension r2@),
+-- or an error will be raised.
+--
+-- The answer is precise if both relations are 'exact'.
+-- If either relation is inexact, this function returns @False@.
+equal :: Rel -> Rel -> Bool
+equal = useRel2 L.equal
+
+-- | Union of two relations.
+-- The relations must have the same dimension
+-- (@inputDimension r1 == inputDimension r2 && outputDimension r1 == outputDimension r2@),
+-- or an error will be raised.
+union :: Rel -> Rel -> Rel
+union s1 s2 = useRel2Rel L.union (relInpDim s1) (relOutDim s1) s1 s2
+
+-- | Intersection of two relations.
+-- The relations must have the same dimension
+-- (@inputDimension r1 == inputDimension r2 && outputDimension r1 == outputDimension r2@),
+-- or an error will be raised.
+intersection :: Rel -> Rel -> Rel
+intersection s1 s2 =
+    useRel2Rel L.intersection (relInpDim s1) (relOutDim s1) s1 s2
+
+-- | Composition of two relations.
+-- The second relation's output must be the same size as the first's input
+-- (@outputDimension r2 == inputDimension r1@),
+-- or an error will be raised.
+composition :: Rel -> Rel -> Rel
+composition s1 s2 =
+    useRel2Rel L.composition (relInpDim s2) (relOutDim s1) s1 s2
+
+-- | Same as 'composition', with the arguments swapped.
+join :: Rel -> Rel -> Rel
+join r1 r2 = composition r2 r1
+
+restrictDomain :: Rel -> Set -> Rel
+restrictDomain r s = unsafePerformIO $
+  omegaRelToRel (relInpDim r) (relOutDim r) =<<
+  L.restrictDomain (relOmegaRel r) (Set.toOmegaSet s)
+
+restrictRange :: Rel -> Set -> Rel
+restrictRange r s = unsafePerformIO $
+  omegaRelToRel (relInpDim r) (relOutDim r) =<<
+  L.restrictRange (relOmegaRel r) (Set.toOmegaSet s)
+
+-- | Difference of two relations.
+-- The relations must have the same dimension
+-- (@inputDimension r1 == inputDimension r2 && outputDimension r1 == outputDimension r2@),
+-- or an error will be raised.
+difference :: Rel -> Rel -> Rel
+difference s1 s2 =
+    useRel2Rel L.difference (relInpDim s1) (relOutDim s1) s1 s2
+
+-- | Cross product of two sets.
+crossProduct :: Set -> Set -> Rel
+crossProduct s1 s2 = unsafePerformIO $
+  omegaRelToRel (Set.dimension s1) (Set.dimension s2) =<<
+  L.crossProduct (Set.toOmegaSet s1) (Set.toOmegaSet s2)
+
+-- | Get the gist of a relation, given some background truth.  The
+-- gist operator uses heuristics to simplify the relation while
+-- retaining sufficient information to regenerate the original by
+-- re-introducing the background truth.  The relations must have the
+-- same input dimensions and the same output dimensions.
+--
+-- Given @x@ computed by
+--
+-- > x <- intersection given =<< gist effort r given
+--
+-- we have @x == r@.
+gist :: Effort -> Rel -> Rel -> Rel
+gist effort r1 r2 =
+    useRel2Rel (L.gist effort) (relInpDim r1) (relOutDim r1) r1 r2
+
+transitiveClosure :: Rel -> Rel
+transitiveClosure r =
+    useRelRel L.transitiveClosure (relInpDim r) (relOutDim r) r
+
+inverse :: Rel -> Rel
+inverse s = useRelRel L.inverse (relOutDim s) (relInpDim s) s
+
+complement :: Rel -> Rel
+complement s = useRelRel L.complement (relInpDim s) (relOutDim s) s
+
+deltas :: Rel -> Set
+deltas = useRel (\wrel -> Set.fromOmegaSet =<< L.deltas wrel)
+
+approximate :: Rel -> Rel
+approximate s = useRelRel L.approximate (relInpDim s) (relOutDim s) s
diff --git a/Data/Presburger/Omega/Set.hs b/Data/Presburger/Omega/Set.hs
new file mode 100644
--- /dev/null
+++ b/Data/Presburger/Omega/Set.hs
@@ -0,0 +1,225 @@
+
+-- | Sets whose members are represented compactly using a
+-- Presburger arithmetic formula.  This is a high-level interface to
+-- 'OmegaSet'.
+--
+-- This module is intended to be imported qualified, e.g.
+--
+-- > import qualified Data.Presburger.Omega.Set as WSet
+
+module Data.Presburger.Omega.Set
+    (Set,
+
+     -- * Building sets
+     set, fromOmegaSet,
+
+     -- * Operations on sets
+     toOmegaSet,
+
+     -- ** Inspecting
+     dimension, predicate,
+     lowerBoundSatisfiable,
+     upperBoundSatisfiable,
+     obviousTautology,
+     definiteTautology,
+     exact,
+     inexact,
+     unknown,
+     equal,
+
+     -- ** Bounds
+     upperBound, lowerBound,
+
+     -- ** Binary operations
+     union, intersection, difference,
+     Effort(..),
+     gist,
+
+     -- ** Unary operations
+     complement,
+     approximate
+    )
+where
+
+import System.IO.Unsafe
+
+import Data.Presburger.Omega.Expr
+import qualified Data.Presburger.Omega.LowLevel as L
+import Data.Presburger.Omega.LowLevel(OmegaSet, Effort(..))
+import Data.Presburger.Omega.SetRel
+
+-- | Sets of points in Z^n defined by a formula.
+data Set = Set
+    { setDim      :: !Int      -- ^ the number of variables
+    , setExp      :: BoolExp   -- ^ a predicate defining the set
+    , setOmegaSet :: OmegaSet  -- ^ low-level representation of this set
+    }
+
+instance Show Set where
+    -- Generate a call to 'set'
+    showsPrec n s = showParen (n >= 10) $
+                    showString "set " .
+                    shows (setDim s) .
+                    showChar ' ' .
+                    showsPrec 10 (setExp s)
+
+-- | Create a set whose members are defined by a predicate.
+--
+-- The expression should have one free variable for each dimension.
+--
+-- For example, the set of all points on the plane is
+-- 
+-- >  set 2 trueE
+-- 
+-- The set of all points (x, y, z) where x > y + z is
+-- 
+-- >  set 3 (case takeFreeVariables' 3 of [x,y,z] -> x |>| y |+| z)
+--
+set :: Int                      -- ^ Number of dimensions
+    -> BoolExp                  -- ^ Predicate defining the set
+    -> Set
+set dim expr
+    | variablesWithinRange dim expr =
+        Set
+        { setDim      = dim
+        , setExp      = expr
+        , setOmegaSet = unsafePerformIO $ mkOmegaSet dim expr
+        }
+    | otherwise = error "set: Variables out of range"
+
+mkOmegaSet :: Int -> BoolExp -> IO OmegaSet
+mkOmegaSet dim expr = L.newOmegaSet dim (\vars -> expToFormula vars expr)
+
+-------------------------------------------------------------------------------
+-- Creating sets from Omega sets
+
+-- | Convert an 'OmegaSet' to a 'Set'.
+fromOmegaSet :: OmegaSet -> IO Set
+fromOmegaSet oset = do
+  (dim, expr) <- setToExpression oset
+  return $ Set
+             { setDim      = dim
+             , setExp      = expr
+             , setOmegaSet = oset
+             }
+
+-- | Internal function to convert an 'OmegaSet' to a 'Set', when we know
+-- the set's dimension.  This can avoid actually building the expression
+-- when all we want is the dimension.
+omegaSetToSet :: Int -> OmegaSet -> IO Set
+omegaSetToSet dim oset = return $
+    Set
+    { setDim      = dim
+    , setExp      = unsafePerformIO $ do (_, expr) <- setToExpression oset
+                                         return expr
+    , setOmegaSet = oset
+    }
+
+-------------------------------------------------------------------------------
+-- Using sets
+
+-- First, some helper functions for applying OmegaSet functions to Sets
+
+useSet :: (OmegaSet -> IO a) -> Set -> a
+useSet f s = unsafePerformIO $ f (setOmegaSet s)
+
+useSetSet :: (OmegaSet -> IO OmegaSet) -> Int -> Set -> Set
+useSetSet f dim s = unsafePerformIO $ do
+  omegaSetToSet dim =<< f (setOmegaSet s)
+
+useSet2 :: (OmegaSet -> OmegaSet -> IO a) -> Set -> Set -> a
+useSet2 f s1 s2 = unsafePerformIO $ f (setOmegaSet s1) (setOmegaSet s2)
+
+useSet2Set :: (OmegaSet -> OmegaSet -> IO OmegaSet)
+           -> Int
+           -> Set
+           -> Set
+           -> Set
+useSet2Set f dim s1 s2 = unsafePerformIO $ do
+  omegaSetToSet dim =<< f (setOmegaSet s1) (setOmegaSet s2)
+
+-- | Get the dimensionality of the space a set inhabits
+dimension :: Set -> Int
+dimension = setDim
+
+-- | Get the predicate defining a set's members
+predicate :: Set -> BoolExp
+predicate = setExp
+
+-- | Convert a 'Set' to an 'OmegaSet'.
+toOmegaSet :: Set -> OmegaSet
+toOmegaSet = setOmegaSet
+
+upperBound :: Set -> Set
+upperBound s = useSetSet L.upperBound (setDim s) s
+
+lowerBound :: Set -> Set
+lowerBound s = useSetSet L.lowerBound (setDim s) s
+
+lowerBoundSatisfiable :: Set -> Bool
+lowerBoundSatisfiable = useSet L.lowerBoundSatisfiable
+
+upperBoundSatisfiable :: Set -> Bool
+upperBoundSatisfiable = useSet L.upperBoundSatisfiable
+
+obviousTautology :: Set -> Bool
+obviousTautology = useSet L.obviousTautology
+
+definiteTautology :: Set -> Bool
+definiteTautology = useSet L.definiteTautology
+
+exact :: Set -> Bool
+exact = useSet L.exact
+
+inexact :: Set -> Bool
+inexact = useSet L.inexact
+
+unknown :: Set -> Bool
+unknown = useSet L.unknown
+
+-- | Test whether two sets are equal.
+-- The sets must have the same dimension
+-- (@dimension s1 == dimension s2@), or an error will be raised.
+--
+-- The answer is precise if both relations are 'exact'.
+-- If either relation is inexact, this function returns @False@.
+equal :: Set -> Set -> Bool
+equal = useSet2 L.equal
+
+-- | Union of two sets.
+-- The sets must have the same dimension
+-- (@dimension s1 == dimension s2@), or an error will be raised.
+union :: Set -> Set -> Set
+union s1 s2 = useSet2Set L.union (setDim s1) s1 s2
+
+-- | Intersection of two sets.
+-- The sets must have the same dimension
+-- (@dimension s1 == dimension s2@), or an error will be raised.
+intersection :: Set -> Set -> Set
+intersection s1 s2 = useSet2Set L.intersection (setDim s1) s1 s2
+
+-- | Difference of two sets.
+-- The sets must have the same dimension
+-- (@dimension s1 == dimension s2@), or an error will be raised.
+difference :: Set -> Set -> Set
+difference s1 s2 = useSet2Set L.difference (setDim s1) s1 s2
+
+-- | Get the gist of a set, given some background truth.  The
+-- gist operator uses heuristics to simplify the set while
+-- retaining sufficient information to regenerate the original by
+-- re-introducing the background truth.  The sets must have the
+-- same dimension.
+--
+-- Given @x@ computed by
+--
+-- > x <- intersection given =<< gist effort r given
+--
+-- we have @x == r@.
+gist :: Effort -> Set -> Set -> Set
+gist effort s1 s2 = useSet2Set (L.gist effort) (setDim s1) s1 s2
+
+complement :: Set -> Set
+complement s = useSetSet L.complement (setDim s) s
+
+approximate :: Set -> Set
+approximate s = useSetSet L.approximate (setDim s) s
diff --git a/Data/Presburger/Omega/SetRel.hs b/Data/Presburger/Omega/SetRel.hs
new file mode 100644
--- /dev/null
+++ b/Data/Presburger/Omega/SetRel.hs
@@ -0,0 +1,104 @@
+
+-- | Internal routines used by both "Data.Presburger.Omega.Set" and
+-- "Data.Presburger.Omega.Rel"
+
+module Data.Presburger.Omega.SetRel where
+
+import Data.Presburger.Omega.LowLevel
+import Data.Presburger.Omega.Expr
+
+-- Make a lookup function for translating 'VarHandle's to 'Var's.
+-- The position of a handle determines what 'Var' it translates to. 
+makeLookupFunction :: [VarHandle] -> (VarHandle -> Var)
+makeLookupFunction lowLevelVars =
+    let expVars = takeFreeVariables (length lowLevelVars)
+        varLookupTable = zip lowLevelVars expVars
+
+        findVar v = case lookup v varLookupTable
+                    of Just v' -> v'
+                       Nothing -> error "Cannot find Omega variable"
+    in findVar
+
+-- Create an expression fom some low-level data.
+--
+-- The boolean parameter is true for equality constraints,
+-- false for inequality constraints.
+constraintToExpr :: Bool          -- ^ Is equality
+                 -> [VarHandle]   -- ^ Bound variables
+                 -> [Coefficient] -- ^ Terms
+                 -> Int           -- ^ Constant part
+                 -> BoolExpr      -- ^ Expression
+constraintToExpr isEquality boundVars terms constant =
+    let -- The existential variables are innermost
+        findVar = makeLookupFunction boundVars
+
+        -- Sum of all products and the constant term
+        sumTerm = sumOfProductsExpr constant $ map productTerm terms
+            where
+              productTerm (Coefficient v n) = (n, [findVar v])
+
+        -- Test whether is equal to zero/nonnegative
+        boolTerm = if isEquality
+                   then testExpr IsZero sumTerm
+                   else testExpr IsGEZ sumTerm
+    in boolTerm
+
+-- Get the set as a function.
+-- We pass list-building routines to the low-level queryDNFSet function.
+setToExpression :: OmegaSet -> IO (Int, BoolExp)
+setToExpression s = do
+  (setVars, conjuncts) <- queryDNFSet addEq [] addGeq [] addConjunct [] s
+  return (length setVars, wrapSimplifiedExpr $ disjExpr conjuncts)
+    where
+      -- Call constraintToExpr with the existential variables bound first,
+      -- then the set variables.
+
+      addEq setVars exVars terms constant =
+          (constraintToExpr True (exVars ++ setVars) terms constant :)
+
+      addGeq setVars exVars terms constant =
+          (constraintToExpr False (exVars ++ setVars) terms constant :)
+
+      addConjunct _ exVars eqs geqs =
+          wrapExistentialVars exVars eqs geqs
+
+-- Get a relation as a boolean expression.
+-- In the formula, we expect to see only the output variable whose index
+-- is given by 'index'.
+--
+-- In the result expression, the chosen output variable is bound innermost,
+--  A mysterious error will occur otherwise.
+
+-- We pass list-building routines to the low-level queryDNFRelation function.
+relToExpression :: OmegaRel -> IO (Int, Int, BoolExp)
+relToExpression s = do
+  (inVars, outVars, cs) <- queryDNFRelation addEq [] addGeq [] addConjunct [] s
+  return (length inVars, length outVars, wrapSimplifiedExpr $ disjExpr cs)
+    where
+      addEq inVars outVars exVars terms constant =
+          let vars = exVars ++ inVars ++ outVars
+          in (constraintToExpr True vars terms constant :)
+
+      addGeq inVars outVars exVars terms constant =
+          let vars = exVars ++ inVars ++ outVars
+          in (constraintToExpr False vars terms constant :)
+
+      addConjunct _ _ exVars eqs geqs =
+          wrapExistentialVars exVars eqs geqs
+
+      hasExistentialVars = error "relToExpression: cannot create expression"
+
+wrapExistentialVars exVars eqs geqs = (conjunct :)
+    where
+      conjunct =
+          -- Create a conjunction of constraints, with one quantifier for each
+          -- existential variable
+          iterateN existsExpr (length exVars) $ conjExpr (geqs ++ eqs)
+    
+
+-- Apply a function n times
+iterateN f n x = go n x
+    where go 0 x = x
+          go n x = go (n-1) (f x)
+
+
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,22 @@
+Copyright (c) 2009 Christopher Rodrigues
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+ * Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+ * Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Makefile.in b/Makefile.in
new file mode 100644
--- /dev/null
+++ b/Makefile.in
@@ -0,0 +1,21 @@
+
+CXX=@CXX@
+
+CPPFLAGS=@CPPFLAGS@
+CXXFLAGS=@CXXFLAGS@
+LDFLAGS=@LDFLAGS@
+
+LIBS=@LIBS@
+
+INCLUDES=src/C_omega.h
+
+.PHONY: all clean
+
+all : build/C_omega.o
+
+build :
+	mkdir $@
+
+build/C_omega.o : src/C_omega.cc build $(INCLUDES)
+	$(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $< -o $@
+
diff --git a/Omega.cabal b/Omega.cabal
new file mode 100644
--- /dev/null
+++ b/Omega.cabal
@@ -0,0 +1,41 @@
+Name:			Omega
+Version:		0.1.1
+Cabal-Version:		>= 1.2.3
+Build-Type:		Custom
+License:		BSD3
+License-File:		LICENSE
+Author:			Christopher Rodrigues
+Maintainer:		cirodrig@illinois.edu
+Stability:		Alpha
+Synopsis:		Operations on Presburger arithmetic formulae
+Description:
+        This package provides tools for manipulating sets and relations
+        whose members can be represented compactly as a Presburger
+        arithmetic formula.  The primary interface can be found
+        in "Data.Presburger.Omega.Set" and "Data.Presburger.Omega.Rel".
+
+        The Omega library
+        (<http://github.com/davewathaverford/the-omega-project>) must
+        be installed to build this package.
+Category:		Data
+Extra-Source-Files:
+	README
+	configure.ac
+	Makefile.in
+	src/C_omega.h
+	src/C_omega.cc
+Extra-Tmp-Files:	build/C_omega.o
+
+Library
+  Build-Depends:	base >= 3 && < 4, containers
+  Exposed-Modules:
+        Data.Presburger.Omega.Expr
+        Data.Presburger.Omega.LowLevel
+        Data.Presburger.Omega.Set
+        Data.Presburger.Omega.Rel
+  Other-Modules:
+        Data.Presburger.Omega.SetRel
+  Extensions:		GADTs ScopedTypeVariables
+  Build-Tools:		hsc2hs
+  Include-Dirs:		src
+  Extra-Libraries:	omega stdc++
diff --git a/README b/README
new file mode 100644
--- /dev/null
+++ b/README
@@ -0,0 +1,38 @@
+Omega -- Operations on Presburger arithmetic formulae
+
+BUILDING INSTRUCTIONS
+---------------------
+
+This is a Cabal package.  The typical build process is:
+
+	runhaskell Setup.hs configure
+	runhaskell Setup.hs build
+	runhaskell Setup.hs install
+
+This package requires the C++ Omega library to be installed
+(http://github.com/davewathaverford/the-omega-project).  Because this package
+contains C++ source code, Cabal will probably need help finding the required
+headers and libraries.
+
+You will probably need to provide the paths to the C++ include directory
+(contains STL headers such as "vector") and library directory (contains the
+C runtime library, called "libstdc++.so" on GNU Linux systems).  If the C++
+Omega library is not installed in a standard place, you will also need to
+provide paths to it.
+
+A configuration might look something like this:
+
+	runhaskell Setup.hs configure \
+		--extra-include-dirs=$(YOUR_CXX_INCLUDE_PATH) \
+		--extra-lib-dirs=$(YOUR_CXX_LIB_PATH) \
+		--extra-include-dirs=$(YOUR_OMEGA_PATH)/basic/include \
+		--extra-include-dirs=$(YOUR_OMEGA_PATH)/omega_lib/include \
+		--extra-lib-dirs=$(YOUR_OMEGA_PATH)/omega_lib/obj
+
+DOCUMENTATION
+-------------
+
+The C++ Omega library includes documentation of its exported interface.
+You may wish to look there if the Haddock documentation for a set operation
+or relation operation is lacking.
+
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,166 @@
+
+import Control.Applicative
+import Control.Monad
+import Data.Char
+import Data.Maybe
+import Distribution.PackageDescription
+import Distribution.Simple
+import Distribution.Simple.BuildPaths
+import Distribution.Simple.LocalBuildInfo
+import Distribution.Simple.Program
+import Distribution.Simple.Setup
+import Distribution.Simple.Utils
+import qualified Distribution.Verbosity as Verbosity
+import System.Cmd
+import System.Directory
+import System.Exit(ExitCode(..))
+import System.IO
+import System.FilePath((</>))
+import System.Process
+
+-- Mimic the && command of 'sh'
+(>&&>) :: IO ExitCode -> IO ExitCode -> IO ExitCode
+cmd1 >&&> cmd2 = do
+  rc <- cmd1
+  case rc of
+    ExitSuccess   -> cmd2
+    ExitFailure _ -> return rc
+
+-- We will call 'autoconf' and 'make'
+autoconfProgram = simpleProgram "autoconf"
+makeProgram = simpleProgram "make"
+
+-- Our single C++ source file is here
+cppSourceName = "src" </> "C_omega.cc"
+
+-- It becomes this object file
+cppObjectName = "build" </> "C_omega.o"
+
+-------------------------------------------------------------------------------
+-- Configuration
+
+configureOmega pkgDesc flags = do
+  -- Run Cabal configure
+  lbi <- confHook simpleUserHooks pkgDesc flags
+
+  let verb = fromFlagOrDefault Verbosity.normal $ configVerbosity flags
+      cfg = withPrograms lbi 
+
+      runAutoconf = do rawSystemProgramConf verb autoconfProgram cfg []
+                       return ExitSuccess
+      
+  -- Run autoconf configure
+  runAutoconf >&&> runConfigure lbi
+
+  return lbi
+
+    where
+      -- Run 'configure' with the extra arguments that were passed to
+      -- Setup.hs
+      runConfigure lbi = do
+        currentDir <- getCurrentDirectory
+
+        let opts = autoConfigureOptions lbi
+            configProgramName = currentDir </> "configure"
+
+        rawSystem configProgramName opts
+
+-- Configuration: extract options to pass to 'configure'
+autoConfigureOptions :: LocalBuildInfo -> [String]
+autoConfigureOptions localBuildInfo = [libdirs, includedirs]
+    where
+      libraryDescr = case library $ localPkgDescr localBuildInfo
+                     of Nothing -> error "Library description is missing"
+                        Just l -> l
+
+      buildinfo = libBuildInfo libraryDescr
+
+      -- Create a string "-L/usr/foo -L/usr/bar"
+      ldflagsString =
+          intercalate " " ["-L" ++ dir | dir <- extraLibDirs buildinfo]
+
+      libdirs = "LDFLAGS=" ++ ldflagsString
+
+      -- Create a string "-I/usr/foo -I/usr/bar"
+      cppflagsString =
+          intercalate " " ["-I" ++ dir | dir <- includeDirs buildinfo]
+
+      includedirs = "CPPFLAGS=" ++ cppflagsString
+
+-------------------------------------------------------------------------------
+-- Building
+
+buildOmega pkgDesc lbi userhooks flags = do
+  -- Do default build procedure for hs files
+  buildHook simpleUserHooks pkgDesc lbi userhooks flags
+
+  -- Get 'ar' program
+  let verb = fromFlagOrDefault Verbosity.normal $ buildVerbosity flags
+  (arPgm, _) <- requireProgram verb arProgram AnyVersion (withPrograms lbi)
+
+  -- Build the C++ source file
+  rawSystemProgramConf verb makeProgram (withPrograms lbi) ["all"]
+
+  -- Add the object file to libraries
+  let pkgId   = package $ localPkgDescr lbi
+
+  let addStaticObjectFile objName libName =
+          rawSystemProgram verb arPgm ["r", libName, objName]
+
+  when (withVanillaLib lbi) $
+       let libName = buildDir lbi </> mkLibName pkgId
+       in addStaticObjectFile cppObjectName libName
+
+  when (withProfLib lbi) $
+       let libName = buildDir lbi </> mkProfLibName pkgId
+       in addStaticObjectFile cppObjectName libName
+
+  when (withSharedLib lbi) $
+       die "Sorry, this package is not set up to build shared libraries"
+
+  return ()
+
+-------------------------------------------------------------------------------
+-- Cleaning
+
+cleanOmega pkgDesc mlbi userhooks flags = do
+  let verb = fromFlagOrDefault Verbosity.normal $ cleanVerbosity flags
+
+  -- Clean extra files if we don't need to save configuration
+  -- (Other temp files are automatically cleaned)
+  unless (fromFlag $ cleanSaveConf flags) $ do
+    lenientRemoveFiles configFiles
+    lenientRemoveDirectory "autom4te.cache"
+
+  -- Do default clean procedure
+  cleanHook simpleUserHooks pkgDesc mlbi userhooks flags
+
+    where
+      -- Attempt to remove a file, ignoring errors
+      lenientRemoveFile f =
+          removeFile f `catch` \_ -> return ()
+
+      lenientRemoveFiles = mapM_ lenientRemoveFile
+
+      -- Attempt to remove a directory and its contents
+      -- (one level of recursion only), ignoring errors
+      lenientRemoveDirectory f = do
+        b <- doesDirectoryExist f
+        when b $ do lenientRemoveFiles =<< getDirectoryContents f
+                    removeDirectory f `catch` \_ -> return ()
+
+      -- Extra files produced by configuration
+      configFiles = ["configure", "config.log", "config.status", "Makefile"]
+
+-------------------------------------------------------------------------------
+-- Hooks
+
+hooks =
+    simpleUserHooks
+    { hookedPrograms = [arProgram, autoconfProgram, makeProgram]
+    , confHook = configureOmega
+    , buildHook = buildOmega
+    , cleanHook = cleanOmega
+    }
+
+main = defaultMainWithHooks hooks
diff --git a/configure.ac b/configure.ac
new file mode 100644
--- /dev/null
+++ b/configure.ac
@@ -0,0 +1,39 @@
+
+# Initialization
+AC_INIT(Omega, 0.1)
+AC_LANG(C++)
+
+
+# Check for programs
+AC_PROG_CXX
+
+# Check the omega library
+AC_MSG_CHECKING([whether we can include basic/bool.h])
+AC_COMPILE_IFELSE(
+	[AC_LANG_SOURCE([[#include <basic/bool.h>
+		]])],
+	[AC_MSG_RESULT([ok])],
+	[AC_MSG_FAILURE([cannot include basic/bool.h])])
+
+AC_MSG_CHECKING([whether we can include omega.h])
+AC_COMPILE_IFELSE(
+	[AC_LANG_SOURCE([[#include <omega.h>
+		]])],
+	[AC_MSG_RESULT([ok])],
+	[AC_MSG_FAILURE([cannot include omega.h])])
+
+AC_MSG_CHECKING([whether we can link with omega library])
+{
+ LIBS="${LIBS} -lomega"
+ AC_LINK_IFELSE(
+	[AC_LANG_PROGRAM(
+		[[#include <omega.h>]],
+		[[omega::Relation::Null();]])],
+	[AC_MSG_RESULT([yes])],
+	[AC_MSG_FAILURE([cannot link with the omega library])])
+}
+
+
+# Output
+AC_CONFIG_FILES([Makefile])
+AC_OUTPUT
diff --git a/src/C_omega.cc b/src/C_omega.cc
new file mode 100644
--- /dev/null
+++ b/src/C_omega.cc
@@ -0,0 +1,535 @@
+
+#include <omega.h>
+#include <string.h>
+
+#include "C_omega.h"
+
+extern "C"
+Relation *hsw_new_relation(int n_input, int n_output)
+{
+  return new Relation(n_input, n_output);
+}
+
+extern "C"
+Relation *hsw_new_set(int n)
+{
+  return new Relation(n);
+}
+
+extern "C"
+void hsw_free_relation(Relation *rel)
+{
+  delete rel;
+}
+
+extern "C"
+char *hsw_relation_show(Relation *rel)
+{
+  return strdup((const char *)rel->print_with_subs_to_string());
+}
+
+extern "C"
+int hsw_num_input_vars(Relation *rel)
+{
+  return rel->n_inp();
+}
+
+extern "C"
+int hsw_num_output_vars(Relation *rel)
+{
+  return rel->n_out();
+}
+
+extern "C"
+int hsw_num_set_vars(Relation *rel)
+{
+  return rel->n_set();
+}
+
+extern "C"
+Var_Decl *hsw_input_var(Relation *rel, int n)
+{
+  return rel->input_var(n);
+}
+
+extern "C"
+Var_Decl *hsw_output_var(Relation *rel, int n)
+{
+  return rel->output_var(n);
+}
+extern "C"
+Var_Decl *hsw_set_var(Relation *rel, int n)
+{
+  return rel->set_var(n);
+}
+
+extern "C"
+int hsw_is_lower_bound_satisfiable(Relation *rel)
+{
+  return rel->is_lower_bound_satisfiable();
+}
+
+extern "C"
+int hsw_is_upper_bound_satisfiable(Relation *rel)
+{
+  return rel->is_upper_bound_satisfiable();
+}
+
+extern "C"
+int hsw_is_obvious_tautology(Relation *rel)
+{
+  return rel->is_obvious_tautology();
+}
+extern "C"
+int hsw_is_definite_tautology(Relation *rel)
+{
+  return rel->is_tautology();
+}
+
+extern "C"
+int hsw_is_exact(Relation *rel)
+{
+  return rel->is_exact();
+}
+
+extern "C"
+int hsw_is_inexact(Relation *rel)
+{
+  return rel->is_inexact();
+}
+
+extern "C"
+int hsw_is_unknown(Relation *rel)
+{
+  return rel->is_unknown();
+}
+
+extern "C"
+Relation *hsw_upper_bound(Relation *rel)
+{
+  return new Relation(Upper_Bound(copy(*rel)));
+}
+
+extern "C"
+Relation *hsw_lower_bound(Relation *rel)
+{
+  return new Relation(Lower_Bound(copy(*rel)));
+}
+
+extern "C"
+int hsw_equal(Relation *r, Relation *s)
+{
+  /*   r == s
+   * iff
+   *    r `intersection` not s == False
+   * && r `union` not s        == True
+   */
+  Relation com_s = Complement(copy(*s));
+
+  /* If intersection is satisfiable, unequal */
+  if (Intersection(copy(*r), copy(com_s)).is_upper_bound_satisfiable())
+    return 0;
+
+  /* If union is tautology, equal; else unequal */
+  return Union(copy(*r), com_s).is_tautology();
+}
+
+extern "C"
+Relation *hsw_union(Relation *r, Relation *s)
+{
+  return new Relation(Union(copy(*r), copy(*s)));
+}
+
+extern "C"
+Relation *hsw_intersection(Relation *r, Relation *s)
+{
+  return new Relation(Intersection(copy(*r), copy(*s)));
+}
+
+extern "C"
+Relation *hsw_composition(Relation *r, Relation *s)
+{
+  return new Relation(Composition(copy(*r), copy(*s)));
+}
+
+extern "C"
+Relation *hsw_restrict_domain(Relation *r, Relation *s)
+{
+  return new Relation(Restrict_Domain(copy(*r), copy(*s)));
+}
+
+extern "C"
+Relation *hsw_restrict_range(Relation *r, Relation *s)
+{
+  return new Relation(Restrict_Range(copy(*r), copy(*s)));
+}
+
+extern "C"
+Relation *hsw_difference(Relation *r, Relation *s)
+{
+  return new Relation(Difference(copy(*r), copy(*s)));
+}
+
+extern "C"
+Relation *hsw_cross_product(Relation *r, Relation *s)
+{
+  return new Relation(Cross_Product(copy(*r), copy(*s)));
+}
+
+extern "C"
+Relation *hsw_gist(Relation *r, Relation *s, int effort)
+{
+  return new Relation(Gist(copy(*r), copy(*s), effort));
+}
+
+extern "C"
+Relation *hsw_transitive_closure(Relation *rel)
+{
+  return new Relation(TransitiveClosure(copy(*rel)));
+}
+
+extern "C"
+Relation *hsw_domain(Relation *rel)
+{
+  return new Relation(Domain(copy(*rel)));
+}
+
+extern "C"
+Relation *hsw_range(Relation *rel)
+{
+  return new Relation(Range(copy(*rel)));
+}
+
+extern "C"
+Relation *hsw_inverse(Relation *rel)
+{
+  return new Relation(Inverse(copy(*rel)));
+}
+
+extern "C"
+Relation *hsw_complement(Relation *rel)
+{
+  return new Relation(Complement(copy(*rel)));
+}
+
+extern "C"
+Relation *hsw_deltas(Relation *rel)
+{
+  return new Relation(Deltas(copy(*rel)));
+}
+
+extern "C"
+Relation *hsw_approximate(Relation *rel)
+{
+  return new Relation(Approximate(copy(*rel)));
+}
+
+extern "C"
+F_And *hsw_relation_add_and(Relation *rel)
+{
+  return rel->add_and();
+}
+
+extern "C"
+Formula *hsw_relation_add_or(Relation *rel)
+{
+  return rel->add_or();
+}
+
+extern "C"
+Formula *hsw_relation_add_not(Relation *rel)
+{
+  return rel->add_not();
+}
+
+extern "C"
+F_Declaration *hsw_relation_add_forall(Relation *rel)
+{
+  return rel->add_forall();
+}
+
+extern "C"
+F_Declaration *hsw_relation_add_exists(Relation *rel)
+{
+  return rel->add_exists();
+}
+
+extern "C"
+void hsw_relation_finalize(Relation *rel)
+{
+  rel->finalize();
+}
+
+extern "C"
+Var_Decl *hsw_declaration_declare(F_Declaration *rel)
+{
+  return rel->declare();
+}
+
+extern "C"
+F_And *hsw_formula_to_and(Formula *rel)
+{
+  F_And *and_formula = dynamic_cast<F_And *>(rel);
+
+  /* If the parameter is already an 'and', return it */
+  if (and_formula) return and_formula;
+
+  /* Otherwise add an 'and' */
+  return rel->add_and();
+}
+
+extern "C"
+F_And *hsw_formula_add_and(Formula *rel)
+{
+  return rel->add_and();
+}
+
+extern "C"
+Formula *hsw_formula_add_or(Formula *rel)
+{
+  return rel->add_or();
+}
+
+extern "C"
+Formula *hsw_formula_add_not(Formula *rel)
+{
+  return rel->add_not();
+}
+
+extern "C"
+F_Declaration *hsw_formula_add_forall(Formula *rel)
+{
+  return rel->add_forall();
+}
+
+extern "C"
+F_Declaration *hsw_formula_add_exists(Formula *rel)
+{
+  return rel->add_exists();
+}
+
+extern "C"
+void hsw_formula_finalize(Formula *rel)
+{
+  rel->finalize();
+}
+
+/* hsw_add_constraint creates an equality or inequality constraint,
+ * fills in the coefficients for each variable, and fills in the
+ * constant term. */
+extern "C"
+void hsw_add_constraint(F_And *formula,
+		    int is_eq,
+		    int num_vars,
+		    int *coefficients,
+		    Var_Decl **vars,
+		    int constant)
+{
+  Constraint_Handle *hdl = is_eq
+    ? (Constraint_Handle *)new EQ_Handle(formula->add_EQ())
+    : (Constraint_Handle *)new GEQ_Handle(formula->add_GEQ());
+
+  /* Update each coefficient in the array */
+  for (; num_vars; num_vars--)
+    {
+      int index = num_vars - 1;
+      hdl->update_coef(vars[index], coefficients[index]);
+    }
+
+  /* Update the constant part of the constraint */
+  hdl->update_const(constant);
+
+  hdl->finalize();
+  free(hdl);
+}
+
+/* These are all for inspecting a DNF formula */
+
+extern "C"
+DNF_Iterator *hsw_query_dnf(Relation *rel)
+{
+  return new DNF_Iterator(rel->query_DNF());
+}
+
+extern "C"
+Conjunct *hsw_dnf_iterator_next(DNF_Iterator *iter)
+{
+  if (!iter->live()) return NULL;
+
+  Conjunct *c = **iter;
+  ++*iter;
+  return c;
+}
+
+extern "C"
+void hsw_dnf_iterator_free(DNF_Iterator *iter)
+{
+  delete iter;
+}
+
+/* Use to iterate over the tuple of the variables that are used in the
+ * conjunct.  The variables obtained should not be freed. */
+extern "C"
+struct Tuple_Iter *hsw_get_conjunct_variables(Conjunct *conj)
+{
+  Tuple_Iterator<void *> *ti =
+    reinterpret_cast<Tuple_Iterator<void *> *>
+    (new Tuple_Iterator<Variable_ID>(*conj->variables()));
+  return (struct Tuple_Iter *)ti;
+}
+
+extern "C"
+void *
+hsw_tuple_iterator_next(struct Tuple_Iter *iter)
+{
+  Tuple_Iterator<void *> *ti = (Tuple_Iterator<void *> *)iter;
+
+  if (!ti->live()) return NULL;	// Exhausted?
+
+  void *ret = (void *)**ti;
+  ++*ti;
+  return ret;
+}
+
+extern "C"
+void
+hsw_tuple_iterator_free(struct Tuple_Iter *iter)
+{
+  delete (Tuple_Iterator<void *> *)iter;
+}
+
+/* Use to iterate over the EQ constraints in a conjunct.  The constraints
+ * obtained should be freed once you're done with them. */
+extern "C"
+struct EQ_Iterator *
+hsw_get_eqs(Conjunct *conj)
+{
+  return new EQ_Iterator(conj->EQs());
+}
+
+extern "C"
+struct EQ_Handle *
+hsw_eqs_next(struct EQ_Iterator *g)
+{
+  if (!g->live()) return NULL;	// Exhausted?
+
+  EQ_Handle *hdl = new EQ_Handle(**g);
+  ++*g;
+  return hdl;
+}
+
+extern "C"
+void
+hsw_eqs_free(struct EQ_Iterator *g)
+{
+  delete g;
+}
+
+extern "C"
+void
+hsw_eq_handle_free(struct EQ_Handle *hdl)
+{
+  delete hdl;
+}
+
+/* Use to iterate over the GEQ constraints in a conjunct.  Works like
+ * hsw_get_eqs. */
+extern "C"
+struct GEQ_Iterator *hsw_get_geqs(Conjunct *conj)
+{
+  return new GEQ_Iterator(conj->GEQs());
+}
+
+extern "C"
+struct GEQ_Handle *
+hsw_geqs_next(struct GEQ_Iterator *g)
+{
+  if (!g->live()) return NULL;	// Exhausted?
+
+  GEQ_Handle *hdl = new GEQ_Handle(**g);
+  ++*g;
+  return hdl;
+}
+
+extern "C"
+void
+hsw_geqs_free(struct GEQ_Iterator *g)
+{
+  delete g;
+}
+
+extern "C"
+void
+hsw_geq_handle_free(struct GEQ_Handle *hdl)
+{
+  delete hdl;
+}
+
+extern "C"
+coefficient_t
+hsw_constraint_get_const(struct Constraint_Handle_ *hdl)
+{
+  return ((struct Constraint_Handle *)hdl)->get_const();
+}
+
+extern "C"
+Constr_Vars_Iter *
+hsw_constraint_get_coefficients(struct Constraint_Handle_ *hdl)
+{
+  return new Constr_Vars_Iter(*(Constraint_Handle *)hdl);  
+}
+
+extern "C"
+int
+hsw_constr_vars_next(Variable_Info_struct *out, Constr_Vars_Iter *iter)
+{
+  if (!iter->live()) return 0;
+
+  Variable_Info info(**iter);
+  ++*iter;
+
+  out->var = info.var;
+  out->coef = info.coef;
+
+  return 1;
+}
+
+extern "C"
+void
+hsw_constr_vars_free(Constr_Vars_Iter *iter)
+{
+  delete iter;
+}
+
+/* For debugging */
+
+extern "C"
+void
+hsw_debug_print_eq(struct EQ_Handle *hdl)
+{
+  String s(hdl->print_to_string());
+  puts(s);
+}
+
+extern "C"
+void
+hsw_debug_print_geq(struct GEQ_Handle *hdl)
+{
+  String s(hdl->print_to_string());
+  puts(s);
+}
+
+#if 0 /* Not used? */
+
+/* Find an array element equal to v.  Return the element index,
+ * or -1 if no element matches. */
+static int
+find_variable_index(Var_Decl *v, int num_vars, Var_Decl **vars)
+{
+  int n;
+  for (n = 0; n < num_vars; n++) {
+    if (v == vars[n]) return n;
+  }
+  return -1;
+}
+#endif
diff --git a/src/C_omega.h b/src/C_omega.h
new file mode 100644
--- /dev/null
+++ b/src/C_omega.h
@@ -0,0 +1,115 @@
+
+#ifndef C_OMEGA_H
+#define C_OMEGA_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* This is a copy of 'coef_t'.  Can't use the original because it's in
+ * a C++ header file. */
+typedef long long coefficient_t;
+
+/* This is a copy of struct Variable_Info.  Can't use the original because
+ * it's in a C++ header file. */
+typedef struct Variable_Info_struct {
+  struct Var_Decl *var;
+  coefficient_t    coef;
+} Variable_Info_struct;
+
+struct Relation *hsw_new_relation(int n_input, int n_output);
+struct Relation *hsw_new_set(int n);
+void hsw_free_relation(struct Relation *rel);
+char *hsw_relation_show(struct Relation *rel);
+int hsw_num_input_vars(struct Relation *rel);
+int hsw_num_output_vars(struct Relation *rel);
+int hsw_num_set_vars(struct Relation *rel);
+struct Var_Decl *hsw_input_var(struct Relation *rel, int n);
+struct Var_Decl *hsw_output_var(struct Relation *rel, int n);
+struct Var_Decl *hsw_set_var(struct Relation *rel, int n);
+int hsw_is_lower_bound_satisfiable(struct Relation *rel);
+int hsw_is_upper_bound_satisfiable(struct Relation *rel);
+int hsw_is_obvious_tautology(struct Relation *rel);
+int hsw_is_definite_tautology(struct Relation *rel);
+int hsw_is_exact(struct Relation *rel);
+int hsw_is_inexact(struct Relation *rel);
+int hsw_is_unknown(struct Relation *rel);
+struct Relation *hsw_upper_bound(struct Relation *);
+struct Relation *hsw_lower_bound(struct Relation *);
+int hsw_equal(struct Relation *, struct Relation *);
+struct Relation *hsw_union(struct Relation *, struct Relation *);
+struct Relation *hsw_intersection(struct Relation *, struct Relation *);
+struct Relation *hsw_composition(struct Relation *, struct Relation *);
+struct Relation *hsw_restrict_domain(struct Relation *, struct Relation *);
+struct Relation *hsw_restrict_range(struct Relation *, struct Relation *);
+struct Relation *hsw_difference(struct Relation *, struct Relation *);
+struct Relation *hsw_cross_product(struct Relation *, struct Relation *);
+struct Relation *hsw_gist(struct Relation *, struct Relation *, int);
+struct Relation *hsw_transitive_closure(struct Relation *);
+struct Relation *hsw_domain(struct Relation *);
+struct Relation *hsw_range(struct Relation *);
+struct Relation *hsw_inverse(struct Relation *);
+struct Relation *hsw_complement(struct Relation *);
+struct Relation *hsw_deltas(struct Relation *);
+struct Relation *hsw_approximate(struct Relation *);
+
+struct F_And *hsw_relation_add_and(struct Relation *rel);
+struct Formula *hsw_relation_add_or(struct Relation *rel);
+struct Formula *hsw_relation_add_not(struct Relation *rel);
+struct F_Declaration *hsw_relation_add_forall(struct Relation *rel);
+struct F_Declaration *hsw_relation_add_exists(struct Relation *rel);
+void hsw_relation_finalize(struct Relation *rel);
+
+struct F_And *hsw_formula_add_and(struct Formula *rel);
+struct Formula *hsw_formula_add_or(struct Formula *rel);
+struct Formula *hsw_formula_add_not(struct Formula *rel);
+struct F_Declaration *hsw_formula_add_forall(struct Formula *rel);
+struct F_Declaration *hsw_formula_add_exists(struct Formula *rel);
+void hsw_formula_finalize(struct Formula *rel);
+
+struct Var_Decl *hsw_declaration_declare(struct F_Declaration *rel);
+
+struct F_And *hsw_formula_to_and(struct Formula *rel);
+
+void hsw_add_constraint(struct F_And *formula,
+		    int is_eq,
+		    int num_vars,
+		    int *coefficients,
+		    struct Var_Decl **vars,
+		    int constant);
+
+struct DNF_Iterator *hsw_query_dnf(struct Relation *rel);
+struct Conjunct *hsw_dnf_iterator_next(struct DNF_Iterator *iter);
+void hsw_dnf_iterator_free(struct DNF_Iterator *iter);
+
+struct Tuple_Iter *hsw_get_conjunct_variables(struct Conjunct *conj);
+void *hsw_tuple_iterator_next(struct Tuple_Iter *iter);
+void hsw_tuple_iterator_free(struct Tuple_Iter *iter);
+
+struct EQ_Iterator *hsw_get_eqs(struct Conjunct *conj);
+struct EQ_Handle *hsw_eqs_next(struct EQ_Iterator *g);
+void hsw_eqs_free(struct EQ_Iterator *g);
+void hsw_eq_handle_free(struct EQ_Handle *hdl);
+
+struct GEQ_Iterator *hsw_get_geqs(struct Conjunct *conj);
+struct GEQ_Handle *hsw_geqs_next(struct GEQ_Iterator *g);
+void hsw_geqs_free(struct GEQ_Iterator *g);
+void hsw_geq_handle_free(struct GEQ_Handle *hdl);
+
+struct Constraint_Handle_;	/* Use a different name to get rid of C++ warning */
+coefficient_t hsw_constraint_get_const(struct Constraint_Handle_ *hdl);
+struct Constr_Vars_Iter *hsw_constraint_get_coefficients(struct Constraint_Handle_ *hdl);
+int hsw_constr_vars_next(Variable_Info_struct *out, struct Constr_Vars_Iter *iter);
+void hsw_constr_vars_free(struct Constr_Vars_Iter *iter);
+
+
+
+void hsw_debug_print_eq(struct EQ_Handle *hdl);
+void hsw_debug_print_geq(struct GEQ_Handle *hdl);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif