NAME

Math::Prime::Util - Utilities related to prime numbers, including fast sieves and factoring

VERSION

Version 0.35

SYNOPSIS

# Normally you would just import the functions you are using.
# Nothing is exported by default.  List the functions, or use :all.
use Math::Prime::Util ':all';


# Get a big array reference of many primes
my $aref = primes( 100_000_000 );

# All the primes between 5k and 10k inclusive
my $aref = primes( 5_000, 10_000 );

# If you want them in an array instead
my @primes = @{primes( 500 )};

# You can do something for every prime in a range.  Twin primes to 10k:
forprimes { say if is_prime($_+2) } 10000;

# For non-bigints, is_prime and is_prob_prime will always be 0 or 2.
# They return 0 (composite), 2 (prime), or 1 (probably prime)
say "$n is prime"  if is_prime($n);
say "$n is ", (qw(composite maybe_prime? prime))[is_prob_prime($n)];

# Strong pseudoprime test with multiple bases, using Miller-Rabin
say "$n is a prime or 2/7/61-psp" if is_strong_pseudoprime($n, 2, 7, 61);

# Standard and strong Lucas-Selfridge, and extra strong Lucas tests
say "$n is a prime or lpsp"   if is_lucas_pseudoprime($n);
say "$n is a prime or slpsp"  if is_strong_lucas_pseudoprime($n);
say "$n is a prime or eslpsp" if is_extra_strong_lucas_pseudoprime($n);

# step to the next prime (returns 0 if not using bigints and we'd overflow)
$n = next_prime($n);

# step back (returns 0 if given input less than 2)
$n = prev_prime($n);


# Return Pi(n) -- the number of primes E<lt>= n.
$primepi = prime_count( 1_000_000 );
$primepi = prime_count( 10**14, 10**14+1000 );  # also does ranges

# Quickly return an approximation to Pi(n)
my $approx_number_of_primes = prime_count_approx( 10**17 );

# Lower and upper bounds.  lower <= Pi(n) <= upper for all n
die unless prime_count_lower($n) <= prime_count($n);
die unless prime_count_upper($n) >= prime_count($n);


# Return p_n, the nth prime
say "The ten thousandth prime is ", nth_prime(10_000);

# Return a quick approximation to the nth prime
say "The one trillionth prime is ~ ", nth_prime_approx(10**12);

# Lower and upper bounds.   lower <= nth_prime(n) <= upper for all n
die unless nth_prime_lower($n) <= nth_prime($n);
die unless nth_prime_upper($n) >= nth_prime($n);


# Get the prime factors of a number
@prime_factors = factor( $n );

# Return ([p1,e1],[p2,e2], ...) for $n = p1^e1 * p2*e2 * ...
@pe = factor_exp( $n );

# Get all divisors other than 1 and n
@divisors = all_factors( $n );

# Euler phi (Euler's totient) on a large number
use bigint;  say euler_phi( 801294088771394680000412 );
say jordan_totient(5, 1234);  # Jordan's totient

# Moebius function used to calculate Mertens
$sum += moebius($_) for (1..200); say "Mertens(200) = $sum";
# Mertens function directly (more efficient for large values)
say mertens(10_000_000);
# Exponential of Mangoldt function
say "lamba(49) = ", log(exp_mangoldt(49));
# Some more number theoretical functions
say liouville(4292384);
say chebyshev_psi(234984);
say chebyshev_theta(92384234);
say partitions(1000);

# divisor sum
$sigma  = divisor_sum( $n );       # sum of divisors
$sigma0 = divisor_sum( $n, 0 );    # count of divisors
$sigmak = divisor_sum( $n, $k );
$sigmaf = divisor_sum( $n, sub { log($_[0]) } ); # arbitrary func

# primorial n#, primorial p(n)#, and lcm
say "The product of primes below 47 is ",     primorial(47);
say "The product of the first 47 primes is ", pn_primorial(47);
say "lcm(1..1000) is ", consecutive_integer_lcm(1000);

# Ei, li, and Riemann R functions
my $ei   = ExponentialIntegral($x);   # $x a real: $x != 0
my $li   = LogarithmicIntegral($x);   # $x a real: $x >= 0
my $R    = RiemannR($x)               # $x a real: $x > 0
my $Zeta = RiemannZeta($x)            # $x a real: $x >= 0


# Precalculate a sieve, possibly speeding up later work.
prime_precalc( 1_000_000_000 );

# Free any memory used by the module.
prime_memfree;

# Alternate way to free.  When this leaves scope, memory is freed.
my $mf = Math::Prime::Util::MemFree->new;


# Random primes
my $small_prime = random_prime(1000);      # random prime <= limit
my $rand_prime = random_prime(100, 10000); # random prime within a range
my $rand_prime = random_ndigit_prime(6);   # random 6-digit prime
my $rand_prime = random_nbit_prime(128);   # random 128-bit prime
my $rand_prime = random_strong_prime(256); # random 256-bit strong prime
my $rand_prime = random_maurer_prime(256); # random 256-bit provable prime

DESCRIPTION

A set of utilities related to prime numbers. These include multiple sieving methods, is_prime, prime_count, nth_prime, approximations and bounds for the prime_count and nth prime, next_prime and prev_prime, factoring utilities, and more.

The default sieving and factoring are intended to be (and currently are) the fastest on CPAN, including Math::Prime::XS, Math::Prime::FastSieve, Math::Factor::XS, Math::Prime::TiedArray, Math::Big::Factors, Math::Factoring, and Math::Primality (when the GMP module is available). For numbers in the 10-20 digit range, it is often orders of magnitude faster. Typically it is faster than Math::Pari for 64-bit operations.

All operations support both Perl UV's (32-bit or 64-bit) and bignums. If you want high performance with big numbers (larger than Perl's native 32-bit or 64-bit size), you should install Math::Prime::Util::GMP and Math::BigInt::GMP. This will be a recurring theme throughout this documentation -- while all bignum operations are supported in pure Perl, most methods will be much slower than the C+GMP alternative.

The module is thread-safe and allows concurrency between Perl threads while still sharing a prime cache. It is not itself multi-threaded. See the Limitations section if you are using Win32 and threads in your program.

Two scripts are also included and installed by default:

primes.pl displays primes between start and end values or expressions, with many options for filtering (e.g. twin, safe, circular, good, lucky, etc.). Use --help to see all the options.
factor.pl operates similar to the GNU factor program. It supports bigint and expression inputs.

BIGNUM SUPPORT

By default all functions support bignums. For performance, you should install and use Math::BigInt::GMP or Math::BigInt::Pari, and Math::Prime::Util::GMP.

Some of the functions, including:

factor
factor_exp
is_pseudoprime
is_strong_pseudoprime
nth_prime
moebius
mertens
euler_phi
chebyshev_theta
chebyshev_psi
is_prime
is_prob_prime
next_prime
prev_prime

work very fast (under 1 microsecond) on small inputs, but the wrappers for input validation and bigint support take more time than the function itself. Using the flag '-bigint', e.g.:

use Math::Prime::Util qw(-bigint);

will turn off bigint support for those functions. Those functions will then go directly to the XS versions, which will speed up very small inputs a lot. This is useful if you're using the functions in a loop, but since the difference is less than a millisecond, it's really not important in general. The last five functions have shortcuts by default so will only skip validation.

If you are using bigints, here are some performance suggestions:

Install Math::Prime::Util::GMP, as that will vastly increase the speed of many of the functions. This does require the GMP library be installed on your system, but this increasingly comes pre-installed or easily available using the OS vendor package installation tool.
Install and use Math::BigInt::GMP or Math::BigInt::Pari, then use use bigint try => 'GMP,Pari' in your script, or on the command line -Mbigint=lib,GMP. Large modular exponentiation is much faster using the GMP or Pari backends, as are the math and approximation functions when called with very large inputs.
Install Math::MPFR if you use the Ei, li, Zeta, or R functions. If that module can be loaded, these functions will run much faster on bignum inputs, and are able to provide higher accuracy.
I have run these functions on many versions of Perl, and my experience is that if you're using anything older than Perl 5.14, I would recommend you upgrade if you are using bignums a lot. There are some brittle behaviors on 5.12.4 and earlier with bignums.

PRIMALITY TESTING

This module provides three functions for general primality testing, as well as numerous specialized functions. The three main functions are: "is_prob_prime" and "is_prime" for general use, and "is_provable_prime" for proofs. For inputs below 2^64 the functions are identical and fast deterministic testing is performed. That is, the results will always be correct and should take at most a few microseconds for any input. This is hundreds to thousands of times faster than other CPAN modules. For inputs larger than 2^64, an extra-strong BPSW test is used. See the "PRIMALITY TESTING NOTES" section for more discussion.

FUNCTIONS

is_prime

print "$n is prime" if is_prime($n);

Returns 2 if the number is prime, 0 if not. For numbers larger than 2^64 it will return 0 for composite and 1 for probably prime, using an extra-strong BPSW test. If Math::Prime::Util::GMP is installed, some additional primality tests are also performed on large inputs, and a quick attempt is made to perform a primality proof, so it will return 2 for many other inputs.

Also see the "is_prob_prime" function, which will never do additional tests, and the "is_provable_prime" function which will construct a proof that the input is number prime and returns 2 for almost all primes (at the expense of speed).

