diff --git a/ChangeLog.md b/ChangeLog.md
new file mode 100644
--- /dev/null
+++ b/ChangeLog.md
@@ -0,0 +1,5 @@
+# Revision history for hafar
+
+## 0.1.0.0 -- 2020-02-13
+
+* First version. Released on an unsuspecting world.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2019, Joosep Jääger
+
+All rights reserved.
+
+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.
+
+    * Neither the name of Joosep Jääger nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+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,37 @@
+# Hafar
+
+Hafar is an implementation of affine arithmetic in haskell.
+
+## Building
+
+To build the library, simply run 
+```
+# stack build
+```
+or if using cabal 
+```
+# cabal install --only-dependencies
+# cabal build
+```
+
+## Example
+
+All operations with affine forms must be done inside the AFM monad.
+
+```
+import Numeric.Interval hiding (interval)
+
+x1 = do
+  a <- newFromInterval $ 4...6
+  b <- newFromInterval $ 4...6
+  return . interval $ a - b
+
+evalAFM x1 -- evaluates to approximately -2 ... 2
+
+x2 = do
+  a <- newFromInterval $ 4...6
+  return . interval $ a - a
+
+evalAFM x2 -- evaluates to approximately 0 ... 0
+
+```
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/hafar.cabal b/hafar.cabal
new file mode 100644
--- /dev/null
+++ b/hafar.cabal
@@ -0,0 +1,63 @@
+cabal-version: 1.12
+
+-- This file has been generated from package.yaml by hpack version 0.31.2.
+--
+-- see: https://github.com/sol/hpack
+--
+-- hash: 11f7d969e8f00b5e392eda2b0897debe20363e72178cbd6a8a87a7ab29e2b819
+
+name:           hafar
+version:        0.1.0.0
+synopsis:       Affine arithmetic library for Haskell
+description:    Hafar is an affine arithmetic library for Haskell. It is an efficient way to work with ranges of values or imprecise values.
+category:       Numeric
+homepage:       https://github.com/Soupstraw/hafar#readme
+bug-reports:    https://github.com/Soupstraw/hafar/issues
+author:         Joosep Jääger
+maintainer:     Joosep Jääger
+copyright:      2019 Joosep Jääger
+license:        BSD3
+license-file:   LICENSE
+build-type:     Simple
+extra-source-files:
+    README.md
+    ChangeLog.md
+
+source-repository head
+  type: git
+  location: https://github.com/Soupstraw/hafar
+
+library
+  exposed-modules:
+      Numeric.AffineForm
+      Numeric.AffineForm.ExplicitRounding
+  other-modules:
+      Numeric.AffineForm.Utils
+      Numeric.AffineForm.Internal
+  hs-source-dirs:
+      src
+  build-depends:
+      base >=4.12 && <4.14
+    , intervals >=0.8 && <0.9
+    , mtl >=2.2 && <2.3
+  default-language: Haskell2010
+
+test-suite hafar-test
+  type: exitcode-stdio-1.0
+  main-is: Spec.hs
+  other-modules:
+      Numeric.AffineForm
+      Numeric.AffineForm.Internal
+      Numeric.AffineForm.Utils
+      Numeric.AffineForm.ExplicitRounding
+  hs-source-dirs:
+      test
+      src
+  ghc-options: -threaded -rtsopts -with-rtsopts=-N
+  build-depends:
+      QuickCheck >=2.13 && <2.14
+    , base >=4.12 && <4.14
+    , hafar
+    , intervals >=0.8 && <0.9
+    , mtl >=2.2 && <2.3
+  default-language: Haskell2010
diff --git a/src/Numeric/AffineForm.hs b/src/Numeric/AffineForm.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AffineForm.hs
@@ -0,0 +1,17 @@
+module Numeric.AffineForm (AFM, AF, newEps,
+                           newFromInterval,
+                           singleton,
+                           evalAFM,
+                           radius,
+                           midpoint,
+                           inf, sup,
+                           interval,
+                           member,
+                           epscount_,
+                           setMidpoint,
+                           fix,
+                           addError,
+                           (.+), (.*)
+                          ) where
+
+import Numeric.AffineForm.Internal
diff --git a/src/Numeric/AffineForm/ExplicitRounding.hs b/src/Numeric/AffineForm/ExplicitRounding.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AffineForm/ExplicitRounding.hs
@@ -0,0 +1,57 @@
+-- | ExplicitRounding defines the ExplicitRounding class and
+-- instances for some more common numeric types
+module Numeric.AffineForm.ExplicitRounding (
+                                ExplicitRounding,
+                                eps, prev, next,
+                                (+/), (+\),
+                                (-/), (-\),
+                                (*/), (*\)
+                                ) where
+
+import Numeric.Interval as IA
+import Data.Ratio
+
+-- | The class of numeric values that can be rounded explicitly
+class (Ord a, Num a) => ExplicitRounding a where
+  -- | Return some number so that all the values that could be rounded to the parameter of this function
+  -- would be at most that distance away from that parameter.
+  eps :: a -> a
+  -- | Returns parameter plus its epsilon
+  prev :: a -> a
+  -- | Returns parameter minus its epsilon
+  next :: a -> a
+  -- | Add the two values, rounding the result up
+  (+/) :: a -> a -> a
+  -- | Add the two values, rounding the result down
+  (+\) :: a -> a -> a
+  -- | Subtract the two values, rounding the result up
+  (-/) :: a -> a -> a
+  -- | Subtract the two values, rounding the result down
+  (-\) :: a -> a -> a
+  -- | Multiply the two values, rounding the result up
+  (*/) :: a -> a -> a
+  -- | Multiply the two values, rounding the result down
+  (*\) :: a -> a -> a
+
+  prev x     = x - eps x
+  next x     = x + eps x
+  x +/ y     = next $ x + y
+  x +\ y     = prev $ x + y
+  x -/ y     = next $ x - y
+  x -\ y     = prev $ x - y
+  x */ y     = next $ x * y
+  x *\ y     = prev $ x * y
+
+instance ExplicitRounding Int where
+  eps = const 0
+
+instance (Integral a, ExplicitRounding a) => ExplicitRounding (Ratio a) where
+  eps x = (eps $ numerator x) % (abs (denominator x) - (eps $ denominator x))
+
+instance ExplicitRounding Float where
+  eps 0 = eps $ 2e-36
+  eps x = encodeFloat 2 (snd $ decodeFloat x)
+
+instance ExplicitRounding Double where
+  eps 0 = eps $ 1e-300
+  eps x = encodeFloat 2 (snd $ decodeFloat x)
diff --git a/src/Numeric/AffineForm/Internal.hs b/src/Numeric/AffineForm/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AffineForm/Internal.hs
@@ -0,0 +1,270 @@
+{-# LANGUAGE RankNTypes#-}
+{-# LANGUAGE ExistentialQuantification #-}
+{-# LANGUAGE KindSignatures #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE DataKinds, ScopedTypeVariables #-}
+
+-- | This module defines the affine form, the AFM monad
+-- and many operations for affine forms.
+module Numeric.AffineForm.Internal where
+
+import Control.Monad.State hiding (fix)
+import Control.Monad.Identity hiding (fix)
+import Control.Exception (Exception, throw, evaluate, try)
+
+import Numeric.AffineForm.Utils
+import Numeric.AffineForm.ExplicitRounding
+import qualified Numeric.Interval as IA
+import Numeric.Interval ((...))
+import Data.Fixed (mod')
+import Data.Ratio (approxRational, (%))
+import Data.Either (fromLeft, fromRight)
+
+-- | An affine form is defined by its midpoint, list of epsilon coefficients and an error coefficient
+data AF s a
+  = AF a [a] a
+  deriving (Show)
+
+data Curvature = Convex | Concave
+
+data AFException
+  = DivisionByZero
+  | LogFromNegative
+  | AddingNegativeError
+
+instance Show AFException where
+  show DivisionByZero = "division by zero"
+  show LogFromNegative = "logarithm from a negative number"
+  show AddingNegativeError = "cannot add a negative error to an affine form"
+
+instance Exception AFException
+
+instance (Fractional a, ExplicitRounding a, Ord a) => Num (AF s a) where
+  (+) = add
+  (*) = multiply
+  abs = absAF
+  signum = signumAF
+  fromInteger = singleton . fromInteger
+  negate = negateAF
+
+instance (Fractional a, ExplicitRounding a, Ord a) => Fractional (AF s a) where
+  recip = recipAF
+  fromRational = singleton . fromRational
+
+instance (Floating a, RealFrac a, ExplicitRounding a, Ord a) => Floating (AF s a) where
+  pi = approxSingleton pi
+  exp = minrange exp exp Convex
+  log x
+    | inf x > 0  = minrange log recip Concave x
+    | otherwise = throw LogFromNegative
+  sin = sinAF
+  cos = cosAF
+  asin = minrange asin (\x -> 1/sqrt (1-x^2)) undefined
+  acos = minrange acos (\x -> -1/sqrt (1-x^2)) undefined
+  atan = minrange atan (\x -> 1/(x^2+1)) undefined
+  sinh = minrange sinh cosh undefined
+  cosh = minrange cosh sinh Convex
+  asinh = minrange asinh (\x -> 1/sqrt (x^2+1)) undefined
+  acosh = minrange acosh (\x -> 1/((sqrt (x-1))*(sqrt (x+1)))) Concave
+  atanh = minrange atanh (\x -> 1/(1-x^2)) undefined
+
+type AFIndex = Int
+
+-- | AFM is a state monad that ensures that any new noise symbols have not been used by any previous affine form.
+-- All affine arithmetic calculations should be done inside the AFM monad. Affine forms do not make sense outside of their monad context.
+newtype AFMT t s m a = AFMT {runAFMT :: s -> m (a, s)}
+type AFM t a = AFMT t AFIndex Identity a
+
+instance (Monad m) => Functor (AFMT t s m) where
+  fmap = liftM
+
+instance (Monad m) => Applicative (AFMT t s m) where
+  pure = return
+  (<*>) = ap
+
+instance (Monad m) => Monad (AFMT t s m) where
+  return a = AFMT $ \s -> return (a, s)
+  (AFMT x) >>= f = AFMT $ \s -> do
+    (v, s') <- x s
+    (AFMT x') <- return $ f v
+    x' s'
+
+instance (Monad m) => MonadState s (AFMT t s m) where
+  get   = AFMT $ \s -> return (s, s)
+  put s = AFMT $ \_ -> return ((), s)
+
+instance MonadTrans (AFMT t s) where
+  lift c = AFMT $ \s -> c >>= (\x -> return (x, s))
+
+-- | This gives an affine form with midpoint 0 and radius 1.
+-- This affine form does not share epsilons with any affine forms created before it.
+-- It can be used to instantiate new affine forms.
+newEps :: Num a => AFM t (AF t a)
+newEps = do
+  idx <- get
+  put $ idx + 1
+  return $ AF 0 (replicate idx 0 ++ [1]) 0
+
+-- | Creates a new affine form that covers the interval.
+-- This affine form does not share epsilons with any affine forms created before it.
+newFromInterval :: (Eq a, Fractional a, ExplicitRounding a) => IA.Interval a -> AFM t (AF t a)
+newFromInterval i = do
+  eps <- newEps
+  let mult = ((IA.width i) / 2) .* eps
+  return $ (IA.midpoint i) .+ mult
+
+-- | Creates a new affine form that represents some exact value
+singleton :: (Num a) => a -> AF s a
+singleton x = AF x [] 0
+
+-- | Creates a new affine form that approximately represents some value.
+-- This function adds a small error to account for the 'wobble' in the computer representation of the value.
+approxSingleton :: (ExplicitRounding a) => a -> AF s a
+approxSingleton x = AF x [] $ eps x
+
+-- | Evaluates the AFM monad. It is not possible to get an AF out of an AFM monad.
+evalAFM :: forall a b. (forall t. AFM t b) -> b
+evalAFM (AFMT x) = fst . runIdentity $ x 0
+
+-- | Gives the radius of the affine form
+radius :: (Num a, ExplicitRounding a) => AF s a -> a
+radius (AF _ xs xe) = v
+  where v = xe +/ (sumup $ abs <$> xs)
+
+-- | Gives the midpoint of the affine form (the first term of the affine form).
+midpoint :: AF s a -> a
+midpoint (AF x _ _) = x
+
+-- | Gives the minimal possible value of the affine form
+inf :: (Num a, ExplicitRounding a) => AF s a -> a
+inf af = x - eps x
+  where x = (midpoint af) - (radius af)
+
+-- | Gives the maximal possible value of the affine form
+sup :: (Num a, ExplicitRounding a) => AF s a -> a
+sup af = x + eps x
+  where x = (midpoint af) + (radius af)
+
+-- | Gives the corresponding interval of the affine form
+interval :: (Num a, Ord a, ExplicitRounding a) => AF s a -> IA.Interval a
+interval af = (inf af)...(sup af)
+
+-- | Returns whether the element is representable by the affine form
+member :: (Num a, Ord a, ExplicitRounding a) => a -> AF s a -> Bool
+member x af = x `IA.member` (interval af)
+
+-- | Returns the number of noise symbols in the affine form.
+epscount_ :: AF s a -> Int
+epscount_ (AF _ xs _) = length xs
+
+-- Affine arithmetic operations
+
+-- | Sets the midpoint of the affine form
+setMidpoint :: (Num a, ExplicitRounding a) => a -> AF s a -> AF s a
+setMidpoint m (AF x xs xe) = AF m xs $ xe + eps m
+
+-- | Adds the value to the error term of the affine form
+addError :: (Num a, Ord a) => AF s a -> a -> AF s a
+addError (AF x xs xe) e
+  | e >= 0 = AF x xs (xe + e)
+  | otherwise = throw AddingNegativeError
+
+-- | Adds a scalar value to the affine form
+(.+) :: (Num a, ExplicitRounding a) => a -> AF s a -> AF s a
+a .+ (AF x xs xe) = AF m xs (xe + rnd)
+  where m = x + a
+        rnd = eps $ a + xe
+
+add :: (ExplicitRounding a, Num a, Ord a) => AF s a -> AF s a -> AF s a
+(AF x xs xe) `add` (AF y ys ye) = addError af rnd
+  where zs  = (uncurry (+)) <$> embed xs ys
+        af  = AF (x + y) zs (xe +/ ye)
+        rnd = sumup $ (uncurry (+/)) <$> embed (eps <$> xs ++ [x]) (eps <$> ys ++ [y])
+
+negateAF :: (Num a) => AF s a -> AF s a
+negateAF (AF x xs xe) = AF (-x) (negate <$> xs) xe
+
+multiply :: (Num a, Ord a, ExplicitRounding a) => AF s a -> AF s a -> AF s a
+af1@(AF x xs xe) `multiply` af2@(AF y ys ye) = addError af rnd
+  where zs = uncurry (+) <$> embed ((y*) <$> xs) ((x*) <$> ys)
+        ze1 = sum $ liftM2 (*/) (abs <$> xs ++ [xe]) (abs <$> ys ++ [ye])
+        ze2 = (abs x */ ye) +/ (abs y */ xe)
+        af = AF (x * y) zs (ze1 +/ ze2)
+        -- fig-sto-97:74
+        rnd = sumup $ (uncurry (*/)) <$> liftM2 (,) (eps <$> xs ++ [x, xe]) (eps <$> ys ++ [y, ye])
+
+-- | Multiplies the affine form by a scalar
+(.*) :: (Eq a, Num a, Ord a, ExplicitRounding a) => a -> AF s a -> AF s a
+a .* (AF x xs xe) = addError af rnd
+  where af = AF (a*x) ((a*) <$> xs) $ (a * xe)
+        rnd = sumup $ eps . (a */) <$> xs ++ [x, xe]
+
+recipAF :: (Ord a, Fractional a, ExplicitRounding a) => AF s a -> AF s a
+recipAF af
+  | low > 0   = minrange recip (\x -> -1/x^2) Convex af
+  | high < 0  = negateAF . recipAF $ negateAF af
+  | otherwise = throw DivisionByZero
+  where high = sup af
+        low  = inf af
+
+cosAF :: (Ord a, RealFrac a, Floating a, ExplicitRounding a) => AF s a -> AF s a
+cosAF af
+  | radius af < pi = f af
+  | otherwise = AF 0 [] 1
+  where a = inf af `pmod'` (2*pi)
+        b = sup af `pmod'` (2*pi)
+        f x
+          -- function never reaches extremum
+          | a < pi && b < pi || a > pi && b > pi = minrange cos (negate . sin) undefined af
+          -- function reaches extremum exactly once
+          | a < b = AF (rl - 1) [] rl
+          -- function reaches extremum more than once
+          | otherwise = AF (1 - rh) [] rh
+          where rl = abs (1 + (max (cos a) (cos b)))/2
+                rh = abs (1 - (min (cos a) (cos b)))/2
+
+sinAF :: (Ord a, RealFrac a, Floating a, ExplicitRounding a) => AF s a -> AF s a
+sinAF af = cosAF ((-pi/2) .+ af)
+
+absAF :: (Ord a, ExplicitRounding a, Fractional a) => AF s a -> AF s a
+absAF af
+  | inf af >= 0 = af
+  | sup af <= 0 = -af
+  | otherwise = AF x [] x
+    where x = (max (abs . sup $ af) (abs . inf $ af))/2
+
+signumAF :: (Ord a, Num a, ExplicitRounding a) => AF s a -> AF s a
+signumAF af
+  | inf af >= 0 = AF 1 [] 0
+  | sup af <= 0 = AF (-1) [] 0
+  | otherwise = AF 0 [] 1
+
+--
+-- Helper functions
+--
+
+-- | Fixes the epsilons of the affine form to the values in the list.
+-- The list will be padded with zeroes to match the number of coefficients.
+fix :: (Num a, Ord a, ExplicitRounding a) => AF s a -> [a] -> IA.Interval a
+fix (AF x xs xe) vals = (l)...(h)
+  where em = embed xs vals
+        s = sum $ uncurry (*) <$> em
+        m = x + s
+        l = m - xe
+        h = m + xe
+
+-- | Returns a min-range approximation function for given function and its derivative.
+minrange :: (Fractional a, Ord a, ExplicitRounding a) => (a -> a) -> (a -> a) -> Curvature -> (AF s a -> AF s a)
+minrange f f' curv = \af ->
+  let a   = sup af
+      b   = inf af
+      p   = case curv of
+              Convex  -> f' a
+              Concave -> f' b
+      q   = ((f a)+(f b)-p*(a+b))/2
+      d   = abs ((f a)-(f b)+p*(a-b))/2
+      rnd = eps $ (eps $ q +/ a */ p) + (eps $ q +/ b */ p)
+      af1 = q .+ (p .* af) `addError` (d + rnd)
+  in
+    addError af1 rnd
diff --git a/src/Numeric/AffineForm/Utils.hs b/src/Numeric/AffineForm/Utils.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AffineForm/Utils.hs
@@ -0,0 +1,33 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE DataKinds #-}
+
+-- | Provides some useful functions
+module Numeric.AffineForm.Utils (
+                                embed, sumup,
+                                pmod', clamp,
+                                ) where
+
+import Data.Fixed
+import Numeric.AffineForm.ExplicitRounding
+
+-- | Zips the two lists together, padding the shorter list with zeroes
+embed :: (Num a, Num b) => [a] -> [b] -> [(a,b)]
+embed x y = take (max (length x) (length y)) $ zip infx infy
+  where infx = x ++ repeat 0
+        infy = y ++ repeat 0
+
+
+-- | Sawtooth function with period `b`
+pmod' :: (Ord a, RealFrac a) => a -> a -> a
+pmod' a b
+  | a < 0 = pmod' (a + b*(fromIntegral . ceiling . abs $ a/b)) b
+  | otherwise = a `mod'` b
+
+-- | Clamps `x` between values `a` and `b`
+clamp :: (Ord a) => a -> a -> a -> a
+clamp x a b = min (max a x) b
+
+-- | Like sum but rounds all values up
+sumup :: (ExplicitRounding a) => [a] -> a
+sumup = foldl (+/) 0
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -0,0 +1,190 @@
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE RankNTypes #-}
+
+import Test.QuickCheck
+
+import Control.Monad
+import Text.Printf
+import Data.Fixed (mod')
+
+import qualified Numeric.Interval as IA (member, inf, sup, contains, inflate, Interval, midpoint)
+import Numeric.AffineForm.Internal
+import Numeric.AffineForm.Utils
+import Numeric.AffineForm.ExplicitRounding
+
+--
+-- Generators and modifiers
+--
+
+data EpsV a = EpsV [a]
+  deriving (Show)
+
+instance (Real a, Arbitrary a) => Arbitrary (EpsV a) where
+  arbitrary = do
+    l <- listOf $ arbitrary
+    let ls = (\x -> (x `mod'` 2) -1) <$> l
+    return $ EpsV ls
+  shrink (EpsV l) = filter validEV $ EpsV <$> (shrink l)
+
+instance (Num a, Ord a, Arbitrary a) => Arbitrary (AF s a) where
+  arbitrary = do
+                x <- arbitrary
+                xs <- arbitrary
+                (Positive xe) <- arbitrary
+                return $ AF x xs xe
+  shrink (AF x xs xe) =
+    [AF x' xs' xe' | (x', xs', xe') <- shrink (x, xs, xe)]
+
+newtype SmallExponent a = SmallExponent a
+  deriving (Show)
+
+instance (Enum a, Num a, Arbitrary a) => Arbitrary (SmallExponent a) where
+  arbitrary = SmallExponent <$> elements [1..4]
+  shrink (SmallExponent x) = SmallExponent <$> shrink x
+
+newtype ZerolessAF s a = ZerolessAF (AF s a)
+
+instance (Show a) => Show (ZerolessAF s a) where
+  show (ZerolessAF x) = show x
+
+instance (Fractional a, Ord a, Arbitrary a, ExplicitRounding a) => Arbitrary (ZerolessAF s a) where
+  arbitrary = do
+    af <- arbitrary
+    let mh = max (midpoint af) (1 + radius af)
+        ml = min (midpoint af) ((negate $ radius af) - 1)
+        res
+          | midpoint af >= 0 = ZerolessAF $ setMidpoint mh af
+          | otherwise = ZerolessAF $ setMidpoint ml af
+    return res
+
+newtype PositiveAF s a = PositiveAF (AF s a)
+
+instance (Show a) => Show (PositiveAF s a) where
+  show (PositiveAF x) = show x
+
+instance (Fractional a, Ord a, Arbitrary a, ExplicitRounding a) => Arbitrary (PositiveAF s a) where
+  arbitrary = do
+    af <- arbitrary
+    let m = 1/100000 + max (midpoint af) (radius af)
+    return . PositiveAF $ setMidpoint m af
+
+newtype SmallAF s a = SmallAF (AF s a)
+
+instance (Show a) => Show (SmallAF s a) where
+  show (SmallAF x) = show x
+
+instance (Floating a, Ord a, Arbitrary a, ExplicitRounding a) => Arbitrary (SmallAF s a) where
+  arbitrary = do
+    size <- getSize
+    af <- arbitrary
+    let s = log . fromIntegral $ size + 1
+        k = s / (radius af)
+        m = clamp (midpoint af) (-s) s
+    return . SmallAF $ setMidpoint m (k .* af)
+
+validEV :: (Ord a, Num a) => EpsV a -> Bool
+validEV (EpsV l) = all (\x -> -1 <= x && x <= 1) l
+
+--
+-- Properties
+--
+
+-- Generalized
+
+correctnessPropUnary :: (Fractional a, Ord a, Show a, ExplicitRounding a)
+  => (AF s a -> AF s a)
+  -> (a -> a)
+  -> [a]
+  -> AF s a
+  -> Property
+correctnessPropUnary f g e x = withMaxSuccess 5000 $ counterexample str res
+  where af = f x
+        rhs = g (IA.midpoint $ fix x e)
+        rhs_lo = g (IA.inf $ fix x e)
+        rhs_hi = g (IA.sup $ fix x e)
+        res = rhs `IA.member` interval af .&&.
+              rhs_lo `IA.member` interval af .&&.
+              rhs_hi `IA.member` interval af
+        str = "-- RESULTS --\n"
+           ++ "- LHS -\n"
+           ++ "AF: " ++ (show af) ++ "\n"
+           ++ "INTERVAL: " ++ (show $ interval af) ++ "\n"
+           ++ "- RHS -\n"
+           ++ "MID: " ++ (show rhs) ++ "\n"
+           ++ "HI: " ++ (show rhs_hi) ++ "\n"
+           ++ "LO: " ++ (show rhs_lo) ++ "\n"
+
+correctnessPropBinary :: (Fractional a, Ord a, Show a, ExplicitRounding a)
+  => (AF s a -> AF s a -> AF s a)
+  -> (a -> a -> a)
+  -> [a]
+  -> AF s a
+  -> AF s a
+  -> Property
+correctnessPropBinary f g e x y = withMaxSuccess 5000 $ counterexample str res
+  where af = f x y
+        rhs = g (IA.midpoint $ fix x e) (IA.midpoint $ fix y e)
+        rhs_lo = g (IA.inf $ fix x e) (IA.inf $ fix y e)
+        rhs_hi = g (IA.sup $ fix x e) (IA.sup $ fix y e)
+
+        res = rhs `IA.member` interval af .&&.
+              rhs_lo `IA.member` interval af .&&.
+              rhs_hi `IA.member` interval af
+        str = "-- RESULTS --\n"
+           ++ "- LHS -\n"
+           ++ "AF: " ++ (show af) ++ "\n"
+           ++ "INTERVAL: " ++ (show $ interval af) ++ "\n"
+           ++ "- RHS -\n"
+           ++ "MID: " ++ (show rhs) ++ "\n"
+           ++ "HI: " ++ (show rhs_hi) ++ "\n"
+           ++ "LO: " ++ (show rhs_lo) ++ "\n"
+
+-- RuKaS14.pdf [1102:2]
+-- prop_addition :: EpsV Double -> AF Double -> AF Double -> Property
+-- prop_addition (EpsV e) x y = counterexample str res
+--   where lhs = (x + y) `fix` e
+--         rhs = x `fix` e + y `fix` e
+--         res = lhs `IA.contains` rhs
+--         str = "AA: " ++ (show lhs) ++ "\nIA: " ++ (show rhs)
+
+prop_sound_addition :: EpsV Double -> AF s Double -> AF s Double -> Property
+prop_sound_addition (EpsV e) x y = correctnessPropBinary (+) (+) e x y
+
+prop_sound_subtraction :: EpsV Double -> AF s Double -> AF s Double -> Property
+prop_sound_subtraction (EpsV e) x y = correctnessPropBinary (-) (-) e x y
+
+prop_sound_multiplication :: EpsV Double -> AF s Double -> AF s Double -> Property
+prop_sound_multiplication (EpsV e) x y = correctnessPropBinary (*) (*) e x y
+
+prop_sound_power :: EpsV Double -> AF s Double -> SmallExponent Integer -> Property
+prop_sound_power (EpsV e) x (SmallExponent n) = correctnessPropUnary (^n) (^n) e x
+
+prop_sound_recip :: EpsV Double -> ZerolessAF s Double -> Property
+prop_sound_recip (EpsV e) (ZerolessAF x) = correctnessPropUnary recip recip e x
+
+prop_sound_log :: EpsV Double -> PositiveAF s Double -> Property
+prop_sound_log (EpsV e) (PositiveAF x) = correctnessPropUnary log log e x
+
+prop_sound_exp :: EpsV Double -> SmallAF s Double -> Property
+prop_sound_exp (EpsV e) (SmallAF x) = correctnessPropUnary exp exp e x
+
+prop_sound_abs :: EpsV Double -> AF s Double -> Property
+prop_sound_abs (EpsV e) x = correctnessPropUnary abs abs e x
+
+-- prop_sin :: EpsV Double -> AF s Double -> Property
+-- prop_sin (EpsV e) x = correctnessPropUnary sin sin e x
+
+-- prop_cos :: EpsV Double -> AF s Double -> Property
+-- prop_cos (EpsV e) x = correctnessPropUnary cos cos e x
+
+--
+-- Testing boilerplate
+--
+
+return [] -- This is a hack to make the quickCheckAll template work correctly
+
+main :: IO ()
+main = do
+  $quickCheckAll
+  return ()
+