|
Packit |
785658 |
-- |
|
|
Packit |
785658 |
-- Module: Math.NumberTheory.Logarithms
|
|
Packit |
785658 |
-- Copyright: (c) 2011 Daniel Fischer
|
|
Packit |
785658 |
-- Licence: MIT
|
|
Packit |
785658 |
-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>
|
|
Packit |
785658 |
-- Stability: Provisional
|
|
Packit |
785658 |
-- Portability: Non-portable (GHC extensions)
|
|
Packit |
785658 |
--
|
|
Packit |
785658 |
-- Integer Logarithms. For efficiency, the internal representation of 'Integer's
|
|
Packit |
785658 |
-- from integer-gmp is used.
|
|
Packit |
785658 |
--
|
|
Packit |
785658 |
{-# LANGUAGE CPP #-}
|
|
Packit |
785658 |
{-# LANGUAGE MagicHash #-}
|
|
Packit |
785658 |
module Math.NumberTheory.Logarithms
|
|
Packit |
785658 |
( -- * Integer logarithms with input checks
|
|
Packit |
785658 |
integerLogBase
|
|
Packit |
785658 |
, integerLog2
|
|
Packit |
785658 |
, integerLog10
|
|
Packit |
785658 |
|
|
Packit |
785658 |
, naturalLogBase
|
|
Packit |
785658 |
, naturalLog2
|
|
Packit |
785658 |
, naturalLog10
|
|
Packit |
785658 |
|
|
Packit |
785658 |
, intLog2
|
|
Packit |
785658 |
, wordLog2
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- * Integer logarithms without input checks
|
|
Packit |
785658 |
--
|
|
Packit |
785658 |
-- | These functions are total, however, don't rely on the values with out-of-domain arguments.
|
|
Packit |
785658 |
, integerLogBase'
|
|
Packit |
785658 |
, integerLog2'
|
|
Packit |
785658 |
, integerLog10'
|
|
Packit |
785658 |
|
|
Packit |
785658 |
, intLog2'
|
|
Packit |
785658 |
, wordLog2'
|
|
Packit |
785658 |
) where
|
|
Packit |
785658 |
|
|
Packit |
785658 |
import GHC.Exts
|
|
Packit |
785658 |
|
|
Packit |
785658 |
import Data.Bits
|
|
Packit |
785658 |
import Data.Array.Unboxed
|
|
Packit |
785658 |
import Numeric.Natural
|
|
Packit |
785658 |
|
|
Packit |
785658 |
import GHC.Integer.Logarithms.Compat
|
|
Packit |
785658 |
#if Base48 && defined(MIN_VERSION_integer_gmp)
|
|
Packit |
785658 |
import GHC.Integer.GMP.Internals (Integer (..))
|
|
Packit |
785658 |
import GHC.Natural
|
|
Packit |
785658 |
#endif
|
|
Packit |
785658 |
|
|
Packit |
785658 |
#if CheckBounds
|
|
Packit |
785658 |
import Data.Array.IArray (IArray, (!))
|
|
Packit |
785658 |
#else
|
|
Packit |
785658 |
import Data.Array.Base (unsafeAt)
|
|
Packit |
785658 |
#endif
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Calculate the integer logarithm for an arbitrary base.
|
|
Packit |
785658 |
-- The base must be greater than 1, the second argument, the number
|
|
Packit |
785658 |
-- whose logarithm is sought, must be positive, otherwise an error is thrown.
|
|
Packit |
785658 |
-- If @base == 2@, the specialised version is called, which is more
|
|
Packit |
785658 |
-- efficient than the general algorithm.
|
|
Packit |
785658 |
--
|
|
Packit |
785658 |
-- Satisfies:
|
|
Packit |
785658 |
--
|
|
Packit |
785658 |
-- > base ^ integerLogBase base m <= m < base ^ (integerLogBase base m + 1)
|
|
Packit |
785658 |
--
|
|
Packit |
785658 |
-- for @base > 1@ and @m > 0@.
|
|
Packit |
785658 |
integerLogBase :: Integer -> Integer -> Int
|
|
Packit |
785658 |
integerLogBase b n
|
|
Packit |
785658 |
| n < 1 = error "Math.NumberTheory.Logarithms.integerLogBase: argument must be positive."
|
|
Packit |
785658 |
| n < b = 0
|
|
Packit |
785658 |
| b == 2 = integerLog2' n
|
|
Packit |
785658 |
| b < 2 = error "Math.NumberTheory.Logarithms.integerLogBase: base must be greater than one."
|
|
Packit |
785658 |
| otherwise = integerLogBase' b n
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Calculate the integer logarithm of an 'Integer' to base 2.
|
|
Packit |
785658 |
-- The argument must be positive, otherwise an error is thrown.
|
|
Packit |
785658 |
integerLog2 :: Integer -> Int
|
|
Packit |
785658 |
integerLog2 n
|
|
Packit |
785658 |
| n < 1 = error "Math.NumberTheory.Logarithms.integerLog2: argument must be positive"
|
|
Packit |
785658 |
| otherwise = I# (integerLog2# n)
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Cacluate the integer logarithm for an arbitrary base.
|
|
Packit |
785658 |
-- The base must be greater than 1, the second argument, the number
|
|
Packit |
785658 |
-- whose logarithm is sought, must be positive, otherwise an error is thrown.
|
|
Packit |
785658 |
-- If @base == 2@, the specialised version is called, which is more
|
|
Packit |
785658 |
-- efficient than the general algorithm.
|
|
Packit |
785658 |
--
|
|
Packit |
785658 |
-- Satisfies:
|
|
Packit |
785658 |
--
|
|
Packit |
785658 |
-- > base ^ integerLogBase base m <= m < base ^ (integerLogBase base m + 1)
|
|
Packit |
785658 |
--
|
|
Packit |
785658 |
-- for @base > 1@ and @m > 0@.
|
|
Packit |
785658 |
naturalLogBase :: Natural -> Natural -> Int
|
|
Packit |
785658 |
naturalLogBase b n
|
|
Packit |
785658 |
| n < 1 = error "Math.NumberTheory.Logarithms.naturalLogBase: argument must be positive."
|
|
Packit |
785658 |
| n < b = 0
|
|
Packit |
785658 |
| b == 2 = naturalLog2' n
|
|
Packit |
785658 |
| b < 2 = error "Math.NumberTheory.Logarithms.naturalLogBase: base must be greater than one."
|
|
Packit |
785658 |
| otherwise = naturalLogBase' b n
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Calculate the natural logarithm of an 'Natural' to base 2.
|
|
Packit |
785658 |
-- The argument must be non-zero, otherwise an error is thrown.
|
|
Packit |
785658 |
naturalLog2 :: Natural -> Int
|
|
Packit |
785658 |
naturalLog2 n
|
|
Packit |
785658 |
| n < 1 = error "Math.NumberTheory.Logarithms.naturalLog2: argument must be non-zero"
|
|
Packit |
785658 |
| otherwise = I# (naturalLog2# n)
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Calculate the integer logarithm of an 'Int' to base 2.
|
|
Packit |
785658 |
-- The argument must be positive, otherwise an error is thrown.
|
|
Packit |
785658 |
intLog2 :: Int -> Int
|
|
Packit |
785658 |
intLog2 (I# i#)
|
|
Packit |
785658 |
| isTrue# (i# <# 1#) = error "Math.NumberTheory.Logarithms.intLog2: argument must be positive"
|
|
Packit |
785658 |
| otherwise = I# (wordLog2# (int2Word# i#))
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Calculate the integer logarithm of a 'Word' to base 2.
|
|
Packit |
785658 |
-- The argument must be positive, otherwise an error is thrown.
|
|
Packit |
785658 |
wordLog2 :: Word -> Int
|
|
Packit |
785658 |
wordLog2 (W# w#)
|
|
Packit |
785658 |
| isTrue# (w# `eqWord#` 0##) = error "Math.NumberTheory.Logarithms.wordLog2: argument must not be 0."
|
|
Packit |
785658 |
| otherwise = I# (wordLog2# w#)
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Same as 'integerLog2', but without checks, saves a little time when
|
|
Packit |
785658 |
-- called often for known good input.
|
|
Packit |
785658 |
integerLog2' :: Integer -> Int
|
|
Packit |
785658 |
integerLog2' n = I# (integerLog2# n)
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Same as 'naturalLog2', but without checks, saves a little time when
|
|
Packit |
785658 |
-- called often for known good input.
|
|
Packit |
785658 |
naturalLog2' :: Natural -> Int
|
|
Packit |
785658 |
naturalLog2' n = I# (naturalLog2# n)
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Same as 'intLog2', but without checks, saves a little time when
|
|
Packit |
785658 |
-- called often for known good input.
|
|
Packit |
785658 |
intLog2' :: Int -> Int
|
|
Packit |
785658 |
intLog2' (I# i#) = I# (wordLog2# (int2Word# i#))
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Same as 'wordLog2', but without checks, saves a little time when
|
|
Packit |
785658 |
-- called often for known good input.
|
|
Packit |
785658 |
wordLog2' :: Word -> Int
|
|
Packit |
785658 |
wordLog2' (W# w#) = I# (wordLog2# w#)
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Calculate the integer logarithm of an 'Integer' to base 10.
|
|
Packit |
785658 |
-- The argument must be positive, otherwise an error is thrown.
|
|
Packit |
785658 |
integerLog10 :: Integer -> Int
|
|
Packit |
785658 |
integerLog10 n
|
|
Packit |
785658 |
| n < 1 = error "Math.NumberTheory.Logarithms.integerLog10: argument must be positive"
|
|
Packit |
785658 |
| otherwise = integerLog10' n
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Calculate the integer logarithm of an 'Integer' to base 10.
|
|
Packit |
785658 |
-- The argument must be not zero, otherwise an error is thrown.
|
|
Packit |
785658 |
naturalLog10 :: Natural -> Int
|
|
Packit |
785658 |
naturalLog10 n
|
|
Packit |
785658 |
| n < 1 = error "Math.NumberTheory.Logarithms.naturalaLog10: argument must be non-zero"
|
|
Packit |
785658 |
| otherwise = naturalLog10' n
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Same as 'integerLog10', but without a check for a positive
|
|
Packit |
785658 |
-- argument. Saves a little time when called often for known good
|
|
Packit |
785658 |
-- input.
|
|
Packit |
785658 |
integerLog10' :: Integer -> Int
|
|
Packit |
785658 |
integerLog10' n
|
|
Packit |
785658 |
| n < 10 = 0
|
|
Packit |
785658 |
| n < 100 = 1
|
|
Packit |
785658 |
| otherwise = ex + integerLog10' (n `quot` 10 ^ ex)
|
|
Packit |
785658 |
where
|
|
Packit |
785658 |
ln = I# (integerLog2# n)
|
|
Packit |
785658 |
-- u/v is a good approximation of log 2/log 10
|
|
Packit |
785658 |
u = 1936274
|
|
Packit |
785658 |
v = 6432163
|
|
Packit |
785658 |
-- so ex is a good approximation to integerLogBase 10 n
|
|
Packit |
785658 |
ex = fromInteger ((u * fromIntegral ln) `quot` v)
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Same as 'naturalLog10', but without a check for a positive
|
|
Packit |
785658 |
-- argument. Saves a little time when called often for known good
|
|
Packit |
785658 |
-- input.
|
|
Packit |
785658 |
naturalLog10' :: Natural -> Int
|
|
Packit |
785658 |
naturalLog10' n
|
|
Packit |
785658 |
| n < 10 = 0
|
|
Packit |
785658 |
| n < 100 = 1
|
|
Packit |
785658 |
| otherwise = ex + naturalLog10' (n `quot` 10 ^ ex)
|
|
Packit |
785658 |
where
|
|
Packit |
785658 |
ln = I# (naturalLog2# n)
|
|
Packit |
785658 |
-- u/v is a good approximation of log 2/log 10
|
|
Packit |
785658 |
u = 1936274
|
|
Packit |
785658 |
v = 6432163
|
|
Packit |
785658 |
-- so ex is a good approximation to naturalLogBase 10 n
|
|
Packit |
785658 |
ex = fromInteger ((u * fromIntegral ln) `quot` v)
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Same as 'integerLogBase', but without checks, saves a little time when
|
|
Packit |
785658 |
-- called often for known good input.
|
|
Packit |
785658 |
integerLogBase' :: Integer -> Integer -> Int
|
|
Packit |
785658 |
integerLogBase' b n
|
|
Packit |
785658 |
| n < b = 0
|
|
Packit |
785658 |
| ln-lb < lb = 1 -- overflow safe version of ln < 2*lb, implies n < b*b
|
|
Packit |
785658 |
| b < 33 = let bi = fromInteger b
|
|
Packit |
785658 |
ix = 2*bi-4
|
|
Packit |
785658 |
-- u/v is a good approximation of log 2/log b
|
|
Packit |
785658 |
u = logArr `unsafeAt` ix
|
|
Packit |
785658 |
v = logArr `unsafeAt` (ix+1)
|
|
Packit |
785658 |
-- hence ex is a rather good approximation of integerLogBase b n
|
|
Packit |
785658 |
-- most of the time, it will already be exact
|
|
Packit |
785658 |
ex = fromInteger ((fromIntegral u * fromIntegral ln) `quot` fromIntegral v)
|
|
Packit |
785658 |
in case u of
|
|
Packit |
785658 |
1 -> ln `quot` v -- a power of 2, easy
|
|
Packit |
785658 |
_ -> ex + integerLogBase' b (n `quot` b ^ ex)
|
|
Packit |
785658 |
| otherwise = let -- shift b so that 16 <= bi < 32
|
|
Packit |
785658 |
bi = fromInteger (b `shiftR` (lb-4))
|
|
Packit |
785658 |
-- we choose an approximation of log 2 / log (bi+1) to
|
|
Packit |
785658 |
-- be sure we underestimate
|
|
Packit |
785658 |
ix = 2*bi-2
|
|
Packit |
785658 |
-- u/w is a reasonably good approximation to log 2/log b
|
|
Packit |
785658 |
-- it is too small, but not by much, so the recursive call
|
|
Packit |
785658 |
-- should most of the time be caught by one of the first
|
|
Packit |
785658 |
-- two guards unless n is huge, but then it'd still be
|
|
Packit |
785658 |
-- a call with a much smaller second argument.
|
|
Packit |
785658 |
u = fromIntegral $ logArr `unsafeAt` ix
|
|
Packit |
785658 |
v = fromIntegral $ logArr `unsafeAt` (ix+1)
|
|
Packit |
785658 |
w = v + u*fromIntegral (lb-4)
|
|
Packit |
785658 |
ex = fromInteger ((u * fromIntegral ln) `quot` w)
|
|
Packit |
785658 |
in ex + integerLogBase' b (n `quot` b ^ ex)
|
|
Packit |
785658 |
where
|
|
Packit |
785658 |
lb = integerLog2' b
|
|
Packit |
785658 |
ln = integerLog2' n
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- | Same as 'naturalLogBase', but without checks, saves a little time when
|
|
Packit |
785658 |
-- called often for known good input.
|
|
Packit |
785658 |
naturalLogBase' :: Natural -> Natural -> Int
|
|
Packit |
785658 |
naturalLogBase' b n
|
|
Packit |
785658 |
| n < b = 0
|
|
Packit |
785658 |
| ln-lb < lb = 1 -- overflow safe version of ln < 2*lb, implies n < b*b
|
|
Packit |
785658 |
| b < 33 = let bi = fromIntegral b
|
|
Packit |
785658 |
ix = 2*bi-4
|
|
Packit |
785658 |
-- u/v is a good approximation of log 2/log b
|
|
Packit |
785658 |
u = logArr `unsafeAt` ix
|
|
Packit |
785658 |
v = logArr `unsafeAt` (ix+1)
|
|
Packit |
785658 |
-- hence ex is a rather good approximation of integerLogBase b n
|
|
Packit |
785658 |
-- most of the time, it will already be exact
|
|
Packit |
785658 |
ex = fromNatural ((fromIntegral u * fromIntegral ln) `quot` fromIntegral v)
|
|
Packit |
785658 |
in case u of
|
|
Packit |
785658 |
1 -> ln `quot` v -- a power of 2, easy
|
|
Packit |
785658 |
_ -> ex + naturalLogBase' b (n `quot` b ^ ex)
|
|
Packit |
785658 |
| otherwise = let -- shift b so that 16 <= bi < 32
|
|
Packit |
785658 |
bi = fromNatural (b `shiftR` (lb-4))
|
|
Packit |
785658 |
-- we choose an approximation of log 2 / log (bi+1) to
|
|
Packit |
785658 |
-- be sure we underestimate
|
|
Packit |
785658 |
ix = 2*bi-2
|
|
Packit |
785658 |
-- u/w is a reasonably good approximation to log 2/log b
|
|
Packit |
785658 |
-- it is too small, but not by much, so the recursive call
|
|
Packit |
785658 |
-- should most of the time be caught by one of the first
|
|
Packit |
785658 |
-- two guards unless n is huge, but then it'd still be
|
|
Packit |
785658 |
-- a call with a much smaller second argument.
|
|
Packit |
785658 |
u = fromIntegral $ logArr `unsafeAt` ix
|
|
Packit |
785658 |
v = fromIntegral $ logArr `unsafeAt` (ix+1)
|
|
Packit |
785658 |
w = v + u*fromIntegral (lb-4)
|
|
Packit |
785658 |
ex = fromNatural ((u * fromIntegral ln) `quot` w)
|
|
Packit |
785658 |
in ex + naturalLogBase' b (n `quot` b ^ ex)
|
|
Packit |
785658 |
where
|
|
Packit |
785658 |
lb = naturalLog2' b
|
|
Packit |
785658 |
ln = naturalLog2' n
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-- Lookup table for logarithms of 2 <= k <= 32
|
|
Packit |
785658 |
-- In each row "x , y", x/y is a good rational approximation of log 2 / log k.
|
|
Packit |
785658 |
-- For the powers of 2, it is exact, otherwise x/y < log 2/log k, since we don't
|
|
Packit |
785658 |
-- want to overestimate integerLogBase b n = floor $ (log 2/log b)*logBase 2 n.
|
|
Packit |
785658 |
logArr :: UArray Int Int
|
|
Packit |
785658 |
logArr = listArray (0, 61)
|
|
Packit |
785658 |
[ 1 , 1,
|
|
Packit |
785658 |
190537 , 301994,
|
|
Packit |
785658 |
1 , 2,
|
|
Packit |
785658 |
1936274 , 4495889,
|
|
Packit |
785658 |
190537 , 492531,
|
|
Packit |
785658 |
91313 , 256348,
|
|
Packit |
785658 |
1 , 3,
|
|
Packit |
785658 |
190537 , 603988,
|
|
Packit |
785658 |
1936274 , 6432163,
|
|
Packit |
785658 |
1686227 , 5833387,
|
|
Packit |
785658 |
190537 , 683068,
|
|
Packit |
785658 |
5458 , 20197,
|
|
Packit |
785658 |
91313 , 347661,
|
|
Packit |
785658 |
416263 , 1626294,
|
|
Packit |
785658 |
1 , 4,
|
|
Packit |
785658 |
32631 , 133378,
|
|
Packit |
785658 |
190537 , 794525,
|
|
Packit |
785658 |
163451 , 694328,
|
|
Packit |
785658 |
1936274 , 8368437,
|
|
Packit |
785658 |
1454590 , 6389021,
|
|
Packit |
785658 |
1686227 , 7519614,
|
|
Packit |
785658 |
785355 , 3552602,
|
|
Packit |
785658 |
190537 , 873605,
|
|
Packit |
785658 |
968137 , 4495889,
|
|
Packit |
785658 |
5458 , 25655,
|
|
Packit |
785658 |
190537 , 905982,
|
|
Packit |
785658 |
91313 , 438974,
|
|
Packit |
785658 |
390321 , 1896172,
|
|
Packit |
785658 |
416263 , 2042557,
|
|
Packit |
785658 |
709397 , 3514492,
|
|
Packit |
785658 |
1 , 5
|
|
Packit |
785658 |
]
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-------------------------------------------------------------------------------
|
|
Packit |
785658 |
-- Unsafe
|
|
Packit |
785658 |
-------------------------------------------------------------------------------
|
|
Packit |
785658 |
|
|
Packit |
785658 |
#if CheckBounds
|
|
Packit |
785658 |
unsafeAt :: (IArray a e, Ix i) => a i e -> i -> e
|
|
Packit |
785658 |
unsafeAt = (!)`
|
|
Packit |
785658 |
#endif
|
|
Packit |
785658 |
|
|
Packit |
785658 |
-------------------------------------------------------------------------------
|
|
Packit |
785658 |
-- Natural helpers
|
|
Packit |
785658 |
-------------------------------------------------------------------------------
|
|
Packit |
785658 |
|
|
Packit |
785658 |
fromNatural :: Num a => Natural -> a
|
|
Packit |
785658 |
fromNatural = fromIntegral
|
|
Packit |
785658 |
|
|
Packit |
785658 |
naturalLog2# :: Natural -> Int#
|
|
Packit |
785658 |
#if Base48 && defined(MIN_VERSION_integer_gmp)
|
|
Packit |
785658 |
naturalLog2# (NatS# b) = wordLog2# b
|
|
Packit |
785658 |
naturalLog2# (NatJ# n) = integerLog2# (Jp# n)
|
|
Packit |
785658 |
#else
|
|
Packit |
785658 |
naturalLog2# n = integerLog2# (toInteger n)
|
|
Packit |
785658 |
#endif
|
|
Packit |
785658 |
|
|
Packit |
785658 |
#if __GLASGOW_HASKELL__ < 707
|
|
Packit |
785658 |
-- The times they are a-changing. The types of primops too :(
|
|
Packit |
785658 |
isTrue# :: Bool -> Bool
|
|
Packit |
785658 |
isTrue# = id
|
|
Packit |
785658 |
#endif
|