Adjust behaviour of gcd
[ghc-base.git] / GHC / Real.lhs
index f7e2eb2..51d7db4 100644 (file)
@@ -1,12 +1,12 @@
 \begin{code}
-{-# OPTIONS_GHC -fno-implicit-prelude #-}
+{-# LANGUAGE CPP, NoImplicitPrelude, MagicHash, UnboxedTuples #-}
 {-# OPTIONS_HADDOCK hide #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  GHC.Real
--- Copyright   :  (c) The FFI Task Force, 1994-2002
+-- Copyright   :  (c) The University of Glasgow, 1994-2002
 -- License     :  see libraries/base/LICENSE
--- 
+--
 -- Maintainer  :  cvs-ghc@haskell.org
 -- Stability   :  internal
 -- Portability :  non-portable (GHC Extensions)
@@ -24,12 +24,13 @@ import GHC.Num
 import GHC.List
 import GHC.Enum
 import GHC.Show
+import GHC.Err
 
 infixr 8  ^, ^^
 infixl 7  /, `quot`, `rem`, `div`, `mod`
 infixl 7  %
 
-default ()              -- Double isn't available yet, 
+default ()              -- Double isn't available yet,
                         -- and we shouldn't be using defaults anyway
 \end{code}
 
@@ -42,7 +43,7 @@ default ()              -- Double isn't available yet,
 
 \begin{code}
 -- | Rational numbers, with numerator and denominator of some 'Integral' type.
-data  (Integral a)      => Ratio a = !a :% !a  deriving (Eq)
+data  Ratio a = !a :% !a  deriving (Eq)
 
 -- | Arbitrary-precision rational numbers, represented as a ratio of
 -- two 'Integer' values.  A rational number may be constructed using
@@ -57,8 +58,8 @@ infinity, notANumber :: Rational
 infinity   = 1 :% 0
 notANumber = 0 :% 0
 
--- Use :%, not % for Inf/NaN; the latter would 
--- immediately lead to a runtime error, because it normalises. 
+-- Use :%, not % for Inf/NaN; the latter would
+-- immediately lead to a runtime error, because it normalises.
 \end{code}
 
 
@@ -132,10 +133,15 @@ class  (Real a, Enum a) => Integral a  where
     -- | conversion to 'Integer'
     toInteger           :: a -> Integer
 
+    {-# INLINE quot #-}
+    {-# INLINE rem #-}
+    {-# INLINE div #-}
+    {-# INLINE mod #-}
     n `quot` d          =  q  where (q,_) = quotRem n d
     n `rem` d           =  r  where (_,r) = quotRem n d
     n `div` d           =  q  where (q,_) = divMod n d
     n `mod` d           =  r  where (_,r) = divMod n d
+
     divMod n d          =  if signum r == negate (signum d) then (q-1, r+d) else qr
                            where qr@(q,r) = quotRem n d
 
@@ -153,6 +159,8 @@ class  (Num a) => Fractional a  where
     -- @('Fractional' a) => a@.
     fromRational        :: Rational -> a
 
+    {-# INLINE recip #-}
+    {-# INLINE (/) #-}
     recip x             =  1 / x
     x / y               = x * recip y
 
@@ -173,25 +181,28 @@ class  (Real a, Fractional a) => RealFrac a  where
     properFraction      :: (Integral b) => a -> (b,a)
     -- | @'truncate' x@ returns the integer nearest @x@ between zero and @x@
     truncate            :: (Integral b) => a -> b
-    -- | @'round' x@ returns the nearest integer to @x@
+    -- | @'round' x@ returns the nearest integer to @x@;
+    --   the even integer if @x@ is equidistant between two integers
     round               :: (Integral b) => a -> b
     -- | @'ceiling' x@ returns the least integer not less than @x@
     ceiling             :: (Integral b) => a -> b
     -- | @'floor' x@ returns the greatest integer not greater than @x@
     floor               :: (Integral b) => a -> b
 
+    {-# INLINE truncate #-}
     truncate x          =  m  where (m,_) = properFraction x
-    
+
     round x             =  let (n,r) = properFraction x
                                m     = if r < 0 then n - 1 else n + 1
                            in case signum (abs r - 0.5) of
                                 -1 -> n
                                 0  -> if even n then n else m
                                 1  -> m
-    
+                                _  -> error "round default defn: Bad value"
+
     ceiling x           =  if r > 0 then n + 1 else n
                            where (n,r) = properFraction x
-    
+
     floor x             =  if r < 0 then n - 1 else n
                            where (n,r) = properFraction x
 \end{code}
@@ -201,20 +212,21 @@ These 'numeric' enumerations come straight from the Report
 
 \begin{code}
 numericEnumFrom         :: (Fractional a) => a -> [a]
-numericEnumFrom         =  iterate (+1)
+numericEnumFrom n      =  n `seq` (n : numericEnumFrom (n + 1))
 
 numericEnumFromThen     :: (Fractional a) => a -> a -> [a]
-numericEnumFromThen n m =  iterate (+(m-n)) n
+numericEnumFromThen n m        = n `seq` m `seq` (n : numericEnumFromThen m (m+m-n))
 
 numericEnumFromTo       :: (Ord a, Fractional a) => a -> a -> [a]
 numericEnumFromTo n m   = takeWhile (<= m + 1/2) (numericEnumFrom n)
 
 numericEnumFromThenTo   :: (Ord a, Fractional a) => a -> a -> a -> [a]
-numericEnumFromThenTo e1 e2 e3 = takeWhile pred (numericEnumFromThen e1 e2)
+numericEnumFromThenTo e1 e2 e3
+    = takeWhile predicate (numericEnumFromThen e1 e2)
                                 where
                                  mid = (e2 - e1) / 2
-                                 pred | e2 >= e1  = (<= e3 + mid)
-                                      | otherwise = (>= e3 + mid)
+                                 predicate | e2 >= e1  = (<= e3 + mid)
+                                           | otherwise = (>= e3 + mid)
 \end{code}
 
 
@@ -229,36 +241,42 @@ instance  Real Int  where
     toRational x        =  toInteger x % 1
 
 instance  Integral Int  where
-    toInteger i = int2Integer i  -- give back a full-blown Integer
+    toInteger (I# i) = smallInteger i
 
     a `quot` b
      | b == 0                     = divZeroError
-     | a == minBound && b == (-1) = overflowError
+     | b == (-1) && a == minBound = overflowError -- Note [Order of tests]
+                                                  -- in GHC.Int
      | otherwise                  =  a `quotInt` b
 
     a `rem` b
      | b == 0                     = divZeroError
-     | a == minBound && b == (-1) = overflowError
+     | b == (-1) && a == minBound = overflowError -- Note [Order of tests]
+                                                  -- in GHC.Int
      | otherwise                  =  a `remInt` b
 
     a `div` b
      | b == 0                     = divZeroError
-     | a == minBound && b == (-1) = overflowError
+     | b == (-1) && a == minBound = overflowError -- Note [Order of tests]
+                                                  -- in GHC.Int
      | otherwise                  =  a `divInt` b
 
     a `mod` b
      | b == 0                     = divZeroError
-     | a == minBound && b == (-1) = overflowError
+     | b == (-1) && a == minBound = overflowError -- Note [Order of tests]
+                                                  -- in GHC.Int
      | otherwise                  =  a `modInt` b
 
     a `quotRem` b
      | b == 0                     = divZeroError
-     | a == minBound && b == (-1) = overflowError
+     | b == (-1) && a == minBound = overflowError -- Note [Order of tests]
+                                                  -- in GHC.Int
      | otherwise                  =  a `quotRemInt` b
 
     a `divMod` b
      | b == 0                     = divZeroError
-     | a == minBound && b == (-1) = overflowError
+     | b == (-1) && a == minBound = overflowError -- Note [Order of tests]
+                                                  -- in GHC.Int
      | otherwise                  =  a `divModInt` b
 \end{code}
 
@@ -276,17 +294,19 @@ instance  Real Integer  where
 instance  Integral Integer where
     toInteger n      = n
 
-    a `quot` 0 = divZeroError
+    _ `quot` 0 = divZeroError
     n `quot` d = n `quotInteger` d
 
-    a `rem` 0 = divZeroError
+    _ `rem` 0 = divZeroError
     n `rem`  d = n `remInteger`  d
 
-    a `divMod` 0 = divZeroError
-    a `divMod` b = a `divModInteger` b
+    _ `divMod` 0 = divZeroError
+    a `divMod` b = case a `divModInteger` b of
+                   (# x, y #) -> (x, y)
 
-    a `quotRem` 0 = divZeroError
-    a `quotRem` b = a `quotRemInteger` b
+    _ `quotRem` 0 = divZeroError
+    a `quotRem` b = case a `quotRemInteger` b of
+                    (# q, r #) -> (q, r)
 
     -- use the defaults for div & mod
 \end{code}
@@ -314,11 +334,15 @@ instance  (Integral a)  => Num (Ratio a)  where
     signum (x:%_)       =  signum x :% 1
     fromInteger x       =  fromInteger x :% 1
 
+{-# RULES "fromRational/id" fromRational = id :: Rational -> Rational #-}
 instance  (Integral a)  => Fractional (Ratio a)  where
     {-# SPECIALIZE instance Fractional Rational #-}
     (x:%y) / (x':%y')   =  (x*y') % (y*x')
-    recip (x:%y)        =  y % x
-    fromRational (x:%y) =  fromInteger x :% fromInteger y
+    recip (0:%_)        = error "Ratio.%: zero denominator"
+    recip (x:%y)
+        | x < 0         = negate y :% negate x
+        | otherwise     = y :% x
+    fromRational (x:%y) =  fromInteger x % fromInteger y
 
 instance  (Integral a)  => Real (Ratio a)  where
     {-# SPECIALIZE instance Real Rational #-}
@@ -332,9 +356,12 @@ instance  (Integral a)  => RealFrac (Ratio a)  where
 instance  (Integral a)  => Show (Ratio a)  where
     {-# SPECIALIZE instance Show Rational #-}
     showsPrec p (x:%y)  =  showParen (p > ratioPrec) $
-                           showsPrec ratioPrec1 x . 
-                           showString "%" .     -- H98 report has spaces round the %
-                                                -- but we removed them [May 04]
+                           showsPrec ratioPrec1 x .
+                           showString " % " .
+                           -- H98 report has spaces round the %
+                           -- but we removed them [May 04]
+                           -- and added them again for consistency with
+                           -- Haskell 98 [Sep 08, #1920]
                            showsPrec ratioPrec1 y
 
 instance  (Integral a)  => Enum (Ratio a)  where
@@ -342,7 +369,7 @@ instance  (Integral a)  => Enum (Ratio a)  where
     succ x              =  x + 1
     pred x              =  x - 1
 
-    toEnum n            =  fromInteger (int2Integer n) :% 1
+    toEnum n            =  fromIntegral n :% 1
     fromEnum            =  fromInteger . truncate
 
     enumFrom            =  numericEnumFrom
@@ -389,7 +416,7 @@ showSigned :: (Real a)
   -> Int                -- ^ the precedence of the enclosing context
   -> a                  -- ^ the value to show
   -> ShowS
-showSigned showPos p x 
+showSigned showPos p x
    | x < 0     = showParen (p > 6) (showChar '-' . showPos (-x))
    | otherwise = showPos x
 
@@ -403,28 +430,100 @@ odd             =  not . even
         Integer -> Integer -> Integer,
         Integer -> Int -> Integer,
         Int -> Int -> Int #-}
-(^)             :: (Num a, Integral b) => a -> b -> a
-_ ^ 0           =  1
-x ^ n | n > 0   =  f x (n-1) x
-                   where f _ 0 y = y
-                         f a d y = g a d  where
-                                   g b i | even i  = g (b*b) (i `quot` 2)
-                                         | otherwise = f b (i-1) (b*y)
-_ ^ _           = error "Prelude.^: negative exponent"
+{-# INLINABLE (^) #-}    -- See Note [Inlining (^)]
+(^) :: (Num a, Integral b) => a -> b -> a
+x0 ^ y0 | y0 < 0    = error "Negative exponent"
+        | y0 == 0   = 1
+        | otherwise = f x0 y0
+    where -- f : x0 ^ y0 = x ^ y
+          f x y | even y    = f (x * x) (y `quot` 2)
+                | y == 1    = x
+                | otherwise = g (x * x) ((y - 1) `quot` 2) x
+          -- g : x0 ^ y0 = (x ^ y) * z
+          g x y z | even y = g (x * x) (y `quot` 2) z
+                  | y == 1 = x * z
+                  | otherwise = g (x * x) ((y - 1) `quot` 2) (x * z)
 
 -- | raise a number to an integral power
-{-# SPECIALISE (^^) ::
-        Rational -> Int -> Rational #-}
 (^^)            :: (Fractional a, Integral b) => a -> b -> a
+{-# INLINABLE (^^) #-}         -- See Note [Inlining (^)
 x ^^ n          =  if n >= 0 then x^n else recip (x^(negate n))
 
+{- Note [Inlining (^)
+   ~~~~~~~~~~~~~~~~~~~~~
+   The INLINABLE pragma allows (^) to be specialised at its call sites.
+   If it is called repeatedly at the same type, that can make a huge
+   difference, because of those constants which can be repeatedly
+   calculated.
+
+   Currently the fromInteger calls are not floated because we get
+             \d1 d2 x y -> blah
+   after the gentle round of simplification. -}
 
 -------------------------------------------------------
--- | @'gcd' x y@ is the greatest (positive) integer that divides both @x@
+-- Special power functions for Rational
+--
+-- see #4337
+--
+-- Rationale:
+-- For a legitimate Rational (n :% d), the numerator and denominator are
+-- coprime, i.e. they have no common prime factor.
+-- Therefore all powers (n ^ a) and (d ^ b) are also coprime, so it is
+-- not necessary to compute the greatest common divisor, which would be
+-- done in the default implementation at each multiplication step.
+-- Since exponentiation quickly leads to very large numbers and
+-- calculation of gcds is generally very slow for large numbers,
+-- avoiding the gcd leads to an order of magnitude speedup relatively
+-- soon (and an asymptotic improvement overall).
+--
+-- Note:
+-- We cannot use these functions for general Ratio a because that would
+-- change results in a multitude of cases.
+-- The cause is that if a and b are coprime, their remainders by any
+-- positive modulus generally aren't, so in the default implementation
+-- reduction occurs.
+--
+-- Example:
+-- (17 % 3) ^ 3 :: Ratio Word8
+-- Default:
+-- (17 % 3) ^ 3 = ((17 % 3) ^ 2) * (17 % 3)
+--              = ((289 `mod` 256) % 9) * (17 % 3)
+--              = (33 % 9) * (17 % 3)
+--              = (11 % 3) * (17 % 3)
+--              = (187 % 9)
+-- But:
+-- ((17^3) `mod` 256) % (3^3)   = (4913 `mod` 256) % 27
+--                              = 49 % 27
+--
+-- TODO:
+-- Find out whether special-casing for numerator, denominator or
+-- exponent = 1 (or -1, where that may apply) gains something.
+
+-- Special version of (^) for Rational base
+{-# RULES "(^)/Rational"    (^) = (^%^) #-}
+(^%^)           :: Integral a => Rational -> a -> Rational
+(n :% d) ^%^ e
+    | e < 0     = error "Negative exponent"
+    | e == 0    = 1 :% 1
+    | otherwise = (n ^ e) :% (d ^ e)
+
+-- Special version of (^^) for Rational base
+{-# RULES "(^^)/Rational"   (^^) = (^^%^^) #-}
+(^^%^^)         :: Integral a => Rational -> a -> Rational
+(n :% d) ^^%^^ e
+    | e > 0     = (n ^ e) :% (d ^ e)
+    | e == 0    = 1 :% 1
+    | n > 0     = (d ^ (negate e)) :% (n ^ (negate e))
+    | n == 0    = error "Ratio.%: zero denominator"
+    | otherwise = let nn = d ^ (negate e)
+                      dd = (negate n) ^ (negate e)
+                  in if even e then (nn :% dd) else (negate nn :% dd)
+
+-------------------------------------------------------
+-- | @'gcd' x y@ is the greatest (nonnegative) integer that divides both @x@
 -- and @y@; for example @'gcd' (-3) 6@ = @3@, @'gcd' (-3) (-6)@ = @3@,
--- @'gcd' 0 4@ = @4@.  @'gcd' 0 0@ raises a runtime error.
+-- @'gcd' 0 4@ = @4@.  @'gcd' 0 0@ = @0@.
 gcd             :: (Integral a) => a -> a -> a
-gcd 0 0         =  error "Prelude.gcd: gcd 0 0 is undefined"
 gcd x y         =  gcd' (abs x) (abs y)
                    where gcd' a 0  =  a
                          gcd' a b  =  gcd' b (a `rem` b)
@@ -436,13 +535,17 @@ lcm _ 0         =  0
 lcm 0 _         =  0
 lcm x y         =  abs ((x `quot` (gcd x y)) * y)
 
-
+#ifdef OPTIMISE_INTEGER_GCD_LCM
 {-# RULES
 "gcd/Int->Int->Int"             gcd = gcdInt
 "gcd/Integer->Integer->Integer" gcd = gcdInteger
 "lcm/Integer->Integer->Integer" lcm = lcmInteger
  #-}
 
+gcdInt :: Int -> Int -> Int
+gcdInt a b = fromIntegral (gcdInteger (fromIntegral a) (fromIntegral b))
+#endif
+
 integralEnumFrom :: (Integral a, Bounded a) => a -> [a]
 integralEnumFrom n = map fromInteger [toInteger n .. toInteger (maxBound `asTypeOf` n)]