diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,1 @@
+
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,23 @@
+Copyright (c) 2023, QBayLogic B.V.
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+2. 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 OWNER 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/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,11 @@
+# SoP Satisfier
+
+~~Kind of SMT solver with only Non-linear Natural Arithmetic.~~
+
+Library with an SMTlib like interface to declare and assert SoP (kind of) expressions.
+
+Interface:
+- `declare` - to declare expression to the state
+- `assert` - to assert that expression holds in the state
+- `unify` - to get a list of expressions that need to hold for the given expression to hold (only equalities are supported)
+- `range` - to get a range (a pair of minimum and maximum possible values) of expression
diff --git a/sop-satisfier.cabal b/sop-satisfier.cabal
new file mode 100644
--- /dev/null
+++ b/sop-satisfier.cabal
@@ -0,0 +1,53 @@
+cabal-version:      3.0
+name:               sop-satisfier
+version:            0.3.4.5
+synopsis:           Check satisfiability of expressions on natural numbers
+description:
+    Expression satisfier on natural numbers.
+    .
+    It can reason about expressions contatining
+    addition and multiplication.
+    It also provides limited support of exponentiations and subtraction.
+license:            BSD-2-Clause
+license-file:       LICENSE
+author:             Aleksandr Pokatilov <pokatilov0802@gmail.com>
+maintainer:         QBayLogic B.V. <devops@qbaylogic.com>
+build-type:         Simple
+category:           Solver, Symbolic Arithmetic
+extra-doc-files:    CHANGELOG.md
+                    README.md
+tested-with:        GHC == 8.8.4, GHC == 8.10.7, GHC == 9.0.2, GHC == 9.2.8,
+                    GHC == 9.4.8, GHC == 9.6.6, GHC == 9.8.4, GHC == 9.10.1,
+                    GHC == 9.12.1
+
+source-repository head
+    type: git
+    location: https://github.com/clash-lang/sop-satisfier
+
+common warnings
+    ghc-options: -Wall
+
+library
+    import:           warnings
+    exposed-modules:  SoPSat.Satisfier,
+                      SoPSat.SoP,
+                      SoPSat.Internal.Unify,
+                      SoPSat.Internal.Range,
+                      SoPSat.Internal.SoP,
+                      SoPSat.Internal.NewtonsMethod,
+                      SoPSat.Internal.SolverMonad
+    build-depends:    base >=4.13 && <5,
+                      containers >=0.6.2.1 && <0.8,
+                      transformers >=0.5 && <0.7
+    hs-source-dirs:   src
+    default-language: Haskell2010
+
+test-suite system-tests
+    type:             exitcode-stdio-1.0
+    main-is:          SystemTests.hs
+    build-depends:    base >=4.13 && <5,
+                      sop-satisfier,
+                      tasty ^>= 1.5,
+                      tasty-hunit ^>=0.9
+    hs-source-dirs:   tests
+    default-language: Haskell2010
diff --git a/src/SoPSat/Internal/NewtonsMethod.hs b/src/SoPSat/Internal/NewtonsMethod.hs
new file mode 100644
--- /dev/null
+++ b/src/SoPSat/Internal/NewtonsMethod.hs
@@ -0,0 +1,157 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module SoPSat.Internal.NewtonsMethod
+where
+
+import Data.Map (Map)
+import qualified Data.Map as M
+import Data.Maybe (fromJust)
+
+import SoPSat.Internal.SoP (
+  Atom (..),
+  Product (..),
+  SoP (..),
+  Symbol (..),
+ )
+import SoPSat.SoP (atoms)
+
+-- | Evaluates SoP given atom bindings
+evalSoP ::
+  (Ord f, Ord c, Floating n) =>
+  -- | Expression to evaluate
+  SoP f c ->
+  -- | Bindings from atoms to values
+  Map (Atom f c) n ->
+  -- | Evaluation result
+  n
+evalSoP (S []) _ = 0
+evalSoP (S ps) binds = sum $ map (`evalProduct` binds) ps
+
+{- | Evaluates product given atom bindings
+
+Used by @evalSoP@
+-}
+evalProduct ::
+  (Ord f, Ord c, Floating n) =>
+  -- | Product to evalute
+  Product f c ->
+  -- | Atom bindings
+  Map (Atom f c) n ->
+  -- | Evaluation results
+  n
+evalProduct (P ss) binds = product $ map (`evalSymbol` binds) ss
+
+{- | Evaluates symbol given atom bindings
+
+Used by @evalProduct@
+-}
+evalSymbol ::
+  (Ord f, Ord c, Floating n) =>
+  -- | Symbol to evaluate
+  Symbol f c ->
+  -- | Atom bindings
+  Map (Atom f c) n ->
+  -- | Evaluation result
+  n
+evalSymbol (I i) _ = fromInteger i
+evalSymbol (A a) binds = f $ M.lookup a binds
+ where
+  f (Just n) = n
+  f Nothing = 0
+evalSymbol (E b p) binds = exp (evalProduct p binds * log (evalSoP b binds))
+
+{- | Analitically computes derivative of an expression
+with respect to an atom
+
+Returns function similar to @evalSoP@
+-}
+derivative ::
+  (Ord f, Ord c, Floating n) =>
+  -- | Expression to take a derivative of
+  SoP f c ->
+  -- | Atom to take a derivetive with respect to
+  Atom f c ->
+  -- | Function from bindings, representing point,
+  -- to value of the derivative at that point
+  (Map (Atom f c) n -> n)
+derivative sop symb = \binds -> sum $ d <*> [binds]
+ where
+  d = map (`derivativeProduct` symb) $ unS sop
+
+{- | Analitically computes derivative of a product
+with respect to an atom
+
+Used by @derivative@
+-}
+derivativeProduct ::
+  (Ord f, Ord c, Floating n) =>
+  -- | Product to take a derivative of
+  Product f c ->
+  -- | Atom to take a derivative with respect to
+  Atom f c ->
+  -- | Function from bindings to a value
+  (Map (Atom f c) n -> n)
+derivativeProduct (P []) _ = const 0
+derivativeProduct (P (s : ss)) symb = \binds ->
+  derivativeSymbol s symb binds * evalProduct ps binds
+    + evalSymbol s binds * derivativeProduct ps symb binds
+ where
+  ps = P ss
+
+{- | Analitically computes derivative of a symbol
+with respect to an atom
+
+Used by @derivativeProduct@
+-}
+derivativeSymbol ::
+  (Ord f, Ord c, Floating n) =>
+  -- | Symbol to take a derivate of
+  Symbol f c ->
+  -- | Atom to take a derivate with respect to
+  Atom f c ->
+  -- | Function from bindings to a value
+  (Map (Atom f c) n -> n)
+derivativeSymbol (I _) _ = const 0
+derivativeSymbol (A a) atom
+  | a == atom = const 1
+  | otherwise = const 0
+derivativeSymbol e@(E b p) atom = \binds ->
+  expExpr binds
+    * ( derivative b atom binds
+          * evalProduct p binds
+          / evalSoP b binds
+          + logExpr binds
+            * derivativeProduct p atom binds
+      )
+ where
+  expExpr = evalSymbol e
+  logExpr = log . evalSoP b
+
+-- | Finds if an expression can be equal to zero
+newtonMethod ::
+  forall f c n.
+  (Ord f, Ord c, Ord n, Floating n) =>
+  -- | Expression to check
+  SoP f c ->
+  -- | @Right binds@ - Atom bindings when expression is equal to zero
+  --   @Left binds@ - Last checked bindings
+  Either (Map (Atom f c) n) (Map (Atom f c) n)
+newtonMethod sop = go init_guess steps
+ where
+  consts = atoms sop
+  derivs = M.fromSet (derivative sop) consts
+  init_guess = M.fromSet (const 10) consts
+  steps = 40
+
+  go :: Map (Atom f c) n -> Word -> Either (Map (Atom f c) n) (Map (Atom f c) n)
+  go guess 0 = Left guess
+  go guess n
+    | val <= 0.1 = Right guess
+    | otherwise =
+        let
+          new_guess = foldl (\binds (c, x) -> M.insert c (x - val / dsdc c) binds) guess $ M.toList guess
+         in
+          go new_guess (n - 1)
+   where
+    val = evalSoP sop guess
+    dsdc c = fromJust (M.lookup c derivs) guess
diff --git a/src/SoPSat/Internal/Range.hs b/src/SoPSat/Internal/Range.hs
new file mode 100644
--- /dev/null
+++ b/src/SoPSat/Internal/Range.hs
@@ -0,0 +1,97 @@
+module SoPSat.Internal.Range (
+  Range (..),
+  Bound (..),
+  boundSoP,
+  rangeAdd,
+  rangeMul,
+  rangeExp,
+)
+where
+
+import SoPSat.SoP
+
+data Bound f c
+  = Bound (SoP f c)
+  | Inf
+  deriving (Eq, Show)
+
+boundSoP :: Bound f c -> Maybe (SoP f c)
+boundSoP (Bound s) = Just s
+boundSoP Inf = Nothing
+
+data Range f c
+  = Range
+  { lower :: Bound f c
+  , upper :: Bound f c
+  }
+  deriving (Eq, Show)
+
+boundAdd :: (Ord f, Ord c) => Bound f c -> Bound f c -> Bound f c
+boundAdd Inf _ = Inf
+boundAdd _ Inf = Inf
+boundAdd (Bound a) (Bound b) = Bound (a |+| b)
+
+boundMul :: (Ord f, Ord c) => Bound f c -> Bound f c -> Bound f c
+boundMul Inf _ = Inf
+boundMul _ Inf = Inf
+boundMul (Bound a) (Bound b) = Bound (a |*| b)
+
+boundExp :: (Ord f, Ord c) => Bound f c -> Bound f c -> Bound f c
+boundExp Inf _ = Inf
+boundExp _ Inf = Inf
+boundExp (Bound a) (Bound b) = Bound (a |^| b)
+
+rangeAdd :: (Ord f, Ord c) => Range f c -> Range f c -> Maybe (Range f c)
+-- Subtraction of unbounded functions
+rangeAdd (Range _ Inf) (Range Inf _) = Nothing
+rangeAdd (Range Inf _) (Range _ Inf) = Nothing
+rangeAdd (Range low1 up1) (Range low2 up2) =
+  Just $
+    Range (boundAdd low1 low2) (boundAdd up1 up2)
+
+rangeMul :: (Ord f, Ord c) => Range f c -> Range f c -> Maybe (Range f c)
+-- Multiplication of unbounded functions
+rangeMul (Range Inf Inf) _ = Nothing
+rangeMul _ (Range Inf Inf) = Nothing
+-- Multiplication with infinitely increasing/decresing functions
+rangeMul (Range low1 Inf) (Range low2 _) =
+  Just $
+    Range (boundMul low1 low2) Inf
+rangeMul (Range low1 _) (Range low2 Inf) =
+  Just $
+    Range (boundMul low1 low2) Inf
+rangeMul (Range Inf up1) (Range low2 _) =
+  Just $
+    Range Inf (boundMul up1 low2)
+rangeMul (Range low1 _) (Range Inf up2) =
+  Just $
+    Range Inf (boundMul up2 low1)
+rangeMul (Range low1 up1) (Range low2 up2) =
+  Just $
+    Range (boundMul low1 low2) (boundMul up1 up2)
+
+-- rangeMul (Range low1 up1) (Range low2 up2)
+-- --   | sopSign low1 == sopSign low2
+-- --   = Range
+-- rangeMul (Range low1 up1) (Range low2 up2) = let
+--     low1Sign = sopSign =<< boundSoP low1
+--     low2Sign = sopSign =<< boundSoP low2
+--   in case (low1Sign,low2Sign) of
+--        (Just Positive, Just Positive) -> Just $
+--          Range (boundMul low1 low2) (boundMul up1 up2)
+--        (Just Negative, Just Positive) -> Just $
+--          Range (boundMul low1 up2) (boundMul up1 low2)
+--        (Just Positive, Just Negative) -> Just $
+--          Range (boundMul up1 low2) (boundMul low1 up2)
+--        (Just Negative, Just Negative) -> Just $
+--          Range (boundMul
+-- rangeMul _ _ = Nothing
+
+rangeExp :: (Ord f, Ord c) => Range f c -> Range f c -> Maybe (Range f c)
+rangeExp (Range Inf _) (Range Inf _) = Nothing
+rangeExp (Range _ up1) (Range Inf up2) =
+  Just $
+    Range (Bound (int 0)) (boundExp up1 up2)
+rangeExp (Range low1 up1) (Range low2 up2) =
+  Just $
+    Range (boundExp low1 low2) (boundExp up1 up2)
diff --git a/src/SoPSat/Internal/SoP.hs b/src/SoPSat/Internal/SoP.hs
new file mode 100644
--- /dev/null
+++ b/src/SoPSat/Internal/SoP.hs
@@ -0,0 +1,209 @@
+module SoPSat.Internal.SoP
+where
+
+-- External
+import Data.Either (partitionEithers)
+import Data.List (intercalate, sort)
+
+{- | Atomic part of a @SoP@
+like constants and unknown functions
+-}
+data Atom f c
+  = -- | Constant
+    C c
+  | -- | Unknown function
+    F f [SoP f c]
+  deriving (Eq, Ord)
+
+instance (Show f, Show c) => Show (Atom f c) where
+  show (C c) = show c
+  show (F f args) = show f ++ "(" ++ intercalate ", " (map show args) ++ ")"
+
+{- | The most basic part used during reasoning:
+- Numbers
+- Atoms
+- Exponents
+-}
+data Symbol f c
+  = -- | Number in an expression
+    I Integer
+  | -- | Atom in an expression
+    A (Atom f c)
+  | -- | Exponentiation
+    E (SoP f c) (Product f c)
+  deriving (Eq, Ord)
+
+instance (Show f, Show c) => Show (Symbol f c) where
+  show (E s p) = show s ++ "^" ++ show p
+  show (I i) = show i
+  show (A a) = show a
+
+-- | Product of symbols
+newtype Product f c = P {unP :: [Symbol f c]}
+  deriving (Eq)
+
+instance (Ord f, Ord c) => Ord (Product f c) where
+  compare (P [x]) (P [y]) = compare x y
+  compare (P [_]) (P (_ : _)) = LT
+  compare (P (_ : _)) (P [_]) = GT
+  compare (P xs) (P ys) = compare xs ys
+
+instance (Show f, Show c) => Show (Product f c) where
+  show (P [s]) = show s
+  show (P ss) = "(" ++ intercalate " * " (map show ss) ++ ")"
+
+-- | Sum of Products
+newtype SoP f c = S {unS :: [Product f c]}
+  deriving (Ord)
+
+instance (Eq f, Eq c) => Eq (SoP f c) where
+  (S []) == (S [P [I 0]]) = True
+  (S [P [I 0]]) == (S []) = True
+  (S ps1) == (S ps2) = ps1 == ps2
+
+instance (Show f, Show c) => Show (SoP f c) where
+  show (S [p]) = show p
+  show (S ps) = "(" ++ intercalate " + " (map show ps) ++ ")"
+
+mergeWith :: (a -> a -> Either a a) -> [a] -> [a]
+mergeWith _ [] = []
+mergeWith op (f : fs) = case partitionEithers $ map (`op` f) fs of
+  ([], _) -> f : mergeWith op fs
+  (updated, untouched) -> mergeWith op (updated ++ untouched)
+
+reduceExp :: (Ord f, Ord c) => Symbol f c -> Symbol f c
+reduceExp (E _ (P [I 0])) = I 1
+reduceExp (E (S [P [I 0]]) _) = I 0
+reduceExp (E (S [P [I i]]) (P [I j]))
+  | j >= 0 = I (i ^ j)
+reduceExp (E (S [P [E k i]]) j) =
+  case normaliseExp k (S [e]) of
+    (S [P [s]]) -> s
+    _ -> E k e
+ where
+  e = P . sort . map reduceExp $ mergeWith mergeS (unP i ++ unP j)
+reduceExp s = s
+
+mergeS ::
+  (Ord f, Ord c) =>
+  Symbol f c ->
+  Symbol f c ->
+  Either (Symbol f c) (Symbol f c)
+mergeS (I i) (I j) = Left (I (i * j))
+mergeS (I 1) r = Left r
+mergeS l (I 1) = Left l
+mergeS (I 0) _ = Left (I 0)
+mergeS _ (I 0) = Left (I 0)
+-- x * x^4 ==> x^5
+mergeS s (E (S [P [s']]) (P [I i]))
+  | s == s' =
+      Left (E (S [P [s']]) (P [I (i + 1)]))
+-- x^4 * x ==> x^5
+mergeS (E (S [P [s']]) (P [I i])) s
+  | s == s' =
+      Left (E (S [P [s']]) (P [I (i + 1)]))
+-- 4^x * 2^x ==> 8^x
+mergeS (E (S [P [I i]]) p) (E (S [P [I j]]) p')
+  | p == p' =
+      Left (E (S [P [I (i * j)]]) p)
+-- y*y ==> y^2
+mergeS l r
+  | l == r =
+      case normaliseExp (S [P [l]]) (S [P [I 2]]) of
+        (S [P [e]]) -> Left e
+        _ -> Right l
+-- x^y * x^(-y) ==> 1
+mergeS (E s1 (P p1)) (E s2 (P (I i : p2)))
+  | i == (-1)
+  , s1 == s2
+  , p1 == p2 =
+      Left (I 1)
+-- x^(-y) * x^y ==> 1
+mergeS (E s1 (P (I i : p1))) (E s2 (P p2))
+  | i == (-1)
+  , s1 == s2
+  , p1 == p2 =
+      Left (I 1)
+mergeS l _ = Right l
+
+mergeP ::
+  (Eq f, Eq c) =>
+  Product f c ->
+  Product f c ->
+  Either (Product f c) (Product f c)
+-- 2xy + 3xy ==> 5xy
+mergeP (P ((I i) : is)) (P ((I j) : js))
+  | is == js = Left . P $ I (i + j) : is
+-- 2xy + xy  ==> 3xy
+mergeP (P ((I i) : is)) (P js)
+  | is == js = Left . P $ I (i + 1) : is
+-- xy + 2xy  ==> 3xy
+mergeP (P is) (P ((I j) : js))
+  | is == js = Left . P $ I (j + 1) : is
+-- xy + xy ==> 2xy
+mergeP (P is) (P js)
+  | is == js = Left . P $ I 2 : is
+  | otherwise = Right $ P is
+
+normaliseExp :: (Ord f, Ord c) => SoP f c -> SoP f c -> SoP f c
+-- b^1 ==> b
+normaliseExp b (S [P [I 1]]) = b
+-- x^(2xy) ==> x^(2xy)
+normaliseExp b@(S [P [A _]]) (S [e]) = S [P [E b e]]
+-- 2^(y^2) ==> 4^y
+normaliseExp b@(S [P [_]]) (S [e@(P [_])]) = S [P [reduceExp (E b e)]]
+-- (x + 2)^2 ==> x^2 + 4xy + 4
+normaliseExp b (S [P [I i]])
+  | i > 0 =
+      foldr1 mergeSoPMul (replicate (fromInteger i) b)
+-- (x + 2)^(2x) ==> (x^2 + 4xy + 4)^x
+normaliseExp b (S [P (e@(I i) : es)])
+  | i >= 0 =
+      -- Without the "| i >= 0" guard, normaliseExp can loop with itself
+      -- for exponentials such as: 2^(n-k)
+      normaliseExp (normaliseExp b (S [P [e]])) (S [P es])
+-- (x + 2)^(xy) ==> (x+2)^(xy)
+normaliseExp b (S [e]) = S [P [reduceExp (E b e)]]
+-- (x + 2)^(y + 2) ==> 4x(2 + x)^y + 4(2 + x)^y + (2 + x)^yx^2
+normaliseExp b (S e) = foldr1 mergeSoPMul (map (normaliseExp b . S . (: [])) e)
+
+zeroP :: Product f c -> Bool
+zeroP (P ((I 0) : _)) = True
+zeroP _ = False
+
+mkNonEmpty :: (Ord f, Ord c) => SoP f c -> SoP f c
+mkNonEmpty (S []) = S [P [I 0]]
+mkNonEmpty s = s
+
+simplifySoP :: (Ord f, Ord c) => SoP f c -> SoP f c
+simplifySoP = repeatF go
+ where
+  go =
+    mkNonEmpty
+      . S
+      . sort
+      . filter (not . zeroP)
+      . mergeWith mergeP
+      . map (P . sort . map reduceExp . mergeWith mergeS . unP)
+      . unS
+
+  repeatF f x =
+    let x' = f x
+     in if x' == x
+          then x
+          else repeatF f x'
+{-# INLINEABLE simplifySoP #-}
+
+mergeSoPAdd :: (Ord f, Ord c) => SoP f c -> SoP f c -> SoP f c
+mergeSoPAdd (S ps1) (S ps2) = simplifySoP $ S (ps1 ++ ps2)
+
+mergeSoPMul :: (Ord f, Ord c) => SoP f c -> SoP f c -> SoP f c
+mergeSoPMul (S ps1) (S ps2) =
+  simplifySoP . S $
+    concatMap (zipWith (\p1 p2 -> P (unP p1 ++ unP p2)) ps1 . repeat) ps2
+
+mergeSoPSub :: (Ord f, Ord c) => SoP f c -> SoP f c -> SoP f c
+mergeSoPSub a b = mergeSoPAdd a (mergeSoPMul (S [P [I (-1)]]) b)
+
+mergeSoPDiv :: (Ord f, Ord c) => SoP f c -> SoP f c -> (SoP f c, SoP f c)
+mergeSoPDiv (S _ps1) (S _ps2) = undefined
diff --git a/src/SoPSat/Internal/SolverMonad.hs b/src/SoPSat/Internal/SolverMonad.hs
new file mode 100644
--- /dev/null
+++ b/src/SoPSat/Internal/SolverMonad.hs
@@ -0,0 +1,102 @@
+{-# LANGUAGE RecordWildCards #-}
+
+module SoPSat.Internal.SolverMonad
+where
+
+import Control.Monad.Trans.State.Strict (
+  StateT (..),
+  evalStateT,
+  get,
+  gets,
+  put,
+ )
+
+import Data.Map (Map)
+import qualified Data.Map as M
+
+import SoPSat.Internal.Range
+import SoPSat.Internal.SoP (
+  Product (..),
+  SoP (..),
+  Symbol (..),
+ )
+import SoPSat.Internal.Unify
+import SoPSat.SoP
+import qualified SoPSat.SoP as SoP
+
+data (Ord f, Ord c) => State f c
+  = State (Map (Atom f c) (Range f c)) [Unifier f c]
+  deriving (Show)
+
+instance (Ord f, Ord c) => Semigroup (State f c) where
+  (State r1 u1) <> (State r2 u2) = State (M.union r1 r2) (u1 ++ u2)
+
+instance (Ord f, Ord c) => Monoid (State f c) where
+  mempty = State M.empty []
+
+-- TODO: Change Maybe to some MonadError for better error indication
+type SolverState f c = StateT (State f c) Maybe
+
+maybeFail :: (MonadFail m) => Maybe a -> m a
+maybeFail (Just a) = return a
+maybeFail Nothing = fail ""
+
+getRanges :: (Ord f, Ord c) => SolverState f c (Map (Atom f c) (Range f c))
+getRanges = gets (\(State rangeS _) -> rangeS)
+
+getRange :: (Ord f, Ord c) => Atom f c -> SolverState f c (Range f c)
+getRange c = maybeFail . M.lookup c =<< getRanges
+
+getRangeSymbol :: (Ord f, Ord c) => Symbol f c -> SolverState f c (Range f c)
+getRangeSymbol (E b p) = maybeFail =<< rangeExp <$> getRangeSoP b <*> getRangeProduct p
+getRangeSymbol i@(I _) = return range
+ where
+  bound = Bound (toSoP i)
+  range = Range bound bound
+getRangeSymbol (A a) = getRange a
+
+getRangeProduct :: (Ord f, Ord c) => Product f c -> SolverState f c (Range f c)
+getRangeProduct p = maybeFail . foldl rm oneRange =<< mapM getRangeSymbol (unP p)
+ where
+  one = Bound $ SoP.int 1
+  oneRange = Just (Range one one)
+  rm Nothing _ = Nothing
+  rm (Just a) b = rangeMul a b
+
+getRangeSoP :: (Ord f, Ord c) => SoP f c -> SolverState f c (Range f c)
+getRangeSoP s = maybeFail . foldl ra zeroRange =<< mapM getRangeProduct (unS s)
+ where
+  zero = Bound $ SoP.int 0
+  zeroRange = Just (Range zero zero)
+  ra Nothing _ = Nothing
+  ra (Just a) b = rangeAdd a b
+
+putRange :: (Ord f, Ord c) => Atom f c -> Range f c -> SolverState f c ()
+putRange symb range@Range{..} = do
+  -- Anti-symmetry: 5 <= x ^ x <= 5 => x = 5
+  case (lower == upper, upper) of
+    (True, Bound bound) -> putUnifiers [Subst symb (toSoP bound)]
+    _ -> return ()
+  (State rangeS unifyS) <- get
+  let rangeSn = M.insert symb range rangeS
+  put (State rangeSn unifyS)
+
+getUnifiers :: (Ord f, Ord c) => SolverState f c [Unifier f c]
+getUnifiers = gets (\(State _ unifyS) -> unifyS)
+
+putUnifiers :: (Ord f, Ord c) => [Unifier f c] -> SolverState f c ()
+putUnifiers us = do
+  (State rangeS unifyS) <- get
+  put (State rangeS (substsSubst us unifyS ++ us))
+
+-- | Puts a state to use during computations
+withState :: (Ord f, Ord c) => State f c -> SolverState f c ()
+withState = put
+
+-- | Runs computation returning result and resulting state
+runStatements :: (Ord f, Ord c) => SolverState f c a -> Maybe (a, State f c)
+runStatements stmts = runStateT stmts mempty
+
+-- | Similar to @runStatements@ but does not return final state
+evalStatements :: (Ord f, Ord c) => SolverState f c a -> Maybe a
+evalStatements stmts = evalStateT stmts mempty
diff --git a/src/SoPSat/Internal/Unify.hs b/src/SoPSat/Internal/Unify.hs
new file mode 100644
--- /dev/null
+++ b/src/SoPSat/Internal/Unify.hs
@@ -0,0 +1,225 @@
+{-# LANGUAGE RecordWildCards #-}
+
+module SoPSat.Internal.Unify
+where
+
+import Data.Function (on)
+import Data.List (find, intersect, nub, partition, (\\))
+
+import SoPSat.Internal.SoP (
+  Atom (..),
+  Product (..),
+  SoP (..),
+  Symbol (..),
+ )
+import SoPSat.SoP (
+  toSoP,
+  (|*|),
+  (|+|),
+  (|-|),
+  (|^|),
+ )
+import qualified SoPSat.SoP as SoP
+
+data Unifier f c
+  = Subst
+  { sConst :: Atom f c
+  , sSoP :: SoP f c
+  }
+  deriving (Eq, Show)
+
+substsSoP :: (Ord f, Ord c) => [Unifier f c] -> SoP f c -> SoP f c
+substsSoP [] u = u
+substsSoP ((Subst{..}) : s) u = substsSoP s (substSoP sConst sSoP u)
+
+substSoP :: (Ord f, Ord c) => Atom f c -> SoP f c -> SoP f c -> SoP f c
+substSoP cons subs = foldr1 (|+|) . map (substProduct cons subs) . unS
+
+substProduct :: (Ord f, Ord c) => Atom f c -> SoP f c -> Product f c -> SoP f c
+substProduct cons subs = foldr1 (|*|) . map (substSymbol cons subs) . unP
+
+substSymbol :: (Ord f, Ord c) => Atom f c -> SoP f c -> Symbol f c -> SoP f c
+substSymbol _ _ s@(I _) = toSoP s
+substSymbol cons subst s@(A a)
+  | cons == a = subst
+  | otherwise = S [P [s]]
+substSymbol cons subst (E b p) = substSoP cons subst b |^| substProduct cons subst p
+
+substsSubst :: (Ord f, Ord c) => [Unifier f c] -> [Unifier f c] -> [Unifier f c]
+substsSubst s = map subst
+ where
+  subst sub@(Subst{..}) = sub{sSoP = substsSoP s sSoP}
+
+unifiers :: (Ord f, Ord c) => SoP f c -> SoP f c -> [Unifier f c]
+unifiers (S [P [A a]]) (S []) = [Subst a (S [P [I 0]])]
+unifiers (S []) (S [P [A a]]) = [Subst a (S [P [I 0]])]
+unifiers (S [P [I _]]) (S [P [I _]]) = []
+-- (z ^ a) ~ (z ^ b) ==> [a := b]
+unifiers (S [P [E s1 p1]]) (S [P [E s2 p2]])
+  | s1 == s2 = unifiers (toSoP p1) (toSoP p2)
+-- (2*e ^ d) ~ (2*e*a*c) ==> [a*c := 2*e ^ (d-1)]
+unifiers (S [P [E (S [P s1]) p1]]) (S [P p2])
+  | all (`elem` p2) s1 =
+      let base = s1 `intersect` p2
+          diff = p2 \\ s1
+       in unifiers (S [P diff]) (S [P [E (S [P base]) (P [I (-1)]), E (S [P base]) p1]])
+unifiers (S [P p2]) (S [P [E (S [P s1]) p1]])
+  | all (`elem` p2) s1 =
+      let base = s1 `intersect` p2
+          diff = p2 \\ s1
+       in unifiers (S [P diff]) (S [P [E (S [P base]) (P [I (-1)]), E (S [P base]) p1]])
+-- (i ^ a) ~ j ==> [a := round (logBase i j)], when `i` and `j` are integers,
+unifiers (S [P [E (S [P [I i]]) p]]) (S [P [I j]]) =
+  case integerLogBase i j of
+    Just k -> unifiers (S [p]) (S [P [I k]])
+    Nothing -> []
+unifiers (S [P [I j]]) (S [P [E (S [P [I i]]) p]]) =
+  case integerLogBase i j of
+    Just k -> unifiers (S [p]) (S [P [I k]])
+    Nothing -> []
+-- x ^ i = j => [x := root i j]
+unifiers (S [P [E s (P [I p])]]) (S [P [I j]]) =
+  case integerRt p j of
+    Just k -> unifiers s (S [P [I k]])
+    Nothing -> []
+unifiers (S [P [I j]]) (S [P [E s (P [I p])]]) =
+  case integerRt p j of
+    Just k -> unifiers s (S [P [I k]])
+    Nothing -> []
+-- a^d * a^e ~ a^c ==> [c := d + e]
+unifiers (S [P [E s1 p1]]) (S [p2]) = case collectBases p2 of
+  Just (b : bs, ps)
+    | all (== s1) (b : bs) ->
+        unifiers (S [p1]) (S ps)
+  _ -> []
+unifiers (S [p2]) (S [P [E s1 p1]]) = case collectBases p2 of
+  Just (b : bs, ps)
+    | all (== s1) (b : bs) ->
+        unifiers (S ps) (S [p1])
+  _ -> []
+-- (i * a) ~ (j * b) ==> [a := (div j i) * b]
+-- Where 'a' and 'b' are variables, 'i' and 'j' are integer literals, and j `mod` i == 0
+unifiers (S [P ((I i) : ps)]) (S [P ((I j) : ps1)])
+  | (Just k) <- safeDiv j i = unifiers (S [P ps]) (SoP.int k |*| S [P ps1])
+  | (Just k) <- safeDiv i j = unifiers (SoP.int k |*| S [P ps]) (S [P ps1])
+  | otherwise = []
+-- (2*a) ~ (2*b) ==> [a := b]
+-- unifiers (S [P (p:ps1)]) (S [P (p':ps2)])
+--     | p == p'   = unifiers' ct (S [P ps1]) (S [P ps2])
+--     | otherwise = []
+unifiers (S [P ps1@(_ : _ : _)]) (S [P ps2])
+  | null psx = []
+  | otherwise = unifiers (S [P ps1'']) (S [P ps2''])
+ where
+  ps1' = ps1 \\ psx
+  ps2' = ps2 \\ psx
+  ps1''
+    | null ps1' = [I 1]
+    | otherwise = ps1'
+  ps2''
+    | null ps2' = [I 1]
+    | otherwise = ps2'
+  psx = ps1 `intersect` ps2
+unifiers (S [P ps1]) (S [P ps2@(_ : _ : _)])
+  | null psx = []
+  | otherwise = unifiers (S [P ps1'']) (S [P ps2''])
+ where
+  ps1' = ps1 \\ psx
+  ps2' = ps2 \\ psx
+  ps1''
+    | null ps1' = [I 1]
+    | otherwise = ps1'
+  ps2''
+    | null ps2' = [I 1]
+    | otherwise = ps2'
+  psx = ps1 `intersect` ps2
+unifiers (S [P [A a]]) s = [Subst a s]
+unifiers s (S [P [A a]]) = [Subst a s]
+-- (2 + a) ~ 5 ==> [a := 3]
+unifiers (S ((P [I i]) : ps1)) (S ((P [I j]) : ps2))
+  | i < j = unifiers (S ps1) (S (P [I (j - i)] : ps2))
+  | i > j = unifiers (S (P [I (i - j)] : ps1)) (S ps2)
+-- (a + c) ~ (b + c) ==> [a := b]
+unifiers s1@(S ps1) s2@(S ps2) = case splitSoP s1 s2 of
+  (s1', s2')
+    | s1' /= s1 || s2' /= s2 ->
+        unifiers s1' s2'
+  _
+    | null psx
+    , length ps1 == length ps2 ->
+        case nub (concat (zipWith (\x y -> unifiers (S [x]) (S [y])) ps1 ps2)) of
+          [] -> unifiers' s1 s2
+          [k] -> [k]
+          _ -> []
+    | null psx ->
+        unifiers' s1 s2
+  _ -> unifiers' (S ps1'') (S ps2'')
+ where
+  ps1' = ps1 \\ psx
+  ps2' = ps2 \\ psx
+  ps1''
+    | null ps1' = [P [I 0]]
+    | otherwise = ps1'
+  ps2''
+    | null ps2' = [P [I 0]]
+    | otherwise = ps2'
+  psx = ps1 `intersect` ps2
+
+unifiers' :: (Ord f, Ord c) => SoP f c -> SoP f c -> [Unifier f c]
+unifiers' (S [P [I i], P [A a]]) s2 =
+  [Subst a (s2 |+| S [P [I (negate i)]])]
+unifiers' s1 (S [P [I i], P [A a]]) =
+  [Subst a (s1 |+| S [P [I (negate i)]])]
+unifiers' _ _ = []
+
+splitSoP :: (Ord f, Ord c) => SoP f c -> SoP f c -> (SoP f c, SoP f c)
+splitSoP u v = (lhs, rhs)
+ where
+  reduced = v |-| u
+  (lhs', rhs') = partition neg (unS reduced)
+  lhs
+    | null lhs' = SoP.int 0
+    | otherwise = ((|*|) `on` S) lhs' [P [I (-1)]]
+  rhs
+    | null rhs' = SoP.int 0
+    | otherwise = S rhs'
+
+  neg (P ((I i) : _)) = i < 0
+  neg _ = False
+
+collectBases :: Product f c -> Maybe ([SoP f c], [Product f c])
+collectBases = fmap unzip . traverse go . unP
+ where
+  go (E s1 p1) = Just (s1, p1)
+  go _ = Nothing
+
+safeDiv :: Integer -> Integer -> Maybe Integer
+safeDiv i j
+  | j == 0 = Just 0
+  | otherwise = case divMod i j of
+      (k, 0) -> Just k
+      _ -> Nothing
+
+integerLogBase :: Integer -> Integer -> Maybe Integer
+integerLogBase x y
+  | x > 1 && y > 0 =
+      let z1 = integerLogBase' x y
+          z2 = integerLogBase' x (y - 1)
+       in if z1 == z2
+            then Nothing
+            else Just z1
+integerLogBase _ _ = Nothing
+
+integerLogBase' :: Integer -> Integer -> Integer
+integerLogBase' b m = snd (go b)
+ where
+  go :: Integer -> (Integer, Integer)
+  go pw | m < pw = (m, 0)
+  go pw = case go (pw ^ (2 :: Int)) of
+    (q, e) | q < pw -> (q, 2 * e)
+    (q, e) -> (q `quot` pw, 2 * e + 1)
+
+-- Naive implementation of exact integer root
+integerRt :: Integer -> Integer -> Maybe Integer
+integerRt 1 y = Just y
+integerRt x y = find ((== y) . (^ x)) [1 .. y]
diff --git a/src/SoPSat/Satisfier.hs b/src/SoPSat/Satisfier.hs
new file mode 100644
--- /dev/null
+++ b/src/SoPSat/Satisfier.hs
@@ -0,0 +1,443 @@
+{-# LANGUAGE RecordWildCards #-}
+
+module SoPSat.Satisfier (
+  -- * State
+  SolverState,
+
+  -- * State manipulation
+  declare,
+  assert,
+  unify,
+
+  -- * State information
+  range,
+  ranges,
+
+  -- * State execution
+  withState,
+  runStatements,
+  evalStatements,
+
+  -- * Expressions
+  evalSoP,
+)
+where
+
+import Control.Applicative ((<|>))
+import Control.Arrow (second)
+import Control.Monad (unless, when, (>=>))
+
+import Data.Map (Map)
+import qualified Data.Map as M
+import Data.Maybe (isNothing)
+
+import SoPSat.Internal.NewtonsMethod
+import SoPSat.Internal.Range
+import SoPSat.Internal.SoP (
+  Atom (..),
+  Product (..),
+  SoP (..),
+  Symbol (..),
+ )
+import SoPSat.Internal.SolverMonad
+import SoPSat.Internal.Unify
+import SoPSat.SoP
+
+parts :: [a] -> [[a]]
+parts [] = []
+parts (x : xs) = xs : map (x :) (parts xs)
+
+{- | Declares atom in the state
+ignores constants and only declare function arguments
+-}
+declareAtom :: (Ord f, Ord c) => Atom f c -> SolverState f c Bool
+declareAtom (C _) = return True
+declareAtom (F _ args) = and <$> mapM declareSoP args
+
+{- | Declares symbol in the state with the default interval
+If symbol exists preserves the old interval
+-}
+declareSymbol :: (Ord f, Ord c) => Symbol f c -> SolverState f c Bool
+declareSymbol (I _) = return True
+declareSymbol (A a) = do
+  existing <- getRanges
+  when (isNothing (M.lookup a existing)) (putRange a rangeNatural)
+  declareAtom a
+ where
+  rangeNatural = Range (Bound (int 0)) Inf
+declareSymbol (E b p) = (&&) <$> declareSoP b <*> declareProduct p
+
+-- | Similar to @declareSoP@ but for @Product@
+declareProduct :: (Ord f, Ord c) => Product f c -> SolverState f c Bool
+declareProduct = fmap and . mapM declareSymbol . unP
+
+{- | Declare SoP in the state with default values
+Creates range for free-variables
+-}
+declareSoP :: (Ord f, Ord c) => SoP f c -> SolverState f c Bool
+declareSoP s@(S ps)
+  | left /= int 0 =
+      (&&) <$> (and <$> mapM declareProduct ps) <*> assert (SoPE left right LeR)
+  | otherwise = and <$> mapM declareProduct ps
+ where
+  (left, right) = splitSoP (int 0) s
+
+{- | Declare expression to state, returns normalised expression
+
+Common for @declare@, @assert@, and @unify@
+-}
+declareToState :: (Ord f, Ord c) => SoPE f c -> SolverState f c (SoPE f c)
+declareToState SoPE{..} = do
+  r1 <- declareSoP lhs
+  r2 <- declareSoP rhs
+  us <- getUnifiers
+  let
+    lhs' = substsSoP us lhs
+    rhs' = substsSoP us rhs
+  unless (r1 && r2) (fail "")
+  return (SoPE lhs' rhs' op)
+
+{- | Declare equality of two expressions
+Adds new unifiers to the state
+-}
+declareEq ::
+  (Ord f, Ord c) =>
+  -- | First expression
+  SoP f c ->
+  -- | Second expression
+  SoP f c ->
+  -- | Similar to @declare@ but handles only equalities
+  SolverState f c Bool
+declareEq u v =
+  do
+    (Range low1 up1) <- getRangeSoP u
+    (Range low2 up2) <- getRangeSoP v
+    lowRes <- boundComp low1 low2
+    upRes <- boundComp up1 up2
+
+    -- Declaration and assertions of expression is done on the whole domain
+    -- if two expressions are equal, their domains will intersect
+    --
+    -- g(x) in [1,5] and forall x  g(x) = f(x) then f(x) in [1,5]
+    lowerUpdate <-
+      case (lowRes, low1, low2) of
+        (True, _, Bound lowB2) -> propagateInEqSoP u GeR lowB2
+        (False, Bound lowB1, _) -> propagateInEqSoP v GeR lowB1
+        (_, _, _) -> return True
+
+    upperUpdate <-
+      case (upRes, up1, up2) of
+        (True, _, Bound upB2) -> propagateInEqSoP u LeR upB2
+        (False, Bound upB1, _) -> propagateInEqSoP v LeR upB1
+        (_, _, _) -> return True
+
+    declareEq' u v
+    return (lowerUpdate && upperUpdate)
+ where
+  boundComp Inf _ = return False
+  boundComp _ Inf = return True
+  boundComp (Bound a) (Bound b) = assert (SoPE a b LeR)
+
+declareEq' :: (Ord f, Ord c) => SoP f c -> SoP f c -> SolverState f c ()
+declareEq' (S [P [A a]]) v = putUnifiers [Subst a v]
+declareEq' u (S [P [A a]]) = putUnifiers [Subst a u]
+declareEq' u v = putUnifiers $ unifiers u v
+
+-- | Updates interval information for a symbol
+propagateInEqSymbol ::
+  (Ord f, Ord c) =>
+  -- | Updated symbol
+  Symbol f c ->
+  -- | Relationship between the symbol and target
+  OrdRel ->
+  -- | Target Boundary
+  SoP f c ->
+  -- | Similat to @declareInEq@
+  SolverState f c Bool
+propagateInEqSymbol (I _) _ _ =
+  return True -- No need to update numbers
+propagateInEqSymbol (A a) rel bound = do
+  (Range low up) <- getRange a
+  -- New bound is less/greater than the old one
+  -- The check is done before propagation
+  -- This assumption is potentially wrong
+  case rel of
+    LeR ->
+      putRange a (Range low rangeBound)
+    GeR ->
+      putRange a (Range rangeBound up)
+    EqR -> error "propagateInEqSymbol:EqR: unreachable"
+  return True
+ where
+  rangeBound = Bound bound
+propagateInEqSymbol (E b (P [I i])) rel (S [P [I j]])
+  | (Just p) <- integerRt i j =
+      propagateInEqSoP b rel (int p)
+propagateInEqSymbol (E (S [P [I i]]) p) rel (S [P [I j]])
+  | (Just e) <- integerLogBase i j =
+      propagateInEqProduct p rel (int e)
+propagateInEqSymbol _ _ _ = fail ""
+
+-- | Propagates interval information down the Product
+propagateInEqProduct ::
+  (Ord f, Ord c) =>
+  -- | Updates expression
+  Product f c ->
+  -- | Relationship between the expression and target
+  OrdRel ->
+  -- | Target boundary
+  SoP f c ->
+  -- | Similar to @declareInEq@
+  SolverState f c Bool
+propagateInEqProduct (P [symb]) rel target_bound = propagateInEqSymbol symb rel target_bound
+propagateInEqProduct (P ss) rel target_bound =
+  and <$> mapM (uncurry propagate) (zipWith (curry (second P)) ss (parts ss))
+ where
+  -- a <= x * y => a/y <= x and a/x <= y
+  -- Currently simply propagating the bound further
+  -- a <= x * y => a <= x and a <= y
+  propagate symb _prod =
+    propagateInEqSymbol
+      symb
+      rel
+      target_bound
+
+-- (target_bound |/| prod)
+
+-- | Propagates interval information down the SoP
+propagateInEqSoP ::
+  (Ord f, Ord c) =>
+  -- | Updated expression
+  SoP f c ->
+  -- | Relationship between the expression and target
+  OrdRel ->
+  -- | Target boundary
+  SoP f c ->
+  -- | Similar to @declareInEq@
+  SolverState f c Bool
+propagateInEqSoP (S [P [symb]]) rel target_bound = propagateInEqSymbol symb rel target_bound
+propagateInEqSoP (S ps) rel target_bound =
+  and <$> mapM (uncurry propagate) (zipWith (curry (second S)) ps (parts ps))
+ where
+  -- a <= x + y => a - y <= x and a - x <= y
+  propagate prod sm =
+    propagateInEqProduct
+      prod
+      rel
+      (target_bound |-| sm)
+
+{- | Declare inequality of two expressions
+Updates interval information in the state
+-}
+declareInEq ::
+  (Ord f, Ord c) =>
+  -- | Relationship between expressions
+  OrdRel ->
+  -- | Left-hand side expression
+  SoP f c ->
+  -- | Right-hand side expression
+  SoP f c ->
+  -- | Similar to @declare@ but handles only inequalities
+  SolverState f c Bool
+declareInEq EqR u v = declareEq u v >> return True
+declareInEq op u v =
+  let
+    (u', v') = splitSoP u v
+   in
+    -- If inequality holds with current interval information
+    -- then no need to update it
+    do
+      res <- assert (SoPE u' v' op)
+      if res
+        then return True
+        else case op of
+          LeR -> do
+            a1 <- propagateInEqSoP u' LeR v'
+            a2 <- propagateInEqSoP v' GeR u'
+            return (a1 && a2)
+          GeR -> do
+            a1 <- propagateInEqSoP u' GeR v'
+            a2 <- propagateInEqSoP v' LeR u'
+            return (a1 && a2)
+
+-- | Declare expression to the state
+declare ::
+  (Ord f, Ord c) =>
+  -- | Expression to declare
+  SoPE f c ->
+  -- | - True - if expression was declared
+  --   - False - if expression contradicts current state
+  --
+  -- State will become @Nothing@ if it cannot reason about these kind of expressions
+  SolverState f c Bool
+declare =
+  declareToState >=> \SoPE{..} ->
+    case op of
+      EqR -> declareEq lhs rhs
+      _ -> declareInEq op lhs rhs
+
+-- | Assert that two expressions are equal using unifiers from the state
+assertEq ::
+  (Ord f, Ord c) =>
+  -- | Left-hand side expression
+  SoP f c ->
+  -- | Right-hand size expression
+  SoP f c ->
+  -- | Similar to assert but only checks for equality @lhs = rhs@
+  SolverState f c Bool
+assertEq lhs rhs = return (lhs == rhs)
+
+-- | Assert using only ranges stores in the state
+assertRange ::
+  (Ord f, Ord c) =>
+  -- | Left-hand side expression
+  SoP f c ->
+  -- | Right-hand size expression
+  SoP f c ->
+  -- | Similar to @assert@ but uses only intervals from the state to check @lhs <= rhs@
+  SolverState f c Bool
+assertRange lhs rhs = uncurry assertRange' $ splitSoP lhs rhs
+
+assertRange' :: (Ord f, Ord c) => SoP f c -> SoP f c -> SolverState f c Bool
+assertRange' (S [P [I i]]) (S [P [I j]]) = return (i <= j)
+assertRange' lhs rhs = do
+  (Range _ up1) <- getRangeSoP lhs
+  (Range low2 up2) <- getRangeSoP rhs
+  -- If both sides increase infinitely, fail to use Newton's method
+  -- Information about rate of growth is required
+  -- to check inequality on the whole domain
+  if up1 == up2 && up2 == Inf
+    then fail ""
+    else case (up1, low2) of
+      (Inf, _) -> return False
+      (_, Inf) -> return False
+      (Bound ub1, Bound lb2) ->
+        -- Orders of recursive checks matters
+        -- @runLemma2@ in the tests loops indefinitely
+        -- possibly other test cases too
+        do
+          r1 <-
+            if ub1 /= lhs
+              then assert (SoPE ub1 rhs LeR)
+              else return False
+          r2 <-
+            if lb2 /= rhs
+              then assert (SoPE lhs lb2 LeR)
+              else return False
+          return (r1 || r2)
+
+-- | Assert using only Newton's method
+assertNewton ::
+  (Ord f, Ord c) =>
+  -- | Left-hand side expression
+  SoP f c ->
+  -- | Right-hand side expression
+  SoP f c ->
+  -- | Similar to @assert@ but uses only Newton's method to check @lhs <= rhs@
+  SolverState f c Bool
+assertNewton lhs rhs =
+  let
+    expr = rhs |-| lhs |+| int 1
+   in
+    checkExpr expr
+ where
+  -- hasFunction :: (Ord f, Ord c) => SoP f c -> Bool
+  -- hasFunction = any isFunction . atoms
+
+  checkExpr :: (Ord f, Ord c) => SoP f c -> SolverState f c Bool
+  checkExpr expr
+    | (Right binds) <- newtonMethod expr =
+        not <$> checkBinds binds
+    | otherwise =
+        return True
+
+  checkBinds :: (Ord f, Ord c) => Map (Atom f c) Double -> SolverState f c Bool
+  checkBinds binds = and <$> mapM (uncurry (checkBind binds)) (M.toList binds)
+
+  checkBind ::
+    (Ord f, Ord c, Ord n, Floating n) =>
+    Map (Atom f c) n -> Atom f c -> n -> SolverState f c Bool
+  checkBind binds c v = do
+    (Range left right) <- getRange c
+    return (checkLeft binds v left && checkRight binds v right)
+
+  checkLeft ::
+    (Ord f, Ord c, Ord n, Floating n) => Map (Atom f c) n -> n -> Bound f c -> Bool
+  checkLeft _ _ Inf = True
+  checkLeft binds v (Bound sop) = evalSoP sop binds <= v
+
+  checkRight ::
+    (Ord f, Ord c, Ord n, Floating n) => Map (Atom f c) n -> n -> Bound f c -> Bool
+  checkRight _ _ Inf = True
+  checkRight binds v (Bound sop) = v <= evalSoP sop binds
+
+-- | Assert if given expression holds in the current environment
+assert ::
+  (Ord f, Ord c) =>
+  -- | Asserted expression
+  SoPE f c ->
+  -- | - True - if expressions holds
+  --   - False - otherwise
+  --
+  -- State will become @Nothing@ if it cannot reason about these kind of expressions
+  SolverState f c Bool
+assert =
+  declareToState >=> \SoPE{..} ->
+    case op of
+      EqR -> assertEq lhs rhs
+      LeR -> do
+        r1 <- assertEq lhs rhs
+        if r1
+          then return True
+          else do
+            assertRange lhs rhs <|> assertNewton lhs rhs
+      GeR -> do
+        r1 <- assertEq lhs rhs
+        if r1
+          then return True
+          else do
+            assertRange rhs lhs <|> assertNewton rhs lhs
+
+{- | Get unifiers for an expression
+minimal set of expressions that should hold for the expression to hold
+-}
+unify ::
+  (Ord f, Ord c) =>
+  -- | Unified expression
+  SoPE f c ->
+  -- | List of unifiers - Minimal list of unifiers for the expression to hold.
+  -- The list is empty, if it never holds
+  --
+  -- State will always be valid after a call
+  SolverState f c [SoPE f c]
+unify =
+  declareToState >=> \expr@SoPE{..} ->
+    case op of
+      EqR -> return (filter (/= expr) $ map unifier2SoPE (unifiers lhs rhs))
+      _ -> return []
+ where
+  unifier2SoPE Subst{..} = SoPE (symbol sConst) sSoP EqR
+
+-- | Get range of possible values for an expression
+range ::
+  (Ord f, Ord c) =>
+  -- | Expression
+  SoP f c ->
+  -- | (lower bound, upper bound) - Range for an expression
+  --
+  -- @Nothing@ means that the expression is unbounded
+  -- from that side
+  SolverState f c (Maybe (SoP f c), Maybe (SoP f c))
+range sop = do
+  _ <- declareSoP sop
+  (Range low up) <- getRangeSoP sop
+  return (boundSoP low, boundSoP up)
+
+-- | Get list of all ranges stored in a state
+ranges ::
+  (Ord f, Ord c) =>
+  -- | (lower bound, symbol, upper bound) - Similar to @range@
+  --     but also provides expression
+  SolverState f c [(Maybe (SoP f c), SoP f c, Maybe (SoP f c))]
+ranges =
+  map (\(a, Range low up) -> (boundSoP low, symbol a, boundSoP up)) . M.toList <$> getRanges
diff --git a/src/SoPSat/SoP.hs b/src/SoPSat/SoP.hs
new file mode 100644
--- /dev/null
+++ b/src/SoPSat/SoP.hs
@@ -0,0 +1,226 @@
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE RecordWildCards #-}
+
+module SoPSat.SoP (
+  -- * SoP Types
+  Atom,
+  Symbol,
+  Product,
+  SoP,
+  SoPE (..),
+  ToSoP (..),
+
+  -- * Operators
+  (|+|),
+  (|-|),
+  (|*|),
+  (|/|),
+  (|^|),
+
+  -- * Relations
+  OrdRel (..),
+
+  -- * Related
+  constants,
+  atoms,
+  int,
+  cons,
+  symbol,
+  func,
+
+  -- * Predicates
+  isConst,
+  isFunction,
+)
+where
+
+import Data.Set (Set, union)
+import qualified Data.Set as S
+
+import SoPSat.Internal.SoP
+
+{- | Convertable to a sum of products
+with `f` being type to represent functions
+and `c` being type to represent constants
+-}
+class (Ord f, Ord c) => ToSoP f c a where
+  toSoP :: a -> SoP f c
+
+-- | Predicate for constant @Atom@s
+isConst :: Atom f c -> Bool
+isConst (C _) = True
+isConst _ = False
+
+-- | Predicate for function @Atom@s
+isFunction :: Atom f c -> Bool
+isFunction (F _ _) = True
+isFunction _ = False
+
+instance (Ord f, Ord c) => ToSoP f c (Symbol f c) where
+  toSoP s = simplifySoP $ S [P [s]]
+
+instance (Ord f, Ord c) => ToSoP f c (Product f c) where
+  toSoP p = simplifySoP $ S [p]
+
+instance (Ord f, Ord c) => ToSoP f c (SoP f c) where
+  toSoP = simplifySoP
+
+-- | Order relationship
+data OrdRel
+  = -- | Less than or equal relationship
+    LeR
+  | -- | Equality relationship
+    EqR
+  | -- | Greater than or equal relationship
+    GeR
+  deriving (Eq, Ord)
+
+instance Show OrdRel where
+  show LeR = "<="
+  show EqR = "="
+  show GeR = ">="
+
+-- | Expression
+data SoPE f c
+  = SoPE
+  { lhs :: SoP f c
+  -- ^ Left hand side of the expression
+  , rhs :: SoP f c
+  -- ^ Right hand side of the expression
+  , op :: OrdRel
+  -- ^ Relationship between sides
+  }
+
+instance (Eq f, Eq c) => Eq (SoPE f c) where
+  (SoPE l1 r1 op1) == (SoPE l2 r2 op2)
+    | op1 == op2
+    , op1 == EqR =
+        -- a = b is the same as b = a
+        (l1 == l2) && (r1 == r2) || (l1 == r2) && (r1 == l2)
+    | op1 == op2 =
+        -- (a <= b) is itself
+        (l1 == l2) && (r1 == r2)
+    | EqR `notElem` [op1, op2] =
+        -- (a <= b) is the same as (b >= a)
+        (l1 == r2) && (r1 == l2)
+    | otherwise =
+        False
+
+instance (Show f, Show c) => Show (SoPE f c) where
+  show SoPE{..} = unwords [show lhs, show op, show rhs]
+
+-- | Creates an integer expression
+int :: Integer -> SoP f c
+int i = S [P [I i]]
+
+-- | Creates expression from an atom
+symbol :: Atom f c -> SoP f c
+symbol a = S [P [A a]]
+
+-- | Creates a constant expression
+cons :: c -> SoP f c
+cons c = S [P [A (C c)]]
+
+-- | Creates a function expression
+func :: (Ord f, Ord c) => f -> [SoP f c] -> SoP f c
+func f args = S [P [A (F f (map simplifySoP args))]]
+
+infixr 8 |^|
+
+-- | Exponentiation of @SoP@s
+(|^|) :: (Ord f, Ord c) => SoP f c -> SoP f c -> SoP f c
+-- It's a B2 combinator,
+(|^|) = (. simplifySoP) . normaliseExp
+
+infixl 6 |+|
+
+-- | Addition of @SoP@s
+(|+|) :: (Ord f, Ord c) => SoP f c -> SoP f c -> SoP f c
+(|+|) = mergeSoPAdd
+
+infixl 7 |*|
+
+-- | Multiplication of @SoP@s
+(|*|) :: (Ord f, Ord c) => SoP f c -> SoP f c -> SoP f c
+(|*|) = mergeSoPMul
+
+infixl 6 |-|
+
+-- | Subtraction of @SoP@s
+(|-|) :: (Ord f, Ord c) => SoP f c -> SoP f c -> SoP f c
+(|-|) = mergeSoPSub
+
+infixl 7 |/|
+
+{- | Division of @SoP@s
+
+Produces a tuple of a quotient and a remainder
+NB. Not implemented
+-}
+(|/|) :: (Ord f, Ord c) => SoP f c -> SoP f c -> (SoP f c, SoP f c)
+(|/|) = mergeSoPDiv
+
+-- | Collects @Atom@s used in a @SoP@
+atoms :: (Ord f, Ord c) => SoP f c -> Set (Atom f c)
+atoms = S.unions . map atomsProduct . unS
+
+{- | Collects @Atom@s used in a @Product@
+
+Used by @atoms@
+-}
+atomsProduct :: (Ord f, Ord c) => Product f c -> Set (Atom f c)
+atomsProduct = S.unions . map atomsSymbol . unP
+
+{- | Collect @Atom@s used in @Symbol@s
+
+Used by @atomsProduct@
+-}
+atomsSymbol ::
+  (Ord f, Ord c) =>
+  Symbol f c ->
+  -- | - Empty - if the symbol is an integer
+  --   - Singleton - if the symbol is an atom
+  --   - Set of symbols - if the symbol is an exponentiation
+  Set (Atom f c)
+atomsSymbol (I _) = S.empty
+atomsSymbol (A a) = S.singleton a
+atomsSymbol (E b p) = atoms b `union` atomsProduct p
+
+{- | Collects constants used in @SoP@
+
+Almost equivalent to
+@Data.Set.filter isConst . atoms@
+, but also collects constants used in functions
+-}
+constants :: (Ord f, Ord c) => SoP f c -> Set c
+constants = S.unions . map constsProduct . unS
+
+{- | Collects constants used in @Product@
+
+Used by @constants@
+-}
+constsProduct :: (Ord f, Ord c) => Product f c -> Set c
+constsProduct = S.unions . map constsSymbol . unP
+
+{- | Collects constants used in @Symbol@
+
+Used by @constsProduct@
+-}
+constsSymbol :: (Ord f, Ord c) => Symbol f c -> Set c
+constsSymbol (I _) = S.empty
+constsSymbol (A a) = constsAtom a
+constsSymbol (E b p) = constants b `union` constsProduct p
+
+{- | Collects constants used in @Atom@
+
+Used by @constsSymbol@
+-}
+constsAtom ::
+  (Ord f, Ord c) =>
+  Atom f c ->
+  -- | Singleton - if the atom is a constant
+  --   Set of constants - if the atom is a function
+  Set c
+constsAtom (C c) = S.singleton c
+constsAtom (F _ args) = S.unions $ map constants args
diff --git a/tests/SystemTests.hs b/tests/SystemTests.hs
new file mode 100644
--- /dev/null
+++ b/tests/SystemTests.hs
@@ -0,0 +1,597 @@
+import Test.Tasty
+import Test.Tasty.HUnit (
+  testCase,
+  (@=?),
+ )
+
+import Data.Monoid (Any)
+import SoPSat.Satisfier (
+  SolverState,
+  assert,
+  declare,
+  evalStatements,
+  unify,
+ )
+import SoPSat.SoP (
+  OrdRel (..),
+  SoPE (..),
+  (|*|),
+  (|+|),
+  (|-|),
+  (|^|),
+ )
+import qualified SoPSat.SoP as SoP
+
+type SolveTestCase = SolverState String String Bool
+type SolveTestResult = Maybe Bool
+
+type UnifyTestCase = SolverState String String [SoPE String String]
+type UnifyTestResult = Maybe [SoPE String String]
+
+equalityGiven1 :: SolveTestCase
+equalityGiven1 =
+  let
+    one = SoP.int 1
+    m = SoP.cons "m"
+    n = SoP.cons "n"
+    n1 = SoP.cons "n1"
+   in
+    do
+      declare (SoPE m (n1 |+| one) EqR)
+      assert (SoPE (m |+| n) (n |+| n1 |+| one) EqR)
+
+runEqualityGiven1 :: SolveTestResult
+runEqualityGiven1 = evalStatements equalityGiven1
+
+equalityGiven2 :: SolveTestCase
+equalityGiven2 =
+  let
+    one = SoP.int 1
+    m = SoP.cons "m"
+    n = SoP.cons "n"
+    n1 = SoP.cons "n1"
+   in
+    do
+      declare (SoPE m (n1 |+| one) EqR)
+      assert (SoPE (m |*| n) (n |+| n |*| n1) EqR)
+
+runEqualityGiven2 :: SolveTestResult
+runEqualityGiven2 = evalStatements equalityGiven2
+
+equalityGiven3 :: SolveTestCase
+equalityGiven3 =
+  let
+    one = SoP.int 1
+    m = SoP.cons "m"
+    n = SoP.cons "n"
+    n1 = SoP.cons "n1"
+   in
+    do
+      declare (SoPE m (n1 |+| one) EqR)
+      assert (SoPE (n |^| m) (n |*| n |^| n1) EqR)
+
+runEqualityGiven3 :: SolveTestResult
+runEqualityGiven3 = evalStatements equalityGiven3
+
+transitivity :: SolveTestCase
+transitivity =
+  let
+    i = SoP.cons "i"
+    j = SoP.cons "j"
+    k = SoP.cons "k"
+   in
+    do
+      declare (SoPE i j LeR)
+      declare (SoPE j k LeR)
+      assert (SoPE i k LeR)
+
+runTransitivity :: SolveTestResult
+runTransitivity = evalStatements transitivity
+
+antisymmetryZero :: SolveTestCase
+antisymmetryZero =
+  let
+    z = SoP.int 0
+    x = SoP.cons "x"
+   in
+    do
+      declare (SoPE x z LeR)
+      assert (SoPE x z EqR)
+
+runAntisymmetryZero :: SolveTestResult
+runAntisymmetryZero = evalStatements antisymmetryZero
+
+antisymmetryNonZero :: SolveTestCase
+antisymmetryNonZero =
+  let
+    z = SoP.int 42
+    x = SoP.cons "x"
+   in
+    do
+      declare (SoPE x z LeR)
+      declare (SoPE z x LeR)
+      assert (SoPE x z EqR)
+
+runAntisymmetryNonZero :: SolveTestResult
+runAntisymmetryNonZero = evalStatements antisymmetryNonZero
+
+lemma2 :: SolveTestCase
+lemma2 =
+  let
+    o = SoP.int 1
+    j = SoP.cons "j"
+    n = SoP.cons "n"
+   in
+    do
+      declare (SoPE j n LeR)
+      declare (SoPE o (n |-| j) LeR)
+      assert (SoPE (o |+| j) n LeR)
+
+runLemma2 :: SolveTestResult
+runLemma2 = evalStatements lemma2
+
+trueInEq :: SolveTestCase
+trueInEq =
+  let
+    two = SoP.int 2
+    three = SoP.int 3
+    four = SoP.int 4
+    x = SoP.cons "x"
+    inEq1 = two |^| x |+| three |*| x |^| two |+| three
+    inEq2 = x |^| three |-| two |*| x |^| two |+| four
+   in
+    assert (SoPE inEq2 inEq1 LeR)
+
+runTrueInEq :: SolveTestResult
+runTrueInEq = evalStatements trueInEq
+
+falseInEq :: SolveTestCase
+falseInEq =
+  let
+    two = SoP.int 2
+    three = SoP.int 3
+    four = SoP.int 4
+    x = SoP.cons "x"
+    inEq1 = two |^| x |+| x |^| two |+| three
+    inEq2 = x |^| three |-| two |*| x |^| two |+| four
+   in
+    assert (SoPE inEq1 inEq2 GeR)
+
+runFalseInEq :: SolveTestResult
+runFalseInEq = evalStatements falseInEq
+
+falseInEq2 :: SolveTestCase
+falseInEq2 =
+  let
+    one = SoP.int 1
+    m = SoP.cons "m"
+    rp = SoP.cons "rp"
+   in
+    do
+      declare (SoPE one m LeR)
+      declare (SoPE m rp LeR)
+      assert (SoPE one (rp |-| m) LeR)
+
+runFalseInEq2 :: SolveTestResult
+runFalseInEq2 = evalStatements falseInEq2
+
+overlapInEq :: SolveTestCase
+overlapInEq =
+  let
+    t = SoP.int 2
+    f = SoP.int 4
+    x = SoP.cons "x"
+   in
+    do
+      declare (SoPE f x LeR)
+      declare (SoPE t x LeR)
+      assert (SoPE t x LeR)
+
+runOverlapInEq :: SolveTestResult
+runOverlapInEq = evalStatements overlapInEq
+
+eqSubst :: SolveTestCase
+eqSubst =
+  let
+    o = SoP.int 1
+    x = SoP.cons "x"
+    m = SoP.cons "m"
+    m1 = SoP.cons "m1"
+    n1 = SoP.cons "n1"
+    n2 = SoP.cons "n2"
+   in
+    do
+      declare (SoPE (x |+| o) (n1 |+| m |+| o) EqR)
+      declare (SoPE m n1 EqR)
+      declare (SoPE n1 (n2 |+| m1 |+| o) EqR)
+      assert (SoPE (o |+| n2 |+| m1) n1 EqR)
+
+runEqSubst :: SolveTestResult
+runEqSubst = evalStatements eqSubst
+
+eqSubst2 :: SolveTestCase
+eqSubst2 =
+  let
+    o = SoP.int 1
+    t = SoP.int 2
+    y = SoP.cons "y"
+    x = SoP.cons "x"
+   in
+    do
+      declare (SoPE o y LeR)
+      declare (SoPE (o |+| x) (t |*| y) EqR)
+      assert (SoPE (t |*| y |-| o) x EqR)
+
+runEqSubst2 :: SolveTestResult
+runEqSubst2 = evalStatements eqSubst2
+
+multistep :: SolveTestCase
+multistep =
+  let
+    o = SoP.int 1
+    t = SoP.int 2
+    n1 = SoP.cons "n1"
+    n2 = SoP.cons "n2"
+   in
+    do
+      declare (SoPE o n1 LeR)
+      declare (SoPE (t |*| n2) n1 EqR)
+      r1 <- assert (SoPE (t |*| n2 |-| o) (n1 |-| o) EqR)
+      r2 <- assert (SoPE o n2 LeR)
+      return (r1 && r2)
+
+runMultistep :: SolveTestResult
+runMultistep = evalStatements multistep
+
+multistep2 :: SolveTestCase
+multistep2 =
+  let
+    o = SoP.int 1
+    m = SoP.cons "m"
+    n = SoP.cons "n"
+    n1 = SoP.cons "n1"
+    n2 = SoP.cons "n2"
+   in
+    do
+      declare (SoPE m (n1 |+| o) EqR)
+      declare (SoPE (m |+| n) (n2 |+| o) EqR)
+      assert (SoPE (n1 |+| n) n2 EqR)
+
+runMultistep2 :: SolveTestResult
+runMultistep2 = evalStatements multistep2
+
+step3 :: SolveTestCase
+step3 =
+  let
+    o = SoP.int 1
+    t = SoP.int 2
+    m = SoP.cons "m"
+    n = SoP.cons "n"
+    n1 = SoP.cons "n1"
+   in
+    do
+      declare (SoPE (o |+| n) m EqR)
+      declare (SoPE (o |+| n1) m EqR)
+      r1 <- assert (SoPE n n1 EqR)
+      r2 <- assert (SoPE (t |+| t |*| n) (t |*| m) EqR)
+      return (r1 && r2)
+
+runStep3 :: SolveTestResult
+runStep3 = evalStatements step3
+
+implication :: SolveTestCase
+implication =
+  let
+    t = SoP.int 2
+    m = SoP.cons "m"
+    n = SoP.cons "n"
+   in
+    do
+      declare (SoPE t (t |^| (n |+| m)) LeR)
+      assert (SoPE t (t |^| (m |+| n)) LeR)
+
+runImplication :: SolveTestResult
+runImplication = evalStatements implication
+
+func1 :: SolveTestCase
+func1 =
+  let
+    a = SoP.cons "a"
+    b = SoP.cons "b"
+    c = SoP.cons "c"
+    eq1 = a |+| SoP.func "max" [a |+| b, c]
+    eq2 = SoP.func "max" [b |+| a, c] |+| a
+   in
+    assert (SoPE eq1 eq2 EqR)
+
+runFunc1 :: SolveTestResult
+runFunc1 = evalStatements func1
+
+func2 :: SolveTestCase
+func2 =
+  let
+    a = SoP.cons "a"
+    b = SoP.cons "b"
+    c = SoP.cons "c"
+    eq1 = a |+| SoP.func "bar" [a |+| b, c]
+    eq2 = SoP.func "bar" [c, b |+| a] |+| a
+   in
+    assert (SoPE eq1 eq2 EqR)
+
+runFunc2 :: SolveTestResult
+runFunc2 = evalStatements func2
+
+func3 :: SolveTestCase
+func3 =
+  let
+    a = SoP.cons "a"
+    b = SoP.cons "b"
+    c = SoP.cons "c"
+   in
+    assert (SoPE (SoP.func "foo" [a, b, a |^| c]) a GeR)
+
+runFunc3 :: SolveTestResult
+runFunc3 = evalStatements func3
+
+func4 :: SolveTestCase
+func4 =
+  let
+    x = SoP.cons "x"
+    g = SoP.func "g" [x]
+    f = SoP.func "f" [x]
+   in
+    do
+      declare (SoPE (SoP.int 1) g LeR)
+      assert (SoPE f (f |*| g) LeR)
+
+runFunc4 :: SolveTestResult
+runFunc4 = evalStatements func4
+
+unifyExp :: UnifyTestCase
+unifyExp =
+  let
+    t = SoP.int 2
+    x1 = SoP.cons "x1"
+    x2 = SoP.cons "x2"
+   in
+    do
+      unify (SoPE ((t |^| x1) |*| (t |^| (x1 |+| x1))) ((t |^| x2) |*| (t |^| (x2 |+| x2))) EqR)
+
+runUnifyExp :: UnifyTestResult
+runUnifyExp = evalStatements unifyExp
+
+main :: IO ()
+main = defaultMain tests
+
+tests :: TestTree
+tests =
+  testGroup
+    "lib-tests"
+    [ testGroup
+        "Equality tests"
+        [ testGroup
+            "True"
+            [ testCase "m = n1 + 1 implies n + m = n + n1 + 1" $
+                Just True @=? runEqualityGiven1
+            , testCase "m = n1 + 1 implies n * m = n + n * n1" $
+                Just True @=? runEqualityGiven2
+            , testCase "m = n1 + 1 implies n^m = n*n^n1" $
+                Just True @=? runEqualityGiven3
+            , testCase "n + 1 = n1 + m + 1 and m = n1 and n1 = n2 + m1 + 1 implies 1 + n2 + m1 = n1" $
+                Just True @=? runEqSubst
+            , testCase "1 <= y and x + 1 = 2 * y implies 2 * y - 1 = x" $
+                Just True @=? runEqSubst2
+            , testCase "Combined: 1 <= m and 2 * n = m implies 2 * n - 1 = m - 2 and 1 <= m" $
+                Just True @=? runMultistep
+            , testCase "Multistep: m = n1 + 1 and m + n = n2 + 1 implies n1 + n = n2" $
+                Just True @=? runMultistep2
+            , testCase "1 + a = c and 1 + b = c implies a = b and 2 + 2 * a = 2 * c" $
+                Just True @=? runStep3
+            , testGroup
+                "Functions"
+                [ testCase "a + max(a + b, c) = max(b + a, c) + a" $
+                    Just True @=? runFunc1
+                ]
+            , testGroup
+                "Natural numbers"
+                [ testCase "m = n - 1 implies m + 1 = n" $
+                    Nothing
+                      @=? evalStatements
+                        ( declare (SoPE (SoP.cons "m") (SoP.cons "n" |-| SoP.int 1) EqR)
+                            >> assert (SoPE (SoP.cons "m" |+| SoP.int 1) (SoP.cons "n") EqR) ::
+                            SolveTestCase
+                        )
+                ]
+            ]
+        , testGroup
+            "False"
+            [ testCase "x + 2 /= x + 3" $
+                Just False
+                  @=? evalStatements
+                    ( assert
+                        (SoPE (SoP.cons "x" |+| SoP.int 2) (SoP.cons "x" |+| SoP.int 3) EqR) ::
+                        SolveTestCase
+                    )
+            , testCase "8 /= x + x + x" $
+                Just False
+                  @=? evalStatements
+                    ( assert
+                        (SoPE (SoP.int 3 |*| SoP.cons "x") (SoP.int 8) EqR) ::
+                        SolveTestCase
+                    )
+            , testCase "7 /= 2*y+4" $
+                Just False
+                  @=? evalStatements
+                    ( assert
+                        (SoPE (SoP.int 2 |*| SoP.cons "y" |+| SoP.int 4) (SoP.int 7) EqR) ::
+                        SolveTestCase
+                    )
+            , testGroup
+                "Functions"
+                [ testCase "a + bar(a + b, c) /= bar(c, b + a) + a" $
+                    Just False @=? runFunc2
+                ]
+            ]
+        ]
+    , testGroup
+        "Inequality tests"
+        [ testGroup
+            "True"
+            [ testCase "Transitivity: i <= j and j <= k implies i <= k" $
+                Just True @=? runTransitivity
+            , testCase "Antisymmetry with zero: x is Natural and x <= 0 implies x = 0" $
+                Just True @=? runAntisymmetryZero
+            , testCase "Antisymmetry with non-zero: x <= 5 and x >= 5 implies x = 5" $
+                Just True @=? runAntisymmetryNonZero
+            , testCase "Strongly greater: j <= n and 1 <= n - j imples 1 + j <= n" $
+                Just True @=? runLemma2
+            , testCase "Composite function: x^3-2x^2+4<=2^x+3x^2+3" $
+                Just True @=? runTrueInEq
+            , testCase "Overlapping ranges: 4 <= x implies 2 <= x" $
+                Just True @=? runOverlapInEq
+            , testGroup
+                "Functions"
+                [ testCase "g(x) >= 1 implies f(x) <= g(x) * f(x)" $
+                    Just True @=? runFunc4
+                ]
+            , testGroup
+                "Trivial"
+                [ testCase "a <= a + 1" $
+                    Just True
+                      @=? evalStatements
+                        ( assert
+                            (SoPE (SoP.cons "a") (SoP.cons "a" |+| SoP.int 1) LeR) ::
+                            SolveTestCase
+                        )
+                , testCase "1 <= 2^a" $
+                    Just True
+                      @=? evalStatements
+                        ( assert
+                            (SoPE (SoP.int 1) (SoP.int 2 |^| SoP.cons "a") LeR) ::
+                            SolveTestCase
+                        )
+                , testCase "2 <= 2^(n + m) implies 2 <= 2^(m + n)" $
+                    Just True @=? runImplication
+                ]
+            , testGroup
+                "Implications"
+                [ testCase "a = b - 1 implies b >= 1" $
+                    Nothing
+                      @=? evalStatements
+                        ( declare
+                            (SoPE (SoP.cons "a") (SoP.cons "b" |-| SoP.int 1) EqR)
+                            >> assert (SoPE (SoP.cons "b") (SoP.int 1) GeR) ::
+                            SolveTestCase
+                        )
+                ]
+            ]
+        , testGroup
+            "False"
+            [ testCase "Composite function x^3-2x^2+4<=2^x+x^2+3" $
+                Just False @=? runFalseInEq
+            , testCase "1 <= m and m <= rp implies 1 <= rp - m" $
+                Just False @=? runFalseInEq2
+            , testCase "4a <= 2a" $
+                Just False
+                  @=? evalStatements
+                    ( assert
+                        (SoPE (SoP.int 4 |*| SoP.cons "a") (SoP.int 2 |*| SoP.cons "a") LeR) ::
+                        SolveTestCase
+                    )
+            , testCase "foo(a, b, a ^ c) >= a" $
+                Just False @=? runFunc3
+            ]
+        ]
+    , testGroup
+        "Ranges"
+        []
+    , -- TODO: Add test cases for range narrowing consistency
+
+      testGroup
+        "Unifiers"
+        [ testCase "x = x always holds" $
+            Just []
+              @=? evalStatements (unify (SoPE (SoP.cons "x") (SoP.cons "x") EqR) :: UnifyTestCase)
+        , testCase "t = a + b does not produce unifiers" $
+            Just []
+              @=? evalStatements
+                ( unify (SoPE (SoP.cons "t") (SoP.cons "a" |+| SoP.cons "b") EqR) ::
+                    UnifyTestCase
+                )
+        , testCase "a + b = a + c if b = c" $
+            Just [SoPE (SoP.cons "b") (SoP.cons "c") EqR]
+              @=? evalStatements
+                ( unify
+                    ( SoPE
+                        (SoP.cons "a" |+| SoP.cons "b")
+                        (SoP.cons "a" |+| SoP.cons "c")
+                        EqR
+                    ) ::
+                    UnifyTestCase
+                )
+        , testCase "n = n + d" $
+            Just []
+              @=? evalStatements
+                ( unify
+                    ( SoPE
+                        (SoP.cons "n")
+                        (SoP.cons "n" |+| SoP.cons "d")
+                        EqR
+                    ) ::
+                    UnifyTestCase
+                )
+        , testCase "c = n implies n = n + d never holds" $
+            Just []
+              @=? evalStatements
+                ( declare (SoPE (SoP.cons "c") (SoP.cons "n") EqR)
+                    >> unify
+                      ( SoPE
+                          (SoP.cons "n" |+| SoP.cons "d")
+                          (SoP.cons "n")
+                          EqR
+                      ) ::
+                    UnifyTestCase
+                )
+        , testCase "9 = x + x + x if x = 3" $
+            Just [SoPE (SoP.cons "x") (SoP.int 3) EqR]
+              @=? evalStatements
+                ( unify
+                    (SoPE (SoP.int 3 |*| SoP.cons "x") (SoP.int 9) EqR) ::
+                    UnifyTestCase
+                )
+        , testCase "6 = 2 * y + 4 if x = 1" $
+            Just [SoPE (SoP.cons "y") (SoP.int 1) EqR]
+              @=? evalStatements
+                ( unify
+                    (SoPE (SoP.int 2 |*| SoP.cons "y" |+| SoP.int 4) (SoP.int 6) EqR) ::
+                    UnifyTestCase
+                )
+        , testCase "8 /= x + x + x never holds" $
+            Just []
+              @=? evalStatements
+                ( unify
+                    (SoPE (SoP.int 3 |*| SoP.cons "x") (SoP.int 8) EqR) ::
+                    UnifyTestCase
+                )
+        , testCase "7 /= 2*y+4 never holds" $
+            Just []
+              @=? evalStatements
+                ( unify
+                    (SoPE (SoP.int 2 |*| SoP.cons "y" |+| SoP.int 4) (SoP.int 7) EqR) ::
+                    UnifyTestCase
+                )
+        , testCase "a^b = a^c if b = c" $
+            Just [SoPE (SoP.cons "b") (SoP.cons "c") EqR]
+              @=? evalStatements
+                ( unify
+                    ( SoPE
+                        (SoP.cons "a" |^| SoP.cons "b")
+                        (SoP.cons "a" |^| SoP.cons "c")
+                        EqR
+                    ) ::
+                    UnifyTestCase
+                )
+        , testCase "2^x1 * 2^(x1 + x1) = 2^x2 * 2^(x2 + x2) holds if x1 = x2" $
+            Just [SoPE (SoP.cons "x1") (SoP.cons "x2") EqR] @=? runUnifyExp
+        ]
+    ]