Copyright | (c) 2011 Daniel Fischer |
---|---|

License | MIT |

Maintainer | Daniel Fischer <daniel.is.fischer@googlemail.com> |

Stability | Provisional |

Portability | Non-portable (GHC extensions) |

Safe Haskell | None |

Language | Haskell2010 |

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.

- factorise :: Integer -> [(Integer, Int)]
- defaultStdGenFactorisation :: StdGen -> Integer -> [(Integer, Int)]
- stepFactorisation :: Integer -> [(Integer, Int)]
- factorise' :: Integer -> [(Integer, Int)]
- defaultStdGenFactorisation' :: StdGen -> Integer -> [(Integer, Int)]
- data FactorSieve
- factorSieve :: Integer -> FactorSieve
- sieveFactor :: FactorSieve -> Integer -> [(Integer, Int)]
- trialDivisionTo :: Integer -> Integer -> [(Integer, Int)]
- smallFactors :: Integer -> Integer -> ([(Integer, Int)], Maybe Integer)
- stdGenFactorisation :: Maybe Integer -> StdGen -> Maybe Int -> Integer -> [(Integer, Int)]
- curveFactorisation :: Maybe Integer -> (Integer -> Bool) -> (Integer -> g -> (Integer, g)) -> g -> Maybe Int -> Integer -> [(Integer, Int)]
- montgomeryFactorisation :: Integer -> Word -> Word -> Integer -> Maybe Integer
- data TotientSieve
- totientSieve :: Integer -> TotientSieve
- sieveTotient :: TotientSieve -> Integer -> Integer
- data CarmichaelSieve
- carmichaelSieve :: Integer -> CarmichaelSieve
- sieveCarmichael :: CarmichaelSieve -> Integer -> Integer

# Factorisation functions

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.

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

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 `defaultStdGenFactorisation`

`-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)] Source #

is like `stepFactorisation`

`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)] Source #

Like `factorise`

, but without input checking, hence `n > 1`

is required.

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

Like `defaultStdGenFactorisation`

, but without input checking, so
`n`

must be larger than `1`

.

### Factor sieves

data FactorSieve Source #

A compact store of smallest prime factors.

factorSieve :: Integer -> FactorSieve Source #

creates a store of smallest prime factors of the numbers not exceeding `factorSieve`

n`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.
The prime factors are stored as `Word16`

for compactness, so `n`

must be
smaller than `2^32`

.

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

finds the prime factorisation of `sieveFactor`

fs n`n`

using the `FactorSieve`

`fs`

.
For negative `n`

, a factor of `-1`

is included with multiplicity `1`

.
After stripping any present factors `2`

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

### Trial division

trialDivisionTo :: Integer -> Integer -> [(Integer, Int)] Source #

produces a factorisation of `trialDivisionTo`

bound n`n`

using the
primes `<= bound`

. If `n`

has prime divisors `> bound`

, the last entry
in the list is the product of all these. If `n <= bound^2`

, this is a
full factorisation, but very slow if `n`

has large prime divisors.

## Partial factorisation

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

finds all prime divisors of `smallFactors`

bound n`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.

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

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

is the driver for the factorisation. Its performance (and success)
can be influenced by passing appropriate arguments. If you know that `curveFactorisation`

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

requires that small prime factors have been stripped before. Also, it is
unlikely to succeed if `curveFactorisation`

`n`

has more than one (really) large prime factor.

### Single curve worker

montgomeryFactorisation :: Integer -> Word -> Word -> Integer -> Maybe Integer Source #

tries to find a factor of `montgomeryFactorisation`

n b1 b2 s`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

data TotientSieve Source #

A compact store of totients.

totientSieve :: Integer -> TotientSieve Source #

creates a store of the totients of the numbers not exceeding `totientSieve`

n`n`

.
A `TotientSieve`

only stores values for numbers coprime to `30`

to reduce space usage.
The maximal admissible value for `n`

is

.`fromIntegral`

(`maxBound`

:: `Word`

)

sieveTotient :: TotientSieve -> Integer -> Integer Source #

finds the totient `sieveTotient`

ts n`π(n)`

, i.e. the number of integers `k`

with
`1 <= k <= n`

and

, in other words, the order of the group of units
in `gcd`

n k == 1`ℤ/(n)`

, using the `TotientSieve`

`ts`

.
First, factors of `2, 3`

and `5`

are handled individually, if the remaining
cofactor of `n`

is within the sieve range, its totient is looked up, otherwise
the calculation falls back on factorisation, first by trial division and
if necessary, elliptic curves.

# Carmichael function

data CarmichaelSieve Source #

A compact store of values of the Carmichael function.

carmichaelSieve :: Integer -> CarmichaelSieve Source #

creates a store of values of the Carmichael function
for numbers not exceeding `carmichaelSieve`

n`n`

.
Like a `TotientSieve`

, a `CarmichaelSieve`

only stores values for numbers coprime to `30`

to reduce space usage. The maximal admissible value for `n`

is

.`fromIntegral`

(`maxBound`

:: `Word`

)

sieveCarmichael :: CarmichaelSieve -> Integer -> Integer Source #

finds the value of `sieveCarmichael`

cs n`λ(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 `sieveTotient`

.