integer-roots 1.0.1.0 → 1.0.2.0
raw patch · 7 files changed
+148/−27 lines, 7 filesdep +ghc-bignumPVP ok
version bump matches the API change (PVP)
Dependencies added: ghc-bignum
API changes (from Hackage documentation)
Files
- Math/NumberTheory/Roots/Cubes.hs +25/−6
- Math/NumberTheory/Roots/Fourth.hs +24/−5
- Math/NumberTheory/Roots/General.hs +46/−5
- Math/NumberTheory/Roots/Squares/Internal.hs +26/−4
- changelog.md +4/−0
- integer-roots.cabal +8/−7
- test-suite/Math/NumberTheory/Roots/GeneralTests.hs +15/−0
Math/NumberTheory/Roots/Cubes.hs view
@@ -8,6 +8,7 @@ -- cube roots and testing for cubeness. {-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP #-} {-# LANGUAGE MagicHash #-} module Math.NumberTheory.Roots.Cubes@@ -20,10 +21,21 @@ ) where import Data.Bits (finiteBitSize, (.&.))-import GHC.Exts (Int#, Ptr(..), quotInt#, int2Double#, double2Int#, isTrue#, (/##), (**##), (<#), (*#), (-#))+import GHC.Exts (Int#, Ptr(..), int2Double#, double2Int#, isTrue#, (/##), (**##), (<#))+import Numeric.Natural (Natural)++#ifdef MIN_VERSION_integer_gmp+import GHC.Exts (quotInt#, (*#), (-#)) import GHC.Integer.GMP.Internals (Integer(..), shiftLInteger, shiftRInteger, sizeofBigNat#) import GHC.Integer.Logarithms (integerLog2#)-import Numeric.Natural (Natural)+#define IS S#+#define IP Jp#+#define bigNatSize sizeofBigNat+#else+import GHC.Exts (minusWord#, timesWord#, quotWord#)+import GHC.Num.BigNat (bigNatSize#)+import GHC.Num.Integer (Integer(..), integerLog2#, integerShiftR#, integerShiftL#)+#endif import Math.NumberTheory.Utils.BitMask (indexBitSet) @@ -197,16 +209,23 @@ -- | approximate cube root, about 50 bits should be correct for large numbers appCuRt :: Integer -> Integer-appCuRt (S# i#) = case double2Int# (int2Double# i# **## (1.0## /## 3.0##)) of- r# -> S# r#-appCuRt n@(Jp# bn#)- | isTrue# ((sizeofBigNat# bn#) <# thresh#) =+appCuRt (IS i#) = case double2Int# (int2Double# i# **## (1.0## /## 3.0##)) of+ r# -> IS r#+appCuRt n@(IP bn#)+ | isTrue# ((bigNatSize# bn#) <# thresh#) = floor (fromInteger n ** (1.0/3.0) :: Double) | otherwise = case integerLog2# n of+#ifdef MIN_VERSION_integer_gmp l# -> case (l# `quotInt#` 3#) -# 51# of h# -> case shiftRInteger n (3# *# h#) of m -> case floor (fromInteger m ** (1.0/3.0) :: Double) of r -> shiftLInteger r h#+#else+ l# -> case (l# `quotWord#` 3##) `minusWord#` 51## of+ h# -> case integerShiftR# n (3## `timesWord#` h#) of+ m -> case floor (fromInteger m ** (1.0/3.0) :: Double) of+ r -> integerShiftL# r h#+#endif where -- threshold for shifting vs. direct fromInteger -- we shift when we expect more than 256 bits
Math/NumberTheory/Roots/Fourth.hs view
@@ -7,6 +7,7 @@ -- Functions dealing with fourth powers. Efficient calculation of integer fourth -- roots and efficient testing for being a square's square. +{-# LANGUAGE CPP #-} {-# LANGUAGE MagicHash #-} module Math.NumberTheory.Roots.Fourth@@ -19,10 +20,21 @@ ) where import Data.Bits (finiteBitSize, (.&.))-import GHC.Exts (Int#, Ptr(..), uncheckedIShiftRA#, int2Double#, double2Int#, isTrue#, sqrtDouble#, (<#), (*#), (-#))+import GHC.Exts (Int#, Ptr(..), int2Double#, double2Int#, isTrue#, sqrtDouble#, (<#))+import Numeric.Natural (Natural)++#ifdef MIN_VERSION_integer_gmp+import GHC.Exts (uncheckedIShiftRA#, (*#), (-#)) import GHC.Integer.GMP.Internals (Integer(..), shiftLInteger, shiftRInteger, sizeofBigNat#) import GHC.Integer.Logarithms (integerLog2#)-import Numeric.Natural (Natural)+#define IS S#+#define IP Jp#+#define bigNatSize sizeofBigNat+#else+import GHC.Exts (uncheckedShiftRL#, minusWord#, timesWord#)+import GHC.Num.BigNat (bigNatSize#)+import GHC.Num.Integer (Integer(..), integerLog2#, integerShiftR#, integerShiftL#)+#endif import Math.NumberTheory.Utils.BitMask (indexBitSet) @@ -130,15 +142,22 @@ -- Find a fairly good approximation to the fourth root. -- About 48 bits should be correct for large Integers. appBiSqrt :: Integer -> Integer-appBiSqrt (S# i#) = S# (double2Int# (sqrtDouble# (sqrtDouble# (int2Double# i#))))-appBiSqrt n@(Jp# bn#)- | isTrue# ((sizeofBigNat# bn#) <# thresh#) =+appBiSqrt (IS i#) = IS (double2Int# (sqrtDouble# (sqrtDouble# (int2Double# i#))))+appBiSqrt n@(IP bn#)+ | isTrue# ((bigNatSize# bn#) <# thresh#) = floor (sqrt . sqrt $ fromInteger n :: Double) | otherwise = case integerLog2# n of+#ifdef MIN_VERSION_integer_gmp l# -> case uncheckedIShiftRA# l# 2# -# 47# of h# -> case shiftRInteger n (4# *# h#) of m -> case floor (sqrt $ sqrt $ fromInteger m :: Double) of r -> shiftLInteger r h#+#else+ l# -> case uncheckedShiftRL# l# 2# `minusWord#` 47## of+ h# -> case integerShiftR# n (4## `timesWord#` h#) of+ m -> case floor (sqrt $ sqrt $ fromInteger m :: Double) of+ r -> integerShiftL# r h#+#endif where -- threshold for shifting vs. direct fromInteger -- we shift when we expect more than 256 bits
Math/NumberTheory/Roots/General.hs view
@@ -11,6 +11,7 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE CPP #-} {-# LANGUAGE MagicHash #-}+{-# LANGUAGE ViewPatterns #-} module Math.NumberTheory.Roots.General ( integerRoot@@ -20,19 +21,29 @@ , highestPower ) where +#include "MachDeps.h"+ import Data.Bits (countTrailingZeros, shiftL, shiftR) import Data.List (foldl', sortBy) import Data.Maybe (isJust)-import GHC.Exts (Int(..), Word(..), quotInt#, int2Word#, word2Int#, int2Double#, double2Int#, isTrue#, Ptr(..), indexWord16OffAddr#, (/##), (**##), (<#), (*#), (-#), (+#))+import GHC.Exts (Int(..), Word(..), word2Int#, int2Double#, double2Int#, isTrue#, Ptr(..), indexWord16OffAddr#, (/##), (**##)) #if MIN_VERSION_base(4,16,0) import GHC.Exts (word16ToWord#) #endif #ifdef WORDS_BIGENDIAN-import GHC.Exts (byteSwap#)+import GHC.Exts (byteSwap16#) #endif+import Numeric.Natural (Natural)++#ifdef MIN_VERSION_integer_gmp+import GHC.Exts (int2Word#, quotInt#, (<#), (*#), (-#), (+#)) import GHC.Integer.GMP.Internals (Integer(..), shiftLInteger, shiftRInteger) import GHC.Integer.Logarithms (integerLog2#)-import Numeric.Natural (Natural)+#define IS S#+#else+import GHC.Exts (plusWord#, minusWord#, timesWord#, quotWord#, ltWord#)+import GHC.Num.Integer (Integer(..), integerLog2#, integerShiftR#, integerShiftL#)+#endif import qualified Math.NumberTheory.Roots.Squares as P2 import qualified Math.NumberTheory.Roots.Cubes as P3@@ -85,7 +96,11 @@ a = appKthRoot (fromIntegral k) (toInteger n) kTooLarge = (toInteger k /= toInteger (fromIntegral k `asTypeOf` n)) -- k doesn't fit in n's type || (toInteger k > toInteger (maxBound :: Int)) -- 2^k doesn't fit in Integer+#ifdef MIN_VERSION_integer_gmp || (I# (integerLog2# (toInteger n)) < fromIntegral k) -- n < 2^k+#else+ || (W# (integerLog2# (toInteger n)) < fromIntegral k) -- n < 2^k+#endif -- | For a positive exponent \( k \) -- calculate the exact integer \( k \)-th root of the second argument if it exists,@@ -139,7 +154,11 @@ ok = r^k == n kTooLarge = (toInteger k /= toInteger (fromIntegral k `asTypeOf` n)) -- k doesn't fit in n's type || (toInteger k > toInteger (maxBound :: Int)) -- 2^k doesn't fit in Integer+#ifdef MIN_VERSION_integer_gmp || (I# (integerLog2# (toInteger n)) < fromIntegral k) -- n < 2^k+#else+ || (W# (integerLog2# (toInteger n)) < fromIntegral k) -- n < 2^k+#endif -- | For a positive exponent \( k \) test whether the second argument -- is a perfect \( k \)-th power.@@ -208,7 +227,8 @@ -- find an approximation to the k-th root -- here, k > 4 and n > 31 appKthRoot :: Int -> Integer -> Integer-appKthRoot (I# k#) (S# n#) = S# (double2Int# (int2Double# n# **## (1.0## /## int2Double# k#)))+appKthRoot (I# k#) (IS n#) = IS (double2Int# (int2Double# n# **## (1.0## /## int2Double# k#)))+#ifdef MIN_VERSION_integer_gmp appKthRoot k@(I# k#) n | k >= 256 = 1 `shiftLInteger` (integerLog2# n `quotInt#` k# +# 1#) | otherwise =@@ -224,6 +244,23 @@ | otherwise -> floor (scaleFloat 400 (fromInteger (n `shiftRInteger` (h# *# k#)) ** (1/fromIntegral k) :: Double)) `shiftLInteger` (h# -# 400#)+#else+appKthRoot k@(fromIntegral -> W# k#) n+ | k >= 256 = 1 `integerShiftL#` (integerLog2# n `quotWord#` k# `plusWord#` 1##)+ | otherwise =+ case integerLog2# n of+ l# -> case l# `quotWord#` k# of+ 0## -> 1+ 1## -> 3+ 2## -> 5+ 3## -> 11+ h# | isTrue# (h# `ltWord#` 500##) ->+ floor (scaleFloat (I# (word2Int# h#))+ (fromInteger (n `integerShiftR#` (h# `timesWord#` k#)) ** (1/fromIntegral k) :: Double))+ | otherwise ->+ floor (scaleFloat 400 (fromInteger (n `integerShiftR#` (h# `timesWord#` k#)) ** (1/fromIntegral k) :: Double))+ `integerShiftL#` (h# `minusWord#` 400##)+#endif -- assumption: argument is > 1 integerHighPower :: Integer -> (Integer, Word)@@ -258,7 +295,7 @@ smallOddPrimes :: [Integer] smallOddPrimes = takeWhile (< spBound)- $ map (\(I# k#) -> S# (word2Int# (+ $ map (\(I# k#) -> IS (word2Int# ( #if MIN_VERSION_base(4,16,0) #ifdef WORDS_BIGENDIAN byteSwap16# (word16ToWord# (indexWord16OffAddr# smallPrimesAddr# k#))@@ -290,7 +327,11 @@ | e == 0 = rawPower maxExp n | otherwise = go divs where+#ifdef MIN_VERSION_integer_gmp maxExp = (W# (int2Word# (integerLog2# n))) `quot` spBEx+#else+ maxExp = (W# (integerLog2# n)) `quot` spBEx+#endif divs = divisorsTo maxExp e go [] = (foldl' (*) n [p^ex | (p,ex) <- pws], 1) go (d:ds) = case exactRoot d n of
Math/NumberTheory/Roots/Squares/Internal.hs view
@@ -7,6 +7,7 @@ -- Internal functions dealing with square roots. End-users should not import this module. {-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP #-} {-# LANGUAGE MagicHash #-} module Math.NumberTheory.Roots.Squares.Internal@@ -16,9 +17,19 @@ import Data.Bits (finiteBitSize, unsafeShiftL, unsafeShiftR, (.&.), (.|.)) -import GHC.Exts (Int(..), Int#, uncheckedIShiftRA#, isTrue#, int2Double#, sqrtDouble#, double2Int#, (<#), (*#), (-#))+import GHC.Exts (Int(..), Int#, isTrue#, int2Double#, sqrtDouble#, double2Int#, (<#))+#ifdef MIN_VERSION_integer_gmp+import GHC.Exts (uncheckedIShiftRA#, (*#), (-#)) import GHC.Integer.GMP.Internals (Integer(..), shiftLInteger, shiftRInteger, sizeofBigNat#) import GHC.Integer.Logarithms (integerLog2#)+#define IS S#+#define IP Jp#+#define bigNatSize sizeofBigNat+#else+import GHC.Exts (uncheckedShiftRL#, word2Int#, minusWord#, timesWord#)+import GHC.Num.BigNat (bigNatSize#)+import GHC.Num.Integer (Integer(..), integerLog2#, integerShiftR#, integerShiftL#)+#endif -- Find approximation to square root in 'Integer', then -- find the integer square root by the integer variant@@ -47,15 +58,22 @@ -- At most one off for small Integers, about 48 bits should be correct -- for large Integers. appSqrt :: Integer -> Integer-appSqrt (S# i#) = S# (double2Int# (sqrtDouble# (int2Double# i#)))-appSqrt n@(Jp# bn#)- | isTrue# ((sizeofBigNat# bn#) <# thresh#) =+appSqrt (IS i#) = IS (double2Int# (sqrtDouble# (int2Double# i#)))+appSqrt n@(IP bn#)+ | isTrue# ((bigNatSize# bn#) <# thresh#) = floor (sqrt $ fromInteger n :: Double) | otherwise = case integerLog2# n of+#ifdef MIN_VERSION_integer_gmp l# -> case uncheckedIShiftRA# l# 1# -# 47# of h# -> case shiftRInteger n (2# *# h#) of m -> case floor (sqrt $ fromInteger m :: Double) of r -> shiftLInteger r h#+#else+ l# -> case uncheckedShiftRL# l# 1# `minusWord#` 47## of+ h# -> case integerShiftR# n (2## `timesWord#` h#) of+ m -> case floor (sqrt $ fromInteger m :: Double) of+ r -> integerShiftL# r h#+#endif where -- threshold for shifting vs. direct fromInteger -- we shift when we expect more than 256 bits@@ -90,7 +108,11 @@ in (s `unsafeShiftR` 1, r' `unsafeShiftR` 2) where k = lgN `unsafeShiftR` 2 + 1+#ifdef MIN_VERSION_integer_gmp lgN = I# (integerLog2# n)+#else+ lgN = I# (word2Int# (integerLog2# n))+#endif karatsubaStep :: Int -> (Integer, Integer, Integer, Integer) -> (Integer, Integer) karatsubaStep k (a3, a2, a1, a0)
changelog.md view
@@ -1,3 +1,7 @@+# 1.0.2.0++* More fixes for big-endian architectures.+ # 1.0.1.0 * Fixes for big-endian architectures.
integer-roots.cabal view
@@ -1,10 +1,10 @@ name: integer-roots-version: 1.0.1.0+version: 1.0.2.0 cabal-version: >=1.10 build-type: Simple license: MIT license-file: LICENSE-copyright: (c) 2011 Daniel Fischer, 2016-2020 Andrew Lelechenko.+copyright: (c) 2011 Daniel Fischer, 2016-2021 Andrew Lelechenko. maintainer: Andrew Lelechenko andrew dot lelechenko at gmail dot com homepage: https://github.com/Bodigrim/integer-roots bug-reports: https://github.com/Bodigrim/integer-roots/issues@@ -23,8 +23,11 @@ library build-depends:- base >=4.9 && <5,- integer-gmp <1.2+ base >=4.9 && <5+ if impl(ghc < 9.0)+ build-depends: integer-gmp <1.2+ else+ build-depends: ghc-bignum < 1.3 exposed-modules: Math.NumberTheory.Roots other-modules:@@ -37,7 +40,7 @@ Math.NumberTheory.Utils.BitMask Math.NumberTheory.Utils.FromIntegral default-language: Haskell2010- ghc-options: -Wall -Widentities -Wno-deprecations -Wcompat+ ghc-options: -Wall -Widentities -Wcompat test-suite integer-roots-tests build-depends:@@ -62,8 +65,6 @@ ghc-options: -Wall -Widentities -Wcompat test-suite integer-roots-doctests- if impl(ghc >= 9.2)- buildable: False type: exitcode-stdio-1.0 main-is: Doctest.hs hs-source-dirs: test-suite
test-suite/Math/NumberTheory/Roots/GeneralTests.hs view
@@ -65,6 +65,16 @@ (b, k) = highestPower n (b', k') = highestPower b +highestPowerProperty2 :: Integral a => AnySign a -> Bool+highestPowerProperty2 (AnySign n) = case k of+ 1 -> not (isSquare n) && not (isCube n)+ 2 -> isSquare n && not (isCube n)+ 3 -> n + 1 `elem` [0, 1, 2] || (not (isSquare n) && isCube n)+ 4 -> isSquare n && not (isCube n)+ _ -> all (\l -> isKthPower l n == (k `mod` l == 0)) [1..k]+ where+ (_, k) = highestPower n+ highestPowerSpecialCases :: [Assertion] highestPowerSpecialCases = -- Freezes before d44a13b.@@ -76,6 +86,10 @@ , 1013582159576576 , 7) + , a ( 9 :: Int+ , 3+ , 2)+ , a ( -2 ^ 63 :: Int , -2 :: Int , 63)@@ -125,6 +139,7 @@ , testIntegralProperty "isPerfectPower" isPerfectPowerProperty , testGroup "highestPower" ( testIntegralProperty "highestPower" highestPowerProperty+ : testIntegralProperty "highestPower 2" highestPowerProperty2 : zipWith (\i a -> testCase ("special case " ++ show i) a) [1..] highestPowerSpecialCases ) ]