4 files changed, 1345 insertions, 0 deletions

diff --git a/lib/CommandLine.hs b/lib/CommandLine.hs
new file mode 100644
index 0000000..dfc16f9
--- /dev/null
+++ b/lib/CommandLine.hs
@@ -0,0 +1,559 @@
+{-# LANGUAGE DeriveFunctor #-}
+{-# LANGUAGE CPP #-}
+-- {-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE StandaloneDeriving #-}
+{-# LANGUAGE PatternGuards #-}
+module CommandLine
+    ( Args
+    , UsageError(..)
+    , usageErrorMessage
+    , parseInvocation
+    , runArgs
+    , arg
+    , param
+    , params
+    , label
+    ) where
+
+import Control.Applicative
+import Control.Arrow
+import Control.Monad
+import Data.Bits
+import Data.Either
+import Data.Function
+import Data.List
+import Data.Maybe
+import Data.Ord
+import Data.Map.Strict (Map)
+import qualified Data.Map.Strict as Map
+import Data.IntMap.Strict (IntMap)
+import qualified Data.IntMap.Strict as IntMap
+import Debug.Trace
+import Numeric.Interval (Interval(..), singleton, (...), inf, sup, hull)
+import qualified Numeric.Interval as I
+import Numeric.Interval.Bounded
+import SuperOrd
+
+-- trace :: String -> a -> a
+-- trace _ x = x
+
+-- type CompF a = [String] -> [String] -> a
+
+type MergeData = [(Int,Ordering)]
+
+-- | Expr a
+--
+data Expr a where
+    -- | Prim
+    --
+    --   Takes a function from the option arguments and unamed arguments repsectively to
+    --   a value of type a, usually IO (), and gives you an expression tree. As one
+    --   traverses down the tree only the 'interesting' option arguments are passed
+    --   to this function, but all of the unnamed arguments are passed regardless of
+    --   where we are in the tree.
+    --
+    Prim :: ([[String]] -> [String] -> a) -> Interval (SuperOrd Int) -> Expr a
+    -- | Star
+    --   Applicative '<*>'
+    Star :: MergeData -> Expr (b -> a) -> (Expr b) -> Expr a
+    -- | Or
+    --   Alternative '<|>'
+    Or :: MergeData -> Expr a -> Expr a -> Expr a
+    -- | Empty
+    --   Alternative empty
+    Empty :: Expr a
+
+deriving instance Functor Expr
+
+-- | Args
+--
+--   Applicative Functor for interpretting command line arguments.
+data Args a = Args
+    { expr :: Expr a
+      -- ^ Expression tree
+    , accepts :: [String]
+      -- ^ sorted list of acceptable short and long option names (non positional arguments)
+      --   The names include hyphens.
+    }
+ deriving Functor
+
+instance Applicative Args where
+    pure x = Args { expr = Prim (\_ _ -> x) (singleton $ exactly 0), accepts = [] }
+    f <*> b = Args
+        { expr = Star d (expr f) (expr b)
+        , accepts = m
+        }
+        where d = mergeData compare (accepts f) (accepts b)
+              m = mergeLists d const (accepts f) (accepts b)
+
+instance Alternative Args where
+    empty = Args Empty []
+    f <|> g = Args
+        { expr = Or d (expr f) (expr g)
+        , accepts = m
+        }
+        where d = mergeData compare (accepts f) (accepts g)
+              m = mergeLists d const (accepts f) (accepts g)
+
+
+{- dead code?
+unpackBits :: Integer -> [Bool]
+unpackBits 0 = [False]
+unpackBits 1 = [True]
+unpackBits n = ( r /= 0 ) : unpackBits q
+ where
+    (q,r) = divMod n 2
+
+-- requires finite list
+packBits :: [Bool] -> Integer
+packBits bs = sum $ zipWith (\b n -> if b then n else 0) bs $ iterate (*2) 1
+ -}
+
+
+-- | mergeData
+--
+--   > mergeData compare [1,3,5] [2,2,4,6]  ==> [(1,LT),(2,GT),(1,LT),(1,GT),(1,LT),(1,GT)]
+--
+--   Given a comparison function and two sorted lists, 'mergeData' will return
+--   a RLE compressed (run-length encoded) list of the comparison results
+--   encountered while merging the lists.
+--
+--   This data is enough information to perform the merge without doing the
+--   comparisons or to reverse a merged list back to two sorted lists.
+--
+--   When one list is exausted, the length of the remaining list is returned as
+--   a run-length for LT or GT depending on whether the left list or the right
+--   list has elements.
+mergeData :: (a -> a -> Ordering) -> [a] -> [a] -> [(Int,Ordering)]
+mergeData comp (x:xs) (y:ys)
+  | comp x y == LT = case mergeData comp xs (y:ys) of
+                        (n,LT):ys -> let n'=n+1 in n' `seq` (n',LT):ys
+                        ys        -> (1,LT):ys
+  | comp x y == EQ = case mergeData comp xs ys of
+                        (n,EQ):ys -> let n'=n+1 in n' `seq` (n',EQ):ys
+                        ys        -> (1,EQ):ys
+  | comp x y == GT = case mergeData comp (x:xs) ys of
+                        (n,GT):ys -> let n'=n+1 in n' `seq` (n',GT):ys
+                        ys        -> (1,GT):ys
+mergeData comp [] [] = []
+mergeData comp [] ys = (length ys, GT) : []
+mergeData comp xs [] = (length xs, LT) : []
+
+mergeLists :: [(Int,Ordering)] -> (a -> a -> a) -> [a] -> [a] -> [a]
+mergeLists ((n,LT):os) f xs ys = ls ++ mergeLists os f xs' ys
+ where
+    (ls,xs') = splitAt n xs
+mergeLists ((n,EQ):os) f xs ys = es ++ mergeLists os f xs' ys'
+ where
+    (les,xs') = splitAt n xs
+    (res,ys') = splitAt n ys
+    es = zipWith f les res
+mergeLists ((n,GT):os) f xs ys = gs ++ mergeLists os f xs ys'
+ where
+    (gs,ys') = splitAt n ys
+mergeLists [] f [] ys = ys
+mergeLists [] f xs [] = xs
+mergeLists [] f xs ys = error "xs ++ ys"
+
+{-
+computeMask :: Int -> Ordering -> Ordering -> [(Int,Ordering)] -> Integer
+computeMask k w t [] = 0
+computeMask k w t ((n,v):os)
+   = if w==v then r .|. shiftL (bit n - 1) k
+             else r
+ where r = computeMask (k+n') w t os
+       n' | v==t      = n
+          | otherwise = 0
+
+-- WRONG, one-blocks are not spaced the same in input and output, need shifts
+mergeIntegers :: [(Int,Ordering)] -> (Integer -> Integer -> Integer) -> Integer -> Integer -> Integer
+mergeIntegers os f x y = (leftmask .&. x) .|. (rightmask .&. y) .|. (f (leqmask .&. x) (reqmask .&. y))
+ where
+    leftmask  = computeMask 0 LT EQ os
+    leqmask   = computeMask 0 EQ LT os
+    reqmask   = computeMask 0 EQ GT os
+    rightmask = computeMask 0 GT EQ os
+-}
+{- kinda dead code
+mergeIntegers ((n,LT):os) f !x !y = v `seq` m `seq` m .|. v
+ where
+    m = x .&. (bit n - 1)
+    r = mergeIntegers os f (shiftR x n) y
+    v = r `seq` shiftL r n
+mergeIntegers ((n,EQ):os) f !x !y = mm `seq` v `seq` mm .|. v
+ where
+    mm = f mx my
+    mx = x .&. (bit n - 1)
+    my = y .&. (bit n - 1)
+    r = mergeIntegers os f (shiftR x n) (shiftR y n)
+    v = r `seq` shiftL r n
+mergeIntegers ((n,GT):os) f !x !y = v `seq` m `seq` m .|. v
+ where
+    m = y .&. (bit n - 1)
+    r = mergeIntegers os f x (shiftR y n)
+    v = r `seq` shiftL r n
+mergeIntegers [] f !0 !y = y
+mergeIntegers [] f !x !0 = x
+mergeIntegers [] f !x !y = error "x .|. y"
+-}
+
+splitLists :: [(Int,Ordering)] -> [a] -> ([a],[a])
+splitLists ((n,LT):os) xs = (ls ++ lls, rrs)
+ where
+    (ls,xs') = splitAt n xs
+    (lls,rrs) = splitLists os xs'
+splitLists ((n,EQ):os) xs = (es ++ lls, es ++ rrs)
+ where
+    (es,xs') = splitAt n xs
+    (lls,rrs) = splitLists os xs'
+splitLists ((n,GT):os) xs = (lls, rs ++ rrs)
+ where
+    (rs,xs') = splitAt n xs
+    (lls,rrs) = splitLists os xs'
+splitLists [] xs = (xs,xs)
+
+{-
+mergeBy :: Show a => (a -> a -> Ordering) -> [a] -> [a]
+    -> ( (Integer -> Integer -> Integer) -> Integer -> Integer -> Integer
+       , (b -> b -> b) -> [b] -> [b] -> [b]
+       , [b] -> ([b], [b]))
+mergeBy comp xs ys = trace (unlines ["xs="++show xs,"ys="++show ys,"mergeData="++show d]) (mergeIntegers d, mergeLists d, splitLists d)
+ where
+    d = mergeData comp xs ys
+-}
+
+
+param :: Int -> Args String
+param n = Args (Prim (\_ us -> us!!n) (singleton $ exactly (n+1))) []
+
+arg :: String -> Args String
+arg optname = Args (Prim (\opts _ -> head $ concat $ take 1 opts)
+                   (singleton $ exactly 0))
+                   [optname]
+
+params :: Args [String]
+params = Args (Prim (\_ args -> args) (exactly 0 ... PositiveInfinity)) []
+
+
+label :: String -> Args a -> Args a
+label _ = id
+
+data ArgsStyle = ArgsStyle
+    { parseInvocation :: String -> [String] -> ([(String,[String])], [String])
+    }
+
+-- | Phase 1.  This function accepts a list of command line arguments in its
+-- second argument that will be parsed to obtain a pair of lists: named
+-- argument-value pairs and unnamed arguments values.
+--
+-- The first argument indicates which short options will be treated as on/off
+-- flags and given a dummy value of \"\".  Hyphen-prefixed options not in this
+-- list are given their imeediately following string as a value.
+--
+-- The \"--\" convention is implemented to indicate all remaining arguments are
+-- unnamed argument values.
+--
+-- The results of this function are intended to be used as input to 'runArgs'.
+vanilla :: ArgsStyle
+vanilla = ArgsStyle { parseInvocation = parse }
+ where
+  parse flags cli = (opts, concat nons ++ nondashed ++ drop 1 trailing)
+   where
+    (as, trailing)  = span (/= "--") cli
+    (nons,bs)       = span ((==[True]) . map (not . isPrefixOf "-") . take 1) $ groupBy (const $ not . isPrefixOf "-") as
+    (ds, nondashed) = second concat $ unzip $ map splitParams bs
+    opts            = map ((first concat) . splitAt 1) (ds :: [[String]])
+
+    splitParams (('-':[x]):xs) | x `elem` flags = (['-':[x]],xs)
+    splitParams xs = splitAt 2 xs
+
+-- | Information about how the user failed to satisfy a specified usage.
+data UsageError
+    = TooManyParameters Int
+      -- ^ The given number of excessive unnamed arguments occured.
+    | InsufficientParameters Int
+      -- ^ Not enough unnamed arguments.  The number indicates how many are
+      -- total are expected.
+    | TooManyOf String [String]
+      -- ^ An option was supplied too many times.  The list is a set of values
+      -- associated with the repeated instances.
+    | Missing [String]
+      -- ^ A list of required options that the user failed to specify.
+    | ChooseOne [[String]]
+      -- ^ The user failed to choose one of the given set of option combinations.
+    | Misunderstood [String]
+      -- ^ A list of unrecognized options.
+    | Incompatible [[String]]
+      -- ^ A list of supplied options that may not be used together.
+
+    | NamedFailure String UsageError
+      -- ^ Extra context provided via the 'label' primitive.
+
+ deriving (Eq,Show)
+
+-- | Obtain a description of a usage error that can be reported to the user.
+usageErrorMessage :: UsageError -> String
+usageErrorMessage (NamedFailure _ e)         = usageErrorMessage e
+usageErrorMessage (TooManyParameters _)      = "too many arguments"
+usageErrorMessage (InsufficientParameters c) = "insufficient arguments (need "++show c++")"
+usageErrorMessage (TooManyOf n xs)           = n ++" can be specified only once"
+usageErrorMessage (Missing ns)               = "missing: "++intercalate ", " ns
+usageErrorMessage (ChooseOne nss)            = "choose one of: "++intercalate ", " (map (intercalate " ") nss)
+usageErrorMessage (Misunderstood ns)         = "unrecognized: "++intercalate ", " ns
+usageErrorMessage (Incompatible nss)         = intercalate " and " (map (intercalate " ") nss) ++ " cannot be used together"
+
+{-
+rankError :: UsageError -> Int
+rankError (NamedFailure _ e)         = rankError e
+rankError (TooManyParameters _)      = 0
+rankError (InsufficientParameters _) = 1
+rankError (TooManyOf _ xs)           = 1
+rankError (Missing _)                = 2
+rankError (ChooseOne _)              = 2
+rankError (Misunderstood xs)         = 2 + length xs
+rankError (Incompatible ys)          = 2 + length ys
+
+tagError :: UsageError -> Int
+tagError (NamedFailure _ _)         = 0
+tagError (TooManyParameters _)      = 1
+tagError (InsufficientParameters _) = 2
+tagError (TooManyOf _ _)            = 3
+tagError (Missing _)                = 4
+tagError (ChooseOne _)              = 5
+tagError (Misunderstood _)          = 6
+tagError (Incompatible _)           = 7
+
+missingWhat :: UsageError -> [[String]]
+missingWhat (Missing xs)       = [xs]
+missingWhat (ChooseOne ys)     = ys
+missingWhat (NamedFailure _ e) = missingWhat e
+missingWhat _                  = []
+
+misunderstoodWhat :: UsageError -> [String]
+misunderstoodWhat (Misunderstood xs) = xs
+misunderstoodWhat (Incompatible yss) = concatMap (take 1) yss
+misunderstoodWhat (NamedFailure _ e) = misunderstoodWhat e
+misunderstoodWhat _                  = []
+-}
+
+{- dead code
+tryCompute :: [(String,String)] -> [String] -> Computation a -> Either UsageError a
+tryCompute os us c@(Computation { compLabel = lbl })
+    | null lbl = tryCompute' os us c
+    | otherwise = either (Left . NamedFailure lbl) Right $ tryCompute' os us c
+ where
+    tryCompute' os us c
+        | not (null unused_os) = Left $ Misunderstood $ map fst unused_os
+        | not (null missing)   = Left $ Missing missing
+        | not (null repss)     = Left $ TooManyOf (fst $ head $ head repss) (map snd $ head repss)
+        | ulen < clen          = Left $ InsufficientParameters clen
+        | ulen > clen          = Left $ TooManyParameters (ulen - clen)
+        | otherwise            = Right $ compute c os us
+     where
+        (found, missing) = partition (\k -> k `elem` map fst os) $ consumedOptions c
+        (used_os, unused_os) = partition (\(k,v) -> k `elem` consumedOptions c) os
+        ulen = length us
+        repss = filter (not . null . tail) $ groupBy ((==) `on` fst) $ sortBy (comparing fst) used_os
+        clen = case consumedParameters c of
+                -1 -> ulen
+                num -> num
+-}
+
+#if defined(VERSION_base)
+#if !MIN_VERSION_base(4,8,0)
+sortOn :: Ord b => (a -> b) -> [a] -> [a]
+sortOn f =
+  map snd . sortBy (comparing fst) . map (\x -> let y = f x in y `seq` (y, x))
+#endif
+#endif
+
+removeIntersection (x:xs) (y:ys)
+    | x == y = removeIntersection xs ys
+    | x < y  = first (x :) $ removeIntersection xs (y:ys)
+    | x > y  = second (y :) $ removeIntersection (x:xs) ys
+removeIntersection [] ys = ([],ys)
+removeIntersection xs [] = (xs,[])
+
+
+-- ordinary sorted list merge.
+mergeL :: Ord a => [a] -> [a] -> [a]
+mergeL as bs = mergeLists (mergeData compare as bs) const as bs
+
+-- | runArgs
+--
+--    (os,us)   - named arguments(options, name-value pairs), and unnamed arguments
+--    c         - expression tree (applicative functor)
+--
+-- returns either a UsageError or a computed result (perhaps IO ())
+--
+-- Evaluate the given set of named and unnamed values and return
+-- the computed result or else indicate bad usage.
+--
+-- 'parseInvocation' may be used on the results of 'System.Environment.getArgs'
+-- to obtain suitable input for this function.
+runArgs :: ([(String,[String])], [String]) -> Args a -> Either UsageError a
+runArgs (os,us) c
+    | not (null bads) = Left $ Misunderstood $ map fst bads
+    | not (null dups) = Left $ TooManyOf (fst $ head dups) (concat $ take 1 $ snd $ head dups) -- only reports first dup.
+    | otherwise = makeError $ compute (expr c) (zipWith const [0..] (accepts c)) us
+ where
+    os' = sortOn fst os
+    dups = mapMaybe notSingle $ groupBy ((==) `on` fst) (os' :: [(String,[String])])
+        where notSingle [x] = Nothing
+              notSingle ((k,v):xs)  = Just (k,v : map snd xs)
+    getbit = let r = Map.fromList $ zip (accepts c) [0..] in trace ("getbit = "++show r) r
+    goods :: [(Int,[String])]
+    (bads,goods) = let r = partitionEithers $ map f os' in trace ("(bads,goods)="++show r) r
+        where f (k,v) = case Map.lookup k getbit of
+                            Just b  -> Right (b,v)
+                            Nothing -> Left  (k,v)
+
+    valmap = IntMap.fromList goods
+    namemap = IntMap.fromList $ zip [0..] (accepts c)
+
+    vals = map snd goods
+    ulen = length us
+
+    makeError (_,Left e) = Left $ makeError' e
+     where
+        makeError' (Left xss)   = Incompatible $ map (map (namemap IntMap.!)) xss
+        makeError' (Right [xs]) = Missing $ map (namemap IntMap.!) xs
+        makeError' (Right xss)  = ChooseOne $ map (map (namemap IntMap.!)) xss
+    makeError (i,Right v)
+        | exactly ulen > sup i = Left $ TooManyParameters (ulen - superApprox (sup i))
+        | exactly ulen < inf i = Left $ InsufficientParameters (superApprox (inf i))
+        | otherwise            = Right v
+
+
+    -- On success, returns Right, otherwise:
+    --
+    --   * @ Left (Right xss) @ - xss indicates unspecified required named-options.
+    --
+    --   * @ Left (Left xss) @ - xss is a list of mutually-exclusive sets of specified options.
+    --
+    compute :: Expr a -> [Int] -> [String] -> (Interval (SuperOrd Int), Either (Either [[Int]] [[Int]]) a)
+    compute (Prim f i) opts us
+        | null es   = (i, Right $ f vals us )
+        | otherwise = (i, Left  $ Right [es])
+     where
+        (es,vals) = partitionEithers
+                    $ map (\k -> maybe (Left k) Right (k `IntMap.lookup` valmap)) opts
+    compute (Star d f b) opts us = (max (inf fi) (inf bi) ... max (sup fi) (sup bi), r)
+     where
+        r = case (fres,bres) of
+                (Right ff        , Right bb)        -> Right $ ff bb
+                (Left e          , Right _)         -> Left  e
+                (Right _         , Left e)          -> Left  e
+                (Left (Right ls) , Left (Right rs)) -> Left  $ Right [ mergeL l r | l <- ls, r <- rs ]
+                (Left (Left ls)  , Left (Left rs))  -> Left  $ Left (ls ++ rs)
+                (Left e          , Left (Right _))  -> Left  e
+                (Left (Right _)  , Left e        )  -> Left  e
+        (fopts,bopts) = splitLists d opts
+        (fi,fres) = compute f fopts us
+        (bi,bres) = compute b bopts us
+    compute (Or d f g) opts us
+      = case () of
+            () | null fonly
+               , null gonly
+               , Left (Right fms) <- fr
+               , Left (Right gms) <- gr -> (hi, Left $ Right $ fms ++ gms)
+
+            () | Left (Left fss) <- fr
+               , Left (Left gss) <- gr -> (hi, Left (Left (fss ++ gss)))
+
+            () | null gonly, Left _ <- gr -> (fi,fr)
+            () | null fonly, Left _ <- fr -> (gi,gr)
+
+            () | null gonly, Right _ <- fr -> (fi,fr)
+            () | null fonly, Right _ <- gr -> (gi,gr)
+
+            () | Left (Left fss) <- fr -> (hi, Left (Left ( filter (not . null) (gonly : map (filter (not . (`elem` gopts))) fss))))
+            () | Left (Left gss) <- gr -> (hi, Left (Left ( filter (not . null) (fonly : map (filter (not . (`elem` fopts))) gss))))
+            () -> (hi, Left (Left [fonly,gonly]))
+
+     where
+        (fopts,gopts) = splitLists d opts
+        (fonly,gonly) = (filterPresent *** filterPresent) $ removeIntersection fopts gopts
+        filterPresent = filter (`IntMap.member` valmap)
+        (fi,fr) = compute f fopts us
+        (gi,gr) = compute g gopts us
+        hi = hull fi gi
+    compute Empty _ _ = error "CommandLine: empty evaluation"
+
+{-
+-- | Phase 2.  Evaluate the given set of named and unnamed values and return
+-- the computed result or else indicate bad usage.
+--
+-- 'parseInvocation' may be used on the results of 'System.Environment.getArgs'
+-- to obtain suitable input for this function.
+runArgsOlder :: ([(String,String)], [String]) -> ArgsOlder a -> Either UsageError a
+runArgsOlder (os,us) (ArgsOlder alts)
+    | not (null rs) = Right $ head rs
+    | not (null ls) = Left $ chooseError ls
+    | otherwise     = Right $ error $ show (length alts,ls)
+ where
+    recs = map (tryCompute os us) alts
+    rs = rights recs
+    ls = lefts recs
+-}
+
+{-
+chooseError :: [UsageError] -> UsageError
+chooseError ls = case span ((==2) . rankError) $ sortOn rankError ls of
+                        ([e],_)    -> e
+                        (e:es,_)
+                           | overlap -> em
+                           | otherwise -> -- trace ("ms="++show ms) $
+                              case findPartition ms of
+                                Just (xs@(_:_:_)) -> ChooseOne $ map return xs
+                                _                 -> em
+                            where
+                                em:ems = sortBy (comparing (maximum . map length . missingWhat)) (e:es)
+                                ms = concatMap missingWhat (em:ems)
+                                mi = foldr1 intersect ms
+                                overlap = any null $ map (\\ mi) ms
+                        (_,e:es)   -> case takeWhile ((>2) . rankError) (e:es) of
+                                        [f]  -> f
+                                        f:fs -> -- trace ("ws="++show (w:ws))
+                                                case u of
+                                                    [_] -> f
+                                                    _   -> Incompatible u
+                                                 where u = foldr1 union $ w : takeWhile ((==wlen) . length) ws
+                                                       w:ws = map misunderstoodWhat (f:fs)
+                                                       wlen = length w
+                                        []   -> e
+-}
+
+
+{-
+-- Given a collection of sets, return a list of unique reprasentative members.
+findPartition :: Eq x => [[x]] -> Maybe [x]
+findPartition yss =
+    case sortBy (comparing length) yss of
+        []:_ -> Nothing
+        zss | not (null ds) -> Nothing
+            | otherwise     -> _findPartition ps es xss3
+          where
+            (pss,xss0) = span isSingle zss
+            isSingle [x] = True
+            isSingle _   = False
+            ps = foldr union [] pss
+            xss1 = map (partition (`elem` ps)) xss0
+            (xss2,bs) = partition (null . fst) xss1
+            (cs,ds) = partition (null . drop 1 . fst) bs
+            es = foldr union [] $ map snd cs
+            xss3 = map snd xss2
+
+
+_findPartition :: Eq x => [x] -> [x] -> [[x]] -> Maybe [x]
+_findPartition ps qs [] = Just ps
+_findPartition ps qs (xs:xss)
+    | null cs = Nothing
+    | otherwise = listToMaybe ss
+ where
+    cs = filter (not . flip elem qs) xs
+    ts = init $ zipWith (\as (b:bs) -> (b,as++bs)) (inits cs) (tails cs)
+    ss = mapMaybe (\(t,tqs) -> _findPartition (t:ps) (tqs++qs) (filter (not . elem t) xss)) ts
+-}
diff --git a/lib/Numeric/Interval.hs b/lib/Numeric/Interval.hs
new file mode 100644
index 0000000..df4bc33
--- /dev/null
+++ b/lib/Numeric/Interval.hs
@@ -0,0 +1,754 @@
+{-# 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
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  DeriveDataTypeable
+-- Version : intervals-0.4.2 (minus distributive instance)
+--
+-- Interval arithmetic
+--
+-----------------------------------------------------------------------------
+
+module Numeric.Interval
+  ( Interval(..)
+  , (...)
+  , whole
+  , empty
+  , null
+  , singleton
+  , elem
+  , notElem
+  , inf
+  , sup
+  , singular
+  , width
+  , midpoint
+  , intersection
+  , hull
+  , bisection
+  , magnitude
+  , mignitude
+  , contains
+  , isSubsetOf
+  , certainly, (<!), (<=!), (==!), (>=!), (>!)
+  , possibly, (<?), (<=?), (==?), (>=?), (>?)
+  , clamp
+  , idouble
+  , ifloat
+  ) where
+
+import Control.Applicative hiding (empty)
+import Data.Data
+#ifdef VERSION_distributive
+import Data.Distributive
+#endif
+import Data.Foldable hiding (minimum, maximum, elem, notElem, null)
+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 (>>=) #-}
+
+#ifdef VERSION_distributive
+instance Distributive Interval where
+  distribute f = fmap inf f ... fmap sup f
+  {-# INLINE distribute #-}
+#endif
+
+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)
diff --git a/lib/Numeric/Interval/Bounded.hs b/lib/Numeric/Interval/Bounded.hs
new file mode 100644
index 0000000..2dd4d7b
--- /dev/null
+++ b/lib/Numeric/Interval/Bounded.hs
@@ -0,0 +1,9 @@
+module Numeric.Interval.Bounded where
+
+import Numeric.Interval
+
+whole' :: Bounded a => Interval a
+whole' = ( minBound ... maxBound )
+
+empty' :: Bounded a => Interval a
+empty' = ( maxBound ... minBound )
diff --git a/lib/SuperOrd.hs b/lib/SuperOrd.hs
new file mode 100644
index 0000000..258a823
--- /dev/null
+++ b/lib/SuperOrd.hs
@@ -0,0 +1,23 @@
+module SuperOrd where
+
+data SuperOrd a
+        = NegativeInfinity
+        | SuperOrd { superApprox        :: !a
+                   , superCompareApprox :: !Ordering
+                   }
+        | PositiveInfinity
+  deriving (Eq, Ord, Show)
+
+instance Bounded (SuperOrd a) where
+    minBound = NegativeInfinity
+    maxBound = PositiveInfinity
+
+exactly :: a -> SuperOrd a
+exactly a = SuperOrd a EQ
+
+lessThan :: a -> SuperOrd a
+lessThan a = SuperOrd a LT
+
+greaterThan :: a -> SuperOrd a
+greaterThan a = SuperOrd a GT
+