Skip site navigation (1)Skip section navigation (2)

FreeBSD Manual Pages

  
 
  

home | help
Math::Prime::Util(3)  User Contributed Perl Documentation Math::Prime::Util(3)

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:

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

       o   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:

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

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

       o   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),
       BA1/4the	(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 BA(C)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! 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)	a! 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 I1/4(n),	the MA<paragraph>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 MA<paragraph>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
       DelA(C)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 DelA(C)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
       I1/4(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.  DelA(C)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 I(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(I(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 I>>(n), the Liouville function for a non-negative integer
       input.  This is -1 raised to I(C)(n) (the total number of prime
       factors).

   chebyshev_theta
	 say chebyshev_theta(10000);

       Returns I,(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 I(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 I>>(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 a! 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

       o   p is	large.	 This function requires	at least 128 bits.

       o   p-1 has a large prime factor	r.

       o   p+1 has a large prime factor	s

       o   r-1 has a large prime factor	t

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

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

   random_proven_prime
	 my $bigprime =	random_proven_prime(512);

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

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

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

   random_maurer_prime
	 my $bigprime =	random_maurer_prime(512);

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

       The performance with Math::Prime::Util::GMP installed is	hundreds of
       times faster, so	it is highly recommended.

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

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

       If you don't need absolutely proven results, then it is somewhat	faster
       to use "random_nbit_prime" either by itself or with some	additional
       tests, e.g.  "miller_rabin_random" and/or
       "is_frobenius_underwood_pseudoprime".  One could	also run
       is_provable_prime on the	result,	but this will be slow.

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

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

   random_shawe_taylor_prime
	 my $bigprime =	random_shawe_taylor_prime(8192);

       Construct an n-bit provable prime, using	the Shawe-Taylor algorithm in
       section C.6 of FIPS 186-4.  This	uses 512 bits of randomness and
       SHA-256 as the hash.  This is a slightly	simpler	and older (1986)
       method than Maurer's 1999 construction.	It is a	bit faster than
       Maurer's	method,	and uses less system entropy for large sizes.  The
       primary reason to use this rather than Maurer's method is to use	the
       FIPS 186-4 algorithm.

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

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

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

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

   random_semiprime
       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 I(C)(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 I(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 I<paragraph>(s)-1, where I<paragraph>(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 I.3(x) function of	Planat and SolA(C) (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]A.$_->[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:

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

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

       o   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.  TomA!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)

       o   MPU is consistently the fastest solution, and performs the most
	   stringent probable prime tests on bigints.

       o   Math::Primality has a lot of	overhead that makes it quite slow for
	   native size integers.  With bigints we finally see it work well.

       o   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 TomA!s Oliveira e Silva.  I also
       want to thank Kim Walisch for the many discussions about	prime
       counting.

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

       o   Christian Axler, "Aber die Primzahl-ZAxhlfunktion, 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>

       o   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/>

       o   Manuel Benito and Juan L. Varona, "Recursive	formulas related to
	   the summation of the	MA<paragraph>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 MA<paragraph>bius values (not as fast as DelA(C)glise and
	   Rivat, but really simple).
	   <http://www.unirioja.es/cu/jvarona/downloads/Benito-Varona-TOMATJ-Mertens.pdf>

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

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

       o   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/>

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

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

       o   Marc	DelA(C)glise and JoA<paragraph>l Rivat,	"Computing the
	   summation of	the MA<paragraph>bius function", Experimental
	   Mathematics,	v5, n4,	pp 291-295, 1996.  Enhances the
	   MA<paragraph>bius computation in Lioen/van de Lune, and gives a
	   very	efficient way to compute the Mertens function.
	   <http://projecteuclid.org/euclid.em/1047565447>

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

       o   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/>

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

       o   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 MA<paragraph>bius values.
	   <http://walter.lioen.com/papers/LL94.pdf>

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

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

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

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

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

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

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

       o   Douglas A. Stoll and	Patrick	Demichel , "The	impact of
	   I<paragraph>(s) complex zeros on I(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.

perl v5.32.1			  2018-11-15		  Math::Prime::Util(3)

NAME | VERSION | SYNOPSIS | DESCRIPTION | ENVIRONMENT VARIABLES | BIGNUM SUPPORT | PRIMALITY TESTING | FUNCTIONS | RANDOM NUMBERS | RANDOM PRIMES | UTILITY FUNCTIONS | FACTORING FUNCTIONS | MATHEMATICAL FUNCTIONS | EXAMPLES | PRIMALITY TESTING NOTES | LIMITATIONS | SEE ALSO | PERFORMANCE | AUTHORS | ACKNOWLEDGEMENTS | REFERENCES | COPYRIGHT

Want to link to this manual page? Use this URL:
<https://www.freebsd.org/cgi/man.cgi?query=Math::Prime::Util&sektion=3&manpath=FreeBSD+13.0-RELEASE+and+Ports>

home | help