diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -1,3 +1,7 @@
+0.4
+---
+* Distributive Interval
+
 0.3
 ---
 * Removed dependency on `numeric-extras`
diff --git a/Setup.lhs b/Setup.lhs
--- a/Setup.lhs
+++ b/Setup.lhs
@@ -1,3 +1,55 @@
-#!/usr/bin/env runhaskell
-> import Distribution.Simple
-> main = defaultMainWithHooks simpleUserHooks
+#!/usr/bin/runhaskell
+\begin{code}
+{-# OPTIONS_GHC -Wall #-}
+module Main (main) where
+
+import Data.List ( nub )
+import Data.Version ( showVersion )
+import Distribution.Package ( PackageName(PackageName), Package, PackageId, InstalledPackageId, packageVersion, packageName )
+import Distribution.PackageDescription ( PackageDescription(), TestSuite(..) )
+import Distribution.Simple ( defaultMainWithHooks, UserHooks(..), simpleUserHooks )
+import Distribution.Simple.Utils ( rewriteFile, createDirectoryIfMissingVerbose, copyFiles )
+import Distribution.Simple.BuildPaths ( autogenModulesDir )
+import Distribution.Simple.Setup ( BuildFlags(buildVerbosity), Flag(..), fromFlag, HaddockFlags(haddockDistPref))
+import Distribution.Simple.LocalBuildInfo ( withLibLBI, withTestLBI, LocalBuildInfo(), ComponentLocalBuildInfo(componentPackageDeps) )
+import Distribution.Text ( display )
+import Distribution.Verbosity ( Verbosity, normal )
+import System.FilePath ( (</>) )
+
+main :: IO ()
+main = defaultMainWithHooks simpleUserHooks
+  { buildHook = \pkg lbi hooks flags -> do
+     generateBuildModule (fromFlag (buildVerbosity flags)) pkg lbi
+     buildHook simpleUserHooks pkg lbi hooks flags
+  , postHaddock = \args flags pkg lbi -> do
+     copyFiles normal (haddockOutputDir flags pkg) []
+     postHaddock simpleUserHooks args flags pkg lbi
+  }
+
+haddockOutputDir :: Package p => HaddockFlags -> p -> FilePath
+haddockOutputDir flags pkg = destDir where
+  baseDir = case haddockDistPref flags of
+    NoFlag -> "."
+    Flag x -> x
+  destDir = baseDir </> "doc" </> "html" </> display (packageName pkg)
+
+generateBuildModule :: Verbosity -> PackageDescription -> LocalBuildInfo -> IO ()
+generateBuildModule verbosity pkg lbi = do
+  let dir = autogenModulesDir lbi
+  createDirectoryIfMissingVerbose verbosity True dir
+  withLibLBI pkg lbi $ \_ libcfg -> do
+    withTestLBI pkg lbi $ \suite suitecfg -> do
+      rewriteFile (dir </> "Build_" ++ testName suite ++ ".hs") $ unlines
+        [ "module Build_" ++ testName suite ++ " where"
+        , "deps :: [String]"
+        , "deps = " ++ (show $ formatdeps (testDeps libcfg suitecfg))
+        ]
+  where
+    formatdeps = map (formatone . snd)
+    formatone p = case packageName p of
+      PackageName n -> n ++ "-" ++ showVersion (packageVersion p)
+
+testDeps :: ComponentLocalBuildInfo -> ComponentLocalBuildInfo -> [(InstalledPackageId, PackageId)]
+testDeps xs ys = nub $ componentPackageDeps xs ++ componentPackageDeps ys
+
+\end{code}
diff --git a/intervals.cabal b/intervals.cabal
--- a/intervals.cabal
+++ b/intervals.cabal
@@ -1,5 +1,5 @@
 name:              intervals
-version:           0.3
+version:           0.4
 synopsis:          Interval Arithmetic
 description:
   A 'Numeric.Interval.Interval' is a closed, convex set of floating point values.
@@ -15,15 +15,24 @@
 author:             Edward Kmett
 maintainer:         ekmett@gmail.com
 category:           Math
-build-type:         Simple
-cabal-version:      >=1.6
-extra-source-files: .travis.yml CHANGELOG.markdown README.markdown HLint.hs
+build-type:         Custom
+cabal-version:      >=1.8
 tested-with:        GHC == 7.4.2, GHC == 7.6.1, GHC == 7.6.3
+extra-source-files:
+  .travis.yml
+  CHANGELOG.markdown
+  README.markdown
+  HLint.hs
 
 source-repository head
   type: git
   location: git://github.com/ekmett/intervals.git
 
+-- You can disable the doctests test suite with -f-test-doctests
+flag test-doctests
+  default: True
+  manual: True
+
 library
   hs-source-dirs: src
 
@@ -31,9 +40,29 @@
 
   build-depends:
     array          >= 0.3   && < 0.6,
-    base           >= 4     && < 5
+    base           >= 4     && < 5,
+    distributive   >= 0.2   && < 1
 
   if impl(ghc >=7.4)
     build-depends: ghc-prim
 
   ghc-options:     -Wall -O2
+
+
+test-suite doctests
+  type:           exitcode-stdio-1.0
+  main-is:        doctests.hs
+  ghc-options:    -Wall -threaded
+  hs-source-dirs: tests
+
+  if !flag(test-doctests)
+    buildable: False
+  else
+    build-depends:
+      base,
+      directory      >= 1.0,
+      doctest        >= 0.9.1,
+      filepath
+
+  if impl(ghc<7.6.1)
+    ghc-options: -Werror
diff --git a/src/Numeric/Interval.hs b/src/Numeric/Interval.hs
--- a/src/Numeric/Interval.hs
+++ b/src/Numeric/Interval.hs
@@ -46,6 +46,7 @@
 
 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
@@ -55,6 +56,8 @@
 #endif
 import Prelude hiding (null, elem, notElem)
 
+-- $setup
+
 data Interval a = I !a !a deriving
   ( Data
   , Typeable
@@ -92,6 +95,10 @@
     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
@@ -117,31 +124,55 @@
 {-# 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 #-}
@@ -149,6 +180,12 @@
 -- | 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 #-}
@@ -165,16 +202,43 @@
       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 #-}
@@ -202,27 +266,72 @@
   {-# 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 #-}
 
--- | This means that realToFrac will use the midpoint
-
--- | What moron put an Ord instance requirement on Real!
+-- | 'realToFrac' will use the midpoint
 instance Real a => Real (Interval a) where
   toRational x
     | null x   = nan
@@ -389,7 +498,7 @@
 increasing :: (a -> b) -> Interval a -> Interval b
 increasing f (I a b) = f a ... f b
 
--- | lift a monotone increasing function over a given interval
+-- | 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
 
@@ -426,6 +535,9 @@
 -- 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
@@ -433,6 +545,12 @@
 {-# 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
@@ -441,36 +559,82 @@
 {-# 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
@@ -487,11 +651,25 @@
         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 #-}
@@ -543,11 +721,22 @@
         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/tests/doctests.hsc b/tests/doctests.hsc
new file mode 100644
--- /dev/null
+++ b/tests/doctests.hsc
@@ -0,0 +1,73 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE ForeignFunctionInterface #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Main (doctests)
+-- Copyright   :  (C) 2012-13 Edward Kmett
+-- License     :  BSD-style (see the file LICENSE)
+-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- This module provides doctests for a project based on the actual versions
+-- of the packages it was built with. It requires a corresponding Setup.lhs
+-- to be added to the project
+-----------------------------------------------------------------------------
+module Main where
+
+import Build_doctests (deps)
+import Control.Applicative
+import Control.Monad
+import Data.List
+import System.Directory
+import System.FilePath
+import Test.DocTest
+
+##if defined(mingw32_HOST_OS)
+##if defined(i386_HOST_ARCH)
+##define USE_CP
+import Control.Applicative
+import Control.Exception
+import Foreign.C.Types
+foreign import stdcall "windows.h SetConsoleCP" c_SetConsoleCP :: CUInt -> IO Bool
+foreign import stdcall "windows.h GetConsoleCP" c_GetConsoleCP :: IO CUInt
+##elif defined(x86_64_HOST_ARCH)
+##define USE_CP
+import Control.Applicative
+import Control.Exception
+import Foreign.C.Types
+foreign import ccall "windows.h SetConsoleCP" c_SetConsoleCP :: CUInt -> IO Bool
+foreign import ccall "windows.h GetConsoleCP" c_GetConsoleCP :: IO CUInt
+##endif
+##endif
+
+-- | Run in a modified codepage where we can print UTF-8 values on Windows.
+withUnicode :: IO a -> IO a
+##ifdef USE_CP
+withUnicode m = do
+  cp <- c_GetConsoleCP
+  (c_SetConsoleCP 65001 >> m) `finally` c_SetConsoleCP cp
+##else
+withUnicode m = m
+##endif
+
+main :: IO ()
+main = withUnicode $ getSources >>= \sources -> doctest $
+    "-isrc"
+  : "-idist/build/autogen"
+  : "-optP-include"
+  : "-optPdist/build/autogen/cabal_macros.h"
+  : "-hide-all-packages"
+  : map ("-package="++) deps ++ sources
+
+getSources :: IO [FilePath]
+getSources = filter (isSuffixOf ".hs") <$> go "src"
+  where
+    go dir = do
+      (dirs, files) <- getFilesAndDirectories dir
+      (files ++) . concat <$> mapM go dirs
+
+getFilesAndDirectories :: FilePath -> IO ([FilePath], [FilePath])
+getFilesAndDirectories dir = do
+  c <- map (dir </>) . filter (`notElem` ["..", "."]) <$> getDirectoryContents dir
+  (,) <$> filterM doesDirectoryExist c <*> filterM doesFileExist c