-- Hoogle documentation, generated by Haddock
-- See Hoogle, http://www.haskell.org/hoogle/
-- | Efficient basic number-theoretic functions.
Primes, powers, integer logarithms.
--
-- 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).
@package arithmoi
@version 0.2.0.2
-- | 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.
module Math.NumberTheory.Primes.Heap
-- | 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.
primes :: Integral a => [a]
-- | sieveFrom n generates the list of primes >=
-- n. The remarks about overflow and performance in the
-- documentation of primes apply here too.
sieveFrom :: Integral a => a -> [a]
-- | Efficient calculation of Lucas sequences.
module Math.NumberTheory.Lucas
-- | 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
-- | 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)
-- | lucas k computes the k-th Lucas number. Very
-- similar to fibonacci.
lucas :: Int -> Integer
-- | 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)
-- | 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)
-- | 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.
module Math.NumberTheory.GCD.LowLevel
-- | Greatest common divisor of two Ints, calculated with the binary
-- gcd algorithm.
gcdInt :: Int -> Int -> Int
-- | Greatest common divisor of two Words, calculated with the
-- binary gcd algorithm.
gcdWord :: Word -> Word -> Word
-- | Greatest common divisor of two Int#s, calculated with the
-- binary gcd algorithm.
gcdInt# :: Int# -> Int# -> Int#
-- | Greatest common divisor of two Word#s, calculated with the
-- binary gcd algorithm.
gcdWord# :: Word# -> Word# -> Word#
-- | Test whether two Ints are coprime, using an abbreviated binary
-- gcd algorithm.
coprimeInt :: Int -> Int -> Bool
-- | Test whether two Words are coprime, using an abbreviated binary
-- gcd algorithm.
coprimeWord :: Word -> Word -> Bool
-- | Test whether two Int#s are coprime.
coprimeInt# :: Int# -> Int# -> Bool
-- | Test whether two Word#s are coprime.
coprimeWord# :: Word# -> Word# -> Bool
-- | 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.
module Math.NumberTheory.GCD
-- | Calculate the greatest common divisor using the binary gcd algorithm.
-- Depending on type and hardware, that can be considerably faster than
-- 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 <= 32, 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
-- | 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)
-- | 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
-- | Miscellaneous functions related to modular arithmetic.
module Math.NumberTheory.Moduli
-- | 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.
jacobi :: (Integral a, Bits a) => a -> a -> Int
-- | 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
-- | 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.
powerMod :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer
-- | Specialised version of powerMod for Integer exponents.
-- Reduces the number of shifts of the exponent since shifting large
-- Integers is expensive. Call this function directly if you don't
-- want or can't rely on rewrite rules.
powerModInteger :: Integer -> Integer -> Integer -> Integer
-- | Jacobi symbol of two numbers without validity check of the
-- "denominator".
jacobi' :: (Integral a, Bits a) => a -> a -> Int
-- | Modular power worker without input checking. Assumes all arguments
-- strictly positive and modulus greater than 1.
powerMod' :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer
-- | Specialised worker without input checks. Makes the same assumptions as
-- the general version powerMod'.
powerModInteger' :: Integer -> Integer -> Integer -> Integer
-- | Functions dealing with squares. Efficient calculation of integer
-- square roots and efficient testing for squareness.
module Math.NumberTheory.Powers.Squares
-- | 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.
integerSquareRoot :: Integral a => a -> a
-- | 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.
integerSquareRoot' :: Integral a => a -> a
-- | 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.
exactSquareRoot :: Integral a => a -> Maybe a
-- | 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.
isSquare :: Integral a => a -> Bool
-- | 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.
isSquare' :: Integral a => a -> Bool
-- | 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.
isPossibleSquare :: Integral a => a -> Bool
-- | 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).
isPossibleSquare2 :: Integral a => a -> Bool
-- | 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
-- | List of primes. Since the sieve uses unboxed arrays, overflow occurs
-- at some point. On 64-bit systems, that point is beyond the memory
-- limits, on 32-bit systems, it is at about 1.7*10^18.
primes :: [Integer]
-- | sieveFrom n creates the list of primes not less than
-- n.
sieveFrom :: Integer -> [Integer]
-- | Compact store of primality flags.
data PrimeSieve
-- | 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
-- | List of primes in the form of a list of PrimeSieves, 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]
-- | psieveFrom n creates the list of PrimeSieves
-- 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]
-- | Generate a list of primes for consumption from a PrimeSieve.
primeList :: PrimeSieve -> [Integer]
-- | Functions dealing with cubes. Moderately efficient calculation of
-- integer cube roots and testing for cubeness.
module Math.NumberTheory.Powers.Cubes
-- | 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.
integerCubeRoot :: Integral a => a -> a
-- | 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.
integerCubeRoot' :: Integral a => a -> a
-- | Returns Nothing if the argument is not a cube, Just
-- r if n == r^3.
exactCubeRoot :: Integral a => a -> Maybe a
-- | Test whether an integer is a cube.
isCube :: Integral a => a -> Bool
-- | 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.
isCube' :: Integral a => a -> Bool
-- | 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.
isPossibleCube :: Integral a => a -> Bool
-- | Functions dealing with fourth powers. Efficient calculation of integer
-- fourth roots and efficient testing for being a square's square.
module Math.NumberTheory.Powers.Fourth
-- | 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
-- | Calculate the integer fourth root of a nonnegative number, that is,
-- the largest integer r with r^4 <= n. The
-- condition is not checked.
integerFourthRoot' :: Integral a => a -> a
-- | Returns Nothing if n is not a fourth power, Just
-- r if n == r^4 and r >= 0.
exactFourthRoot :: Integral a => a -> Maybe a
-- | Test whether an integer is a fourth power. First nonnegativity is
-- checked, then the unchecked test is called.
isFourthPower :: Integral a => a -> Bool
-- | 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.
isFourthPower' :: Integral a => a -> Bool
-- | Test whether a nonnegative number is a possible fourth power. The
-- condition is not checked. This eliminates about 99.958% of
-- numbers.
isPossibleFourthPower :: Integral a => a -> Bool
-- | Calculating integer roots and determining perfect powers. The
-- algorithms are moderately efficient.
module Math.NumberTheory.Powers.General
-- | 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.
integerRoot :: (Integral a, Integral b) => b -> a -> a
-- | 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
-- | isKthPower k n checks whether n is a
-- k-th power.
isKthPower :: (Integral a, Integral b) => b -> a -> Bool
-- | isPerfectPower n checks whether n == r^k for
-- some k > 1.
isPerfectPower :: Integral a => a -> Bool
-- | 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)
-- | 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
-- | 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.
integerSquareRoot :: Integral a => a -> a
-- | 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.
isSquare :: Integral a => a -> Bool
-- | 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.
exactSquareRoot :: Integral a => a -> Maybe a
-- | 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.
integerCubeRoot :: Integral a => a -> a
-- | Test whether an integer is a cube.
isCube :: Integral a => a -> Bool
-- | Returns Nothing if the argument is not a cube, Just
-- r if n == r^3.
exactCubeRoot :: Integral a => a -> Maybe a
-- | 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
-- | Test whether an integer is a fourth power. First nonnegativity is
-- checked, then the unchecked test is called.
isFourthPower :: Integral a => a -> Bool
-- | Returns Nothing if n is not a fourth power, Just
-- r if n == r^4 and r >= 0.
exactFourthRoot :: Integral a => a -> Maybe a
-- | 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.
integerRoot :: (Integral a, Integral b) => b -> a -> a
-- | isKthPower k n checks whether n is a
-- k-th power.
isKthPower :: (Integral a, Integral b) => b -> a -> Bool
-- | 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
-- | isPerfectPower n checks whether n == r^k for
-- some k > 1.
isPerfectPower :: Integral a => a -> Bool
-- | 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)
-- | 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.
powerMod :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer
-- | Integer Logarithms. For efficiency, the internal representation of
-- Integers from integer-gmp is used.
module Math.NumberTheory.Logarithms
-- | 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
-- | Calculate the integer logarithm of an Integer to base 2. The
-- argument must be positive, otherwise an error is thrown.
integerLog2 :: Integer -> Int
-- | Calculate the integer logarithm of an Int to base 2. The
-- argument must be positive, otherwise an error is thrown.
intLog2 :: Int -> Int
-- | Calculate the integer logarithm of a Word to base 2. The
-- argument must be positive, otherwise an error is thrown.
wordLog2 :: Word -> Int
-- | Same as integerLogBase, but without checks, saves a little time
-- when called often for known good input.
integerLogBase' :: Integer -> Integer -> Int
-- | Same as integerLog2, but without checks, saves a little time
-- when called often for known good input.
integerLog2' :: Integer -> Int
-- | Same as intLog2, but without checks, saves a little time when
-- called often for known good input.
intLog2' :: Int -> Int
-- | Same as wordLog2, but without checks, saves a little time when
-- called often for known good input.
wordLog2' :: Word -> Int
-- | 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
-- | 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)]
-- | 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)]
-- | 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)]
-- | Like factorise, but without input checking, hence n >
-- 1 is required.
factorise' :: Integer -> [(Integer, Int)]
-- | Like defaultStdGenFactorisation, but without input checking, so
-- n must be larger than 1.
defaultStdGenFactorisation' :: StdGen -> Integer -> [(Integer, Int)]
-- | A compact store of smallest prime factors.
data 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.
factorSieve :: Integer -> FactorSieve
-- | 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)]
-- | trialDivisionTo bound n produces a factorisation of
-- 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.
trialDivisionTo :: Integer -> Integer -> [(Integer, Int)]
-- | 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)
-- | 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 -> StdGen -> Maybe Int -> Integer -> [(Integer, Int)]
-- | 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 -> (Integer -> Bool) -> (Integer -> g -> (Integer, g)) -> g -> Maybe Int -> Integer -> [(Integer, Int)]
-- | 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
-- | 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
-- | Alias of totient for people who prefer Greek letters.
φ :: Integer -> Integer
-- | A compact store of totients.
data 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).
totientSieve :: Integer -> TotientSieve
-- | 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
-- | Calculate the totient from the canonical factorisation.
totientFromCanonical :: [(Integer, Int)] -> Integer
-- | Calculates the Carmichael function for a positive integer, that is,
-- the (smallest) exponent of the group of units in
-- &8484;/(n).
carmichael :: Integer -> Integer
-- | Alias of carmichael for people who prefer Greek letters.
λ :: Integer -> Integer
-- | A compact store of values of the Carmichael function.
data 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).
carmichaelSieve :: Integer -> CarmichaelSieve
-- | 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
-- | Calculate the Carmichael function from the factorisation. Requires
-- that the list of prime factors is strictly ascending.
carmichaelFromCanonical :: [(Integer, Int)] -> 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.
divisors :: Integer -> Set Integer
-- | tau n is the number of (positive) divisors of
-- n. tau 0 is an error because 0 has
-- infinitely many divisors.
tau :: Integer -> Integer
-- | Alias for tau for people preferring Greek letters.
τ :: Integer -> Integer
-- | Alias for tau.
divisorCount :: Integer -> Integer
-- | The sum of all (positive) divisors of a positive number n,
-- calculated from its prime factorisation.
divisorSum :: 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).
sigma :: Int -> Integer -> Integer
-- | Alias for sigma for people preferring Greek letters.
σ :: Int -> Integer -> Integer
-- | Alias for sigma.
divisorPowerSum :: Int -> Integer -> Integer
-- | The set of divisors, efficiently calculated from the canonical
-- factorisation.
divisorsFromCanonical :: [(Integer, Int)] -> Set Integer
-- | The number of divisors, efficiently calculated from the canonical
-- factorisation.
tauFromCanonical :: [(a, Int)] -> Integer
-- | The sum of all divisors, efficiently calculated from the canonical
-- factorisation.
divisorSumFromCanonical :: [(Integer, Int)] -> Integer
-- | The sum of the powers (with fixed exponent) of all divisors,
-- efficiently calculated from the canonical factorisation.
sigmaFromCanonical :: Int -> [(Integer, Int)] -> Integer
-- | Factorisation proving the primality of the found factors.
--
-- For large numbers, this will be very slow in general. Use only if
-- you're paranoid or must be really sure.
module Math.NumberTheory.Primes.Factorisation.Certified
-- | certifiedFactorisation n produces the prime
-- factorisation of n, proving the primality of the factors, but
-- doesn't report the proofs.
certifiedFactorisation :: Integer -> [(Integer, Int)]
-- | certificateFactorisation n produces a
-- provenFactorisation with a default bound of 100000.
certificateFactorisation :: Integer -> [((Integer, Int), PrimalityProof)]
-- | provenFactorisation bound n constructs a the prime
-- factorisation of n (which must be positive) together with
-- proofs of primality of the factors, using trial division up to
-- bound (which is arbitrarily replaced by 2000 if the
-- supplied value is smaller) and elliptic curve factorisation for the
-- remaining factors if necessary.
--
-- Construction of primality proofs can take a very long time, so
-- this will usually be slow (but should be faster than using
-- factorise and proving the primality of the factors from
-- scratch).
provenFactorisation :: Integer -> Integer -> [((Integer, Int), PrimalityProof)]
-- | Certificates for primality or compositeness.
module Math.NumberTheory.Primes.Testing.Certificates
-- | A certificate of either compositeness or primality of an
-- Integer. Only numbers > 1 can be certified, trying
-- to create a certificate for other numbers raises an error.
data Certificate
Composite :: !CompositenessProof -> Certificate
Prime :: !PrimalityProof -> Certificate
argueCertificate :: Certificate -> Either CompositenessArgument PrimalityArgument
-- | A proof of compositeness of a positive number. The type is abstract to
-- ensure the validity of proofs.
data CompositenessProof
-- | The number whose compositeness is proved.
composite :: CompositenessProof -> Integer
-- | A proof of primality of a positive number. The type is abstract to
-- ensure the validity of proofs.
data PrimalityProof
-- | The number whose primality is proved.
cprime :: PrimalityProof -> Integer
-- | An argument for compositeness of a number (which must be >
-- 1). CompositenessProofs translate directly to
-- CompositenessArguments, correct arguments can be transformed
-- into proofs. This type allows the manipulation of proofs while
-- maintaining their correctness. The only way to access components of a
-- CompositenessProof except the composite is through this type.
data CompositenessArgument
-- | compo == firstDiv*secondDiv, where all are > 1
Divisors :: Integer -> Integer -> Integer -> CompositenessArgument
compo :: CompositenessArgument -> Integer
firstDivisor :: CompositenessArgument -> Integer
secondDivisor :: CompositenessArgument -> Integer
-- | compo fails the strong Fermat test for fermatBase
Fermat :: Integer -> Integer -> CompositenessArgument
compo :: CompositenessArgument -> Integer
fermatBase :: CompositenessArgument -> Integer
-- | compo fails the Lucas-Selfridge test
Lucas :: Integer -> CompositenessArgument
compo :: CompositenessArgument -> Integer
-- | No particular reason given
Belief :: Integer -> CompositenessArgument
compo :: CompositenessArgument -> Integer
-- | An argument for primality of a number (which must be > 1).
-- PrimalityProofs translate directly to
-- PrimalityArguments, correct arguments can be transformed into
-- proofs. This type allows the manipulation of proofs while maintaining
-- their correctness. The only way to access components of a
-- PrimalityProof except the prime is through this type.
data PrimalityArgument
-- | A suggested Pocklington certificate
Pock :: Integer -> Integer -> Integer -> [(Integer, Int, Integer, PrimalityArgument)] -> PrimalityArgument
aprime :: PrimalityArgument -> Integer
largeFactor :: PrimalityArgument -> Integer
smallFactor :: PrimalityArgument -> Integer
factorList :: PrimalityArgument -> [(Integer, Int, Integer, PrimalityArgument)]
-- | Primality should be provable by trial division to alimit
Division :: Integer -> Integer -> PrimalityArgument
aprime :: PrimalityArgument -> Integer
alimit :: PrimalityArgument -> Integer
-- | aprime is said to be obviously prime, that holds for primes
-- < 30
Obvious :: Integer -> PrimalityArgument
aprime :: PrimalityArgument -> Integer
-- | Primality assumed
Assumption :: Integer -> PrimalityArgument
aprime :: PrimalityArgument -> Integer
-- | arguePrimality transforms a proof of primality into an
-- argument for primality.
arguePrimality :: PrimalityProof -> PrimalityArgument
-- | argueCompositeness transforms a proof of compositeness
-- into an argument for compositeness.
argueCompositeness :: CompositenessProof -> CompositenessArgument
-- | verifyPrimalityArgument checks the given argument and
-- constructs a proof from it, if it is valid. For the explicit
-- arguments, this is simple and resonably fast, for an
-- Assumption, the verification uses certify and hence may
-- take a long time.
verifyPrimalityArgument :: PrimalityArgument -> Maybe PrimalityProof
-- | verifyCompositenessArgument checks the given argument
-- and constructs a proof from it, if it is valid. For the explicit
-- arguments, this is simple and resonably fast, for a Belief, the
-- verification uses certify and hence may take a long time.
verifyCompositenessArgument :: CompositenessArgument -> Maybe CompositenessProof
-- | certify n constructs, for n > 1, a proof
-- of either primality or compositeness of n. This may take a
-- very long time if the number has no small(ish) prime divisors
certify :: Integer -> Certificate
-- | Check the validity of a Certificate. Since it should be
-- impossible to construct invalid certificates by the public interface,
-- this should never return False.
checkCertificate :: Certificate -> Bool
-- | Check the validity of a CompositenessProof. Since it should be
-- impossible to create invalid proofs by the public interface, this
-- should never return False.
checkCompositenessProof :: CompositenessProof -> Bool
-- | Check the validity of a PrimalityProof. Since it should be
-- impossible to create invalid proofs by the public interface, this
-- should never return False.
checkPrimalityProof :: PrimalityProof -> Bool
-- | Primality tests.
module Math.NumberTheory.Primes.Testing
-- | 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
-- | isCertifiedPrime n tests primality of n,
-- first trial division by small primes is performed, then a Baillie PSW
-- test and finally a prime certificate is constructed and verified,
-- provided no step before found n to be composite. Constructing
-- prime certificates can take a very long time, so use this with
-- care.
isCertifiedPrime :: Integer -> Bool
-- | 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
-- | 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
-- | 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
-- | 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
-- | A compact store of smallest prime factors.
data FactorSieve
-- | Test primality using a FactorSieve. If n is out of
-- bounds of the sieve, fall back to isPrime.
fsIsPrime :: FactorSieve -> Integer -> Bool
-- | trialDivisionPrimeTo bound n tests whether n
-- is coprime to all primes <= bound. If n <=
-- bound^2, this is a full prime test, but very slow if n
-- has no small prime divisors.
trialDivisionPrimeTo :: Integer -> Integer -> Bool
-- | Number of primes not exceeding n, π(n), and
-- n-th prime; also fast, but reasonably accurate approximations
-- to these.
module Math.NumberTheory.Primes.Counting
-- | 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
-- | 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
-- | 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
-- | 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
module Math.NumberTheory.Primes