Math.NumberTheory.Primes.Factorisation

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.

Factorisation functions

Factorisation of Integers 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.

factorise :: Integer -> [(Integer, Int)]

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.

defaultStdGenFactorisation :: StdGen -> Integer -> [(Integer, Int)]

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.

stepFactorisation :: Integer -> [(Integer, Int)]

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.

factorise' :: Integer -> [(Integer, Int)]

Like factorise, but without input checking, hence n > 1 is required.

defaultStdGenFactorisation' :: StdGen -> Integer -> [(Integer, Int)]

Like defaultStdGenFactorisation, but without input checking, so n must be larger than 1.

Factor sieves

data FactorSieve

A compact store of smallest prime factors.

factorSieve :: Integer -> FactorSieve

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 Word16s. The maximal admissible value for n is therefore 2^32 - 1. Since φ(30) = 8, the sieve uses only 16 bytes per 30 numbers.

sieveFactor :: FactorSieve -> Integer -> [(Integer, Int)]

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.

Partial factorisation

smallFactors :: Integer -> Integer -> ([(Integer, Int)], Maybe Integer)

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.

stdGenFactorisation

Arguments

:: 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

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.

curveFactorisation

Arguments

:: 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 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.

Single curve worker

montgomeryFactorisation :: Integer -> Word -> Word -> Integer -> Maybe Integer

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.

Totients

totient :: Integer -> Integer

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).

φ :: Integer -> Integer

Alias of totient for people who prefer Greek letters.

data TotientSieve

A compact store of totients.

totientSieve :: Integer -> TotientSieve

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 Words, thus the space usage is 2 or 4 times as high. The maximal admissible value for n is fromIntegral (maxBound :: Word).

sieveTotient :: TotientSieve -> Integer -> Integer

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.

totientFromCanonical :: [(Integer, Int)] -> Integer

Calculate the totient from the canonical factorisation.

Carmichael function

carmichael :: Integer -> Integer

Calculates the Carmichael function for a positive integer, that is, the (smallest) exponent of the group of units in &8484;/(n).

λ :: Integer -> Integer

Alias of carmichael for people who prefer Greek letters.

data CarmichaelSieve

A compact store of values of the Carmichael function.

carmichaelSieve :: Integer -> CarmichaelSieve

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 Words, thus the space usage is 2 or 4 times as high. The maximal admissible value for n is fromIntegral (maxBound :: Word).

sieveCarmichael :: CarmichaelSieve -> Integer -> Integer

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.

carmichaelFromCanonical :: [(Integer, Int)] -> Integer

Calculate the Carmichael function from the factorisation. Requires that the list of prime factors is strictly ascending.

Divisors

divisors :: Integer -> Set Integer

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 Integers.

tau :: Integer -> Integer

tau n is the number of (positive) divisors of n. tau 0 is an error because 0 has infinitely many divisors.

τ :: Integer -> Integer

Alias for tau for people preferring Greek letters.

divisorCount :: Integer -> Integer

Alias for tau.

divisorSum :: Integer -> Integer

The sum of all (positive) divisors of a positive number n, calculated from its prime factorisation.

sigma :: Int -> Integer -> Integer

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).

σ :: Int -> Integer -> Integer

Alias for sigma for people preferring Greek letters.

divisorPowerSum :: Int -> Integer -> Integer

Alias for sigma.

divisorsFromCanonical :: [(Integer, Int)] -> Set Integer

The set of divisors, efficiently calculated from the canonical factorisation.

tauFromCanonical :: [(a, Int)] -> Integer

The number of divisors, efficiently calculated from the canonical factorisation.

divisorSumFromCanonical :: [(Integer, Int)] -> Integer

The sum of all divisors, efficiently calculated from the canonical factorisation.

sigmaFromCanonical :: Int -> [(Integer, Int)] -> Integer

The sum of the powers (with fixed exponent) of all divisors, efficiently calculated from the canonical factorisation.