MPU_NO_XS

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.

MPU_NO_GMP

If this variable is not set or set to anything other than 1, the module operates normally.

is_prime

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

Returns 0 is the number is composite, 1 if it is probably prime, and 2 if it is definitely prime. For numbers smaller than "2^64" it will only return 0 (composite) or 2 (definitely prime), as this range has been exhaustively tested and has no counterexamples. For larger numbers, an extra-strong BPSW test is used. If Math::Prime::Util::GMP is installed, some additional primality tests are also performed, and a quick attempt is made to perform a primality proof, so it will return 2 for many other inputs.

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

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

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

primes

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

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

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

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

Sieving will be done if required. The algorithm used will depend on the range and whether a sieve result already exists. Possibilities include primality testing (for very small ranges), a Sieve of Eratosthenes using wheel factorization, or a segmented sieve.

next_prime

  $n = next_prime($n);

Returns the next prime greater than the input number. The result will be a bigint if it can not be exactly represented in the native int type (larger than "4,294,967,291" in 32-bit Perl; larger than "18,446,744,073,709,551,557" in 64-bit).

prev_prime

  $n = prev_prime($n);

Returns the prime preceding the input number (i.e. the largest prime that is strictly less than the input). "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

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

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

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

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

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

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

forcomposites

  forcomposites { say } 1000;
  forcomposites { say } 2000,2020;

Given a block and either an end number or a start and end pair, calls the block for each composite in the inclusive range. The composites, OEIS A002808 <http://oeis.org/A002808>, are the numbers greater than 1 which are not prime: "4, 6, 8, 9, 10, 12, 14, 15, ...".

foroddcomposites

forsemiprimes

This is essentially equivalent to:

  forcomposites { if (is_semiprime($_)) { ... } }

forfactored

  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 read-only. Any attempt to modify them will result in undefined behavior.

This corresponds to the Pari/GP 2.10 "forfactored" function.

forsquarefree

This corresponds to the Pari/GP 2.10 "forsquarefree" function.

fordivisors

  fordivisors { $prod *= $_ } $n;

Given a block and a non-negative number "n", the block is called with $_ set to each divisor in sorted order. Also see "divisor_sum".

forpart

  forpart { say "@_" } 25;           # unrestricted partitions
  forpart { say "@_" } 25,{n=>5}     # ... with exactly 5 values
  forpart { say "@_" } 25,{nmax=>5}  # ... with <=5 values

Given a non-negative 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 non-negative 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 (non-zero) or non-prime (zero)

Like forcomb and forperm, the partition return values are read-only. Any attempt to modify them will result in undefined behavior.

forcomp

The number of unrestricted compositions of "n" is "2^(n-1)".

forcomb

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 read-only. Any attempt to modify them will result in undefined behavior.

If the second argument "k" is not supplied, then all k-subsets are returned starting with the smallest set "k=0" and continuing to "k=n". Each k-subset is in lexicographical order. This is the power set of "n".

This corresponds to the Pari/GP 2.10 "forsubset" function.

forperm

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 read-only. Any attempt to modify them will result in undefined behavior.

forderange

formultiperm

  # 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

  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 de-duplication 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 read-only. Hence modifying the array inside the subroutine will cause side-effects.

As with other iterators, the "lastfor" function will cause an early exit.

lastfor

  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.

prime_iterator

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

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

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

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

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

prime_iterator_object

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

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

prime_count

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

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

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

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

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

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

prime_count_upper

prime_count_lower

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

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

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

These routines use verified tight limits below a range at least "2^35". 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.

prime_count_approx

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

Returns an approximation to the "prime_count" function, without having to generate any primes. For values under "10^36" this uses the Riemann R function, which is quite accurate: an error of less than "0.0005%" is typical for input values over "2^32", and decreases as the input gets larger.

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

twin_primes

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.

twin_prime_count

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.

twin_prime_count_approx

semi_primes

This works just like the "primes" function. Like that function, an array reference is returned.

semiprime_count

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.

semiprime_count_approx

ramanujan_primes

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.

ramanujan_prime_count

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.

ramanujan_prime_count_approx

