Math::Prime::Util  Utilities related to prime numbers, including fast sieves
and factoring
Version 0.73
# Nothing is exported by default. List the functions, or use :all.
use Math::Prime::Util ':all'; # import all functions
# The ':rand' tag replaces srand and rand (not done by default)
use Math::Prime::Util ':rand'; # import srand, rand, irand, irand64
# Get a big array reference of many primes
my $aref = primes( 100_000_000 );
# All the primes between 5k and 10k inclusive
$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 nonbigints, is_prime and is_prob_prime will always be 0 or 2.
# They return 0 (composite), 2 (prime), or 1 (probably prime)
my $n = 1000003; # for example
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 MillerRabin
say "$n is a prime or 2/7/61psp" if is_strong_pseudoprime($n, 2, 7, 61);
# Standard and strong LucasSelfridge, 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 undef if given input 2 or less)
$n = prev_prime($n);
# Return Pi(n)  the number of primes E<lt>= n.
my $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
my @prime_factors = factor( $n );
# Return ([p1,e1],[p2,e2], ...) for $n = p1^e1 * p2*e2 * ...
my @pe = factor_exp( $n );
# Get all divisors other than 1 and n
my @divisors = divisors( $n );
# Or just apply a block for each one
my $sum = 0; 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);
# Show all prime partitions of 25
forpart { say "@_" unless scalar grep { !is_prime($_) } @_ } 25;
# List all 3way combinations of an array
my @cdata = qw/apple bread curry donut eagle/;
forcomb { say "@cdata[@_]" } @cdata, 3;
# or all permutations
forperm { say "@cdata[@_]" } @cdata;
# divisor sum
my $sigma = divisor_sum( $n ); # sum of divisors
my $sigma0 = divisor_sum( $n, 0 ); # count of divisors
my $sigmak = divisor_sum( $n, $k );
my $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($rand_prime);
$rand_prime = random_prime(1000); # random prime <= limit
$rand_prime = random_prime(100, 10000); # random prime within a range
$rand_prime = random_ndigit_prime(6); # random 6digit prime
$rand_prime = random_nbit_prime(128); # random 128bit prime
$rand_prime = random_strong_prime(256); # random 256bit strong prime
$rand_prime = random_maurer_prime(256); # random 256bit provable prime
$rand_prime = random_shawe_taylor_prime(256); # as above
A module for number theory in Perl. This includes prime sieving, primality
tests, primality proofs, integer factoring, counts / bounds / approximations
for primes, nth primes, and twin primes, random prime generation, and much
more.
This module is the fastest on CPAN for almost all operations it supports. This
includes 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 1020
digit range, it is often orders of magnitude faster. Typically it is faster
than Math::Pari for 64bit operations.
All operations support both Perl UV's (32bit or 64bit) and bignums. If you
want high performance with big numbers (larger than Perl's native 32bit or
64bit 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 threadsafe and allows concurrency between Perl threads while
still sharing a prime cache. It is not itself multithreaded. See the
Limitations section if you are using Win32 and threads in your program. Also
note that Math::Pari is not threadsafe (and will crash as soon as it is
loaded in threads), so if you use Math::BigInt::Pari rather than
Math::BigInt::GMP or the default backend, things will go pearshaped.
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.
There are two environment variables that affect operation. These are typically
used for validation of the different methods or to simulate systems that have
different support.
If set to 1 then everything is run in pure Perl. No C functions are loaded or
used, as XSLoader is not even called. All toplevel XS functions are replaced
by a pure Perl layer (the PPFE.pm module that supplies a "Pure Perl Front
End").
Caveat: This does not change whether the GMP backend is used. For as much pure
Perl as possible, you will need to set both variables.
If this variable is not set or set to anything other than 1, the module operates
normally.
If set to 1 then the Math::Prime::Util::GMP backend is not loaded, and operation
will be exactly as if it was not installed.
If this variable is not set or set to anything other than 1, the module operates
normally.
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 a recent version of Math::Prime::Util::GMP, as that will vastly
increase the speed of many of the functions. This does require the GMP
<http://gmplib.org> library be installed on your system, but this
increasingly comes preinstalled 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.
 •
 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.
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 extrastrong
BPSW test <http://en.wikipedia.org/wiki/BailliePSW_primality_test> is
used. See the "PRIMALITY TESTING NOTES" section for more discussion.
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
extrastrong 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
MillerRabin bases or BPSW). For larger inputs both "is_prob_prime"
and "is_prime" return probable prime results using the extrastrong
BailliePSW 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 MR tests). This can be done using a
different test (e.g. "is_frobenius_underwood_pseudoprime") or using
additional MR 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 or "random_shawe_taylor_prime" in
Math::Prime::Util which construct random provable 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.
$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 32bit Perl; larger than
"18,446,744,073,709,551,557" in 64bit).
$n = prev_prime($n);
Returns the prime preceding the input number (i.e. the largest prime that is
strictly less than the input). "undef" is returned if the input is 2
or lower.
The behavior in various programs of the
previous prime function is
varied. Pari/GP and Math::Pari returns the input if it is prime, as does
"nearest_le" in Math::Prime::FastSieve. When given an input such
that the return value will be the first prime less than 2,
Math::Prime::FastSieve, Math::Pari, Pari/GP, and older versions of MPU will
return 0. Math::Primality and the current MPU will return "undef".
WolframAlpha returns "2". Maple gives a range error.
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 { 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, OEIS A002808
<http://oeis.org/A002808>, are the numbers greater than 1 which are not
prime: "4, 6, 8, 9, 10, 12, 14, 15, ...".
Similar to "forcomposites", but skipping all even numbers. The odd
composites, OEIS A071904 <http://oeis.org/A071904>, are the numbers
greater than 1 which are not prime and not divisible by two: "9, 15, 21,
25, 27, 33, 35, ...".
Similar to "forcomposites", but only giving composites with exactly
two factors. The semiprimes, OEIS A001358 <http://oeis.org/A001358>, are
the products of two primes: "4, 6, 9, 10, 14, 15, 21, 22, 25, ...".
This is essentially equivalent to:
forcomposites { if (is_semiprime($_)) { ... } }
forfactored { say "$_: @_"; } 100;
Given a block and either an end number or start/end pair, calls the block for
each number in the inclusive range. $_ is set to the number while @_ holds the
factors. Especially for small inputs or large ranges, This can be faster than
calling "factor" on each sequential value.
Similar to the arrays returned by similar functions such as "forpart",
the values in @_ are readonly. Any attempt to modify them will result in
undefined behavior.
This corresponds to the Pari/GP 2.10 "forfactored" function.
Similar to "forfactored", but skipping numbers in the range that have
a repeated factor. Inside the block, the moebius function can be cheaply
computed as "((scalar(@_) & 1) ? 1 : 1)" or similar.
This corresponds to the Pari/GP 2.10 "forsquarefree" function.
fordivisors { $prod *= $_ } $n;
Given a block and a nonnegative number "n", the block is called with
$_ set to each divisor in sorted order. Also see "divisor_sum".
forpart { say "@_" } 25; # unrestricted partitions
forpart { say "@_" } 25,{n=>5} # ... with exactly 5 values
forpart { say "@_" } 25,{nmax=>5} # ... with <=5 values
Given a nonnegative number "n", the block is called with @_ set to
the array of additive integer partitions. The operation is very similar to the
"forpart" function in Pari/GP 2.6.x, though the ordering is
different. The ordering is lexicographic. Use "partitions" to get
just the count of unrestricted partitions.
An optional hash reference may be given to produce restricted partitions. Each
value must be a nonnegative integer. The allowable keys are:
n restrict to exactly this many values
amin all elements must be at least this value
amax all elements must be at most this value
nmin the array must have at least this many values
nmax the array must have at most this many values
prime all elements must be prime (nonzero) or nonprime (zero)
Like forcomb and forperm, the partition return values are readonly. Any attempt
to modify them will result in undefined behavior.
Similar to "forpart", but iterates over integer compositions rather
than partitions. This can be thought of as all ordering of partitions, or
alternately partitions may be viewed as an ordered subset of compositions. The
ordering is lexicographic. All options from "forpart" may be used.
The number of unrestricted compositions of "n" is "2^(n1)".
Given nonnegative arguments "n" and "k", the block is
called with @_ set to the "k" element array of values from 0 to
"n1" representing the combinations in lexicographical order. While
the binomial function gives the total number, this function can be used to
enumerate the choices.
Rather than give a data array as input, an integer is used for "n". A
convenient way to map to array elements is:
forcomb { say "@data[@_]" } @data, 3;
where the block maps the combination array @_ to array values, the argument for
"n" is given the array since it will be evaluated as a scalar and
hence give the size, and the argument for "k" is the desired size of
the combinations.
Like forpart and forperm, the index return values are readonly. Any attempt to
modify them will result in undefined behavior.
If the second argument "k" is not supplied, then all ksubsets are
returned starting with the smallest set "k=0" and continuing to
"k=n". Each ksubset is in lexicographical order. This is the power
set of "n".
This corresponds to the Pari/GP 2.10 "forsubset" function.
Given nonnegative argument "n", the block is called with @_ set to
the "k" element array of values from 0 to "n1"
representing permutations in lexicographical order. The total number of calls
will be "n!".
Rather than give a data array as input, an integer is used for "n". A
convenient way to map to array elements is:
forperm { say "@data[@_]" } @data;
where the block maps the permutation array @_ to array values, and the argument
for "n" is given the array since it will be evaluated as a scalar
and hence give the size.
Like forpart and forcomb, the index return values are readonly. Any attempt to
modify them will result in undefined behavior.
Similar to forperm, but iterates over derangements. This is the set of
permutations skipping any which maps an element to its original position.
# Show all anagrams of 'serpent':
formultiperm { say join("",@_) } [split(//,"serpent")];
Similar to "forperm" but takes an array reference as an argument. This
is treated as a multiset, and the block will be called with each multiset
permutation. While the standard permutation iterator takes a scalar and
returns index permutations, this takes the set itself.
If all values are unique, then the results will be the same as a standard
permutation. Otherwise, the results will be similar to a standard permutation
removing duplicate entries. While generating all permutations and filtering
out duplicates works, it is very slow for large sets. This iterator will be
much more efficient.
There is no ordering requirement for the input array reference. The results will
be in lexicographic order.
forsetproduct { say "@_" } [1,2,3],[qw/a b c/],[qw/@ $ !/];
Takes zero or more array references as arguments and iterates over the set
product (i.e. Cartesian product or cross product) of the lists. The given
subroutine is repeatedly called with @_ set to the current list. Since no
deduplication is done, this is not literally a "set" product.
While zero or one array references are valid, the result is not very
interesting. If any array reference is empty, the product is empty, so no
subroutine calls are performed.
The subroutine is given an array whose values are aliased to the inputs, and are
not set to readonly. Hence modifying the array inside the subroutine
will cause sideeffects.
As with other iterators, the "lastfor" function will cause an early
exit.
forprimes { lastfor,return if $_ > 1000; $sum += $_; } 1e9;
Calling lastfor requests that the current for... loop stop after this call.
Ideally this would act exactly like a "last" inside a loop, but
technical reasons mean it does not exit the block early, hence one typically
adds a "return" if needed.
my $it = prime_iterator;
$sum += $it>() for 1..100000;
Returns a closurestyle 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).
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 bidirectional), and also handles iterating
across bigints.
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) + (ba))",
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.
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". For larger inputs various methods are used including Dusart
(2010), Büthe (2014,2015), and Axler (2014). 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 can be used
for very large values.
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.
A slightly faster but much less accurate answer can be obtained by averaging the
upper and lower bounds.
Returns the lesser of twin primes between the lower and upper limits
(inclusive), with a lower limit of 2 if none is given. This is OEIS A001359
<http://oeis.org/A001359>. Given a twin prime pair "(p,q)"
with "q = p + 2", "p prime", and <q prime>, this
function uses "p" to represent the pair. Hence the bounds need to
include "p", and the returned list will have "p" but not
"q".
This works just like the "primes" function, though only the first
primes of twin prime pairs are returned. Like that function, an array
reference is returned.
Similar to prime count, but returns the count of twin primes (primes
"p" where "p+2" is also prime). Takes either a single
number indicating a count from 2 to the argument, or two numbers indicating a
range.
The primes being counted are the first value, so a range of "(3,5)"
will return a count of two, because both 3 and 5 are counted as twin primes. A
range of "(12,13)" will return a count of zero, because neither
"12+2" nor "13+2" are prime. In contrast,
"primesieve" requires all elements of a constellation to be within
the range to be counted, so would return one for the first example (5 is not
counted because its pair 7 is not in the range).
There is no useful formula known for this, unlike prime counts. We sieve for the
answer, using some small table acceleration.
Returns an approximation to the twin prime count of "n". This returns
quickly and has a very small error for large values. The method used is
conjecture B of Hardy and Littlewood 1922, as stated in Sebah and Gourdon
2002. For inputs under 10M, a correction factor is additionally applied to
reduce the mean squared error.
Returns an array reference to semiprimes between the lower and upper limits
(inclusive), with a lower limit of 4 if none is given. This is OEIS A001358
<http://oeis.org/A001358>. The semiprimes are composite integers which
are products of exactly two primes.
This works just like the "primes" function. Like that function, an
array reference is returned.
Similar to prime count, but returns the count of semiprimes (composites with
exactly two factors). Takes either a single number indicating a count from 2
to the argument, or two numbers indicating a range.
A fast method that requires computation only to the square root of the range end
is used, unless the range is so small that walking it is faster.
Returns an approximation to the semiprime count of "n". This returns
quickly and is typically square root accurate.
Returns the Ramanujan primes R_n between the upper and lower limits (inclusive),
with a lower limit of 2 if none is given. This is OEIS A104272
<http://oeis.org/A104272>. These are the Rn such that if "x >
Rn" then "prime_count"(n)  "prime_count"(n/2) >=
"n".
This has a similar API to the "primes" and "twin_primes"
functions, and like them, returns an array reference.
Generating Ramanujan primes takes some effort, including overhead to cover a
range. This will be substantially slower than generating standard primes.
Similar to prime count, but returns the count of Ramanujan primes. Takes either
a single number indicating a count from 2 to the argument, or two numbers
indicating a range.
While not nearly as efficient as prime_count, this does use a number of speedups
that result it in being much more efficient than generating all the Ramanujan
primes.
A fast approximation of the count of Ramanujan primes under "n".
A fast lower limit on the count of Ramanujan primes under "n".
A fast upper limit on the count of Ramanujan primes under "n".
my @candidates = sieve_range(2**1000, 10000, 40000);
Given a start value "n", and native unsigned integers
"width" and "depth", a sieve of maximum depth
"depth" is done for the "width" consecutive numbers
beginning with "n". An array of offsets from the start is returned.
The returned list contains those offsets in the range "n" to
"n+width1" where "n + offset" has no prime factors less
than "depth".
my @s = sieve_prime_cluster(1, 1e9, 2,6,8,12,18,20);
Efficiently finds prime clusters between the first two arguments "low"
and "high". The remaining arguments describe the cluster. The
cluster values must be even, less than 31 bits, and strictly increasing. Given
a cluster set "C", the returned values are all primes in the range
where "p+c" is prime for each "c" in the cluster set
"C". For returned values under "2^64", all cluster values
are definitely prime. Above this range, all cluster values are BPSW probable
primes (no counterexamples known).
This function returns an array rather than an array reference. Typically the
number of returned values is much lower than for other primes functions, so
this uses the more convenient array return. This function has an identical
signature to the function of the same name in Math::Prime::Util:GMP.
The cluster is described as offsets from 0, with the implicit prime at 0. Hence
an empty list is asking for all primes (the cluster "p+0"). A list
with the single value 2 will find all twin primes (the cluster where
"p+0" and "p+2" are prime). The list "2,6,8"
will find prime quadruplets. Note that there is no requirement that the list
denote a constellation (a cluster with minimal distance)  the list
"42,92,606" is just fine.
Returns the summation of primes between the lower and upper limits (inclusive),
with a lower limit of 2 if none is given. This is essentially similar to
either of:
$sum = 0; forprimes { $sum += $_ } $low,$high; $sum;
# or
vecsum( @{ primes($low,$high) } );
but is much more efficient.
The current implementation is a smalltableenhanced sieve count for sums that
fit in a UV, an efficient sieve count for small ranges, and a Legendre sum
method for larger values.
While this is fairly efficient, the state of the art is Kim Walisch's primesum
<https://github.com/kimwalisch/primesum>. It is recommended for very
large values, as it can be hundreds of times faster.
print_primes(1_000_000); # print the first 1 million primes
print_primes(1000, 2000); # print primes in range
print_primes(2,1000,fileno(STDERR)) # print to a different descriptor
With a single argument this prints all primes from 2 to "n" to
standard out. With two arguments it prints primes between "low" and
"high" to standard output. With three arguments it prints primes
between "low" and "high" to the file descriptor given. If
the file descriptor cannot be written to, this will croak with
"print_primes write error". It will produce identical output to:
forprimes { say } $low,$high;
The point of this function is just efficiency. It is over 10x faster than using
"say", "print", or "printf", though much more
limited in functionality. A later version may allow a file handle as the third
argument.
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".
Like most programs with similar functionality, this is onebased. nth_prime(0)
returns "undef", nth_prime(1) returns 2.
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 lowbiased 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 30 seconds.
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 (32bit or 64bit), the
result will take a very long time. A later version of Math::Prime::Util::GMP
may include this functionality which would help for 32bit machines.
my $lower_limit = nth_prime_lower($n);
my $upper_limit = nth_prime_upper($n);
# For all $n: $lower_limit <= nth_prime($n) <= $upper_limit
Returns an analytical upper or lower bound on the Nth prime. No sieving is done,
so these are fast even for large inputs.
For tiny values of "n". exact answers are returned. For small inputs,
an inverse of the opposite prime count bound is used. For larger values, the
Dusart (2010) and Axler (2013) bounds are used.
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. For values where the nth prime is smaller than
"2^64", the inverse Riemann R function is used. For larger values,
the inverse logarithmic integral is used.
The value returned will not necessarily be prime. This applies to all the
following nth prime approximations, where the returned value is close to the
real value, but no effort is made to coerce the result to the nearest set
element.
Returns the Nth twin prime. This is done via sieving and counting, so is not
very fast for large values.
Returns an approximation to the Nth twin prime. A curve fit is used for small
inputs (under 1200), while for larger inputs a binary search is done on the
approximate twin prime count.
Returns the Nth semiprime, similar to where a "forsemiprimes" loop
would end after "N" iterations, but much more efficiently.
Returns an approximation to the Nth semiprime. Curve fitting is used to get a
fairly close approximation that is orders of magnitude better than the simple
"n log n / log log n" approximation for large "n".
Returns the Nth Ramanujan prime. For reasonable size values of "n",
e.g. under "10^8" or so, this is relatively efficient for single
calls. If multiple calls are being made, it will be much more efficient to get
the list once.
A fast approximation of the Nth Ramanujan prime.
A fast lower limit on the Nth Ramanujan prime.
A fast upper limit on the Nth Ramanujan prime.
Takes a positive number "n" and one or more nonzero positive bases as
input. Returns 1 if the input is a probable prime to each base, 0 if not. This
is the simple Fermat primality test. Removing primes, given base 2 this
produces the sequence OEIS A001567 <http://oeis.org/A001567>.
For practical use, "is_strong_pseudoprime" is a much stronger test
with similar or better performance.
Note that there is a set of composites (the Carmichael numbers) that will pass
this test for all bases. This downside is not shared by the Euler and strong
probable prime tests (also called the SolovayStrassen and MillerRabin
tests).
Takes a positive number "n" and one or more nonzero positive bases as
input. Returns 1 if the input is an Euler probable prime to each base, 0 if
not. This is the Euler test, sometimes called the EulerJacobi test. Removing
primes, given base 2 this produces the sequence OEIS A047713
<http://oeis.org/A047713>.
If 0 is returned, then the number really is a composite. If 1 is returned, then
it is either a prime or an Euler pseudoprime to all the given bases. Given
enough distinct bases, the chances become very high that the number is
actually prime.
This test forms the basis of the SolovayStrassen test, which is a precursor to
the MillerRabin test (which uses the strong probable prime test). There are
no analogies to the Carmichael numbers for this test. For the Euler test, at
most 1/2 of witnesses pass for a composite, while at most 1/4 pass for
the strong pseudoprime test.
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 "n" and one or more nonzero positive bases as
input. Returns 1 if the input is a strong probable prime to each base, 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 high 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
32bit 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 non2 even input makes the function match most other
implementations including Math::Primality's "is_strong_pseudoprime"
function.
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 LucasSelfridge pseudoprime). Removing primes, this
produces the sequence OEIS A217120 <http://oeis.org/A217120>.
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 LucasSelfridge pseudoprime). This is one half of
the BPSW primality test (the MillerRabin strong pseudoprime test with base 2
being the other half). Removing primes, this produces the sequence OEIS
A217255 <http://oeis.org/A217255>.
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/mathscinetgetitem?mr=1680879>). This test has more
stringent conditions than the strong Lucas test, and produces about 60% fewer
pseudoprimes. Performance is typically 2030%
faster than the strong
Lucas test.
The parameters are selected using the BaillieOEIS 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>.
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 extrastrong 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.ubordeaux.fr/faq.html#primetest>.
Because the "U = 0" condition is ignored, this produces about 5% more
pseudoprimes than the extrastrong Lucas test. However this is still only 66%
of the number produced by the strong LucasSelfridge test. No BPSW
counterexamples have been found with any of the Lucas tests described.
Takes a positive number "n" as input and returns 1 if "n"
passes Colin Plumb's Euler Criterion primality test. Pseudoprimes to this test
are a subset of the base 2 Fermat and Euler tests, but a superset of the base
2 strong pseudoprime (MillerRabin) test.
The main reason for this test is that is a bit more efficient than other
probable prime tests.
Takes a positive number "n" as input and returns 1 if "n"
divides P(n) where P(n) is the Perrin number of "n". The Perrin
sequence is defined by "P(n) = P(n2) + P(n3)" with "P(0) = 3,
P(1) = 0, P(2) = 2".
While pseudoprimes are relatively rare (the first two are 271441 and 904631),
infinitely many exist. They have significant overlap with the base2
pseudoprimes and strong pseudoprimes, making the test inferior to the Lucas or
Frobenius tests for combined testing. The pseudoprime sequence is OEIS A013998
<http://oeis.org/A013998>.
The implementation uses modular prefilters, Montgomery math, and the
Adams/Shanks doubling method. This is significantly more efficient than other
known implementations.
An optional second argument "r" indicates whether to run additional
tests. With "r=1", "P(n) = 1 mod n" is also verified,
creating the "minimal restricted" test. With "r=2", the
full signature is also tested using the Adams and Shanks (1982) rules (without
the quadratic form test). With "r=3", the full signature is testing
using the Grantham (2000) test, which additionally does not allow pseudoprimes
to be divisible by 2 or 23. The minimal restricted pseudoprime sequence is
OEIS A018187 <http://oeis.org/A018187>.
Takes a positive number "n" as input and returns 1 if
"1^((n1/2)) C_((n1/2)" is congruent to 2 mod "n", where
"C_n" is the nth Catalan number. The nth Catalan number is equal to
"binomial(2n,n)/(n+1)". All odd primes satisfy this condition, and
only three known composites.
The pseudoprime sequence is OEIS A163209 <http://oeis.org/A163209>.
There is no known efficient method to perform the Catalan primality test, so it
is a curiosity rather than a practical test. The implementation uses a method
from Charles Greathouse IV (2015) and results from Aebi and Cairns (2008) to
produce results many orders of magnitude faster than other known
implementations, but it is still vastly slower than other compositeness tests.
Takes a positive number "n" as input, and two optional parameters
"a" and "b", and returns 1 if the "n" is a
Frobenius probable prime with respect to the polynomial "x^2  ax +
b". Without the parameters, "b = 2" and "a" is the
least positive odd number such that "(a^24bn) = 1". This
selection has no pseudoprimes below "2^64" and none known. In any
case, the discriminant "a^24b" must not be a perfect square.
Some authors use the Fibonacci polynomial "x^2x1" corresponding to
"(1,1)" as the default method for a Frobenius probable prime test.
This creates a weaker test than most other parameter choices (e.g. over twenty
times more pseudoprimes than "(3,5)"), so is not used as the
default here. With the "(1,1)" parameters the pseudoprime sequence
is OEIS A212424 <http://oeis.org/A212424>.
The Frobenius test is a stronger test than the Lucas test. Any Frobenius
"(a,b)" pseudoprime is also a Lucas "(a,b)" pseudoprime
but the converse is not true, as any Frobenius "(a,b)" pseudoprime
is also a Fermat pseudoprime to the base "b". We can see that with
the default parameters this is similar to, but somewhat weaker than, the BPSW
test used by this module (which uses the strong and extrastrong versions of
the probable prime and Lucas tests respectively).
Also see the more efficient "is_frobenius_khashin_pseudoprime" and
"is_frobenius_underwood_pseudoprime" which have no known
counterexamples and run quite a bit faster.
Takes a positive number as input, and returns 1 if the input passes the
efficient Frobenius test of Paul Underwood. This selects a parameter
"a" as the least nonnegative integer such that
"(a^24n)=1", then verifies that "(x+2)^(n+1) = 2a + 5 mod
(x^2ax+1,n)". This combines a Fermat and Lucas test with a cost of only
slightly more than 2 strong pseudoprime tests. This makes it similar to, but
faster than, a regular Frobenius test.
There are no known pseudoprimes to this test and extensive computation has shown
no counterexamples under "2^50". This test also has no overlap with
the BPSW test, making it a very effective method for adding additional
certainty. Performance at 1e12 is about 60% slower than BPSW.
Takes a positive number as input, and returns 1 if the input passes the
Frobenius test of Sergey Khashin. This ensures "n" is not a perfect
square, selects the parameter "c" as the smallest odd prime such
that "(cn)=1", then verifies that "(1+D)^n = (1D) mod
n" where "D = sqrt(c) mod n".
There are no known pseudoprimes to this test and Khashin (2018) shows there are
no counterexamples under "2^64". Performance at 1e12 is about 40%
slower than BPSW.
Takes a positive number ("n") as input and a positive number
("k") of bases to use. Performs "k" MillerRabin tests
using uniform random bases between 2 and "n2".
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" MillerRabin
tests. This function can be used if one wants more assurances for nonproven
primes, such as for cryptographic uses where the size is large enough that
proven primes are not desired.
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 64bit input (native or bignum), this uses either a deterministic set of
MillerRabin tests (1, 2, or 3 tests) or a strong BPSW test consisting of a
single base2 strong probable prime test followed by a strong Lucas test. This
has been verified with Jan Feitsma's 2PSP database to produce no false
results for 64bit inputs. Hence the result will always be 0 (composite) or 2
(prime).
For inputs larger than "2^64", an extrastrong BailliePSW 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.
Given a positive number input, returns 0 (composite), 2 (definitely prime), or 1
(probably prime), using the BPSW primality test (extrastrong 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 pretests to find easy composites. If you know this is not
necessary, then calling "is_bpsw_prime" may save a small amount of
time.
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
MillerRabin bases or the BPSW test. Both of these have been tested for all
small (64bit) 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 PocklingtonLehmer test. This
requires "n1" 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
"n1" proof before using ECPP. Certificates are available with
either method. This results in proofs of 200digit 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.
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 multiline 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.
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.
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 selfcontained. 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 MR).
 "BLS3"

