intervals (empty) → 0.1.0
raw patch · 4 files changed
+538/−0 lines, 4 filesdep +arraydep +basedep +numeric-extrassetup-changed
Dependencies added: array, base, numeric-extras, rounding
Files
- LICENSE +31/−0
- Numeric/Interval.hs +485/−0
- Setup.lhs +3/−0
- intervals.cabal +19/−0
+ LICENSE view
@@ -0,0 +1,31 @@+Copyright (c) 2010, Edward Kmett++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are+met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Edward Kmett nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Numeric/Interval.hs view
@@ -0,0 +1,485 @@+{-# LANGUAGE Rank2Types, TypeFamilies #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Interval+-- Copyright : (c) Edward Kmett 2010+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- Interval arithmetic+--+-----------------------------------------------------------------------------++module Numeric.Interval + ( whole+ , empty+ , null+ , singleton+ , elem+ , notElem+ , inf+ , sup+ , singular+ , width+ , intersection+ , hull+ , bisection+ , magnitude+ , mignitude+ , contains+ , isSubsetOf+ , certainly, (<!), (<=!), (==!), (>=!), (>!)+ , possibly, (<?), (<=?), (==?), (>=?), (>?)+ , idouble + , ifloat + ) where++import Prelude hiding (null, elem, notElem)+import Numeric.Extras+import Numeric.Rounding+import Data.Function (on)++data Interval a = I !(Round Down a) !(Round Up a)++infix 3 ...++negInfinity :: Fractional a => a+negInfinity = (-1)/0 +{-# INLINE negInfinity #-}++posInfinity :: Fractional a => a+posInfinity = 1/0+{-# INLINE posInfinity #-}++nan :: Fractional a => a +nan = 0/0++-- | The rule of thumb is you should only use this to construct using values+-- that you took out of the interval. Otherwise, use I, to force rounding+(...) :: a -> a -> Interval a +a ... b = I (Round a) (Round b)+{-# INLINE (...) #-}++-- | The whole real number line+whole :: Precision a => Interval a +whole = negInfinity ... posInfinity+{-# INLINE whole #-}++-- | An empty interval+empty :: Precision a => Interval a +empty = nan ... nan+{-# INLINE empty #-}++-- | negation handles NaN properly+null :: Ord a => Interval a -> Bool+null x = not (inf x <= sup x)+{-# INLINE null #-}++-- | A singleton point+singleton :: a -> Interval a +singleton a = a ... a+{-# INLINE singleton #-}++-- | The infinumum (lower bound) of an interval+inf :: Interval a -> a+inf (I (Round a) _) = a+{-# INLINE inf #-}++-- | The supremum (upper bound) of an interval+sup :: Interval a -> a+sup (I _ (Round b)) = b+{-# INLINE sup #-}++-- | Is the interval a singleton point? +-- N.B. This is fairly fragile and likely will not hold after+-- even a few operations that only involve singletons+singular :: Ord a => Interval a -> Bool+singular x = not (null x) && inf x == sup x+{-# INLINE singular #-}++instance Precision a => Eq (Interval a) where+ (==) = (==) `on` midpoint++instance Show a => Show (Interval a) where+ showsPrec n (I (Round a) (Round b)) = + showParen (n > 3) $+ showsPrec 3 a .+ showString " ... " . + showsPrec 3 b++-- flip the rounding mode up+u :: Round Down a -> Round Up a+u (Round a) = Round a+{-# INLINE u #-}++-- flip the rounding mode down+d :: Round Up a -> Round Down a+d (Round a) = Round a+{-# INLINE d #-}++-- | Calculate the width of an interval.+-- N.B. the width of an interval is an interval itself due to rounding+width :: Precision a => Interval a -> Interval a+width (I a b) = I (d b - a) (b - u a)+{-# INLINE width #-}++-- | magnitude +magnitude :: Precision a => Interval a -> a +magnitude x = (max `on` abs) (inf x) (sup x)+{-# INLINE magnitude #-}++-- | "mignitude"+mignitude :: Precision a => Interval a -> a +mignitude x = (min `on` abs) (inf x) (sup x)+{-# INLINE mignitude #-}++instance Precision a => Num (Interval a) where+ I a b + I a' b' = I (a + a') (b + b')+ I a b - I a' b' = I (a - d b') (b - u a')+ I a b * I a' b' = minimum [a * a',a * d b',d b * a',d b * d b'] + `I` + maximum [u a * u a',u a * b',b * u a',b * b']+ abs x@(I a b) + | a >= 0 = x + | b <= 0 = negate x+ | otherwise = max (- a) (d b) `I` b++ signum (I a b) = signum a `I` signum b++ fromInteger i = fromInteger i `I` fromInteger i++-- | Bisect an interval at its midpoint.+bisection :: Precision a => Interval a -> (Interval a, Interval a)+bisection (I a b) = (I a (u a + (b - u a) / 2), I (a + (d b - a) / 2) b)+{-# INLINE bisection #-}++-- | Nearest point to the midpoint of the interval.+midpoint :: Precision a => Interval a -> a+midpoint x = inf x + (sup x - inf x) / 2+{-# INLINE midpoint #-}++elem :: Precision a => a -> Interval a -> Bool+elem x xs = x >= inf xs && x <= sup xs+{-# INLINE elem #-}++notElem :: Precision a => a -> Interval a -> Bool+notElem x xs = not (elem x xs)+{-# INLINE notElem #-}++-- | This means that realToFrac will use the midpoint+instance Precision a => Real (Interval a) where+ toRational x = toRational (midpoint x)++-- @'divNonZero' X Y@ assumes @0 `'notElem'` Y@+divNonZero :: Precision a => Interval a -> Interval a -> Interval a +divNonZero (I a b) (I a' b') = + minimum [a / a',a / d b',d b / a',d b / d b'] + `I`+ maximum [u a / u a',u a / b',b / u a',b / b']++-- @'divPositive' X y@ assumes y > 0, and divides @X@ by [0 ... y]+divPositive :: Precision a => Interval a -> a -> Interval a +divPositive x@(I a b) y+ | a == 0 && b == 0 = x+ | b < 0 || isNegativeZero b = negInfinity `I` ( b / up y)+ | a < 0 = whole + -- | isNegativeZero a = whole+ | otherwise = (a / down y) `I` posInfinity++-- divNegative assumes y < 0 and divides the interval @X@ by [y ... 0]+divNegative :: Precision a => Interval a -> a -> Interval a+divNegative x@(I a b) y+ | a == 0 && b == 0 = - x -- flips negative zeros+ | b < 0 || isNegativeZero b = (d b / down y) `I` posInfinity+ | a < 0 = whole+ -- | isNegativeZero a = whole+ | otherwise = negInfinity `I` (u a / up y)++divZero :: Precision a => Interval a -> Interval a+divZero x | inf x == 0 && sup x == 0 = x+ | otherwise = whole++instance Precision a => Fractional (Interval a) where+ -- TODO: check isNegativeZero properly+ x / y+ | 0 `notElem` y = divNonZero x y + | iz && sz = empty -- division by 0+ | iz = divPositive x (inf y)+ | sz = divNegative x (sup y)+ | otherwise = divZero x+ where + iz = inf y == 0+ sz = sup y == 0+ recip (I a b) = on min recip a (d b) `I` on max recip (u a) b+ fromRational r = fromRational r `I` fromRational r++instance Precision a => RealFrac (Interval a) where+ properFraction x = (b, x - fromIntegral b)+ where + b = truncate (midpoint x)+ ceiling x = ceiling (sup x)+ floor x = floor (inf x)+ round x = round (midpoint x)+ truncate x = truncate (midpoint x)++instance Precision a => Floating (Interval a) where+ pi = pi `I` pi + exp = increasing exp+ log (I a b) = (if a > 0 then log a else negInfinity) `I` log b+ cos x + | null x = empty+ | inf (width t) >= inf pi = (-1) ... 1+ | tl >= d pih = - cos (t - pi)+ | th <= u pil = cos (d th) `I` cos (u tl)+ | th <= u pi2l = (-1) `I` cos (u (min (pi2l - d th) tl))+ | otherwise = (-1) ... 1+ where + I pil pih = pi+ pi2@(I pi2l _) = pi * 2+ t@(I tl th) = x `fmod` pi2+ sin x + | null x = empty+ | otherwise = cos (x - pi / 2)+ tan x + | null x = empty+ | inf t' <= -hpil || sup t' >= hpil = whole+ | otherwise = increasing tan x+ where+ t = x `fmod` pi + t' | inf t >= hpil = t - pi+ | otherwise = t+ hpil = inf (pi / 2)+ asin x@(I a b)+ | null x || b < -1 || a > 1 = empty+ | otherwise = + (if a <= - 1 then - d hpis else asin a)+ `I`+ (if b >= 1 then hpis else asin b)+ where+ I _ hpis = pi / 2+ acos x@(I a b)+ | null x || b < -1 || a > 1 = empty+ | otherwise = + (if b >= 1 then 0 else acos (d b))+ `I`+ (if a < -1 then pis else acos (u a))+ where+ I _ pis = pi+ atan = increasing atan+ sinh = increasing sinh+ cosh x@(I a b)+ | null x = empty+ | b < 0 = decreasing cosh x+ | a >= 0 = increasing cosh x+ | otherwise = I 0 $ cosh $ if - u a > b+ then u a + else b+ tanh = increasing tanh+ asinh = increasing asinh+ acosh x@(I a b)+ | null x || b < 1 = empty -- acosh is only defined on [1..1/0)+ | otherwise = I lo $ acosh b+ where lo | a <= 1 = 0 + | otherwise = acosh a+ atanh x@(I a b)+ | null x || b < -1 || a > 1 = empty+ | otherwise =+ (if a <= - 1 then negInfinity else atanh a)+ `I`+ (if b >= 1 then posInfinity else atanh b)+ +-- | lift a monotone increasing function over a given interval +increasing :: Precision a => + (forall d. Rounding d => Round d a -> Round d a) -> + Interval a -> Interval a+increasing f (I a b) = I (f a) (f b)++-- | lift a monotone increasing function over a given interval +decreasing :: Precision a => + (forall d. Rounding d => Round d a -> Round d a) -> + Interval a -> Interval a+decreasing f (I a b) = I (f (d b)) (f (u a))+++-- | We have to play some semantic games to make these methods make sense.+-- Most compute with the midpoint of the interval.+instance Precision a => RealFloat (Interval a) where+ floatRadix = floatRadix . midpoint+ floatDigits = floatDigits . midpoint+ floatRange = floatRange . midpoint+ decodeFloat = decodeFloat . midpoint+ encodeFloat m e = singleton (encodeFloat m e)+ exponent = exponent . midpoint+ significand x = min a b ... max a b+ where+ (_ ,em) = decodeFloat (midpoint x)+ (mi,ei) = decodeFloat (inf x)+ (ms,es) = decodeFloat (sup x)+ a = encodeFloat mi (ei - em - floatDigits x) + b = encodeFloat ms (es - em - floatDigits x)+ scaleFloat n x = scaleFloat n (inf x) ... scaleFloat n (sup x)+ isNaN x = isNaN (inf x) || isNaN (sup x)+ isInfinite x = isInfinite (inf x) || isInfinite (sup x)+ isDenormalized x = isDenormalized (inf x) || isDenormalized (sup x)+ -- contains negative zero+ isNegativeZero x = not (inf x > 0) + && not (sup x < 0)+ && ( (sup x == 0 && (inf x < 0 || isNegativeZero (inf x)))+ || (inf x == 0 && isNegativeZero (inf x)) + || (inf x < 0 && sup x >= 0))+ isIEEE x = isIEEE (inf x) && isIEEE (sup x)+ atan2 = error "unimplemented"++-- TODO: (^), (^^) to give tighter bounds+ +-- | Calculate the intersection of two intervals.+intersection :: Precision a => Interval a -> Interval a -> Interval a+intersection x@(I a b) y@(I a' b')+ | x /=! y = empty+ | otherwise = I (max a a') (min b b')++-- | Calculate the convex hull of two intervals+hull :: Ord a => Interval a -> Interval a -> Interval a+hull x@(I a b) y@(I a' b') + | null x = y+ | null y = x+ | otherwise = I (min a a') (max b b')+ +instance Precision a => RealExtras (Interval a) where+ type C (Interval a) = C a+ -- output always lies within the interval y if y >=! 0+ fmod x y | null y = empty + | otherwise = r -- `intersection` bounds+ where + n :: Integer+ n = floor (inf x / if inf x < 0 then inf y else sup y)+ r = x - fromIntegral n * y + -- bounds | inf y >= 0 = y+ -- | otherwise = y `hull` negate y+ expm1 = increasing expm1+ log1p (I a b) = (if a > (-1) then log1p a else negInfinity) `I` log1p b+ hypot x y = hypot a a' `I` hypot b b'+ where+ I a b = abs x+ I a' b' = abs y+ cbrt = increasing cbrt+ erf = increasing erf++-- | For all @x@ in @X@, @y@ in @Y@. @x '<' y@+(<!) :: Ord a => Interval a -> Interval a -> Bool+x <! y = sup x < inf y+{-# INLINE (<!) #-}++-- | For all @x@ in @X@, @y@ in @Y@. @x '<=' y@+(<=!) :: Ord a => Interval a -> Interval a -> Bool+x <=! y = sup x <= inf y+{-# INLINE (<=!) #-}++-- | For all @x@ in @X@, @y@ in @Y@. @x '==' y@+(==!) :: Eq a => Interval a -> Interval a -> Bool+x ==! y = sup x == inf y && inf x == sup y+{-# INLINE (==!) #-}++-- | For all @x@ in @X@, @y@ in @Y@. @x '/=' y@+(/=!) :: Ord a => Interval a -> Interval a -> Bool+x /=! y = sup x < inf y || inf x > sup y+{-# INLINE (/=!) #-}++-- | For all @x@ in @X@, @y@ in @Y@. @x '>' y@+(>!) :: Ord a => Interval a -> Interval a -> Bool+x >! y = inf x > sup y+{-# INLINE (>!) #-}++-- | For all @x@ in @X@, @y@ in @Y@. @x '>=' y@+(>=!) :: Ord a => Interval a -> Interval a -> Bool+x >=! y = inf x >= sup y+{-# INLINE (>=!) #-}++-- | For all @x@ in @X@, @y@ in @Y@. @x `op` y@+certainly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool+certainly cmp l r + | lt && eq && gt = True+ | lt && eq = l <=! r+ | lt && gt = l /=! r+ | lt = l <! r + | eq && gt = l >=! r + | eq = l ==! r+ | gt = l >! r+ | otherwise = False+ where + lt = cmp LT EQ+ eq = cmp EQ EQ+ gt = cmp GT EQ+{-# INLINE certainly #-}++contains :: Ord a => Interval a -> Interval a -> Bool+contains x y = null y + || (not (null x) && inf x <= inf y && sup y <= sup x)+{-# INLINE contains #-}++isSubsetOf :: Ord a => Interval a -> Interval a -> Bool+isSubsetOf = flip contains++-- | Comparisons are made on the midpoint+instance Precision a => Ord (Interval a) where+ compare = compare `on` midpoint+ max (I a b) (I a' b') = I (max a a') (max b b')+ min (I a b) (I a' b') = I (min a a') (min b b')++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<' y@?+(<?) :: Ord a => Interval a -> Interval a -> Bool+x <? y = inf x < sup y+{-# INLINE (<?) #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<=' y@?+(<=?) :: Ord a => Interval a -> Interval a -> Bool+x <=? y = inf x <= sup y+{-# INLINE (<=?) #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '==' y@?+(==?) :: Ord a => Interval a -> Interval a -> Bool+x ==? y = inf x <= sup y && sup x >= inf y+{-# INLINE (==?) #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '/=' y@?+(/=?) :: Eq a => Interval a -> Interval a -> Bool+x /=? y = inf x /= sup y || sup x /= inf y+{-# INLINE (/=?) #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>' y@?+(>?) :: Ord a => Interval a -> Interval a -> Bool+x >? y = sup x > inf y+{-# INLINE (>?) #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>=' y@?+(>=?) :: Ord a => Interval a -> Interval a -> Bool+x >=? y = sup x >= inf y+{-# INLINE (>=?) #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x `op` y@?+possibly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool+possibly cmp l r + | lt && eq && gt = True+ | lt && eq = l <=? r+ | lt && gt = l /=? r+ | lt = l <? r + | eq && gt = l >=? r + | eq = l ==? r+ | gt = l >? r+ | otherwise = False+ where + lt = cmp LT EQ+ eq = cmp EQ EQ+ gt = cmp GT EQ+{-# INLINE possibly #-}++idouble :: Interval Double -> Interval Double+idouble = id++ifloat :: Interval Float -> Interval Float+ifloat = id++-- Bugs:+-- sin 1 :: Interval Double
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMainWithHooks simpleUserHooks
+ intervals.cabal view
@@ -0,0 +1,19 @@+Name: intervals+Version: 0.1.0+Synopsis: Interval Arithmetic+Homepage: http://patch-tag.com/r/ekmett/intervals+License: BSD3+License-file: LICENSE+Author: Edward Kmett+Maintainer: ekmett@gmail.com+Category: Math+Build-type: Simple+Cabal-version: >=1.6++Library+ Exposed-modules: Numeric.Interval+ Build-depends: base >= 4 && < 5,+ array >= 0.3.0 && < 0.4,+ numeric-extras >= 0.0.1 && < 0.1,+ rounding >= 0.3.0 && < 0.4+ GHC-Options: -Wall