ramanujan_prime_count_lower

ramanujan_prime_count_upper

sieve_range

  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+width-1" where "n + offset" has no prime factors less than "depth".

sieve_prime_cluster

  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.

sum_primes

    $sum = 0; forprimes { $sum += $_ } $low,$high;  $sum;
    # or
    vecsum( @{ primes($low,$high) } );

but is much more efficient.

The current implementation is a small-table-enhanced 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

  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.

nth_prime

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

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

Like most programs with similar functionality, this is one-based. 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 low-biased estimate using the inverse logarithmic integral, use a fast prime count, then sieve in the small difference.

While this method is thousands of times faster than generating primes, and doesn't involve big tables of precomputed values, it still can take a fair amount of time for large inputs. Calculating the "10^12th" prime takes about 1 second, the "10^13th" prime takes under 10 seconds, and the "10^14th" prime (3475385758524527) takes under 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 (32-bit or 64-bit), the result will take a very long time. A later version of Math::Prime::Util::GMP may include this functionality which would help for 32-bit machines.

nth_prime_upper

nth_prime_lower

  my $lower_limit = nth_prime_lower($n);
  my $upper_limit = nth_prime_upper($n);
  # 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.

nth_prime_approx

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

Returns an approximation to the "nth_prime" function, without having to generate any primes. 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.

nth_twin_prime

nth_twin_prime_approx

nth_semiprime

nth_semiprime_approx

nth_ramanujan_prime

nth_ramanujan_prime_approx

nth_ramanujan_prime_lower

nth_ramanujan_prime_upper

is_pseudoprime

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 Solovay-Strassen and Miller-Rabin tests).

is_euler_pseudoprime

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 Solovay-Strassen test, which is a precursor to the Miller-Rabin 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.

is_strong_pseudoprime

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

Takes a positive number "n" and one or more non-zero 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 32-bit number). Given the small chances of passing multiple bases, there are some math packages that just use multiple MR tests for primality testing.

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

is_lucas_pseudoprime

is_strong_lucas_pseudoprime

is_extra_strong_lucas_pseudoprime

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

is_almost_extra_strong_lucas_pseudoprime

An optional second argument (an integer between 1 and 256) indicates the increment amount for "P" parameter selection. The default value of 1 yields the parameter selection described in "is_extra_strong_lucas_pseudoprime", creating a pseudoprime sequence which is a superset of the latter's pseudoprime sequence OEIS A217719 <http://oeis.org/A217719>. A value of 2 yields the method used by Pari <http://pari.math.u-bordeaux.fr/faq.html#primetest>.

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

is_euler_plumb_pseudoprime

The main reason for this test is that is a bit more efficient than other probable prime tests.

is_perrin_pseudoprime

While pseudoprimes are relatively rare (the first two are 271441 and 904631), infinitely many exist. They have significant overlap with the base-2 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 pre-filters, 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>.

is_catalan_pseudoprime

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.

is_frobenius_pseudoprime

Some authors use the Fibonacci polynomial "x^2-x-1" 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 extra-strong 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.

is_frobenius_underwood_pseudoprime

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.

is_frobenius_khashin_pseudoprime

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.

miller_rabin_random

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

is_prob_prime

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

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

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

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

is_bpsw_prime

is_provable_prime

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

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

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

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

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

prime_certificate

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

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

is_provable_prime_with_cert

verify_prime

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

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

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

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

  [MPU - Primality Certificate]

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

  Version 1.0

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

  Base 16

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

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

  Proof for:
  N 100003

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

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

  If Q is prime, then N is prime.

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

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

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

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

"Small"

  Type Small
  N 5791
    

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

"BLS3"

  Type BLS3
  N  2297612322987260054928384863
  Q  16501461106821092981
  A  5
    

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

"Pocklington"

  Type Pocklington
  N  2297612322987260054928384863
  Q  16501461106821092981
  A  5
    

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

"BLS15"

  Type BLS15
  N  8087094497428743437627091507362881
  Q  175806402118016161687545467551367
  LP 1
  LQ 22
    

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

"BLS5"

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

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

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

"ECPP"

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

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

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

is_aks_prime

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

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

