diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -1,3 +1,11 @@
+2.1.0
+-----
+
+* introduce operations for Allen's interval algebra (#18, thanks to marcosh)
+* make `recip` precise when 0 is not an interior point (#21)
+* add `instance Storable` for `Interval` (#25)
+* add `instance Floating` for `Interval` (#26)
+
 2.0.0
 -----
 * change internal representation of `Interval` and `IntegerInterval` to
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -1,7 +1,8 @@
 data-interval
 =============
 
-[![Build Status](https://travis-ci.org/msakai/data-interval.svg?branch=master)](https://travis-ci.org/msakai/data-interval)
+[![Build Status (Travis-CI)](https://travis-ci.org/msakai/data-interval.svg?branch=master)](https://travis-ci.org/msakai/data-interval)
+[![Build Status (GitHub Actions)](https://github.com/msakai/data-interval/actions/workflows/build.yaml/badge.svg)](https://github.com/msakai/data-interval/actions/workflows/build.yaml)
 [![Hackage](https://img.shields.io/hackage/v/data-interval.svg)](https://hackage.haskell.org/package/data-interval)
 [![Hackage Deps](https://img.shields.io/hackage-deps/v/data-interval.svg)](https://packdeps.haskellers.com/feed?needle=data-interval)
 [![Coverage Status](https://coveralls.io/repos/msakai/data-interval/badge.svg)](https://coveralls.io/r/msakai/data-interval)
diff --git a/data-interval.cabal b/data-interval.cabal
--- a/data-interval.cabal
+++ b/data-interval.cabal
@@ -1,5 +1,5 @@
 Name:		data-interval
-Version:	2.0.1
+Version:	2.1.0
 License:	BSD3
 License-File:	COPYING
 Author:		Masahiro Sakai (masahiro.sakai@gmail.com)
@@ -24,24 +24,32 @@
    GHC ==8.0.2
    GHC ==8.2.2
    GHC ==8.4.4
-   GHC ==8.6.4
+   GHC ==8.6.5
+   GHC ==8.8.4
+   GHC ==8.10.2
 
 source-repository head
   type:     git
   location: git://github.com/msakai/data-interval.git
 
+flag lattices
+  description: Derive lattice instances
+  default: True
+
 Library
   Hs-source-dirs: src
   Build-Depends:
        base >=4 && <5
      , containers
-     , lattices >=1.2.1.1 && <2.1
      , deepseq
      , hashable >=1.1.2.5 && <1.4
      , extended-reals >=0.2 && <1.0
   if impl(ghc <8.0)
     Build-depends:
       semigroups
+  if flag(lattices)
+    build-depends:
+     lattices >=1.2.1.1 && <2.1
   Default-Language: Haskell2010
   Other-Extensions:
      CPP
@@ -56,6 +64,7 @@
      Data.Interval
      Data.IntervalMap.Lazy
      Data.IntervalMap.Strict
+     Data.IntervalRelation
      Data.IntervalSet
      Data.IntegerInterval
   Other-Modules:
@@ -67,7 +76,13 @@
   Type:              exitcode-stdio-1.0
   HS-Source-Dirs:    test
   Main-is:           TestSuite.hs
-  Other-Modules:     TestInterval, TestIntervalMap, TestIntervalSet, TestIntegerInterval
+  Other-Modules:
+     TestInterval
+     TestIntervalMap
+     TestIntervalRelation
+     TestIntervalSet
+     TestIntegerInterval
+     TestInstances
   Build-depends:
        base >=4 && <5
      , ChasingBottoms
@@ -83,6 +98,9 @@
      , tasty-th
      , HUnit
      , QuickCheck >=2.5 && <3
+  if impl(ghc >=8.0)
+    Build-depends:
+      quickcheck-classes-base
   if impl(ghc <8.0)
     Build-depends:
       semigroups
diff --git a/src/Data/IntegerInterval.hs b/src/Data/IntegerInterval.hs
--- a/src/Data/IntegerInterval.hs
+++ b/src/Data/IntegerInterval.hs
@@ -44,6 +44,7 @@
   , notMember
   , isSubsetOf
   , isProperSubsetOf
+  , isConnected
   , lowerBound
   , upperBound
   , lowerBound'
@@ -77,18 +78,24 @@
   , fromInterval
   , fromIntervalOver
   , fromIntervalUnder
+
+  -- * Intervals relation
+  , relate
   ) where
 
+#if MIN_VERSION_lattices
 import Algebra.Lattice
+#endif
 import Control.Exception (assert)
 import Control.Monad hiding (join)
 import Data.ExtendedReal
-import Data.List hiding (null)
+import Data.List (foldl')
 import Data.Maybe
 import Prelude hiding (null)
 import Data.IntegerInterval.Internal
-import Data.Interval (Boundary(..))
-import qualified Data.Interval as Interval
+import Data.Interval.Internal (Boundary(..))
+import qualified Data.Interval.Internal as Interval
+import Data.IntervalRelation
 
 infix 5 <..<=
 infix 5 <=..<
@@ -128,6 +135,7 @@
     ub@(Finite _) -> (ub, Closed)
     ub@_ -> (ub, Open)
 
+#if MIN_VERSION_lattices
 #if MIN_VERSION_lattices(2,0,0)
 
 instance Lattice IntegerInterval where
@@ -159,6 +167,7 @@
 instance BoundedLattice IntegerInterval
 
 #endif
+#endif
 
 instance Show IntegerInterval where
   showsPrec _ x | null x = showString "empty"
@@ -274,6 +283,23 @@
 isProperSubsetOf :: IntegerInterval -> IntegerInterval -> Bool
 isProperSubsetOf i1 i2 = i1 /= i2 && i1 `isSubsetOf` i2
 
+-- | Does the union of two range form a set which is the intersection between the integers and a connected real interval?
+isConnected :: IntegerInterval -> IntegerInterval -> Bool
+isConnected x y = null x || null y || x ==? y || lb1nearUb2 || ub1nearLb2
+  where
+    lb1 = lowerBound x
+    lb2 = lowerBound y
+    ub1 = upperBound x
+    ub2 = upperBound y
+
+    lb1nearUb2 = case (lb1, ub2) of
+      (Finite lb1Int, Finite ub2Int) -> lb1Int == ub2Int + 1
+      _                              -> False
+
+    ub1nearLb2 = case (ub1, lb2) of
+      (Finite ub1Int, Finite lb2Int) -> ub1Int + 1 == lb2Int
+      _                              -> False
+
 -- | Width of a interval. Width of an unbounded interval is @undefined@.
 width :: IntegerInterval -> Integer
 width x
@@ -462,7 +488,9 @@
 
 -- | Convert the interval to 'Interval.Interval' data type.
 toInterval :: Real r => IntegerInterval -> Interval.Interval r
-toInterval x = fmap fromInteger (lowerBound x) Interval.<=..<= fmap fromInteger (upperBound x)
+toInterval x = Interval.interval
+  (fmap fromInteger (lowerBound x), Closed)
+  (fmap fromInteger (upperBound x), Closed)
 
 -- | Conversion from 'Interval.Interval' data type.
 fromInterval :: Interval.Interval Integer -> IntegerInterval
@@ -481,18 +509,46 @@
 fromIntervalOver :: RealFrac r => Interval.Interval r -> IntegerInterval
 fromIntervalOver i = fmap floor lb <=..<= fmap ceiling ub
   where
-    lb = Interval.lowerBound i
-    ub = Interval.upperBound i
+    (lb, _) = Interval.lowerBound' i
+    (ub, _) = Interval.upperBound' i
 
 -- | Given a 'Interval.Interval' @I@ over R, compute the largest 'IntegerInterval' @J@ such that @J ⊆ I@.
 fromIntervalUnder :: RealFrac r => Interval.Interval r -> IntegerInterval
-fromIntervalUnder i = fmap f lb <=..<= fmap g ub
+fromIntervalUnder i = lb <=..<= ub
   where
-    lb = Interval.lowerBound i
-    ub = Interval.upperBound i
-    f x = if fromIntegral y `Interval.member` i then y else y+1
-      where
-        y = ceiling x
-    g x = if fromIntegral y `Interval.member` i then y else y-1
-      where
-        y = floor x
+    lb = case Interval.lowerBound' i of
+      (Finite x, Open)
+        | fromInteger (ceiling x) == x
+        -> Finite (ceiling x + 1)
+      (x, _) -> fmap ceiling x
+    ub = case Interval.upperBound' i of
+      (Finite x, Open)
+        | fromInteger (floor x) == x
+        -> Finite (floor x - 1)
+      (x, _) -> fmap floor x
+
+-- | Computes how two intervals are related according to the @`Data.IntervalRelation.Relation`@ classification
+relate :: IntegerInterval -> IntegerInterval -> Relation
+relate i1 i2 =
+  case (i1 `isSubsetOf` i2, i2 `isSubsetOf` i1) of
+    -- 'i1' ad 'i2' are equal
+    (True , True ) -> Equal
+    -- 'i1' is strictly contained in `i2`
+    (True , False) | lowerBound i1 == lowerBound i2 -> Starts
+                   | upperBound i1 == upperBound i2 -> Finishes
+                   | otherwise                      -> During
+    -- 'i2' is strictly contained in `i1`
+    (False, True ) | lowerBound i1 == lowerBound i2 -> StartedBy
+                   | upperBound i1 == upperBound i2 -> FinishedBy
+                   | otherwise                      -> Contains
+    -- neither `i1` nor `i2` is contained in the other
+    (False, False) -> case ( null (i1 `intersection` i2)
+                           , lowerBound i1 <= lowerBound i2
+                           , i1 `isConnected` i2
+                           ) of
+      (True , True , True ) -> JustBefore
+      (True , True , False) -> Before
+      (True , False, True ) -> JustAfter
+      (True , False, False) -> After
+      (False, True , _    ) -> Overlaps
+      (False, False, _    ) -> OverlappedBy
diff --git a/src/Data/Interval.hs b/src/Data/Interval.hs
--- a/src/Data/Interval.hs
+++ b/src/Data/Interval.hs
@@ -7,7 +7,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Interval
--- Copyright   :  (c) Masahiro Sakai 2011-2013
+-- Copyright   :  (c) Masahiro Sakai 2011-2013, Andrew Lelechenko 2020
 -- License     :  BSD-style
 --
 -- Maintainer  :  masahiro.sakai@gmail.com
@@ -45,6 +45,7 @@
   -- * Query
   , null
   , isSingleton
+  , extractSingleton
   , member
   , notMember
   , isSubsetOf
@@ -77,14 +78,20 @@
   -- * Operations
   , pickup
   , simplestRationalWithin
+
+  -- * Intervals relation
+  , relate
   ) where
 
+#if MIN_VERSION_lattices
 import Algebra.Lattice
+#endif
 import Control.Exception (assert)
 import Control.Monad hiding (join)
 import Data.ExtendedReal
 import Data.Interval.Internal
-import Data.List hiding (null)
+import Data.IntervalRelation
+import Data.List (foldl', maximumBy, minimumBy)
 import Data.Maybe
 import Data.Monoid
 import Data.Ratio
@@ -113,6 +120,7 @@
 infix 4 >??
 infix 4 /=??
 
+#if MIN_VERSION_lattices
 #if MIN_VERSION_lattices(2,0,0)
 
 instance (Ord r) => Lattice (Interval r) where
@@ -144,6 +152,7 @@
 instance (Ord r) => BoundedLattice (Interval r)
 
 #endif
+#endif
 
 instance (Ord r, Show r) => Show (Interval r) where
   showsPrec _ x | null x = showString "empty"
@@ -320,10 +329,17 @@
 --
 -- @since 2.0.0
 isSingleton :: Ord r => Interval r -> Bool
-isSingleton i = case (lowerBound' i, upperBound' i) of
-  ((Finite l, Closed), (Finite u, Closed)) -> l==u
-  _ -> False
+isSingleton = isJust . extractSingleton
 
+-- | If the interval is a single point, return this point.
+--
+-- @since 2.1.0
+extractSingleton :: Ord r => Interval r -> Maybe r
+extractSingleton i = case (lowerBound' i, upperBound' i) of
+  ((Finite l, Closed), (Finite u, Closed))
+    | l == u -> Just l
+  _ -> Nothing
+
 -- | Is the element in the interval?
 member :: Ord r => r -> Interval r -> Bool
 member x i = condLB && condUB
@@ -426,13 +442,22 @@
         Finite lb = lowerBound j
         lb_floor  = floor lb
 
--- | @mapMonotonic f i@ is the image of @i@ under @f@, where @f@ must be a strict monotone function.
+-- | @mapMonotonic f i@ is the image of @i@ under @f@, where @f@ must be a strict monotone function,
+-- preserving negative and positive infinities.
 mapMonotonic :: (Ord a, Ord b) => (a -> b) -> Interval a -> Interval b
 mapMonotonic f i = interval (fmap f lb, in1) (fmap f ub, in2)
   where
     (lb, in1) = lowerBound' i
     (ub, in2) = upperBound' i
 
+mapAntiMonotonic :: (Ord a, Ord b) => (a -> b) -> Interval a -> Interval b
+mapAntiMonotonic f i
+  | null i = empty
+  | otherwise = interval (fmap f ub, in2) (fmap f lb, in1)
+  where
+    (lb, in1) = lowerBound' i
+    (ub, in2) = upperBound' i
+
 -- | For all @x@ in @X@, @y@ in @Y@. @x '<' y@?
 (<!) :: Ord r => Interval r -> Interval r -> Bool
 a <! b =
@@ -598,6 +623,8 @@
     lb = lowerBound' x
     ub = upperBound' x
 
+-- | When results of 'abs' or 'signum' do not form a connected interval,
+-- a convex hull is returned instead.
 instance (Num r, Ord r) => Num (Interval r) where
   a + b
     | null a || null b = empty
@@ -639,17 +666,142 @@
       ub3 = maximumBy cmpUB xs
       lb3 = minimumBy cmpLB xs
 
+-- | 'recip' returns 'whole' when 0 is an interior point.
+-- Otherwise @recip (recip xs)@ equals to @xs@ without 0.
 instance forall r. (Real r, Fractional r) => Fractional (Interval r) where
   fromRational r = singleton (fromRational r)
   recip a
     | null a = empty
-    | 0 `member` a = whole -- should be error?
+    | a == 0 = empty
+    | 0 `member` a && 0 /= lowerBound a && 0 /= upperBound a = whole
     | otherwise = interval lb3 ub3
     where
       ub3 = maximumBy cmpUB xs
       lb3 = minimumBy cmpLB xs
       xs = [recipLB (lowerBound' a), recipUB (upperBound' a)]
 
+-- | When results of 'tan' or '**' do not form a connected interval,
+-- a convex hull is returned instead.
+instance (RealFrac r, Floating r) => Floating (Interval r) where
+  pi = singleton pi
+
+  exp = intersection (0 <..< PosInf) . mapMonotonic exp
+  log a = interval (logB (lowerBound' b)) (logB (upperBound' b))
+    where
+      b = intersection (0 <..< PosInf) a
+
+  sqrt = mapMonotonic sqrt . intersection (0 <=..< PosInf)
+
+  a ** b = hulls (posBase : negBasePosPower : negBaseNegPower : zeroPower ++ zeroBase)
+    where
+      posBase = exp (log a * b)
+      zeroPower = [ 1 | 0 `member` b, not (null a) ]
+      zeroBase  = [ 0 | 0 `member` a, not (null (b `intersection` (0 <..< PosInf))) ]
+      negBasePosPower = positiveIntegralPowersOfNegativeValues
+        (a `intersection` (NegInf <..< 0))
+        (b `intersection` (0 <..< PosInf))
+      negBaseNegPower = positiveIntegralPowersOfNegativeValues
+        (recip  (a `intersection` (NegInf <..< 0)))
+        (negate (b `intersection` (NegInf <..< 0)))
+
+  cos a = case lowerBound' a of
+    (NegInf, _) -> -1 <=..<= 1
+    (PosInf, _) -> empty
+    (Finite lb, in1) -> case upperBound' a of
+      (NegInf, _) -> empty
+      (PosInf, _) -> -1 <=..<= 1
+      (Finite ub, in2)
+        | ub - lb > 2 * pi                                             -> -1 <=..<= 1
+        | clb == -1 && ub - lb == 2 * pi && in1 == Open && in2 == Open -> -1 <..<= 1
+        | clb ==  1 && ub - lb == 2 * pi && in1 == Open && in2 == Open -> -1 <=..< 1
+        | ub - lb == 2 * pi                                            -> -1 <=..<= 1
+
+        | lbNorth, ubNorth, clb >= cub -> interval (cub, in2) (clb, in1)
+        | lbNorth, ubNorth -> -1 <=..<= 1
+        | lbNorth -> interval (-1, Closed) $ case clb `compare` cub of
+          LT -> (cub, in2)
+          EQ -> (cub, in1 `max` in2)
+          GT -> (clb, in1)
+        | ubNorth -> (`interval` (1, Closed)) $ case clb `compare` cub of
+          LT -> (clb, in1)
+          EQ -> (clb, in1 `max` in2)
+          GT -> (cub, in2)
+        | clb > cub -> -1 <=..<= 1
+        | otherwise -> interval (clb, in1) (cub, in2)
+        where
+          mod2pi x = let y = x / (2 * pi) in y - fromInteger (floor y)
+          -- is lower bound in the northern half-plane [0,pi)?
+          lbNorth = (mod2pi lb, in1) < (1 / 2, Closed)
+          -- is upper bound in the northern half-plane [0,pi)?
+          ubNorth = (mod2pi ub, in2) < (1 / 2, Closed)
+          clb = Finite (cos lb)
+          cub = Finite (cos ub)
+
+  acos = mapAntiMonotonic acos . intersection (-1 <=..<= 1)
+
+  sin a = cos (pi / 2 - a)
+  asin = mapMonotonic asin . intersection (-1 <=..<= 1)
+
+  tan a = case lowerBound' a of
+    (NegInf, _) -> whole
+    (PosInf, _) -> empty
+    (Finite lb, in1) -> case upperBound' a of
+      (NegInf, _) -> empty
+      (PosInf, _) -> whole
+      (Finite ub, in2)
+        | ub - lb > pi -> whole
+        -- the next case corresponds to (tan lb, +inf) + (-inf, tan ub)
+        -- with tan lb == tan ub, but a convex hull is returned instead
+        | ub - lb == pi && in1 == Open && in2 == Open && modpi lb /= 1/2 -> whole
+        | ub - lb == pi -> whole
+        | tan lb <= tan ub -> interval (Finite $ tan lb, in1) (Finite $ tan ub, in2)
+        -- the next case corresponds to (tan lb, +inf) + (-inf, tan ub),
+        -- but a convex hull is returned instead
+        | otherwise -> whole
+        where
+          modpi x = let y = x / pi in y - fromInteger (floor y)
+
+  atan = intersection (Finite (-pi / 2) <=..<= Finite (pi / 2)) . mapMonotonic atan
+
+  sinh  = mapMonotonic sinh
+  asinh = mapMonotonic asinh
+
+  cosh  = mapMonotonic cosh . abs
+  acosh = mapMonotonic acosh . intersection (1 <=..< PosInf)
+
+  tanh  = intersection (-1 <..< 1) . mapMonotonic tanh
+  atanh a = interval (atanhB (lowerBound' b)) (atanhB (upperBound' b))
+    where
+      b = intersection (-1 <..< 1) a
+
+positiveIntegralPowersOfNegativeValues
+  :: RealFrac r => Interval r -> Interval r -> Interval r
+positiveIntegralPowersOfNegativeValues a b
+  | null a || null b         = empty
+  | Just ub <- mub, lb > ub  = empty
+  | Just ub <- mub, lb == ub = a ^ lb
+  -- cases below connects two intervals (a ^ k, 0) + (0, a ^ k'))
+  -- into a single convex hull
+  | lowerBound a >= -1       = hull (a ^ lb) (a ^ (lb + 1))
+  | Just ub <- mub           = hull (a ^ ub) (a ^ (ub - 1))
+  | Nothing <- mub           = whole
+  where
+    -- Similar to Data.IntegerInterval.fromIntervalUnder
+    lb :: Integer
+    lb = case lowerBound' b of
+      (Finite x, Open)
+        | fromInteger (ceiling x) == x
+        -> ceiling x + 1
+      (Finite x, _) -> ceiling x
+      _ -> 0 -- PosInf is not expected, because b is not null
+    mub :: Maybe Integer
+    mub = case upperBound' b of
+      (Finite x, Open)
+        | fromInteger (floor x) == x
+        -> Just $ floor x - 1
+      (Finite x, _) -> Just $ floor x
+      _ -> Nothing -- NegInf is not expected, because b is not null
+
 cmpUB, cmpLB :: Ord r => (Extended r, Boundary) -> (Extended r, Boundary) -> Ordering
 cmpUB (x1,in1) (x2,in2) = compare x1 x2 `mappend` compare in1 in2
 cmpLB (x1,in1) (x2,in2) = compare x1 x2 `mappend` compare in2 in1
@@ -684,3 +836,42 @@
 recipUB :: (Fractional r, Ord r) => (Extended r, Boundary) -> (Extended r, Boundary)
 recipUB (0, _) = (NegInf, Open)
 recipUB (x1, in1) = (recip x1, in1)
+
+logB :: (Floating r, Ord r) => (Extended r, Boundary) -> (Extended r, Boundary)
+logB (NegInf, in1) = (Finite $ log (log 0), in1)
+logB (Finite 0, _) = (NegInf, Open)
+logB (Finite x1, in1) = (Finite $ log x1, in1)
+logB (PosInf, in1) = (PosInf, in1)
+
+atanhB :: (Floating r, Ord r) => (Extended r, Boundary) -> (Extended r, Boundary)
+atanhB (NegInf, in1) = (Finite $ atanh (-1/0), in1)
+atanhB (Finite (-1), _) = (NegInf, Open)
+atanhB (Finite 1, _) = (PosInf, Open)
+atanhB (Finite x1, in1) = (Finite $ atanh x1, in1)
+atanhB (PosInf, in1) = (Finite $ atanh (1/0), in1)
+
+-- | Computes how two intervals are related according to the @`Data.IntervalRelation.Relation`@ classification
+relate :: Ord r => Interval r -> Interval r -> Relation
+relate i1 i2 =
+  case (i1 `isSubsetOf` i2, i2 `isSubsetOf` i1) of
+    -- 'i1' ad 'i2' are equal
+    (True , True ) -> Equal
+    -- 'i1' is strictly contained in `i2`
+    (True , False) | lowerBound i1 == lowerBound i2 -> Starts
+                   | upperBound i1 == upperBound i2 -> Finishes
+                   | otherwise                      -> During
+    -- 'i2' is strictly contained in `i1`
+    (False, True ) | lowerBound i1 == lowerBound i2 -> StartedBy
+                   | upperBound i1 == upperBound i2 -> FinishedBy
+                   | otherwise                      -> Contains
+    -- neither `i1` nor `i2` is contained in the other
+    (False, False) -> case ( null (i1 `intersection` i2)
+                           , upperBound' i1 <= upperBound' i2
+                           , i1 `isConnected` i2
+                           ) of
+      (True , True , True ) -> JustBefore
+      (True , True , False) -> Before
+      (True , False, True ) -> JustAfter
+      (True , False, False) -> After
+      (False, True , _    ) -> Overlaps
+      (False, False, _    ) -> OverlappedBy
diff --git a/src/Data/Interval/Internal.hs b/src/Data/Interval/Internal.hs
--- a/src/Data/Interval/Internal.hs
+++ b/src/Data/Interval/Internal.hs
@@ -1,5 +1,5 @@
 {-# OPTIONS_GHC -Wall #-}
-{-# LANGUAGE CPP, DeriveDataTypeable, DeriveGeneric, LambdaCase #-}
+{-# LANGUAGE CPP, DeriveDataTypeable, DeriveGeneric, LambdaCase, ScopedTypeVariables #-}
 {-# LANGUAGE Safe #-}
 #if __GLASGOW_HASKELL__ >= 708
 {-# LANGUAGE RoleAnnotations #-}
@@ -14,10 +14,17 @@
   , empty
   ) where
 
+#if __GLASGOW_HASKELL__ < 710
+import Control.Applicative hiding (empty)
+#endif
 import Control.DeepSeq
 import Data.Data
 import Data.ExtendedReal
 import Data.Hashable
+import Data.Int
+import Foreign.Marshal.Array
+import Foreign.Ptr
+import Foreign.Storable
 import GHC.Generics (Generic)
 
 -- | Boundary of an interval may be
@@ -42,11 +49,62 @@
   | LessOrEqual !r
   | GreaterThan !r
   | GreaterOrEqual !r
+  -- For constructors below
+  -- the first argument is strictly less than the second one
   | BothClosed !r !r
   | LeftOpen !r !r
   | RightOpen !r !r
   | BothOpen !r !r
   deriving (Eq, Typeable)
+
+peekInterval :: (Applicative m, Monad m, Ord r) => m Int8 -> m r -> m r -> m (Interval r)
+peekInterval tagM x y = do
+  tag <- tagM
+  case tag of
+    0 -> pure Whole
+    1 -> pure Empty
+    2 -> Point           <$> x
+    3 -> LessThan        <$> x
+    4 -> LessOrEqual     <$> x
+    5 -> GreaterThan     <$> x
+    6 -> GreaterOrEqual  <$> x
+    7 -> wrap BothClosed <$> x <*> y
+    8 -> wrap LeftOpen   <$> x <*> y
+    9 -> wrap RightOpen  <$> x <*> y
+    _ -> wrap BothOpen   <$> x <*> y
+
+-- | Enforce the internal invariant
+-- of 'BothClosed' / 'LeftOpen' / 'RightOpen' / 'BothOpen'.
+wrap :: Ord r => (r -> r -> Interval r) -> r -> r -> Interval r
+wrap f x y
+  | x < y = f x y
+  | otherwise = Empty
+
+pokeInterval :: Applicative m => (Int8 -> m ()) -> (r -> m ()) -> (r -> m ()) -> Interval r -> m ()
+pokeInterval tag actX actY = \case
+  Whole            -> tag (0 :: Int8)
+  Empty            -> tag (1 :: Int8)
+  Point          x -> tag (2 :: Int8) *> actX x
+  LessThan       x -> tag (3 :: Int8) *> actX x
+  LessOrEqual    x -> tag (4 :: Int8) *> actX x
+  GreaterThan    x -> tag (5 :: Int8) *> actX x
+  GreaterOrEqual x -> tag (6 :: Int8) *> actX x
+  BothClosed   x y -> tag (7 :: Int8) *> actX x *> actY y
+  LeftOpen     x y -> tag (8 :: Int8) *> actX x *> actY y
+  RightOpen    x y -> tag (9 :: Int8) *> actX x *> actY y
+  BothOpen     x y -> tag (10 :: Int8) *> actX x *> actY y
+
+instance (Storable r, Ord r) => Storable (Interval r) where
+  sizeOf _ = 3 * sizeOf (undefined :: r)
+  alignment _ = alignment (undefined :: r)
+  peek ptr = peekInterval
+    (peek $ castPtr ptr)
+    (peek $ castPtr ptr `advancePtr` 1)
+    (peek $ castPtr ptr `advancePtr` 2)
+  poke ptr = pokeInterval
+    (poke $ castPtr ptr)
+    (poke $ castPtr ptr `advancePtr` 1)
+    (poke $ castPtr ptr `advancePtr` 2)
 
 -- | Lower endpoint (/i.e./ greatest lower bound) of the interval,
 -- together with 'Boundary' information.
diff --git a/src/Data/IntervalRelation.hs b/src/Data/IntervalRelation.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/IntervalRelation.hs
@@ -0,0 +1,59 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# LANGUAGE CPP, DeriveDataTypeable, DeriveGeneric #-}
+{-# LANGUAGE Safe #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.IntervalRelation
+-- Copyright   :  (c) Masahiro Sakai 2016
+-- License     :  BSD-style
+--
+-- Maintainer  :  masahiro.sakai@gmail.com
+-- Stability   :  provisional
+-- Portability :  non-portable (CPP, DeriveDataTypeable, DeriveGeneric)
+--
+-- Interval relations and their algebra.
+--
+-----------------------------------------------------------------------------
+module Data.IntervalRelation
+  ( Relation(..)
+  , invert
+  )
+  where
+
+import Data.Data
+import GHC.Generics (Generic)
+
+-- | describes how two intervals @x@ and @y@ can be related.
+-- See [Allen's interval algebra](https://en.wikipedia.org/wiki/Allen%27s_interval_algebra)
+data Relation
+  = Before
+  | JustBefore
+  | Overlaps
+  | Starts
+  | During
+  | Finishes
+  | Equal
+  | FinishedBy
+  | Contains
+  | StartedBy
+  | OverlappedBy
+  | JustAfter
+  | After
+  deriving (Eq, Ord, Enum, Bounded, Show, Read, Generic, Data, Typeable)
+
+-- | inverts a relation, such that @'invert' ('Data.Interval.relate' x y) = 'Data.Interval.relate' y x@
+invert :: Relation -> Relation
+invert relation = case relation of
+  Before       -> After
+  JustBefore   -> JustAfter
+  Overlaps     -> OverlappedBy
+  Starts       -> StartedBy
+  During       -> Contains
+  Finishes     -> FinishedBy
+  Equal        -> Equal
+  FinishedBy   -> Finishes
+  Contains     -> During
+  StartedBy    -> Starts
+  OverlappedBy -> Overlaps
+  JustAfter    -> JustBefore
+  After        -> Before
diff --git a/src/Data/IntervalSet.hs b/src/Data/IntervalSet.hs
--- a/src/Data/IntervalSet.hs
+++ b/src/Data/IntervalSet.hs
@@ -62,7 +62,9 @@
   where
 
 import Prelude hiding (null, span)
+#if MIN_VERSION_lattices
 import Algebra.Lattice
+#endif
 import Control.DeepSeq
 import Data.Data
 import Data.ExtendedReal
@@ -131,6 +133,7 @@
 instance Hashable r => Hashable (IntervalSet r) where
   hashWithSalt s (IntervalSet m) = hashWithSalt s (Map.toList m)
 
+#if MIN_VERSION_lattices
 #if MIN_VERSION_lattices(2,0,0)
 
 instance (Ord r) => Lattice (IntervalSet r) where
@@ -162,6 +165,7 @@
 instance (Ord r) => BoundedLattice (IntervalSet r)
 
 #endif
+#endif
 
 instance Ord r => Monoid (IntervalSet r) where
   mempty = empty
@@ -216,9 +220,10 @@
       ]
     return y
 
+-- | @recip (recip xs) == delete 0 xs@
 instance forall r. (Real r, Fractional r) => Fractional (IntervalSet r) where
   fromRational r = singleton (fromRational r)
-  recip = lift1 recip
+  recip xs = lift1 recip (delete (Interval.singleton 0) xs)
 
 #if __GLASGOW_HASKELL__ >= 708
 instance Ord r => GHCExts.IsList (IntervalSet r) where
diff --git a/test/TestInstances.hs b/test/TestInstances.hs
new file mode 100644
--- /dev/null
+++ b/test/TestInstances.hs
@@ -0,0 +1,47 @@
+module TestInstances where
+
+import Control.Monad
+
+import Test.Tasty.QuickCheck
+
+import Data.Interval
+import Data.IntervalRelation
+
+instance Arbitrary Boundary where
+  arbitrary = arbitraryBoundedEnum
+
+instance Arbitrary r => Arbitrary (Extended r) where
+  arbitrary = frequency
+    [ (1, return NegInf)
+    , (1, return PosInf)
+    , (3, liftM Finite arbitrary)
+    ]
+  shrink NegInf = []
+  shrink (Finite x) = NegInf : PosInf : map Finite (shrink x)
+  shrink PosInf = []
+
+instance (Arbitrary r, Ord r) => Arbitrary (Interval r) where
+  arbitrary = do
+    x <- arbitrary
+    y <- arbitrary
+    frequency
+      [ (1, return $ interval x y)
+      , (3, return $ interval (min x y) (max x y))
+      ]
+  shrink a
+    | isSingleton a = case lowerBound a of
+      Finite x -> map singleton $ shrink x
+      _ -> []
+    | otherwise = mkPoint lb ++ mkPoint ub ++ map (lb `interval`) (shrink ub) ++ map (`interval` ub) (shrink lb)
+    where
+      lb = lowerBound' a
+      ub = upperBound' a
+
+      mkPoint (Finite x, _) = [singleton x]
+      mkPoint _ = []
+
+intervals :: Gen (Interval Rational)
+intervals = arbitrary
+
+instance Arbitrary Relation where
+  arbitrary = arbitraryBoundedEnum
diff --git a/test/TestIntegerInterval.hs b/test/TestIntegerInterval.hs
--- a/test/TestIntegerInterval.hs
+++ b/test/TestIntegerInterval.hs
@@ -1,7 +1,9 @@
-{-# LANGUAGE TemplateHaskell, ScopedTypeVariables #-}
+{-# LANGUAGE CPP, TemplateHaskell, ScopedTypeVariables #-}
 module TestIntegerInterval (integerIntervalTestGroup) where
 
+#if MIN_VERSION_lattices
 import qualified Algebra.Lattice as L
+#endif
 import Control.DeepSeq
 import Control.Monad
 import Data.Generics.Schemes
@@ -24,6 +26,7 @@
 import qualified Data.IntegerInterval as IntegerInterval
 import Data.Interval (Interval)
 import qualified Data.Interval as Interval
+import Data.IntervalRelation
 
 {--------------------------------------------------------------------
   empty
@@ -228,6 +231,19 @@
 case_isProperSubsetOf =
   (0 <=..<= 1) `IntegerInterval.isProperSubsetOf` (0 <=..<= 2) @?= True
 
+{-- -----------------------------------------------------------------
+  isConnected
+----------------------------------------------------------------- --}
+
+prop_isConnected_reflexive =
+  forAll integerIntervals $ \a ->
+    a `IntegerInterval.isConnected` a
+
+prop_isConnected_symmetric =
+  forAll integerIntervals $ \a ->
+    forAll integerIntervals $ \b ->
+      (a `IntegerInterval.isConnected` b) == (b `IntegerInterval.isConnected` a)
+
 {--------------------------------------------------------------------
   simplestIntegerWithin
 --------------------------------------------------------------------}
@@ -284,6 +300,43 @@
     IntegerInterval.pickup (IntegerInterval.singleton x) == Just x
 
 {--------------------------------------------------------------------
+  relate
+--------------------------------------------------------------------}
+
+prop_relate_equals =
+  forAll integerIntervals $ \a ->
+    IntegerInterval.relate a a == Equal
+
+prop_relate_empty_contained_in_non_empty =
+  forAll (integerIntervals `suchThat` (not . IntegerInterval.null)) $ \a ->
+    IntegerInterval.relate a IntegerInterval.empty == Contains
+
+prop_relate_detects_before =
+  forAll (nonEmptyIntegerIntervalPairs (\_ ub1 lb2 _ -> ub1 < lb2 - 1)) $ \(a, b) ->
+    IntegerInterval.relate a b == Before
+
+prop_relate_detects_just_before =
+  forAll (arbitrary `suchThat` \(b1, b2, i) -> b1 <= Finite i &&  Finite (i + 1) <= b2) $
+      \(b1, b2, i) ->
+        IntegerInterval.relate (b1 <=..<= Finite i) (Finite (i + 1) <=..<= b2) == JustBefore
+
+prop_relate_two_intervals_overlap =
+  forAll (nonEmptyIntegerIntervalPairs (\lb1 ub1 lb2 ub2 -> lb1 < lb2 && lb2 < ub1 && ub1 < ub2)) $ \(a, b) ->
+    IntegerInterval.relate a b == Overlaps
+
+prop_relate_interval_starts_another =
+  forAll (nonEmptyIntegerIntervalPairs (\lb1 ub1 lb2 ub2 -> lb1 == lb2 && ub1 < ub2)) $ \(a, b) ->
+    IntegerInterval.relate a b == Starts
+
+prop_relate_interval_finishes_another =
+  forAll (nonEmptyIntegerIntervalPairs (\lb1 ub1 lb2 ub2 -> lb1 > lb2 && ub1 == ub2)) $ \(a, b) ->
+    IntegerInterval.relate a b == Finishes
+
+prop_relate_interval_contains_another =
+  forAll (nonEmptyIntegerIntervalPairs (\lb1 ub1 lb2 ub2 -> lb1 < lb2 && ub1 > ub2)) $ \(a, b) ->
+    IntegerInterval.relate a b == Contains
+
+{--------------------------------------------------------------------
   Comparison
 --------------------------------------------------------------------}
 
@@ -671,6 +724,8 @@
   Lattice
 --------------------------------------------------------------------}
 
+#if MIN_VERSION_lattices
+
 prop_Lattice_Leq_welldefined =
   forAll integerIntervals $ \a b ->
     a `L.meetLeq` b == a `L.joinLeq` b
@@ -683,6 +738,14 @@
   forAll integerIntervals $ \a ->
     L.bottom `L.joinLeq` a
 
+#else
+
+prop_Lattice_Leq_welldefined = True
+prop_top                     = True
+prop_bottom                  = True
+
+#endif
+
 {--------------------------------------------------------------------
   Read
 --------------------------------------------------------------------}
@@ -789,6 +852,24 @@
 
 integerIntervals :: Gen IntegerInterval
 integerIntervals = arbitrary
+
+nonEmptyIntegerIntervalPairs
+  :: ( Extended Integer
+    -> Extended Integer
+    -> Extended Integer
+    -> Extended Integer
+    -> Bool)
+  -> Gen (IntegerInterval, IntegerInterval)
+nonEmptyIntegerIntervalPairs boundariesComparer = ap (fmap (,) integerIntervals) integerIntervals `suchThat`
+  (\(i1, i2) ->
+    (not . IntegerInterval.null $ i1) &&
+    (not . IntegerInterval.null $ i2) &&
+    boundariesComparer
+      (IntegerInterval.lowerBound i1)
+      (IntegerInterval.upperBound i1)
+      (IntegerInterval.lowerBound i2)
+      (IntegerInterval.upperBound i2)
+  )
 
 intervals :: Gen (Interval.Interval Rational)
 intervals = arbitrary
diff --git a/test/TestInterval.hs b/test/TestInterval.hs
--- a/test/TestInterval.hs
+++ b/test/TestInterval.hs
@@ -1,11 +1,15 @@
-{-# LANGUAGE TemplateHaskell, ScopedTypeVariables #-}
+{-# LANGUAGE CPP, TemplateHaskell, RankNTypes, ScopedTypeVariables #-}
 module TestInterval (intervalTestGroup) where
 
+#if MIN_VERSION_lattices
 import qualified Algebra.Lattice as L
+#endif
 import Control.DeepSeq
+import Control.Exception
 import Control.Monad
 import Data.Generics.Schemes
 import Data.Hashable
+import Data.Int
 import Data.Maybe
 import Data.Ratio
 import Data.Typeable
@@ -13,7 +17,11 @@
 import Test.Tasty
 import Test.Tasty.QuickCheck
 import Test.Tasty.HUnit
+import Test.Tasty.Options
 import Test.Tasty.TH
+#ifdef MIN_VERSION_quickcheck_classes_base
+import Test.QuickCheck.Classes.Base
+#endif
 
 import Data.Interval
   ( Interval, Extended (..), (<=..<=), (<=..<), (<..<=), (<..<)
@@ -22,7 +30,10 @@
   , (<??), (<=??), (==??), (>=??), (>??), (/=??)
   )
 import qualified Data.Interval as Interval
+import Data.IntervalRelation
 
+import TestInstances
+
 {--------------------------------------------------------------------
   empty
 --------------------------------------------------------------------}
@@ -198,7 +209,11 @@
   forAll intervals $ \a ->
     Interval.isSubsetOf a a
 
-prop_isSubsetOf_trans =
+test_isSubsetOf_trans :: [TestTree]
+test_isSubsetOf_trans =
+  (: []) $
+  adjustOption (\(QuickCheckMaxRatio r) -> QuickCheckMaxRatio (r * 10)) $
+  testProperty "isSubsetOf trans" $
   forAll intervals $ \a ->
   forAll intervals $ \b ->
   forAll intervals $ \c ->
@@ -319,6 +334,73 @@
   Interval.mapMonotonic (+1) (0 <=..< 10) @?= ((1 <=..<11) :: Interval Rational)
 
 {--------------------------------------------------------------------
+  relate
+--------------------------------------------------------------------}
+
+prop_relate_equals =
+  forAll intervals $ \a ->
+    Interval.relate a a == Equal
+
+prop_relate_empty_contained_in_non_empty =
+  forAll (intervals `suchThat` (not . Interval.null)) $ \a ->
+    Interval.relate a Interval.empty == Contains
+
+prop_relate_detects_before =
+  forAll (nonEmptyIntervalPairs (\_ (ub1, _) (lb2, _) _ -> ub1 < lb2)) $ \(a, b) ->
+    Interval.relate a b == Before
+
+prop_relate_open_intervals_with_common_boundary_are_before =
+  forAll (arbitrary `suchThat` \(b1, b2, i) -> fst b1 < i && i < fst b2) $
+      \(b1 :: (Extended Rational, Interval.Boundary), b2, i :: Extended Rational) ->
+        Interval.relate (Interval.interval b1 (i, Interval.Open)) (Interval.interval (i, Interval.Open) b2) == Before
+
+prop_relate_right_closed_interval_just_before =
+  forAll (arbitrary `suchThat` \(b1, b2, i) -> fst b1 < i && i < fst b2) $
+      \(b1 :: (Extended Rational, Interval.Boundary), b2, i :: Extended Rational) ->
+        Interval.relate (Interval.interval b1 (i, Interval.Closed)) (Interval.interval (i, Interval.Open) b2) == JustBefore
+
+prop_relate_right_open_interval_just_before =
+  forAll (arbitrary `suchThat` \(b1, b2, i) -> fst b1 < i && i < fst b2) $
+      \(b1 :: (Extended Rational, Interval.Boundary), b2, i :: Extended Rational) ->
+        Interval.relate (Interval.interval b1 (i, Interval.Open)) (Interval.interval (i, Interval.Closed) b2) == JustBefore
+
+prop_relate_two_intervals_overlap =
+  forAll (nonEmptyIntervalPairs (\(lb1, _) (ub1, _) (lb2, _) (ub2, _) -> lb1 < lb2 && lb2 < ub1 && ub1 < ub2)) $ \(a, b) ->
+    Interval.relate a b == Overlaps
+
+prop_relate_interval_starts_another =
+  forAll (nonEmptyIntervalPairs (\lb1 (ub1, _) lb2 (ub2, _) -> lb1 == lb2 && ub1 < ub2)) $ \(a, b) ->
+    Interval.relate a b == Starts
+
+prop_relate_interval_finishes_another =
+  forAll (nonEmptyIntervalPairs (\(lb1, _) ub1 (lb2, _) ub2 -> lb1 > lb2 && ub1 == ub2)) $ \(a, b) ->
+    Interval.relate a b == Finishes
+
+prop_relate_interval_contains_another =
+  forAll (nonEmptyIntervalPairs (\(lb1, _) (ub1, _) (lb2, _) (ub2, _) -> lb1 < lb2 && ub1 > ub2)) $ \(a, b) ->
+    Interval.relate a b == Contains
+
+prop_relate_one_singleton_before_another =
+  forAll (arbitrary `suchThat` uncurry (<)) $ \(r1 :: Rational, r2) ->
+    Interval.relate (Interval.singleton r1) (Interval.singleton r2) == Before
+
+prop_relate_singleton_starts_interval =
+  forAll (arbitrary `suchThat` uncurry (<)) $ \(r1 :: Rational, r2) b ->
+    Interval.relate (Interval.singleton r1) (Interval.interval (Finite r1, Interval.Closed) (Finite r2, b)) == Starts
+
+prop_relate_singleton_just_before_interval =
+  forAll (arbitrary `suchThat` uncurry (<)) $ \(r1 :: Rational, r2) b ->
+    Interval.relate (Interval.singleton r1) (Interval.interval (Finite r1, Interval.Open) (Finite r2, b)) == JustBefore
+
+prop_relate_singleton_finishes_interval =
+  forAll (arbitrary `suchThat` uncurry (<)) $ \(r1 :: Rational, r2) b ->
+    Interval.relate (Interval.singleton r2) (Interval.interval (Finite r1, b) (Finite r2, Interval.Closed)) == Finishes
+
+prop_relate_singleton_just_after_interval =
+  forAll (arbitrary `suchThat` uncurry (<)) $ \(r1 :: Rational, r2) b ->
+    Interval.relate (Interval.singleton r2) (Interval.interval (Finite r1, b) (Finite r2, Interval.Open)) == JustAfter
+
+{--------------------------------------------------------------------
   Comparison
 --------------------------------------------------------------------}
 
@@ -749,14 +831,332 @@
     i1 = -10 <=..< 0
     i2 = NegInf <..<= (-1/10)
 
-prop_recip_zero =
+case_recip_test4 = recip i1 @?= i2
+  where
+    i1, i2 :: Interval Rational
+    i1 = 0 <=..<= 10
+    i2 = (1/10) <=..< PosInf
+
+case_recip_test5 = recip i1 @?= i2
+  where
+    i1, i2 :: Interval Rational
+    i1 = -10 <=..<= 0
+    i2 = NegInf <..<= (-1/10)
+
+case_recip_test6 = recip i1 @?= i2
+  where
+    i1, i2 :: Interval Rational
+    i1 = 0 <=..<= 0
+    i2 = Interval.empty
+
+prop_recip =
   forAll intervals $ \a ->
-    0 `Interval.member` a ==> recip a == Interval.whole
+    if 0 `isInteriorPoint` a
+    then recip a === Interval.whole
+    else recip (recip a) === without0 a
 
+isInteriorPoint :: (Ord a, Show a) => a -> Interval a -> Bool
+isInteriorPoint x xs
+  = x `Interval.member` xs
+  && Finite x /= Interval.lowerBound xs
+  && Finite x /= Interval.upperBound xs
+
+without0 :: (Ord a, Num a) => Interval a -> Interval a
+without0 xs = case Interval.lowerBound' xs of
+  (0, Interval.Closed) ->
+    Interval.interval (0, Interval.Open) (Interval.upperBound' xs)
+  _ -> case Interval.upperBound' xs of
+    (0, Interval.Closed) ->
+      Interval.interval (Interval.lowerBound' xs) (0, Interval.Open)
+    _ -> xs
+
 {--------------------------------------------------------------------
+  Floating
+--------------------------------------------------------------------}
+
+prop_exp_singleton = floatingSingleton exp
+
+prop_exp_mid_point = floatingMidPoint exp
+
+case_exp_whole = exp Interval.whole @?= 0 <..< PosInf
+
+case_exp_empty = exp Interval.empty @?= Interval.empty
+
+prop_log_singleton a = a > 0 ==>
+  floatingSingleton log a
+
+prop_log_mid_point = floatingMidPoint log . Interval.intersection (0 <..< PosInf)
+
+case_log_whole = log Interval.whole   @?= Interval.whole
+case_log_half1 = log (0 <=..< PosInf) @?= Interval.whole
+case_log_half2 = log (0 <..< PosInf)  @?= Interval.whole
+case_log_zero  = log (0 :: Interval Double) @?= Interval.empty
+
+case_log_empty = log Interval.empty @?= Interval.empty
+
+prop_log_exp a = log (exp a) =~= a
+
+prop_exp_log a = exp (log a) =~= a `Interval.intersection` (0 <..< PosInf)
+
+-------------------------------------------------------------------------------
+
+prop_sqrt_singleton = floatingSingleton sqrt
+
+prop_sqrt_mid_point = floatingMidPoint sqrt . Interval.intersection (0 <=..< PosInf)
+
+case_sqrt_whole = sqrt Interval.whole @?= 0 <=..< PosInf
+
+case_sqrt_empty = sqrt Interval.empty @?= Interval.empty
+
+prop_sqr_sqrt a = sqrt a * sqrt a =~= a `Interval.intersection` (0 <=..< PosInf)
+
+prop_sqrt_sqr a = sqrt (a * a) =~= abs a
+
+-------------------------------------------------------------------------------
+
+prop_pow_singleton_Double_Double a' b' =
+  not (isInfinite c || isNaN c) ==>
+    Interval.singleton a ** Interval.singleton b =~= Interval.singleton c
+  where
+    a = min 5 $ max (-5) a'
+    b = min 5 $ max (-5) b'
+    c = a ** b
+
+prop_pow_singleton_Double_Integer 0 b'
+  | b' < 0 = discard
+prop_pow_singleton_Double_Integer a' b' =
+  Interval.singleton a ** Interval.singleton b =~= Interval.singleton (a ** b)
+  where
+    a = min 5 $ max (-5) a'
+    b = min 5 $ max (-5) $ fromInteger b'
+
+prop_pow_singleton_Integer_Double a' b =
+  not (isInfinite c || isNaN c) ==>
+    Interval.singleton a ** Interval.singleton b =~= Interval.singleton (a ** b)
+  where
+    a = fromInteger a'
+    c = a ** b
+
+prop_pow_mid_point a' b' = case (Interval.pickup a, Interval.pickup b) of
+  (Nothing, _) -> discard
+  (_, Nothing) -> discard
+  (Just x, Just y) -> let z = x ** y :: Double in not (isInfinite z || isNaN z) ==>
+    ioProperty $ do
+      x <- try (evaluate (a ** b))
+      return $ case x of
+        Left LossOfPrecision -> discard
+        Right c -> distance z c < Finite (1e-10 * (1 `max` abs z))
+  where
+    -- for larger intervals the loss of precision becomes exponentially huge
+    a = Interval.mapMonotonic (min 5 . max (-5)) a'
+    b = Interval.mapMonotonic (min 5 . max (-5)) b'
+
+prop_pow_empty_1 :: Interval Double -> Bool
+prop_pow_empty_1 x = Interval.null (Interval.empty ** x)
+
+prop_pow_empty_2 :: Interval Double -> Bool
+prop_pow_empty_2 x = Interval.null (x ** Interval.empty)
+
+-------------------------------------------------------------------------------
+
+prop_sin_singleton a =
+  distance (sin a :: Double) (sin (Interval.singleton a)) <= 1e-10
+
+prop_sin_mid_point a
+  | Interval.isSingleton a = discard
+  | otherwise = floatingMidPoint sin a
+
+case_sin_whole = sin Interval.whole @?= -1 <=..<= 1
+
+case_sin_empty = sin Interval.empty @?= Interval.empty
+
+prop_asin_singleton a = floatingSingleton asin (if abs a < 1 then a else recip a)
+
+prop_asin_mid_point = floatingMidPoint asin . Interval.intersection (-1 <=..<= 1)
+
+case_asin_whole = asin Interval.whole @?= Finite (-pi / 2) <=..<= Finite (pi / 2)
+
+case_asin_empty = asin Interval.empty @?= Interval.empty
+
+prop_sin_asin a = sin (asin a) =~= a `Interval.intersection` (-1 <=..<= 1)
+
+-------------------------------------------------------------------------------
+
+prop_cos_singleton a =
+  distance (cos a :: Double) (cos (Interval.singleton a)) <= 1e-10
+
+prop_cos_mid_point a
+  | Interval.isSingleton a = discard
+  | otherwise = floatingMidPoint cos a
+
+case_cos_whole = cos Interval.whole @?= -1 <=..<= 1
+
+case_cos_empty = cos Interval.empty @?= Interval.empty
+
+prop_acos_singleton a = floatingSingleton acos (if abs a < 1 then a else recip a)
+
+prop_acos_mid_point = floatingMidPoint acos . Interval.intersection (-1 <=..<= 1)
+
+case_acos_whole = acos Interval.whole @?= 0 <=..<= Finite pi
+
+case_acos_empty = acos Interval.empty @?= Interval.empty
+
+prop_cos_acos a = cos (acos a) =~= a `Interval.intersection` (-1 <=..<= 1)
+
+-------------------------------------------------------------------------------
+
+prop_tan_singleton a =
+  distance (tan a :: Double) (tan (Interval.singleton a)) <= 1e-10
+
+prop_tan_mid_point a = case Interval.pickup a of
+  Nothing -> discard
+  Just x -> let z = tan x :: Double in not (isInfinite z || isNaN z) ==>
+    ioProperty $ do
+      x <- try (evaluate (tan a))
+      return $ case x of
+        Left LossOfPrecision -> discard
+        Right c -> distance z c < Finite (1e-10 * (1 `max` abs z))
+
+case_tan_whole = tan Interval.whole @?= Interval.whole
+
+case_tan_empty = tan Interval.empty @?= Interval.empty
+
+prop_atan_singleton = floatingSingleton atan
+
+prop_atan_mid_point = floatingMidPoint atan
+
+case_atan_whole = atan Interval.whole @?= Finite (-pi / 2) <=..<= Finite (pi / 2)
+
+case_atan_empty = atan Interval.empty @?= Interval.empty
+
+prop_tan_atan a = case (Interval.lowerBound a, Interval.upperBound a) of
+  (Finite{}, Finite{}) -> tan (atan a) =~= a
+  _ -> discard
+
+-------------------------------------------------------------------------------
+
+prop_sinh_singleton = floatingSingleton sinh
+
+prop_sinh_mid_point = floatingMidPoint sinh
+
+case_sinh_whole = sinh Interval.whole @?= Interval.whole
+
+case_sinh_empty = sinh Interval.empty @?= Interval.empty
+
+prop_asinh_singleton = floatingSingleton asinh
+
+prop_asinh_mid_point = floatingMidPoint asinh
+
+case_asinh_whole = asinh Interval.whole @?= Interval.whole
+
+case_asinh_empty = asinh Interval.empty @?= Interval.empty
+
+prop_asinh_sinh a' = asinh (sinh a) =~= a
+  where
+    -- for larger intervals the loss of precision becomes exponentially huge
+    a = Interval.mapMonotonic (min 5 . max (-5)) a'
+
+prop_sinh_asinh a = sinh (asinh a) =~= a
+
+-------------------------------------------------------------------------------
+
+prop_cosh_singleton = floatingSingleton cosh
+
+prop_cosh_mid_point = floatingMidPoint cosh
+
+case_cosh_whole = cosh Interval.whole @?= 1 <=..< PosInf
+
+case_cosh_empty = cosh Interval.empty @?= Interval.empty
+
+prop_acosh_singleton = floatingSingleton acosh
+
+prop_acosh_mid_point = floatingMidPoint acosh . Interval.intersection (1 <=..< PosInf)
+
+case_acosh_whole = acosh Interval.whole @?= 0 <=..< PosInf
+
+case_acosh_empty = acosh Interval.empty @?= Interval.empty
+
+prop_acosh_cosh a' = acosh (cosh a) =~= abs a
+  where
+    -- for larger intervals the loss of precision becomes exponentially huge
+    a = Interval.mapMonotonic (min 5 . max (-5)) a'
+
+prop_cosh_acosh a = cosh (acosh a) =~= a `Interval.intersection` (1 <=..< PosInf)
+
+-------------------------------------------------------------------------------
+
+prop_tanh_singleton a = abs a <= 10 ==>
+  floatingSingleton tanh a
+
+prop_tanh_mid_point = floatingMidPoint tanh . Interval.intersection (-5 <=..<= 5)
+
+case_tanh_whole = tanh Interval.whole @?= -1 <..< 1
+
+case_tanh_empty = tanh Interval.empty @?= Interval.empty
+
+prop_atanh_singleton 1    = atanh 1 === Interval.empty
+prop_atanh_singleton (-1) = atanh (-1) === Interval.empty
+prop_atanh_singleton a    = floatingSingleton atanh (if abs a < 1 then a else recip a)
+
+prop_atanh_mid_point = floatingMidPoint atanh . Interval.intersection (-1 <..< 1)
+
+case_atanh_whole = atanh Interval.whole @?= Interval.whole
+
+case_atanh_empty = atanh Interval.empty @?= Interval.empty
+
+prop_atanh_tanh a' = atanh (tanh a) =~= a
+  where
+    -- for larger intervals the loss of precision becomes exponentially huge
+    a = Interval.mapMonotonic (min 5 . max (-5)) a'
+
+prop_tanh_atanh = uncurry (=~=) . tanhAtanh
+
+case_tanh_atanh_1 = uncurry (@?=) $ tanhAtanh (-1 <=..<= 1)
+case_tanh_atanh_2 = uncurry (@?=) $ tanhAtanh (-1 <=..< 1)
+case_tanh_atanh_3 = uncurry (@?=) $ tanhAtanh (-1 <..<= 1)
+case_tanh_atanh_4 = uncurry (@?=) $ tanhAtanh (-1 <..< 1)
+
+tanhAtanh :: Interval Double -> (Interval Double, Interval Double)
+tanhAtanh a = (tanh (atanh a), a `Interval.intersection` (-1 <..< 1))
+
+-------------------------------------------------------------------------------
+
+floatingSingleton :: (forall a. Floating a => a -> a) -> Double -> Property
+floatingSingleton f a = Interval.singleton (f a) === f (Interval.singleton a)
+
+distance :: (Ord r, Num r) => r -> Interval r -> Extended r
+distance x xs
+  | Interval.member x xs = 0
+  | otherwise
+  = abs (Finite x - Interval.lowerBound xs) `min`
+    abs (Finite x - Interval.upperBound xs)
+
+floatingMidPoint :: (forall a. Floating a => a -> a) -> Interval Double -> Property
+floatingMidPoint f a = case Interval.pickup a of
+  Nothing -> discard
+  Just x  -> property $ f x `Interval.member` f a
+
+infix 4 =~=
+(=~=) :: Interval Double -> Interval Double -> Property
+a =~= b
+  | eqPair (Interval.lowerBound' a) (Interval.lowerBound' b)
+  , eqPair (Interval.upperBound' a) (Interval.upperBound' b)
+  = property True
+  | otherwise
+  = a === b
+  where
+    eqPair (x, a) (y, b) = eqExt x y && a == b
+
+    eqExt (Finite x) (Finite y) =
+      abs (x - y) < 1e-10 * (1 `max` abs x `max` abs y)
+    eqExt x y = x == y
+
+{--------------------------------------------------------------------
   Lattice
 --------------------------------------------------------------------}
 
+#if MIN_VERSION_lattices
+
 prop_Lattice_Leq_welldefined =
   forAll intervals $ \a b ->
     a `L.meetLeq` b == a `L.joinLeq` b
@@ -769,6 +1169,14 @@
   forAll intervals $ \a ->
     L.bottom `L.joinLeq` a
 
+#else
+
+prop_Lattice_Leq_welldefined = True
+prop_top                     = True
+prop_bottom                  = True
+
+#endif
+
 {--------------------------------------------------------------------
   Read
 --------------------------------------------------------------------}
@@ -810,28 +1218,44 @@
       | otherwise = x
 
 {--------------------------------------------------------------------
-  Generators
+  Storable
 --------------------------------------------------------------------}
 
-instance Arbitrary Interval.Boundary where
-  arbitrary = arbitraryBoundedEnum
+#ifdef MIN_VERSION_quickcheck_classes_base
+test_Storable_Int8 = map (uncurry testProperty) $ lawsProperties $
+  storableLaws (Proxy :: Proxy (Interval Int8))
+test_Storable_Int = map (uncurry testProperty) $ lawsProperties $
+  storableLaws (Proxy :: Proxy (Interval Int))
+#else
+test_Storable_Int8 = []
+test_Storable_Int = []
+#endif
 
-instance Arbitrary r => Arbitrary (Extended r) where
-  arbitrary =
-    oneof
-    [ return NegInf
-    , return PosInf
-    , liftM Finite arbitrary
-    ]
+{--------------------------------------------------------------------
+  Generators
+--------------------------------------------------------------------}
 
-instance (Arbitrary r, Ord r) => Arbitrary (Interval r) where
-  arbitrary = do
-    lb <- arbitrary
-    ub <- arbitrary
-    return $ Interval.interval lb ub
+nonEmptyIntervalPairs
+  :: ( (Extended Rational, Interval.Boundary)
+    -> (Extended Rational, Interval.Boundary)
+    -> (Extended Rational, Interval.Boundary)
+    -> (Extended Rational, Interval.Boundary)
+    -> Bool)
+  -> Gen (Interval Rational, Interval Rational)
+nonEmptyIntervalPairs boundariesComparer = ap (fmap (,) intervals) intervals `suchThat`
+  (\(i1, i2) ->
+    (not . Interval.null $ i1) &&
+    (not . Interval.null $ i2) &&
+    boundariesComparer
+      (Interval.lowerBound' i1)
+      (Interval.upperBound' i1)
+      (Interval.lowerBound' i2)
+      (Interval.upperBound' i2)
+  )
 
-intervals :: Gen (Interval Rational)
-intervals = arbitrary
+{--------------------------------------------------------------------
+  Test intervals
+--------------------------------------------------------------------}
 
 pos :: Interval Rational
 pos = 0 <..< PosInf
diff --git a/test/TestIntervalRelation.hs b/test/TestIntervalRelation.hs
new file mode 100644
--- /dev/null
+++ b/test/TestIntervalRelation.hs
@@ -0,0 +1,27 @@
+{-# LANGUAGE ScopedTypeVariables, TemplateHaskell #-}
+module TestIntervalRelation (intervalRelationTestGroup) where
+
+import Test.Tasty.QuickCheck
+import Test.Tasty.TH
+
+import Data.Interval
+import Data.IntervalRelation
+
+import TestInstances
+
+{--------------------------------------------------------------------
+  invert
+--------------------------------------------------------------------}
+
+prop_invert_is_involution a =
+  invert (invert a) == a
+
+prop_invert_inverts_relation =
+  forAll intervals $ \a ->
+  forAll intervals $ \b ->
+    relate a b == invert (relate b a)
+
+------------------------------------------------------------------------
+-- Test harness
+
+intervalRelationTestGroup = $(testGroupGenerator)
diff --git a/test/TestIntervalSet.hs b/test/TestIntervalSet.hs
--- a/test/TestIntervalSet.hs
+++ b/test/TestIntervalSet.hs
@@ -1,7 +1,9 @@
-{-# LANGUAGE TemplateHaskell, ScopedTypeVariables #-}
+{-# LANGUAGE CPP, TemplateHaskell, ScopedTypeVariables #-}
 module TestIntervalSet (intervalSetTestGroup) where
 
+#if MIN_VERSION_lattices
 import qualified Algebra.Lattice as L
+#endif
 import Control.Applicative ((<$>))
 import Control.DeepSeq
 import Control.Monad
@@ -295,6 +297,8 @@
   Lattice
 --------------------------------------------------------------------}
 
+#if MIN_VERSION_lattices
+
 prop_Lattice_Leq_welldefined =
   forAll arbitrary $ \(a :: IntervalSet Rational) (b :: IntervalSet Rational) ->
     a `L.meetLeq` b == a `L.joinLeq` b
@@ -307,6 +311,14 @@
   forAll arbitrary $ \(a :: IntervalSet Rational) ->
     L.bottom `L.joinLeq` a
 
+#else
+
+prop_Lattice_Leq_welldefined = True
+prop_top                     = True
+prop_bottom                  = True
+
+#endif
+
 {--------------------------------------------------------------------
   Show / Read
 --------------------------------------------------------------------}
@@ -443,9 +455,8 @@
         d = fromIntegral (denominator r)
     in fromRational n / fromRational d == (fromRational (r::Rational) :: IntervalSet Rational)
 
-prop_recip_zero =
-  forAll arbitrary $ \(a :: IntervalSet Rational) ->
-    0 `IntervalSet.member` a ==> recip a == IntervalSet.whole
+prop_recip (a :: IntervalSet Rational) =
+  recip (recip a) === IntervalSet.delete (Interval.singleton 0) a
 
 {- ------------------------------------------------------------------
   Data
diff --git a/test/TestSuite.hs b/test/TestSuite.hs
--- a/test/TestSuite.hs
+++ b/test/TestSuite.hs
@@ -2,6 +2,7 @@
 
 import TestInterval
 import TestIntervalMap
+import TestIntervalRelation
 import TestIntervalSet
 import TestIntegerInterval
 import Test.Tasty
@@ -10,6 +11,7 @@
 main = defaultMain $ testGroup "data-interval test suite"
   [ intervalTestGroup
   , intervalMapTestGroup
+  , intervalRelationTestGroup
   , intervalSetTestGroup
   , integerIntervalTestGroup
   ]