For native precision numbers (anything smaller than 2^64, all three functions are identical and use a deterministic set of tests (selected Miller-Rabin bases or BPSW). For larger inputs both "is_prob_prime" and "is_prime" return probable prime results using the extra-strong Baillie-PSW test, which has had no counterexample found since it was published in 1980.

For cryptographic key generation, you may want even more testing for probable primes (NIST recommends some additional M-R tests). This can be done using additional random bases with "is_strong_pseudoprime", or a different test such as "is_frobenius_underwood_pseudoprime". Even better, make sure Math::Prime::Util::GMP is installed and use "is_provable_prime" which should be reasonably fast for sizes under 2048 bits. Another possibility is to use "random_maurer_prime" in Math::Prime::Util which constructs a random provable prime.

primes

Returns all the primes between the lower and upper limits (inclusive), with a lower limit of 2 if none is given.

An array reference is returned (with large lists this is much faster and uses less memory than returning an array directly).

my $aref1 = primes( 1_000_000 );
my $aref2 = primes( 1_000_000_000_000, 1_000_000_001_000 );

my @primes = @{ primes( 500 ) };

print "$_\n" for (@{primes( 20, 100 )});

Sieving will be done if required. The algorithm used will depend on the range and whether a sieve result already exists. Possibilities include trial division (for ranges with only one expected prime), a Sieve of Eratosthenes using wheel factorization, or a segmented sieve.

next_prime

$n = next_prime($n);

Returns the next prime greater than the input number. If the input is not a bigint, then 0 is returned if the next prime is larger than a native integer type (the last representable primes being 4,294,967,291 in 32-bit Perl and 18,446,744,073,709,551,557 in 64-bit).

prev_prime

$n = prev_prime($n);

Returns the prime preceding the input number (i.e. the largest prime that is strictly less than the input). 0 is returned if the input is 2 or lower.

forprimes

forprimes { say } 100,200;                  # print primes from 100 to 200

$sum=0;  forprimes { $sum += $_ } 100000;   # sum primes to 100k

forprimes { say if is_prime($_+2) } 10000;  # print twin primes to 10k

Given a block and either an end count or a start and end pair, calls the block for each prime in the range. Compared to getting a big array of primes and iterating through it, this is more memory efficient and perhaps more convenient. This will almost always be the fastest way to loop over a range of primes. Nesting and using in threads are allowed.

Math::BigInt objects may be used for the range.

For some uses an iterator ("prime_iterator", "prime_iterator_object") or a tied array (Math::Prime::Util::PrimeArray) may be more convenient. Objects can be passed to functions, and allow early loop exits without exceptions. Here is a clumsy "forprimes" exception example:

use bigint;
eval { forprimes { die "$_\n" if $_ % 123 == 1 } 2**100, 2**101 };
my $n = 0+$@;

prime_iterator

my $it = prime_iterator;
$sum += $it->() for 1..100000;

Returns a closure-style iterator. The start value defaults to the first prime (2) but an initial value may be given as an argument, which will result in the first value returned being the next prime greater than or equal to the argument. For example, this:

my $it = prime_iterator(200);  say $it->();  say $it->();

will return 211 followed by 223, as those are the next primes >= 200. On each call, the iterator returns the current value and increments to the next prime.

Other options include "forprimes" (more efficiency, less flexibility), Math::Prime::Util::PrimeIterator (an iterator with more functionality), or Math::Prime::Util::PrimeArray (a tied array).

prime_iterator_object

my $it = prime_iterator_object;
while ($it->value < 100) { say $it->value; $it->next; }
$sum += $it->iterate for 1..100000;

Returns a Math::Prime::Util::PrimeIterator object. A shortcut that loads the package if needed, calls new, and returns the object. See the documentation for that package for details. This object has more features than the simple one above (e.g. the iterator is bi-directional), and also handles iterating across bigints.

prime_count

my $primepi = prime_count( 1_000 );
my $pirange = prime_count( 1_000, 10_000 );

Returns the Prime Count function Pi(n), also called primepi in some math packages. When given two arguments, it returns the inclusive count of primes between the ranges. E.g. (13,17) returns 2, (14,17) and (13,16) return 1, (14,16) returns 0.

The current implementation decides based on the ranges whether to use a segmented sieve with a fast bit count, or the extended LMO algorithm. The former is preferred for small sizes as well as small ranges. The latter is much faster for large ranges.

The segmented sieve is very memory efficient and is quite fast even with large base values. Its complexity is approximately O(sqrt(a) + (b-a)), where the first term is typically negligible below ~ 10^11. Memory use is proportional only to sqrt(a), with total memory use under 1MB for any base under 10^14.

The extended LMO method has complexity approximately O(b^(2/3)) + O(a^(2/3)), and also uses low memory. A calculation of Pi(10^14) completes in a few seconds, Pi(10^15) in well under a minute, and Pi(10^16) in about one minute. In contrast, even parallel primesieve would take over a week on a similar machine to determine Pi(10^16).

Also see the function "prime_count_approx" which gives a very good approximation to the prime count, and "prime_count_lower" and "prime_count_upper" which give tight bounds to the actual prime count. These functions return quickly for any input, including bigints.

prime_count_upper

prime_count_lower

my $lower_limit = prime_count_lower($n);
my $upper_limit = prime_count_upper($n);
#   $lower_limit  <=  prime_count(n)  <=  $upper_limit

Returns an upper or lower bound on the number of primes below the input number. These are analytical routines, so will take a fixed amount of time and no memory. The actual prime_count will always be equal to or between these numbers.

A common place these would be used is sizing an array to hold the first $n primes. It may be desirable to use a bit more memory than is necessary, to avoid calling prime_count.

These routines use verified tight limits below a range at least 2^35, and use the Dusart (2010) bounds of

x/logx * (1 + 1/logx + 2.000/log^2x) <= Pi(x)

x/logx * (1 + 1/logx + 2.334/log^2x) >= Pi(x)

above that range. These bounds do not assume the Riemann Hypothesis. If the configuration option assume_rh has been set (it is off by default), then the Schoenfeld (1976) bounds are used for large values.

prime_count_approx

print "there are about ",
      prime_count_approx( 10 ** 18 ),
      " primes below one quintillion.\n";

Returns an approximation to the prime_count function, without having to generate any primes. For values under 10^36 this uses the Riemann R function, which is quite accurate: an error of less than 0.0005% is typical for input values over 2^32, and decreases as the input gets larger. If Math::MPFR is installed, the Riemann R function is used for all values, and will be very fast. If not, then values of 10^36 and larger will use the approximation li(x) - li(sqrt(x))/2. While not as accurate as the Riemann R function, it still should have error less than 0.00000000000000001%.

A slightly faster but much less accurate answer can be obtained by averaging the upper and lower bounds.

nth_prime

say "The ten thousandth prime is ", nth_prime(10_000);

Returns the prime that lies in index n in the array of prime numbers. Put another way, this returns the smallest p such that Pi(p) >= n.

For relatively small inputs (below 1 million or so), this does a sieve over a range containing the nth prime, then counts up to the number. This is fairly efficient in time and memory. For larger values, create a low-biased estimate using the inverse logarithmic integral, use a fast prime count, then sieve in the small difference.

While this method is thousands of times faster than generating primes, and doesn't involve big tables of precomputed values, it still can take a fair amount of time for large inputs. Calculating the 10^12th prime takes about 1 second, the 10^13th prime takes under 10 seconds, and the 10^14th prime (3475385758524527) takes under one minute. Think about whether a bound or approximation would be acceptable, as they can be computed analytically.

If the result is larger than a native integer size (32-bit or 64-bit), the result will take a very long time. A later version of Math::Prime::Util::GMP may include this functionality which would help for 32-bit machines.

nth_prime_upper

nth_prime_lower

my $lower_limit = nth_prime_lower($n);
my $upper_limit = nth_prime_upper($n);
#   $lower_limit  <=  nth_prime(n)  <=  $upper_limit

Returns an analytical upper or lower bound on the Nth prime. These are very fast as they do not need to sieve or search through primes or tables. An exact answer is returned for tiny values of n. The lower limit uses the Dusart 2010 bound for all n, while the upper bound uses one of the two Dusart 2010 bounds for n >= 178974, a Dusart 1999 bound for n >= 39017, and a simple bound of n * (logn + 0.6 * loglogn) for small n.

nth_prime_approx

say "The one trillionth prime is ~ ", nth_prime_approx(10**12);

Returns an approximation to the nth_prime function, without having to generate any primes. Uses the Cipolla 1902 approximation with two polynomials, plus a correction for small values to reduce the error.

is_pseudoprime

Takes a positive number n and a base a as input, and returns 1 if n is a probable prime to base a. This is the simple Fermat primality test. Removing primes, given base 2 this produces the sequence OEIS A001567.

is_strong_pseudoprime

my $maybe_prime = is_strong_pseudoprime($n, 2);
my $probably_prime = is_strong_pseudoprime($n, 2, 3, 5, 7, 11, 13, 17);

Takes a positive number as input and one or more bases. The bases must be greater than 1. Returns 1 if the input is a strong probable prime to all of the bases, and 0 if not.

If 0 is returned, then the number really is a composite. If 1 is returned, then it is either a prime or a strong pseudoprime to all the given bases. Given enough distinct bases, the chances become very, very strong that the number is actually prime.

This is usually used in combination with other tests to make either stronger tests (e.g. the strong BPSW test) or deterministic results for numbers less than some verified limit (e.g. it has long been known that no more than three selected bases are required to give correct primality test results for any 32-bit number). Given the small chances of passing multiple bases, there are some math packages that just use multiple MR tests for primality testing.

Even inputs other than 2 will always return 0 (composite). While the algorithm does run with even input, most sources define it only on odd input. Returning composite for all non-2 even input makes the function match most other implementations including Math::Primality's is_strong_pseudoprime function.

miller_rabin

An alias for is_strong_pseudoprime. This name is deprecated.

is_lucas_pseudoprime

Takes a positive number as input, and returns 1 if the input is a standard Lucas probable prime using the Selfridge method of choosing D, P, and Q (some sources call this a Lucas-Selfridge pseudoprime). Removing primes, this produces the sequence OEIS A217120.

is_strong_lucas_pseudoprime

Takes a positive number as input, and returns 1 if the input is a strong Lucas probable prime using the Selfridge method of choosing D, P, and Q (some sources call this a strong Lucas-Selfridge pseudoprime). This is one half of the BPSW primality test (the Miller-Rabin strong pseudoprime test with base 2 being the other half). Removing primes, this produces the sequence OEIS A217255.

is_extra_strong_lucas_pseudoprime

Takes a positive number as input, and returns 1 if the input passes the extra strong Lucas test (as defined in Grantham 2000). This test has more stringent conditions than the strong Lucas test, and produces about 60% fewer pseudoprimes. Performance is typically 20-30% faster than the strong Lucas test.

The parameters are selected using the Baillie-OEIS method method: increment P from 3 until jacobi(D,n) = -1. Removing primes, this produces the sequence OEIS A217719.

is_almost_extra_strong_lucas_pseudoprime

This is similar to the "is_extra_strong_lucas_pseudoprime" function, but does not calculate U, so is a little faster, but also weaker. With the current implementations, there is little reason to prefer this unless trying to reproduce specific results. The extra-strong implementation has been optimized to use similar features, removing most of the performance advantage.

An optional second argument (an integer between 1 and 256) indicates the increment amount for P parameter selection. The default value of 1 yields the parameter selection described in "is_extra_strong_lucas_pseudoprime", creating a pseudoprime sequence which is a superset of the latter's pseudoprime sequence OEIS A217719. A value of 2 yields the method used by Pari.

Because the U = 0 condition is ignored, this produces about 5% more pseudoprimes than the extra-strong Lucas test. However this is still only 66% of the number produced by the strong Lucas-Selfridge test. No BPSW counterexamples have been found with any of the Lucas tests described.

is_frobenius_underwood_pseudoprime

Takes a positive number as input, and returns 1 if the input passes the minimal lambda+2 test (see Underwood 2012 "Quadratic Compositeness Tests"), where (L+2)^(n-1) = 5 + 2x mod (n, L^2 - Lx + 1). The computational cost for this is between the cost of 2 and 3 strong pseudoprime tests. There are no known counterexamples, but this is not a well studied test.

miller_rabin_random

Takes a positive number (n) as input and a positive number (k) of bases to use. Performs k Miller-Rabin tests using uniform random bases between 2 and n-2.

This should not be used in place of "is_prob_prime", "is_prime", or "is_provable_prime". Those functions will be faster and provide better results than running k Miller-Rabin tests. This function can be used if one wants more assurances for non-proven primes, such as for cryptographic uses where the size is large enough that proven primes are not desired.

is_prob_prime

my $prob_prime = is_prob_prime($n);
# Returns 0 (composite), 2 (prime), or 1 (probably prime)

Takes a positive number as input and returns back either 0 (composite), 2 (definitely prime), or 1 (probably prime).

For 64-bit input (native or bignum), this uses either a deterministic set of Miller-Rabin tests (1, 2, or 3 tests) or a strong BPSW test consisting of a single base-2 strong probable prime test followed by a strong Lucas test. This has been verified with Jan Feitsma's 2-PSP database to produce no false results for 64-bit inputs. Hence the result will always be 0 (composite) or 2 (prime).

For inputs larger than 2^64, an extra-strong Baillie-PSW primality test is performed (also called BPSW or BSW). This is a probabilistic test, so only 0 (composite) and 1 (probably prime) are returned. There is a possibility that composites may be returned marked prime, but since the test was published in 1980, not a single BPSW pseudoprime has been found, so it is extremely likely to be prime. While we believe (Pomerance 1984) that an infinite number of counterexamples exist, there is a weak conjecture (Martin) that none exist under 10000 digits.

is_provable_prime

say "$n is definitely prime" if is_provable_prime($n) == 2;

Takes a positive number as input and returns back either 0 (composite), 2 (definitely prime), or 1 (probably prime). This gives it the same return values as "is_prime" and "is_prob_prime". Note that numbers below 2^64 are considered proven by the deterministic set of Miller-Rabin bases or the BPSW test. Both of these have been tested for all small (64-bit) composites and do not return false positives.

Using the Math::Prime::Util::GMP module is highly recommended for doing primality proofs, as it is much, much faster. The pure Perl code is just not fast for this type of operation, nor does it have the best algorithms. It should suffice for proofs of up to 40 digit primes, while the latest MPU::GMP works for primes of hundreds of digits (thousands with an optional larger polynomial set).

The pure Perl implementation uses theorem 5 of BLS75 (Brillhart, Lehmer, and Selfridge's 1975 paper), an improvement on the Pocklington-Lehmer test. This requires n-1 to be factored to (n/2)^(1/3)). This is often fast, but as n gets larger, it takes exponentially longer to find factors.

Math::Prime::Util::GMP implements both the BLS75 theorem 5 test as well as ECPP (elliptic curve primality proving). It will typically try a quick n-1 proof before using ECPP. Certificates are available with either method. This results in proofs of 200-digit primes in under 1 second on average, and many hundreds of digits are possible. This makes it significantly faster than Pari 2.1.7's is_prime(n,1) which is the default for Math::Pari.

prime_certificate

my $cert = prime_certificate($n);
say verify_prime($cert) ? "proven prime" : "not prime";

Given a positive integer n as input, returns a primality certificate as a multi-line string. If we could not prove n prime, an empty string is returned (n may or may not be composite). This may be examined or given to "verify_prime" for verification. The latter function contains the description of the format.

is_provable_prime_with_cert

Given a positive integer as input, returns a two element array containing the result of "is_provable_prime": 0 definitely composite 1 probably prime 2 definitely prime and a primality certificate like "prime_certificate". The certificate will be an empty string if the first element is not 2.

verify_prime

my $cert = prime_certificate($n);
say verify_prime($cert) ? "proven prime" : "not prime";

Given a primality certificate, returns either 0 (not verified) or 1 (verified). Most computations are done using pure Perl with Math::BigInt, so you probably want to install and use Math::BigInt::GMP, and ECPP certificates will be faster with Math::Prime::Util::GMP for its elliptic curve computations.

If the certificate is malformed, the routine will carp a warning in addition to returning 0. If the verbose option is set (see "prime_set_config") then if the validation fails, the reason for the failure is printed in addition to returning 0. If the verbose option is set to 2 or higher, then a message indicating success and the certificate type is also printed.

A certificate may have arbitrary text before the beginning (the primality routines from this module will not have any extra text, but this way verbose output from the prover can be safely stored in a certificate). The certificate begins with the line:

[MPU - Primality Certificate]

All lines in the certificate beginning with # are treated as comments and ignored, as are blank lines. A version number may follow, such as:

Version 1.0

For all inputs, base 10 is the default, but at any point this may be changed with a line like:

Base 16

where allowed bases are 10, 16, and 62. This module will only use base 10, so its routines will not output Base commands.

Next, we look for (using "100003" as an example):

Proof for:
N 100003

where the text Proof for: indicates we will read an N value. Skipping comments and blank lines, the next line should be "N " followed by the number.

After this, we read one or more blocks. Each block is a proof of the form:

If Q is prime, then N is prime.

Some of the blocks have more than one Q value associated with them, but most only have one. Each block has its own set of conditions which must be verified, and this can be done completely self-contained. That is, each block is independent of the other blocks and may be processed in any order. To be a complete proof, each block must successfully verify. The block types and their conditions are shown below.

Finally, when all blocks have been read and verified, we must ensure we can construct a proof tree from the set of blocks. The root of the tree is the initial N, and for each node (block), all Q values must either have a block using that value as its N or Q must be less than 2^64 and pass BPSW.

Some other certificate formats (e.g. Primo) use an ordered chain, where the first block must be for the initial N, a single Q is given which is the implied N for the next block, and so on. This simplifies validation implementation somewhat, and removes some redundant information from the certificate, but has no obvious way to add proof types such as Lucas or the various BLS75 theorems that use multiple factors. I decided that the most general solution was to have the certificate contain the set in any order, and let the verifier do the work of constructing the tree.

The blocks begin with the text "Type ..." where ... is the type. One or more values follow. The defined types are:

Small

Type Small
N 5791

N must be less than 2^64 and be prime (use BPSW or deterministic M-R).

BLS3

Type BLS3
N  2297612322987260054928384863
Q  16501461106821092981
A  5

A simple n-1 style proof using BLS75 theorem 3. This block verifies if: a Q is odd b Q > 2 c Q divides N-1 . Let M = (N-1)/Q d MQ+1 = N e M > 0 f 2Q+1 > sqrt(N) g A^((N-1)/2) mod N = N-1 h A^(M/2) mod N != N-1

Pocklington

Type Pocklington
N  2297612322987260054928384863
Q  16501461106821092981
A  5

A simple n-1 style proof using generalized Pocklington. This is more restrictive than BLS3 and much more than BLS5. This is Primo's type 1, and this module does not currently generate these blocks. This block verifies if: a Q divides N-1 . Let M = (N-1)/Q b M > 0 c M < Q d MQ+1 = N e A > 1 f A^(N-1) mod N = 1 g gcd(A^M - 1, N) = 1

BLS15

Type BLS15
N  8087094497428743437627091507362881
Q  175806402118016161687545467551367
LP 1
LQ 22

A simple n+1 style proof using BLS75 theorem 15. This block verifies if: a Q is odd b Q > 2 c Q divides N+1 . Let M = (N+1)/Q d MQ-1 = N e M > 0 f 2Q-1 > sqrt(N) . Let D = LP*LP - 4*LQ g D != 0 h Jacobi(D,N) = -1 . Note: V_{k} indicates the Lucas V sequence with LP,LQ i V_{m/2} mod N != 0 j V_{(N+1)/2} mod N == 0

BLS5

Type BLS5
N  8087094497428743437627091507362881
Q[1]  98277749
Q[2]  3631
A[0]  11
----

A more sophisticated n-1 proof using BLS theorem 5. This requires N-1 to be factored only to (N/2)^(1/3). While this looks much more complicated, it really isn't much more work. The biggest drawback is just that we have multiple Q values to chain rather than a single one. This block verifies if:

a  N > 2
b  N is odd
.  Note: the block terminates on the first line starting with a C<->.
.  Let Q[0] = 2
.  Let A[i] = 2 if Q[i] exists and A[i] does not
c  For each i (0 .. maxi):
c1   Q[i] > 1
c2   Q[i] < N-1
c3   A[i] > 1
c4   A[i] < N
c5   Q[i] divides N-1
. Let F = N-1 divided by each Q[i] as many times as evenly possible
. Let R = (N-1)/F
d  F is even
e  gcd(F, R) = 1
. Let s = integer    part of R / 2F
. Let f = fractional part of R / 2F
. Let P = (F+1) * (2*F*F + (r-1)*F + 1)
f  n < P
g  s = 0  OR  r^2-8s is not a perfect square
h  For each i (0 .. maxi):
h1   A[i]^(N-1) mod N = 1
h2   gcd(A[i]^((N-1)/Q[i])-1, N) = 1

ECPP

Type ECPP
N  175806402118016161687545467551367
A  96642115784172626892568853507766
B  111378324928567743759166231879523
M  175806402118016177622955224562171
Q  2297612322987260054928384863
X  3273750212
Y  82061726986387565872737368000504

An elliptic curve primality block, typically generated with an Atkin/Morain ECPP implementation, but this should be adequate for anything using the Atkin-Goldwasser-Kilian-Morain style certificates. Some basic elliptic curve math is needed for these. This block verifies if:

.  Note: A and B are allowed to be negative, with -1 not uncommon.
.  Let A = A % N
.  Let B = B % N
a  N > 0
b  gcd(N, 6) = 1
c  gcd(4*A^3 + 27*B^2, N) = 1
d  Y^2 mod N = X^3 + A*X + B mod N
e  M >= N - 2*sqrt(N) + 1
f  M <= N + 2*sqrt(N) + 1
g  Q > (N^(1/4)+1)^2
h  Q < N
i  M != Q
j  Q divides M
.  Note: EC(A,B,N,X,Y) is the point (X,Y) on Y^2 = X^3 + A*X + B, mod N
.        All values work in affine coordinates, but in theory other
.        representations work just as well.
.  Let POINT1 = (M/Q) * EC(A,B,N,X,Y)
.  Let POINT2 = M * EC(A,B,N,X,Y)  [ = Q * POINT1 ]
k  POINT1 is not the identity
l  POINT2 is the identity

is_aks_prime

say "$n is definitely prime" if is_aks_prime($n);

Takes a positive number as input, and returns 1 if the input passes the Agrawal-Kayal-Saxena (AKS) primality test. This is a deterministic unconditional primality test which runs in polynomial time for general input.

This function is only included for completeness and as an example. The Perl implementation is fast compared to the only other Perl implementation available (in Math::Primality), and the implementation in Math::Prime::Util::GMP compares favorably to others in the literature. However AKS in general is far too slow to be of practical use. R.P. Brent, 2010: "AKS is not a practical algorithm. ECPP is much faster."

lucas_sequence

my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k)

Computes U_k, V_k, and Q_k for the Lucas sequence defined by P,Q, modulo n. The modular Lucas sequence is used in a number of primality tests and proofs.

The following conditions must hold: - D = P*P - 4*Q != 0 - P > 0 - P < n - Q < n - k >= 0 - n >= 2

moebius

say "$n is square free" if moebius($n) != 0;
$sum += moebius($_) for (1..200); say "Mertens(200) = $sum";

Returns μ(n), the Möbius function (also called the Moebius, Mobius, or MoebiusMu function) for a non-negative integer input. This function is 1 if n = 1, 0 if n is not square free (i.e. n has a repeated factor), and -1^t if n is a product of t distinct primes. This is an important function in prime number theory. Like SAGE, we define moebius(0) = 0 for convenience.

If called with two arguments, they define a range low to high, and the function returns an array with the value of the Möbius function for every n from low to high inclusive. Large values of high will result in a lot of memory use. The algorithm used is Deléglise and Rivat (1996) algorithm 4.1, which is a segmented version of Lioen and van de Lune (1994) algorithm 3.2.

mertens

say "Mertens(10M) = ", mertens(10_000_000);   # = 1037

Returns M(n), the Mertens function for a non-negative integer input. This function is defined as sum(moebius(1..n)), but calculated more efficiently for large inputs. For example, computing Mertens(100M) takes:

time    approx mem
  0.4s      0.1MB   mertens(100_000_000)
 74.8s   7000MB     List::Util::sum(moebius(1,100_000_000))
 88.5s      0MB     $sum += moebius($_) for 1..100_000_000   [-nobigint]
181.8s      0MB     $sum += moebius($_) for 1..100_000_000

The summation of individual terms via factoring is quite expensive in time, though uses O(1) space. This function will generate the equivalent output via a sieving method, which will use some more memory, but be much faster. The current method is a simple n^1/2 version of Deléglise and Rivat (1996), which involves calculating all moebius values to n^1/2, which in turn will require prime sieving to n^1/4.

Various algorithms exist for this, using differing quantities of μ(n). The simplest way is to efficiently sum all n values. Benito and Varona (2008) show a clever and simple method that only requires n/3 values. Deléglise and Rivat (1996) describe a segmented method using only n^1/3 values. The current implementation does a simple non-segmented n^1/2 version of their method. Kuznetsov (2011) gives an alternate method that he indicates is even faster. Lastly, one of the advanced prime count algorithms could be theoretically used to create a faster solution.

euler_phi

say "The Euler totient of $n is ", euler_phi($n);

Returns φ(n), the Euler totient function (also called Euler's phi or phi function) for an integer value. This is an arithmetic function that counts the number of positive integers less than or equal to n that are relatively prime to n. Given the definition used, euler_phi will return 0 for all n < 1. This follows the logic used by SAGE. Mathematic/WolframAlpha also returns 0 for input 0, but returns euler_phi(-n) for n < 0.

If called with two arguments, they define a range low to high, and the function returns an array with the totient of every n from low to high inclusive. Large values of high will result in a lot of memory use.

jordan_totient

say "Jordan's totient J_$k($n) is ", jordan_totient($k, $n);

Returns Jordan's totient function for a given integer value. Jordan's totient is a generalization of Euler's totient, where jordan_totient(1,$n) == euler_totient($n) This counts the number of k-tuples less than or equal to n that form a coprime tuple with n. As with euler_phi, 0 is returned for all n < 1. This function can be used to generate some other useful functions, such as the Dedikind psi function, where psi(n) = J(2,n) / J(1,n).

exp_mangoldt

say "exp(lambda($_)) = ", exp_mangoldt($_) for 1 .. 100;

Returns EXP(Λ(n)), the exponential of the Mangoldt function (also known as von Mangoldt's function) for an integer value. It is equal to log p if n is prime or a power of a prime, and 0 otherwise. We return the exponential so all results are integers. Hence the return value for exp_mangoldt is:

p   if n = p^m for some prime p and integer m >= 1
1   otherwise.

liouville

Returns λ(n), the Liouville function for a non-negative integer input. This is -1 raised to Ω(n) (the total number of prime factors).

chebyshev_theta

say chebyshev_theta(10000);

Returns θ(n), the first Chebyshev function for a non-negative integer input. This is the sum of the logarithm of each prime where p <= n. An alternate computation is as the logarithm of n primorial. Hence these functions:

use List::Util qw/sum/;  use Math::BigFloat;

sub c1a { 0+sum( map { log($_) } @{primes(shift)} ) }
sub c1b { Math::BigFloat->new(primorial(shift))->blog }

yield similar results, albeit slower and using more memory.

chebyshev_psi

say chebyshev_psi(10000);

Returns ψ(n), the second Chebyshev function for a non-negative integer input. This is the sum of the logarithm of each prime where p^k <= n for an integer k. An alternate computation is as the summatory Mangoldt function. Another alternate computation is as the logarithm of LCM(1,2,...,n). Hence these functions:

use List::Util qw/sum/;  use Math::BigFloat;

sub c2a { 0+sum( map { log(exp_mangoldt($_)) } 1 .. shift ) }
sub c2b { Math::BigFloat->new(consecutive_integer_lcm(shift))->blog }

yield similar results, albeit slower and using more memory.

divisor_sum

say "Sum of divisors of $n:", divisor_sum( $n );

This function takes a positive integer as input and returns the sum of the k-th powers of the divisors of the input, including 1 and itself. If the second argument (k) is omitted it is assumed to be 1. This is known as the sigma function (see Hardy and Wright section 16.7, or OEIS A000203). The API is identical to Pari/GP's sigma function.

The second argument can be a code reference, which is called for each divisor and the results are summed. This allows computation of other functions, but will be less efficient than using the numeric second argument.

An example of the 5th Jordan totient (OEIS A059378):

divisor_sum( $n, sub { my $d=shift; $d**5 * moebius($n/$d); } );

though we have a function "jordan_totient" which is more efficient.

This function is useful for calculating things like aliquot sums, abundant numbers, perfect numbers, etc.

For numeric second arguments (sigma computations), the result will be a bigint if necessary. For the code reference case, the user must take care to return bigints if overflow will be a concern.

primorial

$prim = primorial(11); #        11# = 2*3*5*7*11 = 2310

Returns the primorial n# of the positive integer input, defined as the product of the prime numbers less than or equal to n. This is the OEIS series A034386: primorial numbers second definition.

primorial(0)  == 1
primorial($n) == pn_primorial( prime_count($n) )

The result will be a Math::BigInt object if it is larger than the native bit size.

Be careful about which version (primorial or pn_primorial) matches the definition you want to use. Not all sources agree on the terminology, though they should give a clear definition of which of the two versions they mean. OEIS, Wikipedia, and Mathworld are all consistent, and these functions should match that terminology. This function should return the same result as the mpz_primorial_ui function added in GMP 5.1.

pn_primorial

$prim = pn_primorial(5); #      p_5# = 2*3*5*7*11 = 2310

Returns the primorial number p_n# of the positive integer input, defined as the product of the first n prime numbers (compare to the factorial, which is the product of the first n natural numbers). This is the OEIS series A002110: primorial numbers first definition.

pn_primorial(0)  == 1
pn_primorial($n) == primorial( nth_prime($n) )

The result will be a Math::BigInt object if it is larger than the native bit size.

consecutive_integer_lcm

$lcm = consecutive_integer_lcm($n);

Given an unsigned integer argument, returns the least common multiple of all integers from 1 to n. This can be done by manipulation of the primes up to n, resulting in much faster and memory-friendly results than using a factorial.

partitions

Calculates the partition function p(n) for a non-negative integer input. This is the number of ways of writing the integer n as a sum of positive integers, without restrictions. This corresponds to Pari's numbpart function and Mathematica's PartitionsP function. The values produced in order are OEIS series A000041.

This uses a combinatorial calculation, which means it will not be very fast compared to Pari, Mathematica, or FLINT which use the Rademacher formula using multi-precision floating point. In 10 seconds, the pure Perl version can produce partitions(10000) while with Math::Prime::Util::GMP it can do partitions(200000). In contrast, in about 10 seconds Pari can solve numbpart(22000000).

If you want the enumerated partitions, see Integer::Partition. It is very fast and uses an extremely memory efficient iterator. It is not, however, practical for producing the partition number for values over 100 or so.

carmichael_lambda

Returns the Carmichael function (also called the reduced totient function, or Carmichael λ(n)) of a positive integer argument. It is the smallest positive integer m such that a^m = 1 mod n for every integer a coprime to n. This is OEIS series A002322.

znorder

$order = znorder(2, next_prime(10**19)-6);

Given two positive integers a and n, returns the multiplicative order of a modulo <n>. This is the smallest positive integer k such that a^k ≡ 1 mod n. Returns 1 if a = 1. Return undef if a = 0 or if a and n are not coprime, since no value will result in 1 mod n.

random_prime

my $small_prime = random_prime(1000);      # random prime <= limit
my $rand_prime = random_prime(100, 10000); # random prime within a range

Returns a pseudo-randomly selected prime that will be greater than or equal to the lower limit and less than or equal to the upper limit. If no lower limit is given, 2 is implied. Returns undef if no primes exist within the range.

The goal is to return a uniform distribution of the primes in the range, meaning for each prime in the range, the chances are equally likely that it will be seen. This is removes from consideration such algorithms as PRIMEINC, which although efficient, gives very non-random output. This also implies that the numbers will not be evenly distributed, since the primes are not evenly distributed. Stated again, the random prime functions return a uniformly selected prime from the set of primes within the range. Hence given random_prime(1000), the numbers 2, 3, 487, 631, and 997 all have the same probability of being returned.

For small numbers, a random index selection is done, which gives ideal uniformity and is very efficient with small inputs. For ranges larger than this ~16-bit threshold but within the native bit size, a Monte Carlo method is used (multiple calls to irand will be made if necessary). This also gives ideal uniformity and can be very fast for reasonably sized ranges. For even larger numbers, we partition the range, choose a random partition, then select a random prime from the partition. This gives some loss of uniformity but results in many fewer bits of randomness being consumed as well as being much faster.

If an irand function has been set via "prime_set_config", it will be used to construct any ranged random numbers needed. The function should return a uniformly random 32-bit integer, which is how the irand functions exported by Math::Random::Secure, Math::Random::MT, Math::Random::ISAAC, and most other modules behave.

If no irand function was set, then Bytes::Random::Secure is used with a non-blocking seed. This will create good quality random numbers, so there should be little reason to change unless one is generating long-term keys, where using the blocking random source may be preferred.

Examples of various ways to set your own irand function:

# Math::Random::Secure.  Uses ISAAC and strong seed methods.
use Math::Random::Secure;
prime_set_config(irand => \&Math::Random::Secure::irand);

# Bytes::Random::Secure (OO interface with full control of options):
use Bytes::Random::Secure ();
BEGIN {
  my $rng = Bytes::Random::Secure->new( Bits => 512 );
  sub irand { return $rng->irand; }
}
prime_set_config(irand => \&irand);

# Crypt::Random.  Uses Pari and /dev/random.  Very slow.
use Crypt::Random qw/makerandom/;
prime_set_config(irand => sub { makerandom(Size=>32, Uniform=>1); });

# Mersenne Twister.  Very fast, decent RNG, auto seeding.
use Math::Random::MT::Auto;
prime_set_config(irand=>sub {Math::Random::MT::Auto::irand() & 0xFFFFFFFF});

random_ndigit_prime

say "My 4-digit prime number is: ", random_ndigit_prime(4);

Selects a random n-digit prime, where the input is an integer number of digits between 1 and the maximum native type (10 for 32-bit, 20 for 64-bit, 10000 if bigint is active). One of the primes within that range (e.g. 1000 - 9999 for 4-digits) will be uniformly selected using the irand function as described above.

If the number of digits is greater than or equal to the maximum native type, then the result will be returned as a BigInt. However, if the -nobigint tag was used, then numbers larger than the threshold will be flagged as an error, and numbers on the threshold will be restricted to native numbers. For better performance with large bit sizes, install Math::Prime::Util::GMP.

random_nbit_prime

my $bigprime = random_nbit_prime(512);

Selects a random n-bit prime, where the input is an integer number of bits between 2 and the maximum representable bits (32, 64, or 100000 for native 32-bit, native 64-bit, and bigint respectively). A prime with the nth bit set will be uniformly selected, with randomness supplied via calls to the irand function as described above.

For bit sizes of 64 and lower, "random_prime" is used, which gives completely uniform results in this range. For sizes larger than 64, Algorithm 1 of Fouque and Tibouchi (2011) is used, wherein we select a random odd number for the lower bits, then loop selecting random upper bits until the result is prime. This allows a more uniform distribution than the general "random_prime" case while running slightly faster (in contrast, for large bit sizes "random_prime" selects a random upper partition then loops on the values within the partition, which very slightly skews the results towards smaller numbers).

The irand function is used for randomness, so all the discussion in "random_prime" about that applies here. The result will be a BigInt if the number of bits is greater than the native bit size. For better performance with large bit sizes, install Math::Prime::Util::GMP.

random_strong_prime

my $bigprime = random_strong_prime(512);

Constructs an n-bit strong prime using Gordon's algorithm. We consider a strong prime p to be one where

p is large. This function requires at least 128 bits.
p-1 has a large prime factor r.
p+1 has a large prime factor s
r-1 has a large prime factor t

Using a strong prime in cryptography guards against easy factoring with algorithms like Pollard's Rho. Rivest and Silverman (1999) present a case that using strong primes is unnecessary, and most modern cryptographic systems agree. First, the smoothness does not affect more modern factoring methods such as ECM. Second, modern factoring methods like GNFS are far faster than either method so make the point moot. Third, due to key size growth and advances in factoring and attacks, for practical purposes, using large random primes offer security equivalent to strong primes.

Similar to "random_nbit_prime", the result will be a BigInt if the number of bits is greater than the native bit size. For better performance with large bit sizes, install Math::Prime::Util::GMP.

random_proven_prime

my $bigprime = random_proven_prime(512);

Constructs an n-bit random proven prime. Internally this may use "is_provable_prime"("random_nbit_prime") or "random_maurer_prime" depending on the platform and bit size.

random_proven_prime_with_cert

my($n, $cert) = random_proven_prime_with_cert(512)

Similar to "random_proven_prime", but returns a two-element array containing the n-bit provable prime along with a primality certificate. The certificate is the same as produced by "prime_certificate" or "is_provable_prime_with_cert", and can be parsed by "verify_prime" or any other software that understands MPU primality certificates.

random_maurer_prime

my $bigprime = random_maurer_prime(512);

Construct an n-bit provable prime, using the FastPrime algorithm of Ueli Maurer (1995). This is the same algorithm used by Crypt::Primes. Similar to "random_nbit_prime", the result will be a BigInt if the number of bits is greater than the native bit size. For better performance with large bit sizes, install Math::Prime::Util::GMP.

The differences between this function and that in Crypt::Primes are described in the "SEE ALSO" section.

Internally this additionally runs the BPSW probable prime test on every partial result, and constructs a primality certificate for the final result, which is verified. These provide additional checks that the resulting value has been properly constructed.

An alternative to this function is to run "is_provable_prime" on the result of "random_nbit_prime", which will provide more diversity and will be faster up to 512 or so bits. Maurer's method should be much faster for large bit sizes (larger than 2048). If you don't need absolutely proven results, then using "random_nbit_prime" followed by additional tests ("is_strong_pseudoprime" and/or "is_frobenius_underwood_pseudoprime") should be much faster.

random_maurer_prime_with_cert

my($n, $cert) = random_maurer_prime_with_cert(512)

As with "random_maurer_prime", but returns a two-element array containing the n-bit provable prime along with a primality certificate. The certificate is the same as produced by "prime_certificate" or "is_provable_prime_with_cert", and can be parsed by "verify_prime" or any other software that understands MPU primality certificates. The proof construction consists of a single chain of BLS3 types.

UTILITY FUNCTIONS

prime_precalc

prime_precalc( 1_000_000_000 );

Let the module prepare for fast operation up to a specific number. It is not necessary to call this, but it gives you more control over when memory is allocated and gives faster results for multiple calls in some cases. In the current implementation this will calculate a sieve for all numbers up to the specified number.

prime_memfree

prime_memfree;

Frees any extra memory the module may have allocated. Like with prime_precalc, it is not necessary to call this, but if you're done making calls, or want things cleanup up, you can use this. The object method might be a better choice for complicated uses.

Math::Prime::Util::MemFree->new

my $mf = Math::Prime::Util::MemFree->new;
# perform operations.  When $mf goes out of scope, memory will be recovered.

This is a more robust way of making sure any cached memory is freed, as it will be handled by the last MemFree object leaving scope. This means if your routines were inside an eval that died, things will still get cleaned up. If you call another function that uses a MemFree object, the cache will stay in place because you still have an object.

prime_get_config

my $cached_up_to = prime_get_config->{'precalc_to'};

Returns a reference to a hash of the current settings. The hash is copy of the configuration, so changing it has no effect. The settings include:

precalc_to      primes up to this number are calculated
maxbits         the maximum number of bits for native operations
xs              0 or 1, indicating the XS code is available
gmp             0 or 1, indicating GMP code is available
maxparam        the largest value for most functions, without bigint
maxdigits       the max digits in a number, without bigint
maxprime        the largest representable prime, without bigint
maxprimeidx     the index of maxprime, without bigint
assume_rh       whether to assume the Riemann hypothesis (default 0)

prime_set_config

prime_set_config( assume_rh => 1 );

Allows setting of some parameters. Currently the only parameters are:

xs              Allows turning off the XS code, forcing the Pure Perl
                code to be used.  Set to 0 to disable XS, set to 1 to
                re-enable.  You probably will never want to do this.

gmp             Allows turning off the use of L<Math::Prime::Util::GMP>,
                which means using Pure Perl code for big numbers.  Set
                to 0 to disable GMP, set to 1 to re-enable.
                You probably will never want to do this.

assume_rh       Allows functions to assume the Riemann hypothesis is
                true if set to 1.  This defaults to 0.  Currently this
                setting only impacts prime count lower and upper
                bounds, but could later be applied to other areas such
                as primality testing.  A later version may also have a
                way to indicate whether no RH, RH, GRH, or ERH is to
                be assumed.

irand           Takes a code ref to an irand function returning a
                uniform number between 0 and 2**32-1.  This will be
                used for all random number generation in the module.

FACTORING FUNCTIONS

factor

my @factors = factor(3_369_738_766_071_892_021);
# returns (204518747,16476429743)

Produces the prime factors of a positive number input, in numerical order. The special cases of n = 0 and n = 1 will return n, which guarantees multiplying the factors together will always result in the input value, though those are the only cases where the returned factors are not prime.

In scalar context, returns Ω(n), the total number of prime factors (OEIS A001222). This corresponds to Pari's bigomega(n) function and Mathematica's PrimeOmega[n] function. This is same result that we would get if we evaluated the resulting array in scalar context. Do note that the inputs of 0 and 1 will return 1, contrary to the standard definition of Omega.

The current algorithm for non-bigints is a sequence of small trial division, a few rounds of Pollard's Rho, SQUFOF, Pollard's p-1, Hart's OLF, a long run of Pollard's Rho, and finally trial division if anything survives. This process is repeated for each non-prime factor. In practice, it is very rare to require more than the first Rho + SQUFOF to find a factor, and I have not seen anything go to the last step.

Factoring bigints works with pure Perl, and can be very handy on 32-bit machines for numbers just over the 32-bit limit, but it can be very slow for "hard" numbers. Installing the Math::Prime::Util::GMP module will speed up bigint factoring a lot, and all future effort on large number factoring will be in that module. If you do not have that module for some reason, use the GMP or Pari version of bigint if possible (e.g. use bigint try => 'GMP,Pari'), which will run 2-3x faster (though still 100x slower than the real GMP code).

factor_exp

my @factor_exponent_pairs = factor_exp(29513484000);
# returns ([2,5], [3,4], [5,3], [7,2], [11,1], [13,2])
# factor(29513484000)
# returns (2,2,2,2,2,3,3,3,3,5,5,5,7,7,11,13,13)

Produces pairs of prime factors and exponents in numerical factor order. This may be more convenient for some algorithms. This is the same form that Mathematica's FactorInteger[n] function returns.

In scalar context, returns ω(n), the number of unique prime factors (OEIS A001221). This corresponds to Pari's omega(n) function and Mathematica's PrimeNu[n] function. This is same result that we would get if we evaluated the resulting array in scalar context. Do note that the inputs of 0 and 1 will return 1, contrary to the standard definition of omega.

The internals are identical to "factor", so all comments there apply. Just the way the factors are arranged is different.

all_factors

my @divisors = all_factors(30);   # returns (1, 2, 3, 5, 6, 10, 15, 30)

Produces all the divisors of a positive number input, including 1 and the input number. The divisors are a power set of multiplications of the prime factors, returned as a uniqued sorted list. The result is identical to that of Pari's divisors function.

In scalar context this returns the sigma0 function, the sigma function (see Hardy and Wright section 16.7, or OEIS A000203).

trial_factor

my @factors = trial_factor($n);

Produces the prime factors of a positive number input. The factors will be in numerical order. The special cases of n = 0 and n = 1 will return n, while with all other inputs the factors are guaranteed to be prime. For large inputs this will be very slow.

fermat_factor

my @factors = fermat_factor($n);

Produces factors, not necessarily prime, of the positive number input. The particular algorithm is Knuth's algorithm C. For small inputs this will be very fast, but it slows down quite rapidly as the number of digits increases. It is very fast for inputs with a factor close to the midpoint (e.g. a semiprime p*q where p and q are the same number of digits).

holf_factor

my @factors = holf_factor($n);

Produces factors, not necessarily prime, of the positive number input. An optional number of rounds can be given as a second parameter. It is possible the function will be unable to find a factor, in which case a single element, the input, is returned. This uses Hart's One Line Factorization with no premultiplier. It is an interesting alternative to Fermat's algorithm, and there are some inputs it can rapidly factor. In the long run it has the same advantages and disadvantages as Fermat's method.

squfof_factor

rsqufof_factor

my @factors = squfof_factor($n);
my @factors = rsqufof_factor($n);  # racing multiplier version

prho_factor

pbrent_factor

my @factors = prho_factor($n);
my @factors = pbrent_factor($n);

# Use a very small number of rounds
my @factors = prho_factor($n, 1000);

Produces factors, not necessarily prime, of the positive number input. An optional number of rounds can be given as a second parameter. These attempt to find a single factor using Pollard's Rho algorithm, either the original version or Brent's modified version. These are more specialized algorithms usually used for pre-factoring very large inputs, as they are very fast at finding small factors.

pminus1_factor

my @factors = pminus1_factor($n);
my @factors = pminus1_factor($n, 1_000);          # set B1 smoothness
my @factors = pminus1_factor($n, 1_000, 50_000);  # set B1 and B2

Produces factors, not necessarily prime, of the positive number input. This is Pollard's p-1 method, using two stages with default smoothness settings of 1_000_000 for B1, and 10 * B1 for B2. This method can rapidly find a factor p of n where p-1 is smooth (it has no large factors).

pplus1_factor

my @factors = pplus1_factor($n);
my @factors = pplus1_factor($n, 1_000);          # set B1 smoothness

Produces factors, not necessarily prime, of the positive number input. This is Williams' p+1 method, using one stage and two predefined initial points.

MATHEMATICAL FUNCTIONS

ExponentialIntegral

my $Ei = ExponentialIntegral($x);

Given a non-zero floating point input x, this returns the real-valued exponential integral of x, defined as the integral of e^t/t dt from -infinity to x.

If the bignum module has been loaded, all inputs will be treated as if they were Math::BigFloat objects.

For non-BigInt/BigFloat objects, the result should be accurate to at least 14 digits.

For BigInt / BigFloat objects, we first check to see if Math::MPFR is available. If so, then it is used since it is very fast and has high accuracy. Accuracy when using MPFR will be equal to the accuracy() value of the input (or the default BigFloat accuracy, which is 40 by default).

MPFR is used for positive inputs only. If Math::MPFR is not available or the input is negative, then other methods are used: continued fractions (x < -1), rational Chebyshev approximation ( -1 < x < 0), a convergent series (small positive x), or an asymptotic divergent series (large positive x). Accuracy should be at least 14 digits.

LogarithmicIntegral

my $li = LogarithmicIntegral($x)

Given a positive floating point input, returns the floating point logarithmic integral of x, defined as the integral of dt/ln t from 0 to x. If given a negative input, the function will croak. The function returns 0 at x = 0, and -infinity at x = 1.

This is often known as li(x). A related function is the offset logarithmic integral, sometimes known as Li(x) which avoids the singularity at 1. It may be defined as Li(x) = li(x) - li(2). Crandall and Pomerance use the term li0 for this function, and define li(x) = Li0(x) - li0(2). Due to this terminology confusion, it is important to check which exact definition is being used.

If the bignum module has been loaded, all inputs will be treated as if they were Math::BigFloat objects.

For non-BigInt/BigFloat objects, the result should be accurate to at least 14 digits.

For BigInt / BigFloat objects, we first check to see if Math::MPFR is available. If so, then it is used, as it will return results much faster and can be more accurate. Accuracy when using MPFR will be equal to the accuracy() value of the input (or the default BigFloat accuracy, which is 40 by default).

MPFR is used for inputs greater than 1 only. If Math::MPFR is not installed or the input is less than 1, results will be calculated as Ei(ln x).

RiemannZeta

my $z = RiemannZeta($s);

Given a floating point input s where s >= 0, returns the floating point value of ζ(s)-1, where ζ(s) is the Riemann zeta function. One is subtracted to ensure maximum precision for large values of s. The zeta function is the sum from k=1 to infinity of 1 / k^s. This function only uses real arguments, so is basically the Euler Zeta function.

If the bignum module has been loaded, all inputs will be treated as if they were Math::BigFloat objects.

For non-BigInt/BigFloat objects, the result should be accurate to at least 14 digits. The XS code uses a rational Chebyshev approximation between 0.5 and 5, and a series for other values. The PP code uses an identical series for all values.

For BigInt / BigFloat objects, we first check to see if the Math::MPFR module is installed. If so, then it is used, as it will return results much faster and can be more accurate. Accuracy when using MPFR will be equal to the accuracy() value of the input (or the default BigFloat accuracy, which is 40 by default).

If Math::MPFR is not installed, then results are calculated using either Borwein (1991) algorithm 2, or the basic series. Full input accuracy is attempted, but Math::BigFloat RT 43692 produces incorrect high-accuracy computations without the fix. It is also very slow. I highly recommend installing Math::MPFR for BigFloat computations.

RiemannR

my $r = RiemannR($x);

Given a positive non-zero floating point input, returns the floating point value of Riemann's R function. Riemann's R function gives a very close approximation to the prime counting function.

If the bignum module has been loaded, all inputs will be treated as if they were Math::BigFloat objects.

For non-BigInt/BigFloat objects, the result should be accurate to at least 14 digits.

EXAMPLES

Print strong pseudoprimes to base 17 up to 10M:

# Similar to A001262's isStrongPsp function, but over 4x faster
perl -MMath::Prime::Util=:all -E 'my $n=3; while($n <= 10000000) { print "$n " if is_strong_pseudoprime($n,$base) && !is_prime($n); $n+=2; } BEGIN {$|=1; $base=17}'

or, slightly faster, use forprimes and loop over the odds between primes:

perl -MMath::Prime::Util=:all -E '$|=1; $base=17; my $prev = 1; forprimes { $prev += 2; while ($prev < $_) { print "$prev " if is_strong_pseudoprime($prev,$base); $prev += 2; } } 3,10000000'

Print some primes above 64-bit range:

perl -MMath::Prime::Util=:all -Mbigint -E 'my $start=100000000000000000000; say join "\n", @{primes($start,$start+1000)}'
# Similar code using Pari:
# perl -MMath::Pari=:int,PARI,nextprime -E 'my $start = PARI "100000000000000000000"; my $end = $start+1000; my $p=nextprime($start); while ($p <= $end) { say $p; $p = nextprime($p+1); }'

Examining the η3(x) function of Planat and Solé (2011):

sub nu3 {
  my $n = shift;
  my $phix = chebyshev_psi($n);
  my $nu3 = 0;
  foreach my $nu (1..3) {
    $nu3 += (moebius($nu)/$nu)*LogarithmicIntegral($phix**(1/$nu));
  }
  return $nu3;
}
say prime_count(1000000);
say prime_count_approx(1000000);
say nu3(1000000);

Construct and use a Sophie-Germain prime iterator:

sub make_sophie_germain_iterator {
  my $p = shift || 2;
  my $it = prime_iterator($p);
  return sub {
    do { $p = $it->() } while !is_prime(2*$p+1);
    $p;
  };
}
my $sgit = make_sophie_germain_iterator();
for (1 .. 10000) {
  print $sgit->(), "\n";
}

Project Euler, problem 3 (Largest prime factor):

use Math::Prime::Util qw/factor/;
use bigint;  # Only necessary for 32-bit machines.
say "", (factor(600851475143))[-1]

Project Euler, problem 7 (10001st prime):

use Math::Prime::Util qw/nth_prime/;
say nth_prime(10_001);

Project Euler, problem 10 (summation of primes):

use Math::Prime::Util qw/primes/;
my $sum = 0;
$sum += $_ for @{primes(2_000_000)};
say $sum;

Project Euler, problem 21 (Amicable numbers):

use Math::Prime::Util qw/divisor_sum/;
sub dsum { my $n = shift; divisor_sum($n) - $n; }
my $sum = 0;
foreach my $a (1..10000) {
  my $b = dsum($a);
  $sum += $a + $b if $b > $a && dsum($b) == $a;
}
say $sum;

Project Euler, problem 41 (Pandigital prime), brute force command line:

perl -MMath::Prime::Util=:all -E 'my @p = grep { /1/&&/2/&&/3/&&/4/&&/5/&&/6/&&/7/} @{primes(1000000,9999999)}; say $p[-1];'

Project Euler, problem 47 (Distinct primes factors):

use Math::Prime::Util qw/pn_primorial factor/;
use List::MoreUtils qw/distinct/;
sub nfactors { scalar distinct factor(shift); }
my $n = pn_primorial(4);  # Start with the first 4-factor number
$n++ while (nfactors($n) != 4 || nfactors($n+1) != 4 || nfactors($n+2) != 4 || nfactors($n+3) != 4);
say $n;

Project Euler, problem 69, stupid brute force solution (about 2 seconds):

use Math::Prime::Util qw/euler_phi/;
my ($n, $max) = (0,0);
do {
  my $ndivphi = $_ / euler_phi($_);
  ($n, $max) = ($_, $ndivphi) if $ndivphi > $max;
} for 1..1000000;
say "$n  $max";

Here is the right way to do PE problem 69 (under 0.03s):

use Math::Prime::Util qw/pn_primorial/;
my $n = 0;
$n++ while pn_primorial($n+1) < 1000000;
say pn_primorial($n);

Project Euler, problem 187, stupid brute force solution, ~4 minutes:

use Math::Prime::Util qw/factor -nobigint/;
my $nsemis = 0;
do { $nsemis++ if scalar factor($_) == 2; }
   for 1 .. int(10**8)-1;
say $nsemis;

Here is the best way for PE187. Under 30 milliseconds from the command line:

use Math::Prime::Util qw/primes prime_count/;
use List::Util qw/sum/;
my $limit = shift || int(10**8);
my @primes = @{primes(int(sqrt($limit)))};
say sum( map { prime_count(int(($limit-1)/$primes[$_-1])) - $_ + 1 }
             1 .. scalar @primes );

PRIMALITY TESTING NOTES

Above 2^64, "is_prob_prime" performs an extra-strong BPSW test which is fast (a little less than the time to perform 3 Miller-Rabin tests) and has no known counterexamples. If you believe the primality testing done by Pari, Maple, SAGE, FLINT, etc., then this function should be appropriate for you. "is_prime" will do the same BPSW test as well as some additional testing, making it slightly more time consuming but less likely to produce a false result. This is a little more stringent than Mathematica. "is_provable_prime" constructs a primality proof. If a certificate is requested, then either BLS75 theorem 5 or ECPP is performed. Without a certificate, the method is implementation specific (currently it is identical, but later releases may use APRCL). With Math::Prime::Util::GMP installed, this is quite fast through 300 or so digits.

Math systems 30 years ago used to use Miller-Rabin tests with k bases (typically fixed bases, sometimes random) for primality testing, but these have generally been replaced by some form of BPSW as used in this module. See Pinch's 1993 paper for examples of why using k M-R tests leads to poor results. The three exceptions in common contemporary use I am aware of are:

libtommath: Uses the first k prime bases. This is unacceptable for cryptographic use, as there are known methods (e.g. Arnault 1994) for constructing counterexamples. The number of bases required to avoid false results is unreasonably high, hence performance is slow even if one ignores counterexamples. Unfortunately this is the multi-precision math library used for Perl 6 and at least one CPAN Crypto module.
GMP/MPIR: Uses a set of k static-random bases. The bases are randomly chosen using a PRNG that is seeded identically each call (the seed changes with each release). This offers a very slight advantage over using the first k prime bases, but not much. See, for example, Nicely's mpz_probab_prime_p pseudoprimes page.
Math::Pari: Pari 2.1.7 is the default version installed with the Math::Pari module. It uses 10 random M-R bases (the PRNG uses a fixed seed set at compile time). Pari 2.3.0 was released in May 2006 and it, like all later releases through at least 2.6.1, use BPSW / APRCL, after complaints of false results from using M-R tests.

Basically the problem is that it is just too easy to get counterexamples from running k M-R tests, forcing one to use a very large number of tests (at least 20) to avoid frequent false results. Using the BPSW test results in no known counterexamples after 30+ years and runs much faster. It can be enhanced with one or more random bases if one desires, and will still be much faster.

Using k fixed bases has another problem, which is that in any adversarial situation we can assume the inputs will be selected such that they are one of our counterexamples. Now we need absurdly large numbers of tests. This is like playing "pick my number" but the number is fixed forever at the start, the guesser gets to know everyone else's guesses and results who has every played, and can keep playing as long as they like. It's only valid if the players are completely oblivious to what is happening.

LIMITATIONS

Perl versions earlier than 5.8.0 have a rather broken 64-bit implementation, in that the values are accessed as doubles. Hence any value larger than ~ 2^49 will start losing bottom bits. This causes numerous functions to not work properly. The test suite will try to determine if your Perl is broken (this only applies to really old versions of Perl compiled for 64-bit when using numbers larger than ~ 2^49). The best solution is updating to a more recent Perl.

The module is thread-safe and should allow good concurrency on all platforms that support Perl threads except Win32. With Win32, either don't use threads or make sure prime_precalc is called before using primes, prime_count, or nth_prime with large inputs. This is only an issue if you use non-Cygwin Win32 and call these routines from within Perl threads.

PERFORMANCE

First, for those looking for the state of the art non-Perl solutions:

Primality testing

PFGW is the fastest primality testing software I'm aware of once past 2000 or so digits, has fast trial division, and is especially fast on many special forms. It does not have a BPSW test however, and there are quite a few counterexamples for a given base of its PRP test, so for primality testing it is most useful for fast filtering of very large candidates. A test such as the BPSW test in this module is then recommended.

Primality proofs

Primo is the best method for open source primality proving for inputs over 1000 digits. Primo also does well below that size, but other good alternatives are WraithX APRCL, the APRCL from the modern Pari package, or the standalone ECPP from this module with large polynomial set.

Factoring

yafu, msieve, and gmp-ecm are all good choices for large inputs. The factoring code in this module (and all other CPAN modules) is very limited compared to those.

Primes

primesieve is the fastest publically available code I am aware of. It is much faster than any of the alternatives, and even more so when run multi-threaded. Tomás Oliveira e Silva's private code may be faster for very large values, but it isn't available for testing.

Note that the Sieve of Atkin is not faster than the Sieve of Eratosthenes when both are well implemented. The only Sieve of Atkin that is even competitive is Bernstein's super optimized primegen, which runs about 10% faster than the simple SoE in this module, slower than Pari and yafu's SoE implementations, and 2x slower than primesieve.

Prime Counts and Nth Prime

Up to a limit, extensive use of tables plus a good segmented sieve will produce the fastest results, but the number of tables needed to maintain good performance grows exponentially. The code in this module approaches the best publically available results, with the notable exception of Christian Bau's L-M-O implementation. The author of primesieve is also working on an L-M-O implementation. None of these are state of the art compared to private research methods.

PRIME COUNTS

Counting the primes to 10^10 (10 billion), with time in seconds. Pi(10^10) = 455,052,511. The numbers below are for sieving. Calculating Pi(10^10) takes 0.064 seconds using the Lehmer algorithm in version 0.12.

External C programs in C / C++:

    1.9  primesieve 3.6 forced to use only a single thread
    2.2  yafu 1.31
    3.8  primegen (optimized Sieve of Atkin, conf-word 8192)
    5.6  Tomás Oliveira e Silva's unoptimized segmented sieve v2 (Sep 2010)
    6.7  Achim Flammenkamp's prime_sieve (32k segments)
    9.3  http://tverniquet.com/prime/ (mod 2310, single thread)
   11.2  Tomás Oliveira e Silva's unoptimized segmented sieve v1 (May 2003)
   17.0  Pari 2.3.5 (primepi)

Small portable functions suitable for plugging into XS:

    4.1  My segmented SoE used in this module (with unrolled inner loop)
   15.6  My Sieve of Eratosthenes using a mod-30 wheel
   17.2  A slightly modified verion of Terje Mathisen's mod-30 sieve
   35.5  Basic Sieve of Eratosthenes on odd numbers
   33.4  Sieve of Atkin, from Praxis (not correct)
   72.8  Sieve of Atkin, 10-minute fixup of basic algorithm
   91.6  Sieve of Atkin, Wikipedia-like

Perl modules, counting the primes to 800_000_000 (800 million):

Time (s)   Module                      Version  Notes
---------  --------------------------  -------  -----------
     0.002 Math::Prime::Util           0.35     using extended LMO
     0.007 Math::Prime::Util           0.12     using Lehmer's method
     0.27  Math::Prime::Util           0.17     segmented mod-30 sieve
     0.39  Math::Prime::Util::PP       0.24     Perl (Lehmer's method)
     0.9   Math::Prime::Util           0.01     mod-30 sieve
     2.9   Math::Prime::FastSieve      0.12     decent odd-number sieve
    11.7   Math::Prime::XS             0.26     needs some optimization
    15.0   Bit::Vector                 7.2
    48.9   Math::Prime::Util::PP       0.14     Perl (fastest I know of)
   170.0   Faster Perl sieve (net)     2012-01  array of odds
   548.1   RosettaCode sieve (net)     2012-06  simplistic Perl
  3048.1   Math::Primality             0.08     Perl + Math::GMPz
~11000     Math::Primality             0.04     Perl + Math::GMPz
>20000     Math::Big                   1.12     Perl, > 26GB RAM used

Python's standard modules are very slow: MPMATH v0.17 primepi takes 169.5s and 25+ GB of RAM. SymPy 0.7.1 primepi takes 292.2s. However there are very fast solutions written by Robert William Hanks (included in the xt/ directory of this distribution): pure Python in 12.1s and NUMPY in 2.8s.

PRIMALITY TESTING

is_prime: my impressions for various sized inputs:

Module                   1-10 digits  10-20 digits  BigInts
-----------------------  -----------  ------------  --------------
Math::Prime::Util        Very fast    Very fast     Slow / Very Fast (1)
Math::Prime::XS          Very fast    Very slow (2) --
Math::Prime::FastSieve   Very fast    N/A (3)       --
Math::Primality          Very slow    Very slow     Fast
Math::Pari               Slow         OK            OK / Fast (4)

(1) Without / With L<Math::Prime::Util::GMP> installed.
(2) Trial division only.  Very fast if every factor is tiny.
(3) Too much memory to hold the sieve (11dig = 6GB, 12dig = ~50GB)
(4) By default L<Math::Pari> installs Pari 2.1.7, which uses 10 M-R tests
    for isprime and is not fast.  See notes below for 2.3.5.

The differences are in the implementations:

Math::Prime::Util: first does simple divisibility tests to quickly recognize composites, then looks in the sieve for a fast bit lookup if possible (default up to 30,000 but can be expanded via prime_precalc). Next, for relatively small inputs, a deterministic set of Miller-Rabin tests are used, while for larger inputs a strong BPSW test is performed. For native integers, this is faster than any of the other modules. With Bigints, you need the Math::Prime::Util::GMP module installed to get good performance. With that installed, it is about 2x faster than Math::Primality and 10x faster than Math::Pari (default 2.1.7).
Math::Prime::XS: does trial divisions, which is wonderful if the input has a small factor (or is small itself). If given a large prime it can be tens of thousands of times slower than MPU. It does not support bigints.
Math::Prime::FastSieve: only works in a sieved range, which is really fast if you can do it (M::P::U will do the same if you call prime_precalc). Larger inputs just need too much time and memory for the sieve.
Math::Primality: uses GMP (in Perl) for all work. Under ~32-bits it uses 2 or 3 MR tests, while above 4759123141 it performs a BPSW test. This is great for bigints over 2^64, but it is significantly slower than native precision tests. With 64-bit numbers it is generally an order of magnitude or more slower than any of the others. Once bigints are being used, its performance is quite good. It is faster than this module unless Math::Prime::Util::GMP has been installed, in which case Math::Prime::Util is faster.
Math::Pari: has some very effective code, but it has some overhead to get to it from Perl. That means for small numbers it is relatively slow: an order of magnitude slower than M::P::XS and M::P::Util (though arguably this is only important for benchmarking since "slow" is ~2 microseconds). Large numbers transition over to smarter tests so don't slow down much. With the default Pari version, isprime will do M-R tests for 10 randomly chosen bases, but can perform a Pocklington-Lehmer proof if requested using isprime(x,1). Both could fail to identify a composite. If pari 2.3.5 is used instead (this requires hand-building the Math::Pari module) then the options are quite different. ispseudoprime(x,0) performs a strong BPSW test, while isprime now performs a primality proof using a fast implementation of the APRCL method. While the APRCL method is very fast it is still much, much slower than a BPSW probable prime test for large inputs.

FACTORING

Factoring performance depends on the input, and the algorithm choices used are still being tuned. Math::Factor::XS is very fast when given input with only small factors, but it slows down rapidly as the smallest factor increases in size. For numbers larger than 32 bits, Math::Prime::Util can be 100x or more faster (a number with only very small factors will be nearly identical, while a semiprime with large factors will be the extreme end). Math::Pari is much slower with native sized inputs, probably due to calling overhead. For bigints, the Math::Prime::Util::GMP module is needed or performance will be far worse than Math::Pari. With the GMP module, performance is pretty similar from 20 through 70 digits, which the caveat that the current MPU factoring uses more memory for 60+ digit numbers.

This slide presentation has a lot of data on 64-bit and GMP factoring performance I collected in 2009. Assuming you do not know anything about the inputs, trial division and optimized Fermat or Lehman work very well for small numbers (<= 10 digits), while native SQUFOF is typically the method of choice for 11-18 digits (I've seen claims that a lightweight QS can be faster for 15+ digits). Some form of Quadratic Sieve is usually used for inputs in the 19-100 digit range, and beyond that is the General Number Field Sieve. For serious factoring, I recommend looking at yafu, msieve, gmp-ecm, GGNFS, and Pari. The latest yafu should cover most uses, with GGNFS likely only providing a benefit for numbers large enough to warrant distributed processing.

PRIMALITY PROVING

The n-1 proving algorithm in Math::Prime::Util::GMP compares well to the version including in Pari. Both are pretty fast to about 60 digits, and work reasonably well to 80 or so before starting to take over many minutes per number on a fast computer. Version 0.09 and newer of MPU::GMP contain an ECPP implementation that, while not state of the art compared to closed source solutions, works quite well. It averages less than a second for proving 200-digit primes including creating a certificate. Times below 200 digits are faster than Pari 2.3.5's APR-CL proof. For larger inputs the bottleneck is a limited set of discriminants, and time becomes more variable. There is a larger set of discriminants on github that help, with 300-digit primes taking ~5 seconds on average and typically under a minute for 500-digits. For primality proving with very large numbers, I recommend Primo.

RANDOM PRIME GENERATION

Seconds per prime for random prime generation on a circa-2009 workstation, with Math::BigInt::GMP, Math::Prime::Util::GMP, and Math::Random::ISAAC::XS installed.

bits    random   +testing  rand_prov   Maurer   CPMaurer
-----  --------  --------  ---------  --------  --------
   64    0.0001  +0.000008   0.0002     0.0001    0.022
  128    0.0020  +0.00023    0.011      0.063     0.057
  256    0.0034  +0.0004     0.058      0.13      0.16
  512    0.0097  +0.0012     0.28       0.28      0.41
 1024    0.060   +0.0060     0.65       0.65      2.19
 2048    0.57    +0.039      4.8        4.8      10.99
 4096    6.24    +0.25      31.9       31.9      79.71
 8192   58.6     +1.61     234.0      234.0     947.3

random    = random_nbit_prime  (results pass BPSW)
random+   = additional time for 3 M-R and a frobenius test
rand_prov = random_proven_prime
maurer    = random_maurer_prime
CPMaurer  = Crypt::Primes::maurer

"random_nbit_prime" is reasonably fast, and for most purposes should suffice. For cryptographic purposes, one may want additional tests or a proven prime. Additional tests are quite cheap, as shown by the time for three extra M-R and a Frobenius test. At these bit sizes, the chances a composite number passes BPSW, three more M-R tests, and a Frobenius test is extraordinarily small.

"random_proven_prime" provides a randomly selected prime with an optional certificate, without specifying the particular method. Below 512 bits, using "is_provable_prime"("random_nbit_prime") is typically faster than Maurer's algorithm, but becomes quite slow as the bit size increases. This leaves the decision of the exact method of proving the result to the implementation.

"random_maurer_prime" constructs a provable prime. A primality test is run on each intermediate, and it also constructs a complete primality certificate which is verified at the end (and can be returned). While the result is uniformly distributed, only about 10% of the primes in the range are selected for output. This is a result of the FastPrime algorithm and is usually unimportant.

"maurer" in Crypt::Primes times are included for comparison. It is pretty fast for small sizes but gets slow as the size increases. It does not perform any primality checks on the intermediate results or the final result (I highly recommended you run a primality test on the output). Additionally important for servers, "maurer" in Crypt::Primes uses excessive system entropy and can grind to a halt if /dev/random is exhausted (it can take days to return). The times above are on a machine running HAVEGED so never waits for entropy. Without this, the times would be much higher.

AUTHORS

Dana Jacobsen <dana@acm.org>

ACKNOWLEDGEMENTS

Eratosthenes of Cyrene provided the elegant and simple algorithm for finding primes.

Terje Mathisen, A.R. Quesada, and B. Van Pelt all had useful ideas which I used in my wheel sieve.

Tomás Oliveira e Silva has released the source for a very fast segmented sieve. The current implementation does not use these ideas. Future versions might.

The SQUFOF implementation being used is a slight modification to the public domain racing version written by Ben Buhrow. Enhancements with ideas from Ben's later code as well as Jason Papadopoulos's public domain implementations are planned for a later version. The old SQUFOF implementation, still included in the code, is my modifications to Ben Buhrow's modifications to Bob Silverman's code.

REFERENCES

Pierre Dusart, "Estimates of Some Functions Over Primes without R.H.", preprint, 2010. Updates to the best non-RH bounds for prime count and nth prime. http://arxiv.org/abs/1002.0442/
Pierre Dusart, "Autour de la fonction qui compte le nombre de nombres premiers", PhD thesis, 1998. In French. The mathematics is readable and highly recommended reading if you're interesting in prime number bounds. http://www.unilim.fr/laco/theses/1998/T1998_01.html
Gabriel Mincu, "An Asymptotic Expansion", Journal of Inequalities in Pure and Applied Mathematics, v4, n2, 2003. A very readable account of Cipolla's 1902 nth prime approximation. http://www.emis.de/journals/JIPAM/images/153_02_JIPAM/153_02.pdf
Christian Bau, "The Extended Meissel-Lehmer Algorithm", 2003, preprint with example C++ implementation. Very detailed implementation-specific paper which was used for the implementation here. Highly recommended for implementing a sieve-based LMO. http://cs.swan.ac.uk/~csoliver/ok-sat-library/OKplatform/ExternalSources/sources/NumberTheory/ChristianBau/
David M. Smith, "Multiple-Precision Exponential Integral and Related Functions", ACM Transactions on Mathematical Software, v37, n4, 2011. http://myweb.lmu.edu/dmsmith/toms2011.pdf
Vincent Pegoraro and Philipp Slusallek, "On the Evaluation of the Complex-Valued Exponential Integral", Journal of Graphics, GPU, and Game Tools, v15, n3, pp 183-198, 2011. http://www.cs.utah.edu/~vpegorar/research/2011_JGT/paper.pdf
William H. Press et al., "Numerical Recipes", 3rd edition.
W. J. Cody and Henry C. Thacher, Jr., "Chebyshev approximations for the exponential integral Ei(x)", Mathematics of Computation, v23, pp 289-303, 1969. http://www.ams.org/journals/mcom/1969-23-106/S0025-5718-1969-0242349-2/
W. J. Cody and Henry C. Thacher, Jr., "Rational Chebyshev Approximations for the Exponential Integral E_1(x)", Mathematics of Computation, v22, pp 641-649, 1968.
W. J. Cody, K. E. Hillstrom, and Henry C. Thacher Jr., "Chebyshev Approximations for the Riemann Zeta Function", "Mathematics of Computation", v25, n115, pp 537-547, July 1971.
Ueli M. Maurer, "Fast Generation of Prime Numbers and Secure Public-Key Cryptographic Parameters", 1995. Generating random provable primes by building up the prime. http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.26.2151
Pierre-Alain Fouque and Mehdi Tibouchi, "Close to Uniform Prime Number Generation With Fewer Random Bits", pre-print, 2011. Describes random prime distributions, their algorithm for creating random primes using few random bits, and comparisons to other methods. Definitely worth reading for the discussions of uniformity. http://eprint.iacr.org/2011/481
Douglas A. Stoll and Patrick Demichel , "The impact of ζ(s) complex zeros on π(x) for x < 10^{10^{13}}", "Mathematics of Computation", v80, n276, pp 2381-2394, October 2011. http://www.ams.org/journals/mcom/2011-80-276/S0025-5718-2011-02477-4/home.html
OEIS: Primorial
Walter M. Lioen and Jan van de Lune, "Systematic Computations on Mertens' Conjecture and Dirichlet's Divisor Problem by Vectorized Sieving", in From Universal Morphisms to Megabytes, Centrum voor Wiskunde en Informatica, pp. 421-432, 1994. Describes a nice way to compute a range of Möbius values. http://walter.lioen.com/papers/LL94.pdf
Marc Deléglise and Joöl Rivat, "Computing the summation of the Möbius function", Experimental Mathematics, v5, n4, pp 291-295, 1996. Enhances the Möbius computation in Lioen/van de Lune, and gives a very efficient way to compute the Mertens function. http://projecteuclid.org/euclid.em/1047565447
Manuel Benito and Juan L. Varona, "Recursive formulas related to the summation of the Möbius function", The Open Mathematics Journal, v1, pp 25-34, 2007. Among many other things, shows a simple formula for computing the Mertens functions with only n/3 Möbius values (not as fast as Deléglise and Rivat, but really simple). http://www.unirioja.es/cu/jvarona/downloads/Benito-Varona-TOMATJ-Mertens.pdf
John Brillhart, D. H. Lehmer, and J. L. Selfridge, "New Primality Criteria and Factorizations of 2^m +/- 1", Mathematics of Computation, v29, n130, Apr 1975, pp 620-647. http://www.ams.org/journals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf

COPYRIGHT

This program is free software; you can redistribute it and/or modify it under the same terms as Perl itself.

To install Math::Prime::Util, copy and paste the appropriate command in to your terminal.

cpanm

cpanm Math::Prime::Util

CPAN shell

perl -MCPAN -e shell
install Math::Prime::Util

For more information on module installation, please visit the detailed CPAN module installation guide.

	Global
`s`	Focus search bar
`?`	Bring up this help dialog

	GitHub
`g` `p`	Go to pull requests
`g` `i`	go to github issues (only if github is preferred repository)

	POD
`g` `a`	Go to author
`g` `c`	Go to changes
`g` `i`	Go to issues
`g` `d`	Go to dist
`g` `r`	Go to repository/SCM
`g` `s`	Go to source
`g` `b`	Go to file browse

	Search terms
module: (e.g. module:Plugin)
distribution: (e.g. distribution:Dancer auth)
author: (e.g. author:SONGMU Redis)
version: (e.g. version:1.00)