While this is an important theoretical algorithm, and makes an interesting example, it is hard to overstate just how impractically slow it is in practice. It is not used for any purpose in non-theoretical work, as it is literally millions of times slower than other algorithms. From R.P. Brent, 2010: "AKS is not a practical algorithm. ECPP is much faster." We have ECPP, and indeed it is much faster.

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

is_mersenne_prime

  say "2^607-1 (M607) is a Mersenne prime" if is_mersenne_prime(607);

Takes a non-negative number "p" as input and returns 1 if the Mersenne number "2^p-1" 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 Lucas-Lehmer test.

The Lucas-Lehmer 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 general-form 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.

is_ramanujan_prime

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.

is_power

  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 non-negative 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 "k-th" 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 k-th power, then this will be set to the k-th 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.

is_prime_power

If a second argument is present, it must be a scalar reference. If the return value is non-zero, then it will be set to "p".

This corresponds to Pari/GP's "isprimepower" function.

is_square

This corresponds to Pari/GP's "issquare" function.

sqrtint

This corresponds to Pari/GP's "sqrtint" function.

rootint

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.

logint

  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 non-zero 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.

lucasu

  say "Fibonacci($_) = ", lucasu(1,-1,$_) for 0..100;

Given integers "P", "Q", and the non-negative 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.

lucasv

  say "Lucas($_) = ", lucasv(1,-1,$_) for 0..100;

Given integers "P", "Q", and the non-negative 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.

lucas_sequence

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

Computes "U_k", "V_k", and "Q_k" for the Lucas sequence defined by "P","Q", modulo "n". The modular Lucas sequence is used in a number of primality tests and proofs. The following conditions must hold: " |P| < n" ; " |Q| < n" ; " k >= 0" ; " n >= 2".

gcd

lcm

  lcm(0, n) = 0              Any zero in list results in zero return
  lcm(n,-m) = lcm(n, m)      We use the absolute values

gcdext

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.

chinese

  say chinese( [14,643], [254,419], [87,733] );  # 87041638

Solves a system of simultaneous congruences using the Chinese Remainder Theorem (with extension to non-coprime 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, ...}]

vecsum

  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.

vecprod

  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.

vecmin

  say "Smallest Totient 100k-200k: ", 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.

vecmax

  say "Largest Totient 100k-200k: ", 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.

vecreduce

  say "Count of non-zero 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 singular-element 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.

vecany

vecall

vecnone

vecnotall

vecfirst

  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 short-circuit and is less obvious.

vecfirstidx

  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.

vecextract

  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 zero-based indices into the first array. Duplicate values are allowed and the ordering is preserved. Hence these are equivalent:

    vecextract($aref, $iref);
    @$aref[@$iref];

todigits

  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 read-only. 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.

todigitstring

  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.

fromdigits

  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.

sumdigits

  # 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 non-digit characters of "n" are ignored (including negative signs and decimal points). This is similar to the command "vecsum(split(//,$n))" but faster, allows non-positive-integer 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 "base-1" 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))".

invmod

  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.

sqrtmod

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.

addmod

mulmod

powmod

divmod

valuation

  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.

hammingweight

is_square_free

  say "$n has no repeating factors" if is_square_free($n);

Returns 1 if the input "n" has no repeated factor.

is_carmichael

  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^(n-1) ≡ 1 mod n" for all "1 < b < n" relatively prime to "n". Alternately Korselt's theorem says these are composites such that "n" is square-free and "p-1" divides "n-1" 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 square-free-ness and the Korselt condition, while those that do not divide "n" are used as the pseudoprime base. The chances of a non-Carmichael passing this test are less than "2^-150".

This is the OEIS series A002997 <http://oeis.org/A002997>.

is_quasi_carmichael

This is the OEIS series A257750 <http://oeis.org/A257750>.

is_semiprime

The boolean result is the same as "scalar(factor(n)) == 2", but this function performs shortcuts that can greatly speed up the operation.

is_fundamental

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.

is_totient

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.

is_pillai

For n prime, this is the OEIS series A063980 <http://oeis.org/A063980>.

is_polygonal

