diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,16 @@ +Copyright (c) 2011 Daniel Fischer, 2017 Oleg Grenrus + +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and + associated documentation files (the "Software"), to deal in the Software without restriction, + including without limitation the rights to use, copy, modify, merge, publish, distribute, + sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all copies or +substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT +LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,5 @@ +module Main where + +import Distribution.Simple + +main = defaultMain diff --git a/integer-logarithms.cabal b/integer-logarithms.cabal new file mode 100644 --- /dev/null +++ b/integer-logarithms.cabal @@ -0,0 +1,101 @@ +name: integer-logarithms +version: 1 +cabal-version: >= 1.10 +author: Daniel Fischer +copyright: (c) 2011 Daniel Fischer +license: MIT +license-file: LICENSE +maintainer: Oleg Grenrus +build-type: Simple +stability: Provisional +homepage: https://github.com/phadej/integer-logarithms +bug-reports: https://github.com/phadej/integer-logarithms/issues + +synopsis: Integer logarithms. +description: + "Math.NumberTheory.Logarithms" and "Math.NumberTheory.Powers.Integer" + from the arithmoi package. + . + Also provides "GHC.Integer.Logarithms.Compat" and + "Math.NumberTheory.Power.Natural" modules, as well as some + additional functions in migrated modules. + +category: Math, Algorithms, Number Theory + +tested-with : + GHC==7.0.4, + GHC==7.2.2, + GHC==7.4.2, + GHC==7.6.3, + GHC==7.8.4, + GHC==7.10.3, + GHC==8.0.1 + +extra-source-files : readme.md + +flag check-bounds + description: Replace unsafe array operations with safe ones + default: False + manual: True + +library + default-language: Haskell2010 + hs-source-dirs: src + build-depends: + base >= 4.3 && < 4.10, + array >= 0.3 && < 0.6, + ghc-prim < 0.6, + integer-gmp < 1.1 + if impl(ghc >= 7.10) + cpp-options: -DBase48 + else + build-depends: nats >= 1.1 && <1.2 + + exposed-modules: + Math.NumberTheory.Logarithms + Math.NumberTheory.Powers.Integer + Math.NumberTheory.Powers.Natural + GHC.Integer.Logarithms.Compat + other-extensions: + BangPatterns + CPP + MagicHash + + ghc-options: -O2 -Wall + if flag(check-bounds) + cpp-options: -DCheckBounds + +source-repository head + type: git + location: https://github.com/phadej/integer-logarithms + +test-suite spec + type: exitcode-stdio-1.0 + hs-source-dirs: test-suite + ghc-options: -Wall + main-is: Test.hs + default-language: Haskell2010 + other-extensions: + StandaloneDeriving + FlexibleContexts + FlexibleInstances + GeneralizedNewtypeDeriving + MultiParamTypeClasses + build-depends: + base, + integer-logarithms, + tasty >= 0.10 && < 0.12, + tasty-smallcheck >= 0.8 && < 0.9, + tasty-quickcheck >= 0.8 && < 0.9, + tasty-hunit >= 0.9 && < 0.10, + QuickCheck >= 2.9 && < 2.10, + smallcheck >= 1.1 && < 1.2 + if !impl(ghc >= 7.10) + build-depends: nats >= 1.1 && <1.2 + + other-modules: + Math.NumberTheory.LogarithmsTests + Math.NumberTheory.Powers.IntegerTests + Math.NumberTheory.Powers.NaturalTests + other-modules: + Math.NumberTheory.TestUtils diff --git a/readme.md b/readme.md new file mode 100644 --- /dev/null +++ b/readme.md @@ -0,0 +1,3 @@ +# integer-logarithms + +`Math.NumberTheory.Logarithms` splitted out of [`arithmoi`](http://hackage.haskell.org/package/arithmoi) diff --git a/src/GHC/Integer/Logarithms/Compat.hs b/src/GHC/Integer/Logarithms/Compat.hs new file mode 100644 --- /dev/null +++ b/src/GHC/Integer/Logarithms/Compat.hs @@ -0,0 +1,152 @@ +-- | +-- Module: GHC.Integer.Logarithms.Compat +-- Copyright: (c) 2011 Daniel Fischer +-- Licence: MIT +-- Maintainer: Daniel Fischer +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Low level stuff for integer logarithms. +{-# LANGUAGE CPP, MagicHash, UnboxedTuples #-} +module GHC.Integer.Logarithms.Compat + ( -- * Functions + integerLogBase# + , integerLog2# + , wordLog2# + ) where + +#if __GLASGOW_HASKELL__ >= 702 + +-- Stuff is already there +import GHC.Integer.Logarithms + +#else + +-- We have to define it here +#include "MachDeps.h" + +import GHC.Base +import GHC.Integer.GMP.Internals + +#if (WORD_SIZE_IN_BITS != 32) && (WORD_SIZE_IN_BITS != 64) +#error Only word sizes 32 and 64 are supported. +#endif + + +#if WORD_SIZE_IN_BITS == 32 + +#define WSHIFT 5 +#define MMASK 31 + +#else + +#define WSHIFT 6 +#define MMASK 63 + +#endif + +-- Reference implementation only, the algorithm in M.NT.Logarithms is better. + +-- | Calculate the integer logarithm for an arbitrary base. +-- The base must be greater than 1, the second argument, the number +-- whose logarithm is sought; should be positive, otherwise the +-- result is meaningless. +-- +-- > base ^ integerLogBase# base m <= m < base ^ (integerLogBase# base m + 1) +-- +-- for @base > 1@ and @m > 0@. +integerLogBase# :: Integer -> Integer -> Int# +integerLogBase# b m = case step b of + (# _, e #) -> e + where + step pw = + if m < pw + then (# m, 0# #) + else case step (pw * pw) of + (# q, e #) -> + if q < pw + then (# q, 2# *# e #) + else (# q `quot` pw, 2# *# e +# 1# #) + +-- | Calculate the integer base 2 logarithm of an 'Integer'. +-- The calculation is much more efficient than for the general case. +-- +-- The argument must be strictly positive, that condition is /not/ checked. +integerLog2# :: Integer -> Int# +integerLog2# (S# i) = wordLog2# (int2Word# i) +integerLog2# (J# s ba) = check (s -# 1#) + where + check i = case indexWordArray# ba i of + 0## -> check (i -# 1#) + w -> wordLog2# w +# (uncheckedIShiftL# i WSHIFT#) + +-- | This function calculates the integer base 2 logarithm of a 'Word#'. +-- @'wordLog2#' 0## = -1#@. +{-# INLINE wordLog2# #-} +wordLog2# :: Word# -> Int# +wordLog2# w = + case leadingZeros of + BA lz -> + let zeros u = indexInt8Array# lz (word2Int# u) in +#if WORD_SIZE_IN_BITS == 64 + case uncheckedShiftRL# w 56# of + a -> + if a `neWord#` 0## + then 64# -# zeros a + else + case uncheckedShiftRL# w 48# of + b -> + if b `neWord#` 0## + then 56# -# zeros b + else + case uncheckedShiftRL# w 40# of + c -> + if c `neWord#` 0## + then 48# -# zeros c + else + case uncheckedShiftRL# w 32# of + d -> + if d `neWord#` 0## + then 40# -# zeros d + else +#endif + case uncheckedShiftRL# w 24# of + e -> + if e `neWord#` 0## + then 32# -# zeros e + else + case uncheckedShiftRL# w 16# of + f -> + if f `neWord#` 0## + then 24# -# zeros f + else + case uncheckedShiftRL# w 8# of + g -> + if g `neWord#` 0## + then 16# -# zeros g + else 8# -# zeros w + +-- Lookup table +data BA = BA ByteArray# + +leadingZeros :: BA +leadingZeros = + let mkArr s = + case newByteArray# 256# s of + (# s1, mba #) -> + case writeInt8Array# mba 0# 9# s1 of + s2 -> + let fillA lim val idx st = + if idx ==# 256# + then st + else if idx <# lim + then case writeInt8Array# mba idx val st of + nx -> fillA lim val (idx +# 1#) nx + else fillA (2# *# lim) (val -# 1#) idx st + in case fillA 2# 8# 1# s2 of + s3 -> case unsafeFreezeByteArray# mba s3 of + (# _, ba #) -> ba + in case mkArr realWorld# of + b -> BA b + +#endif diff --git a/src/Math/NumberTheory/Logarithms.hs b/src/Math/NumberTheory/Logarithms.hs new file mode 100644 --- /dev/null +++ b/src/Math/NumberTheory/Logarithms.hs @@ -0,0 +1,330 @@ +-- | +-- Module: Math.NumberTheory.Logarithms +-- Copyright: (c) 2011 Daniel Fischer +-- Licence: MIT +-- Maintainer: Daniel Fischer +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Integer Logarithms. For efficiency, the internal representation of 'Integer's +-- from integer-gmp is used. +-- +{-# LANGUAGE CPP #-} +{-# LANGUAGE MagicHash #-} +module Math.NumberTheory.Logarithms + ( -- * Integer logarithms with input checks + integerLogBase + , integerLog2 + , integerLog10 + + , naturalLogBase + , naturalLog2 + , naturalLog10 + + , intLog2 + , wordLog2 + + -- * Integer logarithms without input checks + -- + -- | These functions are total, however, don't rely on the values with out-of-domain arguments. + , integerLogBase' + , integerLog2' + , integerLog10' + + , intLog2' + , wordLog2' + ) where + +import GHC.Exts + +import Data.Bits +import Data.Array.Unboxed +import Numeric.Natural + +import GHC.Integer.Logarithms.Compat +#if Base48 +import GHC.Integer.GMP.Internals (Integer (..)) +import GHC.Natural +#endif + +#if CheckBounds +import Data.Array.IArray (IArray, (!)) +#else +import Data.Array.Base (unsafeAt) +#endif + +import Math.NumberTheory.Powers.Integer +import Math.NumberTheory.Powers.Natural + +-- | Calculate the integer logarithm for an arbitrary base. +-- The base must be greater than 1, the second argument, the number +-- whose logarithm is sought, must be positive, otherwise an error is thrown. +-- If @base == 2@, the specialised version is called, which is more +-- efficient than the general algorithm. +-- +-- Satisfies: +-- +-- > base ^ integerLogBase base m <= m < base ^ (integerLogBase base m + 1) +-- +-- for @base > 1@ and @m > 0@. +integerLogBase :: Integer -> Integer -> Int +integerLogBase b n + | n < 1 = error "Math.NumberTheory.Logarithms.integerLogBase: argument must be positive." + | n < b = 0 + | b == 2 = integerLog2' n + | b < 2 = error "Math.NumberTheory.Logarithms.integerLogBase: base must be greater than one." + | otherwise = integerLogBase' b n + +-- | Calculate the integer logarithm of an 'Integer' to base 2. +-- The argument must be positive, otherwise an error is thrown. +integerLog2 :: Integer -> Int +integerLog2 n + | n < 1 = error "Math.NumberTheory.Logarithms.integerLog2: argument must be positive" + | otherwise = I# (integerLog2# n) + +-- | Cacluate the integer logarithm for an arbitrary base. +-- The base must be greater than 1, the second argument, the number +-- whose logarithm is sought, must be positive, otherwise an error is thrown. +-- If @base == 2@, the specialised version is called, which is more +-- efficient than the general algorithm. +-- +-- Satisfies: +-- +-- > base ^ integerLogBase base m <= m < base ^ (integerLogBase base m + 1) +-- +-- for @base > 1@ and @m > 0@. +naturalLogBase :: Natural -> Natural -> Int +naturalLogBase b n + | n < 1 = error "Math.NumberTheory.Logarithms.naturalLogBase: argument must be positive." + | n < b = 0 + | b == 2 = naturalLog2' n + | b < 2 = error "Math.NumberTheory.Logarithms.naturalLogBase: base must be greater than one." + | otherwise = naturalLogBase' b n + +-- | Calculate the natural logarithm of an 'Natural' to base 2. +-- The argument must be non-zero, otherwise an error is thrown. +naturalLog2 :: Natural -> Int +naturalLog2 n + | n < 1 = error "Math.NumberTheory.Logarithms.naturalLog2: argument must be non-zero" + | otherwise = I# (naturalLog2# n) + +-- | Calculate the integer logarithm of an 'Int' to base 2. +-- The argument must be positive, otherwise an error is thrown. +intLog2 :: Int -> Int +intLog2 (I# i#) + | isTrue# (i# <# 1#) = error "Math.NumberTheory.Logarithms.intLog2: argument must be positive" + | otherwise = I# (wordLog2# (int2Word# i#)) + +-- | Calculate the integer logarithm of a 'Word' to base 2. +-- The argument must be positive, otherwise an error is thrown. +wordLog2 :: Word -> Int +wordLog2 (W# w#) + | isTrue# (w# `eqWord#` 0##) = error "Math.NumberTheory.Logarithms.wordLog2: argument must not be 0." + | otherwise = I# (wordLog2# w#) + +-- | Same as 'integerLog2', but without checks, saves a little time when +-- called often for known good input. +integerLog2' :: Integer -> Int +integerLog2' n = I# (integerLog2# n) + +-- | Same as 'naturalLog2', but without checks, saves a little time when +-- called often for known good input. +naturalLog2' :: Natural -> Int +naturalLog2' n = I# (naturalLog2# n) + +-- | Same as 'intLog2', but without checks, saves a little time when +-- called often for known good input. +intLog2' :: Int -> Int +intLog2' (I# i#) = I# (wordLog2# (int2Word# i#)) + +-- | Same as 'wordLog2', but without checks, saves a little time when +-- called often for known good input. +wordLog2' :: Word -> Int +wordLog2' (W# w#) = I# (wordLog2# w#) + +-- | Calculate the integer logarithm of an 'Integer' to base 10. +-- The argument must be positive, otherwise an error is thrown. +integerLog10 :: Integer -> Int +integerLog10 n + | n < 1 = error "Math.NumberTheory.Logarithms.integerLog10: argument must be positive" + | otherwise = integerLog10' n + +-- | Calculate the integer logarithm of an 'Integer' to base 10. +-- The argument must be not zero, otherwise an error is thrown. +naturalLog10 :: Natural -> Int +naturalLog10 n + | n < 1 = error "Math.NumberTheory.Logarithms.naturalaLog10: argument must be non-zero" + | otherwise = naturalLog10' n + +-- | Same as 'integerLog10', but without a check for a positive +-- argument. Saves a little time when called often for known good +-- input. +integerLog10' :: Integer -> Int +integerLog10' n + | n < 10 = 0 + | n < 100 = 1 + | otherwise = ex + integerLog10' (n `quot` integerPower 10 ex) + where + ln = I# (integerLog2# n) + -- u/v is a good approximation of log 2/log 10 + u = 1936274 + v = 6432163 + -- so ex is a good approximation to integerLogBase 10 n + ex = fromInteger ((u * fromIntegral ln) `quot` v) + +-- | Same as 'naturalLog10', but without a check for a positive +-- argument. Saves a little time when called often for known good +-- input. +naturalLog10' :: Natural -> Int +naturalLog10' n + | n < 10 = 0 + | n < 100 = 1 + | otherwise = ex + naturalLog10' (n `quot` naturalPower 10 ex) + where + ln = I# (naturalLog2# n) + -- u/v is a good approximation of log 2/log 10 + u = 1936274 + v = 6432163 + -- so ex is a good approximation to naturalLogBase 10 n + ex = fromInteger ((u * fromIntegral ln) `quot` v) + +-- | Same as 'integerLogBase', but without checks, saves a little time when +-- called often for known good input. +integerLogBase' :: Integer -> Integer -> Int +integerLogBase' b n + | n < b = 0 + | ln-lb < lb = 1 -- overflow safe version of ln < 2*lb, implies n < b*b + | b < 33 = let bi = fromInteger b + ix = 2*bi-4 + -- u/v is a good approximation of log 2/log b + u = logArr `unsafeAt` ix + v = logArr `unsafeAt` (ix+1) + -- hence ex is a rather good approximation of integerLogBase b n + -- most of the time, it will already be exact + ex = fromInteger ((fromIntegral u * fromIntegral ln) `quot` fromIntegral v) + in case u of + 1 -> ln `quot` v -- a power of 2, easy + _ -> ex + integerLogBase' b (n `quot` integerPower b ex) + | otherwise = let -- shift b so that 16 <= bi < 32 + bi = fromInteger (b `shiftR` (lb-4)) + -- we choose an approximation of log 2 / log (bi+1) to + -- be sure we underestimate + ix = 2*bi-2 + -- u/w is a reasonably good approximation to log 2/log b + -- it is too small, but not by much, so the recursive call + -- should most of the time be caught by one of the first + -- two guards unless n is huge, but then it'd still be + -- a call with a much smaller second argument. + u = fromIntegral $ logArr `unsafeAt` ix + v = fromIntegral $ logArr `unsafeAt` (ix+1) + w = v + u*fromIntegral (lb-4) + ex = fromInteger ((u * fromIntegral ln) `quot` w) + in ex + integerLogBase' b (n `quot` integerPower b ex) + where + lb = integerLog2' b + ln = integerLog2' n + +-- | Same as 'naturalLogBase', but without checks, saves a little time when +-- called often for known good input. +naturalLogBase' :: Natural -> Natural -> Int +naturalLogBase' b n + | n < b = 0 + | ln-lb < lb = 1 -- overflow safe version of ln < 2*lb, implies n < b*b + | b < 33 = let bi = fromIntegral b + ix = 2*bi-4 + -- u/v is a good approximation of log 2/log b + u = logArr `unsafeAt` ix + v = logArr `unsafeAt` (ix+1) + -- hence ex is a rather good approximation of integerLogBase b n + -- most of the time, it will already be exact + ex = fromNatural ((fromIntegral u * fromIntegral ln) `quot` fromIntegral v) + in case u of + 1 -> ln `quot` v -- a power of 2, easy + _ -> ex + naturalLogBase' b (n `quot` naturalPower b ex) + | otherwise = let -- shift b so that 16 <= bi < 32 + bi = fromNatural (b `shiftR` (lb-4)) + -- we choose an approximation of log 2 / log (bi+1) to + -- be sure we underestimate + ix = 2*bi-2 + -- u/w is a reasonably good approximation to log 2/log b + -- it is too small, but not by much, so the recursive call + -- should most of the time be caught by one of the first + -- two guards unless n is huge, but then it'd still be + -- a call with a much smaller second argument. + u = fromIntegral $ logArr `unsafeAt` ix + v = fromIntegral $ logArr `unsafeAt` (ix+1) + w = v + u*fromIntegral (lb-4) + ex = fromNatural ((u * fromIntegral ln) `quot` w) + in ex + naturalLogBase' b (n `quot` naturalPower b ex) + where + lb = naturalLog2' b + ln = naturalLog2' n + +-- Lookup table for logarithms of 2 <= k <= 32 +-- In each row "x , y", x/y is a good rational approximation of log 2 / log k. +-- For the powers of 2, it is exact, otherwise x/y < log 2/log k, since we don't +-- want to overestimate integerLogBase b n = floor $ (log 2/log b)*logBase 2 n. +logArr :: UArray Int Int +logArr = listArray (0, 61) + [ 1 , 1, + 190537 , 301994, + 1 , 2, + 1936274 , 4495889, + 190537 , 492531, + 91313 , 256348, + 1 , 3, + 190537 , 603988, + 1936274 , 6432163, + 1686227 , 5833387, + 190537 , 683068, + 5458 , 20197, + 91313 , 347661, + 416263 , 1626294, + 1 , 4, + 32631 , 133378, + 190537 , 794525, + 163451 , 694328, + 1936274 , 8368437, + 1454590 , 6389021, + 1686227 , 7519614, + 785355 , 3552602, + 190537 , 873605, + 968137 , 4495889, + 5458 , 25655, + 190537 , 905982, + 91313 , 438974, + 390321 , 1896172, + 416263 , 2042557, + 709397 , 3514492, + 1 , 5 + ] + +------------------------------------------------------------------------------- +-- Unsafe +------------------------------------------------------------------------------- + +#if CheckBounds +unsafeAt :: (IArray a e, Ix i) => a i e -> i -> e +unsafeAt = (!)` +#endif + +------------------------------------------------------------------------------- +-- Natural helpers +------------------------------------------------------------------------------- + +fromNatural :: Num a => Natural -> a +fromNatural = fromIntegral + +naturalLog2# :: Natural -> Int# +#if Base48 +naturalLog2# (NatS# b) = wordLog2# b +naturalLog2# (NatJ# n) = integerLog2# (Jp# n) +#else +naturalLog2# n = integerLog2# (toInteger n) +#endif + +#if __GLASGOW_HASKELL__ < 707 +-- The times they are a-changing. The types of primops too :( +isTrue# :: Bool -> Bool +isTrue# = id +#endif diff --git a/src/Math/NumberTheory/Powers/Integer.hs b/src/Math/NumberTheory/Powers/Integer.hs new file mode 100644 --- /dev/null +++ b/src/Math/NumberTheory/Powers/Integer.hs @@ -0,0 +1,62 @@ +-- | +-- Module: Math.NumberTheory.Powers.Integer +-- Copyright: (c) 2011-2014 Daniel Fischer +-- Licence: MIT +-- Maintainer: Daniel Fischer +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Potentially faster power function for 'Integer' base and 'Int' +-- or 'Word' exponent. +-- +{-# LANGUAGE CPP #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE BangPatterns #-} +module Math.NumberTheory.Powers.Integer + ( integerPower + , integerWordPower + ) where + +import GHC.Exts +import GHC.Integer.Logarithms.Compat (wordLog2#) + +-- | Power of an 'Integer' by the left-to-right repeated squaring algorithm. +-- This needs two multiplications in each step while the right-to-left +-- algorithm needs only one multiplication for 0-bits, but here the +-- two factors always have approximately the same size, which on average +-- gains a bit when the result is large. +-- +-- For small results, it is unlikely to be any faster than '(^)', quite +-- possibly slower (though the difference shouldn't be large), and for +-- exponents with few bits set, the same holds. But for exponents with +-- many bits set, the speedup can be significant. +-- +-- /Warning:/ No check for the negativity of the exponent is performed, +-- a negative exponent is interpreted as a large positive exponent. +integerPower :: Integer -> Int -> Integer +integerPower b (I# e#) = power b (int2Word# e#) + +-- | Same as 'integerPower', but for exponents of type 'Word'. +integerWordPower :: Integer -> Word -> Integer +integerWordPower b (W# w#) = power b w# + +power :: Integer -> Word# -> Integer +power b w# + | isTrue# (w# `eqWord#` 0##) = 1 + | isTrue# (w# `eqWord#` 1##) = b + | otherwise = go (wordLog2# w# -# 1#) b (b*b) + where + go 0# l h = if isTrue# ((w# `and#` 1##) `eqWord#` 0##) then l*l else (l*h) + go i# l h + | w# `hasBit#` i# = go (i# -# 1#) (l*h) (h*h) + | otherwise = go (i# -# 1#) (l*l) (l*h) + +-- | A raw version of testBit for 'Word#'. +hasBit# :: Word# -> Int# -> Bool +hasBit# w# i# = isTrue# (((w# `uncheckedShiftRL#` i#) `and#` 1##) `neWord#` 0##) + +#if __GLASGOW_HASKELL__ < 707 +-- The times they are a-changing. The types of primops too :( +isTrue# :: Bool -> Bool +isTrue# = id +#endif diff --git a/src/Math/NumberTheory/Powers/Natural.hs b/src/Math/NumberTheory/Powers/Natural.hs new file mode 100644 --- /dev/null +++ b/src/Math/NumberTheory/Powers/Natural.hs @@ -0,0 +1,63 @@ +-- | +-- Module: Math.NumberTheory.Powers.Natural +-- Copyright: (c) 2011-2014 Daniel Fischer +-- Licence: MIT +-- Maintainer: Daniel Fischer +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Potentially faster power function for 'Natural' base and 'Int' +-- or 'Word' exponent. +-- +{-# LANGUAGE CPP #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE BangPatterns #-} +module Math.NumberTheory.Powers.Natural + ( naturalPower + , naturalWordPower + ) where + +import GHC.Exts +import Numeric.Natural +import GHC.Integer.Logarithms.Compat (wordLog2#) + +-- | Power of an 'Natural' by the left-to-right repeated squaring algorithm. +-- This needs two multiplications in each step while the right-to-left +-- algorithm needs only one multiplication for 0-bits, but here the +-- two factors always have approximately the same size, which on average +-- gains a bit when the result is large. +-- +-- For small results, it is unlikely to be any faster than '(^)', quite +-- possibly slower (though the difference shouldn't be large), and for +-- exponents with few bits set, the same holds. But for exponents with +-- many bits set, the speedup can be significant. +-- +-- /Warning:/ No check for the negativity of the exponent is performed, +-- a negative exponent is interpreted as a large positive exponent. +naturalPower :: Natural -> Int -> Natural +naturalPower b (I# e#) = power b (int2Word# e#) + +-- | Same as 'naturalPower', but for exponents of type 'Word'. +naturalWordPower :: Natural -> Word -> Natural +naturalWordPower b (W# w#) = power b w# + +power :: Natural -> Word# -> Natural +power b w# + | isTrue# (w# `eqWord#` 0##) = 1 + | isTrue# (w# `eqWord#` 1##) = b + | otherwise = go (wordLog2# w# -# 1#) b (b*b) + where + go 0# l h = if isTrue# ((w# `and#` 1##) `eqWord#` 0##) then l*l else (l*h) + go i# l h + | w# `hasBit#` i# = go (i# -# 1#) (l*h) (h*h) + | otherwise = go (i# -# 1#) (l*l) (l*h) + +-- | A raw version of testBit for 'Word#'. +hasBit# :: Word# -> Int# -> Bool +hasBit# w# i# = isTrue# (((w# `uncheckedShiftRL#` i#) `and#` 1##) `neWord#` 0##) + +#if __GLASGOW_HASKELL__ < 707 +-- The times they are a-changing. The types of primops too :( +isTrue# :: Bool -> Bool +isTrue# = id +#endif diff --git a/test-suite/Math/NumberTheory/LogarithmsTests.hs b/test-suite/Math/NumberTheory/LogarithmsTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/LogarithmsTests.hs @@ -0,0 +1,146 @@ +-- | +-- Module: Math.NumberTheory.LogarithmsTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- +-- Tests for Math.NumberTheory.Logarithms +-- + +{-# LANGUAGE CPP #-} + +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.LogarithmsTests + ( testSuite + ) where + +import Test.Tasty + +#if MIN_VERSION_base(4,8,0) +#else +import Data.Word +#endif +import Numeric.Natural + +import Math.NumberTheory.Logarithms +import Math.NumberTheory.TestUtils + +-- | Check that 'integerLogBase' returns the largest integer @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@. +integerLogBaseProperty :: Positive Integer -> Positive Integer -> Bool +integerLogBaseProperty (Positive b) (Positive n) = b == 1 || b ^ l <= n && b ^ (l + 1) > n + where + l = toInteger $ integerLogBase b n + +-- | Check that 'integerLog2' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. +integerLog2Property :: Positive Integer -> Bool +integerLog2Property (Positive n) = 2 ^ l <= n && 2 ^ (l + 1) > n + where + l = toInteger $ integerLog2 n + +integerLog2HugeProperty :: Huge (Positive Integer) -> Bool +integerLog2HugeProperty (Huge (Positive n)) = 2 ^ l <= n && 2 ^ (l + 1) > n + where + l = toInteger $ integerLog2 n + +-- | Check that 'integerLog10' returns the largest integer @l@ such that 10 ^ @l@ <= @n@ and 10 ^ (@l@+1) > @n@. +integerLog10Property :: Positive Integer -> Bool +integerLog10Property (Positive n) = 10 ^ l <= n && 10 ^ (l + 1) > n + where + l = toInteger $ integerLog10 n + +-- | Check that 'naturalLogBase' returns the largest natural @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@. +naturalLogBaseProperty :: Positive Natural -> Positive Natural -> Bool +naturalLogBaseProperty (Positive b) (Positive n) = b == 1 || b ^ l <= n && b ^ (l + 1) > n + where + l = fromIntegral $ naturalLogBase b n + +-- | Check that 'naturalLog2' returns the largest natural @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. +naturalLog2Property :: Positive Natural -> Bool +naturalLog2Property (Positive n) = 2 ^ l <= n && 2 ^ (l + 1) > n + where + l = fromIntegral $ naturalLog2 n + +naturalLog2HugeProperty :: Huge (Positive Natural) -> Bool +naturalLog2HugeProperty (Huge (Positive n)) = 2 ^ l <= n && 2 ^ (l + 1) > n + where + l = fromIntegral $ naturalLog2 n + +-- | Check that 'naturalLog10' returns the largest natural @l@ such that 10 ^ @l@ <= @n@ and 10 ^ (@l@+1) > @n@. +naturalLog10Property :: Positive Natural -> Bool +naturalLog10Property (Positive n) = 10 ^ l <= n && 10 ^ (l + 1) > n + where + l = fromIntegral $ naturalLog10 n + +-- | Check that 'intLog2' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. +intLog2Property :: Positive Int -> Bool +intLog2Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2) + where + l = intLog2 n + +-- | Check that 'wordLog2' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. +wordLog2Property :: Positive Word -> Bool +wordLog2Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2) + where + l = wordLog2 n + +-- | Check that 'integerLogBase'' returns the largest integer @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@. +integerLogBase'Property :: Positive Integer -> Positive Integer -> Bool +integerLogBase'Property (Positive b) (Positive n) = b == 1 || b ^ l <= n && b ^ (l + 1) > n + where + l = toInteger $ integerLogBase' b n + +-- | Check that 'integerLogBase'' returns the largest integer @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@ for @b@ > 32 and @n@ >= @b@ ^ 2. +integerLogBase'Property2 :: Positive Integer -> Positive Integer -> Bool +integerLogBase'Property2 (Positive b') (Positive n') = b ^ l <= n && b ^ (l + 1) > n + where + b = b' + 32 + n = n' + b ^ 2 - 1 + l = toInteger $ integerLogBase' b n + +-- | Check that 'integerLog2'' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. +integerLog2'Property :: Positive Integer -> Bool +integerLog2'Property (Positive n) = 2 ^ l <= n && 2 ^ (l + 1) > n + where + l = toInteger $ integerLog2' n + +-- | Check that 'integerLog10'' returns the largest integer @l@ such that 10 ^ @l@ <= @n@ and 10 ^ (@l@+1) > @n@. +integerLog10'Property :: Positive Integer -> Bool +integerLog10'Property (Positive n) = 10 ^ l <= n && 10 ^ (l + 1) > n + where + l = toInteger $ integerLog10' n + +-- | Check that 'intLog2'' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. +intLog2'Property :: Positive Int -> Bool +intLog2'Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2) + where + l = intLog2' n + +-- | Check that 'wordLog2'' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. +wordLog2'Property :: Positive Word -> Bool +wordLog2'Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2) + where + l = wordLog2' n + +testSuite :: TestTree +testSuite = testGroup "Logarithms" + [ testSmallAndQuick "integerLogBase" integerLogBaseProperty + , testSmallAndQuick "integerLog2" integerLog2Property + , testSmallAndQuick "integerLog2Huge" integerLog2HugeProperty + , testSmallAndQuick "integerLog10" integerLog10Property + , testSmallAndQuick "naturalLogBase" naturalLogBaseProperty + , testSmallAndQuick "naturalLog2" naturalLog2Property + , testSmallAndQuick "naturalLog2Huge" naturalLog2HugeProperty + , testSmallAndQuick "naturalLog10" naturalLog10Property + , testSmallAndQuick "intLog2" intLog2Property + , testSmallAndQuick "wordLog2" wordLog2Property + + , testSmallAndQuick "integerLogBase'" integerLogBase'Property + , testSmallAndQuick "integerLogBase' with base > 32 and n >= base ^ 2" + integerLogBase'Property2 + , testSmallAndQuick "integerLog2'" integerLog2'Property + , testSmallAndQuick "integerLog10'" integerLog10'Property + , testSmallAndQuick "intLog2'" intLog2'Property + , testSmallAndQuick "wordLog2'" wordLog2'Property + ] diff --git a/test-suite/Math/NumberTheory/Powers/IntegerTests.hs b/test-suite/Math/NumberTheory/Powers/IntegerTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Powers/IntegerTests.hs @@ -0,0 +1,41 @@ +-- | +-- Module: Math.NumberTheory.Powers.IntegerTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- +-- Tests for Math.NumberTheory.Powers.Integer +-- + +{-# LANGUAGE CPP #-} + +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.Powers.IntegerTests + ( testSuite + ) where + +import Test.Tasty + +#if MIN_VERSION_base(4,8,0) +#else +import Data.Word +#endif + +import Math.NumberTheory.Powers.Integer +import Math.NumberTheory.TestUtils + +-- | Check that 'integerPower' == '^'. +integerPowerProperty :: Integer -> Power Int -> Bool +integerPowerProperty a (Power b) = integerPower a b == a ^ b + +-- | Check that 'integerWordPower' == '^'. +integerWordPowerProperty :: Integer -> Power Word -> Bool +integerWordPowerProperty a (Power b) = integerWordPower a b == a ^ b + +testSuite :: TestTree +testSuite = testGroup "Integer" + [ testSmallAndQuick "integerPower" integerPowerProperty + , testSmallAndQuick "integerWordPower" integerWordPowerProperty + ] diff --git a/test-suite/Math/NumberTheory/Powers/NaturalTests.hs b/test-suite/Math/NumberTheory/Powers/NaturalTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Powers/NaturalTests.hs @@ -0,0 +1,42 @@ +-- | +-- Module: Math.NumberTheory.Powers.NaturalTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- +-- Tests for Math.NumberTheory.Powers.Natural +-- + +{-# LANGUAGE CPP #-} + +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.Powers.NaturalTests + ( testSuite + ) where + +import Test.Tasty + +#if MIN_VERSION_base(4,8,0) +#else +import Data.Word +#endif + +import Math.NumberTheory.Powers.Natural +import Numeric.Natural +import Math.NumberTheory.TestUtils + +-- | Check that 'naturalPower' == '^'. +naturalPowerProperty :: Natural -> Power Int -> Bool +naturalPowerProperty a (Power b) = naturalPower a b == a ^ b + +-- | Check that 'naturalWordPower' == '^'. +naturalWordPowerProperty :: Natural -> Power Word -> Bool +naturalWordPowerProperty a (Power b) = naturalWordPower a b == a ^ b + +testSuite :: TestTree +testSuite = testGroup "Natural" + [ testSmallAndQuick "naturalPower" naturalPowerProperty + , testSmallAndQuick "naturalWordPower" naturalWordPowerProperty + ] diff --git a/test-suite/Math/NumberTheory/TestUtils.hs b/test-suite/Math/NumberTheory/TestUtils.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/TestUtils.hs @@ -0,0 +1,96 @@ +-- | +-- Module: Math.NumberTheory.TestUtils +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- + +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE StandaloneDeriving #-} + +{-# OPTIONS_GHC -fno-warn-orphans #-} + +module Math.NumberTheory.TestUtils + ( module Test.SmallCheck.Series + , Power (..) + , Huge (..) + , testSmallAndQuick + ) where + +import Test.SmallCheck.Series (cons2) +import Test.Tasty +import Test.Tasty.SmallCheck as SC +import Test.Tasty.QuickCheck as QC hiding (Positive, NonNegative, generate, getNonNegative) +import Test.SmallCheck.Series (Positive(..), NonNegative(..), Serial(..), Series, generate) + +import Control.Applicative +import Data.Word +import Numeric.Natural + +testSmallAndQuick + :: SC.Testable IO a + => QC.Testable a + => String -> a -> TestTree +testSmallAndQuick name f = testGroup name + [ SC.testProperty "smallcheck" f + , QC.testProperty "quickcheck" f + ] + +------------------------------------------------------------------------------- +-- Serial monadic actions + +instance Monad m => Serial m Word where + series = + generate (\d -> if d >= 0 then pure 0 else empty) <|> nats + where + nats = generate $ \d -> if d > 0 then [1 .. fromInteger (toInteger d)] else empty + +instance Monad m => Serial m Natural where + series = + generate (\d -> if d >= 0 then pure 0 else empty) <|> nats + where + nats = generate $ \d -> if d > 0 then [1 .. fromInteger (toInteger d)] else empty + +------------------------------------------------------------------------------- +-- Power + +newtype Power a = Power { getPower :: a } + deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real) + +instance (Monad m, Num a, Ord a, Serial m a) => Serial m (Power a) where + series = Power <$> series `suchThatSerial` (> 0) + +instance (Num a, Ord a, Integral a, Arbitrary a) => Arbitrary (Power a) where + arbitrary = Power <$> (getSmall <$> arbitrary) `suchThat` (> 0) + shrink (Power x) = Power <$> filter (> 0) (shrink x) + +suchThatSerial :: Series m a -> (a -> Bool) -> Series m a +suchThatSerial s p = s >>= \x -> if p x then pure x else empty + +------------------------------------------------------------------------------- +-- Huge + +newtype Huge a = Huge { getHuge :: a } + deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real) + +instance (Num a, Arbitrary a) => Arbitrary (Huge a) where + arbitrary = do + Positive l <- arbitrary + ds <- vector (l :: Int) + return $ Huge $ foldl1 (\acc n -> acc * 2^(63 :: Int) + n) ds + +-- | maps 'Huge' constructor over series +instance Serial m a => Serial m (Huge a) where + series = fmap Huge series + +------------------------------------------------------------------------------- +-- Positive from smallcheck + +instance (Num a, Ord a, Arbitrary a) => Arbitrary (Positive a) where + arbitrary = Positive <$> (arbitrary `suchThat` (> 0)) + shrink (Positive x) = Positive <$> filter (> 0) (shrink x) diff --git a/test-suite/Test.hs b/test-suite/Test.hs new file mode 100644 --- /dev/null +++ b/test-suite/Test.hs @@ -0,0 +1,15 @@ +import Test.Tasty + +import qualified Math.NumberTheory.LogarithmsTests as Logarithms +import qualified Math.NumberTheory.Powers.IntegerTests as PowerInteger +import qualified Math.NumberTheory.Powers.NaturalTests as PowerNatural + +main :: IO () +main = defaultMain tests + +tests :: TestTree +tests = testGroup "All" + [ Logarithms.testSuite + , PowerInteger.testSuite + , PowerNatural.testSuite + ]