Provided by: libmath-prime-util-perl_0.37-1_amd64 
NAME
Math::Prime::Util - Utilities related to prime numbers, including fast sieves and
factoring
VERSION
Version 0.37
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;
# Or for the composites in a range
forcomposites { say if is_strong_pseudoprime($_,2) } 10000, 10**6;
# 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 = divisors( $n );
# Or just apply a block for each one
fordivisors { $sum += $_ + $_*$_ } $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.
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 <gttp://gmplib.org> 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. For example, the default BigInt backend in older versions of Perl will
sometimes convert small results to doubles, resulting in corrupted output.
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 <http://en.wikipedia.org/wiki/Baillie-
PSW_primality_test> is used. See the "PRIMALITY TESTING NOTES" section for more
discussion.
FUNCTIONS
is_prime
print "$n is prime" if is_prime($n);
Returns 0 is the number is composite, 1 if it is probably prime, and 2 if it is definitely
prime. For numbers smaller than "2^64" it will only return 0 (composite) or 2 (definitely
prime), as this range has been exhaustively tested and has no counterexamples. For larger
numbers, an extra-strong BPSW test is used. If Math::Prime::Util::GMP is installed, some
additional primality tests are also performed, 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 a different test (e.g.
"is_frobenius_underwood_pseudoprime") or using additional M-R tests with random bases with
"miller_rabin_random". 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 primality testing (for very small
ranges), 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. The result will be a bigint if it
can not be exactly represented in the native int type (larger than "4,294,967,291" in
32-bit Perl; larger than "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 use 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.
forcomposites
forcomposites { say } 1000;
forcomposites { say } 2000,2020;
Given a block and either an end number or a start and end pair, calls the block for each
composite in the inclusive range. The composites are the numbers greater than 1 which are
not prime: "4, 6, 8, 9, 10, 12, 14, 15, ..."
fordivisors
fordivisors { $prod *= $_ } $n;
Given a block and a non-negative number "n", the block is called with $_ set to each
divisor in sorted order. Also see "divisor_sum".
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 <http://oeis.org/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
<http://oeis.org/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 <http://oeis.org/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 <http://www.ams.org/mathscinet-getitem?mr=1680879>).
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 <http://oeis.org/A217719>
method: increment "P" from 3 until "jacobi(D,n) = -1". Removing primes, this produces the
sequence OEIS A217719 <http://oeis.org/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 <http://oeis.org/A217719>. A
value of 2 yields the method used by Pari <http://pari.math.u-
bordeaux.fr/faq.html#primetest>.
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_bpsw_prime
Given a positive number input, returns 0 (composite), 2 (definitely prime), or 1 (probably
prime), using the BPSW primality test (extra-strong variant). Normally one of the
"is_prime" in Math::Prime::Util or "is_prob_prime" in Math::Prime::Util functions will
suffice, but those functions do pre-tests to find easy composites. If you know this is
not necessary, then calling "is_bpsw_prime" may save a small amount of time.
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.
While this is an important theoretical algorithm, and makes an interesting example, it is
hard to overstate just how impractically slow it is in practice. It is not used for any
purpose in non-theoretical work, as it is literally millions of times slower than other
algorithms. From R.P. Brent, 2010: "AKS is not a practical algorithm. ECPP is much
faster." We have ECPP, and indeed it 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" ; " 0 < P < n" ; " Q < n" ; " k
>= 0" ; " n >= 2".
gcd
Given a list of integers, returns the greatest common divisor. This is often used to test
for coprimality <https://oeis.org/wiki/Coprimality>.
lcm
Given a list of integers, returns the least common multiple. Note that we follow the
semantics of Mathematica, Pari, and Perl 6, re:
lcm(0, n) = 0 Any zero in list results in zero return
lcm(n,-m) = lcm(n, m) We use the absolute values
moebius
say "$n is square free" if moebius($n) != 0;
$sum += moebius($_) for (1..200); say "Mertens(200) = $sum";
Returns X(n), the Moebius function (also known as the Moebius, Mobius, or MoebiusMu
function) for an 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 Moebius function for every n from low to high
inclusive. Large values of high will result in a lot of memory use. The algorithm used
for ranges is Deleglise and Rivat (1996) algorithm 4.1, which is a segmented version of
Lioen and van de Lune (1994) algorithm 3.2.
The return values are read-only constants. This should almost never come up, but it means
trying to modify aliased return values will cause an exception (modifying the returned
scalar or array is fine).
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.3s 0.1MB mertens(100_000_000)
1.2s 890MB List::Util::sum(moebius(1,100_000_000))
77s 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. Using the range version of moebius is much faster, but returns a 100M element
array which is not good for memory with this many items. In comparison, this function
will generate the equivalent output via a sieving method that is relatively sparse memory
and very fast. The current method is a simple "n^1/2" version of Deleglise 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 X(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. Deleglise 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 X(n), the Euler totient function (also called Euler's phi or phi function) for an
integer value. This is an arithmetic function which 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.
Mathematica and Pari return "euler_phi(-n)" for "n < 0". Mathematica returns 0 for "n =
0" while Pari raises an exception.
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.
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(X(n)), the exponential of the Mangoldt function (also known as von Mangoldt's
function) for an integer value. The Mangoldt function 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 X(n), the Liouville function for a non-negative integer input. This is -1 raised
to X(n) (the total number of prime factors).
chebyshev_theta
say chebyshev_theta(10000);
Returns X(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 X(n), the second Chebyshev function for a non-negative integer input. This is the
sum of the logarithm of each prime power 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 );
say "sigma_2($n) = ", divisor_sum($n, 2);
say "Number of divisors: sigma_0($n) = ", divisor_sum($n, 0);
This function takes a positive integer as input and returns the sum of its divisors,
including 1 and itself. An optional second argument "k" may be given, which will result
in the sum of the "k-th" powers of the divisors to be returned.
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. This function is useful for
calculating things like aliquot sums, abundant numbers, perfect numbers, etc.
The second argument may also 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. This corresponds to Pari/GP's "sumdiv"
function.
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.
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
<http://oeis.org/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 <http://oeis.org/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
<http://oeis.org/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:
65 Integer::Partition
10_000 MPU pure Perl partitions
200_000 MPU GMP partitions
22_000_000 Pari's numbpart
500_000_000 Jonathan Bober's partitions_c.cc v0.6
If you want the enumerated partitions, see Integer::Partition. It uses a memory efficient
iterator and is very fast for enumeration. It is not practical for producing large
partition numbers as seen above.
carmichael_lambda
Returns the Carmichael function (also called the reduced totient function, or Carmichael
X(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
<http://oeis.org/A002322>.
kronecker
Returns the Kronecker symbol "(a|n)" for two integers. The possible return values with
their meanings for odd positive "n" are:
0 a = 0 mod n
1 a is a quadratic residue modulo n (a = x^2 mod n for some x)
-1 a is a quadratic non-residue modulo n
The Kronecker symbol is an extension of the Jacobi symbol to all integer values of "n"
from the latter's domain of positive odd values of "n". The Jacobi symbol is itself an
extension of the Legendre symbol, which is only defined for odd prime values of "n". This
corresponds to Pari's "kronecker(a,n)" function and Mathematica's "KroneckerSymbol[n,m]"
function.
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 X 1 mod n". Returns 1 if
"a = 1". Returns undef if "a = 0" or if "a" and "n" are not coprime, since no value will
result in 1 mod n. This corresponds to Pari's "znorder(Mod(a,n))" function and
Mathematica's "MultiplicativeOrder[n]" function.
znprimroot
Given a positive integer "n", returns the smallest primitive root of "(Z/nZ)^*", or
"undef" if no root exists. A root exists when "euler_phi($n) == carmichael_lambda($n)",
which will be true for all prime "n" and some composites.
OEIS A033948 <http://oeis.org/A033948> is a sequence of integers where the primitive root
exists, while OEIS A046145 <http://oeis.org/A046145> is a list of the smallest primitive
roots, which is what this function produces.
znlog
$k = znlog($a, $g, $p)
Returns the integer "k" that solves the equation "a = g^k mod p", or undef if no solution
is found. This is the discrete logarithm problem. The implementation in this version is
not very useful, but may be improved.
legendre_phi
$phi = legendre_phi(1000000000, 41);
Given a non-negative integer "n" and a non-negative prime number "a", returns the Legendre
phi function (also called Legendre's sum). This is the count of positive integers <= "n"
which are not divisible by any of the first "a" primes.
RANDOM PRIMES
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 differently, 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:
# System rand. You probably don't want to do this.
prime_set_config(irand => sub { int(rand(4294967296)) });
# 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});
# Go back to MPU's default configuration
prime_set_config(irand => undef);
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. 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" configuration option is
on, then output will be restricted to native size numbers, and requests for more digits
than natively supported will result in an error. 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. 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 product of
the returned factors will be equal to the input. "n = 1" will return an empty list, and
"n = 0" will return 0. This matches Pari.
In scalar context, returns X(n), the total number of prime factors (OEIS A001222
<http://oeis.org/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.
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 is more
convenient for some algorithms. This is the same form that Mathematica's
"FactorInteger[n]" and Pari/GP's "factorint" functions return. Note that Math::Pari
transposes the Pari result matrix.
In scalar context, returns X(n), the number of unique prime factors (OEIS A001221
<http://oeis.org/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.
The internals are identical to "factor", so all comments there apply. Just the way the
factors are arranged is different.
divisors
all_factors
my @divisors = divisors(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" and
Mathematica's "Divisors[n]" functions.
In scalar context this returns the sigma0 function, the sigma function (see Hardy and
Wright section 16.7, or OEIS A000203). This is the same result as evaluating the array in
scalar context.
Also see the "for_divisors" functions for looping over the divisors.
"all_factors" is the deprecated name for this function.
trial_factor
my @factors = trial_factor($n);
Produces the prime factors of a positive number input. The factors will be in numerical
order. 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
my @factors = squfof_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 function
typically runs very fast.
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
X(s)-1, where X(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 <https://rt.cpan.org/Ticket/Display.html?id=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.
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). Accuracy without MPFR should
be 35 digits.
EXAMPLES
Print strong pseudoprimes to base 17 up to 10M:
# Similar to A001262's isStrongPsp function, but much faster
perl -MMath::Prime::Util=:all -E 'forcomposites { say if is_strong_pseudoprime($_,17) } 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)}'
# Another way
perl -MMath::Prime::Util=:all -E 'forprimes { say } "100000000000000000039", "100000000000000000993"'
# Similar using Math::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 X3(x) function of Planat and Sole (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();
print $sgit->(), "\n" for 1 .. 10000;
Project Euler, problem 3 (Largest prime factor):
use Math::Prime::Util qw/factor/;
use bigint; # Only necessary for 32-bit machines.
say 0+(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/forprimes/;
my $sum = 0;
forprimes { $sum += $_ } 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=primes -MList::Util=first -E 'say first { /1/&&/2/&&/3/&&/4/&&/5/&&/6/&&/7/} reverse @{primes(1000000,9999999)};'
Project Euler, problem 47 (Distinct primes factors):
use Math::Prime::Util qw/pn_primorial factor_exp/;
my $n = pn_primorial(4); # Start with the first 4-factor number
# factor_exp in scalar context returns the number of distinct prime factors
$n++ while (factor_exp($n) != 4 || factor_exp($n+1) != 4 || factor_exp($n+2) != 4 || factor_exp($n+3) != 4);
say $n;
Project Euler, problem 69, stupid brute force solution (about 1 second):
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, ~3 minutes:
use Math::Prime::Util qw/factor/;
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 );
Produce the "matches" result from Math::Factor::XS without skipping:
use Math::Prime::Util qw/divisors/;
use Algorithm::Combinatorics qw/combinations_with_repetition/;
my $n = 139650;
my @matches = grep { $_->[0] * $_->[1] == $n && $_->[0] > 1 }
combinations_with_repetition( [divisors($n)], 2 );
Compute OEIS A054903 <http://oeis.org/A054903> just like CRG4's Pari example:
use Math::Prime::Util qw/forcomposite divisor_sum/;
forcomposites {
say if divisor_sum($_)+6 == divisor_sum($_+6)
} 9,1e7;
Construct the table shown in OEIS A046147 <http://oeis.org/A046147>:
use Math::Prime::Util qw/znorder euler_phi gcd/;
foreach my $n (1..100) {
if (!znprimroot($n)) {
say "$n -";
} else {
my $phi = euler_phi($n);
my @r = grep { gcd($_,$n) == 1 && znorder($_,$n) == $phi } 1..$n-1;
say "$n ", join(" ", @r);
}
}
PRIMALITY TESTING NOTES
Above "2^64", "is_prob_prime" performs an extra-strong BPSW test
<http://en.wikipedia.org/wiki/Baillie-PSW_primality_test> which is fast (a little less
than the time to perform 3 Miller-Rabin tests) and has no known counterexamples. If you
trust 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 typically used Miller-Rabin tests with "k" bases (usually 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 problematic 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
<http://www.trnicely.net/misc/mpzspsp.html> 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, 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 problems doing exact integer math. Some operations
will flip signs, and many operations will convert intermediate or output results to
doubles, which loses precision on 64-bit systems. 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.
SEE ALSO
This section describes other CPAN modules available that have some feature overlap with
this one. Also see the "REFERENCES" section. Please let me know if any of this
information is inaccurate. Also note that just because a module doesn't match what I
believe are the best set of features, doesn't mean it isn't perfect for someone else.
I will use SoE to indicate the Sieve of Eratosthenes, and MPU to denote this module
(Math::Prime::Util). Some quick alternatives I can recommend if you don't want to use
MPU:
• Math::Prime::FastSieve is the alternative module I use for basic functionality with
small integers. It's fast and simple, and has a good set of features.
• Math::Primality is the alternative module I use for primality testing on bigints. The
downside is that it can be slow, and the functions other than primality tests are very
slow.
• Math::Pari if you want the kitchen sink and can install it and handle using it. There
are still some functions it doesn't do well (e.g. prime count and nth_prime).
Math::Prime::XS has "is_prime" and "primes" functionality. There is no bigint support.
The "is_prime" function uses well-written trial division, meaning it is very fast for
small numbers, but terribly slow for large 64-bit numbers. MPU is similarly fast with
small numbers, but becomes faster as the size increases. MPXS's prime sieve is an
unoptimized non-segmented SoE which returns an array. Sieve bases larger than "10^7"
start taking inordinately long and using a lot of memory (gigabytes beyond "10^10"). E.g.
"primes(10**9, 10**9+1000)" takes 36 seconds with MPXS, but only 0.00015 seconds with MPU.
Math::Prime::FastSieve supports "primes", "is_prime", "next_prime", "prev_prime",
"prime_count", and "nth_prime". The caveat is that all functions only work within the
sieved range, so are limited to about "10^10". It uses a fast SoE to generate the main
sieve. The sieve is 2-3x slower than the base sieve for MPU, and is non-segmented so
cannot be used for larger values. Since the functions work with the sieve, they are very
fast. The fast bit-vector-lookup functionality can be replicated in MPU using
"prime_precalc" but is not required.
Bit::Vector supports the "primes" and "prime_count" functionality in a somewhat similar
way to Math::Prime::FastSieve. It is the slowest of all the XS sieves, and has the most
memory use. It is faster than pure Perl code.
Crypt::Primes supports "random_maurer_prime" functionality. MPU has more options for
random primes (n-digit, n-bit, ranged, and strong) in addition to Maurer's algorithm. MPU
does not have the critical bug RT81858 <https://rt.cpan.org/Ticket/Display.html?id=81858>.
MPU should have a more uniform distribution as well as return a larger subset of primes
(RT81871 <https://rt.cpan.org/Ticket/Display.html?id=81871>). MPU does not depend on
Math::Pari though can run slow for bigints unless the Math::BigInt::GMP or
Math::BigInt::Pari modules are installed. Having Math::Prime::Util::GMP installed also
helps performance for MPU. Crypt::Primes is hardcoded to use Crypt::Random, while MPU
uses Bytes::Random::Secure, and also allows plugging in a random function. This is more
flexible, faster, has fewer dependencies, and uses a CSPRNG for security. MPU can return
a primality certificate. What Crypt::Primes has that MPU does not is the ability to
return a generator.
Math::Factor::XS calculates prime factors and factors, which correspond to the "factor"
and "divisors" functions of MPU. These functions do not support bigints. Both are
implemented with trial division, meaning they are very fast for really small values, but
quickly become unusably slow (factoring 19 digit semiprimes is over 700 times slower).
The function "count_prime_factors" can be done in MPU using "scalar factor($n)". MPU has
no equivalent to "matches", but see the "EXAMPLES" section for a way to produce the
results.
Math::Big version 1.12 includes "primes" functionality. The current code is only usable
for very tiny inputs as it is incredibly slow and uses lots of memory. RT81986
<https://rt.cpan.org/Ticket/Display.html?id=81986> has a patch to make it run much faster
and use much less memory. Since it is in pure Perl it will still run quite slow compared
to MPU.
Math::Big::Factors supports factorization using wheel factorization (smart trial
division). It supports bigints. Unfortunately it is extremely slow on any input that
isn't the product of just small factors. Even 7 digit inputs can take hundreds or
thousands of times longer to factor than MPU or Math::Factor::XS. 19-digit semiprimes
will take hours versus MPU's single milliseconds.
Math::Factoring is a placeholder module for bigint factoring. Version 0.02 only supports
trial division (the Pollard-Rho method does not work).
Math::Prime::TiedArray allows random access to a tied primes array, almost identically to
what MPU provides in Math::Prime::Util::PrimeArray. MPU has attempted to fix
Math::Prime::TiedArray's shift bug (RT58151
<https://rt.cpan.org/Ticket/Display.html?id=58151>). MPU is typically much faster and
will use less memory, but there are some cases where MP:TA is faster (MP:TA stores all
entries up to the largest request, while MPU:PA stores only a window around the last
request).
Math::Primality supports "is_prime", "is_pseudoprime", "is_strong_pseudoprime",
"is_strong_lucas_pseudoprime", "next_prime", "prev_prime", "prime_count", and
"is_aks_prime" functionality. This is a great little module that implements primality
functionality. It was the first CPAN module to support the BPSW test. All inputs are
processed using GMP, so it of course supports bigints. In fact, Math::Primality was made
originally with bigints in mind, while MPU was originally targeted to native integers, but
both have added better support for the other. The main differences are extra
functionality (MPU has more functions) and performance. With native integer inputs, MPU
is generally much faster, especially with "prime_count". For bigints, MPU is slower
unless the Math::Prime::Util::GMP module is installed, in which case MPU is ~2x faster.
Math::Primality also installs a "primes.pl" program, but it has much less functionality
than the one included with MPU.
Math::NumSeq does not have a one-to-one mapping between functions in MPU, but it does
offer a way to get many similar results such as primes, twin primes, Sophie-Germain
primes, lucky primes, moebius, divisor count, factor count, Euler totient, primorials,
etc. Math::NumSeq is set up for accessing these values in order rather than for arbitrary
values, though a few sequences support random access. The primary advantage I see is the
uniform access mechanism for a lot of sequences. For those methods that overlap, MPU is
usually much faster. Importantly, most of the sequences in Math::NumSeq are limited to
32-bit indices.
Math::Pari supports a lot of features, with a great deal of overlap. In general, MPU will
be faster for native 64-bit integers, while it's differs for bigints (Pari will always be
faster if Math::Prime::Util::GMP is not installed; with it, it varies by function). Note
that Pari extends many of these functions to other spaces (Gaussian integers, complex
numbers, vectors, matrices, polynomials, etc.) which are beyond the realm of this module.
Some of the highlights:
"isprime"
The default Math::Pari is built with Pari 2.1.7. This uses 10 M-R tests with randomly
chosen bases (fixed seed, but doesn't reset each invocation like GMP's
"is_probab_prime"). This has a greater chance of false positives compared to the BPSW
test. Calling with "isprime($n,1)" will perform a Pocklington-Lehmer "n-1" proof, but
this becomes unreasonably slow past 70 or so digits.
If Math::Pari is built using Pari 2.3.5 (this requires manual configuration) then the
primality tests are completely different. Using "ispseudoprime" will perform a BPSW
test and is quite a bit faster than the older test. "isprime" now does an APR-CL
proof (fast, but no certificate).
Math::Primality uses a strong BPSW test, which is the standard BPSW test based on the
1980 paper. It has no known counterexamples (though like all these tests, we know
some exist). Pari 2.3.5 (and through at least 2.6.2) uses an almost-extra-strong BPSW
test for its "ispseudoprime" function. This is deterministic for native integers, and
should be excellent for bigints, with a slightly lower chance of counterexamples than
the traditional strong test. Math::Prime::Util uses the full extra-strong BPSW test,
which has an even lower chance of counterexample. With Math::Prime::Util::GMP,
"is_prime" adds 1 to 5 extra M-R tests using random bases, which further reduces the
probability of a composite being allowed to pass.
"primepi"
Only available with version 2.3 of Pari. Similar to MPU's "prime_count" function in
API, but uses a naive counting algorithm with its precalculated primes, so is not of
practical use. Incidently, Pari 2.6 (not usable from Perl) has fixed the pre-
calculation requirement so it is more useful, but is still thousands of times slower
than MPU.
"primes"
Doesn't support ranges, requires bumping up the precalculated primes for larger
numbers, which means knowing in advance the upper limit for primes. Support for
numbers larger than 400M requires using Pari version 2.3.5. If that is used, sieving
is about 2x faster than MPU, but doesn't support segmenting.
"factorint"
Similar to MPU's "factor_exp" though with a slightly different return. MPU offers
"factor" for a linear array of prime factors where
n = p1 * p2 * p3 * ... as (p1,p2,p3,...) and "factor_exp" for an array of
factor/exponent pairs where:
n = p1^e1 * p2^e2 * ... as ([p1,e1],[p2,e2],...) Pari/GP returns an array similar
to the latter. Math::Pari returns a transposed matrix like:
n = p1^e1 * p2^e2 * ... as ([p1,p2,...],[e1,e2,...]) Slower than MPU for all
64-bit inputs on an x86_64 platform, it may be faster for large values on other
platforms. With the newer Math::Prime::Util::GMP releases, bigint factoring is
slightly faster on average in MPU.
"divisors"
Similar to MPU's "divisors".
"forprime", "forcomposite", "fordiv", "sumdiv"
Similar to MPU's "forprimes", "forcomposites", "fordivisors", and "divisor_sum".
"eulerphi", "moebius"
Similar to MPU's "euler_phi" and "moebius". MPU is 2-20x faster for native integers.
MPU also supported range inputs, which can be much more efficient. Without
Math::Prime::Util::GMP installed, MPU is very slow with bigints. With it installed,
it is about 2x slower than Math::Pari.
"gcd", "lcm", "kronecker", "znorder", "znprimroot"
Similar to MPU's "gcd", "lcm", "kronecker", "znorder", and "znprimroot". Pari's
"znprimroot" only returns the smallest root for prime powers. The behavior is
undefined when the group is not cyclic (sometimes it throws an exception, sometimes it
returns an incorrect answer). MPU's "znprimroot" will always return the smallest root
if it exists, and "undef" otherwise.
"sigma"
Similar to MPU's "divisor_sum". MPU is ~10x faster for native integers and about 2x
slower for bigints.
"numbpart"
Similar to MPU's "partitions". This function is not in Pari 2.1, which is the default
version used by Math::Pari. With Pari 2.3 or newer, the functions produce identical
results, but Pari is much, much faster.
"eint1"
Similar to MPU's "ExponentialIntegral".
"zeta"
MPU has "RiemannZeta" which takes non-negative real inputs, while Pari's function
supports negative and complex inputs.
Overall, Math::Pari supports a huge variety of functionality and has a sophisticated and
mature code base behind it (noting that the default version of Pari used is about 10 years
old now). For native integers often using Math::Pari will be slower, but bigints are
often superior and it rarely has any performance surprises. Some of the unique features
MPU offers include super fast prime counts, nth_prime, ECPP primality proofs with
certificates, approximations and limits for both, random primes, fast Mertens
calculations, Chebyshev theta and psi functions, and the logarithmic integral and Riemann
R functions. All with fairly minimal installation requirements.
PERFORMANCE
First, for those looking for the state of the art non-Perl solutions:
Primality testing
For general numbers smaller than 2000 or so digits, I believe MPU is the fastest
solution (it is faster than Pari 2.6.2 and PFGW), though FLINT might be a little
faster for native sizes. For large inputs, PFGW
<http://sourceforge.net/projects/openpfgw/> is the fastest primality testing software
I'm aware of. It 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 <http://www.ellipsa.eu/> 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 <http://sourceforge.net/projects/mpzaprcl/>, the APRCL
from the modern Pari <http://pari.math.u-bordeaux.fr/> package, or the standalone ECPP
from this module with large polynomial set.
Factoring
yafu <http://sourceforge.net/projects/yafu/>, msieve
<http://sourceforge.net/projects/msieve/>, and gmp-ecm <http://ecm.gforge.inria.fr/>
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 <http://code.google.com/p/primesieve/> and yafu
<http://sourceforge.net/projects/yafu/> are the fastest publically available code I am
aware of. Primesieve will additionally take advantage of multiple cores with
excellent efficiency. Tomas Oliveira e Silva's private code may be faster for very
large values, but 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 on par with the SoE in this module. The SoE's in
Pari, yafu, and primesieve are all faster.
Prime Counts and Nth Prime
Outside of private research implementations doing prime counts for "n > 2^64", this
module should be close to state of the art in performance, and supports results up to
"2^64". Further performance improvements are planned, as well as expansion to larger
values.
The fastest solution for small inputs is a hybrid table/sieve method. This module
does this for values below 60M. As the inputs get larger, either the tables have to
grow exponentially or speed must be sacrificed. Hence this is not a good general
solution for most uses.
PRIME COUNTS
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
>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
Small inputs: is_prime from 1 to 20M
2.6s Math::Prime::Util (sieve lookup if prime_precalc used)
3.4s Math::Prime::FastSieve (sieve lookup)
4.4s Math::Prime::Util (trial + deterministic M-R)
10.9s Math::Prime::XS (trial)
36.5s Math::Pari w/2.3.5 (BPSW)
78.2s Math::Pari (10 random M-R)
501.3s Math::Primality (deterministic M-R)
Large native inputs: is_prime from 10^16 to 10^16 + 20M
7.0s Math::Prime::Util (BPSW)
42.6s Math::Pari w/2.3.5 (BPSW)
144.3s Math::Pari (10 random M-R)
664.0s Math::Primality (BPSW)
30 HRS Math::Prime::XS (trial)
These inputs are too large for Math::Prime::FastSieve.
bigints: is_prime from 10^100 to 10^100 + 0.2M
2.5s Math::Prime::Util (BPSW + 1 random M-R)
3.0s Math::Pari w/2.3.5 (BPSW)
12.9s Math::Primality (BPSW)
35.3s Math::Pari (10 random M-R)
53.5s Math::Prime::Util w/o GMP (BPSW)
94.4s Math::Prime::Util (n-1 or ECPP proof)
102.7s Math::Pari w/2.3.5 (APR-CL proof)
• MPU is consistently the fastest solution, and performs the most stringent probable
prime tests on bigints.
• Math::Primality has a lot of overhead that makes it quite slow for native size
integers. With bigints we finally see it work well.
• Math::Pari build with 2.3.5 not only has a better primality test, but runs faster. It
still has quite a bit of overhead with native size integers. Pari/gp 2.5.0's takes
11.3s, 16.9s, and 2.9s respectively for the tests above. MPU is still faster, but
clearly the time for native integers is dominated by the calling overhead.
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
<http://math.boisestate.edu/~liljanab/BOISECRYPTFall09/Jacobsen.pdf> 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 <http://sourceforge.net/projects/yafu/>, msieve
<http://sourceforge.net/projects/msieve/>, gmp-ecm <http://ecm.gforge.inria.fr/>, GGNFS
<http://sourceforge.net/projects/ggnfs/>, and Pari <http://pari.math.u-bordeaux.fr/>. 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 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 <http://www.ellipsa.eu/>.
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
<http://www.issihosts.com/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.
Tomas 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 LMO implementation is based on the 2003 preprint from Christian Bau, as well as the
2006 paper from Tomas Oliveira e Silva. I also want to thank Kim Walisch for the many
discussions about prime counting.
REFERENCES
• Henri Cohen, "A Course in Computational Algebraic Number Theory", Springer, 1996.
Practical computational number theory from the team lead of Pari <http://pari.math.u-
bordeaux.fr/>. Lots of explicit algorithms.
• Hans Riesel, "Prime Numbers and Computer Methods for Factorization", Birkh?user, 2nd
edition, 1994. Lots of information, some code, easy to follow.
• 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 X(s) complex zeros on X(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 <http://oeis.org/wiki/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 Moebius values. <http://walter.lioen.com/papers/LL94.pdf>
• Marc Deleglise and Jooel Rivat, "Computing the summation of the Moebius function",
Experimental Mathematics, v5, n4, pp 291-295, 1996. Enhances the Moebius 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
Moebius 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
Moebius values (not as fast as Deleglise 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
Copyright 2011-2014 by Dana Jacobsen <dana@acm.org>
This program is free software; you can redistribute it and/or modify it under the same
terms as Perl itself.