If a third argument is present, it must be a scalar reference. It will be set to n if x is the nth s-gonal number. If the function returns 0, then it will be unchanged.

This corresponds to Pari's "ispolygonal" function.

moebius

  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 square-free (i.e. "n" has a repeated factor), and "-1^t" if "n" is a product of "t" distinct primes. This is an important function in prime number theory. Like SAGE, we define "moebius(0) = 0" for convenience.

If called with two arguments, they define a range "low" to "high", and the function returns an array with the value of the Möbius function for every n from low to high inclusive. Large values of high will result in a lot of memory use. The algorithm used 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 read-only constants. This should almost never come up, but it means trying to modify aliased return values will cause an exception (modifying the returned scalar or array is fine).

mertens

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

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

   time    approx mem
     0.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 non-segmented "n^1/2" version of their method. Kuznetsov (2011) gives an alternate method that he indicates is even faster. Lastly, one of the advanced prime count algorithms could be theoretically used to create a faster solution.

euler_phi

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

Returns φ(n), the Euler totient function (also called Euler's phi or phi function) for an integer value. This is an arithmetic function 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 pre-2.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.

inverse_totient

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

jordan_totient

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

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

ramanujan_sum

exp_mangoldt

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

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

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

liouville

chebyshev_theta

  say chebyshev_theta(10000);

Returns θ(n), the first Chebyshev function for a non-negative integer input. This is the sum of the logarithm of each prime where "p <= n". This is effectively:

  my $s = 0;  forprimes { $s += log($_) } $n;  return $s;

chebyshev_psi

  say chebyshev_psi(10000);

Returns ψ(n), the second Chebyshev function for a non-negative integer input. This is the sum of the logarithm of each prime 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;

divisor_sum

  say "Sum of divisors of $n:", divisor_sum( $n );
  say "sigma_2($n) = ", divisor_sum($n, 2);
  say "Number of divisors: sigma_0($n) = ", divisor_sum($n, 0);

This function takes a positive integer as input and returns the sum of its divisors, including 1 and itself. An optional second argument "k" may be given, which will result in the sum of the "k-th" powers of the divisors to be returned.

This is known as the sigma function (see Hardy and Wright section 16.7). 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.

ramanujan_tau

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.

primorial

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

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

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

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

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

pn_primorial

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

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

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

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

consecutive_integer_lcm

  $lcm = consecutive_integer_lcm($n);

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

partitions

This uses a combinatorial calculation, which means it will not be very fast compared to Pari, Mathematica, or FLINT which use the Rademacher formula using multi-precision floating point. In 10 seconds:

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

carmichael_lambda

kronecker

   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 non-residue 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.

factorial

factorialmod

While very efficient, this is not state of the art. Currently, Fredrik Johansson's fast multi-point 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.

binomial

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.

hclassno

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

bernfrac

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.

bernreal

This corresponds to Pari's "bernreal" function.

stirling

  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 non-negative 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^(n-k)" 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" non-empty subsets. The Lah numbers are the number of ways to split a set of "n" elements into "k" non-empty lists.

harmfrac

Binary splitting (Fredrik Johansson's elegant formulation) is used.

This corresponds to Mathematica's "HarmonicNumber" function.

harmreal

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.

znorder

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

znprimroot

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.

is_primitive_root

znlog

  $k = znlog($a, $g, $p)

Returns the integer "k" that solves the equation "a = g^k mod p", or undef if no solution is found. This is the discrete logarithm problem.

The implementation for native integers first applies Silver-Pohlig-Hellman 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 memory-heavy BSGS being used.

legendre_phi

  $phi = legendre_phi(1000000000, 41);

Given a non-negative integer "n" and a non-negative prime number "a", returns the Legendre phi function (also called Legendre's sum). This is the count of positive integers <= "n" which are not divisible by any of the first "a" primes.

inverse_li

  $approx_prime_count = inverse_li(1000000000);

Given a non-negative 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.

numtoperm

  @p = numtoperm(10,654321);  # @p=(1,8,2,7,6,5,3,4,9,0)

Given a non-negative 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.

permtonum

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

randperm

  @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 "n-1".

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.

shuffle

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

OVERVIEW

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 64-bit 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 512-bit state which is significantly larger than the 48-bit "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 non-CSPRNG 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".

irand

  $n32 = irand;     # random 32-bit integer

Returns a random 32-bit integer using the CSPRNG.

irand64

  $n64 = irand64;   # random 64-bit integer

Returns a random 64-bit integer using the CSPRNG (on 64-bit Perl).

drand

  $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 32-bit 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 not-very-good 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".

random_bytes

  $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 - AES-CTR

    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)

entropy_bytes

The actual performance will be highly system dependent.

urandomb

  $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^b-1" inclusive.

urandomm

  $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 "n-1" inclusive. Care is taken to prevent modulo bias.

csrand

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.

srand

The API is nearly identical to the system function "srand". It uses a UV which can be 64-bit rather than always 32-bit. 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.

rand

random_factored_integer

  my($n, $factors) = random_factored_integer(1000000);

Given a positive non-zero 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.

random_prime

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

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

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

For small numbers, a random index selection is done, which gives ideal uniformity and is very efficient with small inputs. For ranges larger than this ~16-bit threshold but within the native bit size, a Monte Carlo method is used. 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.

random_ndigit_prime

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

Selects a random n-digit prime, where the input is an integer number of digits. One of the primes within that range (e.g. 1000 - 9999 for 4-digits) will be uniformly selected.

If the number of digits is greater than or equal to the maximum native type, then the result will be returned as a BigInt. However, if the "nobigint" configuration option is on, then output will be restricted to native size numbers, and requests for more digits than natively supported will result in an error. For better performance with large bit sizes, install Math::Prime::Util::GMP.

random_nbit_prime

  my $bigprime = random_nbit_prime(512);

Selects a random n-bit prime, where the input is an integer number of bits. A prime with the nth bit set will be uniformly selected.

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.

random_strong_prime

  my $bigprime = random_strong_prime(512);

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

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

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

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

random_proven_prime

  my $bigprime = random_proven_prime(512);

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

random_proven_prime_with_cert

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

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

random_maurer_prime

  my $bigprime = random_maurer_prime(512);

Construct an n-bit provable prime, using the FastPrime algorithm of Ueli Maurer (1995). This is the same algorithm used by Crypt::Primes. Similar to "random_nbit_prime", the result will be a BigInt if the number of bits is greater than the native bit size.

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.

random_maurer_prime_with_cert

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

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

random_shawe_taylor_prime

  my $bigprime = random_shawe_taylor_prime(8192);

Construct an n-bit provable prime, using the Shawe-Taylor algorithm in section C.6 of FIPS 186-4. This uses 512 bits of randomness and SHA-256 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 186-4 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.

random_shawe_taylor_prime_with_cert

  my($n, $cert) = random_shawe_taylor_prime_with_cert(4096)

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

random_semiprime

The factors will be approximately equal size, which is typical for cryptographic use. For example, a 64-bit semiprime of this type is the product of two 32-bit 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.

random_unrestricted_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

  prime_precalc( 1_000_000_000 );

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

prime_memfree

  prime_memfree;

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

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

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

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

prime_get_config

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

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

  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

  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
               re-enable.  You probably will never want to do this.

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

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

  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.

factor

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

Produces the prime factors of a positive number input, in numerical order. The product of the returned factors will be equal to the input. "n = 1" will return an empty list, and "n = 0" will return 0. This matches Pari.

In scalar context, returns Ω(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 p-1. The combination is applied to each non-prime factor found.

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

factor_exp

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

Produces pairs of prime factors and exponents in numerical factor order. This is more convenient for some algorithms. This is the same form that Mathematica's "FactorInteger[n]" and Pari/GP's "factorint" functions return. Note that Math::Pari transposes the Pari result matrix.

In scalar context, returns ω(n), the number of unique prime factors (OEIS A001221 <http://oeis.org/A001221>). This corresponds to Pari's omega(n) function and Mathematica's "PrimeNu[n]" function. This is same result that we would get if we evaluated the resulting array in scalar context.

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

divisors

  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.

trial_factor

  my @factors = trial_factor($n);

Produces the prime factors of a positive number input. The factors will be in numerical order. For large inputs this will be very slow. Like all the specific-algorithm *_factor routines, this is not exported unless explicitly requested.

fermat_factor

  my @factors = fermat_factor($n);

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

holf_factor

  my @factors = holf_factor($n);

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

lehman_factor

  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 42-bit inputs: "n < 8796393022208".

squfof_factor

  my @factors = squfof_factor($n);

Produces factors, not necessarily prime, of the positive number input. An optional number of rounds can be given as a second parameter. It is possible the function will be unable to find a factor, in which case a single element, the input, is returned. This function typically runs very fast.

prho_factor

pbrent_factor

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

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

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

pminus1_factor

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

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

pplus1_factor

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

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

ecm_factor

  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.

ExponentialIntegral

  my $Ei = ExponentialIntegral($x);

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

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

For non-BigInt/BigFloat 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").

LogarithmicIntegral

  my $li = LogarithmicIntegral($x)

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

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

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

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

For BigInt / BigFloat inputs, full accuracy and performance is obtained only if Math::Prime::Util::GMP is installed.

RiemannZeta

  my $z = RiemannZeta($s);

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

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

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

For BigInt / BigFloat 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 high-accuracy computations when GMP is not used.

RiemannR

  my $r = RiemannR($x);

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

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

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

For BigInt / BigFloat 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.

LambertW

This function handles all real value inputs with non-complex 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.

Pi

  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.

PRIME COUNTS

  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 mod-30 sieve
       0.39  Math::Prime::Util::PP       0.24     Perl (Lehmer's method)
       0.9   Math::Prime::Util           0.01     mod-30 sieve
       2.9   Math::Prime::FastSieve      0.12     decent odd-number sieve
      11.7   Math::Prime::XS             0.26     needs some optimization
      15.0   Bit::Vector                 7.2
      48.9   Math::Prime::Util::PP       0.14     Perl (fastest I know of)
     170.0   Faster Perl sieve (net)     2012-01  array of odds
     548.1   RosettaCode sieve (net)     2012-06  simplistic Perl
    3048.1   Math::Primality             0.08     Perl + Math::GMPz
  >20000     Math::Big                   1.12     Perl, > 26GB RAM used

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

PRIMALITY TESTING

Small inputs: is_prime from 1 to 20M

    2.0s  Math::Prime::Util      (sieve lookup if prime_precalc used)
    2.5s  Math::Prime::FastSieve (sieve lookup)
    3.3s  Math::Prime::Util      (trial + deterministic M-R)
   10.4s  Math::Prime::XS        (trial)
   19.1s  Math::Pari w/2.3.5     (BPSW)
   52.4s  Math::Pari             (10 random M-R)
  480s    Math::Primality        (deterministic M-R)
    

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 M-R)
  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 M-R)
    2.7s  Math::Pari w/2.3.5         (BPSW)
   13.0s  Math::Primality            (BPSW)
   35.2s  Math::Pari                 (10 random M-R)
   38.6s  Math::Prime::Util w/o GMP  (BPSW)
   70.7s  Math::Prime::Util          (n-1 or ECPP proof)
  102.9s  Math::Pari w/2.3.5         (APR-CL proof)
    

• MPU is consistently the fastest solution, and performs the most stringent probable prime tests on bigints.
• Math::Primality has a lot of overhead that makes it quite slow for native size integers. With bigints we finally see it work well.
• Math::Pari 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

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

PRIMALITY PROVING

RANDOM PRIME GENERATION

  bits    random   +testing   Maurer   Shw-Tylr  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 M-R and a Frobenius test
  maurer    = random_maurer_prime
  Shw-Tylr  = 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 M-R and a Frobenius test. At these bit sizes, the chances a composite number passes BPSW, three more M-R tests, and a Frobenius test is extraordinarily small.

"random_proven_prime" provides a randomly selected prime with an optional certificate, without specifying the particular method. 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 SHA-256 as specified by FIPS 186-4. The GMP implementation uses the same FIPS 186-4 algorithm but uses its own CSPRNG which may not be SHA-256.

"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).