diff --git a/Math/NumberTheory/Roots/Cubes.hs b/Math/NumberTheory/Roots/Cubes.hs --- a/Math/NumberTheory/Roots/Cubes.hs +++ b/Math/NumberTheory/Roots/Cubes.hs @@ -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 diff --git a/Math/NumberTheory/Roots/Fourth.hs b/Math/NumberTheory/Roots/Fourth.hs --- a/Math/NumberTheory/Roots/Fourth.hs +++ b/Math/NumberTheory/Roots/Fourth.hs @@ -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 diff --git a/Math/NumberTheory/Roots/General.hs b/Math/NumberTheory/Roots/General.hs --- a/Math/NumberTheory/Roots/General.hs +++ b/Math/NumberTheory/Roots/General.hs @@ -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 diff --git a/Math/NumberTheory/Roots/Squares/Internal.hs b/Math/NumberTheory/Roots/Squares/Internal.hs --- a/Math/NumberTheory/Roots/Squares/Internal.hs +++ b/Math/NumberTheory/Roots/Squares/Internal.hs @@ -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) diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,7 @@ +# 1.0.2.0 + +* More fixes for big-endian architectures. + # 1.0.1.0 * Fixes for big-endian architectures. diff --git a/integer-roots.cabal b/integer-roots.cabal --- a/integer-roots.cabal +++ b/integer-roots.cabal @@ -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 diff --git a/test-suite/Math/NumberTheory/Roots/GeneralTests.hs b/test-suite/Math/NumberTheory/Roots/GeneralTests.hs --- a/test-suite/Math/NumberTheory/Roots/GeneralTests.hs +++ b/test-suite/Math/NumberTheory/Roots/GeneralTests.hs @@ -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 ) ]