cb6e9e973e33ff931aca95f16e47b111a98eab4c
[ghc.git] / libraries / base / GHC / Real.hs
1 {-# LANGUAGE Trustworthy #-}
2 {-# LANGUAGE CPP, NoImplicitPrelude, MagicHash, UnboxedTuples, BangPatterns #-}
3 {-# OPTIONS_GHC -Wno-orphans #-}
4 {-# OPTIONS_HADDOCK not-home #-}
5
6 -----------------------------------------------------------------------------
7 -- |
8 -- Module : GHC.Real
9 -- Copyright : (c) The University of Glasgow, 1994-2002
10 -- License : see libraries/base/LICENSE
11 --
12 -- Maintainer : cvs-ghc@haskell.org
13 -- Stability : internal
14 -- Portability : non-portable (GHC Extensions)
15 --
16 -- The types 'Ratio' and 'Rational', and the classes 'Real', 'Fractional',
17 -- 'Integral', and 'RealFrac'.
18 --
19 -----------------------------------------------------------------------------
20
21 module GHC.Real where
22
23 #include "MachDeps.h"
24
25 import GHC.Base
26 import GHC.Num
27 import GHC.List
28 import GHC.Enum
29 import GHC.Show
30 import {-# SOURCE #-} GHC.Exception( divZeroException, overflowException
31 , underflowException
32 , ratioZeroDenomException )
33
34 #if defined(OPTIMISE_INTEGER_GCD_LCM)
35 # if defined(MIN_VERSION_integer_gmp)
36 import GHC.Integer.GMP.Internals
37 # else
38 # error unsupported OPTIMISE_INTEGER_GCD_LCM configuration
39 # endif
40 #endif
41
42 infixr 8 ^, ^^
43 infixl 7 /, `quot`, `rem`, `div`, `mod`
44 infixl 7 %
45
46 default () -- Double isn't available yet,
47 -- and we shouldn't be using defaults anyway
48
49 ------------------------------------------------------------------------
50 -- Divide by zero and arithmetic overflow
51 ------------------------------------------------------------------------
52
53 -- We put them here because they are needed relatively early
54 -- in the libraries before the Exception type has been defined yet.
55
56 {-# NOINLINE divZeroError #-}
57 divZeroError :: a
58 divZeroError = raise# divZeroException
59
60 {-# NOINLINE ratioZeroDenominatorError #-}
61 ratioZeroDenominatorError :: a
62 ratioZeroDenominatorError = raise# ratioZeroDenomException
63
64 {-# NOINLINE overflowError #-}
65 overflowError :: a
66 overflowError = raise# overflowException
67
68 {-# NOINLINE underflowError #-}
69 underflowError :: a
70 underflowError = raise# underflowException
71
72
73 --------------------------------------------------------------
74 -- The Ratio and Rational types
75 --------------------------------------------------------------
76
77 -- | Rational numbers, with numerator and denominator of some 'Integral' type.
78 --
79 -- Note that `Ratio`'s instances inherit the deficiencies from the type
80 -- parameter's. For example, @Ratio Natural@'s 'Num' instance has similar
81 -- problems to `Numeric.Natural.Natural`'s.
82 data Ratio a = !a :% !a deriving Eq -- ^ @since 2.01
83
84 -- | Arbitrary-precision rational numbers, represented as a ratio of
85 -- two 'Integer' values. A rational number may be constructed using
86 -- the '%' operator.
87 type Rational = Ratio Integer
88
89 ratioPrec, ratioPrec1 :: Int
90 ratioPrec = 7 -- Precedence of ':%' constructor
91 ratioPrec1 = ratioPrec + 1
92
93 infinity, notANumber :: Rational
94 infinity = 1 :% 0
95 notANumber = 0 :% 0
96
97 -- Use :%, not % for Inf/NaN; the latter would
98 -- immediately lead to a runtime error, because it normalises.
99
100 -- | Forms the ratio of two integral numbers.
101 {-# SPECIALISE (%) :: Integer -> Integer -> Rational #-}
102 (%) :: (Integral a) => a -> a -> Ratio a
103
104 -- | Extract the numerator of the ratio in reduced form:
105 -- the numerator and denominator have no common factor and the denominator
106 -- is positive.
107 numerator :: Ratio a -> a
108
109 -- | Extract the denominator of the ratio in reduced form:
110 -- the numerator and denominator have no common factor and the denominator
111 -- is positive.
112 denominator :: Ratio a -> a
113
114
115 -- | 'reduce' is a subsidiary function used only in this module.
116 -- It normalises a ratio by dividing both numerator and denominator by
117 -- their greatest common divisor.
118 reduce :: (Integral a) => a -> a -> Ratio a
119 {-# SPECIALISE reduce :: Integer -> Integer -> Rational #-}
120 reduce _ 0 = ratioZeroDenominatorError
121 reduce x y = (x `quot` d) :% (y `quot` d)
122 where d = gcd x y
123
124 x % y = reduce (x * signum y) (abs y)
125
126 numerator (x :% _) = x
127 denominator (_ :% y) = y
128
129 --------------------------------------------------------------
130 -- Standard numeric classes
131 --------------------------------------------------------------
132
133 class (Num a, Ord a) => Real a where
134 -- | the rational equivalent of its real argument with full precision
135 toRational :: a -> Rational
136
137 -- | Integral numbers, supporting integer division.
138 --
139 -- The Haskell Report defines no laws for 'Integral'. However, 'Integral'
140 -- instances are customarily expected to define a Euclidean domain and have the
141 -- following properties for the `div`\/`mod` and `quot`\/`rem` pairs, given
142 -- suitable Euclidean functions @f@ and @g@:
143 --
144 -- * @x@ = @y * quot x y + rem x y@ with @rem x y@ = @fromInteger 0@ or
145 -- @g (rem x y)@ < @g y@
146 -- * @x@ = @y * div x y + mod x y@ with @mod x y@ = @fromInteger 0@ or
147 -- @f (mod x y)@ < @f y@
148 --
149 -- An example of a suitable Euclidean function, for `Integer`'s instance, is
150 -- 'abs'.
151 class (Real a, Enum a) => Integral a where
152 -- | integer division truncated toward zero
153 quot :: a -> a -> a
154 -- | integer remainder, satisfying
155 --
156 -- > (x `quot` y)*y + (x `rem` y) == x
157 rem :: a -> a -> a
158 -- | integer division truncated toward negative infinity
159 div :: a -> a -> a
160 -- | integer modulus, satisfying
161 --
162 -- > (x `div` y)*y + (x `mod` y) == x
163 mod :: a -> a -> a
164 -- | simultaneous 'quot' and 'rem'
165 quotRem :: a -> a -> (a,a)
166 -- | simultaneous 'div' and 'mod'
167 divMod :: a -> a -> (a,a)
168 -- | conversion to 'Integer'
169 toInteger :: a -> Integer
170
171 {-# INLINE quot #-}
172 {-# INLINE rem #-}
173 {-# INLINE div #-}
174 {-# INLINE mod #-}
175 n `quot` d = q where (q,_) = quotRem n d
176 n `rem` d = r where (_,r) = quotRem n d
177 n `div` d = q where (q,_) = divMod n d
178 n `mod` d = r where (_,r) = divMod n d
179
180 divMod n d = if signum r == negate (signum d) then (q-1, r+d) else qr
181 where qr@(q,r) = quotRem n d
182
183 -- | Fractional numbers, supporting real division.
184 --
185 -- The Haskell Report defines no laws for 'Fractional'. However, @('+')@ and
186 -- @('*')@ are customarily expected to define a division ring and have the
187 -- following properties:
188 --
189 -- [__'recip' gives the multiplicative inverse__]:
190 -- @x * recip x@ = @recip x * x@ = @fromInteger 1@
191 --
192 -- Note that it /isn't/ customarily expected that a type instance of
193 -- 'Fractional' implement a field. However, all instances in @base@ do.
194 class (Num a) => Fractional a where
195 {-# MINIMAL fromRational, (recip | (/)) #-}
196
197 -- | Fractional division.
198 (/) :: a -> a -> a
199 -- | Reciprocal fraction.
200 recip :: a -> a
201 -- | Conversion from a 'Rational' (that is @'Ratio' 'Integer'@).
202 -- A floating literal stands for an application of 'fromRational'
203 -- to a value of type 'Rational', so such literals have type
204 -- @('Fractional' a) => a@.
205 fromRational :: Rational -> a
206
207 {-# INLINE recip #-}
208 {-# INLINE (/) #-}
209 recip x = 1 / x
210 x / y = x * recip y
211
212 -- | Extracting components of fractions.
213 class (Real a, Fractional a) => RealFrac a where
214 -- | The function 'properFraction' takes a real fractional number @x@
215 -- and returns a pair @(n,f)@ such that @x = n+f@, and:
216 --
217 -- * @n@ is an integral number with the same sign as @x@; and
218 --
219 -- * @f@ is a fraction with the same type and sign as @x@,
220 -- and with absolute value less than @1@.
221 --
222 -- The default definitions of the 'ceiling', 'floor', 'truncate'
223 -- and 'round' functions are in terms of 'properFraction'.
224 properFraction :: (Integral b) => a -> (b,a)
225 -- | @'truncate' x@ returns the integer nearest @x@ between zero and @x@
226 truncate :: (Integral b) => a -> b
227 -- | @'round' x@ returns the nearest integer to @x@;
228 -- the even integer if @x@ is equidistant between two integers
229 round :: (Integral b) => a -> b
230 -- | @'ceiling' x@ returns the least integer not less than @x@
231 ceiling :: (Integral b) => a -> b
232 -- | @'floor' x@ returns the greatest integer not greater than @x@
233 floor :: (Integral b) => a -> b
234
235 {-# INLINE truncate #-}
236 truncate x = m where (m,_) = properFraction x
237
238 round x = let (n,r) = properFraction x
239 m = if r < 0 then n - 1 else n + 1
240 in case signum (abs r - 0.5) of
241 -1 -> n
242 0 -> if even n then n else m
243 1 -> m
244 _ -> errorWithoutStackTrace "round default defn: Bad value"
245
246 ceiling x = if r > 0 then n + 1 else n
247 where (n,r) = properFraction x
248
249 floor x = if r < 0 then n - 1 else n
250 where (n,r) = properFraction x
251
252 -- These 'numeric' enumerations come straight from the Report
253
254 numericEnumFrom :: (Fractional a) => a -> [a]
255 numericEnumFrom n = go 0
256 where
257 -- See Note [Numeric Stability of Enumerating Floating Numbers]
258 go !k = let !n' = n + k
259 in n' : go (k + 1)
260
261 numericEnumFromThen :: (Fractional a) => a -> a -> [a]
262 numericEnumFromThen n m = go 0
263 where
264 step = m - n
265 -- See Note [Numeric Stability of Enumerating Floating Numbers]
266 go !k = let !n' = n + k * step
267 in n' : go (k + 1)
268
269 numericEnumFromTo :: (Ord a, Fractional a) => a -> a -> [a]
270 numericEnumFromTo n m = takeWhile (<= m + 1/2) (numericEnumFrom n)
271
272 numericEnumFromThenTo :: (Ord a, Fractional a) => a -> a -> a -> [a]
273 numericEnumFromThenTo e1 e2 e3
274 = takeWhile predicate (numericEnumFromThen e1 e2)
275 where
276 mid = (e2 - e1) / 2
277 predicate | e2 >= e1 = (<= e3 + mid)
278 | otherwise = (>= e3 + mid)
279
280 {- Note [Numeric Stability of Enumerating Floating Numbers]
281 -----------------------------------------------------------
282 When enumerate floating numbers, we could add the increment to the last number
283 at every run (as what we did previously):
284
285 numericEnumFrom n = n `seq` (n : numericEnumFrom (n + 1))
286
287 This approach is concise and really fast, only needs an addition operation.
288 However when a floating number is large enough, for `n`, `n` and `n+1` will
289 have the same binary representation. For example (all number has type
290 `Double`):
291
292 9007199254740990 is: 0x433ffffffffffffe
293 9007199254740990 + 1 is: 0x433fffffffffffff
294 (9007199254740990 + 1) + 1 is: 0x4340000000000000
295 ((9007199254740990 + 1) + 1) + 1 is: 0x4340000000000000
296
297 When we evaluate ([9007199254740990..9007199254740991] :: Double), we would
298 never reach the condition in `numericEnumFromTo`
299
300 9007199254740990 + 1 + 1 + ... > 9007199254740991 + 1/2
301
302 We would fall into infinite loop (as reported in Trac #15081).
303
304 To remedy the situation, we record the number of `1` that needed to be added
305 to the start number, rather than increasing `1` at every time. This approach
306 can improvement the numeric stability greatly at the cost of a multiplication.
307
308 Furthermore, we use the type of the enumerated number, `Fractional a => a`,
309 as the type of multiplier. In rare situations, the multiplier could be very
310 large and will lead to the enumeration to infinite loop, too, which should
311 be very rare. Consider the following example:
312
313 [1..9007199254740994]
314
315 We could fix that by using an Integer as multiplier but we don't do that.
316 The benchmark on T7954.hs shows that this approach leads to significant
317 degeneration on performance (33% increase allocation and 300% increase on
318 elapsed time).
319
320 See Trac #15081 and Phab:D4650 for the related discussion about this problem.
321 -}
322
323 --------------------------------------------------------------
324 -- Instances for Int
325 --------------------------------------------------------------
326
327 -- | @since 2.0.1
328 instance Real Int where
329 toRational x = toInteger x :% 1
330
331 -- | @since 2.0.1
332 instance Integral Int where
333 toInteger (I# i) = smallInteger i
334
335 a `quot` b
336 | b == 0 = divZeroError
337 | b == (-1) && a == minBound = overflowError -- Note [Order of tests]
338 -- in GHC.Int
339 | otherwise = a `quotInt` b
340
341 a `rem` b
342 | b == 0 = divZeroError
343 -- The quotRem CPU instruction fails for minBound `quotRem` -1,
344 -- but minBound `rem` -1 is well-defined (0). We therefore
345 -- special-case it.
346 | b == (-1) = 0
347 | otherwise = a `remInt` b
348
349 a `div` b
350 | b == 0 = divZeroError
351 | b == (-1) && a == minBound = overflowError -- Note [Order of tests]
352 -- in GHC.Int
353 | otherwise = a `divInt` b
354
355 a `mod` b
356 | b == 0 = divZeroError
357 -- The divMod CPU instruction fails for minBound `divMod` -1,
358 -- but minBound `mod` -1 is well-defined (0). We therefore
359 -- special-case it.
360 | b == (-1) = 0
361 | otherwise = a `modInt` b
362
363 a `quotRem` b
364 | b == 0 = divZeroError
365 -- Note [Order of tests] in GHC.Int
366 | b == (-1) && a == minBound = (overflowError, 0)
367 | otherwise = a `quotRemInt` b
368
369 a `divMod` b
370 | b == 0 = divZeroError
371 -- Note [Order of tests] in GHC.Int
372 | b == (-1) && a == minBound = (overflowError, 0)
373 | otherwise = a `divModInt` b
374
375 --------------------------------------------------------------
376 -- Instances for @Word@
377 --------------------------------------------------------------
378
379 -- | @since 2.01
380 instance Real Word where
381 toRational x = toInteger x % 1
382
383 -- | @since 2.01
384 instance Integral Word where
385 quot (W# x#) y@(W# y#)
386 | y /= 0 = W# (x# `quotWord#` y#)
387 | otherwise = divZeroError
388 rem (W# x#) y@(W# y#)
389 | y /= 0 = W# (x# `remWord#` y#)
390 | otherwise = divZeroError
391 div (W# x#) y@(W# y#)
392 | y /= 0 = W# (x# `quotWord#` y#)
393 | otherwise = divZeroError
394 mod (W# x#) y@(W# y#)
395 | y /= 0 = W# (x# `remWord#` y#)
396 | otherwise = divZeroError
397 quotRem (W# x#) y@(W# y#)
398 | y /= 0 = case x# `quotRemWord#` y# of
399 (# q, r #) ->
400 (W# q, W# r)
401 | otherwise = divZeroError
402 divMod (W# x#) y@(W# y#)
403 | y /= 0 = (W# (x# `quotWord#` y#), W# (x# `remWord#` y#))
404 | otherwise = divZeroError
405 toInteger (W# x#) = wordToInteger x#
406
407 --------------------------------------------------------------
408 -- Instances for Integer
409 --------------------------------------------------------------
410
411 -- | @since 2.0.1
412 instance Real Integer where
413 toRational x = x :% 1
414
415 #if defined(MIN_VERSION_integer_gmp)
416 -- | @since 4.8.0.0
417 instance Real Natural where
418 toRational (NatS# w) = toRational (W# w)
419 toRational (NatJ# bn) = toRational (Jp# bn)
420 #else
421 -- | @since 4.8.0.0
422 instance Real Natural where
423 toRational (Natural a) = toRational a
424 {-# INLINE toRational #-}
425 #endif
426
427 -- Note [Integer division constant folding]
428 -- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
429 --
430 -- Constant folding of quot, rem, div, mod, divMod and quotRem for
431 -- Integer arguments depends crucially on inlining. Constant folding
432 -- rules defined in compiler/prelude/PrelRules.hs trigger for
433 -- quotInteger, remInteger and so on. So if calls to quot, rem and so on
434 -- were not inlined the rules would not fire. The rules would also not
435 -- fire if calls to quotInteger and so on were inlined, but this does not
436 -- happen because they are all marked with NOINLINE pragma - see documentation
437 -- of integer-gmp or integer-simple.
438
439 -- | @since 2.0.1
440 instance Integral Integer where
441 toInteger n = n
442
443 {-# INLINE quot #-}
444 _ `quot` 0 = divZeroError
445 n `quot` d = n `quotInteger` d
446
447 {-# INLINE rem #-}
448 _ `rem` 0 = divZeroError
449 n `rem` d = n `remInteger` d
450
451 {-# INLINE div #-}
452 _ `div` 0 = divZeroError
453 n `div` d = n `divInteger` d
454
455 {-# INLINE mod #-}
456 _ `mod` 0 = divZeroError
457 n `mod` d = n `modInteger` d
458
459 {-# INLINE divMod #-}
460 _ `divMod` 0 = divZeroError
461 n `divMod` d = case n `divModInteger` d of
462 (# x, y #) -> (x, y)
463
464 {-# INLINE quotRem #-}
465 _ `quotRem` 0 = divZeroError
466 n `quotRem` d = case n `quotRemInteger` d of
467 (# q, r #) -> (q, r)
468
469 #if defined(MIN_VERSION_integer_gmp)
470 -- | @since 4.8.0.0
471 instance Integral Natural where
472 toInteger = naturalToInteger
473
474 divMod = quotRemNatural
475 div = quotNatural
476 mod = remNatural
477
478 quotRem = quotRemNatural
479 quot = quotNatural
480 rem = remNatural
481 #else
482 -- | @since 4.8.0.0
483 instance Integral Natural where
484 {-# INLINE toInteger #-}
485 toInteger (Natural a) = a
486
487 {-# INLINE quot #-}
488 Natural a `quot` Natural b = Natural (a `quot` b)
489
490 {-# INLINE rem #-}
491 Natural a `rem` Natural b = Natural (a `rem` b)
492
493 {-# INLINE div #-}
494 Natural a `div` Natural b = Natural (a `div` b)
495
496 {-# INLINE mod #-}
497 Natural a `mod` Natural b = Natural (a `mod` b)
498
499 {-# INLINE divMod #-}
500 Natural a `divMod` Natural b = coerce (a `divMod` b)
501
502 {-# INLINE quotRem #-}
503 Natural a `quotRem` Natural b = coerce (a `quotRem` b)
504 #endif
505
506 --------------------------------------------------------------
507 -- Instances for @Ratio@
508 --------------------------------------------------------------
509
510 -- | @since 2.0.1
511 instance (Integral a) => Ord (Ratio a) where
512 {-# SPECIALIZE instance Ord Rational #-}
513 (x:%y) <= (x':%y') = x * y' <= x' * y
514 (x:%y) < (x':%y') = x * y' < x' * y
515
516 -- | @since 2.0.1
517 instance (Integral a) => Num (Ratio a) where
518 {-# SPECIALIZE instance Num Rational #-}
519 (x:%y) + (x':%y') = reduce (x*y' + x'*y) (y*y')
520 (x:%y) - (x':%y') = reduce (x*y' - x'*y) (y*y')
521 (x:%y) * (x':%y') = reduce (x * x') (y * y')
522 negate (x:%y) = (-x) :% y
523 abs (x:%y) = abs x :% y
524 signum (x:%_) = signum x :% 1
525 fromInteger x = fromInteger x :% 1
526
527 -- | @since 2.0.1
528 {-# RULES "fromRational/id" fromRational = id :: Rational -> Rational #-}
529 instance (Integral a) => Fractional (Ratio a) where
530 {-# SPECIALIZE instance Fractional Rational #-}
531 (x:%y) / (x':%y') = (x*y') % (y*x')
532 recip (0:%_) = ratioZeroDenominatorError
533 recip (x:%y)
534 | x < 0 = negate y :% negate x
535 | otherwise = y :% x
536 fromRational (x:%y) = fromInteger x % fromInteger y
537
538 -- | @since 2.0.1
539 instance (Integral a) => Real (Ratio a) where
540 {-# SPECIALIZE instance Real Rational #-}
541 toRational (x:%y) = toInteger x :% toInteger y
542
543 -- | @since 2.0.1
544 instance (Integral a) => RealFrac (Ratio a) where
545 {-# SPECIALIZE instance RealFrac Rational #-}
546 properFraction (x:%y) = (fromInteger (toInteger q), r:%y)
547 where (q,r) = quotRem x y
548
549 -- | @since 2.0.1
550 instance (Show a) => Show (Ratio a) where
551 {-# SPECIALIZE instance Show Rational #-}
552 showsPrec p (x:%y) = showParen (p > ratioPrec) $
553 showsPrec ratioPrec1 x .
554 showString " % " .
555 -- H98 report has spaces round the %
556 -- but we removed them [May 04]
557 -- and added them again for consistency with
558 -- Haskell 98 [Sep 08, #1920]
559 showsPrec ratioPrec1 y
560
561 -- | @since 2.0.1
562 instance (Integral a) => Enum (Ratio a) where
563 {-# SPECIALIZE instance Enum Rational #-}
564 succ x = x + 1
565 pred x = x - 1
566
567 toEnum n = fromIntegral n :% 1
568 fromEnum = fromInteger . truncate
569
570 enumFrom = numericEnumFrom
571 enumFromThen = numericEnumFromThen
572 enumFromTo = numericEnumFromTo
573 enumFromThenTo = numericEnumFromThenTo
574
575 --------------------------------------------------------------
576 -- Coercions
577 --------------------------------------------------------------
578
579 -- | general coercion from integral types
580 {-# NOINLINE [1] fromIntegral #-}
581 fromIntegral :: (Integral a, Num b) => a -> b
582 fromIntegral = fromInteger . toInteger
583
584 {-# RULES
585 "fromIntegral/Int->Int" fromIntegral = id :: Int -> Int
586 #-}
587
588 {-# RULES
589 "fromIntegral/Int->Word" fromIntegral = \(I# x#) -> W# (int2Word# x#)
590 "fromIntegral/Word->Int" fromIntegral = \(W# x#) -> I# (word2Int# x#)
591 "fromIntegral/Word->Word" fromIntegral = id :: Word -> Word
592 #-}
593
594 {-# RULES
595 "fromIntegral/Natural->Natural" fromIntegral = id :: Natural -> Natural
596 "fromIntegral/Natural->Integer" fromIntegral = toInteger :: Natural->Integer
597 "fromIntegral/Natural->Word" fromIntegral = naturalToWord
598 #-}
599
600 {-# RULES
601 "fromIntegral/Word->Natural" fromIntegral = wordToNatural
602 "fromIntegral/Int->Natural" fromIntegral = intToNatural
603 #-}
604
605 -- | general coercion to fractional types
606 realToFrac :: (Real a, Fractional b) => a -> b
607 {-# NOINLINE [1] realToFrac #-}
608 realToFrac = fromRational . toRational
609
610 --------------------------------------------------------------
611 -- Overloaded numeric functions
612 --------------------------------------------------------------
613
614 -- | Converts a possibly-negative 'Real' value to a string.
615 showSigned :: (Real a)
616 => (a -> ShowS) -- ^ a function that can show unsigned values
617 -> Int -- ^ the precedence of the enclosing context
618 -> a -- ^ the value to show
619 -> ShowS
620 showSigned showPos p x
621 | x < 0 = showParen (p > 6) (showChar '-' . showPos (-x))
622 | otherwise = showPos x
623
624 even, odd :: (Integral a) => a -> Bool
625 even n = n `rem` 2 == 0
626 odd = not . even
627 {-# INLINABLE even #-}
628 {-# INLINABLE odd #-}
629
630 -------------------------------------------------------
631 -- | raise a number to a non-negative integral power
632 {-# SPECIALISE [1] (^) ::
633 Integer -> Integer -> Integer,
634 Integer -> Int -> Integer,
635 Int -> Int -> Int #-}
636 {-# INLINABLE [1] (^) #-} -- See Note [Inlining (^)]
637 (^) :: (Num a, Integral b) => a -> b -> a
638 x0 ^ y0 | y0 < 0 = errorWithoutStackTrace "Negative exponent"
639 | y0 == 0 = 1
640 | otherwise = f x0 y0
641 where -- f : x0 ^ y0 = x ^ y
642 f x y | even y = f (x * x) (y `quot` 2)
643 | y == 1 = x
644 | otherwise = g (x * x) (y `quot` 2) x -- See Note [Half of y - 1]
645 -- g : x0 ^ y0 = (x ^ y) * z
646 g x y z | even y = g (x * x) (y `quot` 2) z
647 | y == 1 = x * z
648 | otherwise = g (x * x) (y `quot` 2) (x * z) -- See Note [Half of y - 1]
649
650 -- | raise a number to an integral power
651 (^^) :: (Fractional a, Integral b) => a -> b -> a
652 {-# INLINABLE [1] (^^) #-} -- See Note [Inlining (^)
653 x ^^ n = if n >= 0 then x^n else recip (x^(negate n))
654
655 {- Note [Half of y - 1]
656 ~~~~~~~~~~~~~~~~~~~~~
657 Since y is guaranteed to be odd and positive here,
658 half of y - 1 can be computed as y `quot` 2, optimising subtraction away.
659 -}
660
661 {- Note [Inlining (^)
662 ~~~~~~~~~~~~~~~~~~~~~
663 The INLINABLE pragma allows (^) to be specialised at its call sites.
664 If it is called repeatedly at the same type, that can make a huge
665 difference, because of those constants which can be repeatedly
666 calculated.
667
668 Currently the fromInteger calls are not floated because we get
669 \d1 d2 x y -> blah
670 after the gentle round of simplification. -}
671
672 {- Rules for powers with known small exponent
673 see #5237
674 For small exponents, (^) is inefficient compared to manually
675 expanding the multiplication tree.
676 Here, rules for the most common exponent types are given.
677 The range of exponents for which rules are given is quite
678 arbitrary and kept small to not unduly increase the number of rules.
679 0 and 1 are excluded based on the assumption that nobody would
680 write x^0 or x^1 in code and the cases where an exponent could
681 be statically resolved to 0 or 1 are rare.
682
683 It might be desirable to have corresponding rules also for
684 exponents of other types (e. g., Word), but it's doubtful they
685 would fire, since the exponents of other types tend to get
686 floated out before the rule has a chance to fire.
687
688 Also desirable would be rules for (^^), but I haven't managed
689 to get those to fire.
690
691 Note: Trying to save multiplications by sharing the square for
692 exponents 4 and 5 does not save time, indeed, for Double, it is
693 up to twice slower, so the rules contain flat sequences of
694 multiplications.
695 -}
696
697 {-# RULES
698 "^2/Int" forall x. x ^ (2 :: Int) = let u = x in u*u
699 "^3/Int" forall x. x ^ (3 :: Int) = let u = x in u*u*u
700 "^4/Int" forall x. x ^ (4 :: Int) = let u = x in u*u*u*u
701 "^5/Int" forall x. x ^ (5 :: Int) = let u = x in u*u*u*u*u
702 "^2/Integer" forall x. x ^ (2 :: Integer) = let u = x in u*u
703 "^3/Integer" forall x. x ^ (3 :: Integer) = let u = x in u*u*u
704 "^4/Integer" forall x. x ^ (4 :: Integer) = let u = x in u*u*u*u
705 "^5/Integer" forall x. x ^ (5 :: Integer) = let u = x in u*u*u*u*u
706 #-}
707
708 -------------------------------------------------------
709 -- Special power functions for Rational
710 --
711 -- see #4337
712 --
713 -- Rationale:
714 -- For a legitimate Rational (n :% d), the numerator and denominator are
715 -- coprime, i.e. they have no common prime factor.
716 -- Therefore all powers (n ^ a) and (d ^ b) are also coprime, so it is
717 -- not necessary to compute the greatest common divisor, which would be
718 -- done in the default implementation at each multiplication step.
719 -- Since exponentiation quickly leads to very large numbers and
720 -- calculation of gcds is generally very slow for large numbers,
721 -- avoiding the gcd leads to an order of magnitude speedup relatively
722 -- soon (and an asymptotic improvement overall).
723 --
724 -- Note:
725 -- We cannot use these functions for general Ratio a because that would
726 -- change results in a multitude of cases.
727 -- The cause is that if a and b are coprime, their remainders by any
728 -- positive modulus generally aren't, so in the default implementation
729 -- reduction occurs.
730 --
731 -- Example:
732 -- (17 % 3) ^ 3 :: Ratio Word8
733 -- Default:
734 -- (17 % 3) ^ 3 = ((17 % 3) ^ 2) * (17 % 3)
735 -- = ((289 `mod` 256) % 9) * (17 % 3)
736 -- = (33 % 9) * (17 % 3)
737 -- = (11 % 3) * (17 % 3)
738 -- = (187 % 9)
739 -- But:
740 -- ((17^3) `mod` 256) % (3^3) = (4913 `mod` 256) % 27
741 -- = 49 % 27
742 --
743 -- TODO:
744 -- Find out whether special-casing for numerator, denominator or
745 -- exponent = 1 (or -1, where that may apply) gains something.
746
747 -- Special version of (^) for Rational base
748 {-# RULES "(^)/Rational" (^) = (^%^) #-}
749 (^%^) :: Integral a => Rational -> a -> Rational
750 (n :% d) ^%^ e
751 | e < 0 = errorWithoutStackTrace "Negative exponent"
752 | e == 0 = 1 :% 1
753 | otherwise = (n ^ e) :% (d ^ e)
754
755 -- Special version of (^^) for Rational base
756 {-# RULES "(^^)/Rational" (^^) = (^^%^^) #-}
757 (^^%^^) :: Integral a => Rational -> a -> Rational
758 (n :% d) ^^%^^ e
759 | e > 0 = (n ^ e) :% (d ^ e)
760 | e == 0 = 1 :% 1
761 | n > 0 = (d ^ (negate e)) :% (n ^ (negate e))
762 | n == 0 = ratioZeroDenominatorError
763 | otherwise = let nn = d ^ (negate e)
764 dd = (negate n) ^ (negate e)
765 in if even e then (nn :% dd) else (negate nn :% dd)
766
767 -------------------------------------------------------
768 -- | @'gcd' x y@ is the non-negative factor of both @x@ and @y@ of which
769 -- every common factor of @x@ and @y@ is also a factor; for example
770 -- @'gcd' 4 2 = 2@, @'gcd' (-4) 6 = 2@, @'gcd' 0 4@ = @4@. @'gcd' 0 0@ = @0@.
771 -- (That is, the common divisor that is \"greatest\" in the divisibility
772 -- preordering.)
773 --
774 -- Note: Since for signed fixed-width integer types, @'abs' 'minBound' < 0@,
775 -- the result may be negative if one of the arguments is @'minBound'@ (and
776 -- necessarily is if the other is @0@ or @'minBound'@) for such types.
777 gcd :: (Integral a) => a -> a -> a
778 {-# NOINLINE [1] gcd #-}
779 gcd x y = gcd' (abs x) (abs y)
780 where gcd' a 0 = a
781 gcd' a b = gcd' b (a `rem` b)
782
783 -- | @'lcm' x y@ is the smallest positive integer that both @x@ and @y@ divide.
784 lcm :: (Integral a) => a -> a -> a
785 {-# SPECIALISE lcm :: Int -> Int -> Int #-}
786 {-# SPECIALISE lcm :: Word -> Word -> Word #-}
787 {-# NOINLINE [1] lcm #-}
788 lcm _ 0 = 0
789 lcm 0 _ = 0
790 lcm x y = abs ((x `quot` (gcd x y)) * y)
791
792 #if defined(OPTIMISE_INTEGER_GCD_LCM)
793 {-# RULES
794 "gcd/Int->Int->Int" gcd = gcdInt'
795 "gcd/Integer->Integer->Integer" gcd = gcdInteger
796 "lcm/Integer->Integer->Integer" lcm = lcmInteger
797 "gcd/Natural->Natural->Natural" gcd = gcdNatural
798 "lcm/Natural->Natural->Natural" lcm = lcmNatural
799 #-}
800
801 gcdInt' :: Int -> Int -> Int
802 gcdInt' (I# x) (I# y) = I# (gcdInt x y)
803
804 {-# RULES
805 "gcd/Word->Word->Word" gcd = gcdWord'
806 #-}
807
808 gcdWord' :: Word -> Word -> Word
809 gcdWord' (W# x) (W# y) = W# (gcdWord x y)
810 #endif
811
812 integralEnumFrom :: (Integral a, Bounded a) => a -> [a]
813 integralEnumFrom n = map fromInteger [toInteger n .. toInteger (maxBound `asTypeOf` n)]
814
815 integralEnumFromThen :: (Integral a, Bounded a) => a -> a -> [a]
816 integralEnumFromThen n1 n2
817 | i_n2 >= i_n1 = map fromInteger [i_n1, i_n2 .. toInteger (maxBound `asTypeOf` n1)]
818 | otherwise = map fromInteger [i_n1, i_n2 .. toInteger (minBound `asTypeOf` n1)]
819 where
820 i_n1 = toInteger n1
821 i_n2 = toInteger n2
822
823 integralEnumFromTo :: Integral a => a -> a -> [a]
824 integralEnumFromTo n m = map fromInteger [toInteger n .. toInteger m]
825
826 integralEnumFromThenTo :: Integral a => a -> a -> a -> [a]
827 integralEnumFromThenTo n1 n2 m
828 = map fromInteger [toInteger n1, toInteger n2 .. toInteger m]