arithmoi 0.4.0.3 → 0.4.0.4
raw patch · 16 files changed
+178/−86 lines, 16 filesdep ~array
Dependency ranges changed: array
Files
- Changes +4/−0
- Math/NumberTheory/GCD.hs +24/−12
- Math/NumberTheory/GCD/LowLevel.hs +15/−14
- Math/NumberTheory/Logarithms.hs +8/−4
- Math/NumberTheory/Moduli.hs +29/−15
- Math/NumberTheory/Powers/Cubes.hs +10/−2
- Math/NumberTheory/Powers/Fourth.hs +9/−2
- Math/NumberTheory/Powers/General.hs +26/−19
- Math/NumberTheory/Powers/Integer.hs +9/−6
- Math/NumberTheory/Powers/Squares.hs +5/−2
- Math/NumberTheory/Primes/Counting/Impl.hs +4/−1
- Math/NumberTheory/Primes/Factorisation/Montgomery.hs +14/−2
- Math/NumberTheory/Primes/Sieve/Eratosthenes.hs +4/−1
- Math/NumberTheory/Primes/Testing/Probabilistic.hs +2/−1
- Math/NumberTheory/Utils.hs +12/−1
- arithmoi.cabal +3/−4
Changes view
@@ -1,3 +1,7 @@+0.4.0.4:+ Update for GHC-7.8, the type of some primops changed, they return Int# now+ instead of Bool.+ Fixed bugs in modular square roots and factorisation. 0.4.0.3: Relaxed dependencies on mtl and containers Fixed warnings from GHC-7.5, Word(..) moved to GHC.Types
Math/NumberTheory/GCD.hs view
@@ -57,11 +57,8 @@ "binaryGCD/Int64" binaryGCD = gi64 "binaryGCD/Word64" binaryGCD = gw64 #-}-#else-{-# SPECIALISE binaryGCD :: Word64 -> Word64 -> Word64,- Int64 -> Int64 -> Int64 #-} #endif-{-# SPECIALISE binaryGCD :: Integer -> Integer -> Integer #-}+{-# INLINE [1] binaryGCD #-} -- | Calculate the greatest common divisor using the binary gcd algorithm. -- Depending on type and hardware, that can be considerably faster than -- @'Prelude.gcd'@ but it may also be significantly slower.@@ -75,9 +72,17 @@ -- -- Relies on twos complement or sign and magnitude representaion for signed types. binaryGCD :: (Integral a, Bits a) => a -> a -> a-binaryGCD a 0 = abs a-binaryGCD 0 b = abs b-binaryGCD a b =+binaryGCD = binaryGCDImpl++#if WORD_SIZE_IN_BITS < 64+{-# SPECIALISE binaryGCDImpl :: Word64 -> Word64 -> Word64,+ Int64 -> Int64 -> Int64 #-}+#endif+{-# SPECIALISE binaryGCDImpl :: Integer -> Integer -> Integer #-}+binaryGCDImpl :: (Integral a, Bits a) => a -> a -> a+binaryGCDImpl a 0 = abs a+binaryGCDImpl 0 b = abs b+binaryGCDImpl a b = case shiftToOddCount a of (!za, !oa) -> case shiftToOddCount b of@@ -139,11 +144,8 @@ "coprime/Int64" coprime = ci64 "coprime/Word64" coprime = cw64 #-}-#else-{-# SPECIALISE coprime :: Word64 -> Word64 -> Bool,- Int64 -> Int64 -> Bool #-} #endif-{-# SPECIALISE coprime :: Integer -> Integer -> Bool #-}+{-# INLINE [1] coprime #-} -- | Test whether two numbers are coprime using an abbreviated binary gcd algorithm. -- A little bit faster than checking @binaryGCD a b == 1@ if one of the arguments -- is even, much faster if both are even.@@ -153,7 +155,17 @@ -- -- Relies on twos complement or sign and magnitude representaion for signed types. coprime :: (Integral a, Bits a) => a -> a -> Bool-coprime a b =+coprime = coprimeImpl++-- Separate implementation to give the rules a chance to fire by not inlining+-- before phase 1, and yet have a specialisation for the types without rules+#if WORD_SIZE_IN_BITS < 64+{-# SPECIALISE coprimeImpl :: Word64 -> Word64 -> Bool,+ Int64 -> Int64 -> Bool #-}+#endif+{-# SPECIALISE coprimeImpl :: Integer -> Integer -> Bool #-}+coprimeImpl :: (Integral a, Bits a) => a -> a -> Bool+coprimeImpl a b = (a' == 1 || b' == 1) || (a' /= 0 && b' /= 0 && ((a .|. b) .&. 1) == 1 && gcdOdd (abs (shiftToOdd a')) (abs (shiftToOdd b')) == 1)
Math/NumberTheory/GCD/LowLevel.hs view
@@ -27,6 +27,7 @@ ) where import GHC.Base+ #if __GLASGOW_HASKELL__ < 705 import GHC.Word (Word(..)) -- Moved to GHC.Types #endif@@ -66,15 +67,15 @@ case shiftToOddCount# a# of (# za#, oa# #) -> case shiftToOddCount# b# of- (# zb#, ob# #) -> gcdWordOdd# oa# ob# `uncheckedShiftL#` (if za# <# zb# then za# else zb#)+ (# zb#, ob# #) -> gcdWordOdd# oa# ob# `uncheckedShiftL#` (if isTrue# (za# <# zb#) then za# else zb#) -- | Test whether two 'Word#'s are coprime. coprimeWord# :: Word# -> Word# -> Bool coprimeWord# a# b# =- (a# `eqWord#` 1## || b# `eqWord#` 1##)- || ((((a# `or#` b#) `and#` 1##) `eqWord#` 1##) -- not both even- && ((a# `neWord#` 0## && b# `neWord#` 0##) -- neither is zero- && gcdWordOdd# (shiftToOdd# a#) (shiftToOdd# b#) `eqWord#` 1##))+ (isTrue# (a# `eqWord#` 1##) || isTrue# (b# `eqWord#` 1##))+ || (isTrue# (((a# `or#` b#) `and#` 1##) `eqWord#` 1##) -- not both even+ && ((isTrue# (a# `neWord#` 0##) && isTrue# (b# `neWord#` 0##)) -- neither is zero+ && isTrue# (gcdWordOdd# (shiftToOdd# a#) (shiftToOdd# b#) `eqWord#` 1##))) -- Various auxiliary functions @@ -82,10 +83,10 @@ {-# INLINE gcdWordOdd# #-} gcdWordOdd# :: Word# -> Word# -> Word# gcdWordOdd# a# b#- | a# `eqWord#` 1## || b# `eqWord#` 1## = 1##- | a# `eqWord#` b# = a#- | a# `ltWord#` b# = oddGCD# b# a#- | otherwise = oddGCD# a# b#+ | isTrue# (a# `eqWord#` 1##) || isTrue# (b# `eqWord#` 1##) = 1##+ | isTrue# (a# `eqWord#` b#) = a#+ | isTrue# (a# `ltWord#` b#) = oddGCD# b# a#+ | otherwise = oddGCD# a# b# -- calculate the gcd of two odd numbers using the binary gcd algorithm -- Precondition: first argument strictly larger than second (which should be greater than 1)@@ -93,11 +94,11 @@ oddGCD# a# b# = case shiftToOdd# (a# `minusWord#` b#) of 1## -> 1##- c# | c# `ltWord#` b# -> oddGCD# b# c#- | c# `gtWord#` b# -> oddGCD# c# b#- | otherwise -> c#+ c# | isTrue# (c# `ltWord#` b#) -> oddGCD# b# c#+ | isTrue# (c# `gtWord#` b#) -> oddGCD# c# b#+ | otherwise -> c# absInt# :: Int# -> Int# absInt# i#- | i# <# 0# = negateInt# i#- | otherwise = i#+ | isTrue# (i# <# 0#) = negateInt# i#+ | otherwise = i#
Math/NumberTheory/Logarithms.hs view
@@ -24,6 +24,7 @@ ) where import GHC.Base+ #if __GLASGOW_HASKELL__ < 705 import GHC.Word (Word(..)) -- Moved to GHC.Types #endif@@ -34,6 +35,9 @@ import Math.NumberTheory.Logarithms.Internal import Math.NumberTheory.Powers.Integer+#if __GLASGOW_HASKELL__ < 707+import Math.NumberTheory.Utils (isTrue#)+#endif -- | Calculate the integer logarithm for an arbitrary base. -- The base must be greater than 1, the second argument, the number@@ -65,15 +69,15 @@ -- The argument must be positive, otherwise an error is thrown. intLog2 :: Int -> Int intLog2 (I# i#)- | i# <# 1# = error "Math.NumberTheory.Logarithms.intLog2: argument must be positive"- | otherwise = I# (wordLog2# (int2Word# 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#)- | w# `eqWord#` 0## = error "Math.NumberTheory.Logarithms.wordLog2: argument must not be 0."- | otherwise = I# (wordLog2# 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.
Math/NumberTheory/Moduli.hs view
@@ -162,37 +162,46 @@ -- calculates @(base ^ exponent) \`mod\` modulus@ by repeated squaring and reduction. -- If @exponent < 0@ and @base@ is invertible modulo @modulus@, @(inverse ^ |exponent|) \`mod\` modulus@ -- is calculated. This function does some input checking and sanitation before calling the unsafe worker.-{-# SPECIALISE powerMod :: Integer -> Int -> Integer -> Integer,- Integer -> Word -> Integer -> Integer- #-} {-# RULES "powerMod/Integer" powerMod = powerModInteger #-}+{-# INLINE [1] powerMod #-} powerMod :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer-powerMod base expo md+powerMod = powerModImpl++{-# SPECIALISE powerModImpl :: Integer -> Int -> Integer -> Integer,+ Integer -> Word -> Integer -> Integer+ #-}+powerModImpl :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer+powerModImpl base expo md | md == 0 = base ^ expo | md' == 1 = 0 | expo == 0 = 1 | bse' == 1 = 1 | expo < 0 = case invertMod bse' md' of- Just i -> powerMod' i (negate expo) md'+ Just i -> powerMod'Impl i (negate expo) md' Nothing -> error "Math.NumberTheory.Moduli.powerMod: Base isn't invertible with respect to modulus" | bse' == 0 = 0- | otherwise = powerMod' bse' expo md'+ | otherwise = powerMod'Impl bse' expo md' where md' = abs md bse' = if base < 0 || md' <= base then base `mod` md' else base -- | Modular power worker without input checking. -- Assumes all arguments strictly positive and modulus greater than 1.-{-# SPECIALISE powerMod' :: Integer -> Int -> Integer -> Integer,- Integer -> Word -> Integer -> Integer- #-} {-# RULES "powerMod'/Integer" powerMod' = powerModInteger' #-}+{-# INLINE [1] powerMod' #-} powerMod' :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer-powerMod' base expo md = go expo 1 base+powerMod' = powerMod'Impl+++{-# SPECIALISE powerMod'Impl :: Integer -> Int -> Integer -> Integer,+ Integer -> Word -> Integer -> Integer+ #-}+powerMod'Impl :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer+powerMod'Impl base expo md = go expo 1 base where go 1 !a !s = (a*s) `rem` md go e a s@@ -296,7 +305,7 @@ sqrtModP n 2 = Just (n `mod` 2) sqrtModP n prime = case jacobi' n prime of 0 -> Just 0- 1 -> Just (tonelliShanks (n `mod` prime) prime)+ 1 -> Just (sqrtModP' (n `mod` prime) prime) _ -> Nothing -- | @sqrtModPList n prime@ computes the list of all square roots of @n@@@ -358,16 +367,21 @@ Just r -> Just $ fixup r _ -> Nothing where- fixup r = case splitOff prime (r*r-n) of- (e,q) | expo <= e -> r- | otherwise -> hoist (fromJust $ invertMod (2*r) prime) r (q `mod` prime) (prime^e)+ fixup r = let diff' = r*r-n+ in if diff' == 0+ then r+ else case splitOff prime diff' of+ (e,q) | expo <= e -> r+ | otherwise -> hoist (fromJust $ invertMod (2*r) prime) r (q `mod` prime) (prime^e) -- hoist inv root elim pp+ | diff' == 0 = root' | expo <= ex = root' | otherwise = hoist inv root' (nelim `mod` prime) (prime^ex) where root' = (root + (inv*(prime-elim))*pp) `mod` (prime*pp)- (ex, nelim) = splitOff prime (root'*root' - n)+ diff' = root'*root' - n+ (ex, nelim) = splitOff prime diff' -- dirty, dirty sqM2P :: Integer -> Int -> Maybe Integer
Math/NumberTheory/Powers/Cubes.hs view
@@ -32,6 +32,9 @@ import GHC.Integer.GMP.Internals import Math.NumberTheory.Logarithms.Internal (integerLog2#)+#if __GLASGOW_HASKELL__ < 707+import Math.NumberTheory.Utils (isTrue#)+#endif -- | Calculate the integer cube root of an integer @n@, -- that is the largest integer @r@ such that @r^3 <= n@.@@ -58,8 +61,9 @@ {-# RULES "integerCubeRoot'/Int" integerCubeRoot' = cubeRootInt' "integerCubeRoot'/Word" integerCubeRoot' = cubeRootWord+"integerCubeRoot'/Igr" integerCubeRoot' = cubeRootIgr #-}-{-# SPECIALISE integerCubeRoot' :: Integer -> Integer #-}+{-# INLINE [1] integerCubeRoot' #-} integerCubeRoot' :: Integral a => a -> a integerCubeRoot' 0 = 0 integerCubeRoot' n = newton3 n (approxCuRt n)@@ -164,6 +168,10 @@ d = 3*r*(r+1) e = c+d +cubeRootIgr :: Integer -> Integer+cubeRootIgr 0 = 0+cubeRootIgr n = newton3 n (approxCuRt n)+ {-# SPECIALISE newton3 :: Int -> Int -> Int #-} {-# SPECIALISE newton3 :: Integer -> Integer -> Integer #-} newton3 :: Integral a => a -> a -> a@@ -194,7 +202,7 @@ appCuRt (S# i#) = case double2Int# (int2Double# i# **## (1.0## /## 3.0##)) of r# -> S# r# appCuRt n@(J# s# _)- | s# <# THRESH# = floor (fromInteger n ** (1.0/3.0) :: Double)+ | isTrue# (s# <# THRESH#) = floor (fromInteger n ** (1.0/3.0) :: Double) | otherwise = case integerLog2# n of l# -> case (l# `quotInt#` 3#) -# 51# of h# -> case shiftRInteger n (3# *# h#) of
Math/NumberTheory/Powers/Fourth.hs view
@@ -34,6 +34,9 @@ #endif import Math.NumberTheory.Logarithms.Internal (integerLog2#)+#if __GLASGOW_HASKELL__ < 707+import Math.NumberTheory.Utils (isTrue#)+#endif -- | Calculate the integer fourth root of a nonnegative number, -- that is, the largest integer @r@ with @r^4 <= n@.@@ -49,8 +52,9 @@ {-# RULES "integerFourthRoot'/Int" integerFourthRoot' = biSqrtInt "integerFourthRoot'/Word" integerFourthRoot' = biSqrtWord+"integerFourthRoot'/Igr" integerFourthRoot' = biSqrtIgr #-}-{-# SPECIALISE integerFourthRoot' :: Integer -> Integer #-}+{-# INLINE [1] integerFourthRoot' #-} integerFourthRoot' :: Integral a => a -> a integerFourthRoot' 0 = 0 integerFourthRoot' n = newton4 n (approxBiSqrt n)@@ -137,7 +141,7 @@ appBiSqrt :: Integer -> Integer appBiSqrt (S# i#) = S# (double2Int# (sqrtDouble# (sqrtDouble# (int2Double# i#)))) appBiSqrt n@(J# s# _)- | s# <# THRESH# = floor (sqrt . sqrt $ fromInteger n :: Double)+ | isTrue# (s# <# THRESH#) = floor (sqrt . sqrt $ fromInteger n :: Double) | otherwise = case integerLog2# n of l# -> case uncheckedIShiftRA# l# 2# -# 47# of h# -> case shiftRInteger n (4# *# h#) of@@ -203,3 +207,6 @@ r2 = r*r r4 = r2*r2 +biSqrtIgr :: Integer -> Integer+biSqrtIgr 0 = 0+biSqrtIgr n = newton4 n (approxBiSqrt n)
Math/NumberTheory/Powers/General.hs view
@@ -17,7 +17,7 @@ , isKthPower , isPerfectPower , highestPower--- , largePFPower+ , largePFPower ) where #include "MachDeps.h"@@ -33,9 +33,14 @@ import Data.List (foldl') import qualified Data.Set as Set --- import Math.NumberTheory.Logarithms+import Math.NumberTheory.Logarithms (integerLogBase') import Math.NumberTheory.Logarithms.Internal (integerLog2#)-import Math.NumberTheory.Utils (shiftToOddCount, splitOff)+import Math.NumberTheory.Utils (shiftToOddCount+ , splitOff+#if __GLASGOW_HASKELL__ < 707+ , isTrue#+#endif+ ) import qualified Math.NumberTheory.Powers.Squares as P2 import qualified Math.NumberTheory.Powers.Cubes as P3 import qualified Math.NumberTheory.Powers.Fourth as P4@@ -169,20 +174,19 @@ sqr 0 m = m sqr k m = sqr (k-1) (m*m) --- Not used, at least not yet--- -- | @'largePFPower' bd n@ produces the pair @(b,k)@ with the largest--- -- exponent @k@ such that @n == b^k@, where @bd > 1@ (it is expected--- -- that @bd@ is much larger, at least @1000@ or so), @n > bd^2@ and @n@--- -- has no prime factors @p <= bd@, skipping the trial division phase--- -- of @'highestPower'@ when that is a priori known to be superfluous.--- -- It is only present to avoid duplication of work in factorisation--- -- and primality testing, it is not expected to be generally useful.--- -- The assumptions are not checked, if they are not satisfied, wrong--- -- results and wasted work may be the consequence.--- largePFPower :: Integer -> Integer -> (Integer, Int)--- largePFPower bd n = rawPower ln n--- where--- ln = integerLogBase' (bd+1) n+-- | @'largePFPower' bd n@ produces the pair @(b,k)@ with the largest+-- exponent @k@ such that @n == b^k@, where @bd > 1@ (it is expected+-- that @bd@ is much larger, at least @1000@ or so), @n > bd^2@ and @n@+-- has no prime factors @p <= bd@, skipping the trial division phase+-- of @'highestPower'@ when that is a priori known to be superfluous.+-- It is only present to avoid duplication of work in factorisation+-- and primality testing, it is not expected to be generally useful.+-- The assumptions are not checked, if they are not satisfied, wrong+-- results and wasted work may be the consequence.+largePFPower :: Integer -> Integer -> (Integer, Int)+largePFPower bd n = rawPower ln n+ where+ ln = integerLogBase' (bd+1) n ------------------------------------------------------------------------------------------ -- Auxiliary functions --@@ -220,7 +224,7 @@ 1# -> 3 2# -> 5 3# -> 11- h# | h# <# 500# ->+ h# | isTrue# (h# <# 500#) -> floor (scaleFloat (I# (h# -# 1#)) (fromInteger (n `shiftRInteger` (h# *# k# -# k#)) ** (1/fromIntegral k) :: Double)) | otherwise ->@@ -279,7 +283,10 @@ rawPower :: Int -> Integer -> (Integer, Int) rawPower mx n- | mx < 2 = (n,1)+ | mx < 2 = (n,1)+ | mx == 2 = case P2.exactSquareRoot n of+ Just r -> (r,2)+ Nothing -> (n,1) rawPower mx n = case P4.exactFourthRoot n of Just r -> case rawPower (mx `quot` 4) r of (m,e) -> (m, 4*e)
Math/NumberTheory/Powers/Integer.hs view
@@ -8,7 +8,7 @@ -- -- Slightly faster power function for Integer base and Int exponent. ---{-# LANGUAGE MagicHash, BangPatterns #-}+{-# LANGUAGE MagicHash, BangPatterns, CPP #-} {-# OPTIONS_HADDOCK hide #-} module Math.NumberTheory.Powers.Integer ( integerPower@@ -17,6 +17,9 @@ import GHC.Base import Math.NumberTheory.Logarithms.Internal ( wordLog2# )+#if __GLASGOW_HASKELL__ < 707+import Math.NumberTheory.Utils (isTrue#)+#endif -- | Power of an 'Integer' by the left-to-right repeated squaring algorithm. -- This needs two multiplications in each step while the right-to-left@@ -25,16 +28,16 @@ -- gains a bit. integerPower :: Integer -> Int -> Integer integerPower b (I# e#)- | e# ==# 0# = 1- | e# ==# 1# = b- | otherwise = go (wordLog2# w# -# 1#) b (b*b)+ | isTrue# (e# ==# 0#) = 1+ | isTrue# (e# ==# 1#) = b+ | otherwise = go (wordLog2# w# -# 1#) b (b*b) where !w# = int2Word# e#- go 0# l h = if (w# `and#` 1##) `eqWord#` 0## then l*l else (l*h)+ 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# = ((w# `uncheckedShiftRL#` i#) `and#` 1##) `neWord#` 0##+hasBit# w# i# = isTrue# (((w# `uncheckedShiftRL#` i#) `and#` 1##) `neWord#` 0##)
Math/NumberTheory/Powers/Squares.hs view
@@ -37,6 +37,9 @@ #endif import Math.NumberTheory.Logarithms.Internal (integerLog2#)+#if __GLASGOW_HASKELL__ < 707+import Math.NumberTheory.Utils (isTrue#)+#endif -- | Calculate the integer square root of a nonnegative number @n@,@@ -58,7 +61,7 @@ "integerSquareRoot'/Int" integerSquareRoot' = isqrtInt' "integerSquareRoot'/Word" integerSquareRoot' = isqrtWord #-}-{-# SPECIALISE integerSquareRoot' :: Integer -> Integer #-}+{-# INLINE [1] integerSquareRoot' #-} integerSquareRoot' :: Integral a => a -> a integerSquareRoot' = isqrtA @@ -191,7 +194,7 @@ appSqrt :: Integer -> Integer appSqrt (S# i#) = S# (double2Int# (sqrtDouble# (int2Double# i#))) appSqrt n@(J# s# _)- | s# <# THRESH# = floor (sqrt $ fromInteger n :: Double)+ | isTrue# (s# <# THRESH#) = floor (sqrt $ fromInteger n :: Double) | otherwise = case integerLog2# n of l# -> case uncheckedIShiftRA# l# 1# -# 47# of h# -> case shiftRInteger n (2# *# h#) of
Math/NumberTheory/Primes/Counting/Impl.hs view
@@ -25,7 +25,10 @@ import Math.NumberTheory.Logarithms import Data.Array.Base-import Data.Array.ST hiding (unsafeThaw)+import Data.Array.ST+#if !MIN_VERSION_array(0,5,0)+ hiding (unsafeThaw)+#endif import Control.Monad.ST import Data.Bits import Data.Int
Math/NumberTheory/Primes/Factorisation/Montgomery.hs view
@@ -56,6 +56,8 @@ import Math.NumberTheory.Logarithms import Math.NumberTheory.Logarithms.Internal+import Math.NumberTheory.Powers.General (highestPower, largePFPower)+import Math.NumberTheory.Powers.Squares (integerSquareRoot') import Math.NumberTheory.Primes.Sieve.Eratosthenes import Math.NumberTheory.Primes.Sieve.Indexing import Math.NumberTheory.Primes.Testing.Probabilistic@@ -158,6 +160,9 @@ Just bd -> \k -> k <= bd || primeTest k Nothing -> primeTest rndR k = state (\gen -> prng k gen)+ perfPw = case primeBound of+ Nothing -> highestPower+ Just bd -> largePFPower (integerSquareRoot' bd) fact m digs = do let (b1,b2,ct) = findParms digs (pfs,cfs) <- repFact m b1 b2 ct if null cfs@@ -166,12 +171,19 @@ nfs <- forM cfs $ \(k,j) -> mult j <$> fact k (if null pfs then digs+4 else digs) return (mergeAll $ pfs:nfs)- repFact m b1 b2 count+ repFact m b1 b2 count = case perfPw m of+ (_,1) -> workFact m b1 b2 count+ (b,e)+ | ptest b -> return ([(b,e)],[])+ | otherwise -> do+ (as,bs) <- workFact b b1 b2 count+ return $ (mult e as, mult e bs)+ workFact m b1 b2 count | count < 0 = return ([],[(m,1)]) | otherwise = do s <- rndR m case montgomeryFactorisation m b1 b2 s of- Nothing -> repFact m b1 b2 (count-1)+ Nothing -> workFact m b1 b2 (count-1) Just d -> do let !cof = m `quot` d case gcd cof d of
Math/NumberTheory/Primes/Sieve/Eratosthenes.hs view
@@ -36,7 +36,10 @@ import Control.Monad.ST import Data.Array.Base-import Data.Array.ST hiding (unsafeFreeze, unsafeThaw, castSTUArray)+import Data.Array.ST+#if !MIN_VERSION_array(0,5,0)+ hiding (unsafeFreeze, unsafeThaw, castSTUArray)+#endif import Control.Monad (when) import Data.Bits import Data.Word
Math/NumberTheory/Primes/Testing/Probabilistic.hs view
@@ -27,6 +27,7 @@ import Data.Bits import GHC.Base+ #if __GLASGOW_HASKELL__ < 705 import GHC.Word -- Moved to GHC.Types #endif@@ -203,7 +204,7 @@ | testBit w i = go j# w (i - 1) 1 1 1 q | otherwise = look j# w (i-1) go k# w i un un1 vn qn- | i < 0 = if k# <# 0#+ | i < 0 = if isTrue# (k# <# 0#) then (un,vn,qn) else go (k# -# 1#) (W# (indexWordArray# ba# k#)) (WORD_SIZE_IN_BITS - 1) un un1 vn qn | testBit w i = go k# w (i-1) u2n1 u2n2 v2n1 q2n1
Math/NumberTheory/Utils.hs view
@@ -20,6 +20,9 @@ , bitCountWord# , uncheckedShiftR , splitOff+#if __GLASGOW_HASKELL__ < 707+ , isTrue#+#endif ) where #include "MachDeps.h"@@ -59,6 +62,7 @@ "shiftToOddCount/Word" shiftToOddCount = shiftOCWord "shiftToOddCount/Integer" shiftToOddCount = shiftOCInteger #-}+{-# INLINE [1] shiftToOddCount #-} shiftToOddCount :: (Integral a, Bits a) => a -> (Int, a) shiftToOddCount n = case shiftOCInteger (fromIntegral n) of (z, o) -> (z, fromInteger o)@@ -81,7 +85,7 @@ shiftOCInteger n@(S# i#) = case shiftToOddCount# (int2Word# i#) of (# z#, w# #)- | z# ==# 0# -> (0, n)+ | isTrue# (z# ==# 0#) -> (0, n) | otherwise -> (I# z#, S# (word2Int# w#)) shiftOCInteger n@(J# _ ba#) = case count 0# 0# of 0# -> (0, n)@@ -100,6 +104,7 @@ "shiftToOdd/Word" shiftToOdd = shiftOWord "shiftToOdd/Integer" shiftToOdd = shiftOInteger #-}+{-# INLINE [1] shiftToOdd #-} shiftToOdd :: (Integral a, Bits a) => a -> a shiftToOdd n = fromInteger (shiftOInteger (fromIntegral n)) @@ -191,3 +196,9 @@ go !k m = case m `quotRem` p of (q,r) | r == 0 -> go (k+1) q | otherwise -> (k,m)++#if __GLASGOW_HASKELL__ < 707+-- The times they are a-changing. The types of primops too :(+isTrue# :: Bool -> Bool+isTrue# = id+#endif
arithmoi.cabal view
@@ -1,5 +1,5 @@ name : arithmoi-version : 0.4.0.3+version : 0.4.0.4 cabal-version : >= 1.6 author : Daniel Fischer copyright : (c) 2011 Daniel Fischer@@ -27,13 +27,12 @@ category : Math, Algorithms, Number Theory -tested-with : GHC == 6.12.3, GHC == 7.0.2, GHC == 7.0.4, GHC == 7.2.1, GHC == 7.4.1,- GHC == 7.4.2+tested-with : GHC == 7.4.2, GHC==7.6.3, GHC==7.8.1 extra-source-files : Changes, TODO library- build-depends : base >= 4 && < 5, array >= 0.3 && < 0.5, ghc-prim,+ build-depends : base >= 4 && < 5, array >= 0.3 && < 0.6, ghc-prim, integer-gmp, containers >= 0.3 && < 0.6, random >= 1.0 && < 1.1, mtl >= 2.0 && < 2.2 exposed-modules : Math.NumberTheory.Logarithms