summaryrefslogtreecommitdiff
path: root/lib/Numeric/Interval.hs
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/Interval.hs')
-rw-r--r--lib/Numeric/Interval.hs754
1 files changed, 754 insertions, 0 deletions
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 @@
1{-# LANGUAGE CPP #-}
2{-# LANGUAGE Rank2Types #-}
3{-# LANGUAGE DeriveDataTypeable #-}
4#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
5{-# LANGUAGE DeriveGeneric #-}
6#endif
7-----------------------------------------------------------------------------
8-- |
9-- Module : Numeric.Interval
10-- Copyright : (c) Edward Kmett 2010-2013
11-- License : BSD3
12-- Maintainer : ekmett@gmail.com
13-- Stability : experimental
14-- Portability : DeriveDataTypeable
15-- Version : intervals-0.4.2 (minus distributive instance)
16--
17-- Interval arithmetic
18--
19-----------------------------------------------------------------------------
20
21module Numeric.Interval
22 ( Interval(..)
23 , (...)
24 , whole
25 , empty
26 , null
27 , singleton
28 , elem
29 , notElem
30 , inf
31 , sup
32 , singular
33 , width
34 , midpoint
35 , intersection
36 , hull
37 , bisection
38 , magnitude
39 , mignitude
40 , contains
41 , isSubsetOf
42 , certainly, (<!), (<=!), (==!), (>=!), (>!)
43 , possibly, (<?), (<=?), (==?), (>=?), (>?)
44 , clamp
45 , idouble
46 , ifloat
47 ) where
48
49import Control.Applicative hiding (empty)
50import Data.Data
51#ifdef VERSION_distributive
52import Data.Distributive
53#endif
54import Data.Foldable hiding (minimum, maximum, elem, notElem, null)
55import Data.Function (on)
56import Data.Monoid
57import Data.Traversable
58#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
59import GHC.Generics
60#endif
61import Prelude hiding (null, elem, notElem)
62
63-- $setup
64
65data Interval a = I !a !a deriving
66 ( Data
67 , Typeable
68#if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704
69 , Generic
70#if __GLASGOW_HASKELL__ >= 706
71 , Generic1
72#endif
73#endif
74 )
75
76instance Functor Interval where
77 fmap f (I a b) = I (f a) (f b)
78 {-# INLINE fmap #-}
79
80instance Foldable Interval where
81 foldMap f (I a b) = f a `mappend` f b
82 {-# INLINE foldMap #-}
83
84instance Traversable Interval where
85 traverse f (I a b) = I <$> f a <*> f b
86 {-# INLINE traverse #-}
87
88instance Applicative Interval where
89 pure a = I a a
90 {-# INLINE pure #-}
91 I f g <*> I a b = I (f a) (g b)
92 {-# INLINE (<*>) #-}
93
94instance Monad Interval where
95 return a = I a a
96 {-# INLINE return #-}
97 I a b >>= f = I a' b' where
98 I a' _ = f a
99 I _ b' = f b
100 {-# INLINE (>>=) #-}
101
102#ifdef VERSION_distributive
103instance Distributive Interval where
104 distribute f = fmap inf f ... fmap sup f
105 {-# INLINE distribute #-}
106#endif
107
108infix 3 ...
109
110negInfinity :: Fractional a => a
111negInfinity = (-1)/0
112{-# INLINE negInfinity #-}
113
114posInfinity :: Fractional a => a
115posInfinity = 1/0
116{-# INLINE posInfinity #-}
117
118nan :: Fractional a => a
119nan = 0/0
120
121fmod :: RealFrac a => a -> a -> a
122fmod a b = a - q*b where
123 q = realToFrac (truncate $ a / b :: Integer)
124{-# INLINE fmod #-}
125
126-- | The rule of thumb is you should only use this to construct using values
127-- that you took out of the interval. Otherwise, use I, to force rounding
128(...) :: a -> a -> Interval a
129(...) = I
130{-# INLINE (...) #-}
131
132-- | The whole real number line
133--
134-- >>> whole
135-- -Infinity ... Infinity
136whole :: Fractional a => Interval a
137whole = negInfinity ... posInfinity
138{-# INLINE whole #-}
139
140-- | An empty interval
141--
142-- >>> empty
143-- NaN ... NaN
144empty :: Fractional a => Interval a
145empty = nan ... nan
146{-# INLINE empty #-}
147
148-- | negation handles NaN properly
149--
150-- >>> null (1 ... 5)
151-- False
152--
153-- >>> null (1 ... 1)
154-- False
155--
156-- >>> null empty
157-- True
158null :: Ord a => Interval a -> Bool
159null x = not (inf x <= sup x)
160{-# INLINE null #-}
161
162-- | A singleton point
163--
164-- >>> singleton 1
165-- 1 ... 1
166singleton :: a -> Interval a
167singleton a = a ... a
168{-# INLINE singleton #-}
169
170-- | The infinumum (lower bound) of an interval
171--
172-- >>> inf (1 ... 20)
173-- 1
174inf :: Interval a -> a
175inf (I a _) = a
176{-# INLINE inf #-}
177
178-- | The supremum (upper bound) of an interval
179--
180-- >>> sup (1 ... 20)
181-- 20
182sup :: Interval a -> a
183sup (I _ b) = b
184{-# INLINE sup #-}
185
186-- | Is the interval a singleton point?
187-- N.B. This is fairly fragile and likely will not hold after
188-- even a few operations that only involve singletons
189--
190-- >>> singular (singleton 1)
191-- True
192--
193-- >>> singular (1.0 ... 20.0)
194-- False
195singular :: Ord a => Interval a -> Bool
196singular x = not (null x) && inf x == sup x
197{-# INLINE singular #-}
198
199instance Eq a => Eq (Interval a) where
200 (==) = (==!)
201 {-# INLINE (==) #-}
202
203instance Show a => Show (Interval a) where
204 showsPrec n (I a b) =
205 showParen (n > 3) $
206 showsPrec 3 a .
207 showString " ... " .
208 showsPrec 3 b
209
210-- | Calculate the width of an interval.
211--
212-- >>> width (1 ... 20)
213-- 19
214--
215-- >>> width (singleton 1)
216-- 0
217--
218-- >>> width empty
219-- NaN
220width :: Num a => Interval a -> a
221width (I a b) = b - a
222{-# INLINE width #-}
223
224-- | Magnitude
225--
226-- >>> magnitude (1 ... 20)
227-- 20
228--
229-- >>> magnitude (-20 ... 10)
230-- 20
231--
232-- >>> magnitude (singleton 5)
233-- 5
234magnitude :: (Num a, Ord a) => Interval a -> a
235magnitude x = (max `on` abs) (inf x) (sup x)
236{-# INLINE magnitude #-}
237
238-- | \"mignitude\"
239--
240-- >>> mignitude (1 ... 20)
241-- 1
242--
243-- >>> mignitude (-20 ... 10)
244-- 10
245--
246-- >>> mignitude (singleton 5)
247-- 5
248mignitude :: (Num a, Ord a) => Interval a -> a
249mignitude x = (min `on` abs) (inf x) (sup x)
250{-# INLINE mignitude #-}
251
252instance (Num a, Ord a) => Num (Interval a) where
253 I a b + I a' b' = (a + a') ... (b + b')
254 {-# INLINE (+) #-}
255 I a b - I a' b' = (a - b') ... (b - a')
256 {-# INLINE (-) #-}
257 I a b * I a' b' =
258 minimum [a * a', a * b', b * a', b * b']
259 ...
260 maximum [a * a', a * b', b * a', b * b']
261 {-# INLINE (*) #-}
262 abs x@(I a b)
263 | a >= 0 = x
264 | b <= 0 = negate x
265 | otherwise = 0 ... max (- a) b
266 {-# INLINE abs #-}
267
268 signum = increasing signum
269 {-# INLINE signum #-}
270
271 fromInteger i = singleton (fromInteger i)
272 {-# INLINE fromInteger #-}
273
274-- | Bisect an interval at its midpoint.
275--
276-- >>> bisection (10.0 ... 20.0)
277-- (10.0 ... 15.0,15.0 ... 20.0)
278--
279-- >>> bisection (singleton 5.0)
280-- (5.0 ... 5.0,5.0 ... 5.0)
281--
282-- >>> bisection empty
283-- (NaN ... NaN,NaN ... NaN)
284bisection :: Fractional a => Interval a -> (Interval a, Interval a)
285bisection x = (inf x ... m, m ... sup x)
286 where m = midpoint x
287{-# INLINE bisection #-}
288
289-- | Nearest point to the midpoint of the interval.
290--
291-- >>> midpoint (10.0 ... 20.0)
292-- 15.0
293--
294-- >>> midpoint (singleton 5.0)
295-- 5.0
296--
297-- >>> midpoint empty
298-- NaN
299midpoint :: Fractional a => Interval a -> a
300midpoint x = inf x + (sup x - inf x) / 2
301{-# INLINE midpoint #-}
302
303-- | Determine if a point is in the interval.
304--
305-- >>> elem 3.2 (1.0 ... 5.0)
306-- True
307--
308-- >>> elem 5 (1.0 ... 5.0)
309-- True
310--
311-- >>> elem 1 (1.0 ... 5.0)
312-- True
313--
314-- >>> elem 8 (1.0 ... 5.0)
315-- False
316--
317-- >>> elem 5 empty
318-- False
319--
320elem :: Ord a => a -> Interval a -> Bool
321elem x xs = x >= inf xs && x <= sup xs
322{-# INLINE elem #-}
323
324-- | Determine if a point is not included in the interval
325--
326-- >>> notElem 8 (1.0 ... 5.0)
327-- True
328--
329-- >>> notElem 1.4 (1.0 ... 5.0)
330-- False
331--
332-- And of course, nothing is a member of the empty interval.
333--
334-- >>> notElem 5 empty
335-- True
336notElem :: Ord a => a -> Interval a -> Bool
337notElem x xs = not (elem x xs)
338{-# INLINE notElem #-}
339
340-- | 'realToFrac' will use the midpoint
341instance Real a => Real (Interval a) where
342 toRational x
343 | null x = nan
344 | otherwise = a + (b - a) / 2
345 where
346 a = toRational (inf x)
347 b = toRational (sup x)
348 {-# INLINE toRational #-}
349
350instance Ord a => Ord (Interval a) where
351 compare x y
352 | sup x < inf y = LT
353 | inf x > sup y = GT
354 | sup x == inf y && inf x == sup y = EQ
355 | otherwise = error "Numeric.Interval.compare: ambiguous comparison"
356 {-# INLINE compare #-}
357
358 max (I a b) (I a' b') = max a a' ... max b b'
359 {-# INLINE max #-}
360
361 min (I a b) (I a' b') = min a a' ... min b b'
362 {-# INLINE min #-}
363
364-- @'divNonZero' X Y@ assumes @0 `'notElem'` Y@
365divNonZero :: (Fractional a, Ord a) => Interval a -> Interval a -> Interval a
366divNonZero (I a b) (I a' b') =
367 minimum [a / a', a / b', b / a', b / b']
368 ...
369 maximum [a / a', a / b', b / a', b / b']
370
371-- @'divPositive' X y@ assumes y > 0, and divides @X@ by [0 ... y]
372divPositive :: (Fractional a, Ord a) => Interval a -> a -> Interval a
373divPositive x@(I a b) y
374 | a == 0 && b == 0 = x
375 -- b < 0 || isNegativeZero b = negInfinity ... ( b / y)
376 | b < 0 = negInfinity ... ( b / y)
377 | a < 0 = whole
378 | otherwise = (a / y) ... posInfinity
379{-# INLINE divPositive #-}
380
381-- divNegative assumes y < 0 and divides the interval @X@ by [y ... 0]
382divNegative :: (Fractional a, Ord a) => Interval a -> a -> Interval a
383divNegative x@(I a b) y
384 | a == 0 && b == 0 = - x -- flip negative zeros
385 -- b < 0 || isNegativeZero b = (b / y) ... posInfinity
386 | b < 0 = (b / y) ... posInfinity
387 | a < 0 = whole
388 | otherwise = negInfinity ... (a / y)
389{-# INLINE divNegative #-}
390
391divZero :: (Fractional a, Ord a) => Interval a -> Interval a
392divZero x
393 | inf x == 0 && sup x == 0 = x
394 | otherwise = whole
395{-# INLINE divZero #-}
396
397instance (Fractional a, Ord a) => Fractional (Interval a) where
398 -- TODO: check isNegativeZero properly
399 x / y
400 | 0 `notElem` y = divNonZero x y
401 | iz && sz = empty -- division by 0
402 | iz = divPositive x (inf y)
403 | sz = divNegative x (sup y)
404 | otherwise = divZero x
405 where
406 iz = inf y == 0
407 sz = sup y == 0
408 recip (I a b) = on min recip a b ... on max recip a b
409 {-# INLINE recip #-}
410 fromRational r = let r' = fromRational r in r' ... r'
411 {-# INLINE fromRational #-}
412
413instance RealFrac a => RealFrac (Interval a) where
414 properFraction x = (b, x - fromIntegral b)
415 where
416 b = truncate (midpoint x)
417 {-# INLINE properFraction #-}
418 ceiling x = ceiling (sup x)
419 {-# INLINE ceiling #-}
420 floor x = floor (inf x)
421 {-# INLINE floor #-}
422 round x = round (midpoint x)
423 {-# INLINE round #-}
424 truncate x = truncate (midpoint x)
425 {-# INLINE truncate #-}
426
427instance (RealFloat a, Ord a) => Floating (Interval a) where
428 pi = singleton pi
429 {-# INLINE pi #-}
430 exp = increasing exp
431 {-# INLINE exp #-}
432 log (I a b) = (if a > 0 then log a else negInfinity) ... log b
433 {-# INLINE log #-}
434 cos x
435 | null x = empty
436 | width t >= pi = (-1) ... 1
437 | inf t >= pi = - cos (t - pi)
438 | sup t <= pi = decreasing cos t
439 | sup t <= 2 * pi = (-1) ... cos ((pi * 2 - sup t) `min` inf t)
440 | otherwise = (-1) ... 1
441 where
442 t = fmod x (pi * 2)
443 {-# INLINE cos #-}
444 sin x
445 | null x = empty
446 | otherwise = cos (x - pi / 2)
447 {-# INLINE sin #-}
448 tan x
449 | null x = empty
450 | inf t' <= - pi / 2 || sup t' >= pi / 2 = whole
451 | otherwise = increasing tan x
452 where
453 t = x `fmod` pi
454 t' | t >= pi / 2 = t - pi
455 | otherwise = t
456 {-# INLINE tan #-}
457 asin x@(I a b)
458 | null x || b < -1 || a > 1 = empty
459 | otherwise =
460 (if a <= -1 then -halfPi else asin a)
461 ...
462 (if b >= 1 then halfPi else asin b)
463 where
464 halfPi = pi / 2
465 {-# INLINE asin #-}
466 acos x@(I a b)
467 | null x || b < -1 || a > 1 = empty
468 | otherwise =
469 (if b >= 1 then 0 else acos b)
470 ...
471 (if a < -1 then pi else acos a)
472 {-# INLINE acos #-}
473 atan = increasing atan
474 {-# INLINE atan #-}
475 sinh = increasing sinh
476 {-# INLINE sinh #-}
477 cosh x@(I a b)
478 | null x = empty
479 | b < 0 = decreasing cosh x
480 | a >= 0 = increasing cosh x
481 | otherwise = I 0 $ cosh $ if - a > b
482 then a
483 else b
484 {-# INLINE cosh #-}
485 tanh = increasing tanh
486 {-# INLINE tanh #-}
487 asinh = increasing asinh
488 {-# INLINE asinh #-}
489 acosh x@(I a b)
490 | null x || b < 1 = empty
491 | otherwise = I lo $ acosh b
492 where lo | a <= 1 = 0
493 | otherwise = acosh a
494 {-# INLINE acosh #-}
495 atanh x@(I a b)
496 | null x || b < -1 || a > 1 = empty
497 | otherwise =
498 (if a <= - 1 then negInfinity else atanh a)
499 ...
500 (if b >= 1 then posInfinity else atanh b)
501 {-# INLINE atanh #-}
502
503-- | lift a monotone increasing function over a given interval
504increasing :: (a -> b) -> Interval a -> Interval b
505increasing f (I a b) = f a ... f b
506
507-- | lift a monotone decreasing function over a given interval
508decreasing :: (a -> b) -> Interval a -> Interval b
509decreasing f (I a b) = f b ... f a
510
511-- | We have to play some semantic games to make these methods make sense.
512-- Most compute with the midpoint of the interval.
513instance RealFloat a => RealFloat (Interval a) where
514 floatRadix = floatRadix . midpoint
515
516 floatDigits = floatDigits . midpoint
517 floatRange = floatRange . midpoint
518 decodeFloat = decodeFloat . midpoint
519 encodeFloat m e = singleton (encodeFloat m e)
520 exponent = exponent . midpoint
521 significand x = min a b ... max a b
522 where
523 (_ ,em) = decodeFloat (midpoint x)
524 (mi,ei) = decodeFloat (inf x)
525 (ms,es) = decodeFloat (sup x)
526 a = encodeFloat mi (ei - em - floatDigits x)
527 b = encodeFloat ms (es - em - floatDigits x)
528 scaleFloat n x = scaleFloat n (inf x) ... scaleFloat n (sup x)
529 isNaN x = isNaN (inf x) || isNaN (sup x)
530 isInfinite x = isInfinite (inf x) || isInfinite (sup x)
531 isDenormalized x = isDenormalized (inf x) || isDenormalized (sup x)
532 -- contains negative zero
533 isNegativeZero x = not (inf x > 0)
534 && not (sup x < 0)
535 && ( (sup x == 0 && (inf x < 0 || isNegativeZero (inf x)))
536 || (inf x == 0 && isNegativeZero (inf x))
537 || (inf x < 0 && sup x >= 0))
538 isIEEE x = isIEEE (inf x) && isIEEE (sup x)
539 atan2 = error "unimplemented"
540
541-- TODO: (^), (^^) to give tighter bounds
542
543-- | Calculate the intersection of two intervals.
544--
545-- >>> intersection (1 ... 10 :: Interval Double) (5 ... 15 :: Interval Double)
546-- 5.0 ... 10.0
547intersection :: (Fractional a, Ord a) => Interval a -> Interval a -> Interval a
548intersection x@(I a b) y@(I a' b')
549 | x /=! y = empty
550 | otherwise = max a a' ... min b b'
551{-# INLINE intersection #-}
552
553-- | Calculate the convex hull of two intervals
554--
555-- >>> hull (0 ... 10 :: Interval Double) (5 ... 15 :: Interval Double)
556-- 0.0 ... 15.0
557--
558-- >>> hull (15 ... 85 :: Interval Double) (0 ... 10 :: Interval Double)
559-- 0.0 ... 85.0
560hull :: Ord a => Interval a -> Interval a -> Interval a
561hull x@(I a b) y@(I a' b')
562 | null x = y
563 | null y = x
564 | otherwise = min a a' ... max b b'
565{-# INLINE hull #-}
566
567-- | For all @x@ in @X@, @y@ in @Y@. @x '<' y@
568--
569-- >>> (5 ... 10 :: Interval Double) <! (20 ... 30 :: Interval Double)
570-- True
571--
572-- >>> (5 ... 10 :: Interval Double) <! (10 ... 30 :: Interval Double)
573-- False
574--
575-- >>> (20 ... 30 :: Interval Double) <! (5 ... 10 :: Interval Double)
576-- False
577(<!) :: Ord a => Interval a -> Interval a -> Bool
578x <! y = sup x < inf y
579{-# INLINE (<!) #-}
580
581-- | For all @x@ in @X@, @y@ in @Y@. @x '<=' y@
582--
583-- >>> (5 ... 10 :: Interval Double) <=! (20 ... 30 :: Interval Double)
584-- True
585--
586-- >>> (5 ... 10 :: Interval Double) <=! (10 ... 30 :: Interval Double)
587-- True
588--
589-- >>> (20 ... 30 :: Interval Double) <=! (5 ... 10 :: Interval Double)
590-- False
591(<=!) :: Ord a => Interval a -> Interval a -> Bool
592x <=! y = sup x <= inf y
593{-# INLINE (<=!) #-}
594
595-- | For all @x@ in @X@, @y@ in @Y@. @x '==' y@
596--
597-- Only singleton intervals return true
598--
599-- >>> (singleton 5 :: Interval Double) ==! (singleton 5 :: Interval Double)
600-- True
601--
602-- >>> (5 ... 10 :: Interval Double) ==! (5 ... 10 :: Interval Double)
603-- False
604(==!) :: Eq a => Interval a -> Interval a -> Bool
605x ==! y = sup x == inf y && inf x == sup y
606{-# INLINE (==!) #-}
607
608-- | For all @x@ in @X@, @y@ in @Y@. @x '/=' y@
609--
610-- >>> (5 ... 15 :: Interval Double) /=! (20 ... 40 :: Interval Double)
611-- True
612--
613-- >>> (5 ... 15 :: Interval Double) /=! (15 ... 40 :: Interval Double)
614-- False
615(/=!) :: Ord a => Interval a -> Interval a -> Bool
616x /=! y = sup x < inf y || inf x > sup y
617{-# INLINE (/=!) #-}
618
619-- | For all @x@ in @X@, @y@ in @Y@. @x '>' y@
620--
621-- >>> (20 ... 40 :: Interval Double) >! (10 ... 19 :: Interval Double)
622-- True
623--
624-- >>> (5 ... 20 :: Interval Double) >! (15 ... 40 :: Interval Double)
625-- False
626(>!) :: Ord a => Interval a -> Interval a -> Bool
627x >! y = inf x > sup y
628{-# INLINE (>!) #-}
629
630-- | For all @x@ in @X@, @y@ in @Y@. @x '>=' y@
631--
632-- >>> (20 ... 40 :: Interval Double) >=! (10 ... 20 :: Interval Double)
633-- True
634--
635-- >>> (5 ... 20 :: Interval Double) >=! (15 ... 40 :: Interval Double)
636-- False
637(>=!) :: Ord a => Interval a -> Interval a -> Bool
638x >=! y = inf x >= sup y
639{-# INLINE (>=!) #-}
640
641-- | For all @x@ in @X@, @y@ in @Y@. @x `op` y@
642--
643--
644certainly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool
645certainly cmp l r
646 | lt && eq && gt = True
647 | lt && eq = l <=! r
648 | lt && gt = l /=! r
649 | lt = l <! r
650 | eq && gt = l >=! r
651 | eq = l ==! r
652 | gt = l >! r
653 | otherwise = False
654 where
655 lt = cmp LT EQ
656 eq = cmp EQ EQ
657 gt = cmp GT EQ
658{-# INLINE certainly #-}
659
660-- | Check if interval @X@ totally contains interval @Y@
661--
662-- >>> (20 ... 40 :: Interval Double) `contains` (25 ... 35 :: Interval Double)
663-- True
664--
665-- >>> (20 ... 40 :: Interval Double) `contains` (15 ... 35 :: Interval Double)
666-- False
667contains :: Ord a => Interval a -> Interval a -> Bool
668contains x y = null y
669 || (not (null x) && inf x <= inf y && sup y <= sup x)
670{-# INLINE contains #-}
671
672-- | Flipped version of `contains`. Check if interval @X@ a subset of interval @Y@
673--
674-- >>> (25 ... 35 :: Interval Double) `isSubsetOf` (20 ... 40 :: Interval Double)
675-- True
676--
677-- >>> (20 ... 40 :: Interval Double) `isSubsetOf` (15 ... 35 :: Interval Double)
678-- False
679isSubsetOf :: Ord a => Interval a -> Interval a -> Bool
680isSubsetOf = flip contains
681{-# INLINE isSubsetOf #-}
682
683-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<' y@?
684(<?) :: Ord a => Interval a -> Interval a -> Bool
685x <? y = inf x < sup y
686{-# INLINE (<?) #-}
687
688-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<=' y@?
689(<=?) :: Ord a => Interval a -> Interval a -> Bool
690x <=? y = inf x <= sup y
691{-# INLINE (<=?) #-}
692
693-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '==' y@?
694(==?) :: Ord a => Interval a -> Interval a -> Bool
695x ==? y = inf x <= sup y && sup x >= inf y
696{-# INLINE (==?) #-}
697
698-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '/=' y@?
699(/=?) :: Eq a => Interval a -> Interval a -> Bool
700x /=? y = inf x /= sup y || sup x /= inf y
701{-# INLINE (/=?) #-}
702
703-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>' y@?
704(>?) :: Ord a => Interval a -> Interval a -> Bool
705x >? y = sup x > inf y
706{-# INLINE (>?) #-}
707
708-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>=' y@?
709(>=?) :: Ord a => Interval a -> Interval a -> Bool
710x >=? y = sup x >= inf y
711{-# INLINE (>=?) #-}
712
713-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x `op` y@?
714possibly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool
715possibly cmp l r
716 | lt && eq && gt = True
717 | lt && eq = l <=? r
718 | lt && gt = l /=? r
719 | lt = l <? r
720 | eq && gt = l >=? r
721 | eq = l ==? r
722 | gt = l >? r
723 | otherwise = False
724 where
725 lt = cmp LT EQ
726 eq = cmp EQ EQ
727 gt = cmp GT EQ
728{-# INLINE possibly #-}
729
730-- | The nearest value to that supplied which is contained in the interval.
731clamp :: Ord a => Interval a -> a -> a
732clamp (I a b) x | x < a = a
733 | x > b = b
734 | otherwise = x
735
736-- | id function. Useful for type specification
737--
738-- >>> :t idouble (1 ... 3)
739-- idouble (1 ... 3) :: Interval Double
740idouble :: Interval Double -> Interval Double
741idouble = id
742
743-- | id function. Useful for type specification
744--
745-- >>> :t ifloat (1 ... 3)
746-- ifloat (1 ... 3) :: Interval Float
747ifloat :: Interval Float -> Interval Float
748ifloat = id
749
750-- Bugs:
751-- sin 1 :: Interval Double
752
753
754default (Integer,Double)