    {-# LANGUAGE Rank2Types, TypeFamilies #-}
    -----------------------------------------------------------------------------
    -- |
    -- Module      :  Numeric.Interval
    -- Copyright   :  (c) Edward Kmett 2010
    -- License     :  BSD3
    -- Maintainer  :  ekmett@gmail.com
    -- Stability   :  experimental
    -- Portability :  GHC only
    --
    -- Interval arithmetic
    --
    -----------------------------------------------------------------------------

    module Numeric.Interval
        ( Interval(..)
        , (...)
        , whole
        , empty
        , null
        , singleton
        , elem
        , notElem
        , inf
        , sup
        , singular
        , width
        , midpoint
        , intersection
        , hull
        , bisection
        , magnitude
        , mignitude
        , contains
        , isSubsetOf
        , certainly, (<!), (<=!), (==!), (>=!), (>!)
        , possibly, (<?), (<=?), (==?), (>=?), (>?)
        , idouble
        , ifloat
        ) where

    import Prelude hiding (null, elem, notElem)
    import Numeric.Extras
    import Data.Function (on)

    data Interval a = I !a !a

    infix 3 ...

    negInfinity :: Fractional a => a
    negInfinity = (-1)/0
    {-# INLINE negInfinity #-}

    posInfinity :: Fractional a => a
    posInfinity = 1/0
    {-# INLINE posInfinity #-}

    nan :: Fractional a => a
    nan = 0/0

    -- | The rule of thumb is you should only use this to construct using values
    -- that you took out of the interval. Otherwise, use I, to force rounding
    (...) :: a -> a -> Interval a
    a ... b = I a b
    {-# INLINE (...) #-}

    -- | The whole real number line
    whole :: Fractional a => Interval a
    whole = negInfinity ... posInfinity
    {-# INLINE whole #-}

    -- | An empty interval
    empty :: Fractional a => Interval a
    empty = nan ... nan
    {-# INLINE empty #-}

    -- | negation handles NaN properly
    null :: Ord a => Interval a -> Bool
    null x = not (inf x <= sup x)
    {-# INLINE null #-}

    -- | A singleton point
    singleton :: a -> Interval a
    singleton a = a ... a
    {-# INLINE singleton #-}

    -- | The infinumum (lower bound) of an interval
    inf :: Interval a -> a
    inf (I a _) = a
    {-# INLINE inf #-}

    -- | The supremum (upper bound) of an interval
    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 :: Ord a => Interval a -> Bool
    singular x = not (null x) && inf x == sup x
    {-# INLINE singular #-}

    instance Eq a => Eq (Interval a) where
        (==) = (==!)

    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 :: Num a => Interval a -> a
    width (I a b) = b - a
    {-# INLINE width #-}

    -- | magnitude
    magnitude :: (Num a, Ord a) => Interval a -> a
    magnitude x = (max `on` abs) (inf x) (sup x)
    {-# INLINE magnitude #-}

    -- | "mignitude"
    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')
        I a b - I a' b' = (a - b') ... (b - a')
        I a b * I a' b' = minimum [a * a', a * b', b * a', b * b']
                          ...
                          maximum [a * a', a * b', b * a', b * b']
        abs x@(I a b)
            | a >= 0    = x
            | b <= 0    = negate x
            | otherwise = max (- a) b ... b

        signum = increasing signum

        fromInteger i = singleton (fromInteger i)

    -- | Bisect an interval at its midpoint.
    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 :: Fractional a => Interval a -> a
    midpoint x = inf x + (sup x - inf x) / 2
    {-# INLINE midpoint #-}

    elem :: Ord a => a -> Interval a -> Bool
    elem x xs = x >= inf xs && x <= sup xs
    {-# INLINE elem #-}

    notElem :: Ord a => a -> Interval a -> Bool
    notElem x xs = not (elem x xs)
    {-# INLINE notElem #-}

    -- | This means that realToFrac will use the midpoint

    -- | What moron put an Ord instance requirement on Real!
    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)

    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"
        min = minInterval
        max = maxInterval

    -- @'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']
        `I`
        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 `I` ( b / y)
        | b < 0 = negInfinity `I` ( b / y)
        | a < 0 = whole
        | otherwise = (a / y) `I` posInfinity

    -- 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) `I` posInfinity
        | b < 0 = (b / y) `I` posInfinity
        | a < 0 = whole
        | otherwise = negInfinity `I` (a / y)

    divZero :: (Fractional a, Ord a) => Interval a -> Interval a
    divZero x | inf x == 0 && sup x == 0 = x
              | otherwise = whole

    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
        fromRational r  = fromRational r ... fromRational r

    instance RealFloat a => RealFrac (Interval a) where
        properFraction x = (b, x - fromIntegral b)
            where
                b = truncate (midpoint x)
        ceiling x = ceiling (sup x)
        floor x = floor (inf x)
        round x = round (midpoint x)
        truncate x = truncate (midpoint x)

    instance (RealExtras a, Ord a) => Floating (Interval a) where
        pi = singleton pi
        exp = increasing exp
        log (I a b) = (if a > 0 then log a else negInfinity) ... log b
        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)
        sin x
            | null x = empty
            | otherwise = cos (x - pi / 2)
        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
        asin x@(I a b)
            | null x || b < -1 || a > 1 = empty
            | otherwise =
                (if a <= -1 then -halfPi else asin a)
                `I`
                (if b >= 1 then halfPi else asin b)
            where
                halfPi = pi / 2
        acos x@(I a b)
            | null x || b < -1 || a > 1 = empty
            | otherwise =
                (if b >= 1 then 0 else acos b)
                `I`
                (if a < -1 then pi else acos a)
        atan = increasing atan
        sinh = increasing sinh
        cosh x@(I a b)
            | null x = empty
            | b < 0  = decreasing cosh x
            | a >= 0 = increasing cosh x
            | otherwise  = I 0 $ cosh $ if - a > b
                                        then a
                                        else b
        tanh = increasing tanh
        asinh = increasing asinh
        acosh x@(I a b)
            | null x || b < 1 = empty
            | otherwise = I lo $ acosh b
            where lo | a <= 1 = 0
                     | otherwise = acosh a
        atanh x@(I a b)
            | null x || b < -1 || a > 1 = empty
            | otherwise =
                    (if a <= - 1 then negInfinity else atanh a)
                    `I`
                    (if b >= 1 then posInfinity else atanh b)

    -- | lift a monotone increasing function over a given interval
    increasing :: (a -> a) -> Interval a -> Interval a
    increasing f (I a b) = I (f a) (f b)

    -- | lift a monotone increasing function over a given interval
    decreasing :: (a -> a) -> Interval a -> Interval a
    decreasing f (I a b) = I (f b) (f a)

    -- | We have to play some semantic games to make these methods make sense.
    -- Most compute with the midpoint of the interval.
    instance RealExtras 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 :: (Fractional a, Ord a) => Interval a -> Interval a -> Interval a
    intersection x@(I a b) y@(I a' b')
        | x /=! y = empty
        | otherwise = I (max a a') (min b b')
    {-# INLINE intersection #-}

    -- | Calculate the convex hull of two intervals
    hull :: Ord a => Interval a -> Interval a -> Interval a
    hull x@(I a b) y@(I a' b')
        | null x = y
        | null y = x
        | otherwise = I (min a a') (max b b')
    {-# INLINE hull #-}

    instance RealExtras a => RealExtras (Interval a) where
        type C (Interval a) = C a
        fmod x y | null y = empty
                 | otherwise = r -- `intersection` bounds
            where
                n :: Integer
                n = floor (inf x / if inf x < 0 then inf y else sup y)
                r = x - fromIntegral n * y
                -- bounds | inf y >= 0 = y
                --        | otherwise = y `hull` negate y
        expm1 = increasing expm1
        log1p (I a b) = (if a > (-1) then log1p a else negInfinity) `I` log1p b
        hypot x y = hypot a a' `I` hypot b b'
            where
                I a b = abs x
                I a' b' = abs y
        cbrt = increasing cbrt
        erf = increasing erf

    -- | For all @x@ in @X@, @y@ in @Y@. @x '<' y@
    (<!)  :: Ord a => Interval a -> Interval a -> Bool
    x <! y = sup x < inf y
    {-# INLINE (<!) #-}

    -- | For all @x@ in @X@, @y@ in @Y@. @x '<=' y@
    (<=!) :: Ord a => Interval a -> Interval a -> Bool
    x <=! y = sup x <= inf y
    {-# INLINE (<=!) #-}

    -- | For all @x@ in @X@, @y@ in @Y@. @x '==' y@
    (==!) :: Eq a => Interval a -> Interval a -> Bool
    x ==! y = sup x == inf y && inf x == sup y
    {-# INLINE (==!) #-}

    -- | For all @x@ in @X@, @y@ in @Y@. @x '/=' y@
    (/=!) :: Ord a => Interval a -> Interval a -> Bool
    x /=! y = sup x < inf y || inf x > sup y
    {-# INLINE (/=!) #-}

    -- | For all @x@ in @X@, @y@ in @Y@. @x '>' y@
    (>!)  :: Ord a => Interval a -> Interval a -> Bool
    x >! y = inf x > sup y
    {-# INLINE (>!) #-}

    -- | For all @x@ in @X@, @y@ in @Y@. @x '>=' y@
    (>=!) :: Ord a => Interval a -> Interval a -> Bool
    x >=! y = inf x >= sup y
    {-# INLINE (>=!) #-}

    -- | For all @x@ in @X@, @y@ in @Y@. @x `op` y@
    certainly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool
    certainly cmp l r
        | lt && eq && gt = True
        | lt && eq       = l <=! r
        | lt &&       gt = l /=! r
        | lt             = l <! r
        |       eq && gt = l >=! r
        |       eq       = l ==! r
        |             gt = l >! r
        | otherwise      = False
        where
            lt = cmp LT EQ
            eq = cmp EQ EQ
            gt = cmp GT EQ
    {-# INLINE certainly #-}

    contains :: Ord a => Interval a -> Interval a -> Bool
    contains x y = null y
                || (not (null x) && inf x <= inf y && sup y <= sup x)
    {-# INLINE contains #-}

    isSubsetOf :: Ord a => Interval a -> Interval a -> Bool
    isSubsetOf = flip contains

    maxInterval :: Ord a => Interval a -> Interval a -> Interval a
    maxInterval  (I a b) (I a' b') = I (max a a') (max b b')

    minInterval :: Ord a => Interval a -> Interval a -> Interval a
    minInterval  (I a b) (I a' b') = I (min a a') (min b b')

    -- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<' y@?
    (<?) :: Ord a => Interval a -> Interval a -> Bool
    x <? y = inf x < sup y
    {-# INLINE (<?) #-}

    -- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<=' y@?
    (<=?) :: Ord a => Interval a -> Interval a -> Bool
    x <=? y = inf x <= sup y
    {-# INLINE (<=?) #-}

    -- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '==' y@?
    (==?) :: Ord a => Interval a -> Interval a -> Bool
    x ==? y = inf x <= sup y && sup x >= inf y
    {-# INLINE (==?) #-}

    -- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '/=' y@?
    (/=?) :: Eq a => Interval a -> Interval a -> Bool
    x /=? y = inf x /= sup y || sup x /= inf y
    {-# INLINE (/=?) #-}

    -- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>' y@?
    (>?) :: Ord a => Interval a -> Interval a -> Bool
    x >? y = sup x > inf y
    {-# INLINE (>?) #-}

    -- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>=' y@?
    (>=?) :: Ord a => Interval a -> Interval a -> Bool
    x >=? y = sup x >= inf y
    {-# INLINE (>=?) #-}

    -- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x `op` y@?
    possibly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool
    possibly cmp l r
        | lt && eq && gt = True
        | lt && eq       = l <=? r
        | lt &&       gt = l /=? r
        | lt             = l <? r
        |       eq && gt = l >=? r
        |       eq       = l ==? r
        |             gt = l >? r
        | otherwise      = False
        where
            lt = cmp LT EQ
            eq = cmp EQ EQ
            gt = cmp GT EQ
    {-# INLINE possibly #-}

    idouble :: Interval Double -> Interval Double
    idouble = id

    ifloat :: Interval Float -> Interval Float
    ifloat = id

    -- Bugs:
    -- sin 1 :: Interval Double