Provided by: libmath-prime-util-perl_0.73-2build5_amd64 bug

NAME

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

VERSION

       Version 0.73

SYNOPSIS

         # Nothing is exported by default.  List the functions, or use :all.
         use Math::Prime::Util ':all';  # import all functions

         # The ':rand' tag replaces srand and rand (not done by default)
         use Math::Prime::Util ':rand';  # import srand, rand, irand, irand64

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

         # Get all divisors other than 1 and n
         my @divisors = divisors( $n );
         # Or just apply a block for each one
         my $sum = 0; fordivisors  { $sum += $_ + $_*$_ }  $n;

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

         # Moebius function used to calculate Mertens
         $sum += moebius($_) for (1..200); say "Mertens(200) = $sum";
         # Mertens function directly (more efficient for large values)
         say mertens(10_000_000);
         # Exponential of Mangoldt function
         say "lamba(49) = ", log(exp_mangoldt(49));
         # Some more number theoretical functions
         say liouville(4292384);
         say chebyshev_psi(234984);
         say chebyshev_theta(92384234);
         say partitions(1000);
         # Show all prime partitions of 25
         forpart { say "@_" unless scalar grep { !is_prime($_) } @_ } 25;
         # List all 3-way combinations of an array
         my @cdata = qw/apple bread curry donut eagle/;
         forcomb { say "@cdata[@_]" } @cdata, 3;
         # or all permutations
         forperm { say "@cdata[@_]" } @cdata;

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

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

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

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

         # Free any memory used by the module.
         prime_memfree;

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

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

DESCRIPTION

       A module for number theory in Perl.  This includes prime sieving, primality tests,
       primality proofs, integer factoring, counts / bounds / approximations for primes, nth
       primes, and twin primes, random prime generation, and much more.

       This module is the fastest on CPAN for almost all operations it supports.  This includes
       Math::Prime::XS, Math::Prime::FastSieve, Math::Factor::XS, Math::Prime::TiedArray,
       Math::Big::Factors, Math::Factoring, and Math::Primality (when the GMP module is
       available).  For numbers in the 10-20 digit range, it is often orders of magnitude faster.
       Typically it is faster than Math::Pari for 64-bit operations.

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

       The module is thread-safe and allows concurrency between Perl threads while still sharing
       a prime cache.  It is not itself multi-threaded.  See the Limitations section if you are
       using Win32 and threads in your program.  Also note that Math::Pari is not thread-safe
       (and will crash as soon as it is loaded in threads), so if you use Math::BigInt::Pari
       rather than Math::BigInt::GMP or the default backend, things will go pear-shaped.

       Two scripts are also included and installed by default:

       •   primes.pl displays primes between start and end values or expressions, with many
           options for filtering (e.g. twin, safe, circular, good, lucky, etc.).  Use "--help" to
           see all the options.

       •   factor.pl operates similar to the GNU "factor" program.  It supports bigint and
           expression inputs.

ENVIRONMENT VARIABLES

       There are two environment variables that affect operation.  These are typically used for
       validation of the different methods or to simulate systems that have different support.

   MPU_NO_XS
       If set to 1 then everything is run in pure Perl.  No C functions are loaded or used, as
       XSLoader is not even called.  All top-level XS functions are replaced by a pure Perl layer
       (the PPFE.pm module that supplies a "Pure Perl Front End").

       Caveat: This does not change whether the GMP backend is used.  For as much pure Perl as
       possible, you will need to set both variables.

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

   MPU_NO_GMP
       If set to 1 then the Math::Prime::Util::GMP backend is not loaded, and operation will be
       exactly as if it was not installed.

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

BIGNUM SUPPORT

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

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

       •   Install a recent version of Math::Prime::Util::GMP, as that will vastly increase the
           speed of many of the functions.  This does require the GMP <http://gmplib.org> library
           be installed on your system, but this increasingly comes pre-installed or easily
           available using the OS vendor package installation tool.

       •   Install and use Math::BigInt::GMP or Math::BigInt::Pari, then use "use bigint try =>
           'GMP,Pari'" in your script, or on the command line "-Mbigint=lib,GMP".  Large modular
           exponentiation is much faster using the GMP or Pari backends, as are the math and
           approximation functions when called with very large inputs.

       •   I have run these functions on many versions of Perl, and my experience is that if
           you're using anything older than Perl 5.14, I would recommend you upgrade if you are
           using bignums a lot.  There are some brittle behaviors on 5.12.4 and earlier with
           bignums.  For example, the default BigInt backend in older versions of Perl will
           sometimes convert small results to doubles, resulting in corrupted output.

PRIMALITY TESTING

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

