diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,31 @@
+Copyright (c) 2010, Edward Kmett
+
+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 Edward Kmett 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/Numeric/Interval.hs b/Numeric/Interval.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/Interval.hs
@@ -0,0 +1,485 @@
+{-# LANGUAGE Rank2Types, TypeFamilies #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Interval
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- Interval arithmetic
+--
+-----------------------------------------------------------------------------
+
+module Numeric.Interval 
+    ( whole
+    , empty
+    , null
+    , singleton
+    , elem
+    , notElem
+    , inf
+    , sup
+    , singular
+    , width
+    , intersection
+    , hull
+    , bisection
+    , magnitude
+    , mignitude
+    , contains
+    , isSubsetOf
+    , certainly, (<!), (<=!), (==!), (>=!), (>!)
+    , possibly, (<?), (<=?), (==?), (>=?), (>?)
+    , idouble 
+    , ifloat 
+    ) where
+
+import Prelude hiding (null, elem, notElem)
+import Numeric.Extras
+import Numeric.Rounding
+import Data.Function (on)
+
+data Interval a = I !(Round Down a) !(Round Up a)
+
+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
+
+-- | 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 
+a ... b = I (Round a) (Round b)
+{-# INLINE (...) #-}
+
+-- | The whole real number line
+whole :: Precision a => Interval a 
+whole = negInfinity ... posInfinity
+{-# INLINE whole #-}
+
+-- | An empty interval
+empty :: Precision a => Interval a 
+empty = nan ... nan
+{-# INLINE empty #-}
+
+-- | negation handles NaN properly
+null :: Ord a => Interval a -> Bool
+null x = not (inf x <= sup x)
+{-# INLINE null #-}
+
+-- | A singleton point
+singleton :: a -> Interval a 
+singleton a = a ... a
+{-# INLINE singleton #-}
+
+-- | The infinumum (lower bound) of an interval
+inf :: Interval a -> a
+inf (I (Round a) _) = a
+{-# INLINE inf #-}
+
+-- | The supremum (upper bound) of an interval
+sup :: Interval a -> a
+sup (I _ (Round 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 :: Ord a => Interval a -> Bool
+singular x = not (null x) && inf x == sup x
+{-# INLINE singular #-}
+
+instance Precision a => Eq (Interval a) where
+    (==) = (==) `on` midpoint
+
+instance Show a => Show (Interval a) where
+    showsPrec n (I (Round a) (Round b)) =   
+        showParen (n > 3) $
+            showsPrec 3 a .
+            showString " ... " . 
+            showsPrec 3 b
+
+-- flip the rounding mode up
+u :: Round Down a -> Round Up a
+u (Round a) = Round a
+{-# INLINE u #-}
+
+-- flip the rounding mode down
+d :: Round Up a -> Round Down a
+d (Round a) = Round a
+{-# INLINE d #-}
+
+-- | Calculate the width of an interval.
+-- N.B. the width of an interval is an interval itself due to rounding
+width :: Precision a => Interval a -> Interval a
+width (I a b) = I (d b - a) (b - u a)
+{-# INLINE width #-}
+
+-- | magnitude 
+magnitude :: Precision a => Interval a -> a 
+magnitude x = (max `on` abs) (inf x) (sup x)
+{-# INLINE magnitude #-}
+
+-- | "mignitude"
+mignitude :: Precision a => Interval a -> a 
+mignitude x = (min `on` abs) (inf x) (sup x)
+{-# INLINE mignitude #-}
+
+instance Precision a => Num (Interval a) where
+    I a b + I a' b' = I (a + a') (b + b')
+    I a b - I a' b' = I (a - d b') (b - u a')
+    I a b * I a' b' = minimum [a * a',a * d b',d b * a',d b * d b'] 
+                      `I` 
+                      maximum [u a * u a',u a * b',b * u a',b * b']
+    abs x@(I a b) 
+        | a >= 0    = x 
+        | b <= 0    = negate x
+        | otherwise = max (- a) (d b) `I` b
+
+    signum (I a b)  = signum a `I` signum b
+
+    fromInteger i   = fromInteger i `I` fromInteger i
+
+-- | Bisect an interval at its midpoint.
+bisection :: Precision a => Interval a -> (Interval a, Interval a)
+bisection (I a b) = (I a (u a + (b - u a) / 2), I (a + (d b - a) / 2) b)
+{-# INLINE bisection #-}
+
+-- | Nearest point to the midpoint of the interval.
+midpoint :: Precision a => Interval a -> a
+midpoint x = inf x + (sup x - inf x) / 2
+{-# INLINE midpoint #-}
+
+elem :: Precision a => a -> Interval a -> Bool
+elem x xs = x >= inf xs && x <= sup xs
+{-# INLINE elem #-}
+
+notElem :: Precision a => a -> Interval a -> Bool
+notElem x xs = not (elem x xs)
+{-# INLINE notElem #-}
+
+-- | This means that realToFrac will use the midpoint
+instance Precision a => Real (Interval a) where
+    toRational x = toRational (midpoint x)
+
+-- @'divNonZero' X Y@ assumes @0 `'notElem'` Y@
+divNonZero :: Precision a => Interval a -> Interval a -> Interval a 
+divNonZero (I a b) (I a' b') = 
+    minimum [a / a',a / d b',d b / a',d b / d b'] 
+    `I`
+    maximum [u a / u a',u a / b',b / u a',b / b']
+
+-- @'divPositive' X y@ assumes y > 0, and divides @X@ by [0 ... y]
+divPositive :: Precision a => Interval a -> a -> Interval a 
+divPositive x@(I a b) y
+    | a == 0 && b == 0 = x
+    | b < 0 || isNegativeZero b = negInfinity `I` ( b / up y)
+    | a < 0 = whole 
+ -- | isNegativeZero a = whole
+    | otherwise = (a / down y) `I` posInfinity
+
+-- divNegative assumes y < 0 and divides the interval @X@ by [y ... 0]
+divNegative :: Precision a => Interval a -> a -> Interval a
+divNegative x@(I a b) y
+    | a == 0 && b == 0 = - x -- flips negative zeros
+    | b < 0 || isNegativeZero b = (d b / down y) `I` posInfinity
+    | a < 0 = whole
+ -- | isNegativeZero a = whole
+    | otherwise = negInfinity `I` (u a / up y)
+
+divZero :: Precision a => Interval a -> Interval a
+divZero x | inf x == 0 && sup x == 0 = x
+          | otherwise = whole
+
+instance Precision 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 (d b) `I` on max recip (u a) b
+    fromRational r  = fromRational r `I` fromRational r
+
+instance Precision a => RealFrac (Interval a) where
+    properFraction x = (b, x - fromIntegral b)
+        where 
+            b = truncate (midpoint x)
+    ceiling x = ceiling (sup x)
+    floor x = floor (inf x)
+    round x = round (midpoint x)
+    truncate x = truncate (midpoint x)
+
+instance Precision a => Floating (Interval a) where
+    pi = pi `I` pi 
+    exp = increasing exp
+    log (I a b) = (if a > 0 then log a else negInfinity) `I` log b
+    cos x 
+        | null x = empty
+        | inf (width t) >= inf pi = (-1) ... 1
+        | tl >= d pih  = - cos (t - pi)
+        | th <= u pil  = cos (d th) `I` cos (u tl)
+        | th <= u pi2l = (-1) `I` cos (u (min (pi2l - d th) tl))
+        | otherwise  = (-1) ... 1
+        where 
+            I pil pih = pi
+            pi2@(I pi2l _) = pi * 2
+            t@(I tl th) = x `fmod` pi2
+    sin x 
+        | null x = empty
+        | otherwise = cos (x - pi / 2)
+    tan x 
+        | null x = empty
+        | inf t' <= -hpil || sup t' >= hpil = whole
+        | otherwise = increasing tan x
+        where
+            t = x `fmod` pi 
+            t' | inf t >= hpil = t - pi
+               | otherwise = t
+            hpil = inf (pi / 2)
+    asin x@(I a b)
+        | null x || b < -1 || a > 1 = empty
+        | otherwise = 
+            (if a <= - 1 then - d hpis else asin a)
+            `I`
+            (if b >= 1 then hpis else asin b)
+        where
+            I _ hpis = pi / 2
+    acos x@(I a b)
+        | null x || b < -1 || a > 1 = empty
+        | otherwise = 
+            (if b >= 1 then 0 else acos (d b))
+            `I`
+            (if a < -1 then pis else acos (u a))
+        where
+            I _ pis = pi
+    atan = increasing atan
+    sinh = increasing sinh
+    cosh x@(I a b)
+        | null x = empty
+        | b < 0  = decreasing cosh x
+        | a >= 0 = increasing cosh x
+        | otherwise  = I 0 $ cosh $ if - u a > b
+                                    then u a 
+                                    else b
+    tanh = increasing tanh
+    asinh = increasing asinh
+    acosh x@(I a b)
+        | null x || b < 1 = empty -- acosh is only defined on [1..1/0)
+        | otherwise = I lo $ acosh b
+        where lo | a <= 1 = 0 
+                 | otherwise = acosh a
+    atanh x@(I a b)
+        | null x || b < -1 || a > 1 = empty
+        | otherwise =
+                (if a <= - 1 then negInfinity else atanh a)
+                `I`
+                (if b >= 1 then posInfinity else atanh b)
+    
+-- | lift a monotone increasing function over a given interval 
+increasing :: Precision a => 
+         (forall d. Rounding d => Round d a -> Round d a) -> 
+         Interval a -> Interval a
+increasing f (I a b) = I (f a) (f b)
+
+-- | lift a monotone increasing function over a given interval 
+decreasing :: Precision a => 
+         (forall d. Rounding d => Round d a -> Round d a) -> 
+         Interval a -> Interval a
+decreasing f (I a b) = I (f (d b)) (f (u a))
+
+
+-- | We have to play some semantic games to make these methods make sense.
+-- Most compute with the midpoint of the interval.
+instance Precision 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 :: Precision 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')
+
+-- | Calculate the convex hull of two intervals
+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 = I (min a a') (max b b')
+    
+instance Precision a => RealExtras (Interval a) where
+    type C (Interval a) = C a
+    -- output always lies within the interval y if y >=! 0
+    fmod x y | null y = empty 
+             | otherwise = r -- `intersection` bounds
+        where 
+            n :: Integer
+            n = floor (inf x / if inf x < 0 then inf y else sup y)
+            r = x - fromIntegral n * y 
+            -- bounds | inf y >= 0 = y
+            --        | otherwise = y `hull` negate y
+    expm1 = increasing expm1
+    log1p (I a b) = (if a > (-1) then log1p a else negInfinity) `I` log1p b
+    hypot x y = hypot a a' `I` hypot b b'
+        where
+            I a b = abs x
+            I a' b' = abs y
+    cbrt = increasing cbrt
+    erf = increasing erf
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '<' y@
+(<!)  :: Ord a => Interval a -> Interval a -> Bool
+x <! y = sup x < inf y
+{-# INLINE (<!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '<=' y@
+(<=!) :: Ord a => Interval a -> Interval a -> Bool
+x <=! y = sup x <= inf y
+{-# INLINE (<=!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '==' y@
+(==!) :: 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@
+(/=!) :: 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@
+(>!)  :: Ord a => Interval a -> Interval a -> Bool
+x >! y = inf x > sup y
+{-# INLINE (>!) #-}
+
+-- | For all @x@ in @X@, @y@ in @Y@. @x '>=' y@
+(>=!) :: 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 #-}
+
+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 #-}
+
+isSubsetOf :: Ord a => Interval a -> Interval a -> Bool
+isSubsetOf = flip contains
+
+-- | Comparisons are made on the midpoint
+instance Precision a => Ord (Interval a) where
+    compare = compare `on` midpoint
+    max (I a b) (I a' b') = I (max a a') (max b b')
+    min (I a b) (I a' b') = I (min a a') (min b b')
+
+-- | 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 #-}
+
+idouble :: Interval Double -> Interval Double
+idouble = id
+
+ifloat :: Interval Float -> Interval Float
+ifloat = id
+
+-- Bugs:
+-- sin 1 :: Interval Double
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#!/usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMainWithHooks simpleUserHooks
diff --git a/intervals.cabal b/intervals.cabal
new file mode 100644
--- /dev/null
+++ b/intervals.cabal
@@ -0,0 +1,19 @@
+Name:              intervals
+Version:           0.1.0
+Synopsis:          Interval Arithmetic
+Homepage:          http://patch-tag.com/r/ekmett/intervals
+License:           BSD3
+License-file:      LICENSE
+Author:            Edward Kmett
+Maintainer:        ekmett@gmail.com
+Category:          Math
+Build-type:        Simple
+Cabal-version:     >=1.6
+
+Library
+  Exposed-modules: Numeric.Interval
+  Build-depends:   base >= 4 && < 5,
+                   array >= 0.3.0 && < 0.4,
+                   numeric-extras >= 0.0.1 && < 0.1,
+                   rounding >= 0.3.0 && < 0.4
+  GHC-Options:     -Wall