diff options
Diffstat (limited to 'lib/Numeric')
-rw-r--r-- | lib/Numeric/Interval.hs | 754 | ||||
-rw-r--r-- | lib/Numeric/Interval/Bounded.hs | 4 |
2 files changed, 2 insertions, 756 deletions
diff --git a/lib/Numeric/Interval.hs b/lib/Numeric/Interval.hs deleted file mode 100644 index df4bc33..0000000 --- a/lib/Numeric/Interval.hs +++ /dev/null | |||
@@ -1,754 +0,0 @@ | |||
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 | |||
21 | module 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 | |||
49 | import Control.Applicative hiding (empty) | ||
50 | import Data.Data | ||
51 | #ifdef VERSION_distributive | ||
52 | import Data.Distributive | ||
53 | #endif | ||
54 | import Data.Foldable hiding (minimum, maximum, elem, notElem, null) | ||
55 | import Data.Function (on) | ||
56 | import Data.Monoid | ||
57 | import Data.Traversable | ||
58 | #if defined(__GLASGOW_HASKELL) && __GLASGOW_HASKELL__ >= 704 | ||
59 | import GHC.Generics | ||
60 | #endif | ||
61 | import Prelude hiding (null, elem, notElem) | ||
62 | |||
63 | -- $setup | ||
64 | |||
65 | data 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 | |||
76 | instance Functor Interval where | ||
77 | fmap f (I a b) = I (f a) (f b) | ||
78 | {-# INLINE fmap #-} | ||
79 | |||
80 | instance Foldable Interval where | ||
81 | foldMap f (I a b) = f a `mappend` f b | ||
82 | {-# INLINE foldMap #-} | ||
83 | |||
84 | instance Traversable Interval where | ||
85 | traverse f (I a b) = I <$> f a <*> f b | ||
86 | {-# INLINE traverse #-} | ||
87 | |||
88 | instance 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 | |||
94 | instance 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 | ||
103 | instance Distributive Interval where | ||
104 | distribute f = fmap inf f ... fmap sup f | ||
105 | {-# INLINE distribute #-} | ||
106 | #endif | ||
107 | |||
108 | infix 3 ... | ||
109 | |||
110 | negInfinity :: Fractional a => a | ||
111 | negInfinity = (-1)/0 | ||
112 | {-# INLINE negInfinity #-} | ||
113 | |||
114 | posInfinity :: Fractional a => a | ||
115 | posInfinity = 1/0 | ||
116 | {-# INLINE posInfinity #-} | ||
117 | |||
118 | nan :: Fractional a => a | ||
119 | nan = 0/0 | ||
120 | |||
121 | fmod :: RealFrac a => a -> a -> a | ||
122 | fmod 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 | ||
136 | whole :: Fractional a => Interval a | ||
137 | whole = negInfinity ... posInfinity | ||
138 | {-# INLINE whole #-} | ||
139 | |||
140 | -- | An empty interval | ||
141 | -- | ||
142 | -- >>> empty | ||
143 | -- NaN ... NaN | ||
144 | empty :: Fractional a => Interval a | ||
145 | empty = 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 | ||
158 | null :: Ord a => Interval a -> Bool | ||
159 | null x = not (inf x <= sup x) | ||
160 | {-# INLINE null #-} | ||
161 | |||
162 | -- | A singleton point | ||
163 | -- | ||
164 | -- >>> singleton 1 | ||
165 | -- 1 ... 1 | ||
166 | singleton :: a -> Interval a | ||
167 | singleton a = a ... a | ||
168 | {-# INLINE singleton #-} | ||
169 | |||
170 | -- | The infinumum (lower bound) of an interval | ||
171 | -- | ||
172 | -- >>> inf (1 ... 20) | ||
173 | -- 1 | ||
174 | inf :: Interval a -> a | ||
175 | inf (I a _) = a | ||
176 | {-# INLINE inf #-} | ||
177 | |||
178 | -- | The supremum (upper bound) of an interval | ||
179 | -- | ||
180 | -- >>> sup (1 ... 20) | ||
181 | -- 20 | ||
182 | sup :: Interval a -> a | ||
183 | sup (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 | ||
195 | singular :: Ord a => Interval a -> Bool | ||
196 | singular x = not (null x) && inf x == sup x | ||
197 | {-# INLINE singular #-} | ||
198 | |||
199 | instance Eq a => Eq (Interval a) where | ||
200 | (==) = (==!) | ||
201 | {-# INLINE (==) #-} | ||
202 | |||
203 | instance 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 | ||
220 | width :: Num a => Interval a -> a | ||
221 | width (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 | ||
234 | magnitude :: (Num a, Ord a) => Interval a -> a | ||
235 | magnitude 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 | ||
248 | mignitude :: (Num a, Ord a) => Interval a -> a | ||
249 | mignitude x = (min `on` abs) (inf x) (sup x) | ||
250 | {-# INLINE mignitude #-} | ||
251 | |||
252 | instance (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) | ||
284 | bisection :: Fractional a => Interval a -> (Interval a, Interval a) | ||
285 | bisection 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 | ||
299 | midpoint :: Fractional a => Interval a -> a | ||
300 | midpoint 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 | -- | ||
320 | elem :: Ord a => a -> Interval a -> Bool | ||
321 | elem 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 | ||
336 | notElem :: Ord a => a -> Interval a -> Bool | ||
337 | notElem x xs = not (elem x xs) | ||
338 | {-# INLINE notElem #-} | ||
339 | |||
340 | -- | 'realToFrac' will use the midpoint | ||
341 | instance 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 | |||
350 | instance 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@ | ||
365 | divNonZero :: (Fractional a, Ord a) => Interval a -> Interval a -> Interval a | ||
366 | divNonZero (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] | ||
372 | divPositive :: (Fractional a, Ord a) => Interval a -> a -> Interval a | ||
373 | divPositive 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] | ||
382 | divNegative :: (Fractional a, Ord a) => Interval a -> a -> Interval a | ||
383 | divNegative 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 | |||
391 | divZero :: (Fractional a, Ord a) => Interval a -> Interval a | ||
392 | divZero x | ||
393 | | inf x == 0 && sup x == 0 = x | ||
394 | | otherwise = whole | ||
395 | {-# INLINE divZero #-} | ||
396 | |||
397 | instance (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 | |||
413 | instance 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 | |||
427 | instance (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 | ||
504 | increasing :: (a -> b) -> Interval a -> Interval b | ||
505 | increasing f (I a b) = f a ... f b | ||
506 | |||
507 | -- | lift a monotone decreasing function over a given interval | ||
508 | decreasing :: (a -> b) -> Interval a -> Interval b | ||
509 | decreasing 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. | ||
513 | instance 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 | ||
547 | intersection :: (Fractional a, Ord a) => Interval a -> Interval a -> Interval a | ||
548 | intersection 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 | ||
560 | hull :: Ord a => Interval a -> Interval a -> Interval a | ||
561 | hull 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 | ||
578 | x <! 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 | ||
592 | x <=! 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 | ||
605 | x ==! 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 | ||
616 | x /=! 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 | ||
627 | x >! 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 | ||
638 | x >=! y = inf x >= sup y | ||
639 | {-# INLINE (>=!) #-} | ||
640 | |||
641 | -- | For all @x@ in @X@, @y@ in @Y@. @x `op` y@ | ||
642 | -- | ||
643 | -- | ||
644 | certainly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool | ||
645 | certainly 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 | ||
667 | contains :: Ord a => Interval a -> Interval a -> Bool | ||
668 | contains 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 | ||
679 | isSubsetOf :: Ord a => Interval a -> Interval a -> Bool | ||
680 | isSubsetOf = 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 | ||
685 | x <? 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 | ||
690 | x <=? 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 | ||
695 | x ==? 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 | ||
700 | x /=? 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 | ||
705 | x >? 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 | ||
710 | x >=? y = sup x >= inf y | ||
711 | {-# INLINE (>=?) #-} | ||
712 | |||
713 | -- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x `op` y@? | ||
714 | possibly :: Ord a => (forall b. Ord b => b -> b -> Bool) -> Interval a -> Interval a -> Bool | ||
715 | possibly 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. | ||
731 | clamp :: Ord a => Interval a -> a -> a | ||
732 | clamp (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 | ||
740 | idouble :: Interval Double -> Interval Double | ||
741 | idouble = id | ||
742 | |||
743 | -- | id function. Useful for type specification | ||
744 | -- | ||
745 | -- >>> :t ifloat (1 ... 3) | ||
746 | -- ifloat (1 ... 3) :: Interval Float | ||
747 | ifloat :: Interval Float -> Interval Float | ||
748 | ifloat = id | ||
749 | |||
750 | -- Bugs: | ||
751 | -- sin 1 :: Interval Double | ||
752 | |||
753 | |||
754 | default (Integer,Double) | ||
diff --git a/lib/Numeric/Interval/Bounded.hs b/lib/Numeric/Interval/Bounded.hs index 2dd4d7b..a0d5d0d 100644 --- a/lib/Numeric/Interval/Bounded.hs +++ b/lib/Numeric/Interval/Bounded.hs | |||
@@ -2,8 +2,8 @@ module Numeric.Interval.Bounded where | |||
2 | 2 | ||
3 | import Numeric.Interval | 3 | import Numeric.Interval |
4 | 4 | ||
5 | whole' :: Bounded a => Interval a | 5 | whole' :: (Ord a, Bounded a) => Interval a |
6 | whole' = ( minBound ... maxBound ) | 6 | whole' = ( minBound ... maxBound ) |
7 | 7 | ||
8 | empty' :: Bounded a => Interval a | 8 | empty' :: (Ord a, Bounded a) => Interval a |
9 | empty' = ( maxBound ... minBound ) | 9 | empty' = ( maxBound ... minBound ) |