{-| Stability: experimental -}
module Data.SciRatio
(
-- * The 'SciRatio' type
SciRatio
, SciRational
, (.^)
, fracSignificand
, base10Exponent
-- * Specialized functions
, (^!)
, (^^!)
, fromSciRatio
-- * Miscellaneous utilities
, factorizeBase
, ilogBase
-- * Deprecated
, intLog
) where
import Data.Ratio ((%), denominator, numerator)
import Data.Hashable (Hashable(hashWithSalt))
infixr 8 ^!, ^^!
infixl 7 :^, .^
infixl 1 ~~
-- Implementation note: a SciRatio must *always* remain in canonical form, so
-- don't use :^ directly unless you're absolutely sure it will stay canonical
-- (use .^ instead if you're unsure).
-- | Represents a fractional number stored in scientific notation: a product
-- of a fractional significand and an integral power of 10.
--
-- - The significand has type @a@ and should be both @'Fractional'@ and
-- @'Real'@. Although this could be a floating-point type, it is not
-- recommended as floating-points use inexact arithmetic and strange bugs
-- will occur as a result.
--
-- - The exponent has type @b@ and should be @'Integral'@.
--
-- @'SciRatio'@ behaves in the same way as an ordinary @'Data.Ratio.Ratio'@
-- and supports the same operations. The main property is that it is more
-- efficient than @'Data.Ratio.Ratio'@ when the exponent is large:
--
-- >>> 5 .^ 99999999 :: SciRational -- works fine
-- >>> 5e99999999 :: Rational -- takes forever
--
-- Specialized functions are provided in cases where they can be implemented
-- more efficiently than the default implementation.
--
-- The number is always stored in a unique, canonical form: the significand
-- shall never contain factors of 2 and 5 simultaneously, and their
-- multiplicities shall always be nonnegative. (The significand is treated
-- as a rational number factorized into a product of prime numbers with
-- integral exponents that are not necessarily positive.)
--
-- __Note__: If inputs differ greatly in magnitude, @('+')@ and @('-')@ can
-- be quite slow: complexity is linear with the absolute
-- difference of the exponents. Furthermore, the complexity of
-- @'toRational'@ is linear with the magnitude of the exponent.
-- These also apply to any functions that indirectly use these
-- operations (including all operations in @'RealFrac'@ and
-- @'Enum'@).
--
data SciRatio a b = !a :^ !b deriving Eq
-- | A specialization of 'SciRatio'.
type SciRational = SciRatio Rational Integer
instance (Fractional a, Real a, Integral b, Read a, Read b) =>
Read (SciRatio a b) where
{-# SPECIALIZE instance Read SciRational #-}
readsPrec p = readParen (p > prec) $ \ r -> do
(x, s) <- readsPrec (succ prec) r
(".^", t) <- lex s
(y, u) <- readsPrec (succ prec) t
return (x .^ y, u)
instance (Show a, Show b) => Show (SciRatio a b) where
{-# SPECIALIZE instance Show SciRational #-}
showsPrec p (x :^ a) = showParen (p > prec) $
showsPrec (succ prec) x .
showString " .^ " .
showsPrec (succ prec) a
instance (Real a, Integral b, Ord a) => Ord (SciRatio a b) where
{-# SPECIALIZE instance Ord SciRational #-}
compare (x :^ a) (y :^ b) = case compare sgnX sgnY of
EQ -> case sgnX of
EQ -> EQ
LT -> invert order
GT -> order
k -> k
where x' = toRational x
y' = toRational y
sgnX = compare x 0
sgnY = compare y 0
absM = abs (numerator x' * denominator y')
absN = abs (numerator y' * denominator x')
invert GT = LT
invert EQ = EQ
invert LT = GT
order = case (_ilogBase 10 absM, _ilogBase 10 absN) of
-- initial comparison via floored logarithms to handle easy cases
-- where exponents are wildly different
(logM, logN) -> case compare (a + logM) (b + logN) of
-- compare directly if the exponents are too close
EQ | a >= b -> compare (absM * 10 ^ (a - b)) absN
| otherwise -> compare absM (absN * 10 ^ (b - a))
k -> k
instance (Hashable a, Hashable b) => Hashable (SciRatio a b) where
{-# SPECIALIZE instance Hashable SciRational #-}
hashWithSalt s (r :^ e) = hashWithSalt s (magic, r, e)
where magic = 0xaaa80d6b :: Int
instance (Fractional a, Real a, Integral b) => Num (SciRatio a b) where
{-# SPECIALIZE instance Num SciRational #-}
p + q = (x + y) .^ c where (x, y, c) = p ~~ q
p - q = (x - y) .^ c where (x, y, c) = p ~~ q
x :^ a * (y :^ b) = x * y .^ (a + b)
abs (y :^ b) = abs y :^ b
negate (y :^ b) = negate y :^ b
signum = (:^ 0) . signum . fracSignificand
fromInteger = (.^ 0) . fromInteger
instance (Fractional a, Real a, Integral b) => Fractional (SciRatio a b) where
{-# SPECIALIZE instance Fractional SciRational #-}
x :^ a / (y :^ b) = x / y .^ (a - b)
recip (y :^ b) = recip y .^ (-b)
fromRational = (.^ 0) . fromRational
instance (Fractional a, Real a, Integral b) => Real (SciRatio a b) where
{-# SPECIALIZE instance Real SciRational #-}
toRational (y :^ b) = toRational y * 10 ^^ b
instance (Fractional a, Real a, Integral b) => RealFrac (SciRatio a b) where
{-# SPECIALIZE instance RealFrac SciRational #-}
properFraction q = (\ (n, p) -> (n, fromRational p))
. properFraction $ toRational q
instance (Fractional a, Real a, Integral b) => Enum (SciRatio a b) where
{-# SPECIALIZE instance Enum SciRational #-}
succ = fromRational . succ . toRational
pred = fromRational . pred . toRational
toEnum = fromRational . toEnum
fromEnum = fromEnum . toRational
enumFrom = fmap fromRational . enumFrom . toRational
enumFromThen x = fmap fromRational
. enumFromThen (toRational x)
. toRational
enumFromTo x = fmap fromRational
. enumFromTo (toRational x)
. toRational
enumFromThenTo x y = fmap fromRational
. enumFromThenTo (toRational x) (toRational y)
. toRational
-- | Precedence of the @('.^')@ operator.
prec :: Int
prec = 7
-- | Construct a number such that
-- @significand .^ exponent == significand * 10 ^^ exponent@.
{-# SPECIALISE (.^) :: Rational -> Integer -> SciRational #-}
(.^) :: (Fractional a, Real a, Integral b) =>
a -- ^ significand
-> b -- ^ exponent
-> SciRatio a b
x .^ y = canonicalize $ x :^ y
-- | Extract the fractional significand.
fracSignificand :: SciRatio a b -> a
fracSignificand (x :^ _) = x
-- | Extract the base-10 exponent.
base10Exponent :: SciRatio a b -> b
base10Exponent (_ :^ x) = x
-- | Convert into a 'Fractional' number.
--
-- This is similar to 'realToFrac' but much more efficient for large
-- exponents.
{-# RULES "realToFrac/fromSciRational"
forall (x :: SciRational) . realToFrac x = fromSciRatio x;
"realToFrac/SciRational"
realToFrac = id :: SciRational -> SciRational;
"fromSciRatio/SciRational"
fromSciRatio = id #-}
{-# NOINLINE [1] fromSciRatio #-}
fromSciRatio :: (Real a, Integral b, Fractional c) => SciRatio a b -> c
fromSciRatio (x :^ y) = realToFrac x * 10 ^^ y
-- | Specialized, more efficient version of @('^^')@.
{-# RULES "(^^)/SciRational" forall (x :: SciRational) . (^^) x = (^^!) x #-}
(^^!) :: (Fractional a, Real a, Integral b, Integral c) =>
SciRatio a b -> c -> SciRatio a b
(x :^ a) ^^! b = x ^^ b :^ (a * fromIntegral b)
-- | Specialized, more efficient version of @('^')@.
{-# RULES "(^)/SciRational" forall (x :: SciRational) . (^) x = (^!) x #-}
(^!) :: (Real a, Integral b, Integral c) =>
SciRatio a b -> c -> SciRatio a b
(x :^ a) ^! b = x ^ b :^ (a * fromIntegral b)
-- | Matches the exponents.
(~~) :: (Fractional a, Integral b) => SciRatio a b -> SciRatio a b -> (a, a, b)
x :^ a ~~ y :^ b = (x * 10 ^^ (a - c), y * 10 ^^ (b - c), c)
where c = if abs a <= abs b then a else b
-- | Factorize a nonzero integer into a significand and a power of the base
-- such that the exponent is maximized:
--
-- > inputInteger = significand * base ^ exponent
--
-- That is, the significand shall not divisible by the base. The base must
-- be greater than one.
factorizeBase :: (Integral a, Integral b) =>
a -- ^ base
-> a -- ^ input integer
-> (a, b) -- ^ significand and exponent
factorizeBase _ 0 = error ("Data.SciRatio.factorizeBase: " ++
"input integer must be nonzero")
factorizeBase b n | b < 2 = error ("Data.SciRatio.factorizeBase: " ++
"base must be greater than one")
| otherwise = _factorizeBase b n
-- | Same as @'factorizeBase'@ but doesn't validate the arguments, so it might
-- loop forever or return something nonsensical.
_factorizeBase :: (Integral a, Integral b) => a -> a -> (a, b)
_factorizeBase b n | r /= 0 = (n, 0)
| r' /= 0 = (n'', e2 + 1)
| otherwise = (n''', e2 + 2)
where (n', r) = n `quotRem` b
(n'', e) = _factorizeBase (b * b) n'
(n''', r') = n'' `quotRem` b
e2 = e * 2
-- | Alias of @'factorizeBase'@.
--
-- Note: Despite what the name suggests, the function does /not/ compute the
-- floored logarithm.
{-# DEPRECATED intLog "use @'factorizeBase'@ instead." #-}
intLog :: (Integral a, Integral b) => a -> a -> (a, b)
intLog = factorizeBase
-- | Calculate the floored logarithm of a positive integer. The base must be
-- greater than one.
ilogBase :: (Integral a, Integral b) =>
a -- ^ base
-> a -- ^ input integer
-> b -- ^ floor of the logarithm
ilogBase b | b < 2 = error ("Data.SciRatio.ilogBase: " ++
"base must be greater than one")
| otherwise = _ilogBase b
-- | Same as @'ilogBase'@ but doesn't validate the arguments, so it might loop
-- forever or return something nonsensical.
_ilogBase :: (Integral a, Integral b) => a -> a -> b
_ilogBase = (fst .) . ilogr
-- use the trick from section 14.4 of The Haskell 98 Language Report
where ilogr b n | n < b = (0, n)
| n' < b = (l2, n')
| otherwise = (1 + l2, n' `div` b)
where (l, n') = ilogr (b * b) n
l2 = l * 2
-- | Convert into canonical form by removing all factors of 2 and 5 from the
-- denominator and factoring out as many powers of the base as possible.
canonicalize :: (Fractional a, Real a, Integral b) =>
SciRatio a b -> SciRatio a b
canonicalize (0 :^ _) = 0 :^ 0
canonicalize (r :^ e) =
let r' = toRational r in
case (_factorizeBase 10 (numerator r'), _factorizeBase 2 (denominator r')) of
((n, ne), (d, e2)) -> case _factorizeBase 5 d of
(d', e5) -> case compare e2 e5 of
EQ -> fromRational (n % d') :^ (e + ne - e5)
LT -> fromRational ((n * 2 ^ (e5 - e2)) % d') :^ (e + ne - e5)
GT -> fromRational ((n * 5 ^ (e2 - e5)) % d') :^ (e + ne - e2)