arithmoi (empty) → 0.1.0.0
raw patch · 32 files changed
+5190/−0 lines, 32 filesdep +arraydep +basedep +containerssetup-changed
Dependencies added: array, base, containers, ghc-prim, integer-gmp, mtl, random
Files
- Changes +2/−0
- LICENCE +7/−0
- Math/NumberTheory/GCD.hs +233/−0
- Math/NumberTheory/GCD/LowLevel.hs +101/−0
- Math/NumberTheory/Logarithms.hs +164/−0
- Math/NumberTheory/Logarithms/Internal.hs +155/−0
- Math/NumberTheory/Lucas.hs +99/−0
- Math/NumberTheory/Moduli.hs +283/−0
- Math/NumberTheory/Powers.hs +43/−0
- Math/NumberTheory/Powers/Cubes.hs +241/−0
- Math/NumberTheory/Powers/Fourth.hs +203/−0
- Math/NumberTheory/Powers/General.hs +338/−0
- Math/NumberTheory/Powers/Integer.hs +40/−0
- Math/NumberTheory/Powers/Squares.hs +272/−0
- Math/NumberTheory/Primes.hs +19/−0
- Math/NumberTheory/Primes/Counting.hs +21/−0
- Math/NumberTheory/Primes/Counting/Approximate.hs +62/−0
- Math/NumberTheory/Primes/Counting/Impl.hs +383/−0
- Math/NumberTheory/Primes/Factorisation.hs +166/−0
- Math/NumberTheory/Primes/Factorisation/Montgomery.hs +399/−0
- Math/NumberTheory/Primes/Factorisation/Utils.hs +77/−0
- Math/NumberTheory/Primes/Heap.hs +388/−0
- Math/NumberTheory/Primes/Sieve.hs +28/−0
- Math/NumberTheory/Primes/Sieve/Eratosthenes.hs +450/−0
- Math/NumberTheory/Primes/Sieve/Indexing.hs +140/−0
- Math/NumberTheory/Primes/Sieve/Misc.hs +344/−0
- Math/NumberTheory/Primes/Testing.hs +40/−0
- Math/NumberTheory/Primes/Testing/Probabilistic.hs +224/−0
- Math/NumberTheory/Utils.hs +187/−0
- Setup.hs +5/−0
- TODO +5/−0
- arithmoi.cabal +71/−0
+ Changes view
@@ -0,0 +1,2 @@+0.1.0:+First release
+ LICENCE view
@@ -0,0 +1,7 @@+Copyright (c) 2011 Daniel Fischer++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.
+ Math/NumberTheory/GCD.hs view
@@ -0,0 +1,233 @@+-- |+-- Module: Math.NumberTheory.GCD+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- This module exports GCD and coprimality test using the binary gcd algorithm+-- and GCD with the extended Euclidean algorithm.+--+-- Efficiently counting the number of trailing zeros, the binary gcd algorithm+-- can perform considerably faster than the Euclidean algorithm on average.+-- For 'Int', GHC has a rewrite rule to use GMP's fast gcd, depending on+-- hardware and\/or GMP version, that can be faster or slower than the binary+-- algorithm (on my 32-bit box, binary is faster, on my 64-bit box, GMP).+-- For 'Word' and the sized @IntN\/WordN@ types, there is no rewrite rule (yet)+-- in GHC, and the binary algorithm performs consistently (so far as my tests go)+-- much better (if this module's rewrite rules fire).+--+-- When using this module, always compile with optimisations turned on to+-- benefit from GHC's primops and the rewrite rules.+{-# LANGUAGE CPP, BangPatterns, MagicHash #-}+module Math.NumberTheory.GCD+ ( binaryGCD+ , extendedGCD+ , coprime+ ) where++import Data.Bits+import GHC.Word+import GHC.Int++import Math.NumberTheory.GCD.LowLevel+import Math.NumberTheory.Utils++#include "MachDeps.h"++{-# RULES+"binaryGCD/Int" binaryGCD = gcdInt+"binaryGCD/Word" binaryGCD = gcdWord+"binaryGCD/Int8" binaryGCD = gi8+"binaryGCD/Int16" binaryGCD = gi16+"binaryGCD/Int32" binaryGCD = gi32+"binaryGCD/Word8" binaryGCD = gw8+"binaryGCD/Word16" binaryGCD = gw16+"binaryGCD/Word32" binaryGCD = gw32+ #-}+#if WORD_SIZE_IN_BITS == 64+gi64 :: Int64 -> Int64 -> Int64+gi64 (I64# x#) (I64# y#) = I64# (gcdInt# x# y#)++gw64 :: Word64 -> Word64 -> Word64+gw64 (W64# x#) (W64# y#) = W64# (gcdWord# x# y#)++{-# RULES+"binaryGCD/Int64" binaryGCD = gi64+"binaryGCD/Word64" binaryGCD = gw64+ #-}+#else+{-# SPECIALISE binaryGCD :: Word64 -> Word64 -> Word64,+ Int64 -> Int64 -> Int64 #-}+#endif+{-# SPECIALISE binaryGCD :: Integer -> Integer -> Integer #-}+-- | 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.+--+-- There are specialised functions for @'Int'@ and @'Word'@ and rewrite rules+-- for those and @IntN@ and @WordN@, @N <= WORD_SIZE_IN_BITS@, to use the+-- specialised variants. These types are worth benchmarking, others probably not.+--+-- It is very slow for 'Integer' (and probably every type except the abovementioned),+-- I recommend not using it for those.+--+-- 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 =+ case shiftToOddCount a of+ (!za, !oa) ->+ case shiftToOddCount b of+ (!zb, !ob) -> gcdOdd (abs oa) (abs ob) `shiftL` min za zb++{-# SPECIALISE extendedGCD :: Int -> Int -> (Int, Int, Int),+ Word -> Word -> (Word, Word, Word),+ Integer -> Integer -> (Integer, Integer, Integer)+ #-}+-- | Calculate the greatest common divisor of two numbers and coefficients+-- for the linear combination.+--+-- Satisfies:+--+-- > case extendedGCD a b of+-- > (d, u, v) -> u*a + v*b == d+-- >+-- > d == gcd a b+--+-- and, for signed types,+--+-- >+-- > abs u < abs b || abs b <= 1+-- >+-- > abs v < abs a || abs a <= 1+--+-- (except if one of @a@ and @b@ is 'minBound' of a signed type).+extendedGCD :: Integral a => a -> a -> (a, a, a)+extendedGCD a b = (d, u, v)+ where+ (d, x, y) = eGCD 0 1 1 0 (abs a) (abs b)+ u | a < 0 = negate x+ | otherwise = x+ v | b < 0 = negate y+ | otherwise = y+ eGCD !n1 o1 !n2 o2 r s+ | s == 0 = (r, o1, o2)+ | otherwise = case r `quotRem` s of+ (q, t) -> eGCD (o1 - q*n1) n1 (o2 - q*n2) n2 s t++{-# RULES+"coprime/Int" coprime = coprimeInt+"coprime/Word" coprime = coprimeWord+"coprime/Int8" coprime = ci8+"coprime/Int16" coprime = ci16+"coprime/Int32" coprime = ci32+"coprime/Word8" coprime = cw8+"coprime/Word16" coprime = cw16+"coprime/Word32" coprime = cw32+ #-}+#if WORD_SIZE_IN_BITS == 64+ci64 :: Int64 -> Int64 -> Bool+ci64 (I64# x#) (I64# y#) = coprimeInt# x# y#++cw64 :: Word64 -> Word64 -> Bool+cw64 (W64# x#) (W64# y#) = coprimeWord# x# y#++{-# RULES+"coprime/Int64" coprime = ci64+"coprime/Word64" coprime = cw64+ #-}+#else+{-# SPECIALISE coprime :: Word64 -> Word64 -> Bool,+ Int64 -> Int64 -> Bool #-}+#endif+{-# SPECIALISE coprime :: Integer -> Integer -> Bool #-}+-- | 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.+--+-- The remarks about performance at 'binaryGCD' apply here too, use this function+-- only at the types with rewrite rules.+--+-- Relies on twos complement or sign and magnitude representaion for signed types.+coprime :: (Integral a, Bits a) => a -> a -> Bool+coprime a b =+ (a' == 1 || b' == 1)+ || (a' /= 0 && b' /= 0 && ((a .|. b) .&. 1) == 1+ && gcdOdd (abs (shiftToOdd a')) (abs (shiftToOdd b')) == 1)+ where+ a' = abs a+ b' = abs b++-- Auxiliaries++-- gcd of two odd numbers+{-# SPECIALISE gcdOdd :: Integer -> Integer -> Integer #-}+#if WORD_SIZE_IN_BITS < 64+{-# SPECIALISE gcdOdd :: Int64 -> Int64 -> Int64,+ Word64 -> Word64 -> Word64+ #-}+#endif+{-# INLINE gcdOdd #-}+gcdOdd :: (Integral a, Bits a) => a -> a -> a+gcdOdd a b+ | a == 1 || b == 1 = 1+ | a < b = oddGCD b a+ | a > b = oddGCD a b+ | otherwise = a++{-# SPECIALISE oddGCD :: Integer -> Integer -> Integer #-}+#if WORD_SIZE_IN_BITS < 64+{-# SPECIALISE oddGCD :: Int64 -> Int64 -> Int64,+ Word64 -> Word64 -> Word64+ #-}+#endif+oddGCD :: (Integral a, Bits a) => a -> a -> a+oddGCD a b =+ case shiftToOdd (a-b) of+ 1 -> 1+ c | c < b -> oddGCD b c+ | c > b -> oddGCD c b+ | otherwise -> c++-------------------------------------------------------------------------------+-- Blech! Getting the rules to fire isn't easy. --+-------------------------------------------------------------------------------++gi8 :: Int8 -> Int8 -> Int8+gi8 (I8# x#) (I8# y#) = I8# (gcdInt# x# y#)++gi16 :: Int16 -> Int16 -> Int16+gi16 (I16# x#) (I16# y#) = I16# (gcdInt# x# y#)++gi32 :: Int32 -> Int32 -> Int32+gi32 (I32# x#) (I32# y#) = I32# (gcdInt# x# y#)++gw8 :: Word8 -> Word8 -> Word8+gw8 (W8# x#) (W8# y#) = W8# (gcdWord# x# y#)++gw16 :: Word16 -> Word16 -> Word16+gw16 (W16# x#) (W16# y#) = W16# (gcdWord# x# y#)++gw32 :: Word32 -> Word32 -> Word32+gw32 (W32# x#) (W32# y#) = W32# (gcdWord# x# y#)++ci8 :: Int8 -> Int8 -> Bool+ci8 (I8# x#) (I8# y#) = coprimeInt# x# y#++ci16 :: Int16 -> Int16 -> Bool+ci16 (I16# x#) (I16# y#) = coprimeInt# x# y#++ci32 :: Int32 -> Int32 -> Bool+ci32 (I32# x#) (I32# y#) = coprimeInt# x# y#++cw8 :: Word8 -> Word8 -> Bool+cw8 (W8# x#) (W8# y#) = coprimeWord# x# y#++cw16 :: Word16 -> Word16 -> Bool+cw16 (W16# x#) (W16# y#) = coprimeWord# x# y#++cw32 :: Word32 -> Word32 -> Bool+cw32 (W32# x#) (W32# y#) = coprimeWord# x# y#
+ Math/NumberTheory/GCD/LowLevel.hs view
@@ -0,0 +1,101 @@+-- |+-- Module: Math.NumberTheory.GCD.LowLevel+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Low level gcd and coprimality functions using the binary gcd algorithm.+-- Normally, accessing these via the higher level interface of "Math.NumberTheory.GCD"+-- should be sufficient.+--+{-# LANGUAGE MagicHash, UnboxedTuples #-}+module Math.NumberTheory.GCD.LowLevel+ ( -- * Specialised GCDs+ gcdInt+ , gcdWord+ -- ** GCDs for unboxed types+ , gcdInt#+ , gcdWord#+ -- * Specialised tests for coprimality+ , coprimeInt+ , coprimeWord+ -- ** Coprimality tests for unboxed types+ , coprimeInt#+ , coprimeWord#+ ) where++import GHC.Base+import GHC.Word (Word(..))++import Math.NumberTheory.Utils++-- | Greatest common divisor of two 'Int's, calculated with the binary gcd algorithm.+gcdInt :: Int -> Int -> Int+gcdInt (I# a#) (I# b#) = I# (gcdInt# a# b#)++-- | Test whether two 'Int's are coprime, using an abbreviated binary gcd algorithm.+coprimeInt :: Int -> Int -> Bool+coprimeInt (I# a#) (I# b#) = coprimeInt# a# b#++-- | Greatest common divisor of two 'Word's, calculated with the binary gcd algorithm.+gcdWord :: Word -> Word -> Word+gcdWord (W# a#) (W# b#) = W# (gcdWord# a# b#)++-- | Test whether two 'Word's are coprime, using an abbreviated binary gcd algorithm.+coprimeWord :: Word -> Word -> Bool+coprimeWord (W# a#) (W# b#) = coprimeWord# a# b#++-- | Greatest common divisor of two 'Int#'s, calculated with the binary gcd algorithm.+gcdInt# :: Int# -> Int# -> Int#+gcdInt# a# b# = word2Int# (gcdWord# (int2Word# (absInt# a#)) (int2Word# (absInt# b#)))+++-- | Test whether two 'Int#'s are coprime.+coprimeInt# :: Int# -> Int# -> Bool+coprimeInt# a# b# = coprimeWord# (int2Word# (absInt# a#)) (int2Word# (absInt# b#))++-- | Greatest common divisor of two 'Word#'s, calculated with the binary gcd algorithm.+gcdWord# :: Word# -> Word# -> Word#+gcdWord# a# 0## = a#+gcdWord# 0## b# = b#+gcdWord# a# b# =+ case shiftToOddCount# a# of+ (# za#, oa# #) ->+ case shiftToOddCount# b# of+ (# zb#, ob# #) -> gcdWordOdd# oa# ob# `uncheckedShiftL#` (if 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##))++-- Various auxiliary functions++-- calculate the gcd of two odd numbers+{-# 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#++-- 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)+oddGCD# :: Word# -> Word# -> Word#+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#++absInt# :: Int# -> Int#+absInt# i#+ | i# <# 0# = negateInt# i#+ | otherwise = i#
+ Math/NumberTheory/Logarithms.hs view
@@ -0,0 +1,164 @@+-- |+-- Module: Math.NumberTheory.Logarithms+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Integer Logarithms. For efficiency, the internal representation of 'Integer's+-- from integer-gmp is used.+--+{-# LANGUAGE MagicHash, UnboxedTuples #-}+module Math.NumberTheory.Logarithms+ ( -- * Integer logarithms with input checks+ integerLogBase+ , integerLog2+ , intLog2+ , wordLog2+ -- * Integer logarithms without input checks+ , integerLogBase'+ , integerLog2'+ , intLog2'+ , wordLog2'+ ) where++import GHC.Base+import GHC.Word (Word(..))++import Data.Bits+import Data.Array.Unboxed+import Data.Array.Base (unsafeAt)++import Math.NumberTheory.Logarithms.Internal+import Math.NumberTheory.Powers.Integer++-- | 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)++-- | 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#)+ | 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#)++-- | 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 '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#)++-- | 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++-- 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+ ]
+ Math/NumberTheory/Logarithms/Internal.hs view
@@ -0,0 +1,155 @@+-- |+-- Module: Math.NumberTheory.Logarithms.Internal+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Low level stuff for integer logarithms.+{-# LANGUAGE CPP, MagicHash, UnboxedTuples #-}+{-# OPTIONS_HADDOCK hide #-}+module Math.NumberTheory.Logarithms.Internal+ ( -- * Functions+ 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 `ltInteger` pw+ then (# m, 0# #)+ else case step (pw `timesInteger` pw) of+ (# q, e #) ->+ if q `ltInteger` pw+ then (# q, 2# *# e #)+ else (# q `quotInteger` 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
+ Math/NumberTheory/Lucas.hs view
@@ -0,0 +1,99 @@+-- |+-- Module: Math.NumberTheory.Lucas+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Efficient calculation of Lucas sequences.+{-# LANGUAGE CPP #-}+module Math.NumberTheory.Lucas+ ( fibonacci+ , fibonacciPair+ , lucas+ , lucasPair+ , generalLucas+ ) where++#include "MachDeps.h"++import Data.Bits++-- | @'fibonacci' k@ calculates the @k@-th Fibonacci number in+-- /O/(@log (abs k)@) steps. The index may be negative. This+-- is efficient for calculating single Fibonacci numbers (with+-- large index), but for computing many Fibonacci numbers in+-- close proximity, it is better to use the simple addition+-- formula starting from an appropriate pair of successive+-- Fibonacci numbers.+fibonacci :: Int -> Integer+fibonacci = fst . fibonacciPair++-- | @'fibonacciPair' k@ returns the pair @(F(k), F(k+1))@ of the @k@-th+-- Fibonacci number and its successor, thus it can be used to calculate+-- the Fibonacci numbers from some index on without needing to compute+-- the previous. The pair is efficiently calculated+-- in /O/(@log (abs k)@) steps. The index may be negative.+fibonacciPair :: Int -> (Integer, Integer)+fibonacciPair n+ | n < 0 = let (f,g) = fibonacciPair (-(n+1)) in if testBit n 0 then (g, -f) else (-g, f)+ | n == 0 = (0, 1)+ | otherwise = look (WORD_SIZE_IN_BITS - 2)+ where+ look k+ | testBit n k = go (k-1) 0 1+ | otherwise = look (k-1)+ go k g f+ | k < 0 = (f, f+g)+ | testBit n k = go (k-1) (f*(f+shiftL g 1)) ((f+g)*shiftL f 1 + g*g)+ | otherwise = go (k-1) (f*f+g*g) (f*(f+shiftL g 1))++-- | @'lucas' k@ computes the @k@-th Lucas number. Very similar+-- to @'fibonacci'@.+lucas :: Int -> Integer+lucas = fst . lucasPair++-- | @'lucasPair' k@ computes the pair @(L(k), L(k+1))@ of the @k@-th+-- Lucas number and its successor. Very similar to @'fibonacciPair'@.+lucasPair :: Int -> (Integer, Integer)+lucasPair n+ | n < 0 = let (f,g) = lucasPair (-(n+1)) in if testBit n 0 then (-g, f) else (g, -f)+ | n == 0 = (2, 1)+ | otherwise = look (WORD_SIZE_IN_BITS - 2)+ where+ look k+ | testBit n k = go (k-1) 0 1+ | otherwise = look (k-1)+ go k g f+ | k < 0 = (shiftL g 1 + f,g+3*f)+ | otherwise = go (k-1) g' f'+ where+ (f',g')+ | testBit n k = (shiftL (f*(f+g)) 1 + g*g,f*(shiftL g 1 + f))+ | otherwise = (f*(shiftL g 1 + f),f*f+g*g)+++-- | @'generalLucas' p q k@ calculates the quadruple @(U(k), U(k+1), V(k), V(k+1))@+-- where @U(i)@ is the Lucas sequence of the first kind and @V(i)@ the Lucas+-- sequence of the second kind for the parameters @p@ and @q@, where @p^2-4q /= 0@.+-- Both sequences satisfy the recurrence relation @A(j+2) = p*A(j+1) - q*A(j)@,+-- the starting values are @U(0) = 0, U(1) = 1@ and @V(0) = 2, V(1) = p@.+-- The Fibonacci numbers form the Lucas sequence of the first kind for the+-- parameters @p = 1, q = -1@ and the Lucas numbers form the Lucas sequence of+-- the second kind for these parameters.+-- Here, the index must be non-negative, since the terms of the sequence for+-- negative indices are in general not integers.+generalLucas :: Integer -> Integer -> Int -> (Integer, Integer, Integer, Integer)+generalLucas p q k+ | k < 0 = error "generalLucas: negative index"+ | k == 0 = (0,1,2,p)+ | otherwise = look (WORD_SIZE_IN_BITS - 2)+ where+ look i+ | testBit k i = go (i-1) 1 p p q+ | otherwise = look (i-1)+ go i un un1 vn qn+ | i < 0 = (un, un1, vn, p*un1 - shiftL (q*un) 1)+ | testBit k i = go (i-1) (un1*vn-qn) ((p*un1-q*un)*vn - p*qn) ((p*un1 - (2*q)*un)*vn - p*qn) (qn*qn*q)+ | otherwise = go (i-1) (un*vn) (un1*vn-qn) (vn*vn - 2*qn) (qn*qn)
+ Math/NumberTheory/Moduli.hs view
@@ -0,0 +1,283 @@+-- |+-- Module: Math.NumberTheory.Moduli+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Miscellaneous functions related to modular arithmetic.+--+{-# LANGUAGE CPP, BangPatterns #-}+module Math.NumberTheory.Moduli+ ( -- * Functions with input check+ jacobi+ , invertMod+ , powerMod+ , powerModInteger+ -- * Unchecked functions+ , jacobi'+ , powerMod'+ , powerModInteger'+ ) where++#include "MachDeps.h"++import Data.Word+import Data.Bits+import Data.Array.Unboxed+import Data.Array.Base (unsafeAt)++import Math.NumberTheory.GCD (extendedGCD)+import Math.NumberTheory.Utils (shiftToOddCount)++-- | Invert a number relative to a modulus.+-- If @number@ and @modulus@ are coprime, the result is+-- @Just inverse@ where+--+-- > (number * inverse) `mod` (abs modulus) == 1+-- > 0 <= inverse < abs modulus+--+-- unless @modulus == 0@ and @abs number == 1@, in which case the+-- result is @Just number@.+-- If @gcd number modulus > 1@, the result is @Nothing@.+invertMod :: Integer -> Integer -> Maybe Integer+invertMod k 0 = if k == 1 || k == (-1) then Just k else Nothing+invertMod k m = case extendedGCD k' m' of+ (1, u, _) -> Just (if u < 0 then m' + u else u)+ _ -> Nothing+ where+ m' = abs m+ k' | k >= m' || k < 0 = k `mod` m'+ | otherwise = k++-- | Jacobi symbol of two numbers.+-- The \"denominator\" must be odd and positive, this condition is checked.+--+-- If both numbers have a common prime factor, the result+-- is @0@, otherwise it is ±1.+{-# SPECIALISE jacobi :: Integer -> Integer -> Int,+ Int -> Int -> Int,+ Word -> Word -> Int+ #-}+jacobi :: (Integral a, Bits a) => a -> a -> Int+jacobi a b+ | b < 0 = error "Math.NumberTheory.Moduli.jacobi: negative denominator"+ | evenI b = error "Math.NumberTheory.Moduli.jacobi: even denominator"+ | b == 1 = 1+ | a == 0 = 0+ | a == 1 = 1+ | otherwise = jacobi' a b -- b odd, > 1, a neither 0 or 1++-- Invariant: b > 1 and odd+-- | Jacobi symbol of two numbers without validity check of+-- the \"denominator\".+{-# SPECIALISE jacobi' :: Integer -> Integer -> Int,+ Int -> Int -> Int,+ Word -> Word -> Int+ #-}+jacobi' :: (Integral a, Bits a) => a -> a -> Int+jacobi' a b+ | a == 0 = 0+ | a == 1 = 1+ | a < 0 = let n | rem4 b == 1 = 1+ | otherwise = -1+ -- Blech, minBound may pose problems+ (z,o) = shiftToOddCount (abs $ toInteger a)+ s | evenI z || unsafeAt jac2 (rem8 b) == 1 = n+ | otherwise = (-n)+ in s*jacobi' (fromInteger o) b+ | a >= b = case a `rem` b of+ 0 -> 0+ r -> jacPS 1 r b+ | evenI a = case shiftToOddCount a of+ (z,o) -> let r = 2 - (rem4 o .&. rem4 b)+ s | evenI z || unsafeAt jac2 (rem8 b) == 1 = r+ | otherwise = (-r)+ in jacOL s b o+ | otherwise = case rem4 a .&. rem4 b of+ 3 -> jacOL (-1) b a+ _ -> jacOL 1 b a++-- numerator positive and smaller than denominator+{-# SPECIALISE jacPS :: Int -> Integer -> Integer -> Int,+ Int -> Int -> Int -> Int,+ Int -> Word -> Word -> Int+ #-}+jacPS :: (Integral a, Bits a) => Int -> a -> a -> Int+jacPS !j a b+ | evenI a = case shiftToOddCount a of+ (z,o) | evenI z || unsafeAt jac2 (rem8 b) == 1 ->+ jacOL (if rem4 o .&. rem4 b == 3 then (-j) else j) b o+ | otherwise ->+ jacOL (if rem4 o .&. rem4 b == 3 then j else (-j)) b o+ | otherwise = jacOL (if rem4 a .&. rem4 b == 3 then (-j) else j) b a++-- numerator odd, positive and larger than denominator+{-# SPECIALISE jacOL :: Int -> Integer -> Integer -> Int,+ Int -> Int -> Int -> Int,+ Int -> Word -> Word -> Int+ #-}+jacOL :: (Integral a, Bits a) => Int -> a -> a -> Int+jacOL !j a b+ | b == 1 = j+ | otherwise = case a `rem` b of+ 0 -> 0+ r -> jacPS j r b++-- | Modular power.+--+-- > powerMod base exponent modulus+--+-- 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+ #-}+powerMod :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer+powerMod 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'+ Nothing -> error "Math.NumberTheory.Moduli.powerMod: Base isn't invertible with respect to modulus"+ | bse' == 0 = 0+ | otherwise = powerMod' 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'+ #-}+powerMod' :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer+powerMod' base expo md = go expo 1 base+ where+ go 1 !a !s = (a*s) `rem` md+ go e a s+ | testBit e 0 = go (e `shiftR` 1) a ((s*s) `rem` md)+ | otherwise = go (e `shiftR` 1) ((a*s) `rem` md) ((s*s) `rem` md)++-- | Specialised version of 'powerMod' for 'Integer' exponents.+-- Reduces the number of shifts of the exponent since shifting+-- large 'Integer's is expensive. Call this function directly+-- if you don't want or can't rely on rewrite rules.+powerModInteger :: Integer -> Integer -> Integer -> Integer+powerModInteger base ex mdl+ | mdl == 0 = base ^ ex+ | mdl' == 1 = 0+ | ex == 0 = 1+ | ex < 0 = case invertMod bse' mdl' of+ Just i -> powerModInteger' i (negate ex) mdl'+ Nothing -> error "Math.NumberTheory.Moduli.powerMod: Base isn't invertible with respect to modulus"+ | bse' == 0 = 0+ | bse' == 1 = 1+ | otherwise = powerModInteger' bse' ex mdl'+ where+ mdl' = abs mdl+ bse' = if base < 0 || mdl' <= base then base `mod` mdl' else base++-- | Specialised worker without input checks. Makes the same assumptions+-- as the general version.+powerModInteger' :: Integer -> Integer -> Integer -> Integer+powerModInteger' base expo md = go e1 w1 1 base+ where+ w1 = fromInteger expo+ e1 = expo `shiftR` 64+#if WORD_SIZE_IN_BITS == 32+ -- Shifting large Integers is expensive, hence we reduce the+ -- number of shifts by processing in 64-bit chunks.+ -- On 32-bit systems, every testBit on a Word64 would be a C-call,+ -- thus it is faster to split each Word64 into the constituent 32-bit+ -- Words and process those separately.+ -- The code becomes ugly, unfortunately.+ go :: Integer -> Word64 -> Integer -> Integer -> Integer+ go 0 !w !a !s = end w a s+ go e w a s = inner1 0 a s+ where+ wl :: Word+ !wl = fromIntegral w+ wh :: Word+ !wh = fromIntegral (w `shiftR` 32)+ inner1 32 !au !sq = inner2 0 au sq+ inner1 i au sq+ | testBit wl i = inner1 (i+1) ((au*sq) `rem` md) ((sq*sq) `rem` md)+ | otherwise = inner1 (i+1) au ((sq*sq) `rem` md)+ inner2 32 !au !sq = go (e `shiftR` 64) (fromInteger e) au sq+ inner2 i au sq+ | testBit wh i = inner2 (i+1) ((au*sq) `rem` md) ((sq*sq) `rem` md)+ | otherwise = inner2 (i+1) au ((sq*sq) `rem` md)+ end w !a !s+ | wh == 0 = fin wl a s+ | otherwise = innerE 0 a s+ where+ wl :: Word+ !wl = fromIntegral w+ wh :: Word+ !wh = fromIntegral (w `shiftR` 32)+ innerE 32 !au !sq = fin wh au sq+ innerE i au sq+ | testBit wl i = innerE (i+1) ((au*sq) `rem` md) ((sq*sq) `rem` md)+ | otherwise = innerE (i+1) au ((sq*sq) `rem` md)+ fin :: Word -> Integer -> Integer -> Integer+ fin 1 !a !s = (a*s) `rem` md+ fin w a s+ | testBit w 0 = fin (w `shiftR` 1) ((a*s) `rem` md) ((s*s) `rem` md)+ | otherwise = fin (w `shiftR` 1) a ((s*s) `rem` md)++#else+ -- WORD_SIZE_IN_BITS == 64, otherwise things wouldn't compile anyway+ -- Shorter code since we need not split each 64-bit word.+ go :: Integer -> Word -> Integer -> Integer -> Integer+ go 0 !w !a !s = end w a s+ go e w a s = inner 0 a s+ where+ inner 64 !au !sq = go (e `shiftR` 64) (fromInteger e) au sq+ inner i au sq+ | testBit w i = inner (i+1) ((au*sq) `rem` md) ((sq*sq) `rem` md)+ | otherwise = inner (i+1) au ((sq*sq) `rem` md)+ end 1 !a !s = (a*s) `rem` md+ end w a s+ | testBit w 0 = end (w `shiftR` 1) ((a*s) `rem` md) ((s*s) `rem` md)+ | otherwise = end (w `shiftR` 1) a ((s*s) `rem` md)++#endif++-- Utilities++-- For large Integers, going via Int is much faster than bit-fiddling+-- on the Integer, so we do that.+{-# SPECIALISE evenI :: Integer -> Bool,+ Int -> Bool,+ Word -> Bool+ #-}+evenI :: (Integral a, Bits a) => a -> Bool+evenI n = fromIntegral n .&. 1 == (0 :: Int)++{-# SPECIALISE rem4 :: Integer -> Int,+ Int -> Int,+ Word -> Int+ #-}+rem4 :: (Integral a, Bits a) => a -> Int+rem4 n = fromIntegral n .&. 3++{-# SPECIALISE rem8 :: Integer -> Int,+ Int -> Int,+ Word -> Int+ #-}+rem8 :: (Integral a, Bits a) => a -> Int+rem8 n = fromIntegral n .&. 7++jac2 :: UArray Int Int+jac2 = array (0,7) [(0,0),(1,1),(2,0),(3,-1),(4,0),(5,-1),(6,0),(7,1)]
+ Math/NumberTheory/Powers.hs view
@@ -0,0 +1,43 @@+-- |+-- Module: Math.NumberTheory.Powers+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Calculating integer roots, modular powers and related things.+-- This module reexports the most needed functions from the implementation+-- modules. The implementation modules provide some additional functions,+-- in particular some unsafe functions which omit some tests for performance+-- reasons.+--+module Math.NumberTheory.Powers+ ( -- * Integer Roots+ -- ** Square roots+ integerSquareRoot+ , isSquare+ , exactSquareRoot+ -- ** Cube roots+ , integerCubeRoot+ , isCube+ , exactCubeRoot+ -- ** Fourth roots+ , integerFourthRoot+ , isFourthPower+ , exactFourthRoot+ -- ** General roots+ , integerRoot+ , isKthPower+ , exactRoot+ , isPerfectPower+ , highestPower+ -- Modular powers+ , powerMod+ ) where++import Math.NumberTheory.Powers.Squares+import Math.NumberTheory.Powers.Cubes+import Math.NumberTheory.Powers.Fourth+import Math.NumberTheory.Powers.General+import Math.NumberTheory.Moduli
+ Math/NumberTheory/Powers/Cubes.hs view
@@ -0,0 +1,241 @@+-- |+-- Module: Math.NumberTheory.Powers.Cubes+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Functions dealing with cubes. Moderately efficient calculation of integer+-- cube roots and testing for cubeness.+{-# LANGUAGE MagicHash, BangPatterns, CPP #-}+module Math.NumberTheory.Powers.Cubes+ ( integerCubeRoot+ , integerCubeRoot'+ , exactCubeRoot+ , isCube+ , isCube'+ , isPossibleCube+ ) where++#include "MachDeps.h"++import Data.Array.Unboxed+import Data.Array.Base+import Data.Array.ST++import Data.Bits+import Data.Word++import GHC.Base+import GHC.Integer+import GHC.Integer.GMP.Internals++import Math.NumberTheory.Logarithms.Internal (integerLog2#)++-- | Calculate the integer cube root of an integer @n@,+-- that is the largest integer @r@ such that @r^3 <= n@.+-- Note that this is not symmetric about @0@, for example+-- @integerCubeRoot (-2) = (-2)@ while @integerCubeRoot 2 = 1@.+{-# SPECIALISE integerCubeRoot :: Int -> Int,+ Integer -> Integer,+ Word -> Word+ #-}+integerCubeRoot :: Integral a => a -> a+integerCubeRoot 0 = 0+integerCubeRoot n+ | n > 0 = integerCubeRoot' n+ | otherwise =+ let m = negate n+ r = if m < 0+ then negate . fromInteger $ integerCubeRoot' (negate $ fromIntegral n)+ else negate (integerCubeRoot' m)+ in if r*r*r == n then r else (r-1)++-- | Calculate the integer cube root of a nonnegative integer @n@,+-- that is, the largest integer @r@ such that @r^3 <= n@.+-- The precondition @n >= 0@ is not checked.+{-# RULES+"integerCubeRoot'/Int" integerCubeRoot' = cubeRootInt'+"integerCubeRoot'/Word" integerCubeRoot' = cubeRootWord+ #-}+{-# SPECIALISE integerCubeRoot' :: Integer -> Integer #-}+integerCubeRoot' :: Integral a => a -> a+integerCubeRoot' 0 = 0+integerCubeRoot' n = newton3 n (approxCuRt n)++-- | Returns @Nothing@ if the argument is not a cube,+-- @Just r@ if @n == r^3@.+{-# SPECIALISE exactCubeRoot :: Int -> Maybe Int,+ Word -> Maybe Word,+ Integer -> Maybe Integer+ #-}+exactCubeRoot :: Integral a => a -> Maybe a+exactCubeRoot 0 = Just 0+exactCubeRoot n+ | n < 0 =+ if m < 0+ then fmap (negate . fromInteger) $ exactCubeRoot (negate $ fromIntegral n)+ else fmap negate (exactCubeRoot m)+ | isPossibleCube n && r*r*r == n = Just r+ | otherwise = Nothing+ where+ m = negate n+ r = integerCubeRoot' n++-- | Test whether an integer is a cube.+{-# SPECIALISE isCube :: Int -> Bool,+ Integer -> Bool,+ Word -> Bool+ #-}+isCube :: Integral a => a -> Bool+isCube 0 = True+isCube n+ | n > 0 = isCube' n+ | m > 0 = isCube' m+ | otherwise = isCube' (negate (fromIntegral n) :: Integer)+ where+ m = negate n++-- | Test whether a nonnegative integer is a cube.+-- Before 'integerCubeRoot' is calculated, a few tests+-- of remainders modulo small primes weed out most non-cubes.+-- For testing many numbers, most of which aren't cubes,+-- this is much faster than @let r = cubeRoot n in r*r*r == n@.+-- The condition @n >= 0@ is /not/ checked.+{-# SPECIALISE isCube' :: Int -> Bool,+ Integer -> Bool,+ Word -> Bool+ #-}+isCube' :: Integral a => a -> Bool+isCube' !n = isPossibleCube n+ && (r*r*r == n)+ where+ r = integerCubeRoot' n++-- | Test whether a nonnegative number is possibly a cube.+-- Only about 0.08% of all numbers pass this test.+-- The precondition @n >= 0@ is /not/ checked.+{-# SPECIALISE isPossibleCube :: Int -> Bool,+ Integer -> Bool,+ Word -> Bool+ #-}+isPossibleCube :: Integral a => a -> Bool+isPossibleCube !n =+ unsafeAt cr512 (fromIntegral n .&. 511)+ && unsafeAt cubeRes837 (fromIntegral (n `rem` 837))+ && unsafeAt cubeRes637 (fromIntegral (n `rem` 637))+ && unsafeAt cubeRes703 (fromIntegral (n `rem` 703))++----------------------------------------------------------------------+-- Utility Functions --+----------------------------------------------------------------------++-- Special case for 'Int', a little faster.+-- For @n <= 2^64@, the truncated 'Double' is never+-- more than one off. Things might overflow for @n@+-- close to @maxBound@, so check for overflow.+cubeRootInt' :: Int -> Int+cubeRootInt' 0 = 0+cubeRootInt' n+ | n < c || c < 0 = r-1+ | 0 < d && d < n = r+1+ | otherwise = r+ where+ x = fromIntegral n :: Double+ r = truncate (x ** (1/3))+ c = r*r*r+ d = c+3*r*(r+1)++cubeRootWord :: Word -> Word+cubeRootWord 0 = 0+cubeRootWord w+#if WORD_SIZE_IN_BITS == 64+ | r > 2642245 = 2642245+#else+ | r > 1625 = 1625+#endif+ | w < c = r-1+ | c < w && e < w && c < e = r+1+ | otherwise = r+ where+ r = truncate ((fromIntegral w) ** (1/3) :: Double)+ c = r*r*r+ d = 3*r*(r+1)+ e = c+d++{-# SPECIALISE newton3 :: Int -> Int -> Int #-}+{-# SPECIALISE newton3 :: Integer -> Integer -> Integer #-}+newton3 :: Integral a => a -> a -> a+newton3 n a = go (step a)+ where+ step k = (2*k + n `quot` (k*k)) `quot` 3+ go k+ | m < k = go m+ | otherwise = k+ where+ m = step k++{-# SPECIALISE approxCuRt :: Integer -> Integer #-}+approxCuRt :: Integral a => a -> a+approxCuRt 0 = 0+approxCuRt n = fromInteger $ appCuRt (fromIntegral n)++-- threshold for shifting vs. direct fromInteger+-- we shift when we expect more than 256 bits+#if WORD_SIZE_IN_BITS == 64+#define THRESH 5+#else+#define THRESH 9+#endif++-- | 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@(J# s# _)+ | 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+ m -> case floor (fromInteger m ** (1.0/3.0) :: Double) of+ r -> shiftLInteger r h#++-- not very discriminating, but cheap, so it's an overall gain+cr512 :: UArray Int Bool+cr512 = runSTUArray $ do+ ar <- newArray (0,511) True+ let note s i+ | i < 512 = unsafeWrite ar i False >> note s (i+s)+ | otherwise = return ()+ note 4 2+ note 8 4+ note 32 16+ note 64 32+ note 256 128+ unsafeWrite ar 256 False+ return ar++-- Remainders modulo @3^3 * 31@+cubeRes837 :: UArray Int Bool+cubeRes837 = runSTUArray $ do+ ar <- newArray (0,836) False+ let note 837 = return ar+ note k = unsafeWrite ar ((k*k*k) `rem` 837) True >> note (k+1)+ note 0++-- Remainders modulo @7^2 * 13@+cubeRes637 :: UArray Int Bool+cubeRes637 = runSTUArray $ do+ ar <- newArray (0,636) False+ let note 637 = return ar+ note k = unsafeWrite ar ((k*k*k) `rem` 637) True >> note (k+1)+ note 0++-- Remainders modulo @19 * 37@+cubeRes703 :: UArray Int Bool+cubeRes703 = runSTUArray $ do+ ar <- newArray (0,702) False+ let note 703 = return ar+ note k = unsafeWrite ar ((k*k*k) `rem` 703) True >> note (k+1)+ note 0
+ Math/NumberTheory/Powers/Fourth.hs view
@@ -0,0 +1,203 @@+-- |+-- Module: Math.NumberTheory.Powers.Squares+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Functions dealing with fourth powers. Efficient calculation of integer fourth+-- roots and efficient testing for being a square's square.+{-# LANGUAGE MagicHash, BangPatterns, CPP #-}+module Math.NumberTheory.Powers.Fourth+ ( integerFourthRoot+ , integerFourthRoot'+ , exactFourthRoot+ , isFourthPower+ , isFourthPower'+ , isPossibleFourthPower+ ) where++#include "MachDeps.h"++import GHC.Base+import GHC.Integer+import GHC.Integer.GMP.Internals++import Data.Array.Unboxed+import Data.Array.ST+import Data.Array.Base (unsafeAt, unsafeWrite)++import Data.Bits+import Data.Word++import Math.NumberTheory.Logarithms.Internal (integerLog2#)++-- | Calculate the integer fourth root of a nonnegative number,+-- that is, the largest integer @r@ with @r^4 <= n@.+-- Throws an error on negaitve input.+integerFourthRoot :: Integral a => a -> a+integerFourthRoot n+ | n < 0 = error "integerFourthRoot: negative argument"+ | otherwise = integerFourthRoot' n++-- | Calculate the integer fourth root of a nonnegative number,+-- that is, the largest integer @r@ with @r^4 <= n@.+-- The condition is /not/ checked.+{-# RULES+"integerFourthRoot'/Int" integerFourthRoot' = biSqrtInt+"integerFourthRoot'/Word" integerFourthRoot' = biSqrtWord+ #-}+{-# SPECIALISE integerFourthRoot' :: Integer -> Integer #-}+integerFourthRoot' :: Integral a => a -> a+integerFourthRoot' 0 = 0+integerFourthRoot' n = newton4 n (approxBiSqrt n)++-- | Returns @Nothing@ if @n@ is not a fourth power,+-- @Just r@ if @n == r^4@ and @r >= 0@.+{-# SPECIALISE exactFourthRoot :: Int -> Maybe Int,+ Integer -> Maybe Integer,+ Word -> Maybe Word+ #-}+exactFourthRoot :: Integral a => a -> Maybe a+exactFourthRoot 0 = Just 0+exactFourthRoot n+ | n < 0 = Nothing+ | isPossibleFourthPower n && r2*r2 == n = Just r+ | otherwise = Nothing+ where+ r = integerFourthRoot' n+ r2 = r*r++-- | Test whether an integer is a fourth power.+-- First nonnegativity is checked, then the unchecked+-- test is called.+{-# SPECIALISE isFourthPower :: Int -> Bool,+ Integer -> Bool,+ Word -> Bool+ #-}+isFourthPower :: Integral a => a -> Bool+isFourthPower 0 = True+isFourthPower n = n > 0 && isFourthPower' n++-- | Test whether a nonnegative number is a fourth power.+-- The condition is /not/ checked. If a number passes the+-- 'isPossibleFourthPower' test, its integer fourth root+-- is calculated.+{-# SPECIALISE isFourthPower' :: Int -> Bool,+ Integer -> Bool,+ Word -> Bool+ #-}+isFourthPower' :: Integral a => a -> Bool+isFourthPower' n = isPossibleFourthPower n && r2*r2 == n+ where+ r = integerFourthRoot' n+ r2 = r*r++-- | Test whether a nonnegative number is a possible fourth power.+-- The condition is /not/ checked.+-- This eliminates about 99.958% of numbers.+{-# SPECIALISE isPossibleFourthPower :: Int -> Bool,+ Integer -> Bool,+ Word -> Bool+ #-}+isPossibleFourthPower :: Integral a => a -> Bool+isPossibleFourthPower n =+ biSqRes256 `unsafeAt` (fromIntegral n .&. 255)+ && biSqRes425 `unsafeAt` (fromIntegral (n `rem` 425))+ && biSqRes377 `unsafeAt` (fromIntegral (n `rem` 377))++{-# SPECIALISE newton4 :: Integer -> Integer -> Integer #-}+newton4 :: Integral a => a -> a -> a+newton4 n a = go (step a)+ where+ step k = (3*k + n `quot` (k*k*k)) `quot` 4+ go k+ | m < k = go m+ | otherwise = k+ where+ m = step k++{-# SPECIALISE approxBiSqrt :: Integer -> Integer #-}+approxBiSqrt :: Integral a => a -> a+approxBiSqrt = fromInteger . appBiSqrt . fromIntegral++-- threshold for shifting vs. direct fromInteger+-- we shift when we expect more than 384 bits+#if WORD_SIZE_IN_BITS == 64+#define THRESH 7+#else+#define THRESH 13+#endif++-- 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@(J# s# _)+ | 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+ m -> case floor (sqrt $ sqrt $ fromInteger m :: Double) of+ r -> shiftLInteger r h#+++biSqRes256 :: UArray Int Bool+biSqRes256 = runSTUArray $ do+ ar <- newArray (0,255) False+ let note 257 = return ar+ note i = unsafeWrite ar i True >> note (i+16)+ unsafeWrite ar 0 True+ unsafeWrite ar 16 True+ note 1++biSqRes425 :: UArray Int Bool+biSqRes425 = runSTUArray $ do+ ar <- newArray (0,424) False+ let note 154 = return ar+ note i = unsafeWrite ar ((i*i*i*i) `rem` 425) True >> note (i+1)+ note 0++biSqRes377 :: UArray Int Bool+biSqRes377 = runSTUArray $ do+ ar <- newArray (0,376) False+ let note 144 = return ar+ note i = unsafeWrite ar ((i*i*i*i) `rem` 377) True >> note (i+1)+ note 0++biSqrtInt :: Int -> Int+biSqrtInt 0 = 0+biSqrtInt n+#if WORD_SIZE_IN_BITS == 64+ | r > 55108 = 55108+#else+ | r > 215 = 215+#endif+ | n < r4 = r-1+ | otherwise = r+ where+ x :: Double+ x = fromIntegral n+ -- timed faster than x**0.25, never too small+ r = truncate (sqrt (sqrt x))+ r2 = r*r+ r4 = r2*r2++biSqrtWord :: Word -> Word+biSqrtWord 0 = 0+biSqrtWord n+#if WORD_SIZE_IN_BITS == 64+ | r > 65535 = 65535+#else+ | r > 255 = 255+#endif+ | n < r4 = r-1+ | otherwise = r+ where+ x :: Double+ x = fromIntegral n+ r = truncate (sqrt (sqrt x))+ r2 = r*r+ r4 = r2*r2+
+ Math/NumberTheory/Powers/General.hs view
@@ -0,0 +1,338 @@+-- |+-- Module: Math.NumberTheory.Powers.General+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Calculating integer roots and determining perfect powers.+-- The algorithms are moderately efficient.+--+{-# LANGUAGE MagicHash, BangPatterns, CPP #-}+{-# OPTIONS_GHC -O2 -fspec-constr-count=8 #-}+module Math.NumberTheory.Powers.General+ ( integerRoot+ , exactRoot+ , isKthPower+ , isPerfectPower+ , highestPower+-- , largePFPower+ ) where++#include "MachDeps.h"++import GHC.Base+import GHC.Integer+import GHC.Integer.GMP.Internals++import Data.Bits+import Data.Word+import Data.List (foldl')+import qualified Data.Set as Set++-- import Math.NumberTheory.Logarithms+import Math.NumberTheory.Logarithms.Internal (integerLog2#)+import Math.NumberTheory.Utils (shiftToOddCount, splitOff)+import qualified Math.NumberTheory.Powers.Squares as P2+import qualified Math.NumberTheory.Powers.Cubes as P3+import qualified Math.NumberTheory.Powers.Fourth as P4++-- | Calculate an integer root, @'integerRoot' k n@ computes the (floor of) the @k@-th+-- root of @n@, where @k@ must be positive.+-- @r = 'integerRoot' k n@ means @r^k <= n < (r+1)^k@ if that is possible at all.+-- It is impossible if @k@ is even and @n \< 0@, since then @r^k >= 0@ for all @r@,+-- then, and if @k <= 0@, @'integerRoot'@ raises an error. For @k < 5@, a specialised+-- version is called which should be more efficient than the general algorithm.+-- However, it is not guaranteed that the rewrite rules for those fire, so if @k@ is+-- known in advance, it is safer to directly call the specialised versions.+{-# SPECIALISE integerRoot :: Int -> Int -> Int,+ Int -> Word -> Word,+ Int -> Integer -> Integer,+ Word -> Int -> Int,+ Word -> Word -> Word,+ Word -> Integer -> Integer,+ Integer -> Integer -> Integer+ #-}+integerRoot :: (Integral a, Integral b) => b -> a -> a+integerRoot 1 n = n+integerRoot 2 n = P2.integerSquareRoot n+integerRoot 3 n = P3.integerCubeRoot n+integerRoot 4 n = P4.integerFourthRoot n+integerRoot k n+ | k < 1 = error "integerRoot: negative exponent or exponent 0"+ | n < 0 && even k = error "integerRoot: negative radicand for even exponent"+ | n < 0 =+ let r = negate . fromInteger . integerRoot k . negate $ fromIntegral n+ in if r^k == n then r else (r-1)+ | n == 0 = 0+ | n < 31 = 1+ | kTooLarge = 1+ | otherwise = newtonK k' n a+ where+ k' = fromIntegral k+ a = approxKthRoot (fromIntegral 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+ || (I# (integerLog2# (toInteger n)) < fromIntegral k) -- n < 2^k++-- | @'exactRoot' k n@ returns @'Nothing'@ if @n@ is not a @k@-th power,+-- @'Just' r@ if @n == r^k@. If @k@ is divisible by @4, 3@ or @2@, a+-- residue test is performed to avoid the expensive calculation if it+-- can thus be determined that @n@ is not a @k@-th power.+exactRoot :: (Integral a, Integral b) => b -> a -> Maybe a+exactRoot 1 n = Just n+exactRoot 2 n = P2.exactSquareRoot n+exactRoot 3 n = P3.exactCubeRoot n+exactRoot 4 n = P4.exactFourthRoot n+exactRoot k n+ | n == 1 = Just 1+ | k < 1 = Nothing+ | n < 0 && even k = Nothing+ | n < 0 = fmap negate (exactRoot k (-n))+ | n < 2 = Just n+ | n < 31 = Nothing+ | kTooLarge = Nothing+ | otherwise = case k `rem` 12 of+ 0 | c4 && c3 && ok -> Just r+ | otherwise -> Nothing+ 2 | c2 && ok -> Just r+ | otherwise -> Nothing+ 3 | c3 && ok -> Just r+ | otherwise -> Nothing+ 4 | c4 && ok -> Just r+ | otherwise -> Nothing+ 6 | c3 && c2 && ok -> Just r+ | otherwise -> Nothing+ 8 | c4 && ok -> Just r+ | otherwise -> Nothing+ 9 | c3 && ok -> Just r+ | otherwise -> Nothing+ 10 | c2 && ok -> Just r+ | otherwise -> Nothing+ _ | ok -> Just r+ | otherwise -> Nothing++ where+ k' :: Int+ k' = fromIntegral k+ r = integerRoot k' n+ c2 = P2.isPossibleSquare n+ c3 = P3.isPossibleCube n+ c4 = P4.isPossibleFourthPower n+ 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+ || (I# (integerLog2# (toInteger n)) < fromIntegral k) -- n < 2^k++-- | @'isKthPower' k n@ checks whether @n@ is a @k@-th power.+isKthPower :: (Integral a, Integral b) => b -> a -> Bool+isKthPower k n = case exactRoot k n of+ Just _ -> True+ Nothing -> False++-- | @'isPerfectPower' n@ checks whether @n == r^k@ for some @k > 1@.+isPerfectPower :: Integral a => a -> Bool+isPerfectPower n+ | n == 0 || n == 1 = True+ | otherwise = k > 1+ where+ (_,k) = highestPower n++-- | @'highestPower' n@ produces the pair @(b,k)@ with the largest+-- exponent @k@ such that @n == b^k@, except for @'abs' n <= 1@,+-- in which case arbitrarily large exponents exist, and by an+-- arbitrary decision @(n,3)@ is returned.+--+-- First, by trial division with small primes, the range of possible+-- exponents is reduced (if @p^e@ exactly divides @n@, then @k@ must+-- be a divisor of @e@, if several small primes divide @n@, @k@ must+-- divide the greatest common divisor of their exponents, which mostly+-- will be @1@, generally small; if none of the small primes divides+-- @n@, the range of possible exponents is reduced since the base is+-- necessarily large), if that has not yet determined the result, the+-- remaining factor is examined by trying the divisors of the @gcd@+-- of the prime exponents if some have been found, otherwise by trying+-- prime exponents recursively.+highestPower :: Integral a => a -> (a, Int)+highestPower n+ | abs n <= 1 = (n,3)+ | n < 0 = case integerHighPower (toInteger $ negate n) of+ (r,e) -> case shiftToOddCount e of+ (k, o) -> (negate $ fromInteger (sqr k r), o)+ | otherwise = case integerHighPower (toInteger n) of+ (r,e) -> (fromInteger r, e)+ where+ sqr :: Int -> Integer -> Integer+ 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++------------------------------------------------------------------------------------------+-- Auxiliary functions --+------------------------------------------------------------------------------------------++{-# SPECIALISE newtonK :: Int -> Int -> Int -> Int,+ Integer -> Integer -> Integer -> Integer,+ Word -> Word -> Word -> Word+ #-}+newtonK :: Integral a => a -> a -> a -> a+newtonK k n a = go (step a)+ where+ step m = ((k-1)*m + n `quot` (m^(k-1))) `quot` k+ go m+ | l < m = go l+ | otherwise = m+ where+ l = step m++{-# SPECIALISE approxKthRoot :: Int -> Integer -> Integer,+ Int -> Int -> Int,+ Int -> Word -> Word+ #-}+approxKthRoot :: Integral a => Int -> a -> a+approxKthRoot k = fromInteger . appKthRoot k . fromIntegral++-- 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 k@(I# k#) n@(J# _ _) =+ case integerLog2# n of+ l# -> case l# `quotInt#` k# of+ 0# -> 1+ 1# -> 3+ 2# -> 5+ 3# -> 11+ h# | h# <# 500# ->+ floor (scaleFloat (I# (h# -# 1#))+ (fromInteger (n `shiftRInteger` (h# *# k# -# k#)) ** (1/fromIntegral k) :: Double))+ | otherwise ->+ floor (scaleFloat 400 (fromInteger (n `shiftRInteger` (h# *# k# -# k#)) ** (1/fromIntegral k) :: Double))+ `shiftLInteger` (h# -# 401#)++-- assumption: argument is > 1+integerHighPower :: Integer -> (Integer, Int)+integerHighPower n+ | n < 4 = (n,1)+ | otherwise = case shiftToOddCount n of+ (e2,m) | m == 1 -> (2,e2)+ | otherwise -> findHighPower e2 (if e2 == 0 then [] else [(2,e2)]) m r smallOddPrimes+ where+ r = P2.integerSquareRoot m++findHighPower :: Int -> [(Integer,Int)] -> Integer -> Integer -> [Integer] -> (Integer, Int)+findHighPower 1 pws m _ _ = (foldl' (*) m [p^e | (p,e) <- pws], 1)+findHighPower e pws 1 _ _ = (foldl' (*) 1 [p^(ex `quot` e) | (p,ex) <- pws], e)+findHighPower e pws m s (p:ps)+ | s < p = findHighPower 1 pws m s []+ | otherwise =+ case splitOff p m of+ (0,_) -> findHighPower e pws m s ps+ (k,r) -> findHighPower (gcd k e) ((p,k):pws) r (P2.integerSquareRoot r) ps+findHighPower e pws m _ [] = finishPower e pws m++spBEx :: Int+spBEx = 14++spBound :: Integer+spBound = 2^spBEx++smallOddPrimes :: [Integer]+smallOddPrimes = 3:5:primes'+ where+ primes' = 7:11:13:17:19:23:29:filter isPrime (takeWhile (< spBound) $ scanl (+) 31 (cycle [6,4,2,4,2,4,6,2]))+ isPrime n = go primes'+ where+ go (p:ps) = (p*p > n) || (n `rem` p /= 0 && go ps)+ go [] = True++-- n large, has no prime divisors < spBound+finishPower :: Int -> [(Integer, Int)] -> Integer -> (Integer, Int)+finishPower e pws n+ | n < (1 `shiftL` (2*spBEx)) = (foldl' (*) n [p^ex | (p,ex) <- pws], 1) -- n is prime+ | e == 0 = rawPower maxExp n+ | otherwise = go divs+ where+ maxExp = (I# (integerLog2# n)) `quot` spBEx+ divs = divisorsTo maxExp e+ go [] = (foldl' (*) n [p^ex | (p,ex) <- pws], 1)+ go (d:ds) = case exactRoot d n of+ Just r -> (foldl' (*) r [p^(ex `quot` d) | (p,ex) <- pws], d)+ Nothing -> go ds++rawPower :: Int -> Integer -> (Integer, Int)+rawPower mx n+ | mx < 2 = (n,1)+rawPower mx n = case P4.exactFourthRoot n of+ Just r -> case rawPower (mx `quot` 4) r of+ (m,e) -> (m, 4*e)+ Nothing -> case P2.exactSquareRoot n of+ Just r -> case rawOddPower (mx `quot` 2) r of+ (m,e) -> (m, 2*e)+ Nothing -> rawOddPower mx n++rawOddPower :: Int -> Integer -> (Integer, Int)+rawOddPower mx n+ | mx < 3 = (n,1)+rawOddPower mx n = case P3.exactCubeRoot n of+ Just r -> case rawOddPower (mx `quot` 3) r of+ (m,e) -> (m, 3*e)+ Nothing -> badPower mx n++badPower :: Int -> Integer -> (Integer, Int)+badPower mx n+ | mx < 5 = (n,1)+ | otherwise = go 1 mx n (takeWhile (<= mx) $ scanl (+) 5 $ cycle [2,4])+ where+ go !e b m (k:ks)+ | b < k = (m,e)+ | otherwise = case exactRoot k m of+ Just r -> go (e*k) (b `quot` k) r (k:ks)+ Nothing -> go e b m ks+ go e _ m [] = (m,e)++divisorsTo :: Int -> Int -> [Int]+divisorsTo mx n = case shiftToOddCount n of+ (k,o) | k == 0 -> go (Set.singleton 1) n iops+ | otherwise -> go (Set.fromDistinctAscList $ takeWhile (<= mx) $ take (k+1) (iterate (*2) 1)) o iops+ where+ mset k st = fst (Set.split (mx+1) (Set.mapMonotonic (*k) st))+ unP :: Int -> Int -> (Int,Int)+ unP p m = goP 0 m+ where+ goP :: Int -> Int -> (Int,Int)+ goP !i j = case m `quotRem` p of+ (q,r) | r == 0 -> goP (i+1) q+ | otherwise -> (0,j)+ iops :: [Int]+ iops = 3:5:prs+ prs :: [Int]+ prs = 7:filter prm (scanl (+) 11 $ cycle [2,4,2,4,6,2,6,4])+ prm :: Int -> Bool+ prm k = td prs+ where+ td (p:ps) = (p*p > k) || (k `rem` p /= 0 && td ps)+ td [] = True+ go !st m (p:ps)+ | m == 1 = reverse $ Set.toAscList st+ | m < p*p = reverse . Set.toAscList $ Set.union st (mset m st)+ | otherwise =+ case unP p m of+ (0,_) -> go st m ps+ (k,r) -> go (Set.unions (st:take k (iterate (mset p) st))) r ps+ go st m [] = go st m [m+1]
+ Math/NumberTheory/Powers/Integer.hs view
@@ -0,0 +1,40 @@+-- |+-- Module: Math.NumberTheory.Powers.Integer+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Slightly faster power function for Integer base and Int exponent.+--+{-# LANGUAGE MagicHash, BangPatterns #-}+{-# OPTIONS_HADDOCK hide #-}+module Math.NumberTheory.Powers.Integer+ ( integerPower+ ) where++import GHC.Base++import Math.NumberTheory.Logarithms.Internal ( 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.+integerPower :: Integer -> Int -> Integer+integerPower b (I# e#)+ | e# ==# 0# = 1+ | 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 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##
+ Math/NumberTheory/Powers/Squares.hs view
@@ -0,0 +1,272 @@+-- |+-- Module: Math.NumberTheory.Powers.Squares+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Functions dealing with squares. Efficient calculation of integer square roots+-- and efficient testing for squareness.+{-# LANGUAGE MagicHash, BangPatterns, CPP #-}+module Math.NumberTheory.Powers.Squares+ ( -- * Square root calculation+ integerSquareRoot+ , integerSquareRoot'+ , exactSquareRoot+ -- * Tests for squares+ , isSquare+ , isSquare'+ , isPossibleSquare+ , isPossibleSquare2+ ) where++#include "MachDeps.h"++import GHC.Base+import GHC.Integer+import GHC.Integer.GMP.Internals++import Data.Array.Unboxed+import Data.Array.ST+import Data.Array.Base (unsafeAt, unsafeWrite)++import Data.Bits+import Data.Word++import Math.NumberTheory.Logarithms.Internal (integerLog2#)+++-- | Calculate the integer square root of a nonnegative number @n@,+-- that is, the largest integer @r@ with @r*r <= n@.+-- Throws an error on negative input.+{-# SPECIALISE integerSquareRoot :: Int -> Int,+ Word -> Word,+ Integer -> Integer+ #-}+integerSquareRoot :: Integral a => a -> a+integerSquareRoot n+ | n < 0 = error "integerSquareRoot: negative argument"+ | otherwise = integerSquareRoot' n++-- | Calculate the integer square root of a nonnegative number @n@,+-- that is, the largest integer @r@ with @r*r <= n@.+-- The precondition @n >= 0@ is not checked.+{-# RULES+"integerSquareRoot'/Int" integerSquareRoot' = isqrtInt'+"integerSquareRoot'/Word" integerSquareRoot' = isqrtWord+ #-}+{-# SPECIALISE integerSquareRoot' :: Integer -> Integer #-}+integerSquareRoot' :: Integral a => a -> a+integerSquareRoot' = isqrtA++-- | Returns 'Nothing' if the argument is not a square,+-- @'Just' r@ if @r*r == n@ and @r >= 0@. Avoids the expensive calculation+-- of the square root if @n@ is recognized as a non-square+-- before, prevents repeated calculation of the square root+-- if only the roots of perfect squares are needed.+-- Checks for negativity and 'isPossibleSquare'.+{-# SPECIALISE exactSquareRoot :: Int -> Maybe Int,+ Word -> Maybe Word,+ Integer -> Maybe Integer+ #-}+exactSquareRoot :: Integral a => a -> Maybe a+exactSquareRoot n+ | n < 0 = Nothing+ | isPossibleSquare n && r*r == n = Just r+ | otherwise = Nothing+ where+ r = integerSquareRoot' n++-- | Test whether the argument is a square.+-- After a number is found to be positive, first 'isPossibleSquare'+-- is checked, if it is, the integer square root is calculated.+{-# SPECIALISE isSquare :: Int -> Bool,+ Word -> Bool,+ Integer -> Bool+ #-}+isSquare :: Integral a => a -> Bool+isSquare n = n >= 0 && isSquare' n++-- | Test whether the input (a nonnegative number) @n@ is a square.+-- The same as 'isSquare', but without the negativity test.+-- Faster if many known positive numbers are tested.+--+-- The precondition @n >= 0@ is not tested, passing negative+-- arguments may cause any kind of havoc.+{-# SPECIALISE isSquare' :: Int -> Bool,+ Word -> Bool,+ Integer -> Bool+ #-}+isSquare' :: Integral a => a -> Bool+isSquare' n = isPossibleSquare n && let r = integerSquareRoot' n in r*r == n++-- | Test whether a non-negative number may be a square.+-- Non-negativity is not checked, passing negative arguments may+-- cause any kind of havoc.+--+-- First the remainder modulo 256 is checked (that can be calculated+-- easily without division and eliminates about 82% of all numbers).+-- After that, the remainders modulo 9, 25, 7, 11 and 13 are tested+-- to eliminate altogether about 99.436% of all numbers.+--+-- This is the test used by 'exactSquareRoot'. For large numbers,+-- the slower but more discriminating test 'isPossibleSqure2' is+-- faster.+{-# SPECIALISE isPossibleSquare :: Int -> Bool,+ Integer -> Bool,+ Word -> Bool+ #-}+isPossibleSquare :: Integral a => a -> Bool+isPossibleSquare n =+ unsafeAt sr256 ((fromIntegral n) .&. 255)+ && unsafeAt sr693 (fromIntegral (n `rem` 693))+ && unsafeAt sr325 (fromIntegral (n `rem` 325))++-- | Test whether a non-negative number may be a square.+-- Non-negativity is not checked, passing negative arguments may+-- cause any kind of havoc.+--+-- First the remainder modulo 256 is checked (that can be calculated+-- easily without division and eliminates about 82% of all numbers).+-- After that, the remainders modulo several small primes are tested+-- to eliminate altogether about 99.98954% of all numbers.+--+-- For smallish to medium sized numbers, this hardly performs better+-- than 'isPossibleSquare', which uses smaller arrays, but for large+-- numbers, where calculating the square root becomes more expensive,+-- it is much faster (if the vast majority of tested numbers aren't squares).+{-# SPECIALISE isPossibleSquare2 :: Int -> Bool,+ Integer -> Bool,+ Word -> Bool+ #-}+isPossibleSquare2 :: Integral a => a -> Bool+isPossibleSquare2 n =+ unsafeAt sr256 ((fromIntegral n) .&. 255)+ && unsafeAt sr819 (fromIntegral (n `rem` 819))+ && unsafeAt sr1025 (fromIntegral (n `rem` 1025))+ && unsafeAt sr2047 (fromIntegral (n `rem` 2047))+ && unsafeAt sr4097 (fromIntegral (n `rem` 4097))+ && unsafeAt sr341 (fromIntegral (n `rem` 341))++-----------------------------------------------------------------------------+-- Auxiliary Stuff++-- Find approximation to square root in 'Integer', then+-- find the integer square root by the integer variant+-- of Heron's method. Takes only a handful of steps+-- unless the input is really large.+{-# SPECIALISE isqrtA :: Integer -> Integer #-}+isqrtA :: Integral a => a -> a+isqrtA 0 = 0+isqrtA n = heron n (fromInteger . appSqrt . fromIntegral $ n)++-- Heron's method for integers. First make one step to ensure+-- the value we're working on is @>= r@, then we have+-- @k == r@ iff @k <= step k@.+{-# SPECIALISE heron :: Integer -> Integer -> Integer #-}+heron :: Integral a => a -> a -> a+heron n a = go (step a)+ where+ step k = (k + n `quot` k) `quot` 2+ go k+ | m < k = go m+ | otherwise = k+ where+ m = step k++-- threshold for shifting vs. direct fromInteger+-- we shift when we expect more than 256 bits+#if WORD_SIZE_IN_BITS == 64+#define THRESH 5+#else+#define THRESH 9+#endif++-- Find a fairly good approximation to the square root.+-- 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@(J# s# _)+ | 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+ m -> case floor (sqrt $ fromInteger m :: Double) of+ r -> shiftLInteger r h#++-- Auxiliaries++-- Make an array indicating whether a remainder is a square remainder.+sqRemArray :: Int -> UArray Int Bool+sqRemArray md = runSTUArray $ do+ arr <- newArray (0,md-1) False+ let !stop = (md `quot` 2) + 1+ fill k+ | k < stop = unsafeWrite arr ((k*k) `rem` md) True >> fill (k+1)+ | otherwise = return arr+ unsafeWrite arr 0 True+ unsafeWrite arr 1 True+ fill 2++sr256 :: UArray Int Bool+sr256 = sqRemArray 256++sr819 :: UArray Int Bool+sr819 = sqRemArray 819++sr4097 :: UArray Int Bool+sr4097 = sqRemArray 4097++sr341 :: UArray Int Bool+sr341 = sqRemArray 341++sr1025 :: UArray Int Bool+sr1025 = sqRemArray 1025++sr2047 :: UArray Int Bool+sr2047 = sqRemArray 2047++sr693 :: UArray Int Bool+sr693 = sqRemArray 693++sr325 :: UArray Int Bool+sr325 = sqRemArray 325++-- Specialisations for Int and Word++-- For @n <= 2^64@, the result of+--+-- > truncate (sqrt $ fromIntegral n)+--+-- is never too small and never more than one too large.+-- The multiplication doesn't overflow for 32 or 64 bit Ints.+isqrtInt' :: Int -> Int+isqrtInt' n+ | n < r*r = r-1+ | otherwise = r+ where+ !r = (truncate :: Double -> Int) . sqrt $ fromIntegral n+-- With -O2, that should be translated to the below+{-+isqrtInt' n@(I# i#)+ | r# *# r# ># i# = I# (r# -# 1#)+ | otherwise = I# r#+ where+ !r# = double2Int# (sqrtDouble# (int2Double# i#))+-}++-- Same for Word.+isqrtWord :: Word -> Word+isqrtWord n+ | n < (r*r)+#if WORD_SIZE_IN_BITS == 64+ || r == 4294967296+-- Double interprets values near maxBound as 2^64, we don't have that problem for 32 bits+#endif+ = r-1+ | otherwise = r+ where+ !r = (fromIntegral :: Int -> Word) . (truncate :: Double -> Int) . sqrt $ fromIntegral n+
+ Math/NumberTheory/Primes.hs view
@@ -0,0 +1,19 @@+-- |+-- Module: Math.NumberTheory.Primes+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+module Math.NumberTheory.Primes+ ( module Math.NumberTheory.Primes.Sieve+ , module Math.NumberTheory.Primes.Counting+ , module Math.NumberTheory.Primes.Testing+ , module Math.NumberTheory.Primes.Factorisation+ ) where++import Math.NumberTheory.Primes.Sieve+import Math.NumberTheory.Primes.Counting+import Math.NumberTheory.Primes.Testing hiding (FactorSieve)+import Math.NumberTheory.Primes.Factorisation
+ Math/NumberTheory/Primes/Counting.hs view
@@ -0,0 +1,21 @@+-- |+-- Module: Math.NumberTheory.Primes.Counting+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: non-portable+--+-- Number of primes not exceeding @n@, @π(n)@, and @n@-th prime; also fast, but+-- reasonably accurate approximations to these.+module Math.NumberTheory.Primes.Counting+ ( -- * Exact functions+ primeCount+ , nthPrime+ -- * Approximations+ , approxPrimeCount+ , nthPrimeApprox+ ) where++import Math.NumberTheory.Primes.Counting.Impl+import Math.NumberTheory.Primes.Counting.Approximate
+ Math/NumberTheory/Primes/Counting/Approximate.hs view
@@ -0,0 +1,62 @@+-- |+-- Module: Math.NumberTheory.Primes.Counting.Approximate+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: portable+--+-- Approximations to the number of primes below a limit and the+-- n-th prime.+--+{-# OPTIONS_HADDOCK hide #-}+module Math.NumberTheory.Primes.Counting.Approximate+ ( approxPrimeCount+ , nthPrimeApprox+ ) where++-- | @'approxPrimeCount' n@ gives (for @n > 0@) an+-- approximation of the number of primes not exceeding+-- @n@. The approximation is fairly good for @n@ large enough.+-- The number of primes should be slightly overestimated+-- (so it is suitable for allocation of storage) and is+-- never underestimated for @n <= 10^12@.+approxPrimeCount :: Integral a => a -> a+approxPrimeCount = truncate . appi . fromIntegral++-- | @'nthPrimeApprox' n@ gives (for @n > 0@) an+-- approximation to the n-th prime. The approximation+-- is fairly good for @n@ large enough. Dual to+-- @'approxPrimeCount'@, this estimate should err+-- on the low side (and does for @n < 10^12@).+nthPrimeApprox :: Integral a => a -> a+nthPrimeApprox = max 1 . truncate . nthApp . fromIntegral . max 3++-- Basically the approximation of the prime count by Li(x),+-- adjusted to give close but slightly too high estimates+-- in the interesting range. The constants are empirically+-- determined.+appi :: Double -> Double+appi x = y - y/300000 + 7*ll+ where+ y = x*l*(1+l*(1+l*h))+ w = log x+ l = 1/w+ ll = log w+ h | x < 10000000 = 2.5625+ | x < 50000000 = 2.5+ | x < 120000000 = 617/256+ | otherwise = 2.0625 + l*(3+ll*l*(13.25+ll*l*57.75))++-- Basically an approximation to the inverse of Li(x), with+-- empirically determined constants to get close results+-- in the interesting range.+nthApp :: Double -> Double+nthApp x = a+ where+ l = log x+ ll = log l+ li = 1/l+ l2 = ll*ll+ a = x*(l+ll-1+li*(ll-2-li*(ll*(0.3+li*(1+0.02970812*l2*l2*l2*li))+ + 8.725*(ll-2.749)*(ll-3.892)*li))) + l*ll + 35
+ Math/NumberTheory/Primes/Counting/Impl.hs view
@@ -0,0 +1,383 @@+-- |+-- Module: Math.NumberTheory.Primes.Counting.Impl+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: non-portable+--+-- Number of primes not exceeding @n@, @π(n)@, and @n@-th prime.+--+{-# LANGUAGE CPP, BangPatterns #-}+#if __GLASGOW_HASKELL__ >= 700+{-# OPTIONS_GHC -fspec-constr-count=24 #-}+#endif+{-# OPTIONS_HADDOCK hide #-}+module Math.NumberTheory.Primes.Counting.Impl (primeCount, nthPrime) where++#include "MachDeps.h"++import Math.NumberTheory.Primes.Sieve.Eratosthenes+import Math.NumberTheory.Primes.Sieve.Indexing+import Math.NumberTheory.Primes.Counting.Approximate+import Math.NumberTheory.Powers.Squares+import Math.NumberTheory.Powers.Cubes+import Math.NumberTheory.Logarithms++import Data.Array.Base+import Data.Array.ST+import Control.Monad.ST+import Data.Bits+import Data.Int++#if SIZEOF_HSWORD < 8+#define COUNT_T Int64+#else+#define COUNT_T Int+#endif++-- | @'primeCount' n == π(n)@ is the number of (positive) primes not exceeding @n@.+--+-- For efficiency, the calculations are done on 64-bit signed integers, therefore @n@ must+-- not exceed @8 * 10^18@.+--+-- Requires @/O/(n^0.5)@ space, the time complexity is roughly @/O/(n^0.7)@.+-- For small bounds, @'primeCount' n@ simply counts the primes not exceeding @n@,+-- for bounds from @30000@ on, Meissel's algorithm is used in the improved form due to+-- D.H. Lehmer, cf.+-- <http://en.wikipedia.org/wiki/Prime_counting_function#Algorithms_for_evaluating_.CF.80.28x.29>.+primeCount :: Integer -> Integer+primeCount n+ | n > 8000000000000000000 = error $ "primeCount: can't handle bound " ++ show n+ | n < 2 = 0+ | n < 1000 = fromIntegral . length . takeWhile (<= n) . primeList . primeSieve $ max 242 n+ | n < 30000 = runST $ do+ ba <- sieveTo n+ (s,e) <- getBounds ba+ ct <- countFromTo s e ba+ return (fromIntegral $ ct+3)+ | otherwise =+ let !ub = cop $ fromInteger n+ !sr = integerSquareRoot' ub+ !cr = nxtEnd $ integerCubeRoot' ub + 15+ nxtEnd k = k - (k `rem` 30) + 31+ !phn1 = calc ub cr+ !cs = cr+6+ !pdf = sieveCount ub cs sr+ in phn1 - pdf++-- | @'nthPrime' n@ calculates the @n@-th prime. Numbering of primes is+-- @1@-based, so @'nthPrime' 1 == 2@.+--+-- Requires @/O/((n*log n)^0.5)@ space, the time complexity is roughly @/O/((n*log n)^0.7@.+-- The argument must be strictly positive, and must not exceed @1.5 * 10^17@.+nthPrime :: Integer -> Integer+nthPrime n+ | n < 1 = error "Prime indexing starts at 1"+ | n > 150000000000000000 = error $ "nthPrime: can't handle index " ++ show n+ | n < 4000 = nthPrimeCt n+ | ct0 < n = tooLow n p0 (n-ct0) approxGap+ | otherwise = tooHigh n p0 (ct0-n) approxGap+ where+ p0 = nthPrimeApprox n+ approxGap = (7 * fromIntegral (integerLog2' p0)) `quot` 10+ ct0 = primeCount p0++--------------------------------------------------------------------------------+-- The Works --+--------------------------------------------------------------------------------++-- TODO: do something better in case we guess too high.+-- Not too pressing, since I think a) nthPrimeApprox always underestimates+-- in the range we can handle, and b) it's always "goodEnough"++tooLow :: Integer -> Integer -> Integer -> Integer -> Integer+tooLow n a miss gap+ | goodEnough = lowSieve a miss+ | c1 < n = lowSieve p1 (n-c1)+ | otherwise = lowSieve a miss -- a third count wouldn't make it faster, I think+ where+ est = miss*gap+ p1 = a + (est * 19) `quot` 20+ goodEnough = 3*est*est*est < 2*p1*p1 -- a second counting would be more work than sieving+ c1 = primeCount p1++tooHigh :: Integer -> Integer -> Integer -> Integer -> Integer+tooHigh n a surp gap+ | c < n = lowSieve b (n-c)+ | otherwise = tooHigh n b (c-n) gap+ where+ b = a - (surp * gap * 11) `quot` 10+ c = primeCount b++lowSieve :: Integer -> Integer -> Integer+lowSieve a miss = countToNth (miss+rep) psieves+ where+ strt = if (fromInteger a .&. (1 :: Int)) == 1+ then a+2+ else a+1+ psieves@(PS vO ba:_) = psieveFrom strt+ rep | o0 < 0 = 0+ | otherwise = sum [1 | i <- [0 .. r2], ba `unsafeAt` i]+ where+ o0 = a - vO - 7+ r0 = fromInteger o0 `rem` 30+ r1 = r0 `quot` 3+ r2 = min 7 (if r1 > 5 then r1-1 else r1)++-- highSieve :: Integer -> Integer -> Integer -> Integer+-- highSieve a surp gap = error "Oh shit"++sieveCount :: COUNT_T -> COUNT_T -> COUNT_T -> Integer+sieveCount ub cr sr = runST $ do+ let psieves = psieveFrom (fromIntegral cr)+ pisr = approxPrimeCount sr+ picr = approxPrimeCount cr+ diff = pisr - picr+ size = fromIntegral (diff + diff `quot` 50) + 30+ store <- unsafeNewArray_ (0,size-1) :: ST s (STUArray s Int COUNT_T)+ let feed voff !wi !ri uar sves+ | ri == sieveBits = case sves of+ (PS vO ba : more) -> feed (fromInteger vO) wi 0 ba more+ _ -> error "prime stream ended prematurely"+ | pval > sr = do+ stu <- unsafeThaw uar+ eat 0 0 voff (wi-1) ri stu sves+ | uar `unsafeAt` ri = do+ unsafeWrite store wi (ub `quot` pval)+ feed voff (wi+1) (ri+1) uar sves+ | otherwise = feed voff wi (ri+1) uar sves+ where+ pval = voff + toPrim ri+ eat !acc !btw voff !wi !si stu sves+ | si == sieveBits =+ case sves of+ [] -> error "Premature end of prime stream"+ (PS vO ba : more) -> do+ nstu <- unsafeThaw ba+ eat acc btw (fromInteger vO) wi 0 nstu more+ | wi < 0 = return acc+ | otherwise = do+ qb <- unsafeRead store wi+ let dist = qb - voff - 7+ if dist < fromIntegral sieveRange+ then do+ let (b,j) = idxPr (dist+7)+ !li = (b `shiftL` 3) .|. j+ new <- if li < si then return 0 else countFromTo si li stu+ let nbtw = btw + fromIntegral new + 1+ eat (acc+nbtw) nbtw voff (wi-1) (li+1) stu sves+ else do+ let (cpl,fds) = dist `quotRem` fromIntegral sieveRange+ (b,j) = idxPr (fds+7)+ !li = (b `shiftL` 3) .|. j+ ctLoop !lac 0 (PS vO ba : more) = do+ nstu <- unsafeThaw ba+ new <- countFromTo 0 li nstu+ let nbtw = btw + lac + 1 + fromIntegral new+ eat (acc+nbtw) nbtw (fromIntegral vO) (wi-1) (li+1) nstu more+ ctLoop lac s (ps : more) = do+ !new <- countAll ps+ ctLoop (lac + fromIntegral new) (s-1) more+ ctLoop _ _ [] = error "Primes ended"+ new <- countFromTo si (sieveBits-1) stu+ ctLoop (fromIntegral new) (cpl-1) sves+ case psieves of+ (PS vO ba : more) -> feed (fromInteger vO) 0 0 ba more+ _ -> error "No primes sieved"++calc :: COUNT_T -> COUNT_T -> Integer+calc lim plim = runST $ do+ !parr <- sieveTo (fromIntegral plim)+ (plo,phi) <- getBounds parr+ !pct <- countFromTo plo phi parr+ !ar1 <- unsafeNewArray_ (0,end-1)+ unsafeWrite ar1 0 lim+ unsafeWrite ar1 1 1+ !ar2 <- unsafeNewArray_ (0,end-1)+ let go cap pix old new+ | pix == 2 = coll cap old+ | otherwise = do+ isp <- unsafeRead parr pix+ if isp+ then do+ let !n = fromInteger (toPrim pix)+ !ncap <- treat cap n old new+ go ncap (pix-1) new old+ else go cap (pix-1) old new+ coll stop ar =+ let cgo !acc i+ | i < stop = do+ !k <- unsafeRead ar i+ !v <- unsafeRead ar (i+1)+ cgo (acc + fromIntegral v*cp6 k) (i+2)+ | otherwise = return (acc+fromIntegral pct+2)+ in cgo 0 0+ go 2 start ar1 ar2+ where+ (bt,ri) = idxPr plim+ !start = 8*bt + ri+ !size = fromIntegral $ (integerSquareRoot lim) `quot` 4+ !end = 2*size++treat :: Int -> COUNT_T -> STUArray s Int COUNT_T -> STUArray s Int COUNT_T -> ST s Int+treat end n old new = do+ qi0 <- locate n 0 (end `quot` 2 - 1) old+ let collect stop !acc ix+ | ix < end = do+ !k <- unsafeRead old ix+ if k < stop+ then do+ v <- unsafeRead old (ix+1)+ collect stop (acc-v) (ix+2)+ else return (acc,ix)+ | otherwise = return (acc,ix)+ goTreat !wi !ci qi+ | qi < end = do+ !key <- unsafeRead old qi+ !val <- unsafeRead old (qi+1)+ let !q0 = key `quot` n+ !r0 = fromIntegral (q0 `rem` 30030)+ !nkey = q0 - fromIntegral (cpDfAr `unsafeAt` r0)+ nk0 = q0 + fromIntegral (cpGpAr `unsafeAt` (r0+1) + 1)+ !nlim = n*nk0+ (wi1,ci1) <- copyTo end nkey old ci new wi+ ckey <- unsafeRead old ci1+ (!acc, !ci2) <- if ckey == nkey+ then do+ !ov <- unsafeRead old (ci1+1)+ return (ov-val,ci1+2)+ else return (-val,ci1)+ (!tot, !nqi) <- collect nlim acc (qi+2)+ unsafeWrite new wi1 nkey+ unsafeWrite new (wi1+1) tot+ goTreat (wi1+2) ci2 nqi+ | otherwise = copyRem end old ci new wi+ goTreat 0 0 qi0++--------------------------------------------------------------------------------+-- Auxiliaries --+--------------------------------------------------------------------------------++locate :: COUNT_T -> Int -> Int -> STUArray s Int COUNT_T -> ST s Int+locate p low high arr = do+ let go lo hi+ | lo < hi = do+ let !md = (lo+hi) `quot` 2+ v <- unsafeRead arr (2*md)+ case compare p v of+ LT -> go lo md+ EQ -> return (2*md)+ GT -> go (md+1) hi+ | otherwise = return (2*lo)+ go low high++{-# INLINE copyTo #-}+copyTo :: Int -> COUNT_T -> STUArray s Int COUNT_T -> Int+ -> STUArray s Int COUNT_T -> Int -> ST s (Int,Int)+copyTo end lim old oi new ni = do+ let go ri wi+ | ri < end = do+ ok <- unsafeRead old ri+ if ok < lim+ then do+ !ov <- unsafeRead old (ri+1)+ unsafeWrite new wi ok+ unsafeWrite new (wi+1) ov+ go (ri+2) (wi+2)+ else return (wi,ri)+ | otherwise = return (wi,ri)+ go oi ni++{-# INLINE copyRem #-}+copyRem :: Int -> STUArray s Int COUNT_T -> Int -> STUArray s Int COUNT_T -> Int -> ST s Int+copyRem end old oi new ni = do+ let go ri wi+ | ri < end = do+ unsafeRead old ri >>= unsafeWrite new wi+ go (ri+1) (wi+1)+ | otherwise = return wi+ go oi ni++{-# INLINE cp6 #-}+cp6 :: COUNT_T -> Integer+cp6 k =+ case k `quotRem` 30030 of+ (q,r) -> 5760*fromIntegral q ++ fromIntegral (cpCtAr `unsafeAt` fromIntegral r)++cop :: COUNT_T -> COUNT_T+cop m = m - fromIntegral (cpDfAr `unsafeAt` fromIntegral (m `rem` 30030))+++--------------------------------------------------------------------------------+-- Ugly helper arrays --+--------------------------------------------------------------------------------++cpCtAr :: UArray Int Int16+cpCtAr = runSTUArray $ do+ ar <- newArray (0,30029) 1+ let zilch s i+ | i < 30030 = unsafeWrite ar i 0 >> zilch s (i+s)+ | otherwise = return ()+ accumulate ct i+ | i < 30030 = do+ v <- unsafeRead ar i+ let !ct' = ct+v+ unsafeWrite ar i ct'+ accumulate ct' (i+1)+ | otherwise = return ar+ zilch 2 0+ zilch 6 3+ zilch 10 5+ zilch 14 7+ zilch 22 11+ zilch 26 13+ accumulate 1 2++cpDfAr :: UArray Int Int8+cpDfAr = runSTUArray $ do+ ar <- newArray (0,30029) 0+ let note s i+ | i < 30029 = unsafeWrite ar i 1 >> note s (i+s)+ | otherwise = return ()+ accumulate d i+ | i < 30029 = do+ v <- unsafeRead ar i+ if v == 0+ then accumulate 2 (i+2)+ else do unsafeWrite ar i d+ accumulate (d+1) (i+1)+ | otherwise = return ar+ note 2 0+ note 6 3+ note 10 5+ note 14 7+ note 22 11+ note 26 13+ accumulate 2 3++cpGpAr :: UArray Int Int8+cpGpAr = runSTUArray $ do+ ar <- newArray (0,30030) 0+ unsafeWrite ar 30030 1+ let note s i+ | i < 30029 = unsafeWrite ar i 1 >> note s (i+s)+ | otherwise = return ()+ accumulate d i+ | i < 1 = return ar+ | otherwise = do+ v <- unsafeRead ar i+ if v == 0+ then accumulate 2 (i-2)+ else do unsafeWrite ar i d+ accumulate (d+1) (i-1)+ | otherwise = return ar+ note 2 0+ note 6 3+ note 10 5+ note 14 7+ note 22 11+ note 26 13+ accumulate 2 30027+
+ Math/NumberTheory/Primes/Factorisation.hs view
@@ -0,0 +1,166 @@+-- |+-- Module: Math.NumberTheory.Primes.Factorisation+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Various functions related to prime factorisation.+-- Many of these functions use the prime factorisation of an 'Integer'.+-- If several of them are used on the same 'Integer', it would be inefficient+-- to recalculate the factorisation, hence there are also functions working+-- on the canonical factorisation, these require that the number be positive+-- and in the case of the Carmichael function that the list of prime factors+-- with their multiplicities is ascending.+module Math.NumberTheory.Primes.Factorisation+ ( -- * Factorisation functions+ -- $algorithm+ -- ** Complete factorisation+ factorise+ , defaultStdGenFactorisation+ , stepFactorisation+ , factorise'+ , defaultStdGenFactorisation'+ -- *** Factor sieves+ , FactorSieve+ , factorSieve+ , sieveFactor+ -- ** Partial factorisation+ , smallFactors+ , stdGenFactorisation+ , curveFactorisation+ -- *** Single curve worker+ , montgomeryFactorisation+ -- * Totients+ , totient+ , φ+ , TotientSieve+ , totientSieve+ , sieveTotient+ , totientFromCanonical+ -- * Carmichael function+ , carmichael+ , λ+ , CarmichaelSieve+ , carmichaelSieve+ , sieveCarmichael+ , carmichaelFromCanonical+ -- * Divisors+ , divisors+ , tau+ , τ+ , divisorCount+ , divisorSum+ , sigma+ , σ+ , divisorPowerSum+ , divisorsFromCanonical+ , tauFromCanonical+ , divisorSumFromCanonical+ , sigmaFromCanonical+ ) where++import Data.Set (Set, singleton)++import Math.NumberTheory.Primes.Factorisation.Utils+import Math.NumberTheory.Primes.Factorisation.Montgomery+import Math.NumberTheory.Primes.Sieve.Misc++-- $algorithm+--+-- Factorisation of 'Integer's by the elliptic curve algorithm after Montgomery.+-- The algorithm is explained at+-- <http://programmingpraxis.com/2010/04/23/modern-elliptic-curve-factorization-part-1/>+-- and+-- <http://programmingpraxis.com/2010/04/27/modern-elliptic-curve-factorization-part-2/>+--+-- The implementation is not very optimised, so it is not suitable for factorising numbers+-- with several huge prime divisors. However, factors of 20-25 digits are normally found in+-- acceptable time. The time taken depends, however, strongly on how lucky the curve-picking+-- is. With luck, even large factors can be found in seconds; on the other hand, finding small+-- factors (about 12-15 digits) can take minutes when the curve-picking is bad.+--+-- Given enough time, the algorithm should be able to factor numbers of 100-120 digits, but it+-- is best suited for numbers of up to 50-60 digits.++-- | Calculates the totient of a positive number @n@, i.e.+-- the number of @k@ with @1 <= k <= n@ and @'gcd' n k == 1@,+-- in other words, the order of the group of units in @ℤ/(n)@.+totient :: Integer -> Integer+totient n+ | n < 1 = error "Totient only defined for positive numbers"+ | n == 1 = 1+ | otherwise = totientFromCanonical (factorise' n)++-- | Alias of 'totient' for people who prefer Greek letters.+φ :: Integer -> Integer+φ = totient++-- | Calculates the Carmichael function for a positive integer, that is,+-- the (smallest) exponent of the group of units in @&8484;/(n)@.+carmichael :: Integer -> Integer+carmichael n+ | n < 1 = error "Carmichael function only defined for positive numbers"+ | n == 1 = 1+ | otherwise = carmichaelFromCanonical (factorise' n)++-- | Alias of 'carmichael' for people who prefer Greek letters.+λ :: Integer -> Integer+λ = carmichael++-- | @'divisors' n@ is the set of all (positive) divisors of @n@.+-- @'divisors' 0@ is an error because we can't create the set of all 'Integer's.+divisors :: Integer -> Set Integer+divisors n+ | n < 0 = divisors (-n)+ | n == 0 = error "Can't create set of divisors of 0"+ | n == 1 = singleton 1+ | otherwise = divisorsFromCanonical (factorise' n)++-- | @'tau' n@ is the number of (positive) divisors of @n@.+-- @'tau' 0@ is an error because @0@ has infinitely many divisors.+tau :: Integer -> Integer+tau n+ | n < 0 = tau (-n)+ | n == 0 = error "0 has infinitely many divisors"+ | n == 1 = 1+ | otherwise = tauFromCanonical (factorise' n)++-- | Alias for 'tau'.+divisorCount :: Integer -> Integer+divisorCount = tau++-- | The sum of all (positive) divisors of a positive number @n@,+-- calculated from its prime factorisation.+divisorSum :: Integer -> Integer+divisorSum n+ | n < 1 = error "divisor sum only defined for positive numbers"+ | n == 1 = 1+ | otherwise = divisorSumFromCanonical (factorise' n)++-- | Alias for 'sigma'.+divisorPowerSum :: Int -> Integer -> Integer+divisorPowerSum = sigma++-- | @'sigma' k n@ is the sum of the @k@-th powers of the+-- (positive) divisors of @n@. @k@ must be non-negative and @n@ positive.+-- For @k == 0@, it is the divisor count (@d^0 = 1@).+sigma :: Int -> Integer -> Integer+sigma 0 n = tau n+sigma 1 n = divisorSum n+sigma k n+ | k < 0 = error "sigma: exponent must be non-negative"+ | n < 1 = error "sigma: n must be positive"+ | n == 1 = 1+ | otherwise = sigmaFromCanonical k (factorise' n)++-- | Alias for 'sigma' for people preferring Greek letters.+σ :: Int -> Integer -> Integer+σ 0 = divisorCount+σ 1 = divisorSum+σ k = divisorPowerSum k++-- | Alias for 'tau' for people preferring Greek letters.+τ :: Integer -> Integer+τ = tau
+ Math/NumberTheory/Primes/Factorisation/Montgomery.hs view
@@ -0,0 +1,399 @@+-- |+-- Module: Math.NumberTheory.Primes.Factorisation.Montgomery+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Factorisation of 'Integer's by the elliptic curve algorithm after Montgomery.+-- The algorithm is explained at+-- <http://programmingpraxis.com/2010/04/23/modern-elliptic-curve-factorization-part-1/>+-- and+-- <http://programmingpraxis.com/2010/04/27/modern-elliptic-curve-factorization-part-2/>+--+-- The implementation is not very optimised, so it is not suitable for factorising numbers+-- with only huge prime divisors. However, factors of 20-25 digits are normally found in+-- acceptable time. The time taken depends, however, strongly on how lucky the curve-picking+-- is. With luck, even large factors can be found in seconds; on the other hand, finding small+-- factors (about 10 digits) can take minutes when the curve-picking is bad.+--+-- Given enough time, the algorithm should be able to factor numbers of 100-120 digits, but it+-- is best suited for numbers of up to 50-60 digits.+{-# LANGUAGE CPP, BangPatterns, MagicHash #-}+{-# OPTIONS_HADDOCK hide #-}+module Math.NumberTheory.Primes.Factorisation.Montgomery+ ( -- * Complete factorisation functions+ -- ** Functions with input checking+ factorise+ , defaultStdGenFactorisation+ -- ** Functions without input checking+ , factorise'+ , stepFactorisation+ , defaultStdGenFactorisation'+ -- * Partial factorisation+ , smallFactors+ , stdGenFactorisation+ , curveFactorisation+ -- ** Single curve worker+ , montgomeryFactorisation+ ) where++#include "MachDeps.h"++import GHC.Base+import GHC.Word+import Data.Array.Base++import System.Random+import Control.Monad.State.Strict+import Control.Applicative+import Data.Bits+import Data.Maybe++import Math.NumberTheory.Logarithms+import Math.NumberTheory.Logarithms.Internal+import Math.NumberTheory.Primes.Sieve.Eratosthenes+import Math.NumberTheory.Primes.Sieve.Indexing+import Math.NumberTheory.Primes.Testing.Probabilistic+import Math.NumberTheory.Utils++-- | @'factorise' n@ produces the prime factorisation of @n@, including+-- a factor of @(-1)@ if @n < 0@. @'factorise' 0@ is an error and the+-- factorisation of @1@ is empty. Uses a 'StdGen' produced in an arbitrary+-- manner from the bit-pattern of @n@.+factorise :: Integer -> [(Integer,Int)]+factorise n+ | n < 0 = (-1,1):factorise (-n)+ | n == 0 = error "0 has no prime factorisation"+ | n == 1 = []+ | otherwise = factorise' n++-- | Like 'factorise', but without input checking, hence @n > 1@ is required.+factorise' :: Integer -> [(Integer,Int)]+factorise' n = defaultStdGenFactorisation' (mkStdGen $ fromInteger n `xor` 0xdeadbeef) n++-- | @'stepFactorisation'@ is like 'factorise'', except that it doesn't use a+-- pseudo random generator but steps through the curves in order.+-- This strategy turns out to be surprisingly fast, on average it doesn't+-- seem to be slower than the 'StdGen' based variant.+stepFactorisation :: Integer -> [(Integer,Int)]+stepFactorisation n+ = let (sfs,mb) = smallFactors 100000 n+ in sfs ++ case mb of+ Nothing -> []+ Just r -> curveFactorisation (Just 10000000000) bailliePSW+ (\m k -> (if k < (m-1) then k else error "Curves exhausted",k+1)) 6 Nothing r++-- | @'defaultStdGenFactorisation'@ first strips off all small prime factors and then,+-- if the factorisation is not complete, proceeds to curve factorisation.+-- For negative numbers, a factor of @-1@ is included, the factorisation of @1@+-- is empty. Since @0@ has no prime factorisation, a zero argument causes+-- an error.+defaultStdGenFactorisation :: StdGen -> Integer -> [(Integer,Int)]+defaultStdGenFactorisation sg n+ | n == 0 = error "0 has no prime factorisation"+ | n < 0 = (-1,1) : defaultStdGenFactorisation sg (-n)+ | n == 1 = []+ | otherwise = defaultStdGenFactorisation' sg n++-- | Like 'defaultStdGenFactorisation', but without input checking, so+-- @n@ must be larger than @1@.+defaultStdGenFactorisation' :: StdGen -> Integer -> [(Integer,Int)]+defaultStdGenFactorisation' sg n+ = let (sfs,mb) = smallFactors 100000 n+ in sfs ++ case mb of+ Nothing -> []+ Just m -> stdGenFactorisation (Just 10000000000) sg Nothing m++----------------------------------------------------------------------------------------------------+-- Factorisation wrappers --+----------------------------------------------------------------------------------------------------++-- | A wrapper around 'curveFactorisation' providing a few default arguments.+-- The primality test is 'bailliePSW', the @prng@ function - naturally -+-- 'randomR'. This function also requires small prime factors to have been+-- stripped before.+stdGenFactorisation :: Maybe Integer -- ^ Lower bound for composite divisors+ -> StdGen -- ^ Standard PRNG+ -> Maybe Int -- ^ Estimated number of digits of smallest prime factor+ -> Integer -- ^ The number to factorise+ -> [(Integer,Int)] -- ^ List of prime factors and exponents+stdGenFactorisation primeBound sg digits n+ = curveFactorisation primeBound bailliePSW (\m -> randomR (6,m-2)) sg digits n++-- | @'curveFactorisation'@ is the driver for the factorisation. Its performance (and success)+-- can be influenced by passing appropriate arguments. If you know that @n@ has no prime divisors+-- below @b@, any divisor found less than @b*b@ must be prime, thus giving @Just (b*b)@ as the+-- first argument allows skipping the comparatively expensive primality test for those.+-- If @n@ is such that all prime divisors must have a specific easy to test for structure, a+-- custom primality test can improve the performance (normally, it will make very little+-- difference, since @n@ has not many divisors, and many curves have to be tried to find one).+-- More influence has the pseudo random generator (a function @prng@ with @6 <= fst (prng k s) <= k-2@+-- and an initial state for the PRNG) used to generate the curves to try. A lucky choice here can+-- make a huge difference. So, if the default takes too long, try another one; or you can improve your+-- chances for a quick result by running several instances in parallel.+--+-- @'curveFactorisation'@ requires that small prime factors have been stripped before. Also, it is+-- unlikely to succeed if @n@ has more than one (really) large prime factor.+curveFactorisation :: Maybe Integer -- ^ Lower bound for composite divisors+ -> (Integer -> Bool) -- ^ A primality test+ -> (Integer -> g -> (Integer,g)) -- ^ A PRNG+ -> g -- ^ Initial PRNG state+ -> Maybe Int -- ^ Estimated number of digits of the smallest prime factor+ -> Integer -- ^ The number to factorise+ -> [(Integer,Int)] -- ^ List of prime factors and exponents+curveFactorisation primeBound primeTest prng seed mbdigs n+ | ptest n = [(n,1)]+ | otherwise = evalState (fact n digits) seed+ where+ digits = fromMaybe 6 mbdigs+ mult 1 xs = xs+ mult j xs = [(p,j*k) | (p,k) <- xs]+ dbl (u,v) = (mult 2 u, mult 2 v)+ ptest = case primeBound of+ Just bd -> \k -> k <= bd || primeTest k+ Nothing -> primeTest+ rndR k = state (\gen -> prng k gen)+ fact m digs = do let (b1,b2,ct) = findParms digs+ (pfs,cfs) <- repFact m b1 b2 ct+ if null cfs+ then return pfs+ else do+ 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+ | count < 0 = return ([],[(m,1)])+ | otherwise = do+ s <- rndR m+ case montgomeryFactorisation m b1 b2 s of+ Nothing -> repFact m b1 b2 (count-1)+ Just d -> do+ let !cof = m `quot` d+ case gcd cof d of+ 1 -> do+ (dp,dc) <- if ptest d+ then return ([(d,1)],[])+ else repFact d b1 b2 (count-1)+ (cp,cc) <- if ptest cof+ then return ([(cof,1)],[])+ else repFact cof b1 b2 (count-1)+ return (merge dp cp, dc ++ cc)+ g -> do+ let d' = d `quot` g+ c' = cof `quot` g+ (dp,dc) <- if ptest d'+ then return ([(d',1)],[])+ else repFact d' b1 b2 (count-1)+ (cp,cc) <- if ptest c'+ then return ([(c',1)],[])+ else repFact c' b1 b2 (count-1)+ (gp,gc) <- if ptest g+ then return ([(g,2)],[])+ else dbl <$> repFact g b1 b2 (count-1)+ return (mergeAll [dp,cp,gp], dc ++ cc ++ gc)++----------------------------------------------------------------------------------------------------+-- The workhorse --+----------------------------------------------------------------------------------------------------++-- | @'montgomeryFactorisation' n b1 b2 s@ tries to find a factor of @n@ using the+-- curve and point determined by the seed @s@ (@6 <= s < n-1@), multiplying the+-- point by the least common multiple of all numbers @<= b1@ and all primes+-- between @b1@ and @b2@. The idea is that there's a good chance that the order+-- of the point in the curve over one prime factor divides the multiplier, but the+-- order over another factor doesn't, if @b1@ and @b2@ are appropriately chosen.+-- If they are too small, none of the orders will probably divide the multiplier,+-- if they are too large, all probably will, so they should be chosen to fit+-- the expected size of the smallest factor.+--+-- It is assumed that @n@ has no small prime factors.+--+-- The result is maybe a nontrivial divisor of @n@.+montgomeryFactorisation :: Integer -> Word -> Word -> Integer -> Maybe Integer+montgomeryFactorisation n b1 b2 s = go p5 (list primeStore)+ where+ l2 = wordLog2' b1+ b1i = toInteger b1+ (^~) :: Word -> Int -> Word+ w ^~ i = w ^ i+ (e, p0) = montgomeryData n s+ dbl pt = double n e pt+ dbln 0 !pt = pt+ dbln k pt = dbln (k-1) (dbl pt)+ p2 = dbln l2 p0+#if WORD_SIZE_IN_BITS == 64+ mul a b c = (a*b) `quot` c -- can't overflow, work on Int+#else+ mul a b c = fromInteger ((toInteger a * b) `quot` c) -- might overflow if Int is used+#endif+ adjust bd ml w+ | w <= bd = adjust bd ml (w*ml)+ | otherwise = w+ l3 = mul l2 190537 301994+ w3 = 3 ^~ l3+ pw3 = adjust (b1 `quot` 3) 3 w3+ p3 = multiply n e pw3 p2+ l5 = mul l2 1936274 4495889+ w5 = 5 ^~ l5+ pw5 = adjust (b1 `quot` 5) 5 w5+ p5 = multiply n e pw5 p3+ go (P _ 0) _ = Nothing+ go !pt@(P _ z) (pr:prs)+ | pr <= b1 = let !lp = integerLogBase' (fromIntegral pr) b1i+ in go (multiply n e (pr ^~ lp) pt) prs+ | otherwise = case gcd n z of+ 1 -> lgo (multiply n e pr pt) prs+ g -> Just g+ go (P _ z) _ = case gcd n z of+ 1 -> Nothing+ g -> Just g+ lgo (P _ 0) _ = Nothing+ lgo !pt@(P _ z) (pr:prs)+ | pr <= b2 = lgo (multiply n e pr pt) prs+ | otherwise = case gcd n z of+ 1 -> Nothing+ g -> Just g+ lgo (P _ z) _ = case gcd n z of+ 1 -> Nothing+ g -> Just g++----------------------------------------------------------------------------------------------------+-- Helpers, Curves and elliptic arithmetics --+----------------------------------------------------------------------------------------------------++-- A Montgomery curve is given by y^2 = x^3 + (A_n / A_d - 2)*x^2 + x (mod n).+-- We store A_n and 4*A_d, since A_n occurs with the factor 4 in all formulae.+data Curve = C !Integer !Integer++-- Point in the projective plane, will be on the curve+-- A coordinate transformation eliminates the y-coordinate, hence+-- we store only x and z+data Point = P !Integer !Integer++-- Get curve and point to start+-- Input should satisfy 6 <= s < n-1+montgomeryData :: Integer -> Integer -> (Curve, Point)+montgomeryData n s = (C an ad4, P x z)+ where+ u = (s*s-5) `mod` n+ v = (4*s) `mod` n+ d = (v-u)+ x = (u*u*u) `mod` n+ z = (v*v*v) `mod` n+ an = ((d*d)*(d*(3*u+v))) `mod` n+ ad4 = (16*x*v) `mod` n++-- Addition on the curve, given the modulus n and three points,+-- p0, p1 and p2, with p0 = p2 - p1, calculate the point p1 + p2.+-- Note that the addition does not depend on the curve.+add :: Integer -> Point -> Point -> Point -> Point+add n (P x0 z0) (P x1 z1) (P x2 z2) = P x3 z3+ where+ a = (x1-z1)*(x2+z2)+ b = (x1+z1)*(x2-z2)+ c = a+b+ d = a-b+ x3 = (c*c*z0) `rem` n+ z3 = (d*d*x0) `rem` n++-- Double a point on the curve.+double :: Integer -> Curve -> Point -> Point+double n (C an ad4) (P x z) = P x' z'+ where+ r = x+z+ s = x-z+ u = r*r+ v = s*s+ t = u-v+ x' = (ad4*u*v) `rem` n+ z' = ((ad4*v+t*an)*t) `rem` n++-- Multiply a point on the curve by a Word.+-- Within Word range, we can use the faster variant going+-- from high-order bits to low-order.+multiply :: Integer -> Curve -> Word -> Point -> Point+multiply n cve (W# w##) p =+ case wordLog2# w## of+ l# -> go (l# -# 1#) p (double n cve p)+ where+ go 0# !p0 !p1 = case w## `and#` 1## of+ 0## -> double n cve p0+ _ -> add n p p0 p1+ go i# p0 p1 = case (uncheckedShiftRL# w## i#) `and#` 1## of+ 0## -> go (i# -# 1#) (double n cve p0) (add n p p0 p1)+ _ -> go (i# -# 1#) (add n p p0 p1) (double n cve p1)++{- Not (yet) needed+-- Multiply a point on the curve by an Integer.+multIgr :: Integer -> Curve -> Integer -> Point -> Point+multIgr n cve k p = go k+ where+ go 1 = (p, double n cve p)+ go m = case m `quotRem` 2 of+ (q,r) -> let !(!s, l) = go q+ in case r of+ 0 -> (double n cve s, add n p s l)+ _ -> (add n p s l, double n cve l)+-}++-- primes, compactly stored as a bit sieve+primeStore :: [PrimeSieve]+primeStore = psieveFrom 7++-- generate list of primes from arrays+list :: [PrimeSieve] -> [Word]+list sieves = concat [[off + toPrim i | i <- [0 .. li], unsafeAt bs i]+ | PS vO bs <- sieves, let { (_,li) = bounds bs; off = fromInteger vO; }]++-- | @'smallFactors' bound n@ finds all prime divisors of @n > 1@ up to @bound@ by trial division and returns the+-- list of these together with their multiplicities, and a possible remaining factor which may be composite.+smallFactors :: Integer -> Integer -> ([(Integer,Int)], Maybe Integer)+smallFactors bd n = case shiftToOddCount n of+ (0,m) -> go m prms+ (k,m) -> (2,k) <: if m == 1 then ([],Nothing) else go m prms+ where+ prms = tail (primeStore >>= primeList)+ x <: ~(l,b) = (x:l,b)+ go m (p:ps)+ | m < p*p = ([(m,1)], Nothing)+ | bd < p = ([], Just m)+ | otherwise = case splitOff p m of+ (0,_) -> go m ps+ (k,r) | r == 1 -> ([(p,k)], Nothing)+ | otherwise -> (p,k) <: go r ps+ go m [] = ([(m,1)], Nothing)++-- helpers: merge sorted lists+merge :: [(Integer,Int)] -> [(Integer,Int)] -> [(Integer,Int)]+merge xxs@(x@(p,k):xs) yys@(y@(q,m):ys) = case compare p q of+ LT -> x : merge xs yys+ EQ -> (p,k+m) : merge xs ys+ GT -> y : merge xxs ys+merge xs [] = xs+merge _ ys = ys++mergeAll :: [[(Integer,Int)]] -> [(Integer,Int)]+mergeAll [] = []+mergeAll [xs] = xs+mergeAll (xs:ys:zss) = merge (merge xs ys) (mergeAll zss)++-- Parameters for the factorisation, the two b-parameters for montgomery and the number of tries+-- to use these, depending on the size of the factor we are looking for.+-- The numbers are roughly based on the parameters listed on Dario Alpern's ECM site.+testParms :: [(Int,Word,Word,Int)]+testParms = [ (12, 400, 10000, 10), (15, 2000, 50000, 25), (20, 11000, 150000, 90)+ , (25, 50000, 500000, 300), (30, 250000, 1500000, 700)+ , (35, 1000000, 4000000, 1800), (40, 3000000, 12000000, 5100)+ , (45, 11000000, 45000000, 10600), (50, 43000000, 200000000, 19300)+ , (55, 80000000, 400000000,30000), (60, 120000000, 700000000, 50000)+ ]++findParms :: Int -> (Word, Word, Int)+findParms digs = go (100, 1000, 7) testParms+ where+ go p ((d,b1,b2,ct):rest)+ | digs < d = p+ | otherwise = go (b1,b2,ct) rest+ go p [] = p
+ Math/NumberTheory/Primes/Factorisation/Utils.hs view
@@ -0,0 +1,77 @@+-- |+-- Module: Math.NumberTheory.Primes.Factorisation.Utils+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Some utilities related to factorisation, defined here to avoid import cycles.+{-# LANGUAGE BangPatterns #-}+{-# OPTIONS_HADDOCK hide #-}+module Math.NumberTheory.Primes.Factorisation.Utils+ ( ppTotient+ , totientFromCanonical+ , carmichaelFromCanonical+ , divisorsFromCanonical+ , tauFromCanonical+ , divisorSumFromCanonical+ , sigmaFromCanonical+ ) where++import Data.Set (Set)+import qualified Data.Set as Set+import Data.Bits+import Data.List++import Math.NumberTheory.Powers.Integer++-- | Totient of a prime power.+ppTotient :: (Integer,Int) -> Integer+ppTotient (p,1) = p-1+ppTotient (p,k) = (p-1)*(integerPower p (k-1)) -- slightly faster than (^) usually++-- | Calculate the totient from the canonical factorisation.+totientFromCanonical :: [(Integer,Int)] -> Integer+totientFromCanonical = product . map ppTotient++-- | Calculate the Carmichael function from the factorisation.+-- Requires that the list of prime factors is strictly ascending.+carmichaelFromCanonical :: [(Integer,Int)] -> Integer+carmichaelFromCanonical = go2+ where+ go2 ((2,k):ps) = let acc = case k of+ 1 -> 1+ 2 -> 2+ _ -> 1 `shiftL` (k-2)+ in go acc ps+ go2 ps = go 1 ps+ go !acc ((p,1):pps) = go (lcm acc (p-1)) pps+ go acc ((p,k):pps) = go ((lcm acc (p-1))*integerPower p (k-1)) pps+ go acc [] = acc++-- | The set of divisors, efficiently calculated from the canonical factorisation.+divisorsFromCanonical :: [(Integer,Int)] -> Set Integer+divisorsFromCanonical = foldl' step (Set.singleton 1)+ where+ step st (p,k) = Set.unions (st:[Set.mapMonotonic (*pp) st | pp <- take k (iterate (*p) p) ])++-- | The number of divisors, efficiently calculated from the canonical factorisation.+tauFromCanonical :: [(a,Int)] -> Integer+tauFromCanonical pps = product [fromIntegral k + 1 | (_,k) <- pps]++-- | The sum of all divisors, efficiently calculated from the canonical factorisation.+divisorSumFromCanonical :: [(Integer,Int)] -> Integer+divisorSumFromCanonical = product . map ppDivSum++ppDivSum :: (Integer,Int) -> Integer+ppDivSum (p,1) = p+1+ppDivSum (p,k) = (p^(k+1)-1) `quot` (p-1)++-- | The sum of the powers (with fixed exponent) of all divisors,+-- efficiently calculated from the canonical factorisation.+sigmaFromCanonical :: Int -> [(Integer,Int)] -> Integer+sigmaFromCanonical k = product . map (ppDivPowerSum k)++ppDivPowerSum :: Int -> (Integer,Int) -> Integer+ppDivPowerSum k (p,m) = (p^(k*(m+1)) - 1) `quot` (p^k - 1)
+ Math/NumberTheory/Primes/Heap.hs view
@@ -0,0 +1,388 @@+-- |+-- Module: Math.NumberTheory.Primes.Heap+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Prime generation using a priority queue for the composites.+-- The algorithm is basically the one described in+-- <http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf>, but+-- it uses a more efficient heap for the priority queue and a+-- larger wheel, thus it is faster (in particular, the speed+-- penalty for @'Integer'@ is much smaller) and uses less memory.+-- It is nevertheless very slow compared to a bit sieve.+-- This module is mainly intended for comparison and verification.+{-# LANGUAGE BangPatterns, CPP, MonoLocalBinds #-}+{-# OPTIONS_GHC -funbox-strict-fields #-}+#if __GLASGOW_HASKELL__ >= 700+{-# OPTIONS_GHC -fno-float-in -fno-spec-constr -fno-full-laziness #-}+#endif+module Math.NumberTheory.Primes.Heap (primes, sieveFrom) where++import Data.Array.Unboxed+import Data.Array.Base (unsafeAt, unsafeRead, unsafeWrite)+import Data.Array.ST+import Control.Monad.ST+import Data.List (foldl')+import Data.Word++#ifndef SH_SIZE+#define SH_SIZE 31+#endif++-- Composites to eliminate, components are+-- composite, multiple of prime+-- prime+-- index of step to find next multiple of prime+data Del a = D !a !a {-# UNPACK #-} !Int++step :: Integral a => Int -> a+-- {-# INLINE step #-}+step i = fromIntegral (steps `unsafeAt` i)++-- Priority queue as baby heap+-- Invariant: left subheap one larger than right or both+-- have the same size (and of course, heap property)+data Hipp a+ = E+ | H !a !a {-# UNPACK #-} !Int !(Hipp a) !(Hipp a)+++-- push composite-data down the heap+{-# SPECIALISE push :: Int -> Int -> Int -> Hipp Int -> Hipp Int #-}+{-# SPECIALISE push :: Word -> Word -> Int -> Hipp Word -> Hipp Word #-}+{-# SPECIALISE push :: Integer -> Integer -> Int -> Hipp Integer -> Hipp Integer #-}+push :: Integral a => a -> a -> Int -> Hipp a -> Hipp a+-- GHC 7 does not like the old code, so it gets a new implementation.+-- That is faster than what it does with the old code, but still slower+-- than what GHC 6 did with it. :(+#if __GLASGOW_HASKELL__ >= 700+push !c !p !w = go+ where+ less = (< c)+ go (H hc hp hw l r)+ | less hc = H hc hp hw (go r) l+ | otherwise = H c p w (push hc hp hw r) l+ go _ = H c p w E E+#else+push c p w (H hc hp hw l r)+ | c < hc = H c p w (push hc hp hw r) l+ | otherwise = H hc hp hw (push c p w r) l+push c p w E = H c p w E E+#endif++-- bubble down increased top to regain heap invariant+{-# SPECIALISE bubble :: Hipp Int -> Hipp Int #-}+{-# SPECIALISE bubble :: Hipp Word -> Hipp Word #-}+{-# SPECIALISE bubble :: Hipp Integer -> Hipp Integer #-}+bubble :: Integral a => Hipp a -> Hipp a+-- Again, GHC 6 fared better, so new code for GHC 7, still+-- not as good as GHC 6 was.+#if __GLASGOW_HASKELL__ >= 700+bubble h@(H c p w l r) =+ case r of+ E -> case l of+ E -> h+ H lc lp lw ll lr+ | lc < c -> H lc lp lw (H c p w ll lr) r+ | otherwise -> h+ H rc rp rw rl rr ->+ case l of+ H lc lp lw ll lr+ | lc < c -> if lc < rc+ then H lc lp lw (mkHipp c p w ll lr) r+ else H rc rp rw l (mkHipp c p w rl rr)+ | rc < c -> H rc rp rw l (mkHipp c p w rl rr)+ | otherwise -> h+ _ -> error "Heap invariant violated, left smaller than right!"+#else+bubble h@(H c p w l@(H lc lp lw ll lr) r@(H rc rp rw rl rr))+ | c <= lc && c <= rc = h+ | lc < rc = H lc lp lw (bubble (H c p w ll lr)) r+ | otherwise = H rc rp rw l (bubble (H c p w rl rr))+bubble h@(H c p w (H lc lp lw _ _) _)+ | c <= lc = h+ | otherwise = H lc lp lw (H c p w E E) E+#endif+bubble h = h++-- join two heaps and composite-data, GHC 7 doesn't do well on the old bubble.+#if __GLASGOW_HASKELL__ >= 700+{-# SPECIALISE+ mkHipp :: Int -> Int -> Int -> Hipp Int -> Hipp Int -> Hipp Int,+ Integer -> Integer -> Int -> Hipp Integer -> Hipp Integer -> Hipp Integer,+ Word -> Word -> Int -> Hipp Word -> Hipp Word -> Hipp Word+ #-}+mkHipp :: Integral a => a -> a -> Int -> Hipp a -> Hipp a -> Hipp a+mkHipp !c !p !w = go+ where+ less = (< c)+ go l r =+ case r of+ E -> case l of+ E -> H c p w l r+ H lc lp lw _ _+ | less lc -> H lc lp lw (H c p w E E) E+ | otherwise -> H c p w l r+ H rc rp rw rl rr ->+ case l of+ H lc lp lw ll lr+ | less lc -> if lc < rc+ then H lc lp lw (go ll lr) r+ else H rc rp rw l (go rl rr)+ | less rc -> H rc rp rw l (go rl rr)+ | otherwise -> H c p w l r+ _ -> error "Heap invariant violated, left smaller than right!"+-- {-# INLINE mkHipp #-}+#endif++-- increase the top of the heap and re-heap+{-# SPECIALISE inc :: Hipp Int -> Hipp Int #-}+{-# SPECIALISE inc :: Hipp Word -> Hipp Word #-}+{-# SPECIALISE inc :: Hipp Integer -> Hipp Integer #-}+inc :: Integral a => Hipp a -> Hipp a+inc (H c p i l r)+ = {-# SCC "incBubble" #-} bubble (H (c+p*step i) p (nextIndex i) l r)+inc h = h++-- while top of heap equals composite, increase and re-heap+{-# SPECIALISE adjust :: Int -> Hipp Int -> Hipp Int #-}+{-# SPECIALISE adjust :: Word -> Hipp Word -> Hipp Word #-}+{-# SPECIALISE adjust :: Integer -> Hipp Integer -> Hipp Integer #-}+adjust :: Integral a => a -> Hipp a -> Hipp a+adjust cm h@(H v _ _ _ _)+ | cm == v = adjust cm (inc h)+adjust _ h = h++-- build a heap from a sorted list of Del's+{-# SPECIALISE buildH :: [Del Int] -> Hipp Int #-}+{-# SPECIALISE buildH :: [Del Word] -> Hipp Word #-}+{-# SPECIALISE buildH :: [Del Integer] -> Hipp Integer #-}+buildH :: Integral a => [Del a] -> Hipp a+buildH [] = E+buildH (D s p w : tl) = H s p w l r+ where+ (ll,rl) = goSplit [] [] tl+ goSplit xs ys [] = (reverse ys, reverse xs)+ goSplit xs ys (d:ds) = goSplit ys (d:xs) ds+ l = buildH ll+ r = buildH rl++-- Simple sieve pushing each prime immediately onto the heap,+-- feeds the feeder, runs at about fourth root of the main sieve.+{-# SPECIALISE simpleSieve :: Hipp Int -> Int -> Int -> [Del Int] #-}+{-# SPECIALISE simpleSieve :: Hipp Word -> Word -> Int -> [Del Word] #-}+{-# SPECIALISE simpleSieve :: Hipp Integer -> Integer -> Int -> [Del Integer] #-}+simpleSieve :: Integral a => Hipp a -> a -> Int -> [Del a]+simpleSieve h@(H nc _ _ _ _) cd !i+ | cd < nc = D s cd i : simpleSieve ({-# SCC "simplePush" #-} push s cd i h) (cd + step i) (nextIndex i)+ | otherwise = simpleSieve (adjust cd h) (cd + step i) (nextIndex i)+ where+ s = cd*cd+simpleSieve _ _ _ = [] -- would violate an invariant++-- Feeder sieve, produces composites at the rate of the progress of the main sieve,+-- hence primes at about the square root of it, thus needs about fourth root heap+-- space. The two-step feeding makes the feeder produce faster and hence the main+-- sieve (since we have only one O(n^0.5) heap and not two.+-- Using two heaps, one small for multiples of small primes which change often+-- and one for multiples of larger primes which are less frequently updated+-- speeds things up.+{-# SPECIALISE feederSieve :: [Del Int] -> Hipp Int -> Hipp Int -> Int -> Int -> [Del Int] #-}+{-# SPECIALISE feederSieve :: [Del Word] -> Hipp Word -> Hipp Word -> Word -> Int -> [Del Word] #-}+{-# SPECIALISE feederSieve :: [Del Integer] -> Hipp Integer -> Hipp Integer -> Integer -> Int -> [Del Integer] #-}+feederSieve :: Integral a => [Del a] -> Hipp a -> Hipp a -> a -> Int -> [Del a]+feederSieve dls@((D s p u):ds) sh@(H sc _ _ _ _) lh@(H lc _ _ _ _) cd i+ | cd == sc = feederSieve dls (adjust cd (inc sh)) (adjust cd lh) cd' j+ | cd == lc = feederSieve dls sh (adjust cd (inc lh)) cd' j+ | cd == s = feederSieve ds sh (push (s + p*step u) p (nextIndex u) lh) cd' j+ | otherwise = D (cd*cd) cd i : feederSieve dls sh lh cd' j+ where+ !cd' = cd + step i+ !j = nextIndex i+feederSieve _ _ _ _ _ = [] -- invariant violated++-- Build the feeder sieve, arguments are+-- first prime whose multiples have to be eliminated+-- index of step for this prime.+{-# SPECIALISE feeder :: Int -> Int -> [Del Int] #-}+{-# SPECIALISE feeder :: Word -> Int -> [Del Word] #-}+{-# SPECIALISE feeder :: Integer -> Int -> [Del Integer] #-}+feeder :: Integral a => a -> Int -> [Del a]+feeder p i = feederSieve lrg sh lh p i+ where+ (sml,D s lp w : lrg) = splitAt SH_SIZE (D q p i : {-# SCC "simple" #-} simpleSieve (H q p i E E) (p+step i) (nextIndex i))+ sh = buildH sml+ lh = H s lp w E E+ q = p*p++-- The main sieve. Code almost identical to feederSieve, but we don't construct the Del,+-- which gains some performance.+{-# SPECIALISE primeSieve :: [Del Int] -> Hipp Int -> Hipp Int -> Int -> Int -> [Int] #-}+{-# SPECIALISE primeSieve :: [Del Word] -> Hipp Word -> Hipp Word -> Word -> Int -> [Word] #-}+{-# SPECIALISE primeSieve :: [Del Integer] -> Hipp Integer -> Hipp Integer -> Integer -> Int -> [Integer] #-}+primeSieve :: Integral a => [Del a] -> Hipp a -> Hipp a -> a -> Int -> [a]+primeSieve dls@((D s p u):ds) sh@(H sc _ _ _ _) lh@(H lc _ _ _ _) cd i+ | cd == sc = primeSieve dls ({-# SCC "adjSmall" #-} adjust cd (inc sh)) ({-# SCC "adjLarge" #-}adjust cd lh) cd' j+ | cd == lc = primeSieve dls sh (adjust cd (inc lh)) cd' j+ | cd == s = primeSieve ds sh (push (s + p*step u) p (nextIndex u) lh) cd' j+ | otherwise = cd : primeSieve dls sh lh cd' j+ where+ !cd' = cd + step i+ !j = nextIndex i+primeSieve _ _ _ _ _ = [] -- invariant violated++-- | A list of primes. The sieve does not handle overflow, hence for+-- bounded types, garbage occurs near @'maxBound'@. If primes that+-- large are requested, use+--+-- @+-- 'map' 'fromInteger' $ 'takeWhile' (<= 'fromIntegral' 'maxBound') 'primes'+-- @+--+-- instead. Checking for overflow would be slower. The sieve is specialised+-- for @'Int'@, @'Word'@ and @'Integer'@, since these are the most commonly+-- used. For the fixed-width @Int@ or @Word@ types, sieving at one of the+-- specialised types and converting is faster.+-- To ensure sharing of the list of primes, bind it to a monomorphic variable,+-- to make sure that it is not shared, use @'sieveFrom'@ with different+-- arguments.+{-# SPECIALISE primes :: [Int] #-}+{-# SPECIALISE primes :: [Word] #-}+{-# SPECIALISE primes :: [Integer] #-}+primes :: Integral a => [a]+primes = 2:3:5:7:11:13:sieve 17 0++-- | @'sieveFrom' n@ generates the list of primes @>= n@.+-- The remarks about overflow and performance in the documentation+-- of @'primes'@ apply here too.+{-# SPECIALISE sieveFrom :: Int -> [Int] #-}+{-# SPECIALISE sieveFrom :: Word -> [Word] #-}+{-# SPECIALISE sieveFrom :: Integer -> [Integer] #-}+sieveFrom :: Integral a => a -> [a]+sieveFrom from+ | fromIntegral from < (32768 :: Integer)+ = dropWhile (< from) (foldr ((:) . fromIntegral) (sieve sp si) wheelPrimes)+ | otherwise+ = primeSieve dls sh lh start (nextIndex i0)+ where+ -- trick the compiler into not CAFing feeder 17 0+ sp | odd from = 17+ | otherwise = fromIntegral (remainders `unsafeAt` 0)+ si | even from = 0+ | otherwise = (steps `unsafeAt` 0)-2+ (q, r) = (from - 18) `quotRem` 30030+ i0 = findIx (fromIntegral r + 17)+ -- last number coprime to all wheel primes < from+ before = 30030*q + fromIntegral (remainders `unsafeAt` i0)+ -- first candidate+ !start = before + step i0+ (sml, lrg) = splitAt SH_SIZE (feeder sp si)+ !sh = foldl' pushD E [findMulIx p | D _ p _ <- sml]+ (lh, dls) = {-# SCC "munch" #-} munch E lrg+ pushD h (c, p, i) = push c p i h+ findMulIx p = ((p*mp), p, (nextIndex ip))+ where+ fpq = before `quot` p+ (qq, qr) = (fpq-17) `quotRem` 30030+ !ip = findIx (fromIntegral qr + 17)+ !mp = 30030*qq + fromIntegral (remainders `unsafeAt` ip) + step ip+ munch !h dels@(D s p _ : ds)+ | before < s = (h,dels)+ | otherwise = munch h' ds+ where+ !(!c, pr, i) = findMulIx p+ h' = push c pr i h+ munch h [] = (h,[])++-- Build main sieve.+{-# SPECIALISE sieve :: Int -> Int -> [Int] #-}+{-# SPECIALISE sieve :: Word -> Int -> [Word] #-}+{-# SPECIALISE sieve :: Integer -> Int -> [Integer] #-}+sieve :: Integral a => a -> Int -> [a]+sieve p i = primeSieve lrg sh lh p i+ where+ (sml,D s lp j : lrg) = splitAt SH_SIZE (feeder p i)+ !sh = buildH sml+ lh = H s lp j E E++-- next step index, we have 5760 numbers coprime to all wheel+-- primes in [1 .. product wheelPrimes]+{-# INLINE nextIndex #-}+nextIndex :: Int -> Int+nextIndex 5759 = 0+nextIndex i = i+1++-- The six smallest primes, that makes the supporting arrays small enough+-- and avoids enough composites to get acceptable speed (for sufficiently+-- generous values of acceptable).+wheelPrimes :: [Int]+wheelPrimes = 2:3:5:7:11:13:[]++-- index of largest coprime <= r+findIx :: Int -> Int+findIx r+ | 30030 < r = 5759+ | r == m = a+ | r < m = down (a-1)+ | otherwise = up a+ where+ a = min 5758 ((192*r) `quot` 1001 - 1)+ m = remainders `unsafeAt` a+ down k+ | r < (remainders `unsafeAt` k) = down (k-1)+ | otherwise = k+ up k+ | r < (remainders `unsafeAt` (k+1)) = k+ | otherwise = up (k+1)++-- array of numbers coprime to all wheel primes in wheel range+remainders :: UArray Int Int+remainders = runSTUArray $ do+ sar <- newArray (0,30029) True :: ST s (STUArray s Int Bool)+ let n2 30030 = return ()+ n2 i = unsafeWrite sar i False >> n2 (i+2)+ n3 30033 = return ()+ n3 i = unsafeWrite sar i False >> n3 (i+6)+ n5 30035 = return ()+ n5 i = unsafeWrite sar i False+ >> unsafeWrite sar (i+20) False >> n5 (i+30)+ n7 30037 = return ()+ n7 i = unsafeWrite sar i False >> n7 (i+14)+ n11 30041 = return ()+ n11 i = unsafeWrite sar i False >> n11 (i+22)+ n13 30043 = return ()+ n13 i = unsafeWrite sar i False >> n13 (i+26)+ n2 0+ n3 3+ n5 5+ n7 7+ n11 11+ n13 13+ rar <- newArray_ (0,5759) :: ST s (STUArray s Int Int)+ let loop 30031 _ = unsafeWrite rar 5759 30031 >> return rar+ loop i !r = do+ c <- unsafeRead sar i+ if c+ then do+ unsafeWrite rar r i+ loop (i+2) (r+1)+ else loop (i+2) r+ loop 17 0++-- distance from one coprime remainder to the next+steps :: UArray Int Int+steps = runSTUArray $ do+ sar <- newArray_ (0,5759) :: ST s (STUArray s Int Int)+ let loop 5759 p = do+ unsafeWrite sar 5759 (30047-p)+ return sar+ loop i p = do+ let !j = i+1+ !n = remainders `unsafeAt` j+ unsafeWrite sar i (n-p)+ loop j n+ loop 0 17+
+ Math/NumberTheory/Primes/Sieve.hs view
@@ -0,0 +1,28 @@+-- |+-- Module: Math.NumberTheory.Primes.Sieve+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Prime generation using a sieve.+-- Currently, an enhanced sieve of Eratosthenes is used, switching to an+-- Atkin sieve is planned (if I get around to implementing it and it's not slower).+--+-- The sieve used is segmented, with a chunk size chosen to give good (enough)+-- cache locality while still getting something substantial done per chunk.+-- However, that means we must store data for primes up to the square root of+-- where sieving is done, thus sieving primes up to @n@ requires+-- @/O/(sqrt n/log n)@ space.+module Math.NumberTheory.Primes.Sieve+ ( primes+ , sieveFrom+ , PrimeSieve+ , primeSieve+ , psieveList+ , psieveFrom+ , primeList+ ) where++import Math.NumberTheory.Primes.Sieve.Eratosthenes
+ Math/NumberTheory/Primes/Sieve/Eratosthenes.hs view
@@ -0,0 +1,450 @@+-- |+-- Module: Math.NumberTheory.Primes.Sieve.Eratosthenes+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Sieve+--+{-# LANGUAGE CPP, BangPatterns #-}+#if __GLASGOW_HASKELL__ >= 700+{-# OPTIONS_GHC -fspec-constr-count=6 #-}+#endif+{-# OPTIONS_HADDOCK hide #-}+module Math.NumberTheory.Primes.Sieve.Eratosthenes+ ( primes+ , sieveFrom+ , psieveFrom+ , PrimeSieve(..)+ , psieveList+ , primeList+ , primeSieve+ , nthPrimeCt+ , countFromTo+ , countAll+ , countToNth+ , sieveBytes+ , sieveBits+ , sieveWords+ , sieveRange+ , sieveTo+ ) where++#include "MachDeps.h"++import Control.Monad.ST+import Data.Array.Base (unsafeRead, unsafeWrite, unsafeAt, unsafeNewArray_)+import Data.Array.ST+import Data.Array.Unboxed+import Control.Monad (when)+import Data.Bits+import Data.Word++import Math.NumberTheory.Powers.Squares (integerSquareRoot)+import Math.NumberTheory.Utils+import Math.NumberTheory.Primes.Counting.Approximate+import Math.NumberTheory.Primes.Sieve.Indexing++-- Sieve in 128K chunks.+-- Large enough to get something done per chunk+-- and hopefully small enough to fit in the cache.+sieveBytes :: Int+sieveBytes = 128*1024++-- Number of bits per chunk.+sieveBits :: Int+sieveBits = 8*sieveBytes++-- Last index of chunk.+lastIndex :: Int+lastIndex = sieveBits - 1++-- Range of a chunk.+sieveRange :: Int+sieveRange = 30*sieveBytes++sieveWords :: Int+sieveWords = sieveBytes `quot` SIZEOF_HSWORD++#if SIZEOF_HSWORD == 8+type CacheWord = Word+#define RMASK 63+#define WSHFT 6+#define TOPB 32+#define TOPM 0xFFFFFFFF+#else+type CacheWord = Word64+#define RMASK 31+#define WSHFT 5+#define TOPB 16+#define TOPM 0xFFFF+#endif++-- | Compact store of primality flags.+data PrimeSieve = PS !Integer {-# UNPACK #-} !(UArray Int Bool)++-- | Sieve primes up to (and including) a bound.+-- For small enough bounds, this is more efficient than+-- using the segmented sieve.+--+-- Since arrays are 'Int'-indexed, overflow occurs when the sieve+-- size comes near @'maxBound' :: 'Int'@, that corresponds to an+-- upper bound near @15/8*'maxBound'@. On @32@-bit systems, that+-- is often within memory limits, so don't give bounds larger than+-- @8*10^9@ there.+primeSieve :: Integer -> PrimeSieve+primeSieve bound = PS 0 (runSTUArray $ sieveTo bound)++-- | Generate a list of primes for consumption from a+-- 'PrimeSieve'.+primeList :: PrimeSieve -> [Integer]+primeList (PS 0 bs) = 2:3:5:[toPrim i | let (lo,hi) = bounds bs+ , i <- [lo .. hi]+ , unsafeAt bs i+ ]+primeList (PS vO bs) = [vO + toPrim i+ | let (lo,hi) = bounds bs+ , i <- [lo .. hi]+ , unsafeAt bs i+ ]++-- | List of primes.+-- Since the sieve uses unboxed arrays, overflow occurs at some point,+-- but not before @10^6*'fromIntegral' ('maxBound' :: 'Int')@ (I forgot where exactly).+primes :: [Integer]+primes = 2:3:5:concat [[vO + toPrim i | i <- [0 .. li], unsafeAt bs i]+ | PS vO bs <- psieveList, let (_,li) = bounds bs]++-- | List of primes in the form of a list of 'PrimeSieve's, more compact than+-- 'primes', thus it may be better to use @'psieveList' >>= 'primeList'@+-- than keeping the list of primes alive during the entire run.+psieveList :: [PrimeSieve]+psieveList = makeSieves plim sqlim 0 0 cache+ where+ plim = 4801 -- prime #647+ sqlim = plim*plim+ cache = runSTUArray $ do+ sieve <- sieveTo 4801+ new <- unsafeNewArray_ (0,1287) :: ST s (STUArray s Int CacheWord)+ let fill j indx+ | 1279 < indx = return new -- index of 4801 = 159*30 + 31 ~> 159*8+7+ | otherwise = do+ p <- unsafeRead sieve indx+ if p+ then do+ let !i = indx .&. 7+ k :: Integer+ k = fromIntegral (indx `shiftR` 3)+ strt1 = (k*(30*k + fromIntegral (2*rho i))+ + fromIntegral (byte i)) `shiftL` 3+ + fromIntegral (idx i)+ !strt = fromIntegral strt1 .&. 0xFFFFF+ !skip = fromIntegral (strt1 `shiftR` 20)+ !ixes = fromIntegral indx `shiftL` 23 + strt `shiftL` 3 + fromIntegral i+ unsafeWrite new j skip+ unsafeWrite new (j+1) ixes+ fill (j+2) (indx+1)+ else fill j (indx+1)+ fill 0 0++makeSieves :: Integer -> Integer -> Integer -> Integer -> UArray Int CacheWord -> [PrimeSieve]+makeSieves plim sqlim bitOff valOff cache+ | valOff' < sqlim =+ let (nc, bs) = runST $ do+ cch <- unsafeThaw cache :: ST s (STUArray s Int CacheWord)+ bs0 <- slice cch+ fcch <- unsafeFreeze cch+ fbs0 <- unsafeFreeze bs0+ return (fcch, fbs0)+ in PS valOff bs : makeSieves plim sqlim bitOff' valOff' nc+ | otherwise =+ let plim' = plim + 4800+ sqlim' = plim' * plim'+ (nc,bs) = runST $ do+ cch <- growCache bitOff plim cache+ bs0 <- slice cch+ fcch <- unsafeFreeze cch+ fbs0 <- unsafeFreeze bs0+ return (fcch, fbs0)+ in PS valOff bs : makeSieves plim' sqlim' bitOff' valOff' nc+ where+ valOff' = valOff + fromIntegral sieveRange+ bitOff' = bitOff + fromIntegral sieveBits++slice :: STUArray s Int CacheWord -> ST s (STUArray s Int Bool)+slice cache = do+ hi <- snd `fmap` getBounds cache+ sieve <- newArray (0,lastIndex) True+ let treat pr+ | hi < pr = return sieve+ | otherwise = do+ w <- unsafeRead cache pr+ if w /= 0+ then unsafeWrite cache pr (w-1)+ else do+ ixes <- unsafeRead cache (pr+1)+ let !stj = ixes .&. 0x7FFFFF+ !ixw = ixes `shiftR` 23+ !i = fromIntegral (ixw .&. 7)+ !k = fromIntegral ixw - i+ !o = i `shiftL` 3+ !j = fromIntegral stj .&. 7+ !s = fromIntegral stj `shiftR` 3+ (n, u) <- tick k o j s+ let !skip = fromIntegral n `shiftR` 20+ !strt = fromIntegral n .&. 0xFFFFF+ unsafeWrite cache pr skip+ unsafeWrite cache (pr+1) (ixes - stj + strt `shiftL` 3 + fromIntegral u)+ treat (pr+2)+ tick stp off j ix+ | lastIndex < ix = return (ix - sieveBits, j)+ | otherwise = do+ p <- unsafeRead sieve ix+ when p (unsafeWrite sieve ix False)+ tick stp off ((j+1) .&. 7) (ix + stp*delta j + tau (off+j))+ treat 0++-- | Sieve up to bound in one go.+sieveTo :: Integer -> ST s (STUArray s Int Bool)+sieveTo bound = arr+ where+ (bytes,lidx) = idxPr bound+ !mxidx = 8*bytes+lidx+ mxval :: Integer+ mxval = 30*fromIntegral bytes + fromIntegral (rho lidx)+ !mxsve = integerSquareRoot mxval+ (kr,r) = idxPr mxsve+ !svbd = 8*kr+r+ arr = do+ ar <- newArray (0,mxidx) True+ let start k i = 8*(k*(30*k+2*rho i) + byte i) + idx i+ tick stp off j ix+ | mxidx < ix = return ()+ | otherwise = do+ p <- unsafeRead ar ix+ when p (unsafeWrite ar ix False)+ tick stp off ((j+1) .&. 7) (ix + stp*delta j + tau (off+j))+ sift ix+ | svbd < ix = return ar+ | otherwise = do+ p <- unsafeRead ar ix+ when p (do let i = ix .&. 7+ k = ix `shiftR` 3+ !off = i `shiftL` 3+ !stp = ix - i+ tick stp off i (start k i))+ sift (ix+1)+ sift 0++growCache :: Integer -> Integer -> UArray Int CacheWord -> ST s (STUArray s Int CacheWord)+growCache offset plim old = do+ let (_,num) = bounds old+ (bt,ix) = idxPr plim+ !start = 8*bt+ix+1+ !nlim = plim+4800+ sieve <- sieveTo nlim+ (_,hi) <- getBounds sieve+ more <- countFromToWd start hi sieve+ new <- unsafeNewArray_ (0,num+2*more) :: ST s (STUArray s Int CacheWord)+ let copy i+ | num < i = return ()+ | otherwise = do+ unsafeWrite new i (old `unsafeAt` i)+ copy (i+1)+ copy 0+ let fill j indx+ | hi < indx = return new+ | otherwise = do+ p <- unsafeRead sieve indx+ if p+ then do+ let !i = indx .&. 7+ k :: Integer+ k = fromIntegral (indx `shiftR` 3)+ strt0 = ((k*(30*k + fromIntegral (2*rho i))+ + fromIntegral (byte i)) `shiftL` 3)+ + fromIntegral (idx i)+ strt1 = strt0 - offset+ !strt = fromIntegral strt1 .&. 0xFFFFF+ !skip = fromIntegral (strt1 `shiftR` 20)+ !ixes = fromIntegral indx `shiftL` 23 + strt `shiftL` 3 + fromIntegral i+ unsafeWrite new j skip+ unsafeWrite new (j+1) ixes+ fill (j+2) (indx+1)+ else fill j (indx+1)+ fill (num+1) start++-- Danger: relies on start and end being the first resp. last+-- index in a Word+-- Do not use except in growCache and psieveFrom+{-# INLINE countFromToWd #-}+countFromToWd :: Int -> Int -> STUArray s Int Bool -> ST s Int+countFromToWd start end ba = do+ wa <- (castSTUArray :: STUArray s Int Bool -> ST s (STUArray s Int Word)) ba+ let !sb = start `shiftR` WSHFT+ !eb = end `shiftR` WSHFT+ count !acc i+ | eb < i = return acc+ | otherwise = do+ w <- unsafeRead wa i+ count (acc + bitCountWord w) (i+1)+ count 0 sb++-- count set bits between two indices (inclusive)+-- start and end must both be valid indices and start <= end+countFromTo :: Int -> Int -> STUArray s Int Bool -> ST s Int+countFromTo start end ba = do+ wa <- (castSTUArray :: STUArray s Int Bool -> ST s (STUArray s Int Word)) ba+ let !sb = start `shiftR` WSHFT+ !si = start .&. RMASK+ !eb = end `shiftR` WSHFT+ !ei = end .&. RMASK+ count !acc i+ | i == eb = do+ w <- unsafeRead wa i+ return (acc + bitCountWord (w `shiftL` (RMASK - ei)))+ | otherwise = do+ w <- unsafeRead wa i+ count (acc + bitCountWord w) (i+1)+ if sb < eb+ then do+ w <- unsafeRead wa sb+ count (bitCountWord (w `shiftR` si)) (sb+1)+ else do+ w <- unsafeRead wa sb+ let !w1 = w `shiftR` si+ return (bitCountWord (w1 `shiftL` (RMASK - ei + si)))++-- | @'sieveFrom' n@ creates the list of primes not less than @n@.+sieveFrom :: Integer -> [Integer]+sieveFrom n = case psieveFrom n of+ ps -> dropWhile (< n) (ps >>= primeList)++-- | @'psieveFrom' n@ creates the list of 'PrimeSieve's starting roughly+-- at @n@. Due to the organisation of the sieve, the list may contain+-- a few primes less than @n@.+-- This form uses less memory than @['Integer']@, hence it may be preferable+-- to use this if it is to be reused.+psieveFrom :: Integer -> [PrimeSieve]+psieveFrom n = makeSieves plim sqlim bitOff valOff cache+ where+ k0 = max 0 (n-7) `quot` 30+ valOff = 30*k0+ bitOff = 8*k0+ start = valOff+7+ ssr = integerSquareRoot (start-1) + 1+ end1 = start - 6 + fromIntegral sieveRange+ plim0 = integerSquareRoot end1+ plim = plim0 + 4801 - (plim0 `rem` 4800)+ sqlim = plim*plim+ cache = runSTUArray $ do+ sieve <- sieveTo plim+ (lo,hi) <- getBounds sieve+ pct <- countFromToWd lo hi sieve+ new <- unsafeNewArray_ (0,2*pct-1) :: ST s (STUArray s Int CacheWord)+ let fill j indx+ | hi < indx = return new+ | otherwise = do+ isPr <- unsafeRead sieve indx+ if isPr+ then do+ let !i = indx .&. 7+ !moff = i `shiftL` 3+ k :: Integer+ k = fromIntegral (indx `shiftR` 3)+ p = 30*k+fromIntegral (rho i)+ q0 = (start-1) `quot` p+ (skp0,q1) = q0 `quotRem` fromIntegral sieveRange+ (b0,r0) = idxPr (fromIntegral q1 :: Int)+ (b1,r1) | r0 == 7 = (b0+1,0)+ | otherwise = (b0,r0+1)+ b2 = skp0*fromIntegral sieveBytes + fromIntegral b1+ strt0 = ((k*(30*b2 + fromIntegral (rho r1))+ + b2 * fromIntegral (rho i)+ + fromIntegral (mu (moff + r1))) `shiftL` 3)+ + fromIntegral (nu (moff + r1))+ strt1 = ((k*(30*k + fromIntegral (2*rho i))+ + fromIntegral (byte i)) `shiftL` 3)+ + fromIntegral (idx i)+ (strt2,r2)+ | p < ssr = (strt0 - bitOff,r1)+ | otherwise = (strt1 - bitOff, i)+ !strt = fromIntegral strt2 .&. 0xFFFFF+ !skip = fromIntegral (strt2 `shiftR` 20)+ !ixes = fromIntegral indx `shiftL` 23 + strt `shiftL` 3 + fromIntegral r2+ unsafeWrite new j skip+ unsafeWrite new (j+1) ixes+ fill (j+2) (indx+1)+ else fill j (indx+1)+ fill 0 0++-- prime counting++nthPrimeCt :: Integer -> Integer+nthPrimeCt 1 = 2+nthPrimeCt 2 = 3+nthPrimeCt 3 = 5+nthPrimeCt 4 = 7+nthPrimeCt 5 = 11+nthPrimeCt 6 = 13+nthPrimeCt n+ | n < 1 = error "nthPrimeCt: negative argument"+ | n < 200000 = let bd0 = nthPrimeApprox n+ bnd = bd0 + bd0 `quot` 32 + 37+ !sv = primeSieve bnd+ in countToNth (n-3) [sv]+ | otherwise = countToNth (n-3) (psieveList)++-- find the n-th set bit in a list of PrimeSieves,+-- aka find the (n+3)-rd prime+countToNth :: Integer -> [PrimeSieve] -> Integer+countToNth !n ps = runST (countDown n ps)++countDown :: Integer -> [PrimeSieve] -> ST s Integer+countDown !n (ps@(PS v0 bs) : more)+ | n > 278734 || (v0 /= 0 && n > 253000) = do+ ct <- countAll ps+ countDown (n - fromIntegral ct) more+ | otherwise = do+ stu <- unsafeThaw bs+ wa <- (castSTUArray :: STUArray s Int Bool -> ST s (STUArray s Int Word)) stu+ let go !k i+ | i == sieveWords = countDown k more+ | otherwise = do+ w <- unsafeRead wa i+ let !bc = fromIntegral $ bitCountWord w+ if bc < k+ then go (k-bc) (i+1)+ else let !j = fromIntegral (bc - k)+ !px = top w j (fromIntegral bc)+ in return (v0 + toPrim (px+(i `shiftL` WSHFT)))+ go n 0+countDown _ [] = error "Prime stream ended prematurely"++-- count all set bits in a chunk, do it wordwise for speed.+countAll :: PrimeSieve -> ST s Int+countAll (PS _ bs) = do+ stu <- unsafeThaw bs+ wa <- (castSTUArray :: STUArray s Int Bool -> ST s (STUArray s Int Word)) stu+ let go !ct i+ | i == sieveWords = return ct+ | otherwise = do+ w <- unsafeRead wa i+ go (ct + bitCountWord w) (i+1)+ go 0 0++-- Find the j-th highest of bc set bits in the Word w.+top :: Word -> Int -> Int -> Int+top w j bc = go 0 TOPB TOPM bn w+ where+ !bn = bc-j+ go !_ _ !_ !_ 0 = error "Too few bits set"+ go bs 0 _ _ wd = if wd .&. 1 == 0 then error "Too few bits, shift 0" else bs+ go bs a msk ix wd =+ case bitCountWord (wd .&. msk) of+ lc | lc < ix -> go (bs+a) a msk (ix-lc) (wd `uncheckedShiftR` a)+ | otherwise ->+ let !na = a `shiftR` 1+ in go bs na (msk `uncheckedShiftR` na) ix wd
+ Math/NumberTheory/Primes/Sieve/Indexing.hs view
@@ -0,0 +1,140 @@+-- |+-- Module: Math.NumberTheory.Primes.Sieve.Indexing+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+{-# OPTIONS_HADDOCK hide #-}+module Math.NumberTheory.Primes.Sieve.Indexing+ ( idxPr+ , toPrim+ , toIdx+ , rho+ , delta+ , tau+ , byte+ , idx+ , mu+ , nu+ ) where++import Data.Array.Unboxed+import Data.Array.Base (unsafeAt)+import Data.Bits++-- Auxiliary stuff, conversion between number and index,+-- remainders modulo 30 and related things.++-- {-# SPECIALISE idxPr :: Integer -> (Int,Int),+-- Int -> (Int,Int),+-- Word -> (Int,Int)+-- #-}+{-# INLINE idxPr #-}+idxPr :: Integral a => a -> (Int,Int)+idxPr n0 = (fromIntegral bytes0, rm3)+ where+ n = if (fromIntegral n0 .&. 1 == (1 :: Int))+ then n0 else (n0-1)+ (bytes0,rm0) = (n-7) `quotRem` 30+ rm1 = fromIntegral rm0+ rm2 = rm1 `quot` 3+ rm3 = min 7 (if rm2 > 5 then rm2-1 else rm2)++-- {-# SPECIALISE toPrim :: Int -> Integer,+-- Int -> Int,+-- Int -> Word,+-- Int -> Word16+-- #-}+{-# INLINE toPrim #-}+toPrim :: Integral a => Int -> a+toPrim ix = 30*fromIntegral k + fromIntegral (rho i)+ where+ i = ix .&. 7+ k = ix `shiftR` 3++-- Assumes n >= 7, gcd n 30 == 1+{-# INLINE toIdx #-}+toIdx :: Integral a => a -> Int+toIdx n = 8*fromIntegral q+r2+ where+ (q,r) = (n-7) `quotRem` 30+ r1 = fromIntegral r `quot` 3+ r2 = min 7 (if r1 > 5 then r1-1 else r1)++{-# INLINE rho #-}+rho :: Int -> Int+rho i = unsafeAt residues i++residues :: UArray Int Int+residues = listArray (0,7) [7,11,13,17,19,23,29,31]++{-# INLINE delta #-}+delta :: Int -> Int+delta i = unsafeAt deltas i++deltas :: UArray Int Int+deltas = listArray (0,7) [4,2,4,2,4,6,2,6]++{-# INLINE tau #-}+tau :: Int -> Int+tau i = unsafeAt taus i++taus :: UArray Int Int+taus = listArray (0,63)+ [ 7, 4, 7, 4, 7, 12, 3, 12+ , 12, 6, 11, 6, 12, 18, 5, 18+ , 14, 7, 13, 7, 14, 21, 7, 21+ , 18, 9, 19, 9, 18, 27, 9, 27+ , 20, 10, 21, 10, 20, 30, 11, 30+ , 25, 12, 25, 12, 25, 36, 13, 36+ , 31, 15, 31, 15, 31, 47, 15, 47+ , 33, 17, 33, 17, 33, 49, 17, 49+ ]++{-# INLINE byte #-}+byte :: Int -> Int+byte i = unsafeAt startByte i++startByte :: UArray Int Int+startByte = listArray (0,7) [1,3,5,9,11,17,27,31]++{-# INLINE idx #-}+idx :: Int -> Int+idx i = unsafeAt startIdx i++startIdx :: UArray Int Int+startIdx = listArray (0,7) [4,7,4,4,7,4,7,7]++{-# INLINE mu #-}+mu :: Int -> Int+mu i = unsafeAt mArr i++{-# INLINE nu #-}+nu :: Int -> Int+nu i = unsafeAt nArr i++mArr :: UArray Int Int+mArr = listArray (0,63)+ [ 1, 2, 2, 3, 4, 5, 6, 7+ , 2, 3, 4, 6, 6, 8, 10, 11+ , 2, 4, 5, 7, 8, 9, 12, 13+ , 3, 6, 7, 9, 10, 12, 16, 17+ , 4, 6, 8, 10, 11, 14, 18, 19+ , 5, 8, 9, 12, 14, 17, 22, 23+ , 6, 10, 12, 16, 18, 22, 27, 29+ , 7, 11, 13, 17, 19, 23, 29, 31+ ]++nArr :: UArray Int Int+nArr = listArray (0,63)+ [ 4, 3, 7, 6, 2, 1, 5, 0+ , 3, 7, 5, 0, 6, 2, 4, 1+ , 7, 5, 4, 1, 0, 6, 3, 2+ , 6, 0, 1, 4, 5, 7, 2, 3+ , 2, 6, 0, 5, 7, 3, 1, 4+ , 1, 2, 6, 7, 3, 4, 0, 5+ , 5, 4, 3, 2, 1, 0, 7, 6+ , 0, 1, 2, 3, 4, 5, 6, 7+ ]
+ Math/NumberTheory/Primes/Sieve/Misc.hs view
@@ -0,0 +1,344 @@+-- |+-- Module: Math.NumberTheory.Primes.Sieve.Misc+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, MonoLocalBinds, FlexibleContexts #-}+#if __GLASGOW_HASKELL__ >= 700+{-# OPTIONS_GHC -fspec-constr-count=8 #-}+#endif+{-# OPTIONS_HADDOCK hide #-}+module Math.NumberTheory.Primes.Sieve.Misc+ ( -- * Types+ FactorSieve+ , TotientSieve+ , CarmichaelSieve+ -- * Functions+ -- ** Smallest prime factors+ , factorSieve+ , sieveFactor+ , fsBound+ , fsPrimeTest+ -- ** Totients+ , totientSieve+ , sieveTotient+ -- ** Carmichael+ , carmichaelSieve+ , sieveCarmichael+ ) where++import Control.Monad.ST+import Data.Array.Base (unsafeRead, unsafeWrite, unsafeAt)+import Data.Array.ST+import Data.Array.Unboxed+import Control.Monad (when)+import Data.Bits+import GHC.Word++import System.Random++import Math.NumberTheory.Powers.Squares (integerSquareRoot')+import Math.NumberTheory.Primes.Sieve.Indexing+import Math.NumberTheory.Primes.Factorisation.Montgomery+import Math.NumberTheory.Primes.Factorisation.Utils+import Math.NumberTheory.Utils++-- | A compact store of smallest prime factors.+data FactorSieve = FS {-# UNPACK #-} !Word {-# UNPACK #-} !(UArray Int Word16)++-- | A compact store of totients.+data TotientSieve = TS {-# UNPACK #-} !Word {-# UNPACK #-} !(UArray Int Word)++-- | A compact store of values of the Carmichael function.+data CarmichaelSieve = CS {-# UNPACK #-} !Word {-# UNPACK #-} !(UArray Int Word)++-- | @'factorSieve' n@ creates a store of smallest prime factors of the numbers not exceeding @n@.+-- If you need to factorise many smallish numbers, this can give a big speedup since it avoids+-- many superfluous divisions. However, a too large sieve leads to a slowdown due to cache misses.+-- To reduce space usage, only the smallest prime factors of numbers coprime to @30@ are stored,+-- encoded as 'Word16's. The maximal admissible value for @n@ is therefore @2^32 - 1@.+-- Since @φ(30) = 8@, the sieve uses only @16@ bytes per @30@ numbers.+factorSieve :: Integer -> FactorSieve+factorSieve bound+ | 4294967295 < bound = error "factorSieve: overflow"+ | bound < 8 = FS 7 (array (0,0) [(0,0)])+ | otherwise = FS bnd (runSTUArray (spfSieve bnd))+ where+ bnd = fromInteger bound++fsBound :: FactorSieve -> Word+fsBound (FS b _) = b++fsPrimeTest :: FactorSieve -> Integer -> Bool+fsPrimeTest fs@(FS bnd sve) n+ | n < 0 = fsPrimeTest fs (-n)+ | n < 2 = False+ | fromInteger n .&. (1 :: Int) == 0 = n == 2+ | n `rem` 3 == 0 = n == 3+ | n `rem` 5 == 0 = n == 5+ | n <= fromIntegral bnd = sve `unsafeAt` (toIdx n) == 0+ | otherwise = error "Out of bounds"++-- | @'sieveFactor' fs n@ finds the prime factorisation of @n@ using the 'FactorSieve' @fs@.+-- For negative @n@, a factor of @-1@ is included with multiplicity @1@.+-- After stripping any present factors @2, 3@ or @5@, the remaining cofactor @c@ (if larger+-- than @1@) is factorised with @fs@. This is most efficient of course if @c@ does not+-- exceed the bound with which @fs@ was constructed. If it does, trial division is performed+-- until either the cofactor falls below the bound or the sieve is exhausted. In the latter+-- case, the elliptic curve method is used to finish the factorisation.+sieveFactor :: FactorSieve -> Integer -> [(Integer,Int)]+sieveFactor (FS bnd sve) = check+ where+ bound = fromIntegral bnd+ check 0 = error "0 has no prime factorisation"+ check 1 = []+ check n+ | n < 0 = (-1,1) : check (-n)+ | otherwise = go2 n+ go2 n = case shiftToOddCount n of+ (0,_) -> go3 n+ (k,m) -> (2,k) : if m == 1 then [] else go3 m+ go3 n = case splitOff 3 n of+ (0,_) -> go5 n+ (k,m) -> (3,k) : if m == 1 then [] else go5 m+ go5 n = case splitOff 5 n of+ (0,_) -> if n < 49 then [(n,1)] else sieveLoop n+ (k,m) -> (5,k) : case m of+ 1 -> []+ j | j < 49 -> [(j,1)]+ | otherwise -> sieveLoop j+ sieveLoop n+ | bound < n = tdLoop n (integerSquareRoot' n) 0+ | otherwise = intLoop (fromIntegral n)+ intLoop :: Word -> [(Integer,Int)]+ intLoop n = case unsafeAt sve (toIdx n) of+ 0 -> [(fromIntegral n,1)]+ p -> case splitOff (fromIntegral p) n of+ (k,m) -> (fromIntegral p, k) : case m of+ 1 -> []+ _ -> intLoop m+ lstIdx = snd (bounds sve)+ tdLoop n sr ix+ | lstIdx < ix = curve n+ | sr < p = [(n,1)]+ | pix /= 0 = tdLoop n sr (ix+1) -- not a prime+ | otherwise = case splitOff p n of+ (0,_) -> tdLoop n sr (ix+1)+ (k,m) -> (p,k) : case m of+ 1 -> []+ j | j <= bound -> intLoop (fromIntegral j)+ | otherwise -> tdLoop j (integerSquareRoot' j) (ix+1)+ where+ p = toPrim ix+ pix = unsafeAt sve ix+ curve n = stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n++-- | @'totientSieve' n@ creates a store of the totients of the numbers not exceeding @n@.+-- Like a 'FactorSieve', a 'TotientSieve' only stores values for numbers coprime to @30@+-- to reduce space usage. However, totients are stored as 'Word's, thus the space usage is+-- @2@ or @4@ times as high. The maximal admissible value for @n@ is @'fromIntegral' ('maxBound' :: 'Word')@.+totientSieve :: Integer -> TotientSieve+totientSieve bound+ | fromIntegral (maxBound :: Word) < bound = error "totientSieve: overflow"+ | bound < 8 = TS 7 (array (0,0) [(0,6)])+ | otherwise = TS bnd (totSieve bnd)+ where+ bnd = fromInteger bound++-- | @'sieveTotient' ts n@ finds the totient @π(n)@, i.e. the number of integers @k@ with+-- @1 <= k <= n@ and @'gcd' n k == 1@, in other words, the order of the group of units+-- in @ℤ/(n)@, using the 'TotientSieve' @ts@.+-- The strategy is analogous to 'sieveFactor'.+sieveTotient :: TotientSieve -> Integer -> Integer+sieveTotient (TS bnd sve) = check+ where+ bound = fromIntegral bnd+ check n+ | n < 1 = error "Totient only defined for positive numbers"+ | n == 1 = 1+ | otherwise = go2 n+ go2 n = case shiftToOddCount n of+ (0,_) -> go3 1 n+ (k,m) -> let tt = (shiftL 1 (k-1)) in if m == 1 then tt else go3 tt m+ go3 !tt n = case splitOff 3 n of+ (0,_) -> go5 tt n+ (k,m) -> let nt = tt*(2*3^(k-1)) in if m == 1 then nt else go5 nt m+ go5 !tt n = case splitOff 5 n of+ (0,_) -> sieveLoop tt n+ (k,m) -> let nt = tt*(4*5^(k-1)) in if m == 1 then nt else sieveLoop nt m+ sieveLoop !tt n+ | bound < n = tdLoop tt n (integerSquareRoot' n) 0+ | otherwise = case unsafeAt sve (toIdx n) of+ nt -> tt*fromIntegral nt+ lstIdx = snd (bounds sve)+ tdLoop !tt n sr ix+ | lstIdx < ix = curve tt n+ | sr < p' = tt*(n-1) -- n is a prime+ | pix /= p-1 = tdLoop tt n sr (ix+1) -- not a prime, next+ | otherwise = case splitOff p' n of+ (0,_) -> tdLoop tt n sr (ix+1)+ (k,m) -> let nt = tt*ppTotient (p',k)+ in case m of+ 1 -> nt+ j | j <= bound -> nt*fromIntegral (unsafeAt sve (toIdx j))+ | otherwise -> tdLoop nt j (integerSquareRoot' j) (ix+1)+ where+ p = toPrim ix+ p' = fromIntegral p+ pix = unsafeAt sve ix+ curve tt n = tt * totientFromCanonical (stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n)++-- | @'carmichaelSieve' n@ creates a store of values of the Carmichael function+-- for numbers not exceeding @n@.+-- Like a 'FactorSieve', a 'CarmichaelSieve' only stores values for numbers coprime to @30@+-- to reduce space usage. However, values are stored as 'Word's, thus the space usage is+-- @2@ or @4@ times as high. The maximal admissible value for @n@ is @'fromIntegral' ('maxBound' :: 'Word')@.+carmichaelSieve :: Integer -> CarmichaelSieve+carmichaelSieve bound+ | fromIntegral (maxBound :: Word) < bound = error "carmichaelSieve: overflow"+ | bound < 8 = CS 7 (array (0,0) [(0,6)])+ | otherwise = CS bnd (carSieve bnd)+ where+ bnd = fromInteger bound++-- | @'sieveCarmichael' cs n@ finds the value of @λ(n)@ (or @ψ(n)@), the smallest positive+-- integer @e@ such that for all @a@ with @gcd a n == 1@ the congruence @a^e ≡ 1 (mod n)@ holds,+-- in other words, the (smallest) exponent of the group of units in @ℤ/(n)@.+-- The strategy is analogous to 'sieveFactor'.+sieveCarmichael :: CarmichaelSieve -> Integer -> Integer+sieveCarmichael (CS bnd sve) = check+ where+ bound = fromIntegral bnd+ check n+ | n < 1 = error "Carmichael function only defined for positive numbers"+ | n == 1 = 1+ | otherwise = go2 n+ go2 n = case shiftToOddCount n of+ (0,_) -> go3 1 n+ (k,m) -> let tt = case k of+ 1 -> 1+ 2 -> 2+ _ -> (shiftL 1 (k-2))+ in if m == 1 then tt else go3 tt m+ go3 !tt n = case splitOff 3 n of+ (0,_) -> go5 tt n+ (k,1) | tt == 1 -> 2*3^(k-1)+ | otherwise -> tt*3^(k-1)+ (k,m) | tt == 1 -> go5 (2*3^(k-1)) m+ | otherwise -> go5 (tt*3^(k-1)) m+ go5 !tt n = case splitOff 5 n of+ (0,_) -> sieveLoop tt n+ (k,m) -> let tt' = case fromInteger tt .&. (3 :: Int) of+ 0 -> tt+ 2 -> 2*tt+ _ -> 4*tt+ nt = tt'*5^(k-1)+ in if m == 1 then nt else sieveLoop nt m+ sieveLoop !tt n+ | bound < n = tdLoop tt n (integerSquareRoot' n) 0+ | otherwise = case unsafeAt sve (toIdx n) of+ nt -> tt `lcm` fromIntegral nt+ lstIdx = snd (bounds sve)+ tdLoop !tt n sr ix+ | lstIdx < ix = curve tt n+ | sr < p' = tt `lcm` (n-1) -- n is a prime+ | pix /= p-1 = tdLoop tt n sr (ix+1) -- not a prime, next+ | otherwise = case splitOff p' n of+ (0,_) -> tdLoop tt n sr (ix+1)+ (k,m) -> let nt = (lcm tt (p'-1))*p'^(k-1)+ in case m of+ 1 -> nt+ j | j <= bound -> nt `lcm` fromIntegral (unsafeAt sve (toIdx j))+ | otherwise -> tdLoop nt j (integerSquareRoot' j) (ix+1)+ where+ p = toPrim ix+ p' = fromIntegral p+ pix = unsafeAt sve ix+ curve tt n = tt `lcm` carmichaelFromCanonical (stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n)+++{-# SPECIALISE spfSieve :: Word -> ST s (STUArray s Int Word),+ Word -> ST s (STUArray s Int Word16)+ #-}+spfSieve :: forall s w. (Integral w, MArray (STUArray s) w (ST s)) => Word -> ST s (STUArray s Int w)+spfSieve bound = do+ let (octs,lidx) = idxPr bound+ !mxidx = 8*octs+lidx+ mxval :: Word+ mxval = 30*fromIntegral octs + fromIntegral (rho lidx)+ !mxsve = integerSquareRoot' mxval+ (kr,r) = idxPr mxsve+ !svbd = 8*kr+r+ ar <- newArray (0,mxidx) 0 :: ST s (STUArray s Int w)+ let start k i = 8*(k*(30*k+2*rho i) + byte i) + idx i+ tick p stp off j ix+ | mxidx < ix = return ()+ | otherwise = do+ s <- (unsafeRead :: STUArray s Int w -> Int -> ST s w) ar ix+ when (s == 0) ((unsafeWrite :: STUArray s Int w -> Int -> w -> ST s ()) ar ix p)+ tick p stp off ((j+1) .&. 7) (ix + stp*delta j + tau (off+j))+ sift ix+ | svbd < ix = return ar+ | otherwise = do+ e <- unsafeRead ar ix+ when (e == 0) (do let i = ix .&. 7+ k = ix `shiftR` 3+ !off = i `shiftL` 3+ !stp = ix - i+ !p = toPrim ix+ tick p stp off i (start k i))+ sift (ix+1)+ sift 0++totSieve :: Word -> UArray Int Word+totSieve bound = runSTUArray $ do+ ar <- spfSieve bound+ (_,lst) <- getBounds ar+ let tot ix+ | lst < ix = return ar+ | otherwise = do+ p <- unsafeRead ar ix+ if p == 0+ then unsafeWrite ar ix (toPrim ix - 1)+ else do let !n = toPrim ix+ (tp,m) = unFact p (n `quot` p)+ case m of+ 1 -> unsafeWrite ar ix tp+ _ -> do+ tm <- unsafeRead ar (toIdx m)+ unsafeWrite ar ix (tp*tm)+ tot (ix+1)+ tot 0++carSieve :: Word -> UArray Int Word+carSieve bound = runSTUArray $ do+ ar <- spfSieve bound+ (_,lst) <- getBounds ar+ let car ix+ | lst < ix = return ar+ | otherwise = do+ p <- unsafeRead ar ix+ if p == 0+ then unsafeWrite ar ix (toPrim ix - 1)+ else do let !n = toPrim ix+ (tp,m) = unFact p (n `quot` p)+ case m of+ 1 -> unsafeWrite ar ix tp+ _ -> do+ tm <- unsafeRead ar (toIdx m)+ unsafeWrite ar ix (lcm tp tm)+ car (ix+1)+ car 0++-- Find the p-part of the totient of (p*m) and the cofactor+-- of the p-power in m.+{-# INLINE unFact #-}+unFact :: Word -> Word -> (Word,Word)+unFact p m = go (p-1) m+ where+ go !tt k = case k `quotRem` p of+ (q,0) -> go (p*tt) q+ _ -> (tt,k)
+ Math/NumberTheory/Primes/Testing.hs view
@@ -0,0 +1,40 @@+-- |+-- Module: Math.NumberTheory.Primes.Testing+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Primality tests.+module Math.NumberTheory.Primes.Testing+ ( -- * Standard tests+ isPrime+ -- $certificates+ , bailliePSW+ , millerRabinV+ , isStrongFermatPP+ , isFermatPP+ -- * Using a sieve+ , FactorSieve+ , fsIsPrime+ ) where++import Math.NumberTheory.Primes.Testing.Probabilistic+import Math.NumberTheory.Primes.Sieve.Misc++-- | Test primality using a 'FactorSieve'. If @n@ is out of bounds+-- of the sieve, fall back to 'isPrime'.+fsIsPrime :: FactorSieve -> Integer -> Bool+fsIsPrime fs n+ | n < 0 = fsIsPrime fs (-n)+ | n <= fromIntegral (fsBound fs) = fsPrimeTest fs n+ | otherwise = isPrime n++-- $certificates+--+-- The tests in this module may wrongly consider some composite numbers as prime.+-- For the Baillie-PSW test, no pseudoprimes are known, and it is known that none+-- exist below @2^64@, so for most practical purposes it can be regarded as conclusive.+-- Nevertheless, it is desirable to certify numbers passing it as primes (or find that+-- they are composite). The addition of prime certificates is planned for the next release.
+ Math/NumberTheory/Primes/Testing/Probabilistic.hs view
@@ -0,0 +1,224 @@+-- |+-- Module: Math.NumberTheory.Primes.Testing.Probabilistic+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Probabilistic primality tests, Miller-Rabin and Baillie-PSW.+{-# LANGUAGE CPP, MagicHash, BangPatterns #-}+{-# OPTIONS_HADDOCK hide #-}+module Math.NumberTheory.Primes.Testing.Probabilistic+ ( isPrime+ , millerRabinV+ , bailliePSW+ , isStrongFermatPP+ , isFermatPP+ , lucasTest+ ) where++#include "MachDeps.h"++import Math.NumberTheory.Moduli+import Math.NumberTheory.Utils+import Math.NumberTheory.Powers.Squares++import Data.Bits++import GHC.Base+import GHC.Word+import GHC.Integer.GMP.Internals++-- | @'isPrime' n@ tests whether @n@ is a prime (negative or positive).+-- First, trial division by the primes less than @1200@ is performed.+-- If that hasn't determined primality or compositeness, a Baillie PSW+-- test is performed.+--+-- Since the Baillie PSW test may not be perfect, it is possible that+-- some large composites are wrongly deemed prime, however, no composites+-- passing the test are known and none exist below @2^64@.+isPrime :: Integer -> Bool+isPrime n+ | n < 0 = isPrime (-n)+ | n < 2 = False+ | n < 4 = True+ | otherwise = go smallPrimes+ where+ go (p:ps)+ | p*p > n = True+ | otherwise = case n `rem` p of+ 0 -> False+ _ -> go ps+ go [] = bailliePSW n++-- | A Miller-Rabin like probabilistic primality test with preceding+-- trial division. While the classic Miller-Rabin test uses+-- randomly chosen bases, @'millerRabinV' k n@ uses the @k@+-- smallest primes as bases if trial division has not reached+-- a conclusive result. (Only the primes up to @1200@ are+-- available in this module, so the maximal effective @k@ is @196@.)+millerRabinV :: Int -> Integer -> Bool+millerRabinV k n+ | n < 0 = millerRabinV k (-n)+ | n < 2 = False+ | n < 4 = True+ | otherwise = go smallPrimes+ where+ go (p:ps)+ | p*p > n = True+ | otherwise = (n `rem` p /= 0) && go ps+ go [] = all (isStrongFermatPP n) (take k smallPrimes)++-- | @'isStrongFermatPP' n b@ tests whether @n@ is a strong Fermat+-- probable prime for base @b@, where @n > 2@ and @1 < b < n@.+-- The conditions on the arguments are not checked.+--+-- Apart from primes, also some composite numbers have the tested+-- property, but those are rare. Very rare are composite numbers+-- having the property for many bases, so testing a large prime+-- candidate with several bases can identify composite numbers+-- with high probability. An odd number @n > 3@ is prime if and+-- only if @'isStrongFermatPP' n b@ holds for all @b@ with+-- @2 <= b <= (n-1)/2@, but of course checking all those bases+-- would be less efficient than trial division, so one normally+-- checks only a relatively small number of bases, depending on+-- the desired degree of certainty. The probability that a randomly+-- chosen base doesn't identify a composite number @n@ is less than+-- @1/4@, so five to ten tests give a reasonable level of certainty+-- in general.+--+-- Some notes about the choice of bases: @b@ is a strong Fermat base+-- for @n@ if and only if @n-b@ is, hence one needs only test @b <= (n-1)/2@.+-- If @b@ is a strong Fermat base for @n@, then so is @b^k `mod` n@ for+-- all @k > 1@, hence one needs not test perfect powers, since their+-- base yields a stronger condition. Finally, if @a@ and @b@ are strong+-- Fermat bases for @n@, then @a*b@ is in most cases a strong Fermat+-- base for @n@, it can only fail to be so if @n `mod` 4 == 1@ and+-- the strong Fermat condition is reached at the same step for @a@ as for @b@,+-- so primes are the most powerful bases to test.+isStrongFermatPP :: Integer -> Integer -> Bool+isStrongFermatPP n b = a == 1 || go t a+ where+ m = n-1+ (t,u) = shiftToOddCount m+ a = powerModInteger' b u n+ go 0 _ = False+ go k x = x == m || go (k-1) ((x*x) `rem` n)++-- | @'isFermatPP' n b@ tests whether @n@ is a Fermat probable prime+-- for the base @b@, that is, whether @b^(n-1) `mod` n == 1@.+-- This is a weaker but simpler condition. However, more is lost+-- in strength than is gained in simplicity, so for primality testing,+-- the strong check should be used. The remarks about+-- the choice of bases to test from @'isStrongFermatPP'@ apply+-- with the modification that if @a@ and @b@ are Fermat bases+-- for @n@, then @a*b@ /always/ is a Fermat base for @n@ too.+-- A /Charmichael number/ is a composite number @n@ which is a+-- Fermat probable prime for all bases @b@ coprime to @n@. By the+-- above, only primes @p <= n/2@ not dividing @n@ need to be tested+-- to identify Carmichael numbers (however, testing all those+-- primes would be less efficient than determining Carmichaelness+-- from the prime factorisation; but testing an appropriate number+-- of prime bases is reasonable to find out whether it's worth the+-- effort to undertake the prime factorisation).+isFermatPP :: Integer -> Integer -> Bool+isFermatPP n b = powerModInteger' b (n-1) n == 1++-- | Primality test after Baillie, Pomerance, Selfridge and Wagstaff.+-- The Baillie PSW test consists of a strong Fermat probable primality+-- test followed by a (strong) Lucas primality test. This implementation+-- assumes that the number @n@ to test is odd and larger than @3@.+-- Even and small numbers have to be handled before. Also, before+-- applying this test, trial division by small primes should be performed+-- to identify many composites cheaply (although the Baillie PSW test is+-- rather fast, about the same speed as a strong Fermat test for four or+-- five bases usually, it is, for large numbers, much more costly than+-- trial division by small primes, the primes less than @1000@, say, so+-- eliminating numbers with small prime factors beforehand is more efficient).+--+-- The Baillie PSW test is very reliable, so far no composite numbers+-- passing it are known, and it is known (Gilchrist 2010) that no+-- Baillie PSW pseudoprimes exist below @2^64@. However, a heuristic argument+-- by Pomerance indicates that there are likely infinitely many Baillie PSW+-- pseudoprimes. On the other hand, according to+-- <http://mathworld.wolfram.com/Baillie-PSWPrimalityTest.html> there is+-- reason to believe that there are none with less than several+-- thousand digits, so that for most use cases the test can be+-- considered definitive.+bailliePSW :: Integer -> Bool+bailliePSW n = isStrongFermatPP n 2 && lucasTest n++-- precondition: n odd, > 3 (no small prime factors, typically large)+-- | The Lucas-Selfridge test, including square-check, but without+-- the Fermat test. For package-internal use only.+lucasTest :: Integer -> Bool+lucasTest n+ | square || d == 0 = False+ | d == 1 = True+ | otherwise = uo == 0 || go t vo qo+ where+ square = isPossibleSquare2 n && r*r == n+ r = integerSquareRoot n+ d = find True 5+ find !pos cd = case jacobi' (n `rem` cd) cd of+ 0 -> if cd == n then 1 else 0+ 1 -> find (not pos) (cd+2)+ _ -> if pos then cd else (-cd)+ q = (1-d) `quot` 4+ (t,o) = shiftToOddCount (n+1)+ (uo, vo, qo) = testLucas n q o+ go 0 _ _ = False+ go s vn qn = vn == 0 || go (s-1) ((vn*vn-2*qn) `rem` n) ((qn*qn) `rem` n)+++-- n odd positive, n > abs q, index odd+testLucas :: Integer -> Integer -> Integer -> (Integer, Integer, Integer)+testLucas n q (S# i#) = look (WORD_SIZE_IN_BITS - 2)+ where+ j = I# i#+ look k+ | testBit j k = go (k-1) 1 1 1 q+ | otherwise = look (k-1)+ go k un un1 vn qn+ | k < 0 = (un, vn, qn)+ | testBit j k = go (k-1) u2n1 u2n2 v2n1 q2n1+ | otherwise = go (k-1) u2n u2n1 v2n q2n+ where+ u2n = (un*vn) `rem` n+ u2n1 = (un1*vn-qn) `rem` n+ u2n2 = ((un1-q*un)*vn-qn) `rem` n+ v2n = (vn*vn-2*qn) `rem` n+ v2n1 = ((un1 - (2*q)*un)*vn-qn) `rem` n+ q2n = (qn*qn) `rem` n+ q2n1 = (qn*qn*q) `rem` n+testLucas n q (J# s# ba#) = test (s# -# 1#)+ where+ test j# = case indexWordArray# ba# j# of+ 0## -> test (j# -# 1#)+ w# -> look (j# -# 1#) (W# w#) (WORD_SIZE_IN_BITS - 1)+ look j# w i+ | 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#+ 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+ | otherwise = go k# w (i-1) u2n u2n1 v2n q2n+ where+ u2n = (un*vn) `rem` n+ u2n1 = (un1*vn-qn) `rem` n+ u2n2 = ((un1-q*un)*vn-qn) `rem` n+ v2n = (vn*vn-2*qn) `rem` n+ v2n1 = ((un1 - (2*q)*un)*vn-qn) `rem` n+ q2n = (qn*qn) `rem` n+ q2n1 = (qn*qn*q) `rem` n++smallPrimes :: [Integer]+smallPrimes = 2:3:5:prs+ where+ prs = 7:11:filter isPr (takeWhile (< 1200) . scanl (+) 13 $ cycle [4,2,4,6,2,6,4,2])+ isPr n = td n prs+ td n (p:ps) = (p*p > n) || (n `rem` p /= 0 && td n ps)+ td _ [] = True
+ Math/NumberTheory/Utils.hs view
@@ -0,0 +1,187 @@+-- |+-- Module: Math.NumberTheory.Utils+-- Copyright: (c) 2011 Daniel Fischer+-- Licence: MIT+-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability: Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Some utilities for bit twiddling.+--+{-# LANGUAGE CPP, MagicHash, UnboxedTuples, BangPatterns #-}+{-# OPTIONS_HADDOCK hide #-}+module Math.NumberTheory.Utils+ ( shiftToOddCount+ , shiftToOdd+ , shiftToOdd#+ , shiftToOddCount#+ , bitCountWord+ , bitCountInt+ , bitCountWord#+ , uncheckedShiftR+ , splitOff+ ) where++#include "MachDeps.h"++import GHC.Base+import GHC.Word++import GHC.Integer+import GHC.Integer.GMP.Internals++import Data.Bits++#if WORD_SIZE_IN_BITS == 64+#define m5 0x5555555555555555+#define m3 0x3333333333333333+#define mf 0x0F0F0F0F0F0F0F0F+#define m1 0x0101010101010101+#define sd 56+#else+#define m5 0x55555555+#define m3 0x33333333+#define mf 0x0F0F0F0F+#define m1 0x01010101+#define sd 24+#endif++uncheckedShiftR :: Word -> Int -> Word+uncheckedShiftR (W# w#) (I# i#) = W# (uncheckedShiftRL# w# i#)++-- | Remove factors of @2@ and count them. If+-- @n = 2^k*m@ with @m@ odd, the result is @(k, m)@.+-- Precondition: argument not @0@ (not checked).+{-# RULES+"shiftToOddCount/Int" shiftToOddCount = shiftOCInt+"shiftToOddCount/Word" shiftToOddCount = shiftOCWord+"shiftToOddCount/Integer" shiftToOddCount = shiftOCInteger+ #-}+shiftToOddCount :: (Integral a, Bits a) => a -> (Int, a)+shiftToOddCount n = case shiftOCInteger (fromIntegral n) of+ (z, o) -> (z, fromInteger o)++-- | Specialised version for @'Word'@.+-- Precondition: argument strictly positive (not checked).+shiftOCWord :: Word -> (Int, Word)+shiftOCWord (W# w#) = case shiftToOddCount# w# of+ (# z# , u# #) -> (I# z#, W# u#)++-- | Specialised version for @'Int'@.+-- Precondition: argument nonzero (not checked).+shiftOCInt :: Int -> (Int, Int)+shiftOCInt (I# i#) = case shiftToOddCount# (int2Word# i#) of+ (# z#, u# #) -> (I# z#, I# (word2Int# u#))++-- | Specialised version for @'Integer'@.+-- Precondition: argument nonzero (not checked).+shiftOCInteger :: Integer -> (Int, Integer)+shiftOCInteger n@(S# i#) =+ case shiftToOddCount# (int2Word# i#) of+ (# z#, w# #)+ | z# ==# 0# -> (0, n)+ | otherwise -> (I# z#, S# (word2Int# w#))+shiftOCInteger n@(J# _ ba#) = case count 0# 0# of+ 0# -> (0, n)+ z# -> (I# z#, n `shiftRInteger` z#)+ where+ count a# i# =+ case indexWordArray# ba# i# of+ 0## -> count (a# +# WORD_SIZE_IN_BITS#) (i# +# 1#)+ w# -> a# +# trailZeros# w#+++-- | Remove factors of @2@. If @n = 2^k*m@ with @m@ odd, the result is @m@.+-- Precondition: argument not @0@ (not checked).+{-# RULES+"shiftToOdd/Int" shiftToOdd = shiftOInt+"shiftToOdd/Word" shiftToOdd = shiftOWord+"shiftToOdd/Integer" shiftToOdd = shiftOInteger+ #-}+shiftToOdd :: (Integral a, Bits a) => a -> a+shiftToOdd n = fromInteger (shiftOInteger (fromIntegral n))++-- | Specialised version for @'Int'@.+-- Precondition: argument nonzero (not checked).+shiftOInt :: Int -> Int+shiftOInt (I# i#) = I# (word2Int# (shiftToOdd# (int2Word# i#)))++-- | Specialised version for @'Word'@.+-- Precondition: argument nonzero (not checked).+shiftOWord :: Word -> Word+shiftOWord (W# w#) = W# (shiftToOdd# w#)++-- | Specialised version for @'Int'@.+-- Precondition: argument nonzero (not checked).+shiftOInteger :: Integer -> Integer+shiftOInteger (S# i#) = S# (word2Int# (shiftToOdd# (int2Word# i#)))+shiftOInteger n@(J# _ ba#) = case count 0# 0# of+ 0# -> n+ z# -> n `shiftRInteger` z#+ where+ count a# i# =+ case indexWordArray# ba# i# of+ 0## -> count (a# +# WORD_SIZE_IN_BITS#) (i# +# 1#)+ w# -> a# +# trailZeros# w#++-- | Shift argument right until the result is odd.+-- Precondition: argument not @0@, not checked.+shiftToOdd# :: Word# -> Word#+shiftToOdd# w# = case trailZeros# w# of+ k# -> uncheckedShiftRL# w# k#++-- | Like @'shiftToOdd#'@, but count the number of places to shift too.+shiftToOddCount# :: Word# -> (# Int#, Word# #)+shiftToOddCount# w# = case trailZeros# w# of+ k# -> (# k#, uncheckedShiftRL# w# k# #)++-- | Number of 1-bits in a @'Word#'@.+bitCountWord# :: Word# -> Int#+bitCountWord# w# = case bitCountWord (W# w#) of+ I# i# -> i#++-- | Number of 1-bits in a @'Word'@.+bitCountWord :: Word -> Int+#if __GLASGOW_HASKELL__ >= 703+bitCountWord = popCount+-- should yield a machine instruction+#else+bitCountWord w = case w - (shiftR w 1 .&. m5) of+ !w1 -> case (w1 .&. m3) + (shiftR w1 2 .&. m3) of+ !w2 -> case (w2 + shiftR w2 4) .&. mf of+ !w3 -> fromIntegral (shiftR (w3 * m1) sd)+#endif++-- | Number of 1-bits in an @'Int'@.+bitCountInt :: Int -> Int+#if __GLASGOW_HASKELL__ >= 703+bitCountWord = popCount+-- should yield a machine instruction+#else+bitCountInt (I# i#) = bitCountWord (W# (int2Word# i#))+#endif++-- | Number of trailing zeros in a @'Word#'@, wrong for @0@.+{-# INLINE trailZeros# #-}+trailZeros# :: Word# -> Int#+trailZeros# w =+ case xor# w (w `minusWord#` 1##) `uncheckedShiftRL#` 1# of+ v0 ->+ case v0 `minusWord#` (uncheckedShiftRL# v0 1# `and#` m5##) of+ v1 ->+ case (v1 `and#` m3##) `plusWord#` (uncheckedShiftRL# v1 2# `and#` m3##) of+ v2 ->+ case (v2 `plusWord#` uncheckedShiftRL# v2 4#) `and#` mf## of+ v3 -> word2Int# (uncheckedShiftRL# (v3 `timesWord#` m1##) sd#)++-- {-# SPECIALISE splitOff :: Integer -> Integer -> (Int, Integer),+-- Int -> Int -> (Int, Int),+-- Word -> Word -> (Int, Word)+-- #-}+{-# INLINABLE splitOff #-}+splitOff :: Integral a => a -> a -> (Int, a)+splitOff p n = go 0 n+ where+ go !k m = case m `quotRem` p of+ (q,r) | r == 0 -> go (k+1) q+ | otherwise -> (k,m)
+ Setup.hs view
@@ -0,0 +1,5 @@+module Main where++import Distribution.Simple++main = defaultMain
+ TODO view
@@ -0,0 +1,5 @@+- Prime Certificates+- Atkin sieve+- General number field sieve+- Portability+- Check whether bit twiddling can be as fast as the lookup table for leading and trailing zeros
+ arithmoi.cabal view
@@ -0,0 +1,71 @@+name : arithmoi+version : 0.1.0.0+cabal-version : >= 1.6+author : Daniel Fischer+copyright : (c) 2011 Daniel Fischer+license : MIT+license-file : LICENCE+maintainer : Daniel Fischer <daniel.is.fischer@googlemail.com>+build-type : Simple+stability : Provisional+homepage : https://bitbucket.org/dafis/arithmoi+bug-reports : https://bitbucket.org/dafis/arithmoi/issues++synopsis : Efficient basic number-theoretic functions.+ Primes, powers, integer logarithms.+description : A library of basic functionality needed for+ number-theoretic calculations. The aim of this library+ is to provide efficient implementations of the functions.++ Primes and related things (totients, factorisation),+ powers (integer roots and tests, modular exponentiation),+ integer logarithms.++ Note: Requires GHC >= 6.12 with the integer-gmp package+ for efficiency. Portability is on the to-do list (with+ low priority, however).++category : Math, Algorithms, Number Theory++tested-with : GHC == 6.12.3, GHC == 7.0.2, GHC == 7.0.3, GHC == 7.2.1++extra-source-files : Changes, TODO++library+ build-depends : base >= 4 && < 5, array >= 0.3 && < 0.4, ghc-prim,+ integer-gmp, containers >= 0.3 && < 0.5, random >= 1.0 && < 1.1,+ mtl >= 2.0 && < 2.1+ exposed-modules : Math.NumberTheory.Logarithms+ Math.NumberTheory.Moduli+ Math.NumberTheory.Lucas+ Math.NumberTheory.GCD+ Math.NumberTheory.GCD.LowLevel+ Math.NumberTheory.Powers+ Math.NumberTheory.Powers.Squares+ Math.NumberTheory.Powers.Cubes+ Math.NumberTheory.Powers.Fourth+ Math.NumberTheory.Powers.General+ Math.NumberTheory.Primes+ Math.NumberTheory.Primes.Sieve+ Math.NumberTheory.Primes.Factorisation+ Math.NumberTheory.Primes.Counting+ Math.NumberTheory.Primes.Testing+ Math.NumberTheory.Primes.Heap+ other-modules : Math.NumberTheory.Utils+ Math.NumberTheory.Logarithms.Internal+ Math.NumberTheory.Powers.Integer+ Math.NumberTheory.Primes.Counting.Impl+ Math.NumberTheory.Primes.Counting.Approximate+ Math.NumberTheory.Primes.Factorisation.Montgomery+ Math.NumberTheory.Primes.Factorisation.Utils+ Math.NumberTheory.Primes.Sieve.Eratosthenes+ Math.NumberTheory.Primes.Sieve.Indexing+ Math.NumberTheory.Primes.Sieve.Misc+ Math.NumberTheory.Primes.Testing.Probabilistic+ extensions : BangPatterns+ ghc-options : -O2 -Wall+ ghc-prof-options : -O2 -auto++source-repository head+ type: mercurial+ location: https://bitbucket.org/dafis/arithmoi