Type BLS3
N 2297612322987260054928384863
Q 16501461106821092981
A 5
A simple n1 style proof using BLS75 theorem 3. This block verifies if:
a Q is odd
b Q > 2
c Q divides N1
. Let M = (N1)/Q
d MQ+1 = N
e M > 0
f 2Q+1 > sqrt(N)
g A^((N1)/2) mod N = N1
h A^(M/2) mod N != N1
 "Pocklington"

Type Pocklington
N 2297612322987260054928384863
Q 16501461106821092981
A 5
A simple n1 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 N1
. Let M = (N1)/Q
b M > 0
c M < Q
d MQ+1 = N
e A > 1
f A^(N1) 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 MQ1 = N
e M > 0
f 2Q1 > 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 n1 proof using BLS theorem 5. This requires N1 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] < N1
c3 A[i] > 1
c4 A[i] < N
c5 Q[i] divides N1
. Let F = N1 divided by each Q[i] as many times as evenly possible
. Let R = (N1)/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 + (r1)*F + 1)
f n < P
g s = 0 OR r^28s is not a perfect square
h For each i (0 .. maxi):
h1 A[i]^(N1) mod N = 1
h2 gcd(A[i]^((N1)/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
AtkinGoldwasserKilianMorain 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
say "$n is definitely prime" if is_aks_prime($n);
Takes a nonnegative number as input, and returns 1 if the input passes the
AgrawalKayalSaxena (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 nontheoretical 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.
This implementation uses theorem 4.1 from Bernstein (2003). It runs
substantially faster than the original, v6 revised paper with Lenstra
improvements, or the late 2002 improvements of Voloch and Bornemann. The GMP
implementation uses a binary segmentation method for modular polynomial
multiplication (see Bernstein's 2007 Quartic paper), which reduces to a single
scalar multiplication, at which GMP excels. Because of this, the GMP
implementation is likely to be faster once the input is larger than
"2^33".
say "2^6071 (M607) is a Mersenne prime" if is_mersenne_prime(607);
Takes a nonnegative number "p" as input and returns 1 if the Mersenne
number "2^p1" is prime. Since an enormous effort has gone into
testing these, a list of known Mersenne primes is used to accelerate this.
Beyond the highest sequential Mersenne prime (currently 37,156,667) this
performs pretesting followed by the LucasLehmer test.
The LucasLehmer test is a deterministic unconditional test that runs very fast
compared to other primality methods for numbers of comparable size, and vastly
faster than any known generalform primality proof methods. While this test is
fast, the GMP implementation is not nearly as fast as specialized programs
such as "prime95". Additionally, since we use the table for
"small" numbers, testing via this function call will only occur for
numbers with over 9.8 million digits. At this size, tools such as
"prime95" are greatly preferred.
Takes a positive number "n" as input and returns back either 0 or 1,
indicating whether "n" is a Ramanujan prime. Numbers that can be
produced by the functions "ramanujan_primes" and
"nth_ramanujan_prime" will return 1, while all other numbers will
return 0.
There is no simple function for this predicate, so Ramanujan primes through at
least "n" are generated, then a search is performed for
"n". This is not efficient for multiple calls.
say "$n is a perfect square" if is_power($n, 2);
say "$n is a perfect cube" if is_power($n, 3);
say "$n is a ", is_power($n), "th power";
Given a single nonnegative integer input "n", returns k if "n =
r^k" for some integer "r > 1, k > 1", and 0 otherwise.
The k returned is the largest possible. This can be used in a boolean
statement to determine if "n" is a perfect power.
If given two arguments "n" and "k", returns 1 if
"n" is a "kth" power, and 0 otherwise. For example, if
"k=2" then this detects perfect squares. Setting "k=0"
gives behavior like the first case (the largest root is found and its value is
returned).
If a third argument is present, it must be a scalar reference. If "n"
is a kth power, then this will be set to the kth root of "n". For
example:
my $n = 222657534574035968;
if (my $pow = is_power($n, 0, \my $root)) { say "$n = $root^$pow" }
# prints: 222657534574035968 = 2948^5
This corresponds to Pari/GP's "ispower" function with integer
arguments.
Given an integer input "n", returns "k" if "n =
p^k" for some prime p, and zero otherwise.
If a second argument is present, it must be a scalar reference. If the return
value is nonzero, then it will be set to "p".
This corresponds to Pari/GP's "isprimepower" function.
Given a positive integer "n", returns 1 if "n" is a perfect
square, 0 otherwise. This is identical to "is_power(n,2)".
This corresponds to Pari/GP's "issquare" function.
Given a nonnegative integer input "n", returns the integer square
root. For native integers, this is equal to "int(sqrt(n))".
This corresponds to Pari/GP's "sqrtint" function.
Given an nonnegative integer "n" and positive exponent "k",
return the integer kth root of "n". This is the largest integer
"r" such that "r^k <= n".
If a third argument is present, it must be a scalar reference. It will be set to
"r^k".
Technically if "n" is negative and "k" is odd, the root
exists and is equal to "sign(n) * rootint(abs(n),k)". It was
decided to follow the behavior of Pari/GP and Math::BigInt and disallow
negative "n".
This corresponds to Pari/GP's "sqrtnint" function.
say "decimal digits: ", 1+logint($n, 10);
say "digits in base 12: ", 1+logint($n, 12);
my $be; my $e = logint(1000,2, \$be);
say "smallest power of 2 less than 1000: 2^$e = $be";
Given a nonzero positive integer "n" and an integer base
"b" greater than 1, returns the largest integer "e" such
that "b^e <= n".
If a third argument is present, it must be a scalar reference. It will be set to
"b^e".
This corresponds to Pari/GP's "logint" function.
say "Fibonacci($_) = ", lucasu(1,1,$_) for 0..100;
Given integers "P", "Q", and the nonnegative integer
"k", computes "U_k" for the Lucas sequence defined by
"P","Q". These include the Fibonacci numbers
("1,1"), the Pell numbers ("2,1"), the Jacobsthal
numbers ("1,2"), the Mersenne numbers ("3,2"), and more.
This corresponds to OpenPFGW's "lucasU" function and gmpy2's
"lucasu" function.
say "Lucas($_) = ", lucasv(1,1,$_) for 0..100;
Given integers "P", "Q", and the nonnegative integer
"k", computes "V_k" for the Lucas sequence defined by
"P","Q". These include the Lucas numbers
("1,1").
This corresponds to OpenPFGW's "lucasV" function and gmpy2's
"lucasv" function.
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: " P < n" ; " Q <
n" ; " k >= 0" ; " n >= 2".
Given a list of integers, returns the greatest common divisor. This is often
used to test for coprimality <https://oeis.org/wiki/Coprimality>.
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
Given two integers "x" and "y", returns "u,v,d"
such that "d = gcd(x,y)" and "u*x + v*y = d". This uses
the extended Euclidian algorithm to compute the values satisfying
Bézout's Identity.
This corresponds to Pari's "gcdext" function, which was renamed from
"bezout" out Pari 2.6. The results will hence match
"bezout" in Math::Pari.
say chinese( [14,643], [254,419], [87,733] ); # 87041638
Solves a system of simultaneous congruences using the Chinese Remainder Theorem
(with extension to noncoprime moduli). A list of "[a,n]" pairs are
taken as input, each representing an equation "x ≡ a mod n".
If no solution exists, "undef" is returned. If a solution is
returned, the modulus is equal to the lcm of all the given moduli (see
"lcm". In the standard case where all values of "n" are
coprime, this is just the product. The "n" values must be positive
integers, while the "a" values are integers.
Comparison to similar functions in other software:
Math::ModInt::ChineseRemainder:
cr_combine( mod(a1,m1), mod(a2,m2), ... )
Pari/GP:
chinese( [Mod(a1,m1), Mod(a2,m2), ...] )
Mathematica:
ChineseRemainder[{a1, a2, ...}{m1, m2, ...}]
say "Totient sum 500,000: ", vecsum(euler_phi(0,500_000));
Returns the sum of all arguments, each of which must be an integer. This is
similar to List::Util's "sum0" in List::Util function, but has a
very important difference. List::Util turns all inputs into doubles and
returns a double, which will mean incorrect results with large integers.
"vecsum" sums (signed) integers and returns the untruncated result.
Processing is done on native integers while possible.
say "Totient product 5,000: ", vecprod(euler_phi(1,5_000));
Returns the product of all arguments, each of which must be an integer. This is
similar to List::Util's "product" in List::Util function, but keeps
all results as integers and automatically switches to bigints if needed.
say "Smallest Totient 100k200k: ", vecmin(euler_phi(100_000,200_000));
Returns the minimum of all arguments, each of which must be an integer. This is
similar to List::Util's "min" in List::Util function, but has a very
important difference. List::Util turns all inputs into doubles and returns a
double, which gives incorrect results with large integers. "vecmin"
validates and compares all results as integers. The validation step will make
it a little slower than "min" in List::Util but this prevents
accidental and unintentional use of floats.
say "Largest Totient 100k200k: ", vecmax(euler_phi(100_000,200_000));
Returns the maximum of all arguments, each of which must be an integer. This is
similar to List::Util's "max" in List::Util function, but has a very
important difference. List::Util turns all inputs into doubles and returns a
double, which gives incorrect results with large integers. "vecmax"
validates and compares all results as integers. The validation step will make
it a little slower than "max" in List::Util but this prevents
accidental and unintentional use of floats.
say "Count of nonzero elements: ", vecreduce { $a + !!$b } (0,@v);
my $checksum = vecreduce { $a ^ $b } @{twin_primes(1000000)};
Does a reduce operation via left fold. Takes a block and a list as arguments.
The block uses the special local variables "a" and "b"
representing the accumulation and next element respectively, with the result
of the block being used for the new accumulation. No initial element is used,
so "undef" will be returned with an empty list.
The interface is exactly the same as "reduce" in List::Util. This was
done to increase portability and minimize confusion. See chapter 7 of Higher
Order Perl (or many other references) for a discussion of reduce with empty or
singularelement lists. It is often a good idea to give an identity element as
the first list argument.
While operations like vecmin, vecmax, vecsum, vecprod, etc. can be fairly easily
done with this function, it will not be as efficient. There are a wide variety
of other functions that can be easily made with reduce, making it a useful
tool.
say "all values are Carmichael" if vecall { is_carmichael($_) } @n;
Short circuit evaluations of a block over a list. Takes a block and a list as
arguments. The block is called with $_ set to each list element, and
evaluation on list elements is done until either all list values have been
evaluated or the result condition can be determined. For instance, in the
example of "vecall" above, evaluation stops as soon as any value
returns false.
The interface is exactly the same as the "any", "all",
"none", "notall", and "first" functions in
List::Util. This was done to increase portability and minimize confusion.
Unlike other vector functions like "vecmax", "vecmax",
"vecsum", etc. there is no added value to using these versus the
ones from List::Util. They are here for convenience.
These operations can fairly easily be mapped to "scalar(grep {...}
@n)", but that does not shortcircuit and is less obvious.
say "first Carmichael is index ", vecfirstidx { is_carmichael($_) } @n;
Returns the index of the first element in a list that evaluates to true. Just
like vecfirst, but returns the index instead of the value. Returns 1 if the
item could not be found.
This interface matches "firstidx" and "first_index" from
List::MoreUtils.
say "Power set: ", join(" ",vecextract(\@v,$_)) for 0..2**scalar(@v)1;
@word = vecextract(["a".."z"], [15, 17, 8, 12, 4]);
Extracts elements from an array reference based on a mask, with the result
returned as an array. The mask is either an unsigned integer which is treated
as a bit mask, or an array reference containing integer indices.
If the second argument is an integer, each bit set in the mask results in the
corresponding element from the array reference to be returned. Bits are read
from the right, so a mask of 1 returns the first element, while 5 will return
the first and third. The mask may be a bigint.
If the second argument is an array reference, then its elements will be used as
zerobased indices into the first array. Duplicate values are allowed and the
ordering is preserved. Hence these are equivalent:
vecextract($aref, $iref);
@$aref[@$iref];
say "product of digits of n: ", vecprod(todigits($n));
Given an integer "n", return an array of digits of "n". An
optional second integer argument specifies a base (default 10). For example,
given a base of 2, this returns an array of binary digits of "n". An
optional third argument specifies a length for the returned array. The result
will be either have upper digits truncated or have leading zeros added. This
is most often used with base 2, 8, or 16.
The values returned may be readonly. todigits(0) returns an empty array. The
base must be at least 2, and is limited to an int. Length must be at least
zero and is limited to an int.
This corresponds to Pari's "digits" and "binary" functions,
and Mathematica's "IntegerDigits" function.
say "decimal 456 in hex is ", todigitstring(456, 16);
say "last 4 bits of $n are: ", todigitstring($n, 2, 4);
Similar to "todigits" but returns a string. For bases <= 10, this
is equivalent to joining the array returned by "todigits". For bases
between 11 and 36, lower case characters "a" to "z" are
used to represent larger values. This makes "todigitstring($n,16)"
return a usable hex string.
This corresponds to Mathematica's "IntegerString" function.
say "hex 1c8 in decimal is ", fromdigits("1c8", 16);
say "Base 3 array to number is: ", fromdigits([0,1,2,2,2,1,0],3);
This takes either a string or array reference, and an optional base (default
10). With a string, each character will be interpreted as a digit in the given
base, with both upper and lower case denoting values 11 through 36. With an
array reference, the values indicate the entries in that location, and values
larger than the base are allowed (results are carried). The result is a number
(either a native integer or a bigint).
This corresponds to Pari's "fromdigits" function and Mathematica's
"FromDigits" function.
# Sum digits of primes to 1 million.
my $s=0; forprimes { $s += sumdigits($_); } 1e6; say $s;
Given an input "n", return the sum of the digits of "n". Any
nondigit characters of "n" are ignored (including negative signs
and decimal points). This is similar to the command
"vecsum(split(//,$n))" but faster, allows nonpositiveinteger
inputs, and can sum in other bases.
An optional second argument indicates the base of the input number. This
defaults to 10, and must be between 2 and 36. Any character that is outside
the range 0 to "base1" will be ignored.
If no base is given and the input number "n" begins with
"0x" or "0b" then it will be interpreted as a string in
base 16 or 2 respectively.
Regardless of the base, the output sum is a decimal number.
This is similar but not identical to Pari's "sumdigits" function from
version 2.8 and later. The Pari/GP function always takes the input as a
decimal number, uses the optional base as a base to first convert to, then
sums the digits. This can be done with either "vecsum(todigits($n,
$base))" or "sumdigits(todigitstring($n,$base))".
say "The inverse of 42 mod 2017 = ", invmod(42,2017);
Given two integers "a" and "n", return the inverse of
"a" modulo "n". If not defined, undef is returned. If
defined, then the return value multiplied by "a" equals 1 modulo
"n".
The results correspond to the Pari result of "lift(Mod(1/a,n))". The
semantics with respect to negative arguments match Pari. Notably, a negative
"n" is negated, which is different from Math::BigInt, but in both
cases the return value is still congruent to 1 modulo "n" as
expected.
Given two integers "a" and "n", return the square root of
"a" mod "n". If no square root exists, undef is returned.
If defined, the return value "r" will always satisfy "r^2 = a
mod n".
If the modulus is prime, the function will always return "r", the
smaller of the two square roots (the other being "r mod p". If the
modulus is composite, one of possibly many square roots will be returned, and
it will not necessarily be the smallest.
Given three integers "a", "b", and "n" where
"n" is positive, return "(a+b) mod n". This is
particularly useful when dealing with numbers that are larger than a halfword
but still native size. No bigint package is needed and this can be 10200x
faster than using one.
Given three integers "a", "b", and "n" where
"n" is positive, return "(a*b) mod n". This is
particularly useful when "n" fits in a native integer. No bigint
package is needed and this can be 10200x faster than using one.
Given three integers "a", "b", and "n" where
"n" is positive, return "(a ** b) mod n". Typically binary
exponentiation is used, so the process is very efficient. With native size
inputs, no bigint library is needed.
Given three integers "a", "b", and "n" where
"n" is positive, return "(a/b) mod n". This is done as
"(a * (1/b mod n)) mod n". If no inverse of "b" mod
"n" exists then undef if returned.
say "$n is divisible by 2 ", valuation($n,2), " times.";
Given integers "n" and "k", returns the numbers of times
"n" is divisible by "k". This is a very limited version of
the algebraic valuation meaning, just applied to integers. This corresponds to
Pari's "valuation" function. 0 is returned if "n" or
"k" is one of the values "1", 0, or 1.
Given an integer "n", returns the binary Hamming weight of abs(n).
This is also called the population count, and is the number of 1s in the
binary representation. This corresponds to Pari's "hammingweight"
function for "t_INT" arguments.
say "$n has no repeating factors" if is_square_free($n);
Returns 1 if the input "n" has no repeated factor.
for (1..1e6) { say if is_carmichael($_) } # Carmichaels under 1,000,000
Returns 1 if the input "n" is a Carmichael number. These are
composites that satisfy "b^(n1) ≡ 1 mod n" for all "1
< b < n" relatively prime to "n". Alternately Korselt's
theorem says these are composites such that "n" is squarefree and
"p1" divides "n1" for all prime divisors "p"
of "n".
For inputs larger than 50 digits after removing very small factors, this uses a
probabilistic test since factoring the number could take unreasonably long.
The first 150 primes are used for testing. Any that divide "n" are
checked for squarefreeness and the Korselt condition, while those that do
not divide "n" are used as the pseudoprime base. The chances of a
nonCarmichael passing this test are less than "2^150".
This is the OEIS series A002997 <http://oeis.org/A002997>.
Returns 0 if the input "n" is not a quasiCarmichael number, and the
number of bases otherwise. These are squarefree composites that satisfy
"p+b" divides "n+b" for all prime factors "p" or
"n" and for one or more nonzero integer "b".
This is the OEIS series A257750 <http://oeis.org/A257750>.
Given a positive integer "n", returns 1 if "n" is a
semiprime, 0 otherwise. A semiprime is the product of exactly two primes.
The boolean result is the same as "scalar(factor(n)) == 2", but this
function performs shortcuts that can greatly speed up the operation.
Given an integer "d", returns 1 if "d" is a fundamental
discriminant, 0 otherwise. We consider 1 to be a fundamental discriminant.
This is the OEIS series A003658 <http://oeis.org/A003658> (positive) and
OEIS series A003657 <http://oeis.org/A003657> (negative).
This corresponds to Pari's "isfundamental" function.
Given an integer "n", returns 1 if there exists an integer
"x" where "euler_phi(x) == n".
This corresponds to Pari's "istotient" function, though without the
optional second argument to return an "x". Math::NumSeq::Totient
also has a similar function.
Also see "inverse_totient" which gives the count or list of values
that produce a given totient. This function is more efficient than getting the
full count or list.
Given a positive integer "n", if there exists a "v" where
"v! % n == n1" and "n % v != 1", then "v" is
returned. Otherwise 0.
For n prime, this is the OEIS series A063980 <http://oeis.org/A063980>.
Given integers "x" and "s", return 1 if x is an sgonal
number, 0 otherwise. "s" must be greater than 2.
If a third argument is present, it must be a scalar reference. It will be set to
n if x is the nth sgonal number. If the function returns 0, then it will be
unchanged.
This corresponds to Pari's "ispolygonal" function.
say "$n is square free" if moebius($n) != 0;
$sum += moebius($_) for (1..200); say "Mertens(200) = $sum";
say "Mertens(2000) = ", vecsum(moebius(0,2000));
Returns μ(n), the Möbius 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 squarefree (i.e. "n"
has a repeated factor), and "1^t" if "n" is a product of
"t" distinct primes. This is an important function in prime number
theory. Like SAGE, we define "moebius(0) = 0" for convenience.
If called with two arguments, they define a range "low" to
"high", and the function returns an array with the value of the
Möbius function for every n from low to high inclusive. Large values of
high will result in a lot of memory use. The algorithm used for ranges is
Deléglise 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 readonly 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).
say "Mertens(10M) = ", mertens(10_000_000); # = 1037
Returns M(n), the Mertens function for a nonnegative integer input. This
function is defined as "sum(moebius(1..n))", but calculated more
efficiently for large inputs. For example, computing Mertens(100M) takes:
time approx mem
0.4s 0.1MB mertens(100_000_000)
3.0s 880MB vecsum(moebius(1,100_000_000))
56s 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, even though they are shared constants, is
not good for memory at this size. In comparison, this function will generate
the equivalent output via a sieving method that is relatively memory frugal
and very fast. The current method is a simple "n^1/2" version of
Deléglise and Rivat (1996), which involves calculating all moebius
values to "n^1/2", which in turn will require prime sieving to
"n^1/4".
Various algorithms exist for this, using differing quantities of μ(n).
The simplest way is to efficiently sum all "n" values. Benito and
Varona (2008) show a clever and simple method that only requires
"n/3" values. Deléglise and Rivat (1996) describe a segmented
method using only "n^1/3" values. The current implementation does a
simple nonsegmented "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.
say "The Euler totient of $n is ", euler_phi($n);
Returns φ(n), the Euler totient function (also called Euler's phi or phi
function) for an integer value. This is an arithmetic function 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", Pari pre2.6.2 raises an exception, and Pari 2.6.2 and
newer returns 2.
If called with two arguments, they define a range "low" to
"high", and the function returns a list with the totient of every n
from low to high inclusive.
In array context, given a positive integer "n", returns the complete
list of values "x" where "euler_phi(x) = n". This can be a
memory intensive operation if there are many values.
In scalar context, returns just the count of values. This is faster and uses
substantially less memory. The list/scalar distinction is similar to
"factor" and "divisors".
This roughly corresponds to the Maple function "InverseTotient", and
the hidden Mathematica function "EulerPhiInverse". The algorithm
used is from Max Alekseyev (2016).
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
ktuples 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
Dedekind psi function, where "psi(n) = J(2,n) / J(1,n)".
Returns Ramanujan's sum of the two positive variables "k" and
"n". This is the sum of the nth powers of the primitive kth roots
of unity.
say "exp(lambda($_)) = ", exp_mangoldt($_) for 1 .. 100;
Returns EXP(Λ(n)), the exponential of the Mangoldt function (also known
as von Mangoldt's function) for an integer value. 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.
Returns λ(n), the Liouville function for a nonnegative integer input.
This is 1 raised to Ω(n) (the total number of prime factors).
say chebyshev_theta(10000);
Returns θ(n), the first Chebyshev function for a nonnegative integer
input. This is the sum of the logarithm of each prime where "p <=
n". This is effectively:
my $s = 0; forprimes { $s += log($_) } $n; return $s;
say chebyshev_psi(10000);
Returns ψ(n), the second Chebyshev function for a nonnegative integer
input. This is the sum of the logarithm of each prime power where "p^k
<= n" for an integer k. An alternate but slower computation is as the
summatory Mangoldt function, such as:
my $s = 0; for (1..$n) { $s += log(exp_mangoldt($_)) } return $s;
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 "kth" powers of
the divisors to be returned.
This is known as the sigma function (see Hardy and Wright section 16.7). The API
is identical to Pari/GP's "sigma" function, and not dissimilar to
Mathematica's "DivisorSigma[k,n]" function. This function is useful
for calculating things like aliquot sums, abundant numbers, perfect numbers,
etc.
With various "k" values, the results are the OEIS sequences OEIS
series A000005 <http://oeis.org/A000005> ("k=0", number of
divisors), OEIS series A000203 <http://oeis.org/A000203>
("k=1", sum of divisors), OEIS series A001157
<http://oeis.org/A001157> ("k=2", sum of squares of divisors),
OEIS series A001158 <http://oeis.org/A001158> ("k=4", sum of
cubes of divisors), 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.
Takes a positive integer as input and returns the value of Ramanujan's tau
function. The result is a signed integer. This corresponds to Pari v2.8's
"tauramanujan" function and Mathematica's "RamanujanTau"
function.
This currently uses a simple method based on divisor sums, which does not have a
good computational growth rate. Pari's implementation uses Hurwitz class
numbers and is more efficient for large inputs.
$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.
$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.
$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 memoryfriendly
results than using a factorial.
Calculates the partition function p(n) for a nonnegative 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
multiprecision floating point. In 10 seconds:
70 Integer::Partition
90 MPU forpart { $n++ }
10_000 MPU pure Perl partitions
250_000 MPU GMP partitions
35_000_000 Pari's numbpart
500_000_000 Jonathan Bober's partitions_c.cc v0.6
If you want the enumerated partitions, see "forpart".
Returns the Carmichael function (also called the reduced totient function, or
Carmichael λ(n)) of a positive integer argument. It is the smallest
positive integer "m" such that "a^m = 1 mod n" for every
integer "a" coprime to "n". This is OEIS series A002322
<http://oeis.org/A002322>.
Returns the Kronecker symbol "(an)" for two integers. The possible
return values with their meanings for odd prime "n" are:
0 a = 0 mod n
1 a is a quadratic residue mod n (a = x^2 mod n for some x)
1 a is a quadratic nonresidue mod n (no a where a = x^2 mod 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, Mathematica's
"KroneckerSymbol[n,m]" function, and GMP's
"mpz_kronecker(a,n)", "mpz_jacobi(a,n)", and
"mpz_legendre(a,n)" functions.
Given positive integer argument "n", returns the factorial of
"n", defined as the product of the integers 1 to "n" with
the special case of "factorial(0) = 1". This corresponds to Pari's
factorial(n) and Mathematica's "Factorial[n]" functions.
Given two positive integer arguments "n" and "m", returns
"n! mod m". This is much faster than computing the large
factorial(n) followed by a mod operation.
While very efficient, this is not state of the art. Currently, Fredrik
Johansson's fast multipoint polynomial evaluation method as used in FLINT is
the fastest known method. This becomes noticeable for "n" >
"10^8" or so, and the O(n^.5) versus O(n) complexity makes it quite
extreme as the input gets larger.
Given integer arguments "n" and "k", returns the binomial
coefficient "n*(n1)*...*(nk+1)/k!", also known as the choose
function. Negative arguments use the Kronenburg extensions
<http://arxiv.org/abs/1105.3689/>. This corresponds to Pari's
"binomial(n,k)" function, Mathematica's "Binomial[n,k]"
function, and GMP's "mpz_bin_ui" function.
For negative arguments, this matches Mathematica. Pari does not implement the
"n < 0, k <= n" extension and instead returns 0 for this case.
GMP's API does not allow negative "k" but otherwise matches.
Math::BigInt does not implement any extensions and the results for "n
< 0, k " 0> are undefined.
Returns 12 times the HurwitzKronecker class number of the input integer
"n". This will always be an integer due to the premultiplication by
12. The result is 0 for any input less than zero or congruent to 1 or 2 mod 4.
This is related to Pari's qfbhclassno(n) where hclassno(n) for positive
"n" equals "12 * qfbhclassno(n)" in Pari/GP. This is OEIS
A259825 <http://oeis.org/A259825>.
Returns the Bernoulli number "B_n" for an integer argument
"n", as a rational number represented by two Math::BigInt objects.
B_1 = 1/2. This corresponds to Pari's bernfrac(n) and Mathematica's
"BernoulliB" functions.
Having a modern version of Math::Prime::Util::GMP installed will make a big
difference in speed. That module uses a fast Pi/Zeta method. Our pure Perl
backend uses the Seidel method as shown by Peter Luschny. This is faster than
Math::Pari which uses an older algorithm, but quite a bit slower than modern
Pari, Mathematica, or our GMP backend.
This corresponds to Pari's "bernfrac" function and Mathematica's
"BernoulliB" function.
Returns the Bernoulli number "B_n" for an integer argument
"n", as a Math::BigFloat object using the default precision. An
optional second argument may be given specifying the precision to be used.
This corresponds to Pari's "bernreal" function.
say "s(14,2) = ", stirling(14, 2);
say "S(14,2) = ", stirling(14, 2, 2);
say "L(14,2) = ", stirling(14, 2, 3);
Returns the Stirling numbers of either the first kind (default), the second
kind, or the third kind (the unsigned Lah numbers), with the kind selected as
an optional third argument. It takes two nonnegative integer arguments
"n" and "k" plus the optional "type". This
corresponds to Pari's "stirling(n,k,{type})" function and
Mathematica's "StirlingS1" / "StirlingS2" functions.
Stirling numbers of the first kind are "1^(nk)" times the number of
permutations of "n" symbols with exactly "k" cycles.
Stirling numbers of the second kind are the number of ways to partition a set
of "n" elements into "k" nonempty subsets. The Lah
numbers are the number of ways to split a set of "n" elements into
"k" nonempty lists.
Returns the Harmonic number "H_n" for an integer argument
"n", as a rational number represented by two Math::BigInt objects.
The harmonic numbers are the sum of reciprocals of the first "n"
natural numbers: "1 + 1/2 + 1/3 + ... + 1/n".
Binary splitting (Fredrik Johansson's elegant formulation) is used.
This corresponds to Mathematica's "HarmonicNumber" function.
Returns the Harmonic number "H_n" for an integer argument
"n", as a Math::BigFloat object using the default precision. An
optional second argument may be given specifying the precision to be used.
For large "n" values, using a lower precision may result in faster
computation as an asymptotic formula may be used. For precisions of 13 or
less, native floating point is used for even more speed.
$order = znorder(2, next_prime(10**16)6);
Given two positive integers "a" and "n", returns the
multiplicative order of "a" modulo "n". This is the
smallest positive integer "k" such that "a^k ≡ 1 mod
n". Returns 1 if "a = 1". 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[a,n]" function.
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.
Given two nonnegative numbers "a" and "n", returns 1 if
"a" is a primitive root modulo "n", and 0 if not. If
"a" is a primitive root, then euler_phi(n) is the smallest
"e" for which "a^e = 1 mod n".
$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 for native integers first applies SilverPohligHellman on
the group order to possibly reduce the problem to a set of smaller problems.
The solutions are then performed using a mixture of trial, Shanks' BSGS, and
Pollard's DLP Rho.
The PP implementation is less sophisticated, with only a memoryheavy BSGS being
used.
$phi = legendre_phi(1000000000, 41);
Given a nonnegative integer "n" and a nonnegative 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.
$approx_prime_count = inverse_li(1000000000);
Given a nonnegative integer "n", returns the least integer value
"k" such that Li(k) >= n>. Since the logarithmic integral
Li(n) is a good approximation to the number of primes less than "n",
this function is a good simple approximation to the nth prime.
@p = numtoperm(10,654321); # @p=(1,8,2,7,6,5,3,4,9,0)
Given a nonnegative integer "n" and integer "k", return the
rank "k" lexicographic permutation of "n" elements.
"k" will be interpreted as mod "n!".
This will match iteration number "k" (zero based) of
"forperm".
This corresponds to Pari's "numtoperm(n,k)" function, though Pari uses
an implementation specific ordering rather than lexicographic.
$k = permtonum([1,8,2,7,6,5,3,4,9,0]); # $k = 654321
Given an array reference containing integers from 0 to "n", returns
the lexicographic permutation rank of the set. This is the inverse of the
"numtoperm" function. All integers up to "n" must be
present.
This will match iteration number "k" (zero based) of
"forperm". The result will be between 0 and "n!1".
This corresponds to Pari's permtonum(n) function, though Pari uses an
implementation specific ordering rather than lexicographic.
@p = randperm(100); # returns shuffled 0..99
@p = randperm(100,4) # returns 4 elements from shuffled 0..99
@s = @data[randperm(1+$#data)]; # shuffle an array
@p = @data[randperm(1+$#data,2)]; # pick 2 from an array
With a single argument "n", this returns a random permutation of the
values from 0 to "n1".
When given a second argument "k", the returned list will have only
"k" elements. This is more efficient than truncating the full
shuffled list.
The randomness comes from our CSPRNG.
@shuffled = shuffle(@data);
Takes a list as input, and returns a random permutation of the list. Like
randperm, the randomness comes from our CSPRNG.
This function is functionally identical to the "shuffle" function in
List::Util. The only difference is the random source (Chacha20 with better
randomness, a larger period, and a larger state). This does make it slower.
If the entire shuffled array is desired, this is faster than slicing with
"randperm" as shown in its example above. If, however, a
"pick" operation is desired, e.g. pick 2 random elements from a
large array, then the slice technique can be hundreds of times faster.
Prior to version 5.20, Perl's "rand" function used the system rand
function. This meant it varied by system, and was almost always a poor choice.
For 5.20, Perl standardized on "drand48" and includes the source so
there are no system dependencies. While this was an improvement,
"drand48" is not a good PRNG. It really only has 32 bits of random
values, and fails many statistical tests. See
<http://www.pcgrandom.org/statisticaltests.html> for more information.
There are much better choices for standard random number generators, such as the
Mersenne Twister, PCG, or Xoroshiro128+. Someday perhaps Perl will get one of
these to replace drand48. In the mean time, Math::Random::MTwist provides
numerous features and excellent performance, or this module.
Since we often deal with random primes for cryptographic purposes, we have
additional requirements. This module uses a CSPRNG for its random stream. In
particular, ChaCha20, which is the same algorithm used by BSD's
"arc4random" and "/dev/urandom" on BSD and Linux 4.8+.
Seeding is performed at startup using the Win32 Crypto API (on Windows),
"/dev/urandom", "/dev/random", or Crypt::PRNG, whichever
is found first.
We use the original ChaCha definition rather than RFC7539. This means a 64bit
counter, resulting in a period of 2^72 bytes or 2^68 calls to drand or
<irand64>. This compares favorably to the 2^48 period of Perl's
"drand48". It has a 512bit state which is significantly larger than
the 48bit "drand48" state. When seeding, 320 bits (40 bytes) are
used. Among other things, this means all 52! permutations of a shuffled card
deck are possible, which is not true of "shuffle" in List::Util.
One might think that performance would suffer from using a CSPRNG, but
benchmarking shows it is less than one might expect. does not seem to be the
case. The speed of irand, irand64, and drand are within 20% of the fastest
existing modules using nonCSPRNG methods, and 2 to 20 times faster than most.
While a faster underlying RNG is useful, the Perl call interface overhead is a
majority of the time for these calls. Carefully tuning that interface is
critical.
For performance on large amounts of data, see the tables in
"random_bytes".
Each thread uses its own context, meaning seeding in one thread has no impact on
other threads. In addition to improved security, this is better for
performance than a single context with locks. If explicit control of multiple
independent streams are needed then using a more specific module is
recommended. I believe Crypt::PRNG (part of CryptX) and Bytes::Random::Secure
are good alternatives.
Using the ":rand" export option will define "rand" and
"srand" as similar but improved versions of the system functions of
the same name, as well as "irand" and "irand64".
$n32 = irand; # random 32bit integer
Returns a random 32bit integer using the CSPRNG.
$n64 = irand64; # random 64bit integer
Returns a random 64bit integer using the CSPRNG (on 64bit Perl).
$f = drand; # random floating point value in [0,1)
$r = drand(25.33); # random floating point value in [0,25.33)
Returns a random NV (Perl's native floating point) using the CSPRNG. The API is
similar to Perl's "rand" but giving better results.
The number of bits returned is equal to the number of significand bits of the NV
type used in the Perl build. By default Perl uses doubles and the returned
values have 53 bits (even on 32bit Perl). If Perl is built with long double
or quadmath support, each value may have 64 or even 113 bits. On newer Perls,
one can check the Config variable "nvmantbits" to see how many are
filled.
This gives
substantially better quality random numbers than the default
Perl "rand" function. Among other things, on modern Perl's,
"rand" uses drand48, which has 32 bits of notverygood randomness
and 16 more bits of obvious patterns (e.g. the 48th bit alternates, the 47th
has a period of 4, etc.). Output from "rand" fails at least 5 tests
from the TestU01 SmallCrush suite, while our function easily passes.
With the ":rand" tag, this function is additionally exported as
"rand".
$str = random_bytes(32); # 32 random bytes
Given an unsigned number "n" of bytes, returns a string filled with
random data from the CSPRNG. Performance for large quantities:
Module/Method Rate Type
  
Math::Prime::Util::GMP 1067 MB/s CSPRNG  ISAAC
ntheory random_bytes 384 MB/s CSPRNG  ChaCha20
Crypt::PRNG 140 MB/s CSPRNG  Fortuna
Crypt::OpenSSL::Random 32 MB/s CSPRNG  SHA1 counter
Math::Random::ISAAC::XS 15 MB/s CSPRNG  ISAAC
ntheory entropy_bytes 13 MB/s CSPRNG  /dev/urandom
Crypt::Random 12 MB/s CSPRNG  /dev/urandom
Crypt::Urandom 12 MB/s CSPRNG  /dev/urandom
Bytes::Random::Secure 6 MB/s CSPRNG  ISAAC
ntheory pure perl ISAAC 5 MB/s CSPRNG  ISAAC (no XS)
Math::Random::ISAAC::PP 2.5 MB/s CSPRNG  ISAAC (no XS)
ntheory pure perl ChaCha 1.0 MB/s CSPRNG  ChaCha20 (no XS)
Data::Entropy::Algorithms 0.5 MB/s CSPRNG  AESCTR
Math::Random::MTwist 927 MB/s PRNG  Mersenne Twister
Bytes::Random::XS 109 MB/s PRNG  drand48
pack CORE::rand 25 MB/s PRNG  drand48 (no XS)
Bytes::Random 2.6 MB/s PRNG  drand48 (no XS)
Similar to random_bytes, but directly using the entropy source. This is not
normally recommended as it can consume shared system resources and is
typically slow  on the computer that produced the "random_bytes"
chart above, using "dd" generated the same 13 MB/s performance as
our "entropy_bytes" function.
The actual performance will be highly system dependent.
$n32 = urandomb(32); # Classic irand32, returns a UV
$n = urandomb(1024); # Random integer less than 2^1024
Given a number of bits "b", returns a random unsigned integer less
than "2^b". The result will be uniformly distributed between 0 and
"2^b1" inclusive.
$n = urandomm(100); # random integer in [0,99]
$n = urandomm(1024); # random integer in [0,1023]
Given a positive integer "n", returns a random unsigned integer less
than "n". The results will be uniformly distributed between 0 and
"n1" inclusive. Care is taken to prevent modulo bias.
Takes a binary string "data" as input and seeds the internal CSPRNG.
This is not normally needed as system entropy is used as a seed on startup.
For best security this should be 16128 bytes of good entropy. No more than
1024 bytes will be used.
With no argument, reseeds using system entropy, which is preferred.
If the "secure" configuration has been set, then this will croak if
given an argument. This allows for control of reseeding with entropy the
module gets itself, but not user supplied.
Takes a single UV argument and seeds the CSPRNG with it, as well as returning
it. If no argument is given, a new UV seed is constructed. Note that this
creates a very weak seed from a cryptographic standpoint, so it is useful for
testing or simulations but "csrand" is recommended, or keep using
the system entropy default seed.
The API is nearly identical to the system function "srand". It uses a
UV which can be 64bit rather than always 32bit. The behaviour for
"undef", empty string, empty list, etc. is slightly different (we
treat these as 0).
This function is not exported with the ":all" tag, but is with
":rand".
If the "secure" configuration has been set, this function will croak.
Manual seeding using "srand" is not compatible with cryptographic
security.
An alias for "drand", not exported unless the ":rand" tag is
used.
my($n, $factors) = random_factored_integer(1000000);
Given a positive nonzero input "n", returns a uniform random integer
in the range 1 to "n", along with an array reference containing the
factors.
This uses Kalai's algorithm for generating random integers along with their
factorization, and is much faster than the naive method of generating random
integers followed by a factorization. A later implementation may use Bach's
more efficient algorithm.
my $small_prime = random_prime(1000); # random prime <= limit
my $rand_prime = random_prime(100, 10000); # random prime within a range
Returns a pseudorandomly 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 nonrandom 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 ~16bit threshold but within the native bit size, a Monte Carlo method is
used. 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.
say "My 4digit prime number is: ", random_ndigit_prime(4);
Selects a random ndigit prime, where the input is an integer number of digits.
One of the primes within that range (e.g. 1000  9999 for 4digits) will be
uniformly selected.
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.
my $bigprime = random_nbit_prime(512);
Selects a random nbit prime, where the input is an integer number of bits. A
prime with the nth bit set will be uniformly selected.
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 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.
my $bigprime = random_strong_prime(512);
Constructs an nbit 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.
 •
 p1 has a large prime factor r.
 •
 p+1 has a large prime factor s
 •
 r1 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.
my $bigprime = random_proven_prime(512);
Constructs an nbit random proven prime. Internally this may use
"is_provable_prime"("random_nbit_prime") or
"random_maurer_prime" depending on the platform and bit size.
my($n, $cert) = random_proven_prime_with_cert(512)
Similar to "random_proven_prime", but returns a twoelement array
containing the nbit 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.
my $bigprime = random_maurer_prime(512);
Construct an nbit 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.
The performance with Math::Prime::Util::GMP installed is hundreds of times
faster, so it is highly recommended.
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.
If you don't need absolutely proven results, then it is somewhat faster to use
"random_nbit_prime" either by itself or with some additional tests,
e.g. "miller_rabin_random" and/or
"is_frobenius_underwood_pseudoprime". One could also run
is_provable_prime on the result, but this will be slow.
my($n, $cert) = random_maurer_prime_with_cert(512)
As with "random_maurer_prime", but returns a twoelement array
containing the nbit 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.
my $bigprime = random_shawe_taylor_prime(8192);
Construct an nbit provable prime, using the ShaweTaylor algorithm in section
C.6 of FIPS 1864. This uses 512 bits of randomness and SHA256 as the hash.
This is a slightly simpler and older (1986) method than Maurer's 1999
construction. It is a bit faster than Maurer's method, and uses less system
entropy for large sizes. The primary reason to use this rather than Maurer's
method is to use the FIPS 1864 algorithm.
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. Also see
"random_maurer_prime" and "random_proven_prime".
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.
my($n, $cert) = random_shawe_taylor_prime_with_cert(4096)
As with "random_shawe_taylor_prime", but returns a twoelement array
containing the nbit 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
"Pocklington" types.
Takes a positive integer number of bits "bits", returns a random
semiprime of exactly "bits" bits. The result has exactly two prime
factors (hence semiprime).
The factors will be approximately equal size, which is typical for cryptographic
use. For example, a 64bit semiprime of this type is the product of two 32bit
primes. "bits" must be 4 or greater.
Some effort is taken to select uniformly from the universe of
"bits"bit semiprimes. This takes slightly longer than some methods
that do not select uniformly.
Takes a positive integer number of bits "bits", returns a random
semiprime of exactly "bits" bits. The result has exactly two prime
factors (hence semiprime).
The factors are uniformly selected from the universe of all "bits"bit
semiprimes. This means semiprimes with one factor equal to 2 will be most
common, 3 next most common, etc. "bits" must be 3 or greater.
Some effort is taken to select uniformly from the universe of
"bits"bit semiprimes. This takes slightly longer than some methods
that do not select uniformly.
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;
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.
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.
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:
verbose verbose level. 1 or more will result in extra output.
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)
secure disable ability to manually seed the CSPRNG
prime_set_config( assume_rh => 1 );
Allows setting of some parameters. Currently the only parameters are:
verbose The default setting of 0 will generate no extra output.
Setting to 1 or higher results in extra output. For
example, at setting 1 the AKS algorithm will indicate
the chosen r and s values. At setting 2 it will output
a sequence of dots indicating progress. Similarly, for
random_maurer_prime, setting 3 shows real time progress.
Factoring large numbers is another place where verbose
settings can give progress indications.
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
reenable. 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 reenable.
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.
secure The CSPRNG may no longer be manually seeded. Once set,
this option cannot be disabled. L</srand> will croak
if called, and L</csrand> will croak if called with any
arguments. L</csrand> with no arguments is still allowed,
as that will use system entropy without giving anything
to the caller. The point of this option is to ensure that
any called functions do not try to control the RNG.
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 Ω(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 does a little trial division, a check for perfect powers,
followed by combinations of Pollard's Rho, SQUFOF, and Pollard's p1. The
combination is applied to each nonprime factor found.
Factoring bigints works with pure Perl, and can be very handy on 32bit machines
for numbers just over the 32bit 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 23x faster (though still
100x slower than the real GMP code).
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 ω(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.
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 (see Hardy and Wright section
16.7). This is OEIS A000005 <http://oeis.org/A000005>. The results is
identical to evaluating the array in scalar context, but more efficient. This
corresponds to Pari's "numdiv" and Mathematica's
"DivisorSigma[0,n]" functions.
Also see the "for_divisors" functions for looping over the divisors.
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. Like all the
specificalgorithm *_factor routines, this is not exported unless explicitly
requested.
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).
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. Overall it has the same
advantages and disadvantages as Fermat's method.
my @factors = lehman_factor($n);
Produces factors, not necessarily prime, of the positive number input. An
optional argument, defaulting to 0 (false), indicates whether to run trial
division. Without trial division, is possible the function will be unable to
find a factor, in which case a single element, the input, is returned.
This is Warren D. Smith's Lehman core with minor modifications. It is limited to
42bit inputs: "n < 8796393022208".
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.
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 prefactoring very large inputs, as they are very fast at
finding small factors.
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 "p1" 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 "p1" is
smooth (it has no large factors).
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.
my @factors = ecm_factor($n);
my @factors = ecm_factor($n, 100, 400, 10); # B1, B2, # of curves
Produces factors, not necessarily prime, of the positive number input. This is
the elliptic curve method using two stages.
my $Ei = ExponentialIntegral($x);
Given a nonzero floating point input "x", this returns the
realvalued 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 nonBigInt/BigFloat inputs, the result should be accurate to at least 14
digits.
For BigInt / BigFloat inputs, full accuracy and performance is obtained only if
Math::Prime::Util::GMP is installed. If this module is not available, then
other methods are used and give at least 14 digits of accuracy: continued
fractions ("x < 1"), rational Chebyshev approximation (" 1
< x < 0"), a convergent series (small positive "x"), or
an asymptotic divergent series (large positive "x").
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 nonBigInt/BigFloat objects, the result should be accurate to at least 14
digits.
For BigInt / BigFloat inputs, full accuracy and performance is obtained only if
Math::Prime::Util::GMP is installed.
my $z = RiemannZeta($s);
Given a floating point input "s" where "s >= 0", returns
the floating point value of ζ(s)1, where ζ(s) is the Riemann
zeta function. One is subtracted to ensure maximum precision for large values
of "s". The zeta function is the sum from k=1 to infinity of "1
/ k^s". This function only uses real arguments, so is basically the Euler
Zeta function.
If the bignum module has been loaded, all inputs will be treated as if they were
Math::BigFloat objects.
For nonBigInt/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 inputs, full accuracy and performance is obtained only if
Math::Prime::Util::GMP is installed. If this module is not available, then
other methods are used and give at least 14 digits of accuracy: Either Borwein
(1991) algorithm 2, or the basic series. Math::BigFloat RT 43692
<https://rt.cpan.org/Ticket/Display.html?id=43692> can produce incorrect
highaccuracy computations when GMP is not used.
my $r = RiemannR($x);
Given a positive nonzero 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 nonBigInt/BigFloat objects, the result should be accurate to at least 14
digits.
For BigInt / BigFloat inputs, full accuracy and performance is obtained only if
Math::Prime::Util::GMP is installed. If this module are not available,
accuracy should be 35 digits.
Returns the principal branch of the Lambert W function of a real value. Given a
value "k" this solves for "W" in the equation "k =
We^W". The input must not be less than "1/e". This corresponds
to Pari's "lambertw" function and Mathematica's
"ProductLog" / "LambertW" function.
This function handles all real value inputs with noncomplex return values. This
is a superset of Pari's "lambertw" which is similar but only for
positive arguments. Mathematica's function is much more detailed, with both
branches, complex arguments, and complex results.
Calculation will be done with C long doubles if the input is a standard scalar,
but if bignum is in use or if the input is a BigFloat type, then extended
precision results will be used.
Speed of the native code is about half of the fastest native code (Veberic's
C++), and about 30x faster than Pari/GP. However the bignum calculation is
slower than Pari/GP.
my $tau = 2 * Pi; # $tau = 6.28318530717959
my $tau = 2 * Pi(40); # $tau = 6.283185307179586476925286766559005768394
With no arguments, returns the value of Pi as an NV. With a positive integer
argument, returns the value of Pi with the requested number of digits
(including the leading 3). The return value will be an NV if the number of
digits fits in an NV (typically 15 or less), or a Math::BigFloat object
otherwise.
For sizes over 10k digits, having either Math::Prime::Util::GMP or
Math::BigInt::GMP installed will help performance. For sizes over 50k the one
is highly recommended.
Print Fibonacci numbers:
perl Mntheory=:all E 'say lucasu(1,1,$_) for 0..20'
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 64bit 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); }'
Generate Carmichael numbers (OEIS A002997 <http://oeis.org/A002997>):
perl Mntheory=:all E 'foroddcomposites { say if is_carmichael($_) } 1e6;'
# Less efficient, similar to Mathematica or MAGMA:
perl Mntheory=:all E 'foroddcomposites { say if $_ % carmichael_lambda($_) == 1 } 1e6;'
Examining the η3(x) function of Planat and Solé (2011):
sub nu3 {
my $n = shift;
my $phix = chebyshev_psi($n);
my $nu3 = 0;
foreach my $nu (1..3) {
$nu3 += (moebius($nu)/$nu)*LogarithmicIntegral($phix**(1/$nu));
}
return $nu3;
}
say prime_count(1000000);
say prime_count_approx(1000000);
say nu3(1000000);
Construct and use a SophieGermain 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 32bit 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/sum_primes/;
say sum_primes(2_000_000);
# ... or do it a little more manually ...
use Math::Prime::Util qw/forprimes/;
my $sum = 0;
forprimes { $sum += $_ } 2_000_000;
say $sum;
# ... or do it using a big list ...
use Math::Prime::Util qw/vecsum primes/;
say vecsum( @{primes(2_000_000)} );
Project Euler, problem 21 (Amicable numbers):
use Math::Prime::Util qw/divisor_sum/;
my $sum = 0;
foreach my $x (1..10000) {
my $y = divisor_sum($x)$x;
$sum += $x + $y if $y > $x && $x == divisor_sum($y)$y;
}
say $sum;
# Or using a pipeline:
use Math::Prime::Util qw/vecsum divisor_sum/;
say vecsum( map { divisor_sum($_) }
grep { my $y = divisor_sum($_)$_;
$y > $_ && $_==(divisor_sum($y)$y) }
1 .. 10000 );
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 4factor 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 ($maxn, $maxratio) = (0,0);
foreach my $n (1..1000000) {
my $ndivphi = $n / euler_phi($n);
($maxn, $maxratio) = ($n, $ndivphi) if $ndivphi > $maxratio;
}
say "$maxn $maxratio";
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, 1 to 2 minutes:
use Math::Prime::Util qw/forcomposites factor/;
my $nsemis = 0;
forcomposites { $nsemis++ if scalar factor($_) == 2; } int(10**8)1;
say $nsemis;
Here is one of the best ways for PE187: under 20 milliseconds from the command
line. Much faster than Pari, and competitive with Mathematica.
use Math::Prime::Util qw/forprimes prime_count/;
my $limit = shift  int(10**8);
$limit;
my ($sum, $pc) = (0, 1);
forprimes {
$sum += prime_count(int($limit/$_)) + 1  $pc++;
} int(sqrt($limit));
say $sum;
To get the result of "matches" in Math::Factor::XS:
use Math::Prime::Util qw/divisors/;
sub matches {
my @d = divisors(shift);
return map { [$d[$_],$d[$#d$_]] } 1..(@d1)>>1;
}
my $n = 139650;
say "$n = ", join(" = ", map { "$_>[0]·$_>[1]" } matches($n));
or its "matches" function with the "skip_multiples" option:
sub matches {
my @d = divisors(shift);
return map { [$d[$_],$d[$#d$_]] }
grep { my $div=$d[$_]; !scalar(grep {!($div % $d[$_])} 1..$_1) }
1..(@d1)>>1; }
}
Compute OEIS A054903 <http://oeis.org/A054903> just like CRG4s 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..$n1;
say "$n ", join(" ", @r);
}
}
Find the 7digit palindromic primes in the first 20k digits of Pi:
use Math::Prime::Util qw/Pi is_prime/;
my $pi = "".Pi(20000); # make sure we only stringify once
for my $pos (2 .. length($pi)7) {
my $s = substr($pi, $pos, 7);
say "$s at $pos" if $s eq reverse($s) && is_prime($s);
}
# Or we could use the regex engine to find the palindromes:
while ($pi =~ /(([1379])(\d)(\d)\d\4\3\2)/g) {
say "$1 at ",pos($pi)7 if is_prime($1)
}
The Bell numbers <https://en.wikipedia.org/wiki/Bell_number> B_n:
sub B { my $n = shift; vecsum(map { stirling($n,$_,2) } 0..$n) }
say "$_ ",B($_) for 1..50;
Recognizing tetrahedral numbers (OEIS A000292 <http://oeis.org/A000292>):
sub is_tetrahedral {
my $n6 = vecprod(6,shift);
my $k = rootint($n6,3);
vecprod($k,$k+1,$k+2) == $n6;
}
Recognizing powerful numbers (e.g. "ispowerful" from Pari/GP):
sub ispowerful { 0 + vecall { $_>[1] > 1 } factor_exp(shift); }
Convert from binary to hex (3000x faster than Math::BaseConvert):
my $hex_string = todigitstring(fromdigits($bin_string,2),16);
Calculate and print derangements using permutations:
my @data = qw/a b c d/;
forperm { say "@data[@_]" unless vecany { $_[$_]==$_ } 0..$#_ } @data;
# Using forderange directly is faster
Compute the subfactorial of n (OEIS A000166 <http://oeis.org/A000166>):
sub subfactorial { my $n = shift;
vecsum(map{ vecprod((1)**($n$_),binomial($n,$_),factorial($_)) }0..$n);
}
Compute subfactorial (number of derangements) using simple recursion:
sub subfactorial { my $n = shift;
use bigint;
($n < 1) ? 1 : $n * subfactorial($n1) + (1)**$n;
}
Above "2^64", "is_prob_prime" performs an extrastrong BPSW
test <http://en.wikipedia.org/wiki/BailliePSW_primality_test> which is
fast (a little less than the time to perform 3 MillerRabin 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 MillerRabin 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" MR 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 multiprecision math
library used for Perl 6 and at least one CPAN Crypto module.
 GMP/MPIR
 Uses a set of "k" staticrandom 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 (not recent Pari/GP)
 Pari 2.1.7 is the default version installed with the Math::Pari module. It
uses 10 random MR 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 MR tests. For example, it will indicate 9 is prime
about 1 out of every 276k calls.
Basically the problem is that it is just too easy to get counterexamples from
running "k" MR 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.
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 64bit 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 64bit when using numbers larger than "~ 2^49").
The best solution is updating to a more recent Perl.
The module is threadsafe 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 nonCygwin Win32
and call these routines from within Perl threads.
Because the loop functions like "forprimes" use "MULTICALL",
there is some odd behavior with anonymous sub creation inside the block. This
is shared with most XS modules that use "MULTICALL", and is rarely
seen because it is such an unusual use. An example is:
forprimes { my $var = "p is $_"; push @subs, sub {say $var}; } 50;
$_>() for @subs;
This can be worked around by using double braces for the function, e.g.
"forprimes {{ ... }} 50".
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
wellwritten trial division, meaning it is very fast for small numbers, but
terribly slow for large 64bit numbers. MPU is similarly fast with small
numbers, but becomes faster as the size increases. MPXS's prime sieve is an
unoptimized nonsegmented 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.0001 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 23x slower than the base sieve for MPU,
and is nonsegmented so cannot be used for larger values. Since the functions
work with the sieve, they are very fast. The fast bitvectorlookup
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 (ndigit, nbit, ranged, strong, and ST) in
addition to Maurer's algorithm. MPU does not have the critical bug RT81858
<https://rt.cpan.org/Ticket/Display.html?id=81858>. MPU has 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 makes the speed vastly faster. Crypt::Primes is hardcoded to use
Crypt::Random which uses /dev/random (blocking source), while MPU uses its own
ChaCha20 implementation seeded from /dev/urandom or Win32. 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. Its functions do
not support bigints. Both are implemented with trial division, meaning they
are very fast for really small values, but become very slow as the input gets
larger (factoring 19 digit semiprimes is over 1000 times slower). The function
"count_prime_factors" can be done in MPU using "scalar
factor($n)". See the "EXAMPLES" section for a 2line function
replicating "matches".
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. 19digit 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 PollardRho 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).
List::Gen is very interesting and includes a builtin primes iterator as well as
a "is_prime" filter for arbitrary sequences. Unfortunately both are
very slow.
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 24x 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 onetoone mapping between functions in MPU, but it
does offer a way to get many similar results such as primes, twin primes,
SophieGermain 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 32bit indices.
"cr_combine" in Math::ModInt::ChineseRemainder is similar to MPU's
"chinese", and in fact they use the same algorithm. The former
module uses caching of moduli to speed up further operations. MPU does not do
this. This would only be important for cases where the lcm is larger than a
native int (noting that use in cryptography would always have large moduli).
For combinations and permutations there are many alternatives. One difference
with nearly all of them is that MPU's "forcomb" and
"forperm" functions don't operate directly on a user array but on
generic indices. Math::Combinatorics and Algorithm::Combinatorics have more
features, but will be slower. List::Permutor does permutations with an
iterator. Algorithm::FastPermute and Algorithm::Permute are very similar but
can be 210x faster than MPU (they use the same userblock structure but
twiddle the user array each call).
There are numerous modules to perform a set product (also called Cartesian
product or cross product). These include Set::Product,
Math::Cartesian::Product, Set::Scalar, and Set::CrossProduct, as well as a few
others. The Set::CartesianProduct::Lazy module provides random access, albeit
rather slowly. Our "forsetproduct" matches Set::Product in both high
performance and functionality (that module's single function
"product" in Set::Product is essentially identical to ours).
Math::Pari supports a lot of features, with a great deal of overlap. In general,
MPU will be faster for native 64bit integers, while it 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 MR tests
with randomly chosen bases (fixed seed, but doesn't reset each invocation
like GMP's "is_probab_prime"). This has a much greater chance of
false positives compared to the BPSW test  some composites such as 9,
88831, 38503, etc. (OEIS A141768 <http://oeis.org/A141768>) have a
surprisingly high chance of being indicated prime. Using
"isprime($n,1)" will perform an "n1" 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 APRCL 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 almostextrastrong 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
extrastrong BPSW test, which has an even lower chance of counterexample.
With Math::Prime::Util::GMP, "is_prime" adds an extra MR test
using a random base, 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 precalculation
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 64bit 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
220x faster for native integers. MPU also supported range inputs, which
can be much more efficient. With bigint arguments, MPU is slightly faster
than Math::Pari if the GMP backend is available, but very slow
without.
 "gcd", "lcm", "kronecker",
"znorder", "znprimroot", "znlog"
 Similar to MPU's "gcd", "lcm", "kronecker",
"znorder", "znprimroot", and "znlog". 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, sometimes
it hangs). MPU's "znprimroot" will always return the smallest
root if it exists, and "undef" otherwise. Similarly, MPU's
"znlog" will return the smallest "x" and work with
nonprimitiveroot "g", which is similar to Pari/GP 2.6, but not
the older versions in Math::Pari. The performance of "znlog" is
quite good compared to older Pari/GP, but much worse than 2.6's new
methods.
 "sigma"
 Similar to MPU's "divisor_sum". MPU is ~10x faster when the
result fits in a native integer. Once things overflow it is fairly similar
in performance. However, using Math::BigInt can slow things down quite a
bit, so for best performance in these cases using a Math::GMP object is
best.
 "numbpart", "forpart"
 Similar to MPU's "partitions" and "forpart". These
functions were introduced in Pari 2.3 and 2.6, hence are not in
Math::Pari. "numbpart" produce identical results to
"partitions", but Pari is much faster. forpart is very
similar to Pari's function, but produces a different ordering (MPU is the
standard antilexicographical, Pari uses a size sort). Currently Pari is
somewhat faster due to Perl function call overhead. When using
restrictions, Pari has much better optimizations.
 "eint1"
 Similar to MPU's "ExponentialIntegral".
 "zeta"
 MPU has "RiemannZeta" which takes nonnegative 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 Pari library
used is about 10 years old now). For native integers, typically Math::Pari
will be slower than MPU. For bigints, Math::Pari may be 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.
First, for those looking for the state of the art nonPerl solutions:
 Primality testing
 For general numbers smaller than 2000 or so digits, MPU is the fastest
solution I am aware of (it is faster than Pari 2.7, PFGW, and FLINT). For
very 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 it is commonly used for fast filtering of 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 David Cleaver's mpzaprcl
<http://sourceforge.net/projects/mpzaprcl/>, the APRCL from the
modern Pari <http://pari.math.ubordeaux.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 gmpecm
<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. Tomás 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.
Counting the primes to "800_000_000" (800 million):
Time (s) Module Version Notes
   
0.001 Math::Prime::Util 0.37 using extended LMO
0.007 Math::Prime::Util 0.12 using Lehmer's method
0.27 Math::Prime::Util 0.17 segmented mod30 sieve
0.39 Math::Prime::Util::PP 0.24 Perl (Lehmer's method)
0.9 Math::Prime::Util 0.01 mod30 sieve
2.9 Math::Prime::FastSieve 0.12 decent oddnumber 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) 201201 array of odds
548.1 RosettaCode sieve (net) 201206 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.
 Small inputs: is_prime from 1 to 20M

2.0s Math::Prime::Util (sieve lookup if prime_precalc used)
2.5s Math::Prime::FastSieve (sieve lookup)
3.3s Math::Prime::Util (trial + deterministic MR)
10.4s Math::Prime::XS (trial)
19.1s Math::Pari w/2.3.5 (BPSW)
52.4s Math::Pari (10 random MR)
480s Math::Primality (deterministic MR)
 Large native inputs: is_prime from 10^16 to 10^16 + 20M

4.5s Math::Prime::Util (BPSW)
24.9s Math::Pari w/2.3.5 (BPSW)
117.0s Math::Pari (10 random MR)
682s 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.2s Math::Prime::Util (BPSW + 1 random MR)
2.7s Math::Pari w/2.3.5 (BPSW)
13.0s Math::Primality (BPSW)
35.2s Math::Pari (10 random MR)
38.6s Math::Prime::Util w/o GMP (BPSW)
70.7s Math::Prime::Util (n1 or ECPP proof)
102.9s Math::Pari w/2.3.5 (APRCL 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 built with 2.3.5 not only has a better primality test versus
the default 2.1.7, but runs faster. It still has quite a bit of overhead
with native size integers. Pari/GP 2.5.0 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 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 may be 3000x faster). 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 64bit 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 1118 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 19100 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/>, gmpecm
<http://ecm.gforge.inria.fr/>, GGNFS
<http://sourceforge.net/projects/ggnfs/>, and Pari
<http://pari.math.ubordeaux.fr/>. The latest yafu should cover most
uses, with GGNFS likely only providing a benefit for numbers large enough to
warrant distributed processing.
The "n1" proving algorithm in Math::Prime::Util::GMP compares well to
the version included 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
200digit primes including creating a certificate. Times below 200 digits are
faster than Pari 2.3.5's APRCL 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 300digit primes taking
~5 seconds on average and typically under a minute for 500digits. For
primality proving with very large numbers, I recommend Primo
<http://www.ellipsa.eu/>.
Seconds per prime for random prime generation on a early 2015 Macbook Pro (2.7
GHz i5) with Math::BigInt::GMP and Math::Prime::Util::GMP installed.
bits random +testing Maurer ShwTylr CPMaurer
     
64 0.00002 +0.000009 0.00004 0.00004 0.019
128 0.00008 +0.00014 0.00018 0.00012 0.051
256 0.0004 +0.0003 0.00085 0.00058 0.13
512 0.0023 +0.0007 0.0048 0.0030 0.40
1024 0.019 +0.0033 0.034 0.025 1.78
2048 0.26 +0.014 0.41 0.25 8.02
4096 2.82 +0.11 4.4 3.0 66.7
8192 23.7 +0.65 50.8 38.7 929.4
random = random_nbit_prime (results pass BPSW)
random+ = additional time for 3 MR and a Frobenius test
maurer = random_maurer_prime
ShwTylr = random_shawe_taylor_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
MR and a Frobenius test. At these bit sizes, the chances a composite number
passes BPSW, three more MR 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. With GMP
installed this always uses Maurer's algorithm as it is the best compromise
between speed and diversity.
"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.
"random_shawe_taylor_prime" similarly constructs a provable prime. It
uses a simpler construction method. It is slightly faster than Maurer's
algorithm but provides less diversity (even fewer primes in the range are
selected, though for typical cryptographic sizes this is not important). The
Perl implementation uses a single large random seed followed by SHA256 as
specified by FIPS 1864. The GMP implementation uses the same FIPS 1864
algorithm but uses its own CSPRNG which may not be SHA256.
"maurer" in Crypt::Primes times are included for comparison. It is
reasonably fast for small sizes but gets slow as the size increases. It is 10
to 500 times slower than this module's GMP methods. It does not perform any
primality checks on the intermediate results or the final result (I highly
recommended running 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).
Dana Jacobsen <dana@acm.org>
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.
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 Tomás Oliveira e Silva. I also want to thank Kim
Walisch for the many discussions about prime counting.
 •
 Christian Axler, "New bounds for the prime counting function
π(x)", September 2014. For large values, improved limits
versus Dusart 2010. <http://arxiv.org/abs/1409.1780>
 •
 Christian Axler, "Über die PrimzahlZählfunktion, die
nte Primzahl und verallgemeinerte RamanujanPrimzahlen", January
2013. Prime count and nthprime bounds in more detail. Thesis in German,
but first part is easily read.
<http://docserv.uniduesseldorf.de/servlets/DerivateServlet/Derivate28284/pdfa1b.pdf>
 •
 Christian Bau, "The Extended MeisselLehmer Algorithm", 2003,
preprint with example C++ implementation. Very detailed
implementationspecific paper which was used for the implementation here.
Highly recommended for implementing a sievebased LMO.
<http://cs.swan.ac.uk/~csoliver/oksatlibrary/OKplatform/ExternalSources/sources/NumberTheory/ChristianBau/>
 •
 Manuel Benito and Juan L. Varona, "Recursive formulas related to the
summation of the Möbius function", The Open Mathematics
Journal, v1, pp 2534, 2007. Among many other things, shows a simple
formula for computing the Mertens functions with only n/3 Möbius
values (not as fast as Deléglise and Rivat, but really simple).
<http://www.unirioja.es/cu/jvarona/downloads/BenitoVaronaTOMATJMertens.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 620647.
<http://www.ams.org/journals/mcom/197529130/S00255718197503846731/S00255718197503846731.pdf>
 •
 W. J. Cody and Henry C. Thacher, Jr., "Rational Chebyshev
Approximations for the Exponential Integral E_1(x)", Mathematics
of Computation, v22, pp 641649, 1968.
 •
 W. J. Cody and Henry C. Thacher, Jr., "Chebyshev approximations for
the exponential integral Ei(x)", Mathematics of Computation,
v23, pp 289303, 1969.
<http://www.ams.org/journals/mcom/196923106/S00255718196902423492/>
 •
 W. J. Cody, K. E. Hillstrom, and Henry C. Thacher Jr., "Chebyshev
Approximations for the Riemann Zeta Function", "Mathematics of
Computation", v25, n115, pp 537547, July 1971.
 •
 Henri Cohen, "A Course in Computational Algebraic Number
Theory", Springer, 1996. Practical computational number theory from
the team lead of Pari <http://pari.math.ubordeaux.fr/>. Lots of
explicit algorithms.
 •
 Marc Deléglise and Joöl Rivat, "Computing the summation
of the Möbius function", Experimental Mathematics, v5,
n4, pp 291295, 1996. Enhances the Möbius computation in Lioen/van
de Lune, and gives a very efficient way to compute the Mertens function.
<http://projecteuclid.org/euclid.em/1047565447>
 •
 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 interested in prime number
bounds. <http://www.unilim.fr/laco/theses/1998/T1998_01.html>
 •
 Pierre Dusart, "Estimates of Some Functions Over Primes without
R.H.", preprint, 2010. Updates to the best nonRH bounds for prime
count and nth prime. <http://arxiv.org/abs/1002.0442/>
 •
 PierreAlain Fouque and Mehdi Tibouchi, "Close to Uniform Prime
Number Generation With Fewer Random Bits", preprint, 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>
 •
 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. 421432, 1994. Describes a nice way to
compute a range of Möbius values.
<http://walter.lioen.com/papers/LL94.pdf>
 •
 Ueli M. Maurer, "Fast Generation of Prime Numbers and Secure
PublicKey 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>
 •
 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>
 •
 OEIS: Primorial <http://oeis.org/wiki/Primorial>
 •
 Vincent Pegoraro and Philipp Slusallek, "On the Evaluation of the
ComplexValued Exponential Integral", Journal of Graphics, GPU,
and Game Tools, v15, n3, pp 183198, 2011.
<http://www.cs.utah.edu/~vpegorar/research/2011_JGT/paper.pdf>
 •
 William H. Press et al., "Numerical Recipes", 3rd edition.
 •
 Hans Riesel, "Prime Numbers and Computer Methods for
Factorization", Birkh?user, 2nd edition, 1994. Lots of information,
some code, easy to follow.
 •
 David M. Smith, "MultiplePrecision Exponential Integral and Related
Functions", ACM Transactions on Mathematical Software, v37,
n4, 2011. <http://myweb.lmu.edu/dmsmith/toms2011.pdf>
 •
 Douglas A. Stoll and Patrick Demichel , "The impact of ζ(s)
complex zeros on π(x) for x < 10^{10^{13}}",
"Mathematics of Computation", v80, n276, pp 23812394, October
2011.
<http://www.ams.org/journals/mcom/201180276/S002557182011024774/home.html>
Copyright 20112018 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.