diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -1,3 +1,11 @@
+0.5
+---
+* The default `Numeric.Interval` now deals more conventionally with empty intervals.
+* The old "Kaucher directed interval" behavior is available as `Numeric.Interval.Kaucher`.
+* Strictly Non-Empty intervals are now contained in `Numeric.Interval.NonEmpty`
+* Renamed `bisection` to `bisect`.
+* Added `bisectIntegral`.
+
 0.4.2
 -----
 * Added `clamp`
diff --git a/intervals.cabal b/intervals.cabal
--- a/intervals.cabal
+++ b/intervals.cabal
@@ -1,5 +1,5 @@
 name:              intervals
-version:           0.4.2
+version:           0.5
 synopsis:          Interval Arithmetic
 description:
   A 'Numeric.Interval.Interval' is a closed, convex set of floating point values.
@@ -38,7 +38,12 @@
 library
   hs-source-dirs: src
 
-  exposed-modules: Numeric.Interval
+  exposed-modules:
+    Numeric.Interval
+    Numeric.Interval.Internal
+    Numeric.Interval.Kaucher
+    Numeric.Interval.NonEmpty
+    Numeric.Interval.NonEmpty.Internal
 
   build-depends:
     array          >= 0.3   && < 0.6,
diff --git a/src/Numeric/Interval.hs b/src/Numeric/Interval.hs
--- a/src/Numeric/Interval.hs
+++ b/src/Numeric/Interval.hs
@@ -1,24 +1,16 @@
-{-# LANGUAGE CPP #-}
-{-# LANGUAGE Rank2Types #-}
-{-# LANGUAGE DeriveDataTypeable #-}
-#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
-{-# LANGUAGE DeriveGeneric #-}
-#endif
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Interval
--- Copyright   :  (c) Edward Kmett 2010-2013
+-- Copyright   :  (c) Edward Kmett 2010-2014
 -- License     :  BSD3
 -- Maintainer  :  ekmett@gmail.com
 -- Stability   :  experimental
 -- Portability :  DeriveDataTypeable
 --
 -- Interval arithmetic
---
 -----------------------------------------------------------------------------
-
 module Numeric.Interval
-  ( Interval(..)
+  ( Interval
   , (...)
   , whole
   , empty
@@ -33,717 +25,17 @@
   , midpoint
   , intersection
   , hull
-  , bisection
+  , bisect
+  , bisectIntegral
   , magnitude
   , mignitude
   , contains
   , isSubsetOf
   , certainly, (<!), (<=!), (==!), (>=!), (>!)
   , possibly, (<?), (<=?), (==?), (>=?), (>?)
-  , clamp
   , idouble
   , ifloat
   ) where
 
-import Control.Applicative hiding (empty)
-import Data.Data
-import Data.Distributive
-import Data.Foldable hiding (minimum, maximum, elem, notElem)
-import Data.Function (on)
-import Data.Monoid
-import Data.Traversable
-#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
-import GHC.Generics
-#endif
-import Prelude hiding (null, elem, notElem)
-
--- $setup
-
-data Interval a = I !a !a deriving
-  ( Data
-  , Typeable
-#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
-  , Generic
-#if __GLASGOW_HASKELL__ >= 706
-  , Generic1
-#endif
-#endif
-  )
-
-instance Functor Interval where
-  fmap f (I a b) = I (f a) (f b)
-  {-# INLINE fmap #-}
-
-instance Foldable Interval where
-  foldMap f (I a b) = f a `mappend` f b
-  {-# INLINE foldMap #-}
-
-instance Traversable Interval where
-  traverse f (I a b) = I <$> f a <*> f b
-  {-# INLINE traverse #-}
-
-instance Applicative Interval where
-  pure a = I a a
-  {-# INLINE pure #-}
-  I f g <*> I a b = I (f a) (g b)
-  {-# INLINE (<*>) #-}
-
-instance Monad Interval where
-  return a = I a a
-  {-# INLINE return #-}
-  I a b >>= f = I a' b' where
-    I a' _ = f a
-    I _ b' = f b
-  {-# INLINE (>>=) #-}
-
-instance Distributive Interval where
-  distribute f = fmap inf f ... fmap sup f
-  {-# INLINE distribute #-}
-
-infix 3 ...
-
-negInfinity :: Fractional a => a
-negInfinity = (-1)/0
-{-# INLINE negInfinity #-}
-
-posInfinity :: Fractional a => a
-posInfinity = 1/0
-{-# INLINE posInfinity #-}
-
-nan :: Fractional a => a
-nan = 0/0
-
-fmod :: RealFrac a => a -> a -> a
-fmod a b = a - q*b where
-  q = realToFrac (truncate $ a / b :: Integer)
-{-# INLINE fmod #-}
-
--- | The rule of thumb is you should only use this to construct using values
--- that you took out of the interval. Otherwise, use I, to force rounding
-(...) :: a -> a -> Interval a
-(...) = I
-{-# INLINE (...) #-}
-
--- | The whole real number line
---
--- >>> whole
--- -Infinity ... Infinity
-whole :: Fractional a => Interval a
-whole = negInfinity ... posInfinity
-{-# INLINE whole #-}
-
--- | An empty interval
---
--- >>> empty
--- NaN ... NaN
-empty :: Fractional a => Interval a
-empty = nan ... nan
-{-# INLINE empty #-}
-
--- | negation handles NaN properly
---
--- >>> null (1 ... 5)
--- False
---
--- >>> null (1 ... 1)
--- False
---
--- >>> null empty
--- True
-null :: Ord a => Interval a -> Bool
-null x = not (inf x <= sup x)
-{-# INLINE null #-}
-
--- | A singleton point
---
--- >>> singleton 1
--- 1 ... 1
-singleton :: a -> Interval a
-singleton a = a ... a
-{-# INLINE singleton #-}
-
--- | The infinumum (lower bound) of an interval
---
--- >>> inf (1 ... 20)
--- 1
-inf :: Interval a -> a
-inf (I a _) = a
-{-# INLINE inf #-}
-
--- | The supremum (upper bound) of an interval
---
--- >>> sup (1 ... 20)
--- 20
-sup :: Interval a -> a
-sup (I _ b) = b
-{-# INLINE sup #-}
-
--- | Is the interval a singleton point?
--- N.B. This is fairly fragile and likely will not hold after
--- even a few operations that only involve singletons
---
--- >>> singular (singleton 1)
--- True
---
--- >>> singular (1.0 ... 20.0)
--- False
-singular :: Ord a => Interval a -> Bool
-singular x = not (null x) && inf x == sup x
-{-# INLINE singular #-}
-
-instance Eq a => Eq (Interval a) where
-  (==) = (==!)
-  {-# INLINE (==) #-}
-
-instance Show a => Show (Interval a) where
-  showsPrec n (I a b) =
-    showParen (n > 3) $
-      showsPrec 3 a .
-      showString " ... " .
-      showsPrec 3 b
-
--- | Calculate the width of an interval.
---
--- >>> width (1 ... 20)
--- 19
---
--- >>> width (singleton 1)
--- 0
---
--- >>> width empty
--- NaN
-width :: Num a => Interval a -> a
-width (I a b) = b - a
-{-# INLINE width #-}
-
--- | Magnitude
---
--- >>> magnitude (1 ... 20)
--- 20
---
--- >>> magnitude (-20 ... 10)
--- 20
---
--- >>> magnitude (singleton 5)
--- 5
-magnitude :: (Num a, Ord a) => Interval a -> a
-magnitude x = (max `on` abs) (inf x) (sup x)
-{-# INLINE magnitude #-}
-
--- | \"mignitude\"
---
--- >>> mignitude (1 ... 20)
--- 1
---
--- >>> mignitude (-20 ... 10)
--- 10
---
--- >>> mignitude (singleton 5)
--- 5
-mignitude :: (Num a, Ord a) => Interval a -> a
-mignitude x = (min `on` abs) (inf x) (sup x)
-{-# INLINE mignitude #-}
-
-instance (Num a, Ord a) => Num (Interval a) where
-  I a b + I a' b' = (a + a') ... (b + b')
-  {-# INLINE (+) #-}
-  I a b - I a' b' = (a - b') ... (b - a')
-  {-# INLINE (-) #-}
-  I a b * I a' b' =
-    minimum [a * a', a * b', b * a', b * b']
-    ...
-    maximum [a * a', a * b', b * a', b * b']
-  {-# INLINE (*) #-}
-  abs x@(I a b)
-    | a >= 0    = x
-    | b <= 0    = negate x
-    | otherwise = 0 ... max (- a) b
-  {-# INLINE abs #-}
-
-  signum = increasing signum
-  {-# INLINE signum #-}
-
-  fromInteger i = singleton (fromInteger i)
-  {-# INLINE fromInteger #-}
-
--- | Bisect an interval at its midpoint.
---
--- >>> bisection (10.0 ... 20.0)
--- (10.0 ... 15.0,15.0 ... 20.0)
---
--- >>> bisection (singleton 5.0)
--- (5.0 ... 5.0,5.0 ... 5.0)
---
--- >>> bisection empty
--- (NaN ... NaN,NaN ... NaN)
-bisection :: Fractional a => Interval a -> (Interval a, Interval a)
-bisection x = (inf x ... m, m ... sup x)
-  where m = midpoint x
-{-# INLINE bisection #-}
-
--- | Nearest point to the midpoint of the interval.
---
--- >>> midpoint (10.0 ... 20.0)
--- 15.0
---
--- >>> midpoint (singleton 5.0)
--- 5.0
---
--- >>> midpoint empty
--- NaN
-midpoint :: Fractional a => Interval a -> a
-midpoint x = inf x + (sup x - inf x) / 2
-{-# INLINE midpoint #-}
-
--- | Determine if a point is in the interval.
---
--- >>> elem 3.2 (1.0 ... 5.0)
--- True
---
--- >>> elem 5 (1.0 ... 5.0)
--- True
---
--- >>> elem 1 (1.0 ... 5.0)
--- True
---
--- >>> elem 8 (1.0 ... 5.0)
--- False
---
--- >>> elem 5 empty
--- False
---
-elem :: Ord a => a -> Interval a -> Bool
-elem x xs = x >= inf xs && x <= sup xs
-{-# INLINE elem #-}
-
--- | Determine if a point is not included in the interval
---
--- >>> notElem 8 (1.0 ... 5.0)
--- True
---
--- >>> notElem 1.4 (1.0 ... 5.0)
--- False
---
--- And of course, nothing is a member of the empty interval.
---
--- >>> notElem 5 empty
--- True
-notElem :: Ord a => a -> Interval a -> Bool
-notElem x xs = not (elem x xs)
-{-# INLINE notElem #-}
-
--- | 'realToFrac' will use the midpoint
-instance Real a => Real (Interval a) where
-  toRational x
-    | null x   = nan
-    | otherwise = a + (b - a) / 2
-    where
-      a = toRational (inf x)
-      b = toRational (sup x)
-  {-# INLINE toRational #-}
-
-instance Ord a => Ord (Interval a) where
-  compare x y
-    | sup x < inf y = LT
-    | inf x > sup y = GT
-    | sup x == inf y && inf x == sup y = EQ
-    | otherwise = error "Numeric.Interval.compare: ambiguous comparison"
-  {-# INLINE compare #-}
-
-  max (I a b) (I a' b') = max a a' ... max b b'
-  {-# INLINE max #-}
-
-  min (I a b) (I a' b') = min a a' ... min b b'
-  {-# INLINE min #-}
-
--- @'divNonZero' X Y@ assumes @0 `'notElem'` Y@
-divNonZero :: (Fractional a, Ord a) => Interval a -> Interval a -> Interval a
-divNonZero (I a b) (I a' b') =
-  minimum [a / a', a / b', b / a', b / b']
-  ...
-  maximum [a / a', a / b', b / a', b / b']
-
--- @'divPositive' X y@ assumes y > 0, and divides @X@ by [0 ... y]
-divPositive :: (Fractional a, Ord a) => Interval a -> a -> Interval a
-divPositive x@(I a b) y
-  | a == 0 && b == 0 = x
-  -- b < 0 || isNegativeZero b = negInfinity ... ( b / y)
-  | b < 0 = negInfinity ... ( b / y)
-  | a < 0 = whole
-  | otherwise = (a / y) ... posInfinity
-{-# INLINE divPositive #-}
-
--- divNegative assumes y < 0 and divides the interval @X@ by [y ... 0]
-divNegative :: (Fractional a, Ord a) => Interval a -> a -> Interval a
-divNegative x@(I a b) y
-  | a == 0 && b == 0 = - x -- flip negative zeros
-  -- b < 0 || isNegativeZero b = (b / y) ... posInfinity
-  | b < 0 = (b / y) ... posInfinity
-  | a < 0 = whole
-  | otherwise = negInfinity ... (a / y)
-{-# INLINE divNegative #-}
-
-divZero :: (Fractional a, Ord a) => Interval a -> Interval a
-divZero x
-  | inf x == 0 && sup x == 0 = x
-  | otherwise = whole
-{-# INLINE divZero #-}
-
-instance (Fractional a, Ord a) => Fractional (Interval a) where
-  -- TODO: check isNegativeZero properly
-  x / y
-    | 0 `notElem` y = divNonZero x y
-    | iz && sz  = empty -- division by 0
-    | iz        = divPositive x (inf y)
-    |       sz  = divNegative x (sup y)
-    | otherwise = divZero x
-    where
-      iz = inf y == 0
-      sz = sup y == 0
-  recip (I a b)   = on min recip a b ... on max recip a b
-  {-# INLINE recip #-}
-  fromRational r  = let r' = fromRational r in r' ... r'
-  {-# INLINE fromRational #-}
-
-instance RealFrac a => RealFrac (Interval a) where
-  properFraction x = (b, x - fromIntegral b)
-    where
-      b = truncate (midpoint x)
-  {-# INLINE properFraction #-}
-  ceiling x = ceiling (sup x)
-  {-# INLINE ceiling #-}
-  floor x = floor (inf x)
-  {-# INLINE floor #-}
-  round x = round (midpoint x)
-  {-# INLINE round #-}
-  truncate x = truncate (midpoint x)
-  {-# INLINE truncate #-}
-
-instance (RealFloat a, Ord a) => Floating (Interval a) where
-  pi = singleton pi
-  {-# INLINE pi #-}
-  exp = increasing exp
-  {-# INLINE exp #-}
-  log (I a b) = (if a > 0 then log a else negInfinity) ... log b
-  {-# INLINE log #-}
-  cos x
-    | null x = empty
-    | width t >= pi = (-1) ... 1
-    | inf t >= pi = - cos (t - pi)
-    | sup t <= pi = decreasing cos t
-    | sup t <= 2 * pi = (-1) ... cos ((pi * 2 - sup t) `min` inf t)
-    | otherwise = (-1) ... 1
-    where
-      t = fmod x (pi * 2)
-  {-# INLINE cos #-}
-  sin x
-    | null x = empty
-    | otherwise = cos (x - pi / 2)
-  {-# INLINE sin #-}
-  tan x
-    | null x = empty
-    | inf t' <= - pi / 2 || sup t' >= pi / 2 = whole
-    | otherwise = increasing tan x
-    where
-      t = x `fmod` pi
-      t' | t >= pi / 2 = t - pi
-         | otherwise    = t
-  {-# INLINE tan #-}
-  asin x@(I a b)
-    | null x || b < -1 || a > 1 = empty
-    | otherwise =
-      (if a <= -1 then -halfPi else asin a)
-      ...
-      (if b >= 1 then halfPi else asin b)
-    where
-      halfPi = pi / 2
-  {-# INLINE asin #-}
-  acos x@(I a b)
-    | null x || b < -1 || a > 1 = empty
-    | otherwise =
-      (if b >= 1 then 0 else acos b)
-      ...
-      (if a < -1 then pi else acos a)
-  {-# INLINE acos #-}
-  atan = increasing atan
-  {-# INLINE atan #-}
-  sinh = increasing sinh
-  {-# INLINE sinh #-}
-  cosh x@(I a b)
-    | null x = empty
-    | b < 0  = decreasing cosh x
-    | a >= 0 = increasing cosh x
-    | otherwise  = I 0 $ cosh $ if - a > b
-                                then a
-                                else b
-  {-# INLINE cosh #-}
-  tanh = increasing tanh
-  {-# INLINE tanh #-}
-  asinh = increasing asinh
-  {-# INLINE asinh #-}
-  acosh x@(I a b)
-    | null x || b < 1 = empty
-    | otherwise = I lo $ acosh b
-    where lo | a <= 1 = 0
-             | otherwise = acosh a
-  {-# INLINE acosh #-}
-  atanh x@(I a b)
-    | null x || b < -1 || a > 1 = empty
-    | otherwise =
-      (if a <= - 1 then negInfinity else atanh a)
-      ...
-      (if b >= 1 then posInfinity else atanh b)
-  {-# INLINE atanh #-}
-
--- | lift a monotone increasing function over a given interval
-increasing :: (a -> b) -> Interval a -> Interval b
-increasing f (I a b) = f a ... f b
-
--- | lift a monotone decreasing function over a given interval
-decreasing :: (a -> b) -> Interval a -> Interval b
-decreasing f (I a b) = f b ... f a
-
--- | We have to play some semantic games to make these methods make sense.
--- Most compute with the midpoint of the interval.
-instance RealFloat a => RealFloat (Interval a) where
-  floatRadix = floatRadix . midpoint
-
-  floatDigits = floatDigits . midpoint
-  floatRange = floatRange . midpoint
-  decodeFloat = decodeFloat . midpoint
-  encodeFloat m e = singleton (encodeFloat m e)
-  exponent = exponent . midpoint
-  significand x = min a b ... max a b
-    where
-      (_ ,em) = decodeFloat (midpoint x)
-      (mi,ei) = decodeFloat (inf x)
-      (ms,es) = decodeFloat (sup x)
-      a = encodeFloat mi (ei - em - floatDigits x)
-      b = encodeFloat ms (es - em - floatDigits x)
-  scaleFloat n x = scaleFloat n (inf x) ... scaleFloat n (sup x)
-  isNaN x = isNaN (inf x) || isNaN (sup x)
-  isInfinite x = isInfinite (inf x) || isInfinite (sup x)
-  isDenormalized x = isDenormalized (inf x) || isDenormalized (sup x)
-  -- contains negative zero
-  isNegativeZero x = not (inf x > 0)
-                  && not (sup x < 0)
-                  && (  (sup x == 0 && (inf x < 0 || isNegativeZero (inf x)))
-                     || (inf x == 0 && isNegativeZero (inf x))
-                     || (inf x < 0 && sup x >= 0))
-  isIEEE x = isIEEE (inf x) && isIEEE (sup x)
-  atan2 = error "unimplemented"
-
--- TODO: (^), (^^) to give tighter bounds
-
--- | Calculate the intersection of two intervals.
---
--- >>> intersection (1 ... 10 :: Interval Double) (5 ... 15 :: Interval Double)
--- 5.0 ... 10.0
-intersection :: (Fractional a, Ord a) => Interval a -> Interval a -> Interval a
-intersection x@(I a b) y@(I a' b')
-  | x /=! y = empty
-  | otherwise = max a a' ... min b b'
-{-# INLINE intersection #-}
-
--- | Calculate the convex hull of two intervals
---
--- >>> hull (0 ... 10 :: Interval Double) (5 ... 15 :: Interval Double)
--- 0.0 ... 15.0
---
--- >>> hull (15 ... 85 :: Interval Double) (0 ... 10 :: Interval Double)
--- 0.0 ... 85.0
-hull :: Ord a => Interval a -> Interval a -> Interval a
-hull x@(I a b) y@(I a' b')
-  | null x = y
-  | null y = x
-  | otherwise = min a a' ... max b b'
-{-# INLINE hull #-}
-
--- | For all @x@ in @X@, @y@ in @Y@. @x '<' y@
---
--- >>> (5 ... 10 :: Interval Double) <! (20 ... 30 :: Interval Double)
--- True
---
--- >>> (5 ... 10 :: Interval Double) <! (10 ... 30 :: Interval Double)
--- False
---
--- >>> (20 ... 30 :: Interval Double) <! (5 ... 10 :: Interval Double)
--- False
-(<!)  :: Ord a => Interval a -> Interval a -> Bool
-x <! y = sup x < inf y
-{-# INLINE (<!) #-}
-
--- | For all @x@ in @X@, @y@ in @Y@. @x '<=' y@
---
--- >>> (5 ... 10 :: Interval Double) <=! (20 ... 30 :: Interval Double)
--- True
---
--- >>> (5 ... 10 :: Interval Double) <=! (10 ... 30 :: Interval Double)
--- True
---
--- >>> (20 ... 30 :: Interval Double) <=! (5 ... 10 :: Interval Double)
--- False
-(<=!) :: Ord a => Interval a -> Interval a -> Bool
-x <=! y = sup x <= inf y
-{-# INLINE (<=!) #-}
-
--- | For all @x@ in @X@, @y@ in @Y@. @x '==' y@
---
--- Only singleton intervals return true
---
--- >>> (singleton 5 :: Interval Double) ==! (singleton 5 :: Interval Double)
--- True
---
--- >>> (5 ... 10 :: Interval Double) ==! (5 ... 10 :: Interval Double)
--- False
-(==!) :: Eq a => Interval a -> Interval a -> Bool
-x ==! y = sup x == inf y && inf x == sup y
-{-# INLINE (==!) #-}
-
--- | For all @x@ in @X@, @y@ in @Y@. @x '/=' y@
---
--- >>> (5 ... 15 :: Interval Double) /=! (20 ... 40 :: Interval Double)
--- True
---
--- >>> (5 ... 15 :: Interval Double) /=! (15 ... 40 :: Interval Double)
--- False
-(/=!) :: Ord a => Interval a -> Interval a -> Bool
-x /=! y = sup x < inf y || inf x > sup y
-{-# INLINE (/=!) #-}
-
--- | For all @x@ in @X@, @y@ in @Y@. @x '>' y@
---
--- >>> (20 ... 40 :: Interval Double) >! (10 ... 19 :: Interval Double)
--- True
---
--- >>> (5 ... 20 :: Interval Double) >! (15 ... 40 :: Interval Double)
--- False
-(>!)  :: Ord a => Interval a -> Interval a -> Bool
-x >! y = inf x > sup y
-{-# INLINE (>!) #-}
-
--- | For all @x@ in @X@, @y@ in @Y@. @x '>=' y@
---
--- >>> (20 ... 40 :: Interval Double) >=! (10 ... 20 :: Interval Double)
--- True
---
--- >>> (5 ... 20 :: Interval Double) >=! (15 ... 40 :: Interval Double)
--- False
-(>=!) :: Ord a => Interval a -> Interval a -> Bool
-x >=! y = inf x >= sup y
-{-# INLINE (>=!) #-}
-
--- | For all @x@ in @X@, @y@ in @Y@. @x `op` y@
---
---
-certainly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool
-certainly cmp l r
-    | lt && eq && gt = True
-    | lt && eq       = l <=! r
-    | lt &&       gt = l /=! r
-    | lt             = l <!  r
-    |       eq && gt = l >=! r
-    |       eq       = l ==! r
-    |             gt = l >!  r
-    | otherwise      = False
-    where
-        lt = cmp LT EQ
-        eq = cmp EQ EQ
-        gt = cmp GT EQ
-{-# INLINE certainly #-}
-
--- | Check if interval @X@ totally contains interval @Y@
---
--- >>> (20 ... 40 :: Interval Double) `contains` (25 ... 35 :: Interval Double)
--- True
---
--- >>> (20 ... 40 :: Interval Double) `contains` (15 ... 35 :: Interval Double)
--- False
-contains :: Ord a => Interval a -> Interval a -> Bool
-contains x y = null y
-            || (not (null x) && inf x <= inf y && sup y <= sup x)
-{-# INLINE contains #-}
-
--- | Flipped version of `contains`. Check if interval @X@ a subset of interval @Y@
---
--- >>> (25 ... 35 :: Interval Double) `isSubsetOf` (20 ... 40 :: Interval Double)
--- True
---
--- >>> (20 ... 40 :: Interval Double) `isSubsetOf` (15 ... 35 :: Interval Double)
--- False
-isSubsetOf :: Ord a => Interval a -> Interval a -> Bool
-isSubsetOf = flip contains
-{-# INLINE isSubsetOf #-}
-
--- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<' y@?
-(<?) :: Ord a => Interval a -> Interval a -> Bool
-x <? y = inf x < sup y
-{-# INLINE (<?) #-}
-
--- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<=' y@?
-(<=?) :: Ord a => Interval a -> Interval a -> Bool
-x <=? y = inf x <= sup y
-{-# INLINE (<=?) #-}
-
--- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '==' y@?
-(==?) :: Ord a => Interval a -> Interval a -> Bool
-x ==? y = inf x <= sup y && sup x >= inf y
-{-# INLINE (==?) #-}
-
--- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '/=' y@?
-(/=?) :: Eq a => Interval a -> Interval a -> Bool
-x /=? y = inf x /= sup y || sup x /= inf y
-{-# INLINE (/=?) #-}
-
--- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>' y@?
-(>?) :: Ord a => Interval a -> Interval a -> Bool
-x >? y = sup x > inf y
-{-# INLINE (>?) #-}
-
--- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>=' y@?
-(>=?) :: Ord a => Interval a -> Interval a -> Bool
-x >=? y = sup x >= inf y
-{-# INLINE (>=?) #-}
-
--- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x `op` y@?
-possibly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool
-possibly cmp l r
-    | lt && eq && gt = True
-    | lt && eq       = l <=? r
-    | lt &&       gt = l /=? r
-    | lt             = l <? r
-    |       eq && gt = l >=? r
-    |       eq       = l ==? r
-    |             gt = l >? r
-    | otherwise      = False
-    where
-        lt = cmp LT EQ
-        eq = cmp EQ EQ
-        gt = cmp GT EQ
-{-# INLINE possibly #-}
-
--- | The nearest value to that supplied which is contained in the interval.
-clamp :: Ord a => Interval a -> a -> a
-clamp (I a b) x | x < a     = a
-                | x > b     = b
-                | otherwise = x
-
--- | id function. Useful for type specification
---
--- >>> :t idouble (1 ... 3)
--- idouble (1 ... 3) :: Interval Double
-idouble :: Interval Double -> Interval Double
-idouble = id
-
--- | id function. Useful for type specification
---
--- >>> :t ifloat (1 ... 3)
--- ifloat (1 ... 3) :: Interval Float
-ifloat :: Interval Float -> Interval Float
-ifloat = id
-
--- Bugs:
--- sin 1 :: Interval Double
-
-
-default (Integer,Double)
+import Numeric.Interval.Internal
+import Prelude ()
diff --git a/src/Numeric/Interval/Internal.hs b/src/Numeric/Interval/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Interval/Internal.hs
@@ -0,0 +1,783 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
+{-# LANGUAGE DeriveGeneric #-}
+#endif
+{-# OPTIONS_HADDOCK not-home #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Interval.Internal
+-- Copyright   :  (c) Edward Kmett 2010-2013
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  DeriveDataTypeable
+--
+-- Interval arithmetic
+--
+-----------------------------------------------------------------------------
+module Numeric.Interval.Internal
+  ( Interval(..)
+  , (...)
+  , whole
+  , empty
+  , null
+  , singleton
+  , elem
+  , notElem
+  , inf
+  , sup
+  , singular
+  , width
+  , midpoint
+  , intersection
+  , hull
+  , bisect
+  , bisectIntegral
+  , magnitude
+  , mignitude
+  , contains
+  , isSubsetOf
+  , certainly, (<!), (<=!), (==!), (>=!), (>!)
+  , possibly, (<?), (<=?), (==?), (>=?), (>?)
+  , idouble
+  , ifloat
+  ) where
+
+import Data.Data
+import Data.Foldable hiding (minimum, maximum, elem, notElem)
+import Data.Function (on)
+import Data.Monoid
+#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
+import GHC.Generics
+#endif
+import Prelude hiding (null, elem, notElem)
+
+-- $setup
+
+data Interval a = I !a !a | Empty deriving
+  ( Data
+  , Typeable
+#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
+  , Generic
+#if __GLASGOW_HASKELL__ >= 706
+  , Generic1
+#endif
+#endif
+  )
+
+instance Foldable Interval where
+  foldMap f (I a b) = f a `mappend` f b
+  foldMap _ Empty = mempty
+  {-# INLINE foldMap #-}
+
+infix 3 ...
+
+negInfinity :: Fractional a => a
+negInfinity = (-1)/0
+{-# INLINE negInfinity #-}
+
+posInfinity :: Fractional a => a
+posInfinity = 1/0
+{-# INLINE posInfinity #-}
+
+nan :: Fractional a => a
+nan = 0/0
+
+fmod :: RealFrac a => a -> a -> a
+fmod a b = a - q*b where
+  q = realToFrac (truncate $ a / b :: Integer)
+{-# INLINE fmod #-}
+
+(...) :: Ord a => a -> a -> Interval a
+a ... b
+  | a <= b = I a b
+  | otherwise = Empty
+{-# INLINE (...) #-}
+
+-- | The whole real number line
+--
+-- >>> whole
+-- -Infinity ... Infinity
+whole :: Fractional a => Interval a
+whole = I negInfinity posInfinity
+{-# INLINE whole #-}
+
+-- | An empty interval
+--
+-- >>> empty
+-- Empty
+empty :: Fractional a => Interval a
+empty = Empty
+{-# INLINE empty #-}
+
+-- | negation handles NaN properly
+--
+-- >>> null (1 ... 5)
+-- False
+--
+-- >>> null (1 ... 1)
+-- False
+--
+-- >>> null empty
+-- True
+null :: Interval a -> Bool
+null Empty = True
+null _ = False
+{-# INLINE null #-}
+
+-- | A singleton point
+--
+-- >>> singleton 1
+-- 1 ... 1
+singleton :: a -> Interval a
+singleton a = I a a
+{-# INLINE singleton #-}
+
+-- | The infinumum (lower bound) of an interval
+--
+-- >>> inf (1 ... 20)
+-- 1
+inf :: Fractional a => Interval a -> a
+inf (I a _) = a
+inf Empty = nan
+{-# INLINE inf #-}
+
+-- | The supremum (upper bound) of an interval
+--
+-- >>> sup (1 ... 20)
+-- 20
+sup :: Fractional a => Interval a -> a
+sup (I _ b) = b
+sup Empty = nan
+{-# INLINE sup #-}
+
+-- | Is the interval a singleton point?
+-- N.B. This is fairly fragile and likely will not hold after
+-- even a few operations that only involve singletons
+--
+-- >>> singular (singleton 1)
+-- True
+--
+-- >>> singular (1.0 ... 20.0)
+-- False
+singular :: Ord a => Interval a -> Bool
+singular Empty = False
+singular (I a b) = a == b
+{-# INLINE singular #-}
+
+instance Eq a => Eq (Interval a) where
+  (==) = (==!)
+  {-# INLINE (==) #-}
+
+instance Show a => Show (Interval a) where
+  showsPrec _ Empty = showString "Empty"
+  showsPrec n (I a b) =
+    showParen (n > 3) $
+      showsPrec 3 a .
+      showString " ... " .
+      showsPrec 3 b
+
+-- | Calculate the width of an interval.
+--
+-- >>> width (1 ... 20)
+-- 19
+--
+-- >>> width (singleton 1)
+-- 0
+--
+-- >>> width empty
+-- 0
+width :: Num a => Interval a -> a
+width (I a b) = b - a
+width Empty = 0
+{-# INLINE width #-}
+
+-- | Magnitude
+--
+-- >>> magnitude (1 ... 20)
+-- 20
+--
+-- >>> magnitude (-20 ... 10)
+-- 20
+--
+-- >>> magnitude (singleton 5)
+-- 5
+--
+-- >>> magnitude empty
+-- 0
+magnitude :: (Num a, Ord a) => Interval a -> a
+magnitude (I a b) = on max abs a b
+magnitude Empty = 0
+{-# INLINE magnitude #-}
+
+-- | \"mignitude\"
+--
+-- >>> mignitude (1 ... 20)
+-- 1
+--
+-- >>> mignitude (-20 ... 10)
+-- 10
+--
+-- >>> mignitude (singleton 5)
+-- 5
+--
+-- >>> mignitude empty
+-- 0
+mignitude :: (Num a, Ord a) => Interval a -> a
+mignitude (I a b) = on min abs a b
+mignitude Empty = 0
+{-# INLINE mignitude #-}
+
+instance (Num a, Ord a) => Num (Interval a) where
+  I a b + I a' b' = (a + a') ... (b + b')
+  _ + _ = Empty
+  {-# INLINE (+) #-}
+  I a b - I a' b' = (a - b') ... (b - a')
+  _ - _ = Empty
+  {-# INLINE (-) #-}
+  I a b * I a' b' =
+    minimum [a * a', a * b', b * a', b * b']
+    ...
+    maximum [a * a', a * b', b * a', b * b']
+  _ * _ = Empty
+  {-# INLINE (*) #-}
+  abs x@(I a b)
+    | a >= 0    = x
+    | b <= 0    = negate x
+    | otherwise = 0 ... max (- a) b
+  abs Empty = Empty
+  {-# INLINE abs #-}
+
+  signum = increasing signum
+  {-# INLINE signum #-}
+
+  fromInteger i = singleton (fromInteger i)
+  {-# INLINE fromInteger #-}
+
+-- | Bisect an interval at its midpoint.
+--
+-- >>> bisect (10.0 ... 20.0)
+-- (10.0 ... 15.0,15.0 ... 20.0)
+--
+-- >>> bisect (singleton 5.0)
+-- (5.0 ... 5.0,5.0 ... 5.0)
+--
+-- >>> bisect empty
+-- (NaN ... NaN,NaN ... NaN)
+bisect :: Fractional a => Interval a -> (Interval a, Interval a)
+bisect Empty = (Empty,Empty)
+bisect (I a b) = (I a m, I m b) where m = a + (b - a) / 2
+{-# INLINE bisect #-}
+
+bisectIntegral :: Integral a => Interval a -> (Interval a, Interval a)
+bisectIntegral Empty = (Empty, Empty)
+bisectIntegral (I a b)
+  | a == m || b == m = (I a a, I b b)
+  | otherwise        = (I a m, I m b)
+  where m = a + (b - a) `div` 2
+
+-- | Nearest point to the midpoint of the interval.
+--
+-- >>> midpoint (10.0 ... 20.0)
+-- 15.0
+--
+-- >>> midpoint (singleton 5.0)
+-- 5.0
+--
+-- >>> midpoint empty
+-- NaN
+midpoint :: Fractional a => Interval a -> a
+midpoint (I a b) = a + (b - a) / 2
+midpoint Empty = nan
+{-# INLINE midpoint #-}
+
+-- | Determine if a point is in the interval.
+--
+-- >>> elem 3.2 (1.0 ... 5.0)
+-- True
+--
+-- >>> elem 5 (1.0 ... 5.0)
+-- True
+--
+-- >>> elem 1 (1.0 ... 5.0)
+-- True
+--
+-- >>> elem 8 (1.0 ... 5.0)
+-- False
+--
+-- >>> elem 5 empty
+-- False
+--
+elem :: Ord a => a -> Interval a -> Bool
+elem x (I a b) = x >= a && x <= b
+elem _ Empty = False
+{-# INLINE elem #-}
+
+-- | Determine if a point is not included in the interval
+--
+-- >>> notElem 8 (1.0 ... 5.0)
+-- True
+--
+-- >>> notElem 1.4 (1.0 ... 5.0)
+-- False
+--
+-- And of course, nothing is a member of the empty interval.
+--
+-- >>> notElem 5 empty
+-- True
+notElem :: Ord a => a -> Interval a -> Bool
+notElem x xs = not (elem x xs)
+{-# INLINE notElem #-}
+
+-- | 'realToFrac' will use the midpoint
+instance Real a => Real (Interval a) where
+  toRational Empty = nan
+  toRational (I ra rb) = a + (b - a) / 2 where
+    a = toRational ra
+    b = toRational rb
+  {-# INLINE toRational #-}
+
+instance Ord a => Ord (Interval a) where
+  compare Empty Empty = EQ
+  compare Empty _ = LT
+  compare _ Empty = GT
+  compare (I ax bx) (I ay by)
+    | bx < ay = LT
+    | ax > by = GT
+    | bx == ay && ax == by = EQ
+    | otherwise = error "Numeric.Interval.compare: ambiguous comparison"
+  {-# INLINE compare #-}
+
+  max (I a b) (I a' b') = max a a' ... max b b'
+  max Empty i = i
+  max i Empty = i
+  {-# INLINE max #-}
+
+  min (I a b) (I a' b') = min a a' ... min b b'
+  min Empty _ = Empty
+  min _ Empty = Empty
+  {-# INLINE min #-}
+
+-- @'divNonZero' X Y@ assumes @0 `'notElem'` Y@
+divNonZero :: (Fractional a, Ord a) => Interval a -> Interval a -> Interval a
+divNonZero (I a b) (I a' b') =
+  minimum [a / a', a / b', b / a', b / b']
+  ...
+  maximum [a / a', a / b', b / a', b / b']
+divNonZero _ _ = Empty
+
+-- @'divPositive' X y@ assumes y > 0, and divides @X@ by [0 ... y]
+divPositive :: (Fractional a, Ord a) => Interval a -> a -> Interval a
+divPositive Empty _ = Empty
+divPositive x@(I a b) y
+  | a == 0 && b == 0 = x
+  -- b < 0 || isNegativeZero b = negInfinity ... ( b / y)
+  | b < 0 = negInfinity ... (b / y)
+  | a < 0 = whole
+  | otherwise = (a / y) ... posInfinity
+{-# INLINE divPositive #-}
+
+-- divNegative assumes y < 0 and divides the interval @X@ by [y ... 0]
+divNegative :: (Fractional a, Ord a) => Interval a -> a -> Interval a
+divNegative Empty _ = Empty
+divNegative x@(I a b) y
+  | a == 0 && b == 0 = - x -- flip negative zeros
+  -- b < 0 || isNegativeZero b = (b / y) ... posInfinity
+  | b < 0 = (b / y) ... posInfinity
+  | a < 0 = whole
+  | otherwise = negInfinity ... (a / y)
+{-# INLINE divNegative #-}
+
+divZero :: (Fractional a, Ord a) => Interval a -> Interval a
+divZero x@(I a b)
+  | a == 0 && b == 0 = x
+  | otherwise        = whole
+divZero Empty = Empty
+{-# INLINE divZero #-}
+
+instance (Fractional a, Ord a) => Fractional (Interval a) where
+  -- TODO: check isNegativeZero properly
+  _ / Empty = Empty
+  x / y@(I a b)
+    | 0 `notElem` y = divNonZero x y
+    | iz && sz  = empty -- division by 0
+    | iz        = divPositive x a
+    |       sz  = divNegative x b
+    | otherwise = divZero x
+    where
+      iz = a == 0
+      sz = b == 0
+  recip Empty = Empty
+  recip (I a b)   = on min recip a b ... on max recip a b
+  {-# INLINE recip #-}
+  fromRational r  = let r' = fromRational r in I r' r'
+  {-# INLINE fromRational #-}
+
+instance RealFrac a => RealFrac (Interval a) where
+  properFraction x = (b, x - fromIntegral b)
+    where
+      b = truncate (midpoint x)
+  {-# INLINE properFraction #-}
+  ceiling x = ceiling (sup x)
+  {-# INLINE ceiling #-}
+  floor x = floor (inf x)
+  {-# INLINE floor #-}
+  round x = round (midpoint x)
+  {-# INLINE round #-}
+  truncate x = truncate (midpoint x)
+  {-# INLINE truncate #-}
+
+instance (RealFloat a, Ord a) => Floating (Interval a) where
+  pi = singleton pi
+  {-# INLINE pi #-}
+  exp = increasing exp
+  {-# INLINE exp #-}
+  log (I a b) = (if a > 0 then log a else negInfinity) ... log b
+  log Empty = Empty
+  {-# INLINE log #-}
+  cos Empty = Empty
+  cos x
+    | width t >= pi = (-1) ... 1
+    | inf t >= pi = - cos (t - pi)
+    | sup t <= pi = decreasing cos t
+    | sup t <= 2 * pi = (-1) ... cos ((pi * 2 - sup t) `min` inf t)
+    | otherwise = (-1) ... 1
+    where
+      t = fmod x (pi * 2)
+  {-# INLINE cos #-}
+  sin Empty = Empty
+  sin x = cos (x - pi / 2)
+  {-# INLINE sin #-}
+  tan Empty = Empty
+  tan x
+    | inf t' <= - pi / 2 || sup t' >= pi / 2 = whole
+    | otherwise = increasing tan x
+    where
+      t = x `fmod` pi
+      t' | t >= pi / 2 = t - pi
+         | otherwise    = t
+  {-# INLINE tan #-}
+  asin Empty = Empty
+  asin (I a b)
+    | b < -1 || a > 1 = Empty
+    | otherwise =
+      (if a <= -1 then -halfPi else asin a)
+      ...
+      (if b >= 1 then halfPi else asin b)
+    where
+      halfPi = pi / 2
+  {-# INLINE asin #-}
+  acos Empty = Empty
+  acos (I a b)
+    | b < -1 || a > 1 = Empty
+    | otherwise =
+      (if b >= 1 then 0 else acos b)
+      ...
+      (if a < -1 then pi else acos a)
+  {-# INLINE acos #-}
+  atan = increasing atan
+  {-# INLINE atan #-}
+  sinh = increasing sinh
+  {-# INLINE sinh #-}
+  cosh Empty = Empty
+  cosh x@(I a b)
+    | b < 0  = decreasing cosh x
+    | a >= 0 = increasing cosh x
+    | otherwise  = I 0 $ cosh $ if - a > b
+                                then a
+                                else b
+  {-# INLINE cosh #-}
+  tanh = increasing tanh
+  {-# INLINE tanh #-}
+  asinh = increasing asinh
+  {-# INLINE asinh #-}
+  acosh Empty = Empty
+  acosh (I a b)
+    | b < 1 = Empty
+    | otherwise = I lo $ acosh b
+    where lo | a <= 1 = 0
+             | otherwise = acosh a
+  {-# INLINE acosh #-}
+  atanh Empty = Empty
+  atanh (I a b)
+    | b < -1 || a > 1 = Empty
+    | otherwise =
+      (if a <= - 1 then negInfinity else atanh a)
+      ...
+      (if b >= 1 then posInfinity else atanh b)
+  {-# INLINE atanh #-}
+
+-- | lift a monotone increasing function over a given interval
+increasing :: (a -> b) -> Interval a -> Interval b
+increasing f (I a b) = I (f a) (f b)
+increasing _ Empty = Empty
+
+-- | lift a monotone decreasing function over a given interval
+decreasing :: (a -> b) -> Interval a -> Interval b
+decreasing f (I a b) = I (f b) (f a)
+decreasing _ Empty = Empty
+
+-- | We have to play some semantic games to make these methods make sense.
+-- Most compute with the midpoint of the interval.
+instance RealFloat a => RealFloat (Interval a) where
+  floatRadix = floatRadix . midpoint
+
+  floatDigits = floatDigits . midpoint
+  floatRange = floatRange . midpoint
+  decodeFloat = decodeFloat . midpoint
+  encodeFloat m e = singleton (encodeFloat m e)
+  exponent = exponent . midpoint
+  significand x = min a b ... max a b
+    where
+      (_ ,em) = decodeFloat (midpoint x)
+      (mi,ei) = decodeFloat (inf x)
+      (ms,es) = decodeFloat (sup x)
+      a = encodeFloat mi (ei - em - floatDigits x)
+      b = encodeFloat ms (es - em - floatDigits x)
+  scaleFloat _ Empty = Empty
+  scaleFloat n (I a b) = I (scaleFloat n a) (scaleFloat n b)
+  isNaN (I a b) = isNaN a || isNaN b
+  isNaN Empty = True
+  isInfinite (I a b) = isInfinite a || isInfinite b
+  isInfinite Empty = False
+  isDenormalized (I a b) = isDenormalized a || isDenormalized b
+  isDenormalized Empty = False
+  -- contains negative zero
+  isNegativeZero (I a b) = not (a > 0)
+                  && not (b < 0)
+                  && (  (b == 0 && (a < 0 || isNegativeZero a))
+                     || (a == 0 && isNegativeZero a)
+                     || (a < 0 && b >= 0))
+  isNegativeZero Empty = False
+  isIEEE _ = False
+
+  atan2 = error "unimplemented"
+
+-- TODO: (^), (^^) to give tighter bounds
+
+-- | Calculate the intersection of two intervals.
+--
+-- >>> intersection (1 ... 10 :: Interval Double) (5 ... 15 :: Interval Double)
+-- 5.0 ... 10.0
+intersection :: (Fractional a, Ord a) => Interval a -> Interval a -> Interval a
+intersection x@(I a b) y@(I a' b')
+  | x /=! y   = Empty
+  | otherwise = I (max a a') (min b b')
+intersection _ _ = Empty
+{-# INLINE intersection #-}
+
+-- | Calculate the convex hull of two intervals
+--
+-- >>> hull (0 ... 10 :: Interval Double) (5 ... 15 :: Interval Double)
+-- 0.0 ... 15.0
+--
+-- >>> hull (15 ... 85 :: Interval Double) (0 ... 10 :: Interval Double)
+-- 0.0 ... 85.0
+hull :: Ord a => Interval a -> Interval a -> Interval a
+hull (I a b) (I a' b') = I (min a a') (max b b')
+hull Empty x = x
+hull x Empty = x
+{-# INLINE hull #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '<' y@
+--
+-- >>> (5 ... 10 :: Interval Double) <! (20 ... 30 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 10 :: Interval Double) <! (10 ... 30 :: Interval Double)
+-- False
+--
+-- >>> (20 ... 30 :: Interval Double) <! (5 ... 10 :: Interval Double)
+-- False
+(<!)  :: Ord a => Interval a -> Interval a -> Bool
+Empty <! _ = True
+_ <! Empty = True
+I _ bx <! I ay _ = bx < ay
+{-# INLINE (<!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '<=' y@
+--
+-- >>> (5 ... 10 :: Interval Double) <=! (20 ... 30 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 10 :: Interval Double) <=! (10 ... 30 :: Interval Double)
+-- True
+--
+-- >>> (20 ... 30 :: Interval Double) <=! (5 ... 10 :: Interval Double)
+-- False
+(<=!) :: Ord a => Interval a -> Interval a -> Bool
+Empty <=! _ = True
+_ <=! Empty = True
+I _ bx <=! I ay _ = bx <= ay
+{-# INLINE (<=!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '==' y@
+--
+-- Only singleton intervals or empty intervals can return true
+--
+-- >>> (singleton 5 :: Interval Double) ==! (singleton 5 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 10 :: Interval Double) ==! (5 ... 10 :: Interval Double)
+-- False
+(==!) :: Eq a => Interval a -> Interval a -> Bool
+Empty ==! _ = True
+_ ==! Empty = True
+I ax bx ==! I ay by = bx == ay && ax == by
+{-# INLINE (==!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '/=' y@
+--
+-- >>> (5 ... 15 :: Interval Double) /=! (20 ... 40 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 15 :: Interval Double) /=! (15 ... 40 :: Interval Double)
+-- False
+(/=!) :: Ord a => Interval a -> Interval a -> Bool
+Empty /=! _ = True
+_ /=! Empty = True
+I ax bx /=! I ay by = bx < ay || ax > by
+{-# INLINE (/=!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '>' y@
+--
+-- >>> (20 ... 40 :: Interval Double) >! (10 ... 19 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 20 :: Interval Double) >! (15 ... 40 :: Interval Double)
+-- False
+(>!)  :: Ord a => Interval a -> Interval a -> Bool
+Empty >! _ = True
+_ >! Empty = True
+I ax _ >! I _ by = ax > by
+{-# INLINE (>!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '>=' y@
+--
+-- >>> (20 ... 40 :: Interval Double) >=! (10 ... 20 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 20 :: Interval Double) >=! (15 ... 40 :: Interval Double)
+-- False
+(>=!) :: Ord a => Interval a -> Interval a -> Bool
+Empty >=! _ = True
+_ >=! Empty = True
+I ax _ >=! I _ by = ax >= by
+{-# INLINE (>=!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x `op` y@
+certainly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool
+certainly cmp l r
+    | lt && eq && gt = True
+    | lt && eq       = l <=! r
+    | lt &&       gt = l /=! r
+    | lt             = l <!  r
+    |       eq && gt = l >=! r
+    |       eq       = l ==! r
+    |             gt = l >!  r
+    | otherwise      = False
+    where
+        lt = cmp False True
+        eq = cmp True True
+        gt = cmp True False
+{-# INLINE certainly #-}
+
+-- | Check if interval @X@ totally contains interval @Y@
+--
+-- >>> (20 ... 40 :: Interval Double) `contains` (25 ... 35 :: Interval Double)
+-- True
+--
+-- >>> (20 ... 40 :: Interval Double) `contains` (15 ... 35 :: Interval Double)
+-- False
+contains :: Ord a => Interval a -> Interval a -> Bool
+contains _ Empty = True
+contains (I ax bx) (I ay by) = ax <= ay && by <= bx
+contains Empty I{} = False
+{-# INLINE contains #-}
+
+-- | Flipped version of `contains`. Check if interval @X@ a subset of interval @Y@
+--
+-- >>> (25 ... 35 :: Interval Double) `isSubsetOf` (20 ... 40 :: Interval Double)
+-- True
+--
+-- >>> (20 ... 40 :: Interval Double) `isSubsetOf` (15 ... 35 :: Interval Double)
+-- False
+isSubsetOf :: Ord a => Interval a -> Interval a -> Bool
+isSubsetOf = flip contains
+{-# INLINE isSubsetOf #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<' y@?
+(<?) :: Ord a => Interval a -> Interval a -> Bool
+Empty <? _ = False
+_ <? Empty = False
+I ax _ <? I _ by = ax < by
+{-# INLINE (<?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<=' y@?
+(<=?) :: Ord a => Interval a -> Interval a -> Bool
+Empty <=? _ = False
+_ <=? Empty = False
+I ax _ <=? I _ by = ax <= by
+{-# INLINE (<=?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '==' y@?
+(==?) :: Ord a => Interval a -> Interval a -> Bool
+I ax bx ==? I ay by = ax <= by && bx >= ay
+_ ==? _ = False
+{-# INLINE (==?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '/=' y@?
+(/=?) :: Eq a => Interval a -> Interval a -> Bool
+I ax bx /=? I ay by = ax /= by || bx /= ay
+_ /=? _ = False
+{-# INLINE (/=?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>' y@?
+(>?) :: Ord a => Interval a -> Interval a -> Bool
+I _ bx >? I ay _ = bx > ay
+_ >? _ = False
+{-# INLINE (>?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>=' y@?
+(>=?) :: Ord a => Interval a -> Interval a -> Bool
+I _ bx >=? I ay _ = bx >= ay
+_ >=? _ = False
+{-# INLINE (>=?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x `op` y@?
+possibly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool
+possibly cmp l r
+    | lt && eq && gt = True
+    | lt && eq       = l <=? r
+    | lt &&       gt = l /=? r
+    | lt             = l <? r
+    |       eq && gt = l >=? r
+    |       eq       = l ==? r
+    |             gt = l >? r
+    | otherwise      = False
+    where
+        lt = cmp LT EQ
+        eq = cmp EQ EQ
+        gt = cmp GT EQ
+{-# INLINE possibly #-}
+
+-- | id function. Useful for type specification
+--
+-- >>> :t idouble (1 ... 3)
+-- idouble (1 ... 3) :: Interval Double
+idouble :: Interval Double -> Interval Double
+idouble = id
+
+-- | id function. Useful for type specification
+--
+-- >>> :t ifloat (1 ... 3)
+-- ifloat (1 ... 3) :: Interval Float
+ifloat :: Interval Float -> Interval Float
+ifloat = id
+
+-- Bugs:
+-- sin 1 :: Interval Double
+
+default (Integer,Double)
diff --git a/src/Numeric/Interval/Kaucher.hs b/src/Numeric/Interval/Kaucher.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Interval/Kaucher.hs
@@ -0,0 +1,746 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
+{-# LANGUAGE DeriveGeneric #-}
+#endif
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Interval
+-- Copyright   :  (c) Edward Kmett 2010-2014
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  DeriveDataTypeable
+--
+-- "Directed" Interval arithmetic
+--
+-----------------------------------------------------------------------------
+
+module Numeric.Interval.Kaucher
+  ( Interval(..)
+  , (...)
+  , whole
+  , empty
+  , null
+  , singleton
+  , elem
+  , notElem
+  , inf
+  , sup
+  , singular
+  , width
+  , midpoint
+  , intersection
+  , hull
+  , bisect
+  , magnitude
+  , mignitude
+  , contains
+  , isSubsetOf
+  , certainly, (<!), (<=!), (==!), (>=!), (>!)
+  , possibly, (<?), (<=?), (==?), (>=?), (>?)
+  , clamp
+  , idouble
+  , ifloat
+  ) where
+
+import Control.Applicative hiding (empty)
+import Data.Data
+import Data.Distributive
+import Data.Foldable hiding (minimum, maximum, elem, notElem)
+import Data.Function (on)
+import Data.Monoid
+import Data.Traversable
+#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
+import GHC.Generics
+#endif
+import Prelude hiding (null, elem, notElem)
+
+-- $setup
+
+data Interval a = I !a !a deriving
+  ( Data
+  , Typeable
+#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
+  , Generic
+#if __GLASGOW_HASKELL__ >= 706
+  , Generic1
+#endif
+#endif
+  )
+
+instance Functor Interval where
+  fmap f (I a b) = I (f a) (f b)
+  {-# INLINE fmap #-}
+
+instance Foldable Interval where
+  foldMap f (I a b) = f a `mappend` f b
+  {-# INLINE foldMap #-}
+
+instance Traversable Interval where
+  traverse f (I a b) = I <$> f a <*> f b
+  {-# INLINE traverse #-}
+
+instance Applicative Interval where
+  pure a = I a a
+  {-# INLINE pure #-}
+  I f g <*> I a b = I (f a) (g b)
+  {-# INLINE (<*>) #-}
+
+instance Monad Interval where
+  return a = I a a
+  {-# INLINE return #-}
+  I a b >>= f = I a' b' where
+    I a' _ = f a
+    I _ b' = f b
+  {-# INLINE (>>=) #-}
+
+instance Distributive Interval where
+  distribute f = fmap inf f ... fmap sup f
+  {-# INLINE distribute #-}
+
+infix 3 ...
+
+negInfinity :: Fractional a => a
+negInfinity = (-1)/0
+{-# INLINE negInfinity #-}
+
+posInfinity :: Fractional a => a
+posInfinity = 1/0
+{-# INLINE posInfinity #-}
+
+nan :: Fractional a => a
+nan = 0/0
+
+fmod :: RealFrac a => a -> a -> a
+fmod a b = a - q*b where
+  q = realToFrac (truncate $ a / b :: Integer)
+{-# INLINE fmod #-}
+
+(...) :: a -> a -> Interval a
+(...) = I
+{-# INLINE (...) #-}
+
+-- | The whole real number line
+--
+-- >>> whole
+-- -Infinity ... Infinity
+whole :: Fractional a => Interval a
+whole = negInfinity ... posInfinity
+{-# INLINE whole #-}
+
+-- | An empty interval
+--
+-- >>> empty
+-- NaN ... NaN
+empty :: Fractional a => Interval a
+empty = nan ... nan
+{-# INLINE empty #-}
+
+-- | negation handles NaN properly
+--
+-- >>> null (1 ... 5)
+-- False
+--
+-- >>> null (1 ... 1)
+-- False
+--
+-- >>> null empty
+-- True
+null :: Ord a => Interval a -> Bool
+null x = not (inf x <= sup x)
+{-# INLINE null #-}
+
+-- | A singleton point
+--
+-- >>> singleton 1
+-- 1 ... 1
+singleton :: a -> Interval a
+singleton a = a ... a
+{-# INLINE singleton #-}
+
+-- | The infinumum (lower bound) of an interval
+--
+-- >>> inf (1 ... 20)
+-- 1
+inf :: Interval a -> a
+inf (I a _) = a
+{-# INLINE inf #-}
+
+-- | The supremum (upper bound) of an interval
+--
+-- >>> sup (1 ... 20)
+-- 20
+sup :: Interval a -> a
+sup (I _ b) = b
+{-# INLINE sup #-}
+
+-- | Is the interval a singleton point?
+-- N.B. This is fairly fragile and likely will not hold after
+-- even a few operations that only involve singletons
+--
+-- >>> singular (singleton 1)
+-- True
+--
+-- >>> singular (1.0 ... 20.0)
+-- False
+singular :: Ord a => Interval a -> Bool
+singular x = not (null x) && inf x == sup x
+{-# INLINE singular #-}
+
+instance Eq a => Eq (Interval a) where
+  (==) = (==!)
+  {-# INLINE (==) #-}
+
+instance Show a => Show (Interval a) where
+  showsPrec n (I a b) =
+    showParen (n > 3) $
+      showsPrec 3 a .
+      showString " ... " .
+      showsPrec 3 b
+
+-- | Calculate the width of an interval.
+--
+-- >>> width (1 ... 20)
+-- 19
+--
+-- >>> width (singleton 1)
+-- 0
+--
+-- >>> width empty
+-- NaN
+width :: Num a => Interval a -> a
+width (I a b) = b - a
+{-# INLINE width #-}
+
+-- | Magnitude
+--
+-- >>> magnitude (1 ... 20)
+-- 20
+--
+-- >>> magnitude (-20 ... 10)
+-- 20
+--
+-- >>> magnitude (singleton 5)
+-- 5
+magnitude :: (Num a, Ord a) => Interval a -> a
+magnitude x = (max `on` abs) (inf x) (sup x)
+{-# INLINE magnitude #-}
+
+-- | \"mignitude\"
+--
+-- >>> mignitude (1 ... 20)
+-- 1
+--
+-- >>> mignitude (-20 ... 10)
+-- 10
+--
+-- >>> mignitude (singleton 5)
+-- 5
+mignitude :: (Num a, Ord a) => Interval a -> a
+mignitude x = (min `on` abs) (inf x) (sup x)
+{-# INLINE mignitude #-}
+
+instance (Num a, Ord a) => Num (Interval a) where
+  I a b + I a' b' = (a + a') ... (b + b')
+  {-# INLINE (+) #-}
+  I a b - I a' b' = (a - b') ... (b - a')
+  {-# INLINE (-) #-}
+  I a b * I a' b' =
+    minimum [a * a', a * b', b * a', b * b']
+    ...
+    maximum [a * a', a * b', b * a', b * b']
+  {-# INLINE (*) #-}
+  abs x@(I a b)
+    | a >= 0    = x
+    | b <= 0    = negate x
+    | otherwise = 0 ... max (- a) b
+  {-# INLINE abs #-}
+
+  signum = increasing signum
+  {-# INLINE signum #-}
+
+  fromInteger i = singleton (fromInteger i)
+  {-# INLINE fromInteger #-}
+
+-- | Bisect an interval at its midpoint.
+--
+-- >>> bisect (10.0 ... 20.0)
+-- (10.0 ... 15.0,15.0 ... 20.0)
+--
+-- >>> bisect (singleton 5.0)
+-- (5.0 ... 5.0,5.0 ... 5.0)
+--
+-- >>> bisect empty
+-- (NaN ... NaN,NaN ... NaN)
+bisect :: Fractional a => Interval a -> (Interval a, Interval a)
+bisect x = (inf x ... m, m ... sup x) where m = midpoint x
+{-# INLINE bisect #-}
+
+-- | Nearest point to the midpoint of the interval.
+--
+-- >>> midpoint (10.0 ... 20.0)
+-- 15.0
+--
+-- >>> midpoint (singleton 5.0)
+-- 5.0
+--
+-- >>> midpoint empty
+-- NaN
+midpoint :: Fractional a => Interval a -> a
+midpoint x = inf x + (sup x - inf x) / 2
+{-# INLINE midpoint #-}
+
+-- | Determine if a point is in the interval.
+--
+-- >>> elem 3.2 (1.0 ... 5.0)
+-- True
+--
+-- >>> elem 5 (1.0 ... 5.0)
+-- True
+--
+-- >>> elem 1 (1.0 ... 5.0)
+-- True
+--
+-- >>> elem 8 (1.0 ... 5.0)
+-- False
+--
+-- >>> elem 5 empty
+-- False
+--
+elem :: Ord a => a -> Interval a -> Bool
+elem x xs = x >= inf xs && x <= sup xs
+{-# INLINE elem #-}
+
+-- | Determine if a point is not included in the interval
+--
+-- >>> notElem 8 (1.0 ... 5.0)
+-- True
+--
+-- >>> notElem 1.4 (1.0 ... 5.0)
+-- False
+--
+-- And of course, nothing is a member of the empty interval.
+--
+-- >>> notElem 5 empty
+-- True
+notElem :: Ord a => a -> Interval a -> Bool
+notElem x xs = not (elem x xs)
+{-# INLINE notElem #-}
+
+-- | 'realToFrac' will use the midpoint
+instance Real a => Real (Interval a) where
+  toRational x
+    | null x   = nan
+    | otherwise = a + (b - a) / 2
+    where
+      a = toRational (inf x)
+      b = toRational (sup x)
+  {-# INLINE toRational #-}
+
+instance Ord a => Ord (Interval a) where
+  compare x y
+    | sup x < inf y = LT
+    | inf x > sup y = GT
+    | sup x == inf y && inf x == sup y = EQ
+    | otherwise = error "Numeric.Interval.compare: ambiguous comparison"
+  {-# INLINE compare #-}
+
+  max (I a b) (I a' b') = max a a' ... max b b'
+  {-# INLINE max #-}
+
+  min (I a b) (I a' b') = min a a' ... min b b'
+  {-# INLINE min #-}
+
+-- @'divNonZero' X Y@ assumes @0 `'notElem'` Y@
+divNonZero :: (Fractional a, Ord a) => Interval a -> Interval a -> Interval a
+divNonZero (I a b) (I a' b') =
+  minimum [a / a', a / b', b / a', b / b']
+  ...
+  maximum [a / a', a / b', b / a', b / b']
+
+-- @'divPositive' X y@ assumes y > 0, and divides @X@ by [0 ... y]
+divPositive :: (Fractional a, Ord a) => Interval a -> a -> Interval a
+divPositive x@(I a b) y
+  | a == 0 && b == 0 = x
+  -- b < 0 || isNegativeZero b = negInfinity ... ( b / y)
+  | b < 0 = negInfinity ... ( b / y)
+  | a < 0 = whole
+  | otherwise = (a / y) ... posInfinity
+{-# INLINE divPositive #-}
+
+-- divNegative assumes y < 0 and divides the interval @X@ by [y ... 0]
+divNegative :: (Fractional a, Ord a) => Interval a -> a -> Interval a
+divNegative x@(I a b) y
+  | a == 0 && b == 0 = - x -- flip negative zeros
+  -- b < 0 || isNegativeZero b = (b / y) ... posInfinity
+  | b < 0 = (b / y) ... posInfinity
+  | a < 0 = whole
+  | otherwise = negInfinity ... (a / y)
+{-# INLINE divNegative #-}
+
+divZero :: (Fractional a, Ord a) => Interval a -> Interval a
+divZero x
+  | inf x == 0 && sup x == 0 = x
+  | otherwise = whole
+{-# INLINE divZero #-}
+
+instance (Fractional a, Ord a) => Fractional (Interval a) where
+  -- TODO: check isNegativeZero properly
+  x / y
+    | 0 `notElem` y = divNonZero x y
+    | iz && sz  = empty -- division by 0
+    | iz        = divPositive x (inf y)
+    |       sz  = divNegative x (sup y)
+    | otherwise = divZero x
+    where
+      iz = inf y == 0
+      sz = sup y == 0
+  recip (I a b)   = on min recip a b ... on max recip a b
+  {-# INLINE recip #-}
+  fromRational r  = let r' = fromRational r in r' ... r'
+  {-# INLINE fromRational #-}
+
+instance RealFrac a => RealFrac (Interval a) where
+  properFraction x = (b, x - fromIntegral b)
+    where
+      b = truncate (midpoint x)
+  {-# INLINE properFraction #-}
+  ceiling x = ceiling (sup x)
+  {-# INLINE ceiling #-}
+  floor x = floor (inf x)
+  {-# INLINE floor #-}
+  round x = round (midpoint x)
+  {-# INLINE round #-}
+  truncate x = truncate (midpoint x)
+  {-# INLINE truncate #-}
+
+instance (RealFloat a, Ord a) => Floating (Interval a) where
+  pi = singleton pi
+  {-# INLINE pi #-}
+  exp = increasing exp
+  {-# INLINE exp #-}
+  log (I a b) = (if a > 0 then log a else negInfinity) ... log b
+  {-# INLINE log #-}
+  cos x
+    | null x = empty
+    | width t >= pi = (-1) ... 1
+    | inf t >= pi = - cos (t - pi)
+    | sup t <= pi = decreasing cos t
+    | sup t <= 2 * pi = (-1) ... cos ((pi * 2 - sup t) `min` inf t)
+    | otherwise = (-1) ... 1
+    where
+      t = fmod x (pi * 2)
+  {-# INLINE cos #-}
+  sin x
+    | null x = empty
+    | otherwise = cos (x - pi / 2)
+  {-# INLINE sin #-}
+  tan x
+    | null x = empty
+    | inf t' <= - pi / 2 || sup t' >= pi / 2 = whole
+    | otherwise = increasing tan x
+    where
+      t = x `fmod` pi
+      t' | t >= pi / 2 = t - pi
+         | otherwise    = t
+  {-# INLINE tan #-}
+  asin x@(I a b)
+    | null x || b < -1 || a > 1 = empty
+    | otherwise =
+      (if a <= -1 then -halfPi else asin a)
+      ...
+      (if b >= 1 then halfPi else asin b)
+    where
+      halfPi = pi / 2
+  {-# INLINE asin #-}
+  acos x@(I a b)
+    | null x || b < -1 || a > 1 = empty
+    | otherwise =
+      (if b >= 1 then 0 else acos b)
+      ...
+      (if a < -1 then pi else acos a)
+  {-# INLINE acos #-}
+  atan = increasing atan
+  {-# INLINE atan #-}
+  sinh = increasing sinh
+  {-# INLINE sinh #-}
+  cosh x@(I a b)
+    | null x = empty
+    | b < 0  = decreasing cosh x
+    | a >= 0 = increasing cosh x
+    | otherwise  = I 0 $ cosh $ if - a > b
+                                then a
+                                else b
+  {-# INLINE cosh #-}
+  tanh = increasing tanh
+  {-# INLINE tanh #-}
+  asinh = increasing asinh
+  {-# INLINE asinh #-}
+  acosh x@(I a b)
+    | null x || b < 1 = empty
+    | otherwise = I lo $ acosh b
+    where lo | a <= 1 = 0
+             | otherwise = acosh a
+  {-# INLINE acosh #-}
+  atanh x@(I a b)
+    | null x || b < -1 || a > 1 = empty
+    | otherwise =
+      (if a <= - 1 then negInfinity else atanh a)
+      ...
+      (if b >= 1 then posInfinity else atanh b)
+  {-# INLINE atanh #-}
+
+-- | lift a monotone increasing function over a given interval
+increasing :: (a -> b) -> Interval a -> Interval b
+increasing f (I a b) = f a ... f b
+
+-- | lift a monotone decreasing function over a given interval
+decreasing :: (a -> b) -> Interval a -> Interval b
+decreasing f (I a b) = f b ... f a
+
+-- | We have to play some semantic games to make these methods make sense.
+-- Most compute with the midpoint of the interval.
+instance RealFloat a => RealFloat (Interval a) where
+  floatRadix = floatRadix . midpoint
+
+  floatDigits = floatDigits . midpoint
+  floatRange = floatRange . midpoint
+  decodeFloat = decodeFloat . midpoint
+  encodeFloat m e = singleton (encodeFloat m e)
+  exponent = exponent . midpoint
+  significand x = min a b ... max a b
+    where
+      (_ ,em) = decodeFloat (midpoint x)
+      (mi,ei) = decodeFloat (inf x)
+      (ms,es) = decodeFloat (sup x)
+      a = encodeFloat mi (ei - em - floatDigits x)
+      b = encodeFloat ms (es - em - floatDigits x)
+  scaleFloat n x = scaleFloat n (inf x) ... scaleFloat n (sup x)
+  isNaN x = isNaN (inf x) || isNaN (sup x)
+  isInfinite x = isInfinite (inf x) || isInfinite (sup x)
+  isDenormalized x = isDenormalized (inf x) || isDenormalized (sup x)
+  -- contains negative zero
+  isNegativeZero x = not (inf x > 0)
+                  && not (sup x < 0)
+                  && (  (sup x == 0 && (inf x < 0 || isNegativeZero (inf x)))
+                     || (inf x == 0 && isNegativeZero (inf x))
+                     || (inf x < 0 && sup x >= 0))
+  isIEEE x = isIEEE (inf x) && isIEEE (sup x)
+  atan2 = error "unimplemented"
+
+-- TODO: (^), (^^) to give tighter bounds
+
+-- | Calculate the intersection of two intervals.
+--
+-- >>> intersection (1 ... 10 :: Interval Double) (5 ... 15 :: Interval Double)
+-- 5.0 ... 10.0
+intersection :: (Fractional a, Ord a) => Interval a -> Interval a -> Interval a
+intersection x@(I a b) y@(I a' b')
+  | x /=! y = empty
+  | otherwise = max a a' ... min b b'
+{-# INLINE intersection #-}
+
+-- | Calculate the convex hull of two intervals
+--
+-- >>> hull (0 ... 10 :: Interval Double) (5 ... 15 :: Interval Double)
+-- 0.0 ... 15.0
+--
+-- >>> hull (15 ... 85 :: Interval Double) (0 ... 10 :: Interval Double)
+-- 0.0 ... 85.0
+hull :: Ord a => Interval a -> Interval a -> Interval a
+hull x@(I a b) y@(I a' b')
+  | null x = y
+  | null y = x
+  | otherwise = min a a' ... max b b'
+{-# INLINE hull #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '<' y@
+--
+-- >>> (5 ... 10 :: Interval Double) <! (20 ... 30 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 10 :: Interval Double) <! (10 ... 30 :: Interval Double)
+-- False
+--
+-- >>> (20 ... 30 :: Interval Double) <! (5 ... 10 :: Interval Double)
+-- False
+(<!)  :: Ord a => Interval a -> Interval a -> Bool
+x <! y = sup x < inf y
+{-# INLINE (<!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '<=' y@
+--
+-- >>> (5 ... 10 :: Interval Double) <=! (20 ... 30 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 10 :: Interval Double) <=! (10 ... 30 :: Interval Double)
+-- True
+--
+-- >>> (20 ... 30 :: Interval Double) <=! (5 ... 10 :: Interval Double)
+-- False
+(<=!) :: Ord a => Interval a -> Interval a -> Bool
+x <=! y = sup x <= inf y
+{-# INLINE (<=!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '==' y@
+--
+-- Only singleton intervals return true
+--
+-- >>> (singleton 5 :: Interval Double) ==! (singleton 5 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 10 :: Interval Double) ==! (5 ... 10 :: Interval Double)
+-- False
+(==!) :: Eq a => Interval a -> Interval a -> Bool
+x ==! y = sup x == inf y && inf x == sup y
+{-# INLINE (==!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '/=' y@
+--
+-- >>> (5 ... 15 :: Interval Double) /=! (20 ... 40 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 15 :: Interval Double) /=! (15 ... 40 :: Interval Double)
+-- False
+(/=!) :: Ord a => Interval a -> Interval a -> Bool
+x /=! y = sup x < inf y || inf x > sup y
+{-# INLINE (/=!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '>' y@
+--
+-- >>> (20 ... 40 :: Interval Double) >! (10 ... 19 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 20 :: Interval Double) >! (15 ... 40 :: Interval Double)
+-- False
+(>!)  :: Ord a => Interval a -> Interval a -> Bool
+x >! y = inf x > sup y
+{-# INLINE (>!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '>=' y@
+--
+-- >>> (20 ... 40 :: Interval Double) >=! (10 ... 20 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 20 :: Interval Double) >=! (15 ... 40 :: Interval Double)
+-- False
+(>=!) :: Ord a => Interval a -> Interval a -> Bool
+x >=! y = inf x >= sup y
+{-# INLINE (>=!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x `op` y@
+--
+--
+certainly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool
+certainly cmp l r
+    | lt && eq && gt = True
+    | lt && eq       = l <=! r
+    | lt &&       gt = l /=! r
+    | lt             = l <!  r
+    |       eq && gt = l >=! r
+    |       eq       = l ==! r
+    |             gt = l >!  r
+    | otherwise      = False
+    where
+        lt = cmp LT EQ
+        eq = cmp EQ EQ
+        gt = cmp GT EQ
+{-# INLINE certainly #-}
+
+-- | Check if interval @X@ totally contains interval @Y@
+--
+-- >>> (20 ... 40 :: Interval Double) `contains` (25 ... 35 :: Interval Double)
+-- True
+--
+-- >>> (20 ... 40 :: Interval Double) `contains` (15 ... 35 :: Interval Double)
+-- False
+contains :: Ord a => Interval a -> Interval a -> Bool
+contains x y = null y
+            || (not (null x) && inf x <= inf y && sup y <= sup x)
+{-# INLINE contains #-}
+
+-- | Flipped version of `contains`. Check if interval @X@ a subset of interval @Y@
+--
+-- >>> (25 ... 35 :: Interval Double) `isSubsetOf` (20 ... 40 :: Interval Double)
+-- True
+--
+-- >>> (20 ... 40 :: Interval Double) `isSubsetOf` (15 ... 35 :: Interval Double)
+-- False
+isSubsetOf :: Ord a => Interval a -> Interval a -> Bool
+isSubsetOf = flip contains
+{-# INLINE isSubsetOf #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<' y@?
+(<?) :: Ord a => Interval a -> Interval a -> Bool
+x <? y = inf x < sup y
+{-# INLINE (<?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<=' y@?
+(<=?) :: Ord a => Interval a -> Interval a -> Bool
+x <=? y = inf x <= sup y
+{-# INLINE (<=?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '==' y@?
+(==?) :: Ord a => Interval a -> Interval a -> Bool
+x ==? y = inf x <= sup y && sup x >= inf y
+{-# INLINE (==?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '/=' y@?
+(/=?) :: Eq a => Interval a -> Interval a -> Bool
+x /=? y = inf x /= sup y || sup x /= inf y
+{-# INLINE (/=?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>' y@?
+(>?) :: Ord a => Interval a -> Interval a -> Bool
+x >? y = sup x > inf y
+{-# INLINE (>?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>=' y@?
+(>=?) :: Ord a => Interval a -> Interval a -> Bool
+x >=? y = sup x >= inf y
+{-# INLINE (>=?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x `op` y@?
+possibly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool
+possibly cmp l r
+    | lt && eq && gt = True
+    | lt && eq       = l <=? r
+    | lt &&       gt = l /=? r
+    | lt             = l <? r
+    |       eq && gt = l >=? r
+    |       eq       = l ==? r
+    |             gt = l >? r
+    | otherwise      = False
+    where
+        lt = cmp LT EQ
+        eq = cmp EQ EQ
+        gt = cmp GT EQ
+{-# INLINE possibly #-}
+
+-- | The nearest value to that supplied which is contained in the interval.
+clamp :: Ord a => Interval a -> a -> a
+clamp (I a b) x | x < a     = a
+                | x > b     = b
+                | otherwise = x
+
+-- | id function. Useful for type specification
+--
+-- >>> :t idouble (1 ... 3)
+-- idouble (1 ... 3) :: Interval Double
+idouble :: Interval Double -> Interval Double
+idouble = id
+
+-- | id function. Useful for type specification
+--
+-- >>> :t ifloat (1 ... 3)
+-- ifloat (1 ... 3) :: Interval Float
+ifloat :: Interval Float -> Interval Float
+ifloat = id
+
+-- Bugs:
+-- sin 1 :: Interval Double
+
+
+default (Integer,Double)
diff --git a/src/Numeric/Interval/NonEmpty.hs b/src/Numeric/Interval/NonEmpty.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Interval/NonEmpty.hs
@@ -0,0 +1,48 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
+{-# LANGUAGE DeriveGeneric #-}
+#endif
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Interval.NonEmpty
+-- Copyright   :  (c) Edward Kmett 2010-2013
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  DeriveDataTypeable
+--
+-- Interval arithmetic
+--
+-----------------------------------------------------------------------------
+
+module Numeric.Interval.NonEmpty
+  ( Interval
+  , (...)
+  , whole
+  , singleton
+  , elem
+  , notElem
+  , inf
+  , sup
+  , singular
+  , width
+  , midpoint
+  , intersection
+  , hull
+  , bisect
+  , bisectIntegral
+  , magnitude
+  , mignitude
+  , contains
+  , isSubsetOf
+  , certainly, (<!), (<=!), (==!), (>=!), (>!)
+  , possibly, (<?), (<=?), (==?), (>=?), (>?)
+  , clamp
+  , idouble
+  , ifloat
+  ) where
+
+import Numeric.Interval.NonEmpty.Internal
+import Prelude ()
diff --git a/src/Numeric/Interval/NonEmpty/Internal.hs b/src/Numeric/Interval/NonEmpty/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Interval/NonEmpty/Internal.hs
@@ -0,0 +1,650 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
+{-# LANGUAGE DeriveGeneric #-}
+#endif
+{-# OPTIONS_HADDOCK not-home #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Interval.NonEmpty.Internal
+-- Copyright   :  (c) Edward Kmett 2010-2014
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  DeriveDataTypeable
+--
+-- Interval arithmetic
+-----------------------------------------------------------------------------
+module Numeric.Interval.NonEmpty.Internal
+  ( Interval(..)
+  , (...)
+  , whole
+  , singleton
+  , elem
+  , notElem
+  , inf
+  , sup
+  , singular
+  , width
+  , midpoint
+  , intersection
+  , hull
+  , bisect
+  , bisectIntegral
+  , magnitude
+  , mignitude
+  , contains
+  , isSubsetOf
+  , certainly, (<!), (<=!), (==!), (>=!), (>!)
+  , possibly, (<?), (<=?), (==?), (>=?), (>?)
+  , clamp
+  , idouble
+  , ifloat
+  ) where
+
+import Data.Data
+import Data.Foldable hiding (minimum, maximum, elem, notElem)
+import Data.Function (on)
+import Data.Monoid
+#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
+import GHC.Generics
+#endif
+import Prelude hiding (null, elem, notElem)
+
+-- $setup
+
+data Interval a = I !a !a deriving
+  ( Data
+  , Typeable
+#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
+  , Generic
+#if __GLASGOW_HASKELL__ >= 706
+  , Generic1
+#endif
+#endif
+  )
+
+instance Foldable Interval where
+  foldMap f (I a b) = f a `mappend` f b
+  {-# INLINE foldMap #-}
+
+infix 3 ...
+
+negInfinity :: Fractional a => a
+negInfinity = (-1)/0
+{-# INLINE negInfinity #-}
+
+posInfinity :: Fractional a => a
+posInfinity = 1/0
+{-# INLINE posInfinity #-}
+
+fmod :: RealFrac a => a -> a -> a
+fmod a b = a - q*b where
+  q = realToFrac (truncate $ a / b :: Integer)
+{-# INLINE fmod #-}
+
+(...) :: Ord a => a -> a -> Interval a
+a ... b
+  | a <= b = I a b
+  | otherwise = I b a
+{-# INLINE (...) #-}
+
+-- | The whole real number line
+--
+-- >>> whole
+-- -Infinity ... Infinity
+whole :: Fractional a => Interval a
+whole = I negInfinity posInfinity
+{-# INLINE whole #-}
+
+-- | A singleton point
+--
+-- >>> singleton 1
+-- 1 ... 1
+singleton :: a -> Interval a
+singleton a = I a a
+{-# INLINE singleton #-}
+
+-- | The infinumum (lower bound) of an interval
+--
+-- >>> inf (1 ... 20)
+-- 1
+inf :: Interval a -> a
+inf (I a _) = a
+{-# INLINE inf #-}
+
+-- | The supremum (upper bound) of an interval
+--
+-- >>> sup (1 ... 20)
+-- 20
+sup :: Interval a -> a
+sup (I _ b) = b
+{-# INLINE sup #-}
+
+-- | Is the interval a singleton point?
+-- N.B. This is fairly fragile and likely will not hold after
+-- even a few operations that only involve singletons
+--
+-- >>> singular (singleton 1)
+-- True
+--
+-- >>> singular (1.0 ... 20.0)
+-- False
+singular :: Ord a => Interval a -> Bool
+singular (I a b) = a == b
+{-# INLINE singular #-}
+
+instance Eq a => Eq (Interval a) where
+  (==) = (==!)
+  {-# INLINE (==) #-}
+
+instance Show a => Show (Interval a) where
+  showsPrec n (I a b) =
+    showParen (n > 3) $
+      showsPrec 3 a .
+      showString " ... " .
+      showsPrec 3 b
+
+-- | Calculate the width of an interval.
+--
+-- >>> width (1 ... 20)
+-- 19
+--
+-- >>> width (singleton 1)
+-- 0
+width :: Num a => Interval a -> a
+width (I a b) = b - a
+{-# INLINE width #-}
+
+-- | Magnitude
+--
+-- >>> magnitude (1 ... 20)
+-- 20
+--
+-- >>> magnitude (-20 ... 10)
+-- 20
+--
+-- >>> magnitude (singleton 5)
+-- 5
+magnitude :: (Num a, Ord a) => Interval a -> a
+magnitude (I a b) = on max abs a b
+{-# INLINE magnitude #-}
+
+-- | \"mignitude\"
+--
+-- >>> mignitude (1 ... 20)
+-- 1
+--
+-- >>> mignitude (-20 ... 10)
+-- 10
+--
+-- >>> mignitude (singleton 5)
+-- 5
+mignitude :: (Num a, Ord a) => Interval a -> a
+mignitude (I a b) = on min abs a b
+{-# INLINE mignitude #-}
+
+instance (Num a, Ord a) => Num (Interval a) where
+  I a b + I a' b' = (a + a') ... (b + b')
+  {-# INLINE (+) #-}
+  I a b - I a' b' = (a - b') ... (b - a')
+  {-# INLINE (-) #-}
+  I a b * I a' b' =
+    minimum [a * a', a * b', b * a', b * b']
+    ...
+    maximum [a * a', a * b', b * a', b * b']
+  {-# INLINE (*) #-}
+  abs x@(I a b)
+    | a >= 0    = x
+    | b <= 0    = negate x
+    | otherwise = 0 ... max (- a) b
+  {-# INLINE abs #-}
+
+  signum = increasing signum
+  {-# INLINE signum #-}
+
+  fromInteger i = singleton (fromInteger i)
+  {-# INLINE fromInteger #-}
+
+-- | Bisect an interval at its midpoint.
+--
+-- >>> bisect (10.0 ... 20.0)
+-- (10.0 ... 15.0,15.0 ... 20.0)
+--
+-- >>> bisect (singleton 5.0)
+-- (5.0 ... 5.0,5.0 ... 5.0)
+bisect :: Fractional a => Interval a -> (Interval a, Interval a)
+bisect (I a b) = (I a m, I m b) where m = a + (b - a) / 2
+{-# INLINE bisect #-}
+
+bisectIntegral :: Integral a => Interval a -> (Interval a, Interval a)
+bisectIntegral (I a b)
+  | a == m || b == m = (I a a, I b b)
+  | otherwise        = (I a m, I m b)
+  where m = a + (b - a) `div` 2
+{-# INLINE bisectIntegral #-}
+
+-- | Nearest point to the midpoint of the interval.
+--
+-- >>> midpoint (10.0 ... 20.0)
+-- 15.0
+--
+-- >>> midpoint (singleton 5.0)
+-- 5.0
+midpoint :: Fractional a => Interval a -> a
+midpoint (I a b) = a + (b - a) / 2
+{-# INLINE midpoint #-}
+
+-- | Determine if a point is in the interval.
+--
+-- >>> elem 3.2 (1.0 ... 5.0)
+-- True
+--
+-- >>> elem 5 (1.0 ... 5.0)
+-- True
+--
+-- >>> elem 1 (1.0 ... 5.0)
+-- True
+--
+-- >>> elem 8 (1.0 ... 5.0)
+-- False
+elem :: Ord a => a -> Interval a -> Bool
+elem x (I a b) = x >= a && x <= b
+{-# INLINE elem #-}
+
+-- | Determine if a point is not included in the interval
+--
+-- >>> notElem 8 (1.0 ... 5.0)
+-- True
+--
+-- >>> notElem 1.4 (1.0 ... 5.0)
+-- False
+notElem :: Ord a => a -> Interval a -> Bool
+notElem x xs = not (elem x xs)
+{-# INLINE notElem #-}
+
+-- | 'realToFrac' will use the midpoint
+instance Real a => Real (Interval a) where
+  toRational (I ra rb) = a + (b - a) / 2 where
+    a = toRational ra
+    b = toRational rb
+  {-# INLINE toRational #-}
+
+instance Ord a => Ord (Interval a) where
+  compare (I ax bx) (I ay by)
+    | bx < ay = LT
+    | ax > by = GT
+    | bx == ay && ax == by = EQ
+    | otherwise = error "Numeric.Interval.compare: ambiguous comparison"
+  {-# INLINE compare #-}
+
+  max (I a b) (I a' b') = max a a' ... max b b'
+  {-# INLINE max #-}
+
+  min (I a b) (I a' b') = min a a' ... min b b'
+  {-# INLINE min #-}
+
+-- @'divNonZero' X Y@ assumes @0 `'notElem'` Y@
+divNonZero :: (Fractional a, Ord a) => Interval a -> Interval a -> Interval a
+divNonZero (I a b) (I a' b') =
+  minimum [a / a', a / b', b / a', b / b']
+  ...
+  maximum [a / a', a / b', b / a', b / b']
+
+-- @'divPositive' X y@ assumes y > 0, and divides @X@ by [0 ... y]
+divPositive :: (Fractional a, Ord a) => Interval a -> a -> Interval a
+divPositive x@(I a b) y
+  | a == 0 && b == 0 = x
+  -- b < 0 || isNegativeZero b = negInfinity ... ( b / y)
+  | b < 0 = negInfinity ... (b / y)
+  | a < 0 = whole
+  | otherwise = (a / y) ... posInfinity
+{-# INLINE divPositive #-}
+
+-- divNegative assumes y < 0 and divides the interval @X@ by [y ... 0]
+divNegative :: (Fractional a, Ord a) => Interval a -> a -> Interval a
+divNegative x@(I a b) y
+  | a == 0 && b == 0 = - x -- flip negative zeros
+  -- b < 0 || isNegativeZero b = (b / y) ... posInfinity
+  | b < 0 = (b / y) ... posInfinity
+  | a < 0 = whole
+  | otherwise = negInfinity ... (a / y)
+{-# INLINE divNegative #-}
+
+divZero :: (Fractional a, Ord a) => Interval a -> Interval a
+divZero x@(I a b)
+  | a == 0 && b == 0 = x
+  | otherwise        = whole
+{-# INLINE divZero #-}
+
+instance (Fractional a, Ord a) => Fractional (Interval a) where
+  -- TODO: check isNegativeZero properly
+  x / y@(I a b)
+    | 0 `notElem` y = divNonZero x y
+    | iz && sz  = error "division by zero"
+    | iz        = divPositive x a
+    |       sz  = divNegative x b
+    | otherwise = divZero x
+    where
+      iz = a == 0
+      sz = b == 0
+  recip (I a b)   = on min recip a b ... on max recip a b
+  {-# INLINE recip #-}
+  fromRational r  = let r' = fromRational r in I r' r'
+  {-# INLINE fromRational #-}
+
+instance RealFrac a => RealFrac (Interval a) where
+  properFraction x = (b, x - fromIntegral b)
+    where
+      b = truncate (midpoint x)
+  {-# INLINE properFraction #-}
+  ceiling x = ceiling (sup x)
+  {-# INLINE ceiling #-}
+  floor x = floor (inf x)
+  {-# INLINE floor #-}
+  round x = round (midpoint x)
+  {-# INLINE round #-}
+  truncate x = truncate (midpoint x)
+  {-# INLINE truncate #-}
+
+instance (RealFloat a, Ord a) => Floating (Interval a) where
+  pi = singleton pi
+  {-# INLINE pi #-}
+  exp = increasing exp
+  {-# INLINE exp #-}
+  log (I a b) = (if a > 0 then log a else negInfinity) ... log b
+  {-# INLINE log #-}
+  cos x
+    | width t >= pi = (-1) ... 1
+    | inf t >= pi = - cos (t - pi)
+    | sup t <= pi = decreasing cos t
+    | sup t <= 2 * pi = (-1) ... cos ((pi * 2 - sup t) `min` inf t)
+    | otherwise = (-1) ... 1
+    where
+      t = fmod x (pi * 2)
+  {-# INLINE cos #-}
+  sin x = cos (x - pi / 2)
+  {-# INLINE sin #-}
+  tan x
+    | inf t' <= - pi / 2 || sup t' >= pi / 2 = whole
+    | otherwise = increasing tan x
+    where
+      t = x `fmod` pi
+      t' | t >= pi / 2 = t - pi
+         | otherwise    = t
+  {-# INLINE tan #-}
+  asin (I a b) = I (if a <= -1 then -halfPi else asin a) (if b >= 1 then halfPi else asin b)
+    where halfPi = pi / 2
+  {-# INLINE asin #-}
+  acos (I a b) = I (if b >= 1 then 0 else acos b) (if a < -1 then pi else acos a)
+  {-# INLINE acos #-}
+  atan = increasing atan
+  {-# INLINE atan #-}
+  sinh = increasing sinh
+  {-# INLINE sinh #-}
+  cosh x@(I a b)
+    | b < 0  = decreasing cosh x
+    | a >= 0 = increasing cosh x
+    | otherwise  = I 0 $ cosh $ if - a > b
+                                then a
+                                else b
+  {-# INLINE cosh #-}
+  tanh = increasing tanh
+  {-# INLINE tanh #-}
+  asinh = increasing asinh
+  {-# INLINE asinh #-}
+  acosh (I a b) = I lo $ acosh b
+    where lo | a <= 1 = 0
+             | otherwise = acosh a
+  {-# INLINE acosh #-}
+  atanh (I a b) = I (if a <= - 1 then negInfinity else atanh a) (if b >= 1 then posInfinity else atanh b)
+  {-# INLINE atanh #-}
+
+-- | lift a monotone increasing function over a given interval
+increasing :: (a -> b) -> Interval a -> Interval b
+increasing f (I a b) = I (f a) (f b)
+
+-- | lift a monotone decreasing function over a given interval
+decreasing :: (a -> b) -> Interval a -> Interval b
+decreasing f (I a b) = I (f b) (f a)
+
+-- | We have to play some semantic games to make these methods make sense.
+-- Most compute with the midpoint of the interval.
+instance RealFloat a => RealFloat (Interval a) where
+  floatRadix = floatRadix . midpoint
+
+  floatDigits = floatDigits . midpoint
+  floatRange = floatRange . midpoint
+  decodeFloat = decodeFloat . midpoint
+  encodeFloat m e = singleton (encodeFloat m e)
+  exponent = exponent . midpoint
+  significand x = min a b ... max a b
+    where
+      (_ ,em) = decodeFloat (midpoint x)
+      (mi,ei) = decodeFloat (inf x)
+      (ms,es) = decodeFloat (sup x)
+      a = encodeFloat mi (ei - em - floatDigits x)
+      b = encodeFloat ms (es - em - floatDigits x)
+  scaleFloat n (I a b) = I (scaleFloat n a) (scaleFloat n b)
+  isNaN (I a b) = isNaN a || isNaN b
+  isInfinite (I a b) = isInfinite a || isInfinite b
+  isDenormalized (I a b) = isDenormalized a || isDenormalized b
+  -- contains negative zero
+  isNegativeZero (I a b) = not (a > 0)
+                  && not (b < 0)
+                  && (  (b == 0 && (a < 0 || isNegativeZero a))
+                     || (a == 0 && isNegativeZero a)
+                     || (a < 0 && b >= 0))
+  isIEEE _ = False
+
+  atan2 = error "unimplemented"
+
+-- TODO: (^), (^^) to give tighter bounds
+
+-- | Calculate the intersection of two intervals.
+--
+-- >>> intersection (1 ... 10 :: Interval Double) (5 ... 15 :: Interval Double)
+-- Just (5.0 ... 10.0)
+intersection :: (Fractional a, Ord a) => Interval a -> Interval a -> Maybe (Interval a)
+intersection x@(I a b) y@(I a' b')
+  | x /=! y   = Nothing
+  | otherwise = Just $ I (max a a') (min b b')
+{-# INLINE intersection #-}
+
+-- | Calculate the convex hull of two intervals
+--
+-- >>> hull (0 ... 10 :: Interval Double) (5 ... 15 :: Interval Double)
+-- 0.0 ... 15.0
+--
+-- >>> hull (15 ... 85 :: Interval Double) (0 ... 10 :: Interval Double)
+-- 0.0 ... 85.0
+hull :: Ord a => Interval a -> Interval a -> Interval a
+hull (I a b) (I a' b') = I (min a a') (max b b')
+{-# INLINE hull #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '<' y@
+--
+-- >>> (5 ... 10 :: Interval Double) <! (20 ... 30 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 10 :: Interval Double) <! (10 ... 30 :: Interval Double)
+-- False
+--
+-- >>> (20 ... 30 :: Interval Double) <! (5 ... 10 :: Interval Double)
+-- False
+(<!)  :: Ord a => Interval a -> Interval a -> Bool
+I _ bx <! I ay _ = bx < ay
+{-# INLINE (<!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '<=' y@
+--
+-- >>> (5 ... 10 :: Interval Double) <=! (20 ... 30 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 10 :: Interval Double) <=! (10 ... 30 :: Interval Double)
+-- True
+--
+-- >>> (20 ... 30 :: Interval Double) <=! (5 ... 10 :: Interval Double)
+-- False
+(<=!) :: Ord a => Interval a -> Interval a -> Bool
+I _ bx <=! I ay _ = bx <= ay
+{-# INLINE (<=!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '==' y@
+--
+-- Only singleton intervals or empty intervals can return true
+--
+-- >>> (singleton 5 :: Interval Double) ==! (singleton 5 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 10 :: Interval Double) ==! (5 ... 10 :: Interval Double)
+-- False
+(==!) :: Eq a => Interval a -> Interval a -> Bool
+I ax bx ==! I ay by = bx == ay && ax == by
+{-# INLINE (==!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '/=' y@
+--
+-- >>> (5 ... 15 :: Interval Double) /=! (20 ... 40 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 15 :: Interval Double) /=! (15 ... 40 :: Interval Double)
+-- False
+(/=!) :: Ord a => Interval a -> Interval a -> Bool
+I ax bx /=! I ay by = bx < ay || ax > by
+{-# INLINE (/=!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '>' y@
+--
+-- >>> (20 ... 40 :: Interval Double) >! (10 ... 19 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 20 :: Interval Double) >! (15 ... 40 :: Interval Double)
+-- False
+(>!)  :: Ord a => Interval a -> Interval a -> Bool
+I ax _ >! I _ by = ax > by
+{-# INLINE (>!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '>=' y@
+--
+-- >>> (20 ... 40 :: Interval Double) >=! (10 ... 20 :: Interval Double)
+-- True
+--
+-- >>> (5 ... 20 :: Interval Double) >=! (15 ... 40 :: Interval Double)
+-- False
+(>=!) :: Ord a => Interval a -> Interval a -> Bool
+I ax _ >=! I _ by = ax >= by
+{-# INLINE (>=!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x `op` y@
+certainly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool
+certainly cmp l r
+    | lt && eq && gt = True
+    | lt && eq       = l <=! r
+    | lt &&       gt = l /=! r
+    | lt             = l <!  r
+    |       eq && gt = l >=! r
+    |       eq       = l ==! r
+    |             gt = l >!  r
+    | otherwise      = False
+    where
+        lt = cmp False True
+        eq = cmp True True
+        gt = cmp True False
+{-# INLINE certainly #-}
+
+-- | Check if interval @X@ totally contains interval @Y@
+--
+-- >>> (20 ... 40 :: Interval Double) `contains` (25 ... 35 :: Interval Double)
+-- True
+--
+-- >>> (20 ... 40 :: Interval Double) `contains` (15 ... 35 :: Interval Double)
+-- False
+contains :: Ord a => Interval a -> Interval a -> Bool
+contains (I ax bx) (I ay by) = ax <= ay && by <= bx
+{-# INLINE contains #-}
+
+-- | Flipped version of `contains`. Check if interval @X@ a subset of interval @Y@
+--
+-- >>> (25 ... 35 :: Interval Double) `isSubsetOf` (20 ... 40 :: Interval Double)
+-- True
+--
+-- >>> (20 ... 40 :: Interval Double) `isSubsetOf` (15 ... 35 :: Interval Double)
+-- False
+isSubsetOf :: Ord a => Interval a -> Interval a -> Bool
+isSubsetOf = flip contains
+{-# INLINE isSubsetOf #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<' y@?
+(<?) :: Ord a => Interval a -> Interval a -> Bool
+I ax _ <? I _ by = ax < by
+{-# INLINE (<?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<=' y@?
+(<=?) :: Ord a => Interval a -> Interval a -> Bool
+I ax _ <=? I _ by = ax <= by
+{-# INLINE (<=?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '==' y@?
+(==?) :: Ord a => Interval a -> Interval a -> Bool
+I ax bx ==? I ay by = ax <= by && bx >= ay
+{-# INLINE (==?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '/=' y@?
+(/=?) :: Eq a => Interval a -> Interval a -> Bool
+I ax bx /=? I ay by = ax /= by || bx /= ay
+{-# INLINE (/=?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>' y@?
+(>?) :: Ord a => Interval a -> Interval a -> Bool
+I _ bx >? I ay _ = bx > ay
+{-# INLINE (>?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>=' y@?
+(>=?) :: Ord a => Interval a -> Interval a -> Bool
+I _ bx >=? I ay _ = bx >= ay
+{-# INLINE (>=?) #-}
+
+-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x `op` y@?
+possibly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool
+possibly cmp l r
+    | lt && eq && gt = True
+    | lt && eq       = l <=? r
+    | lt &&       gt = l /=? r
+    | lt             = l <? r
+    |       eq && gt = l >=? r
+    |       eq       = l ==? r
+    |             gt = l >? r
+    | otherwise      = False
+    where
+        lt = cmp LT EQ
+        eq = cmp EQ EQ
+        gt = cmp GT EQ
+{-# INLINE possibly #-}
+
+-- | The nearest value to that supplied which is contained in the interval.
+clamp :: Ord a => Interval a -> a -> a
+clamp (I a b) x
+  | x < a     = a
+  | x > b     = b
+  | otherwise = x
+
+-- | id function. Useful for type specification
+--
+-- >>> :t idouble (1 ... 3)
+-- idouble (1 ... 3) :: Interval Double
+idouble :: Interval Double -> Interval Double
+idouble = id
+
+-- | id function. Useful for type specification
+--
+-- >>> :t ifloat (1 ... 3)
+-- ifloat (1 ... 3) :: Interval Float
+ifloat :: Interval Float -> Interval Float
+ifloat = id
+
+-- Bugs:
+-- sin 1 :: Interval Double
+
+default (Integer,Double)