FUNCTIONS

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

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

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

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

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

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

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

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

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

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

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

   next_prime
         $n = next_prime($n);

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

   prev_prime
         $n = prev_prime($n);

       Returns the prime preceding the input number (i.e. the largest prime that is strictly less
       than the input).  "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
       Similar to "forcomposites", but skipping all even numbers.  The odd composites, OEIS
       A071904 <http://oeis.org/A071904>, are the numbers greater than 1 which are not prime and
       not divisible by two: "9, 15, 21, 25, 27, 33, 35, ...".

   forsemiprimes
       Similar to "forcomposites", but only giving composites with exactly two factors.  The
       semiprimes, OEIS A001358 <http://oeis.org/A001358>, are the products of two primes: "4, 6,
       9, 10, 14, 15, 21, 22, 25, ...".

       This is essentially equivalent to:

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

   forfactored
         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
       Similar to "forfactored", but skipping numbers in the range that have a repeated factor.
       Inside the block, the moebius function can be cheaply computed as "((scalar(@_) & 1) ? -1
       : 1)" or similar.

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

   fordivisors
         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
       Similar to "forpart", but iterates over integer compositions rather than partitions.  This
       can be thought of as all ordering of partitions, or alternately partitions may be viewed
       as an ordered subset of compositions.  The ordering is lexicographic.  All options from
       "forpart" may be used.

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

   forcomb
       Given non-negative arguments "n" and "k", the block is called with @_ set to the "k"
       element array of values from 0 to "n-1" representing the combinations in lexicographical
       order.  While the binomial function gives the total number, this function can be used to
       enumerate the choices.

       Rather than give a data array as input, an integer is used for "n".  A convenient way to
       map to array elements is:

         forcomb { say "@data[@_]" } @data, 3;

       where the block maps the combination array @_ to array values, the argument for "n" is
       given the array since it will be evaluated as a scalar and hence give the size, and the
       argument for "k" is the desired size of the combinations.

       Like forpart and forperm, the index return values are 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
       Given non-negative argument "n", the block is called with @_ set to the "k" element array
       of values from 0 to "n-1" representing permutations in lexicographical order.  The total
       number of calls will be "n!".

       Rather than give a data array as input, an integer is used for "n".  A convenient way to
       map to array elements is:

         forperm { say "@data[@_]" } @data;

       where the block maps the permutation array @_ to array values, and the argument for "n" is
       given the array since it will be evaluated as a scalar and hence give the size.

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

   forderange
       Similar to forperm, but iterates over derangements.  This is the set of permutations
       skipping any which maps an element to its original position.

   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
       Returns the lesser of twin primes between the lower and upper limits (inclusive), with a
       lower limit of 2 if none is given.  This is OEIS A001359 <http://oeis.org/A001359>.  Given
       a twin prime pair "(p,q)" with "q = p + 2", "p prime", and <q prime>, this function uses
       "p" to represent the pair.  Hence the bounds need to include "p", and the returned list
       will have "p" but not "q".

       This works just like the "primes" function, though only the first primes of twin prime
       pairs are returned.  Like that function, an array reference is returned.

   twin_prime_count
       Similar to prime count, but returns the count of twin primes (primes "p" where "p+2" is
       also prime).  Takes either a single number indicating a count from 2 to the argument, or
       two numbers indicating a range.

       The primes being counted are the first value, so a range of "(3,5)" will return a count of
       two, because both 3 and 5 are counted as twin primes.  A range of "(12,13)" will return a
       count of zero, because neither "12+2" nor "13+2" are prime.  In contrast, "primesieve"
       requires all elements of a constellation to be within the range to be counted, so would
       return one for the first example (5 is not counted because its pair 7 is not in the
       range).

       There is no useful formula known for this, unlike prime counts.  We sieve for the answer,
       using some small table acceleration.

   twin_prime_count_approx
       Returns an approximation to the twin prime count of "n".  This returns quickly and has a
       very small error for large values.  The method used is conjecture B of Hardy and
       Littlewood 1922, as stated in Sebah and Gourdon 2002.  For inputs under 10M, a correction
       factor is additionally applied to reduce the mean squared error.

   semi_primes
       Returns an array reference to semiprimes between the lower and upper limits (inclusive),
       with a lower limit of 4 if none is given.  This is OEIS A001358 <http://oeis.org/A001358>.
       The semiprimes are composite integers which are products of exactly two primes.

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

   semiprime_count
       Similar to prime count, but returns the count of semiprimes (composites with exactly two
       factors).  Takes either a single number indicating a count from 2 to the argument, or two
       numbers indicating a range.

       A fast method that requires computation only to the square root of the range end is used,
       unless the range is so small that walking it is faster.

   semiprime_count_approx
       Returns an approximation to the semiprime count of "n".  This returns quickly and is
       typically square root accurate.

   ramanujan_primes
       Returns the Ramanujan primes R_n between the upper and lower limits (inclusive), with a
       lower limit of 2 if none is given.  This is OEIS A104272 <http://oeis.org/A104272>.  These
       are the Rn such that if "x > Rn" then "prime_count"(n) - "prime_count"(n/2) >= "n".

       This has a similar API to the "primes" and "twin_primes" functions, and like them, returns
       an array reference.

       Generating Ramanujan primes takes some effort, including overhead to cover a range.  This
       will be substantially slower than generating standard primes.

   ramanujan_prime_count
       Similar to prime count, but returns the count of Ramanujan primes.  Takes either a single
       number indicating a count from 2 to the argument, or two numbers indicating a range.

       While not nearly as efficient as prime_count, this does use a number of speedups that
       result it in being much more efficient than generating all the Ramanujan primes.

   ramanujan_prime_count_approx
       A fast approximation of the count of Ramanujan primes under "n".

   ramanujan_prime_count_lower
       A fast lower limit on the count of Ramanujan primes under "n".

   ramanujan_prime_count_upper
       A fast upper limit on the count of Ramanujan primes under "n".

   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
       Returns the summation of primes between the lower and upper limits (inclusive), with a
       lower limit of 2 if none is given.  This is essentially similar to either of:

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

       but is much more efficient.

       The current implementation is a 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
       Returns the Nth twin prime.  This is done via sieving and counting, so is not very fast
       for large values.

   nth_twin_prime_approx
       Returns an approximation to the Nth twin prime.  A curve fit is used for small inputs
       (under 1200), while for larger inputs a binary search is done on the approximate twin
       prime count.

   nth_semiprime
       Returns the Nth semiprime, similar to where a "forsemiprimes" loop would end after "N"
       iterations, but much more efficiently.

   nth_semiprime_approx
       Returns an approximation to the Nth semiprime.  Curve fitting is used to get a fairly
       close approximation that is orders of magnitude better than the simple "n log n / log log
       n" approximation for large "n".

   nth_ramanujan_prime
       Returns the Nth Ramanujan prime.  For reasonable size values of "n", e.g.  under "10^8" or
       so, this is relatively efficient for single calls.  If multiple calls are being made, it
       will be much more efficient to get the list once.

   nth_ramanujan_prime_approx
       A fast approximation of the Nth Ramanujan prime.

   nth_ramanujan_prime_lower
       A fast lower limit on the Nth Ramanujan prime.

   nth_ramanujan_prime_upper
       A fast upper limit on the Nth Ramanujan prime.

   is_pseudoprime
       Takes a positive number "n" and one or more non-zero positive bases as input.  Returns 1
       if the input is a probable prime to each base, 0 if not.  This is the simple Fermat
       primality test.  Removing primes, given base 2 this produces the sequence OEIS A001567
       <http://oeis.org/A001567>.

       For practical use, "is_strong_pseudoprime" is a much stronger test with similar or better
       performance.

       Note that there is a set of composites (the Carmichael numbers) that will pass this test
       for all bases.  This downside is not shared by the Euler and strong probable prime tests
       (also called the Solovay-Strassen and Miller-Rabin tests).

   is_euler_pseudoprime
       Takes a positive number "n" and one or more non-zero positive bases as input.  Returns 1
       if the input is an Euler probable prime to each base, 0 if not.  This is the Euler test,
       sometimes called the Euler-Jacobi test.  Removing primes, given base 2 this produces the
       sequence OEIS A047713 <http://oeis.org/A047713>.

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

       This test forms the basis of the 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
       Takes a positive number as input, and returns 1 if the input is a standard Lucas probable
       prime using the Selfridge method of choosing D, P, and Q (some sources call this a Lucas-
       Selfridge pseudoprime).  Removing primes, this produces the sequence OEIS A217120
       <http://oeis.org/A217120>.

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

   is_extra_strong_lucas_pseudoprime
       Takes a positive number as input, and returns 1 if the input passes the extra strong Lucas
       test (as defined in Grantham 2000 <http://www.ams.org/mathscinet-getitem?mr=1680879>).
       This test has more stringent conditions than the strong Lucas test, and produces about 60%
       fewer pseudoprimes.  Performance is typically 20-30% faster than the strong Lucas test.

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

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

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

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

   is_euler_plumb_pseudoprime
       Takes a positive number "n" as input and returns 1 if "n" passes Colin Plumb's Euler
       Criterion primality test.  Pseudoprimes to this test are a subset of the base 2 Fermat and
       Euler tests, but a superset of the base 2 strong pseudoprime (Miller-Rabin) test.

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

   is_perrin_pseudoprime
       Takes a positive number "n" as input and returns 1 if "n" divides P(n) where P(n) is the
       Perrin number of "n".  The Perrin sequence is defined by "P(n) = P(n-2) + P(n-3)" with
       "P(0) = 3, P(1) = 0, P(2) = 2".

       While pseudoprimes are relatively rare (the first two are 271441 and 904631), infinitely
       many exist.  They have significant overlap with the 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
       Takes a positive number "n" as input and returns 1 if "-1^((n-1/2)) C_((n-1/2)" is
       congruent to 2 mod "n", where "C_n" is the nth Catalan number.  The nth Catalan number is
       equal to "binomial(2n,n)/(n+1)".  All odd primes satisfy this condition, and only three
       known composites.

       The pseudoprime sequence is OEIS A163209 <http://oeis.org/A163209>.

       There is no known efficient method to perform the Catalan primality test, so it is a
       curiosity rather than a practical test.  The implementation uses a method from Charles
       Greathouse IV (2015) and results from Aebi and Cairns (2008) to produce results many
       orders of magnitude faster than other known implementations, but it is still vastly slower
       than other compositeness tests.

   is_frobenius_pseudoprime
       Takes a positive number "n" as input, and two optional parameters "a" and "b", and returns
       1 if the "n" is a Frobenius probable prime with respect to the polynomial "x^2 - ax + b".
       Without the parameters, "b = 2" and "a" is the least positive odd number such that
       "(a^2-4b|n) = -1".  This selection has no pseudoprimes below "2^64" and none known.  In
       any case, the discriminant "a^2-4b" must not be a perfect square.

       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
       Takes a positive number as input, and returns 1 if the input passes the efficient
       Frobenius test of Paul Underwood.  This selects a parameter "a" as the least non-negative
       integer such that "(a^2-4|n)=-1", then verifies that "(x+2)^(n+1) = 2a + 5 mod
       (x^2-ax+1,n)".  This combines a Fermat and Lucas test with a cost of only slightly more
       than 2 strong pseudoprime tests.  This makes it similar to, but faster than, a regular
       Frobenius test.

       There are no known pseudoprimes to this test and extensive computation has shown no
       counterexamples under "2^50".  This test also has no overlap with the BPSW test, making it
       a very effective method for adding additional certainty.  Performance at 1e12 is about 60%
       slower than BPSW.

   is_frobenius_khashin_pseudoprime
       Takes a positive number as input, and returns 1 if the input passes the Frobenius test of
       Sergey Khashin.  This ensures "n" is not a perfect square, selects the parameter "c" as
       the smallest odd prime such that "(c|n)=-1", then verifies that "(1+D)^n = (1-D) mod n"
       where "D = sqrt(c) mod n".

       There are no known pseudoprimes to this test and Khashin (2018) shows there are no
       counterexamples under "2^64".  Performance at 1e12 is about 40% slower than BPSW.

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

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

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

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

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

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

   is_bpsw_prime
       Given a positive number input, returns 0 (composite), 2 (definitely prime), or 1 (probably
       prime), using the BPSW primality test (extra-strong variant).  Normally one of the
       "is_prime" in Math::Prime::Util or "is_prob_prime" in Math::Prime::Util functions will
       suffice, but those functions do pre-tests to find easy composites.  If you know this is
       not necessary, then calling "is_bpsw_prime" may save a small amount of time.

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

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

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

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

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

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

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

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

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

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

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

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

         [MPU - Primality Certificate]

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

         Version 1.0

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

         Base 16

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

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

         Proof for:
         N 100003

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

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

         If Q is prime, then N is prime.

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

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

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

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

       "Small"
             Type Small
             N 5791

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

       "BLS3"
             Type BLS3
             N  2297612322987260054928384863
             Q  16501461106821092981
             A  5

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

       "Pocklington"
             Type Pocklington
             N  2297612322987260054928384863
             Q  16501461106821092981
             A  5

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

       "BLS15"
             Type BLS15
             N  8087094497428743437627091507362881
             Q  175806402118016161687545467551367
             LP 1
             LQ 22

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

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

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

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

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

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

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

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

       Takes a 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
       Takes a positive number "n" as input and returns back either 0 or 1, indicating whether
       "n" is a Ramanujan prime.  Numbers that can be produced by the functions
       "ramanujan_primes" and "nth_ramanujan_prime" will return 1, while all other numbers will
       return 0.

       There is no simple function for this predicate, so Ramanujan primes through at least "n"
       are generated, then a search is performed for "n".  This is not efficient for multiple
       calls.

   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
       Given an integer input "n", returns "k" if "n = p^k" for some prime p, and zero otherwise.

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

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

   is_square
       Given a positive integer "n", returns 1 if "n" is a perfect square, 0 otherwise.  This is
       identical to "is_power(n,2)".

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

   sqrtint
       Given a non-negative integer input "n", returns the integer square root.  For native
       integers, this is equal to "int(sqrt(n))".

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

   rootint
       Given an non-negative integer "n" and positive exponent "k", return the integer k-th root
       of "n".  This is the largest integer "r" such that "r^k <= n".

       If a third argument is present, it must be a scalar reference.  It will be set to "r^k".

       Technically if "n" is negative and "k" is odd, the root exists and is equal to "sign(n) *
       |rootint(abs(n),k)".  It was decided to follow the behavior of Pari/GP and Math::BigInt
       and disallow negative "n".

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

   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
       Given a list of integers, returns the greatest common divisor.  This is often used to test
       for coprimality <https://oeis.org/wiki/Coprimality>.

   lcm
       Given a list of integers, returns the least common multiple.  Note that we follow the
       semantics of Mathematica, Pari, and Perl 6, re:

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

   gcdext
       Given two integers "x" and "y", returns "u,v,d" such that "d = gcd(x,y)" and "u*x + v*y =
       d".  This uses the extended Euclidian algorithm to compute the values satisfying Bézout's
       Identity.

       This corresponds to Pari's "gcdext" function, which was renamed from "bezout" out Pari
       2.6.  The results will hence match "bezout" in Math::Pari.

   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
       Given two integers "a" and "n", return the square root of "a" mod "n".  If no square root
       exists, undef is returned.  If defined, the return value "r" will always satisfy "r^2 = a
       mod n".

       If the modulus is prime, the function will always return "r", the smaller of the two
       square roots (the other being "-r mod p".  If the modulus is composite, one of possibly
       many square roots will be returned, and it will not necessarily be the smallest.

   addmod
       Given three integers "a", "b", and "n" where "n" is positive, return "(a+b) mod n".  This
       is particularly useful when dealing with numbers that are larger than a half-word but
       still native size.  No bigint package is needed and this can be 10-200x faster than using
       one.

   mulmod
       Given three integers "a", "b", and "n" where "n" is positive, return "(a*b) mod n".  This
       is particularly useful when "n" fits in a native integer.  No bigint package is needed and
       this can be 10-200x faster than using one.

   powmod
       Given three integers "a", "b", and "n" where "n" is positive, return "(a ** b) mod n".
       Typically binary exponentiation is used, so the process is very efficient.  With native
       size inputs, no bigint library is needed.

   divmod
       Given three integers "a", "b", and "n" where "n" is positive, return "(a/b) mod n".  This
       is done as "(a * (1/b mod n)) mod n".  If no inverse of "b" mod "n" exists then undef if
       returned.

   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
       Given an integer "n", returns the binary Hamming weight of abs(n).  This is also called
       the population count, and is the number of 1s in the binary representation.  This
       corresponds to Pari's "hammingweight" function for "t_INT" arguments.

   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
       Returns 0 if the input "n" is not a quasi-Carmichael number, and the number of bases
       otherwise.  These are square-free composites that satisfy "p+b" divides "n+b" for all
       prime factors "p" or "n" and for one or more non-zero integer "b".

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

   is_semiprime
       Given a positive integer "n", returns 1 if "n" is a semiprime, 0 otherwise.  A semiprime
       is the product of exactly two primes.

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

   is_fundamental
       Given an integer "d", returns 1 if "d" is a fundamental discriminant, 0 otherwise.  We
       consider 1 to be a fundamental discriminant.

       This is the OEIS series A003658 <http://oeis.org/A003658> (positive) and OEIS series
       A003657 <http://oeis.org/A003657> (negative).

       This corresponds to Pari's "isfundamental" function.

   is_totient
       Given an integer "n", returns 1 if there exists an integer "x" where "euler_phi(x) == n".

       This corresponds to Pari's "istotient" function, though without the optional second
       argument to return an "x".  Math::NumSeq::Totient also has a similar function.

       Also see "inverse_totient" which gives the count or list of values that produce a given
       totient.  This function is more efficient than getting the full count or list.

   is_pillai
       Given a positive integer "n", if there exists a "v" where "v! % n == n-1" and "n % v !=
       1", then "v" is returned.  Otherwise 0.

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

   is_polygonal
       Given integers "x" and "s", return 1 if x is an s-gonal number, 0 otherwise.  "s" must be
       greater than 2.

       If a third argument is present, it must be a scalar reference.  It will be set to n if x
       is the nth 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 array context, given a positive integer "n", returns the complete list of values "x"
       where "euler_phi(x) = n".  This can be a memory intensive operation if there are many
       values.

       In scalar context, returns just the count of values.  This is faster and uses
       substantially less memory.  The list/scalar distinction is similar to "factor" and
       "divisors".

       This roughly corresponds to the Maple function "InverseTotient", and the hidden
       Mathematica function "EulerPhiInverse".  The algorithm used is from Max Alekseyev (2016).

   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
       Returns Ramanujan's sum of the two positive variables "k" and "n".  This is the sum of the
       nth powers of the primitive k-th roots of unity.

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

   chebyshev_theta
         say chebyshev_theta(10000);

       Returns θ(n), the first Chebyshev function for a non-negative integer input.  This is the
       sum of the logarithm of each prime where "p <= n".  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
       Takes a positive integer as input and returns the value of Ramanujan's tau function.  The
       result is a signed integer.  This corresponds to Pari v2.8's "tauramanujan" function and
       Mathematica's "RamanujanTau" function.

       This currently uses a simple method based on divisor sums, which does not have a good
       computational growth rate.  Pari's implementation uses Hurwitz class numbers and is more
       efficient for large inputs.

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

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

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

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

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

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

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

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

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

   consecutive_integer_lcm
         $lcm = consecutive_integer_lcm($n);

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

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

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

                  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
       Returns the Carmichael function (also called the reduced totient function, or Carmichael
       λ(n)) of a positive integer argument.  It is the smallest positive integer "m" such that
       "a^m = 1 mod n" for every integer "a" coprime to "n".  This is OEIS series A002322
       <http://oeis.org/A002322>.

   kronecker
       Returns the Kronecker symbol "(a|n)" for two integers.  The possible return values with
       their meanings for odd prime "n" are:

          0   a = 0 mod n
          1   a is a quadratic residue mod n       (a = x^2 mod n for some x)
         -1   a is a quadratic 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
       Given positive integer argument "n", returns the factorial of "n", defined as the product
       of the integers 1 to "n" with the special case of "factorial(0) = 1".  This corresponds to
       Pari's factorial(n) and Mathematica's "Factorial[n]" functions.

   factorialmod
       Given two positive integer arguments "n" and "m", returns "n! mod m".  This is much faster
       than computing the large factorial(n) followed by a mod operation.

       While very efficient, this is not state of the art.  Currently, Fredrik Johansson's fast
       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
       Given integer arguments "n" and "k", returns the binomial coefficient
       "n*(n-1)*...*(n-k+1)/k!", also known as the choose function.  Negative arguments use the
       Kronenburg extensions <http://arxiv.org/abs/1105.3689/>.  This corresponds to Pari's
       "binomial(n,k)" function, Mathematica's "Binomial[n,k]" function, and GMP's "mpz_bin_ui"
       function.

       For negative arguments, this matches Mathematica.  Pari does not implement the "n < 0, k
       <= n" extension and instead returns 0 for this case.  GMP's API does not allow negative
       "k" but otherwise matches.  Math::BigInt does not implement any extensions and the results
       for "n < 0, k " 0> are undefined.

   hclassno
       Returns 12 times the Hurwitz-Kronecker class number of the input integer "n".  This will
       always be an integer due to the pre-multiplication by 12.  The result is 0 for any input
       less than zero or congruent to 1 or 2 mod 4.

       This is related to Pari's qfbhclassno(n) where hclassno(n) for positive "n" equals "12 *
       qfbhclassno(n)" in Pari/GP.  This is OEIS A259825 <http://oeis.org/A259825>.

   bernfrac
       Returns the Bernoulli number "B_n" for an integer argument "n", as a rational number
       represented by two Math::BigInt objects.  B_1 = 1/2.  This corresponds to Pari's
       bernfrac(n) and Mathematica's "BernoulliB" functions.

       Having a modern version of Math::Prime::Util::GMP installed will make a big difference in
       speed.  That module uses a fast Pi/Zeta method.  Our pure Perl backend uses the Seidel
       method as shown by Peter Luschny.  This is faster than Math::Pari which uses an older
       algorithm, but quite a bit slower than modern Pari, Mathematica, or our GMP backend.

       This corresponds to Pari's "bernfrac" function and Mathematica's "BernoulliB" function.

   bernreal
       Returns the Bernoulli number "B_n" for an integer argument "n", as a Math::BigFloat object
       using the default precision.  An optional second argument may be given specifying the
       precision to be used.

       This corresponds to Pari's "bernreal" function.

   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
       Returns the Harmonic number "H_n" for an integer argument "n", as a rational number
       represented by two Math::BigInt objects.  The harmonic numbers are the sum of reciprocals
       of the first "n" natural numbers: "1 + 1/2 + 1/3 + ... + 1/n".

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

       This corresponds to Mathematica's "HarmonicNumber" function.

   harmreal
       Returns the Harmonic number "H_n" for an integer argument "n", as a Math::BigFloat object
       using the default precision.  An optional second argument may be given specifying the
       precision to be used.

       For large "n" values, using a lower precision may result in faster computation as an
       asymptotic formula may be used.  For precisions of 13 or less, native floating point is
       used for even more speed.

   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
       Given a positive integer "n", returns the smallest primitive root of "(Z/nZ)^*", or
       "undef" if no root exists.  A root exists when "euler_phi($n) == carmichael_lambda($n)",
       which will be true for all prime "n" and some composites.

       OEIS A033948 <http://oeis.org/A033948> is a sequence of integers where the primitive root
       exists, while OEIS A046145 <http://oeis.org/A046145> is a list of the smallest primitive
       roots, which is what this function produces.

   is_primitive_root
       Given two non-negative numbers "a" and "n", returns 1 if "a" is a primitive root modulo
       "n", and 0 if not.  If "a" is a primitive root, then euler_phi(n) is the smallest "e" for
       which "a^e = 1 mod n".

   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.

RANDOM NUMBERS

   OVERVIEW
       Prior to version 5.20, Perl's "rand" function used the system rand function.  This meant
       it varied by system, and was almost always a poor choice.  For 5.20, Perl standardized on
       "drand48" and includes the source so there are no system dependencies.  While this was an
       improvement, "drand48" is not a good PRNG.  It really only has 32 bits of random values,
       and fails many statistical tests.  See <http://www.pcg-random.org/statistical-tests.html>
       for more information.

       There are much better choices for standard random number generators, such as the Mersenne
       Twister, PCG, or Xoroshiro128+.  Someday perhaps Perl will get one of these to replace
       drand48.  In the mean time, Math::Random::MTwist provides numerous features and excellent
       performance, or this module.

       Since we often deal with random primes for cryptographic purposes, we have additional
       requirements.  This module uses a CSPRNG for its random stream.  In particular, ChaCha20,
       which is the same algorithm used by BSD's "arc4random" and "/dev/urandom" on BSD and Linux
       4.8+.  Seeding is performed at startup using the Win32 Crypto API (on Windows),
       "/dev/urandom", "/dev/random", or Crypt::PRNG, whichever is found first.

       We use the original ChaCha definition rather than RFC7539.  This means a 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
       Similar to random_bytes, but directly using the entropy source.  This is not normally
       recommended as it can consume shared system resources and is typically slow -- on the
       computer that produced the "random_bytes" chart above, using "dd" generated the same 13
       MB/s performance as our "entropy_bytes" function.

       The actual performance will be highly system dependent.

   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
       Takes a binary string "data" as input and seeds the internal CSPRNG.  This is not normally
       needed as system entropy is used as a seed on startup.  For best security this should be
       16-128 bytes of good entropy.  No more than 1024 bytes will be used.

       With no argument, reseeds using system entropy, which is preferred.

       If the "secure" configuration has been set, then this will croak if given an argument.
       This allows for control of reseeding with entropy the module gets itself, but not user
       supplied.

   srand
       Takes a single UV argument and seeds the CSPRNG with it, as well as returning it.  If no
       argument is given, a new UV seed is constructed.  Note that this creates a very weak seed
       from a cryptographic standpoint, so it is useful for testing or simulations but "csrand"
       is recommended, or keep using the system entropy default seed.

       The API is nearly identical to the system function "srand".  It uses a UV which can be
       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
       An alias for "drand", not exported unless the ":rand" tag is used.

   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 PRIMES

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

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

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

       For small numbers, a random index selection is done, which gives ideal uniformity and is
       very efficient with small inputs.  For ranges larger than this ~16-bit threshold but
       within the native bit size, a Monte Carlo method is used.  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 sr-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
       Takes a positive integer number of bits "bits", returns a random semiprime of exactly
       "bits" bits.  The result has exactly two prime factors (hence semiprime).

       The factors will be approximately equal size, which is typical for cryptographic use.  For
       example, a 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
       Takes a positive integer number of bits "bits", returns a random semiprime of exactly
       "bits" bits.  The result has exactly two prime factors (hence semiprime).

       The factors are uniformly selected from the universe of all "bits"-bit semiprimes.  This
       means semiprimes with one factor equal to 2 will be most common, 3 next most common, etc.
       "bits" must be 3 or greater.

       Some effort is taken to select uniformly from the universe of "bits"-bit semiprimes.  This
       takes slightly longer than some methods that do not select uniformly.

UTILITY FUNCTIONS

   prime_precalc
         prime_precalc( 1_000_000_000 );

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

   prime_memfree
         prime_memfree;

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

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

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

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

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

         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.

FACTORING FUNCTIONS

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

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

       In scalar context, returns Ω(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.

MATHEMATICAL FUNCTIONS

   ExponentialIntegral
         my $Ei = ExponentialIntegral($x);

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

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

       For non-BigInt/BigFloat 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
       Returns the principal branch of the Lambert W function of a real value.  Given a value "k"
       this solves for "W" in the equation "k = We^W".  The input must not be less than "-1/e".
       This corresponds to Pari's "lambertw" function and Mathematica's "ProductLog" / "LambertW"
       function.

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

EXAMPLES

       Print Fibonacci numbers:

           perl -Mntheory=:all -E 'say lucasu(1,-1,$_) for 0..20'

       Print strong pseudoprimes to base 17 up to 10M:

           # Similar to A001262's isStrongPsp function, but much faster
           perl -MMath::Prime::Util=:all -E 'forcomposites { say if is_strong_pseudoprime($_,17) } 10000000;'

       Print some primes above 64-bit range:

           perl -MMath::Prime::Util=:all -Mbigint -E 'my $start=100000000000000000000; say join "\n", @{primes($start,$start+1000)}'

           # Another way
           perl -MMath::Prime::Util=:all -E 'forprimes { say } "100000000000000000039", "100000000000000000993"'

           # Similar using Math::Pari:
           # perl -MMath::Pari=:int,PARI,nextprime -E 'my $start = PARI "100000000000000000000"; my $end = $start+1000; my $p=nextprime($start); while ($p <= $end) { say $p; $p = nextprime($p+1); }'

       Generate Carmichael numbers (OEIS A002997 <http://oeis.org/A002997>):

           perl -Mntheory=:all -E 'foroddcomposites { say if is_carmichael($_) } 1e6;'

           # Less efficient, similar to Mathematica or MAGMA:
           perl -Mntheory=:all -E 'foroddcomposites { say if $_ % carmichael_lambda($_) == 1 } 1e6;'

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

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

       Construct and use a Sophie-Germain prime iterator:

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

       Project Euler, problem 3 (Largest prime factor):

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

       Project Euler, problem 7 (10001st prime):

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

       Project Euler, problem 10 (summation of primes):

         use Math::Prime::Util qw/sum_primes/;
         say sum_primes(2_000_000);
         #  ... or do it a little more manually ...
         use Math::Prime::Util qw/forprimes/;
         my $sum = 0;
         forprimes { $sum += $_ } 2_000_000;
         say $sum;
         #  ... or do it using a big list ...
         use Math::Prime::Util qw/vecsum primes/;
         say vecsum( @{primes(2_000_000)} );

       Project Euler, problem 21 (Amicable numbers):

         use Math::Prime::Util qw/divisor_sum/;
         my $sum = 0;
         foreach my $x (1..10000) {
           my $y = divisor_sum($x)-$x;
           $sum += $x + $y if $y > $x && $x == divisor_sum($y)-$y;
         }
         say $sum;
         # Or using a pipeline:
         use Math::Prime::Util qw/vecsum divisor_sum/;
         say vecsum( map { divisor_sum($_) }
                     grep { my $y = divisor_sum($_)-$_;
                            $y > $_ && $_==(divisor_sum($y)-$y) }
                     1 .. 10000 );

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

         perl -MMath::Prime::Util=primes -MList::Util=first -E 'say first { /1/&&/2/&&/3/&&/4/&&/5/&&/6/&&/7/} reverse @{primes(1000000,9999999)};'

       Project Euler, problem 47 (Distinct primes factors):

         use Math::Prime::Util qw/pn_primorial factor_exp/;
         my $n = pn_primorial(4);  # Start with the first 4-factor number
         # factor_exp in scalar context returns the number of distinct prime factors
         $n++ while (factor_exp($n) != 4 || factor_exp($n+1) != 4 || factor_exp($n+2) != 4 || factor_exp($n+3) != 4);
         say $n;

       Project Euler, problem 69, stupid brute force solution (about 1 second):

         use Math::Prime::Util qw/euler_phi/;
         my ($maxn, $maxratio) = (0,0);
         foreach my $n (1..1000000) {
           my $ndivphi = $n / euler_phi($n);
           ($maxn, $maxratio) = ($n, $ndivphi) if $ndivphi > $maxratio;
         }
         say "$maxn  $maxratio";

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

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

       Project Euler, problem 187, stupid brute force solution, 1 to 2 minutes:

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

       Here is one of the best ways for PE187:  under 20 milliseconds from the command line.
       Much faster than Pari, and competitive with Mathematica.

         use Math::Prime::Util qw/forprimes prime_count/;
         my $limit = shift || int(10**8);
         $limit--;
         my ($sum, $pc) = (0, 1);
         forprimes {
           $sum += prime_count(int($limit/$_)) + 1 - $pc++;
         } int(sqrt($limit));
         say $sum;

       To get the result of "matches" in Math::Factor::XS:

         use Math::Prime::Util qw/divisors/;
         sub matches {
           my @d = divisors(shift);
           return map { [$d[$_],$d[$#d-$_]] } 1..(@d-1)>>1;
         }
         my $n = 139650;
         say "$n = ", join(" = ", map { "$_->[0]·$_->[1]" } matches($n));

       or its "matches" function with the "skip_multiples" option:

         sub matches {
           my @d = divisors(shift);
           return map { [$d[$_],$d[$#d-$_]] }
                  grep { my $div=$d[$_]; !scalar(grep {!($div % $d[$_])} 1..$_-1) }
                  1..(@d-1)>>1; }
         }

       Compute OEIS A054903 <http://oeis.org/A054903> just like CRG4s Pari example:

         use Math::Prime::Util qw/forcomposite divisor_sum/;
         forcomposites {
           say if divisor_sum($_)+6 == divisor_sum($_+6)
         } 9,1e7;

       Construct the table shown in OEIS A046147 <http://oeis.org/A046147>:

         use Math::Prime::Util qw/znorder euler_phi gcd/;
         foreach my $n (1..100) {
           if (!znprimroot($n)) {
             say "$n -";
           } else {
             my $phi = euler_phi($n);
             my @r = grep { gcd($_,$n) == 1 && znorder($_,$n) == $phi } 1..$n-1;
             say "$n ", join(" ", @r);
           }
         }

       Find the 7-digit palindromic primes in the first 20k digits of Pi:

         use Math::Prime::Util qw/Pi is_prime/;
         my $pi = "".Pi(20000);  # make sure we only stringify once
         for my $pos (2 .. length($pi)-7) {
           my $s = substr($pi, $pos, 7);
           say "$s at $pos" if $s eq reverse($s) && is_prime($s);
         }

         # Or we could use the regex engine to find the palindromes:
         while ($pi =~ /(([1379])(\d)(\d)\d\4\3\2)/g) {
           say "$1 at ",pos($pi)-7 if is_prime($1)
         }

       The Bell numbers <https://en.wikipedia.org/wiki/Bell_number> B_n:

         sub B { my $n = shift; vecsum(map { stirling($n,$_,2) } 0..$n) }
         say "$_  ",B($_) for 1..50;

       Recognizing tetrahedral numbers (OEIS A000292 <http://oeis.org/A000292>):

         sub is_tetrahedral {
           my $n6 = vecprod(6,shift);
           my $k  = rootint($n6,3);
           vecprod($k,$k+1,$k+2) == $n6;
         }

       Recognizing powerful numbers (e.g. "ispowerful" from Pari/GP):

         sub ispowerful { 0 + vecall { $_->[1] > 1 } factor_exp(shift); }

       Convert from binary to hex (3000x faster than Math::BaseConvert):

         my $hex_string = todigitstring(fromdigits($bin_string,2),16);

       Calculate and print derangements using permutations:

         my @data = qw/a b c d/;
         forperm { say "@data[@_]" unless vecany { $_[$_]==$_ } 0..$#_ } @data;
         # Using forderange directly is faster

       Compute the subfactorial of n (OEIS A000166 <http://oeis.org/A000166>):

         sub subfactorial { my $n = shift;
           vecsum(map{ vecprod((-1)**($n-$_),binomial($n,$_),factorial($_)) }0..$n);
         }

       Compute subfactorial (number of derangements) using simple recursion:

         sub subfactorial { my $n = shift;
           use bigint;
           ($n < 1)  ?  1  :  $n * subfactorial($n-1) + (-1)**$n;
         }

PRIMALITY TESTING NOTES

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

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

       libtommath
           Uses the first "k" prime bases.  This is problematic for cryptographic use, as there
           are known methods (e.g. Arnault 1994) for constructing counterexamples.  The number of
           bases required to avoid false results is unreasonably high, hence performance is slow
           even if one ignores counterexamples.  Unfortunately this is the multi-precision math
           library used for Perl 6 and at least one CPAN Crypto module.

       GMP/MPIR
           Uses a set of "k" static-random bases.  The bases are randomly chosen using a PRNG
           that is seeded identically each call (the seed changes with each release).  This
           offers a very slight advantage over using the first "k" prime bases, but not much.
           See, for example, Nicely's mpz_probab_prime_p pseudoprimes
           <http://www.trnicely.net/misc/mpzspsp.html> page.

       Math::Pari (not recent Pari/GP)
           Pari 2.1.7 is the default version installed with the Math::Pari module.  It uses 10
           random M-R bases (the PRNG uses a fixed seed set at compile time).  Pari 2.3.0 was
           released in May 2006 and it, like all later releases through at least 2.6.1, use BPSW
           / APRCL, after complaints of false results from using M-R tests.  For example, it will
           indicate 9 is prime about 1 out of every 276k calls.

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

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

LIMITATIONS

       Perl versions earlier than 5.8.0 have problems doing exact integer math.  Some operations
       will flip signs, and many operations will convert intermediate or output results to
       doubles, which loses precision on 64-bit systems.  This causes numerous functions to not
       work properly.  The test suite will try to determine if your Perl is broken (this only
       applies to really old versions of Perl compiled for 64-bit when using numbers larger than
       "~ 2^49").  The best solution is updating to a more recent Perl.

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

       Because the loop functions like "forprimes" use "MULTICALL", there is some odd behavior
       with anonymous sub creation inside the block.  This is shared with most XS modules that
       use "MULTICALL", and is rarely seen because it is such an unusual use.  An example is:

         forprimes { my $var = "p is $_"; push @subs, sub {say $var}; } 50;
         $_->() for @subs;

       This can be worked around by using double braces for the function, e.g.  "forprimes {{ ...
       }} 50".

SEE ALSO

       This section describes other CPAN modules available that have some feature overlap with
       this one.  Also see the "REFERENCES" section.  Please let me know if any of this
       information is inaccurate.  Also note that just because a module doesn't match what I
       believe are the best set of features doesn't mean it isn't perfect for someone else.

       I will use SoE to indicate the Sieve of Eratosthenes, and MPU to denote this module
       (Math::Prime::Util).  Some quick alternatives I can recommend if you don't want to use
       MPU:

       •   Math::Prime::FastSieve is the alternative module I use for basic functionality with
           small integers.  It's fast and simple, and has a good set of features.

       •   Math::Primality is the alternative module I use for primality testing on bigints.  The
           downside is that it can be slow, and the functions other than primality tests are very
           slow.

       •   Math::Pari if you want the kitchen sink and can install it and handle using it.  There
           are still some functions it doesn't do well (e.g. prime count and nth_prime).

       Math::Prime::XS has "is_prime" and "primes" functionality.  There is no bigint support.
       The "is_prime" function uses well-written trial division, meaning it is very fast for
       small numbers, but terribly slow for large 64-bit numbers.  MPU is similarly fast with
       small numbers, but becomes faster as the size increases.  MPXS's prime sieve is an
       unoptimized non-segmented SoE which returns an array.  Sieve bases larger than "10^7"
       start taking inordinately long and using a lot of memory (gigabytes beyond "10^10").  E.g.
       "primes(10**9, 10**9+1000)" takes 36 seconds with MPXS, but only 0.0001 seconds with MPU.

       Math::Prime::FastSieve supports "primes", "is_prime", "next_prime", "prev_prime",
       "prime_count", and "nth_prime".  The caveat is that all functions only work within the
       sieved range, so are limited to about "10^10".  It uses a fast SoE to generate the main
       sieve.  The sieve is 2-3x slower than the base sieve for MPU, and is non-segmented so
       cannot be used for larger values.  Since the functions work with the sieve, they are very
       fast.  The fast bit-vector-lookup functionality can be replicated in MPU using
       "prime_precalc" but is not required.

       Bit::Vector supports the "primes" and "prime_count" functionality in a somewhat similar
       way to Math::Prime::FastSieve.  It is the slowest of all the XS sieves, and has the most
       memory use.  It is faster than pure Perl code.

       Crypt::Primes supports "random_maurer_prime" functionality.  MPU has more options for
       random primes (n-digit, n-bit, ranged, strong, and S-T) in addition to Maurer's algorithm.
       MPU does not have the critical bug RT81858
       <https://rt.cpan.org/Ticket/Display.html?id=81858>.  MPU has a more uniform distribution
       as well as return a larger subset of primes (RT81871
       <https://rt.cpan.org/Ticket/Display.html?id=81871>).  MPU does not depend on Math::Pari
       though can run slow for bigints unless the Math::BigInt::GMP or Math::BigInt::Pari modules
       are installed.  Having Math::Prime::Util::GMP installed makes the speed vastly faster.
       Crypt::Primes is hardcoded to use Crypt::Random which uses /dev/random (blocking source),
       while MPU uses its own ChaCha20 implementation seeded from /dev/urandom or Win32.  MPU can
       return a primality certificate.  What Crypt::Primes has that MPU does not is the ability
       to return a generator.

       Math::Factor::XS calculates prime factors and factors, which correspond to the "factor"
       and "divisors" functions of MPU.  Its functions do not support bigints.  Both are
       implemented with trial division, meaning they are very fast for really small values, but
       become very slow as the input gets larger (factoring 19 digit semiprimes is over 1000
       times slower).  The function "count_prime_factors" can be done in MPU using "scalar
       factor($n)".  See the "EXAMPLES" section for a 2-line function replicating "matches".

       Math::Big version 1.12 includes "primes" functionality.  The current code is only usable
       for very tiny inputs as it is incredibly slow and uses lots of memory.  RT81986
       <https://rt.cpan.org/Ticket/Display.html?id=81986> has a patch to make it run much faster
       and use much less memory.  Since it is in pure Perl it will still run quite slow compared
       to MPU.

       Math::Big::Factors supports factorization using wheel factorization (smart trial
       division).  It supports bigints.  Unfortunately it is extremely slow on any input that
       isn't the product of just small factors.  Even 7 digit inputs can take hundreds or
       thousands of times longer to factor than MPU or Math::Factor::XS.  19-digit semiprimes
       will take hours versus MPU's single milliseconds.

       Math::Factoring is a placeholder module for bigint factoring.  Version 0.02 only supports
       trial division (the Pollard-Rho method does not work).

       Math::Prime::TiedArray allows random access to a tied primes array, almost identically to
       what MPU provides in Math::Prime::Util::PrimeArray.  MPU has attempted to fix
       Math::Prime::TiedArray's shift bug (RT58151
       <https://rt.cpan.org/Ticket/Display.html?id=58151>).  MPU is typically much faster and
       will use less memory, but there are some cases where MP:TA is faster (MP:TA stores all
       entries up to the largest request, while MPU:PA stores only a window around the last
       request).

       List::Gen is very interesting and includes a built-in primes iterator as well as a
       "is_prime" filter for arbitrary sequences.  Unfortunately both are very slow.

       Math::Primality supports "is_prime", "is_pseudoprime", "is_strong_pseudoprime",
       "is_strong_lucas_pseudoprime", "next_prime", "prev_prime", "prime_count", and
       "is_aks_prime" functionality.  This is a great little module that implements primality
       functionality.  It was the first CPAN module to support the BPSW test.  All inputs are
       processed using GMP, so it of course supports bigints.  In fact, Math::Primality was made
       originally with bigints in mind, while MPU was originally targeted to native integers, but
       both have added better support for the other.  The main differences are extra
       functionality (MPU has more functions) and performance.  With native integer inputs, MPU
       is generally much faster, especially with "prime_count".  For bigints, MPU is slower
       unless the Math::Prime::Util::GMP module is installed, in which case MPU is 2-4x faster.
       Math::Primality also installs a "primes.pl" program, but it has much less functionality
       than the one included with MPU.

       Math::NumSeq does not have a one-to-one mapping between functions in MPU, but it does
       offer a way to get many similar results such as primes, twin primes, Sophie-Germain
       primes, lucky primes, moebius, divisor count, factor count, Euler totient, primorials,
       etc.  Math::NumSeq is set up for accessing these values in order rather than for arbitrary
       values, though a few sequences support random access.  The primary advantage I see is the
       uniform access mechanism for a lot of sequences.  For those methods that overlap, MPU is
       usually much faster.  Importantly, most of the sequences in Math::NumSeq are limited to
       32-bit indices.

       "cr_combine" in Math::ModInt::ChineseRemainder is similar to MPU's "chinese", and in fact
       they use the same algorithm.  The former module uses caching of moduli to speed up further
       operations.  MPU does not do this.  This would only be important for cases where the lcm
       is larger than a native int (noting that use in cryptography would always have large
       moduli).

       For combinations and permutations there are many alternatives.  One difference with nearly
       all of them is that MPU's "forcomb" and "forperm" functions don't operate directly on a
       user array but on generic indices.  Math::Combinatorics and Algorithm::Combinatorics have
       more features, but will be slower.  List::Permutor does permutations with an iterator.
       Algorithm::FastPermute and Algorithm::Permute are very similar but can be 2-10x faster
       than MPU (they use the same user-block structure but twiddle the user array each call).

       There are numerous modules to perform a set product (also called Cartesian product or
       cross product).  These include Set::Product, Math::Cartesian::Product, Set::Scalar, and
       Set::CrossProduct, as well as a few others.  The Set::CartesianProduct::Lazy module
       provides random access, albeit rather slowly.  Our "forsetproduct" matches Set::Product in
       both high performance and functionality (that module's single function "product" in
       Set::Product is essentially identical to ours).

       Math::Pari supports a lot of features, with a great deal of overlap.  In general, MPU will
       be faster for native 64-bit integers, while it differs for bigints (Pari will always be
       faster if Math::Prime::Util::GMP is not installed; with it, it varies by function).  Note
       that Pari extends many of these functions to other spaces (Gaussian integers, complex
       numbers, vectors, matrices, polynomials, etc.) which are beyond the realm of this module.
       Some of the highlights:

       "isprime"
           The default Math::Pari is built with Pari 2.1.7.  This uses 10 M-R tests with randomly
           chosen bases (fixed seed, but doesn't reset each invocation like GMP's
           "is_probab_prime").  This has a much greater chance of false positives compared to the
           BPSW test -- some composites such as 9, 88831, 38503, etc.  (OEIS A141768
           <http://oeis.org/A141768>) have a surprisingly high chance of being indicated prime.
           Using "isprime($n,1)" will perform an "n-1" proof, but this becomes unreasonably slow
           past 70 or so digits.

           If Math::Pari is built using Pari 2.3.5 (this requires manual configuration) then the
           primality tests are completely different.  Using "ispseudoprime" will perform a BPSW
           test and is quite a bit faster than the older test.  "isprime" now does an APR-CL
           proof (fast, but no certificate).

           Math::Primality uses a strong BPSW test, which is the standard BPSW test based on the
           1980 paper.  It has no known counterexamples (though like all these tests, we know
           some exist).  Pari 2.3.5 (and through at least 2.6.2) uses an almost-extra-strong BPSW
           test for its "ispseudoprime" function.  This is deterministic for native integers, and
           should be excellent for bigints, with a slightly lower chance of counterexamples than
           the traditional strong test.  Math::Prime::Util uses the full extra-strong BPSW test,
           which has an even lower chance of counterexample.  With Math::Prime::Util::GMP,
           "is_prime" adds an extra M-R test using a random base, which further reduces the
           probability of a composite being allowed to pass.

       "primepi"
           Only available with version 2.3 of Pari.  Similar to MPU's "prime_count" function in
           API, but uses a naive counting algorithm with its precalculated primes, so is not of
           practical use.  Incidently, Pari 2.6 (not usable from Perl) has fixed the pre-
           calculation requirement so it is more useful, but is still thousands of times slower
           than MPU.

       "primes"
           Doesn't support ranges, requires bumping up the precalculated primes for larger
           numbers, which means knowing in advance the upper limit for primes.  Support for
           numbers larger than 400M requires using Pari version 2.3.5.  If that is used, sieving
           is about 2x faster than MPU, but doesn't support segmenting.

       "factorint"
           Similar to MPU's "factor_exp" though with a slightly different return.  MPU offers
           "factor" for a linear array of prime factors where
              n = p1 * p2 * p3 * ...   as (p1,p2,p3,...)  and "factor_exp" for an array of
           factor/exponent pairs where:
              n = p1^e1 * p2^e2 * ...  as ([p1,e1],[p2,e2],...)  Pari/GP returns an array similar
           to the latter.  Math::Pari returns a transposed matrix like:
              n = p1^e1 * p2^e2 * ...  as ([p1,p2,...],[e1,e2,...])  Slower than MPU for all
           64-bit inputs on an x86_64 platform, it may be faster for large values on other
           platforms.  With the newer Math::Prime::Util::GMP releases, bigint factoring is
           slightly faster on average in MPU.

       "divisors"
           Similar to MPU's "divisors".

       "forprime", "forcomposite", "fordiv", "sumdiv"
           Similar to MPU's "forprimes", "forcomposites", "fordivisors", and "divisor_sum".

       "eulerphi", "moebius"
           Similar to MPU's "euler_phi" and "moebius".  MPU is 2-20x faster for native integers.
           MPU also supported range inputs, which can be much more efficient.  With bigint
           arguments, MPU is slightly faster than Math::Pari if the GMP backend is available, but
           very slow without.

       "gcd", "lcm", "kronecker", "znorder", "znprimroot", "znlog"
           Similar to MPU's "gcd", "lcm", "kronecker", "znorder", "znprimroot", and "znlog".
           Pari's "znprimroot" only returns the smallest root for prime powers.  The behavior is
           undefined when the group is not cyclic (sometimes it throws an exception, sometimes it
           returns an incorrect answer, sometimes it hangs).  MPU's "znprimroot" will always
           return the smallest root if it exists, and "undef" otherwise.  Similarly, MPU's
           "znlog" will return the smallest "x" and work with non-primitive-root "g", which is
           similar to Pari/GP 2.6, but not the older versions in Math::Pari.  The performance of
           "znlog" is quite good compared to older Pari/GP, but much worse than 2.6's new
           methods.

       "sigma"
           Similar to MPU's "divisor_sum".  MPU is ~10x faster when the result fits in a native
           integer.  Once things overflow it is fairly similar in performance.  However, using
           Math::BigInt can slow things down quite a bit, so for best performance in these cases
           using a Math::GMP object is best.

       "numbpart", "forpart"
           Similar to MPU's "partitions" and "forpart".  These functions were introduced in Pari
           2.3 and 2.6, hence are not in Math::Pari.  "numbpart" produce identical results to
           "partitions", but Pari is much faster.  forpart is very similar to Pari's function,
           but produces a different ordering (MPU is the standard anti-lexicographical, Pari uses
           a size sort).  Currently Pari is somewhat faster due to Perl function call overhead.
           When using restrictions, Pari has much better optimizations.

       "eint1"
           Similar to MPU's "ExponentialIntegral".

       "zeta"
           MPU has "RiemannZeta" which takes non-negative real inputs, while Pari's function
           supports negative and complex inputs.

       Overall, Math::Pari supports a huge variety of functionality and has a sophisticated and
       mature code base behind it (noting that the Pari library used is about 10 years old now).
       For native integers, typically Math::Pari will be slower than MPU.  For bigints,
       Math::Pari may be superior and it rarely has any performance surprises.  Some of the
       unique features MPU offers include super fast prime counts, nth_prime, ECPP primality
       proofs with certificates, approximations and limits for both, random primes, fast Mertens
       calculations, Chebyshev theta and psi functions, and the logarithmic integral and Riemann
       R functions.  All with fairly minimal installation requirements.

PERFORMANCE

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

       Primality testing
           For general numbers smaller than 2000 or so digits, MPU is the fastest solution I am
           aware of (it is faster than Pari 2.7, PFGW, and FLINT).  For very large inputs, PFGW
           <http://sourceforge.net/projects/openpfgw/> is the fastest primality testing software
           I'm aware of.  It has fast trial division, and is especially fast on many special
           forms.  It does not have a BPSW test however, and there are quite a few
           counterexamples for a given base of its PRP test, so it is commonly used for fast
           filtering of large candidates.  A test such as the BPSW test in this module is then
           recommended.

       Primality proofs
           Primo <http://www.ellipsa.eu/> is the best method for open source primality proving
           for inputs over 1000 digits.  Primo also does well below that size, but other good
           alternatives are David Cleaver's mpzaprcl <http://sourceforge.net/projects/mpzaprcl/>,
           the APRCL from the modern Pari <http://pari.math.u-bordeaux.fr/> package, or the
           standalone ECPP from this module with large polynomial set.

       Factoring
           yafu <http://sourceforge.net/projects/yafu/>, msieve
           <http://sourceforge.net/projects/msieve/>, and gmp-ecm <http://ecm.gforge.inria.fr/>
           are all good choices for large inputs.  The factoring code in this module (and all
           other CPAN modules) is very limited compared to those.

       Primes
           primesieve <http://code.google.com/p/primesieve/> and yafu
           <http://sourceforge.net/projects/yafu/> are the fastest publically available code I am
           aware of.  Primesieve will additionally take advantage of multiple cores with
           excellent efficiency.  Tomás Oliveira e Silva's private code may be faster for very
           large values, but isn't available for testing.

           Note that the Sieve of Atkin is not faster than the Sieve of Eratosthenes when both
           are well implemented.  The only Sieve of Atkin that is even competitive is Bernstein's
           super optimized primegen, which runs on par with the SoE in this module.  The SoE's in
           Pari, yafu, and primesieve are all faster.

       Prime Counts and Nth Prime
           Outside of private research implementations doing prime counts for "n > 2^64", this
           module should be close to state of the art in performance, and supports results up to
           "2^64".  Further performance improvements are planned, as well as expansion to larger
           values.

           The fastest solution for small inputs is a hybrid table/sieve method.  This module
           does this for values below 60M.  As the inputs get larger, either the tables have to
           grow exponentially or speed must be sacrificed.  Hence this is not a good general
           solution for most uses.

   PRIME COUNTS
       Counting the primes to "800_000_000" (800 million):

         Time (s)   Module                      Version  Notes
         ---------  --------------------------  -------  -----------
              0.001 Math::Prime::Util           0.37     using extended LMO
              0.007 Math::Prime::Util           0.12     using Lehmer's method
              0.27  Math::Prime::Util           0.17     segmented 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
       Factoring performance depends on the input, and the algorithm choices used are still being
       tuned.  Math::Factor::XS is very fast when given input with only small factors, but it
       slows down rapidly as the smallest factor increases in size.  For numbers larger than 32
       bits, Math::Prime::Util can be 100x or more faster (a number with only very small factors
       will be nearly identical, while a semiprime may be 3000x faster).  Math::Pari is much
       slower with native sized inputs, probably due to calling overhead.  For bigints, the
       Math::Prime::Util::GMP module is needed or performance will be far worse than Math::Pari.
       With the GMP module, performance is pretty similar from 20 through 70 digits, which the
       caveat that the current MPU factoring uses more memory for 60+ digit numbers.

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

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

   RANDOM PRIME GENERATION
       Seconds per prime for random prime generation on a early 2015 Macbook Pro (2.7 GHz i5)
       with Math::BigInt::GMP and Math::Prime::Util::GMP installed.

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

AUTHORS

       Dana Jacobsen <dana@acm.org>

ACKNOWLEDGEMENTS

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

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

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

       The LMO implementation is based on the 2003 preprint from Christian Bau, as well as the
       2006 paper from Tomás Oliveira e Silva.  I also want to thank Kim Walisch for the many
       discussions about prime counting.

REFERENCES

       •   Christian Axler, "New bounds for the prime counting function π(x)", September 2014.
           For large values, improved limits versus Dusart 2010.
           <http://arxiv.org/abs/1409.1780>

       •   Christian Axler, "Über die Primzahl-Zählfunktion, die n-te Primzahl und
           verallgemeinerte Ramanujan-Primzahlen", January 2013.  Prime count and nth-prime
           bounds in more detail.  Thesis in German, but first part is easily read.
           <http://docserv.uni-duesseldorf.de/servlets/DerivateServlet/Derivate-28284/pdfa-1b.pdf>

       •   Christian Bau, "The Extended Meissel-Lehmer Algorithm", 2003, preprint with example
           C++ implementation.  Very detailed implementation-specific paper which was used for
           the implementation here.  Highly recommended for implementing a sieve-based LMO.
           <http://cs.swan.ac.uk/~csoliver/ok-sat-library/OKplatform/ExternalSources/sources/NumberTheory/ChristianBau/>

       •   Manuel Benito and Juan L. Varona, "Recursive formulas related to the summation of the
           Möbius function", The Open Mathematics Journal, v1, pp 25-34, 2007.  Among many other
           things, shows a simple formula for computing the Mertens functions with only n/3
           Möbius values (not as fast as Deléglise and Rivat, but really simple).
           <http://www.unirioja.es/cu/jvarona/downloads/Benito-Varona-TOMATJ-Mertens.pdf>

       •   John Brillhart, D. H. Lehmer, and J. L. Selfridge, "New Primality Criteria and
           Factorizations of 2^m +/- 1", Mathematics of Computation, v29, n130, Apr 1975, pp
           620-647.
           <http://www.ams.org/journals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf>

       •   W. J. Cody and Henry C. Thacher, Jr., "Rational Chebyshev Approximations for the
           Exponential Integral E_1(x)", Mathematics of Computation, v22, pp 641-649, 1968.

       •   W. J. Cody and Henry C. Thacher, Jr., "Chebyshev approximations for the exponential
           integral Ei(x)", Mathematics of Computation, v23, pp 289-303, 1969.
           <http://www.ams.org/journals/mcom/1969-23-106/S0025-5718-1969-0242349-2/>

       •   W. J. Cody, K. E. Hillstrom, and Henry C. Thacher Jr., "Chebyshev Approximations for
           the Riemann Zeta Function", "Mathematics of Computation", v25, n115, pp 537-547, July
           1971.

       •   Henri Cohen, "A Course in Computational Algebraic Number Theory", Springer, 1996.
           Practical computational number theory from the team lead of Pari <http://pari.math.u-
           bordeaux.fr/>.  Lots of explicit algorithms.

       •   Marc Deléglise and Joöl Rivat, "Computing the summation of the Möbius function",
           Experimental Mathematics, v5, n4, pp 291-295, 1996.  Enhances the Möbius computation
           in Lioen/van de Lune, and gives a very efficient way to compute the Mertens function.
           <http://projecteuclid.org/euclid.em/1047565447>

       •   Pierre Dusart, "Autour de la fonction qui compte le nombre de nombres premiers", PhD
           thesis, 1998.  In French.  The mathematics is readable and highly recommended reading
           if you're interested in prime number bounds.
           <http://www.unilim.fr/laco/theses/1998/T1998_01.html>

       •   Pierre Dusart, "Estimates of Some Functions Over Primes without R.H.", preprint, 2010.
           Updates to the best non-RH bounds for prime count and nth prime.
           <http://arxiv.org/abs/1002.0442/>

       •   Pierre-Alain Fouque and Mehdi Tibouchi, "Close to Uniform Prime Number Generation With
           Fewer Random Bits", pre-print, 2011.  Describes random prime distributions, their
           algorithm for creating random primes using few random bits, and comparisons to other
           methods.  Definitely worth reading for the discussions of uniformity.
           <http://eprint.iacr.org/2011/481>

       •   Walter M. Lioen and Jan van de Lune, "Systematic Computations on Mertens' Conjecture
           and Dirichlet's Divisor Problem by Vectorized Sieving", in From Universal Morphisms to
           Megabytes, Centrum voor Wiskunde en Informatica, pp. 421-432, 1994.  Describes a nice
           way to compute a range of Möbius values.  <http://walter.lioen.com/papers/LL94.pdf>

       •   Ueli M. Maurer, "Fast Generation of Prime Numbers and Secure Public-Key Cryptographic
           Parameters", 1995.  Generating random provable primes by building up the prime.
           <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.26.2151>

       •   Gabriel Mincu, "An Asymptotic Expansion", Journal of Inequalities in Pure and Applied
           Mathematics, v4, n2, 2003.  A very readable account of Cipolla's 1902 nth prime
           approximation.  <http://www.emis.de/journals/JIPAM/images/153_02_JIPAM/153_02.pdf>

       •   OEIS: Primorial <http://oeis.org/wiki/Primorial>

       •   Vincent Pegoraro and Philipp Slusallek, "On the Evaluation of the Complex-Valued
           Exponential Integral", Journal of Graphics, GPU, and Game Tools, v15, n3, pp 183-198,
           2011.  <http://www.cs.utah.edu/~vpegorar/research/2011_JGT/paper.pdf>

       •   William H. Press et al., "Numerical Recipes", 3rd edition.

       •   Hans Riesel, "Prime Numbers and Computer Methods for Factorization", Birkh?user, 2nd
           edition, 1994.  Lots of information, some code, easy to follow.

       •   David M. Smith, "Multiple-Precision Exponential Integral and Related Functions", ACM
           Transactions on Mathematical Software, v37, n4, 2011.
           <http://myweb.lmu.edu/dmsmith/toms2011.pdf>

       •   Douglas A. Stoll and Patrick Demichel , "The impact of ζ(s) complex zeros on π(x) for
           x < 10^{10^{13}}", "Mathematics of Computation", v80, n276, pp 2381-2394, October
           2011.
           <http://www.ams.org/journals/mcom/2011-80-276/S0025-5718-2011-02477-4/home.html>

COPYRIGHT

       Copyright 2011-2018 by Dana Jacobsen <dana@acm.org>

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