.\" Automatically generated by Pod::Man 2.28 (Pod::Simple 3.28) .\" .\" Standard preamble: .\" ======================================================================== .de Sp \" Vertical space (when we can't use .PP) .if t .sp .5v .if n .sp .. .de Vb \" Begin verbatim text .ft CW .nf .ne \\$1 .. .de Ve \" End verbatim text .ft R .fi .. .\" Set up some character translations and predefined strings. \*(-- will .\" give an unbreakable dash, \*(PI will give pi, \*(L" will give a left .\" double quote, and \*(R" will give a right double quote. \*(C+ will .\" give a nicer C++. Capital omega is used to do unbreakable dashes and .\" therefore won't be available. \*(C` and \*(C' expand to `' in nroff, .\" nothing in troff, for use with C<>. .tr \(*W- .ds C+ C\v'-.1v'\h'-1p'\s-2+\h'-1p'+\s0\v'.1v'\h'-1p' .ie n \{\ . ds -- \(*W- . ds PI pi . if (\n(.H=4u)&(1m=24u) .ds -- \(*W\h'-12u'\(*W\h'-12u'-\" diablo 10 pitch . if (\n(.H=4u)&(1m=20u) .ds -- \(*W\h'-12u'\(*W\h'-8u'-\" diablo 12 pitch . ds L" "" . ds R" "" . ds C` "" . ds C' "" 'br\} .el\{\ . ds -- \|\(em\| . ds PI \(*p . ds L" `` . ds R" '' . ds C` . ds C' 'br\} .\" .\" Escape single quotes in literal strings from groff's Unicode transform. .ie \n(.g .ds Aq \(aq .el .ds Aq ' .\" .\" If the F register is turned on, we'll generate index entries on stderr for .\" titles (.TH), headers (.SH), subsections (.SS), items (.Ip), and index .\" entries marked with X<> in POD. Of course, you'll have to process the .\" output yourself in some meaningful fashion. .\" .\" Avoid warning from groff about undefined register 'F'. .de IX .. .nr rF 0 .if \n(.g .if rF .nr rF 1 .if (\n(rF:(\n(.g==0)) \{ . if \nF \{ . de IX . tm Index:\\$1\t\\n%\t"\\$2" .. . if !\nF==2 \{ . nr % 0 . nr F 2 . \} . \} .\} .rr rF .\" .\" Accent mark definitions (@(#)ms.acc 1.5 88/02/08 SMI; from UCB 4.2). .\" Fear. Run. Save yourself. No user-serviceable parts. . \" fudge factors for nroff and troff .if n \{\ . ds #H 0 . ds #V .8m . ds #F .3m . ds #[ \f1 . ds #] \fP .\} .if t \{\ . ds #H ((1u-(\\\\n(.fu%2u))*.13m) . ds #V .6m . ds #F 0 . ds #[ \& . ds #] \& .\} . \" simple accents for nroff and troff .if n \{\ . ds ' \& . ds ` \& . ds ^ \& . ds , \& . ds ~ ~ . ds / .\} .if t \{\ . ds ' \\k:\h'-(\\n(.wu*8/10-\*(#H)'\'\h"|\\n:u" . ds ` \\k:\h'-(\\n(.wu*8/10-\*(#H)'\`\h'|\\n:u' . ds ^ \\k:\h'-(\\n(.wu*10/11-\*(#H)'^\h'|\\n:u' . ds , \\k:\h'-(\\n(.wu*8/10)',\h'|\\n:u' . ds ~ \\k:\h'-(\\n(.wu-\*(#H-.1m)'~\h'|\\n:u' . ds / \\k:\h'-(\\n(.wu*8/10-\*(#H)'\z\(sl\h'|\\n:u' .\} . \" troff and (daisy-wheel) nroff accents .ds : \\k:\h'-(\\n(.wu*8/10-\*(#H+.1m+\*(#F)'\v'-\*(#V'\z.\h'.2m+\*(#F'.\h'|\\n:u'\v'\*(#V' .ds 8 \h'\*(#H'\(*b\h'-\*(#H' .ds o \\k:\h'-(\\n(.wu+\w'\(de'u-\*(#H)/2u'\v'-.3n'\*(#[\z\(de\v'.3n'\h'|\\n:u'\*(#] .ds d- \h'\*(#H'\(pd\h'-\w'~'u'\v'-.25m'\f2\(hy\fP\v'.25m'\h'-\*(#H' .ds D- D\\k:\h'-\w'D'u'\v'-.11m'\z\(hy\v'.11m'\h'|\\n:u' .ds th \*(#[\v'.3m'\s+1I\s-1\v'-.3m'\h'-(\w'I'u*2/3)'\s-1o\s+1\*(#] .ds Th \*(#[\s+2I\s-2\h'-\w'I'u*3/5'\v'-.3m'o\v'.3m'\*(#] .ds ae a\h'-(\w'a'u*4/10)'e .ds Ae A\h'-(\w'A'u*4/10)'E . \" corrections for vroff .if v .ds ~ \\k:\h'-(\\n(.wu*9/10-\*(#H)'\s-2\u~\d\s+2\h'|\\n:u' .if v .ds ^ \\k:\h'-(\\n(.wu*10/11-\*(#H)'\v'-.4m'^\v'.4m'\h'|\\n:u' . \" for low resolution devices (crt and lpr) .if \n(.H>23 .if \n(.V>19 \ \{\ . ds : e . ds 8 ss . ds o a . ds d- d\h'-1'\(ga . ds D- D\h'-1'\(hy . ds th \o'bp' . ds Th \o'LP' . ds ae ae . ds Ae AE .\} .rm #[ #] #H #V #F C .\" ======================================================================== .\" .IX Title "Math::Prime::Util 3pm" .TH Math::Prime::Util 3pm "2014-10-21" "perl v5.20.1" "User Contributed Perl Documentation" .\" For nroff, turn off justification. Always turn off hyphenation; it makes .\" way too many mistakes in technical documents. .if n .ad l .nh .SH "NAME" Math::Prime::Util \- Utilities related to prime numbers, including fast sieves and factoring .SH "VERSION" .IX Header "VERSION" Version 0.46 .SH "SYNOPSIS" .IX Header "SYNOPSIS" .Vb 3 \& # Normally you would just import the functions you are using. \& # Nothing is exported by default. List the functions, or use :all. \& use Math::Prime::Util \*(Aq:all\*(Aq; \& \& \& # 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\*(Aqd overflow) \& $n = next_prime($n); \& \& # step back (returns 0 if given input less than 2) \& $n = prev_prime($n); \& \& \& # Return Pi(n) \-\- the number of primes E= 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\*(Aqs totient) on a large number \& use bigint; say euler_phi( 801294088771394680000412 ); \& say jordan_totient(5, 1234); # Jordan\*(Aqs 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 .Ve .SH "DESCRIPTION" .IX Header "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. .PP This module is the fastest on \s-1CPAN\s0 for almost all operations. Only Math::Pari is faster for a few operations. This includes Math::Prime::XS, Math::Prime::FastSieve, Math::Factor::XS, Math::Prime::TiedArray, Math::Big::Factors, Math::Factoring, and Math::Primality (when the \s-1GMP\s0 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. .PP All operations support both Perl \s-1UV\s0'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. .PP 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. .PP Two scripts are also included and installed by default: .IP "\(bu" 4 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 \f(CW\*(C`\-\-help\*(C'\fR to see all the options. .IP "\(bu" 4 factor.pl operates similar to the \s-1GNU \s0\f(CW\*(C`factor\*(C'\fR program. It supports bigint and expression inputs. .SH "BIGNUM SUPPORT" .IX Header "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. .PP If you are using bigints, here are some performance suggestions: .IP "\(bu" 4 Install Math::Prime::Util::GMP, as that will vastly increase the speed of many of the functions. This does require the \s-1GMP\s0 library be installed on your system, but this increasingly comes pre-installed or easily available using the \s-1OS\s0 vendor package installation tool. .IP "\(bu" 4 Install and use Math::BigInt::GMP or Math::BigInt::Pari, then use \&\f(CW\*(C`use bigint try => \*(AqGMP,Pari\*(Aq\*(C'\fR in your script, or on the command line \&\f(CW\*(C`\-Mbigint=lib,GMP\*(C'\fR. Large modular exponentiation is much faster using the \&\s-1GMP\s0 or Pari backends, as are the math and approximation functions when called with very large inputs. .IP "\(bu" 4 Install Math::MPFR if you use the Ei, li, Zeta, or R functions. If that module can be loaded, these functions will run much faster on bignum inputs, and are able to provide higher accuracy. .IP "\(bu" 4 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. .SH "PRIMALITY TESTING" .IX Header "PRIMALITY TESTING" This module provides three functions for general primality testing, as well as numerous specialized functions. The three main functions are: \&\*(L"is_prob_prime\*(R" and \*(L"is_prime\*(R" for general use, and \*(L"is_provable_prime\*(R" for proofs. For inputs below \f(CW\*(C`2^64\*(C'\fR 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 \s-1CPAN\s0 modules. For inputs larger than \f(CW\*(C`2^64\*(C'\fR, an extra-strong \&\s-1BPSW\s0 test is used. See the \*(L"\s-1PRIMALITY TESTING NOTES\*(R"\s0 section for more discussion. .SH "FUNCTIONS" .IX Header "FUNCTIONS" .SS "is_prime" .IX Subsection "is_prime" .Vb 1 \& print "$n is prime" if is_prime($n); .Ve .PP Returns 0 is the number is composite, 1 if it is probably prime, and 2 if it is definitely prime. For numbers smaller than \f(CW\*(C`2^64\*(C'\fR 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 \s-1BPSW\s0 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. .PP Also see the \*(L"is_prob_prime\*(R" function, which will never do additional tests, and the \*(L"is_provable_prime\*(R" function which will construct a proof that the input is number prime and returns 2 for almost all primes (at the expense of speed). .PP For native precision numbers (anything smaller than \f(CW\*(C`2^64\*(C'\fR, all three functions are identical and use a deterministic set of tests (selected Miller-Rabin bases or \s-1BPSW\s0). For larger inputs both \*(L"is_prob_prime\*(R" and \&\*(L"is_prime\*(R" return probable prime results using the extra-strong Baillie-PSW test, which has had no counterexample found since it was published in 1980. .PP For cryptographic key generation, you may want even more testing for probable primes (\s-1NIST\s0 recommends some additional M\-R tests). This can be done using a different test (e.g. \*(L"is_frobenius_underwood_pseudoprime\*(R") or using additional M\-R tests with random bases with \*(L"miller_rabin_random\*(R". Even better, make sure Math::Prime::Util::GMP is installed and use \&\*(L"is_provable_prime\*(R" which should be reasonably fast for sizes under 2048 bits. Another possibility is to use \&\*(L"random_maurer_prime\*(R" in Math::Prime::Util or \&\*(L"random_shawe_taylor_prime\*(R" in Math::Prime::Util which construct random provable primes. .SS "primes" .IX Subsection "primes" Returns all the primes between the lower and upper limits (inclusive), with a lower limit of \f(CW2\fR if none is given. .PP An array reference is returned (with large lists this is much faster and uses less memory than returning an array directly). .PP .Vb 2 \& my $aref1 = primes( 1_000_000 ); \& my $aref2 = primes( 1_000_000_000_000, 1_000_000_001_000 ); \& \& my @primes = @{ primes( 500 ) }; \& \& print "$_\en" for @{primes(20,100)}; .Ve .PP 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. .SS "next_prime" .IX Subsection "next_prime" .Vb 1 \& $n = next_prime($n); .Ve .PP 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 \f(CW\*(C`4,294,967,291\*(C'\fR in 32\-bit Perl; larger than \f(CW\*(C`18,446,744,073,709,551,557\*(C'\fR in 64\-bit). .SS "prev_prime" .IX Subsection "prev_prime" .Vb 1 \& $n = prev_prime($n); .Ve .PP Returns the prime preceding the input number (i.e. the largest prime that is strictly less than the input). 0 is returned if the input is \f(CW2\fR or lower. .SS "forprimes" .IX Subsection "forprimes" .Vb 1 \& 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 .Ve .PP 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. .PP Math::BigInt objects may be used for the range. .PP For some uses an iterator (\*(L"prime_iterator\*(R", \*(L"prime_iterator_object\*(R") or a tied array (Math::Prime::Util::PrimeArray) may be more convenient. Objects can be passed to functions, and allow early loop exits. .SS "forcomposites" .IX Subsection "forcomposites" .Vb 2 \& forcomposites { say } 1000; \& forcomposites { say } 2000,2020; .Ve .PP 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, \&\s-1OEIS A002808\s0 , are the numbers greater than 1 which are not prime: \f(CW\*(C`4, 6, 8, 9, 10, 12, 14, 15, ...\*(C'\fR. .SS "foroddcomposites" .IX Subsection "foroddcomposites" Similar to \*(L"forcomposites\*(R", but skipping all even numbers. The odd composites, \s-1OEIS A071904\s0 , are the numbers greater than 1 which are not prime and not divisible by two: \&\f(CW\*(C`9, 15, 21, 25, 27, 33, 35, ...\*(C'\fR. .SS "fordivisors" .IX Subsection "fordivisors" .Vb 1 \& fordivisors { $prod *= $_ } $n; .Ve .PP Given a block and a non-negative number \f(CW\*(C`n\*(C'\fR, the block is called with \&\f(CW$_\fR set to each divisor in sorted order. Also see \*(L"divisor_sum\*(R". .SS "forpart" .IX Subsection "forpart" .Vb 3 \& forpart { say "@_" } 25; # unrestricted partitions \& forpart { say "@_" } 25,{n=>5} # ... with exactly 5 values \& forpart { say "@_" } 25,{nmax=>5} # ... with <=5 values .Ve .PP Given a non-negative number \f(CW\*(C`n\*(C'\fR, the block is called with \f(CW@_\fR set to the array of additive integer partitions. The operation is very similar to the \f(CW\*(C`forpart\*(C'\fR function in Pari/GP 2.6.x, though the ordering is different. The algorithm is \s-1ZS1\s0 from Zoghbi and StojmenoviX (1998), hence the ordering is identical to that of Integer::Partition. Use \*(L"partitions\*(R" to get just the count of unrestricted partitions. .PP An optional hash reference may be given to produce restricted partitions. Each value must be a non-negative integer. The allowable keys are: .PP .Vb 5 \& 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 .Ve .PP Like forcomb and forperm, the partition return values are read-only. Any attempt to modify them will result in undefined behavior. .SS "forcomb" .IX Subsection "forcomb" Given non-negative arguments \f(CW\*(C`n\*(C'\fR and \f(CW\*(C`k\*(C'\fR, the block is called with \f(CW@_\fR set to the \f(CW\*(C`k\*(C'\fR element array of values from \f(CW0\fR to \f(CW\*(C`n\-1\*(C'\fR representing the combinations in lexicographical order. While the binomial function gives the total number, this function can be used to enumerate the choices. .PP Rather than give a data array as input, an integer is used for \f(CW\*(C`n\*(C'\fR. A convenient way to map to array elements is: .PP .Vb 1 \& forcomb { say "@data[@_]" } @data, 3; .Ve .PP where the block maps the combination array \f(CW@_\fR to array values, the argument for \f(CW\*(C`n\*(C'\fR is given the array since it will be evaluated as a scalar and hence give the size, and the argument for \f(CW\*(C`k\*(C'\fR is the desired size of the combinations. .PP Like forpart and forperm, the index return values are read-only. Any attempt to modify them will result in undefined behavior. .SS "forperm" .IX Subsection "forperm" Given non-negative argument \f(CW\*(C`n\*(C'\fR, the block is called with \f(CW@_\fR set to the \f(CW\*(C`k\*(C'\fR element array of values from \f(CW0\fR to \f(CW\*(C`n\-1\*(C'\fR representing permutations in lexicographical order. The total number of calls will be \f(CW\*(C`n!\*(C'\fR. .PP Rather than give a data array as input, an integer is used for \f(CW\*(C`n\*(C'\fR. A convenient way to map to array elements is: .PP .Vb 1 \& forperm { say "@data[@_]" } @data; .Ve .PP where the block maps the permutation array \f(CW@_\fR to array values, and the argument for \f(CW\*(C`n\*(C'\fR is given the array since it will be evaluated as a scalar and hence give the size. .PP Like forpart and forcomb, the index return values are read-only. Any attempt to modify them will result in undefined behavior. .SS "prime_iterator" .IX Subsection "prime_iterator" .Vb 2 \& my $it = prime_iterator; \& $sum += $it\->() for 1..100000; .Ve .PP 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: .PP .Vb 1 \& my $it = prime_iterator(200); say $it\->(); say $it\->(); .Ve .PP 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. .PP Other options include \*(L"forprimes\*(R" (more efficiency, less flexibility), Math::Prime::Util::PrimeIterator (an iterator with more functionality), or Math::Prime::Util::PrimeArray (a tied array). .SS "prime_iterator_object" .IX Subsection "prime_iterator_object" .Vb 3 \& my $it = prime_iterator_object; \& while ($it\->value < 100) { say $it\->value; $it\->next; } \& $sum += $it\->iterate for 1..100000; .Ve .PP 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. .SS "prime_count" .IX Subsection "prime_count" .Vb 2 \& my $primepi = prime_count( 1_000 ); \& my $pirange = prime_count( 1_000, 10_000 ); .Ve .PP Returns the Prime Count function \f(CWPi(n)\fR, also called \f(CW\*(C`primepi\*(C'\fR in some math packages. When given two arguments, it returns the inclusive count of primes between the ranges. E.g. \f(CW\*(C`(13,17)\*(C'\fR returns 2, \f(CW\*(C`(14,17)\*(C'\fR and \f(CW\*(C`(13,16)\*(C'\fR return 1, \f(CW\*(C`(14,16)\*(C'\fR returns 0. .PP The current implementation decides based on the ranges whether to use a segmented sieve with a fast bit count, or the extended \s-1LMO\s0 algorithm. The former is preferred for small sizes as well as small ranges. The latter is much faster for large ranges. .PP The segmented sieve is very memory efficient and is quite fast even with large base values. Its complexity is approximately \f(CW\*(C`O(sqrt(a) + (b\-a))\*(C'\fR, where the first term is typically negligible below \f(CW\*(C`~ 10^11\*(C'\fR. Memory use is proportional only to \f(CWsqrt(a)\fR, with total memory use under 1MB for any base under \f(CW\*(C`10^14\*(C'\fR. .PP The extended \s-1LMO\s0 method has complexity approximately \&\f(CW\*(C`O(b^(2/3)) + O(a^(2/3))\*(C'\fR, and also uses low memory. A calculation of \f(CW\*(C`Pi(10^14)\*(C'\fR completes in a few seconds, \f(CW\*(C`Pi(10^15)\*(C'\fR in well under a minute, and \f(CW\*(C`Pi(10^16)\*(C'\fR in about one minute. In contrast, even parallel primesieve would take over a week on a similar machine to determine \f(CW\*(C`Pi(10^16)\*(C'\fR. .PP Also see the function \*(L"prime_count_approx\*(R" which gives a very good approximation to the prime count, and \*(L"prime_count_lower\*(R" and \&\*(L"prime_count_upper\*(R" which give tight bounds to the actual prime count. These functions return quickly for any input, including bigints. .SS "prime_count_upper" .IX Subsection "prime_count_upper" .SS "prime_count_lower" .IX Subsection "prime_count_lower" .Vb 3 \& my $lower_limit = prime_count_lower($n); \& my $upper_limit = prime_count_upper($n); \& # $lower_limit <= prime_count(n) <= $upper_limit .Ve .PP 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 \f(CW\*(C`prime_count\*(C'\fR will always be equal to or between these numbers. .PP A common place these would be used is sizing an array to hold the first \f(CW$n\fR primes. It may be desirable to use a bit more memory than is necessary, to avoid calling \f(CW\*(C`prime_count\*(C'\fR. .PP These routines use verified tight limits below a range at least \f(CW\*(C`2^35\*(C'\fR, and use the Dusart (2010) bounds of .PP .Vb 1 \& x/logx * (1 + 1/logx + 2.000/log^2x) <= Pi(x) \& \& x/logx * (1 + 1/logx + 2.334/log^2x) >= Pi(x) .Ve .PP above that range. These bounds do not assume the Riemann Hypothesis. If the configuration option \f(CW\*(C`assume_rh\*(C'\fR has been set (it is off by default), then the Schoenfeld (1976) bounds are used for large values. .SS "prime_count_approx" .IX Subsection "prime_count_approx" .Vb 3 \& print "there are about ", \& prime_count_approx( 10 ** 18 ), \& " primes below one quintillion.\en"; .Ve .PP Returns an approximation to the \f(CW\*(C`prime_count\*(C'\fR function, without having to generate any primes. For values under \f(CW\*(C`10^36\*(C'\fR this uses the Riemann R function, which is quite accurate: an error of less than \f(CW\*(C`0.0005%\*(C'\fR is typical for input values over \f(CW\*(C`2^32\*(C'\fR, and decreases as the input gets larger. If Math::MPFR is installed, the Riemann R function is used for all values, and will be very fast. If not, then values of \f(CW\*(C`10^36\*(C'\fR and larger will use the approximation \f(CW\*(C`li(x) \- li(sqrt(x))/2\*(C'\fR. While not as accurate as the Riemann R function, it still should have error less than \f(CW\*(C`0.00000000000000001%\*(C'\fR. .PP A slightly faster but much less accurate answer can be obtained by averaging the upper and lower bounds. .SS "twin_primes" .IX Subsection "twin_primes" Returns the lesser of twin primes between the lower and upper limits (inclusive), with a lower limit of \f(CW2\fR if none is given. This is \&\s-1OEIS A001359\s0 . Given a twin prime pair \f(CW\*(C`(p,q)\*(C'\fR with \f(CW\*(C`q = p + 2\*(C'\fR, \f(CW\*(C`p prime\*(C'\fR, and , this function uses \f(CW\*(C`p\*(C'\fR to represent the pair. Hence the bounds need to include \f(CW\*(C`p\*(C'\fR, and the returned list will have \f(CW\*(C`p\*(C'\fR but not \f(CW\*(C`q\*(C'\fR. .PP This works just like the \*(L"primes\*(R" function, though only the first primes of twin prime pairs are returned. Like that function, an array reference is returned. .SS "twin_prime_count" .IX Subsection "twin_prime_count" Similar to prime count, but returns the count of twin primes (primes \f(CW\*(C`p\*(C'\fR where \f(CW\*(C`p+2\*(C'\fR is also prime). Takes either a single number indicating a count from 2 to the argument, or two numbers indicating a range. .PP The primes being counted are the first value, so a range of \f(CW\*(C`(3,5)\*(C'\fR will return a count of two, because both \f(CW3\fR and \f(CW5\fR are counted as twin primes. A range of \f(CW\*(C`(12,13)\*(C'\fR will return a count of zero, because neither \f(CW\*(C`12+2\*(C'\fR nor \f(CW\*(C`13+2\*(C'\fR are prime. In contrast, \f(CW\*(C`primesieve\*(C'\fR requires all elements of a constellation to be within the range to be counted, so would return one for the first example (\f(CW5\fR is not counted because its pair \f(CW7\fR is not in the range). .PP There is no useful formula known for this, unlike prime counts. We sieve for the answer, using some small table acceleration. .SS "twin_prime_count_approx" .IX Subsection "twin_prime_count_approx" Returns an approximation to the twin prime count of \f(CW\*(C`n\*(C'\fR. 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. .SS "nth_prime" .IX Subsection "nth_prime" .Vb 1 \& say "The ten thousandth prime is ", nth_prime(10_000); .Ve .PP Returns the prime that lies in index \f(CW\*(C`n\*(C'\fR in the array of prime numbers. Put another way, this returns the smallest \f(CW\*(C`p\*(C'\fR such that \f(CW\*(C`Pi(p) >= n\*(C'\fR. .PP 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. .PP 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 \f(CW\*(C`10^12th\*(C'\fR prime takes about 1 second, the \f(CW\*(C`10^13th\*(C'\fR prime takes under 10 seconds, and the \&\f(CW\*(C`10^14th\*(C'\fR prime (3475385758524527) takes under 30 seconds. Think about whether a bound or approximation would be acceptable, as they can be computed analytically. .PP 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. .SS "nth_prime_upper" .IX Subsection "nth_prime_upper" .SS "nth_prime_lower" .IX Subsection "nth_prime_lower" .Vb 3 \& my $lower_limit = nth_prime_lower($n); \& my $upper_limit = nth_prime_upper($n); \& # $lower_limit <= nth_prime(n) <= $upper_limit .Ve .PP Returns an analytical upper or lower bound on the Nth prime. These are very fast as they do not need to sieve or search through primes or tables. An exact answer is returned for tiny values of \f(CW\*(C`n\*(C'\fR. The lower limit uses the Dusart 2010 bound for all \f(CW\*(C`n\*(C'\fR, while the upper bound uses one of the two Dusart 2010 bounds for \f(CW\*(C`n >= 178974\*(C'\fR, a Dusart 1999 bound for \&\f(CW\*(C`n >= 39017\*(C'\fR, and a simple bound of \f(CW\*(C`n * (logn + 0.6 * loglogn)\*(C'\fR for small \f(CW\*(C`n\*(C'\fR. .SS "nth_prime_approx" .IX Subsection "nth_prime_approx" .Vb 1 \& say "The one trillionth prime is ~ ", nth_prime_approx(10**12); .Ve .PP Returns an approximation to the \f(CW\*(C`nth_prime\*(C'\fR function, without having to generate any primes. For values where the nth prime is smaller than \&\f(CW\*(C`2^64\*(C'\fR, an inverse Riemann R function is used. For larger values, uses the Cipolla 1902 approximation with up to 2nd order terms, plus a third order correction. .SS "nth_twin_prime" .IX Subsection "nth_twin_prime" Returns the Nth twin prime. This is done via sieving and counting, so is not very fast for large values. .SS "nth_twin_prime_approx" .IX Subsection "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. .SS "is_pseudoprime" .IX Subsection "is_pseudoprime" Takes a positive number \f(CW\*(C`n\*(C'\fR and one or more non-zero positive bases as input. Returns \f(CW1\fR if the input is a probable prime to each base, \f(CW0\fR if not. This is the simple Fermat primality test. Removing primes, given base 2 this produces the sequence \s-1OEIS A001567\s0 . .PP For practical use, \*(L"is_strong_pseudoprime\*(R" is a much stronger test with similar or better performance. .SS "is_strong_pseudoprime" .IX Subsection "is_strong_pseudoprime" .Vb 2 \& my $maybe_prime = is_strong_pseudoprime($n, 2); \& my $probably_prime = is_strong_pseudoprime($n, 2, 3, 5, 7, 11, 13, 17); .Ve .PP Takes a positive number \f(CW\*(C`n\*(C'\fR and one or more non-zero positive bases as input. Returns \f(CW1\fR if the input is a strong probable prime to each base, \f(CW0\fR if not. .PP 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. .PP This is usually used in combination with other tests to make either stronger tests (e.g. the strong \s-1BPSW\s0 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 \s-1MR\s0 tests for primality testing. .PP 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 \f(CW\*(C`is_strong_pseudoprime\*(C'\fR function. .SS "is_lucas_pseudoprime" .IX Subsection "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 \&\s-1OEIS A217120\s0 . .SS "is_strong_lucas_pseudoprime" .IX Subsection "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 \s-1BPSW\s0 primality test (the Miller-Rabin strong pseudoprime test with base 2 being the other half). Removing primes, this produces the sequence \&\s-1OEIS A217255\s0 . .SS "is_extra_strong_lucas_pseudoprime" .IX Subsection "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 ). This test has more stringent conditions than the strong Lucas test, and produces about 60% fewer pseudoprimes. Performance is typically 20\-30% \fIfaster\fR than the strong Lucas test. .PP The parameters are selected using the Baillie-OEIS method method: increment \f(CW\*(C`P\*(C'\fR from \f(CW3\fR until \f(CW\*(C`jacobi(D,n) = \-1\*(C'\fR. Removing primes, this produces the sequence \&\s-1OEIS A217719\s0 . .SS "is_almost_extra_strong_lucas_pseudoprime" .IX Subsection "is_almost_extra_strong_lucas_pseudoprime" This is similar to the \*(L"is_extra_strong_lucas_pseudoprime\*(R" function, but does not calculate \f(CW\*(C`U\*(C'\fR, 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. .PP An optional second argument (an integer between 1 and 256) indicates the increment amount for \f(CW\*(C`P\*(C'\fR parameter selection. The default value of 1 yields the parameter selection described in \*(L"is_extra_strong_lucas_pseudoprime\*(R", creating a pseudoprime sequence which is a superset of the latter's pseudoprime sequence \s-1OEIS A217719\s0 . A value of 2 yields the method used by Pari . .PP Because the \f(CW\*(C`U = 0\*(C'\fR 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 \s-1BPSW\s0 counterexamples have been found with any of the Lucas tests described. .SS "is_perrin_pseudoprime" .IX Subsection "is_perrin_pseudoprime" Takes a positive number \f(CW\*(C`n\*(C'\fR as input and returns 1 if \f(CW\*(C`n\*(C'\fR divides \f(CWP(n)\fR where \f(CWP(n)\fR is the Perrin number of \f(CW\*(C`n\*(C'\fR. The Perrin sequence is defined by .PP .Vb 1 \& C .Ve .PP While pseudoprimes are relatively rare (the first two are 271441 and 904631), infinitely many exist. The pseudoprime sequence is \s-1OEIS A013998\s0 . .PP The implementation uses modular 3x3 matrix exponentiation, which is efficient but still quite slow compared to the other probable prime tests. .SS "is_frobenius_pseudoprime" .IX Subsection "is_frobenius_pseudoprime" Takes a positive number \f(CW\*(C`n\*(C'\fR as input, and two optional parameters \f(CW\*(C`a\*(C'\fR and \&\f(CW\*(C`b\*(C'\fR, and returns 1 if the \f(CW\*(C`n\*(C'\fR is a Frobenius probable prime with respect to the polynomial \f(CW\*(C`x^2 \- ax + b\*(C'\fR. Without the parameters, \f(CW\*(C`b = 2\*(C'\fR and \&\f(CW\*(C`a\*(C'\fR is the least positive odd number such that \f(CW\*(C`(a^2\-4b|n) = \-1\*(C'\fR. This selection has no pseudoprimes below \f(CW\*(C`2^64\*(C'\fR and none known. In any case, the discriminant \f(CW\*(C`a^2\-4b\*(C'\fR must not be a perfect square. .PP Some authors use the Fibonacci polynomial \f(CW\*(C`x^2\-x\-1\*(C'\fR corresponding to \&\f(CW\*(C`(1,\-1)\*(C'\fR 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 \f(CW\*(C`(3,\-5)\*(C'\fR), so is not used as the default here. With the \f(CW\*(C`(1,\-1)\*(C'\fR parameters the pseudoprime sequence is \s-1OEIS A212424\s0 . .PP The Frobenius test is a stronger test than the Lucas test. Any Frobenius \&\f(CW\*(C`(a,b)\*(C'\fR pseudoprime is also a Lucas \f(CW\*(C`(a,b)\*(C'\fR pseudoprime but the converse is not true, as any Frobenius \f(CW\*(C`(a,b)\*(C'\fR pseudoprime is also a Fermat pseudoprime to the base \f(CW\*(C`|b|\*(C'\fR. We can see that with the default parameters this is similar to, but somewhat weaker than, the \s-1BPSW\s0 test used by this module (which uses the strong and extra-strong versions of the probable prime and Lucas tests respectively). .PP The performance cost is slightly more than 3 strong pseudoprime tests. Also see \*(L"is_frobenius_underwood_pseudoprime\*(R" which is an extremely efficient construction of a Frobenius test using good parameter selection, allowing it to run 1.5 to 2 times faster than the general Frobenius test. .SS "is_frobenius_underwood_pseudoprime" .IX Subsection "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 \f(CW\*(C`a\*(C'\fR as the least positive integer such that \f(CW\*(C`(a^2\-4|n)=\-1\*(C'\fR, then verifies that \&\f(CW\*(C`(2+2)^(n+1) = 2a + 5 mod (x^2\-ax+1,n)\*(C'\fR. 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 Frobenius test. .PP There are no known pseudoprimes to this test. This test also has no overlap with the \s-1BPSW\s0 test, making it a very effective method for adding additional certainty. .SS "miller_rabin_random" .IX Subsection "miller_rabin_random" Takes a positive number (\f(CW\*(C`n\*(C'\fR) as input and a positive number (\f(CW\*(C`k\*(C'\fR) of bases to use. Performs \f(CW\*(C`k\*(C'\fR Miller-Rabin tests using uniform random bases between 2 and \f(CW\*(C`n\-2\*(C'\fR. .PP This should not be used in place of \*(L"is_prob_prime\*(R", \*(L"is_prime\*(R", or \*(L"is_provable_prime\*(R". Those functions will be faster and provide better results than running \f(CW\*(C`k\*(C'\fR 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. .SS "is_prob_prime" .IX Subsection "is_prob_prime" .Vb 2 \& my $prob_prime = is_prob_prime($n); \& # Returns 0 (composite), 2 (prime), or 1 (probably prime) .Ve .PP Takes a positive number as input and returns back either 0 (composite), 2 (definitely prime), or 1 (probably prime). .PP 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 \s-1BPSW\s0 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). .PP For inputs larger than \f(CW\*(C`2^64\*(C'\fR, an extra-strong Baillie-PSW primality test is performed (also called \s-1BPSW\s0 or \s-1BSW\s0). 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 \s-1BPSW\s0 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. .SS "is_bpsw_prime" .IX Subsection "is_bpsw_prime" Given a positive number input, returns 0 (composite), 2 (definitely prime), or 1 (probably prime), using the \s-1BPSW\s0 primality test (extra-strong variant). Normally one of the \*(L"is_prime\*(R" in Math::Prime::Util or \&\*(L"is_prob_prime\*(R" 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 \*(L"is_bpsw_prime\*(R" may save a small amount of time. .SS "is_provable_prime" .IX Subsection "is_provable_prime" .Vb 1 \& say "$n is definitely prime" if is_provable_prime($n) == 2; .Ve .PP 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 \*(L"is_prime\*(R" and \*(L"is_prob_prime\*(R". Note that numbers below 2^64 are considered proven by the deterministic set of Miller-Rabin bases or the \&\s-1BPSW\s0 test. Both of these have been tested for all small (64\-bit) composites and do not return false positives. .PP Using the Math::Prime::Util::GMP module is \fBhighly recommended\fR 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 \&\s-1MPU::GMP\s0 works for primes of hundreds of digits (thousands with an optional larger polynomial set). .PP The pure Perl implementation uses theorem 5 of \s-1BLS75 \s0(Brillhart, Lehmer, and Selfridge's 1975 paper), an improvement on the Pocklington-Lehmer test. This requires \f(CW\*(C`n\-1\*(C'\fR to be factored to \f(CW\*(C`(n/2)^(1/3))\*(C'\fR. This is often fast, but as \f(CW\*(C`n\*(C'\fR gets larger, it takes exponentially longer to find factors. .PP Math::Prime::Util::GMP implements both the \s-1BLS75\s0 theorem 5 test as well as \s-1ECPP \s0(elliptic curve primality proving). It will typically try a quick \&\f(CW\*(C`n\-1\*(C'\fR proof before using \s-1ECPP. \s0 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 \f(CW\*(C`is_prime(n,1)\*(C'\fR which is the default for Math::Pari. .SS "prime_certificate" .IX Subsection "prime_certificate" .Vb 2 \& my $cert = prime_certificate($n); \& say verify_prime($cert) ? "proven prime" : "not prime"; .Ve .PP Given a positive integer \f(CW\*(C`n\*(C'\fR as input, returns a primality certificate as a multi-line string. If we could not prove \f(CW\*(C`n\*(C'\fR prime, an empty string is returned (\f(CW\*(C`n\*(C'\fR may or may not be composite). This may be examined or given to \*(L"verify_prime\*(R" for verification. The latter function contains the description of the format. .SS "is_provable_prime_with_cert" .IX Subsection "is_provable_prime_with_cert" Given a positive integer as input, returns a two element array containing the result of \*(L"is_provable_prime\*(R": 0 definitely composite 1 probably prime 2 definitely prime and a primality certificate like \*(L"prime_certificate\*(R". The certificate will be an empty string if the first element is not 2. .SS "verify_prime" .IX Subsection "verify_prime" .Vb 2 \& my $cert = prime_certificate($n); \& say verify_prime($cert) ? "proven prime" : "not prime"; .Ve .PP 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 \s-1ECPP\s0 certificates will be faster with Math::Prime::Util::GMP for its elliptic curve computations. .PP If the certificate is malformed, the routine will carp a warning in addition to returning 0. If the \f(CW\*(C`verbose\*(C'\fR option is set (see \*(L"prime_set_config\*(R") then if the validation fails, the reason for the failure is printed in addition to returning 0. If the \f(CW\*(C`verbose\*(C'\fR option is set to 2 or higher, then a message indicating success and the certificate type is also printed. .PP 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: .PP .Vb 1 \& [MPU \- Primality Certificate] .Ve .PP All lines in the certificate beginning with \f(CW\*(C`#\*(C'\fR are treated as comments and ignored, as are blank lines. A version number may follow, such as: .PP .Vb 1 \& Version 1.0 .Ve .PP For all inputs, base 10 is the default, but at any point this may be changed with a line like: .PP .Vb 1 \& Base 16 .Ve .PP where allowed bases are 10, 16, and 62. This module will only use base 10, so its routines will not output Base commands. .PP Next, we look for (using \*(L"100003\*(R" as an example): .PP .Vb 2 \& Proof for: \& N 100003 .Ve .PP where the text \f(CW\*(C`Proof for:\*(C'\fR indicates we will read an \f(CW\*(C`N\*(C'\fR value. Skipping comments and blank lines, the next line should be \*(L"N \*(R" followed by the number. .PP After this, we read one or more blocks. Each block is a proof of the form: .PP .Vb 1 \& If Q is prime, then N is prime. .Ve .PP 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. .PP 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 \f(CW\*(C`N\*(C'\fR, and for each node (block), all \f(CW\*(C`Q\*(C'\fR values must either have a block using that value as its \f(CW\*(C`N\*(C'\fR or \f(CW\*(C`Q\*(C'\fR must be less than \f(CW\*(C`2^64\*(C'\fR and pass \s-1BPSW.\s0 .PP Some other certificate formats (e.g. Primo) use an ordered chain, where the first block must be for the initial \f(CW\*(C`N\*(C'\fR, a single \f(CW\*(C`Q\*(C'\fR is given which is the implied \f(CW\*(C`N\*(C'\fR 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 \s-1BLS75\s0 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. .PP The blocks begin with the text \*(L"Type ...\*(R" where ... is the type. One or more values follow. The defined types are: .ie n .IP """Small""" 4 .el .IP "\f(CWSmall\fR" 4 .IX Item "Small" .Vb 2 \& Type Small \& N 5791 .Ve .Sp N must be less than 2^64 and be prime (use \s-1BPSW\s0 or deterministic M\-R). .ie n .IP """BLS3""" 4 .el .IP "\f(CWBLS3\fR" 4 .IX Item "BLS3" .Vb 4 \& Type BLS3 \& N 2297612322987260054928384863 \& Q 16501461106821092981 \& A 5 .Ve .Sp A simple n\-1 style proof using \s-1BLS75\s0 theorem 3. This block verifies if: a Q is odd b Q > 2 c Q divides N\-1 . Let M = (N\-1)/Q d \s-1MQ+1\s0 = 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 .ie n .IP """Pocklington""" 4 .el .IP "\f(CWPocklington\fR" 4 .IX Item "Pocklington" .Vb 4 \& Type Pocklington \& N 2297612322987260054928384863 \& Q 16501461106821092981 \& A 5 .Ve .Sp A simple n\-1 style proof using generalized Pocklington. This is more restrictive than \s-1BLS3\s0 and much more than \s-1BLS5. \s0 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 \s-1MQ+1\s0 = N e A > 1 f A^(N\-1) mod N = 1 g gcd(A^M \- 1, N) = 1 .ie n .IP """BLS15""" 4 .el .IP "\f(CWBLS15\fR" 4 .IX Item "BLS15" .Vb 5 \& Type BLS15 \& N 8087094497428743437627091507362881 \& Q 175806402118016161687545467551367 \& LP 1 \& LQ 22 .Ve .Sp A simple n+1 style proof using \s-1BLS75\s0 theorem 15. This block verifies if: a Q is odd b Q > 2 c Q divides N+1 . Let M = (N+1)/Q d \s-1MQ\-1\s0 = 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 \s-1LP,LQ \s0 i V_{m/2} mod N != 0 j V_{(N+1)/2} mod N == 0 .ie n .IP """BLS5""" 4 .el .IP "\f(CWBLS5\fR" 4 .IX Item "BLS5" .Vb 6 \& Type BLS5 \& N 8087094497428743437627091507362881 \& Q[1] 98277749 \& Q[2] 3631 \& A[0] 11 \& \-\-\-\- .Ve .Sp A more sophisticated n\-1 proof using \s-1BLS\s0 theorem 5. This requires N\-1 to be factored only to \f(CW\*(C`(N/2)^(1/3)\*(C'\fR. 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: .Sp .Vb 10 \& 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 .Ve .ie n .IP """ECPP""" 4 .el .IP "\f(CWECPP\fR" 4 .IX Item "ECPP" .Vb 8 \& Type ECPP \& N 175806402118016161687545467551367 \& A 96642115784172626892568853507766 \& B 111378324928567743759166231879523 \& M 175806402118016177622955224562171 \& Q 2297612322987260054928384863 \& X 3273750212 \& Y 82061726986387565872737368000504 .Ve .Sp An elliptic curve primality block, typically generated with an Atkin/Morain \&\s-1ECPP\s0 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: .Sp .Vb 10 \& . 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 .Ve .SS "is_aks_prime" .IX Subsection "is_aks_prime" .Vb 1 \& say "$n is definitely prime" if is_aks_prime($n); .Ve .PP Takes a positive number as input, and returns 1 if the input passes the Agrawal-Kayal-Saxena (\s-1AKS\s0) primality test. This is a deterministic unconditional primality test which runs in polynomial time for general input. .PP 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 \fBmillions\fR of times slower than other algorithms. From R.P. Brent, 2010: \*(L"\s-1AKS\s0 is not a practical algorithm. \s-1ECPP\s0 is much faster.\*(R" We have \s-1ECPP,\s0 and indeed it is much faster. .PP This implementation includes the v6 improvements from Lenstra as well as further improvements from Bernstein and Voloch. It runs substantially faster than the original or v6 versions. The \s-1GMP\s0 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 \s-1GMP\s0 excels. Because of this, the \s-1GMP\s0 implementation is likely to be faster once the input is larger than \f(CW\*(C`2^32\*(C'\fR. .SS "is_power" .IX Subsection "is_power" .Vb 3 \& 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"; .Ve .PP Given a single positive integer input \f(CW\*(C`n\*(C'\fR, returns k if \f(CW\*(C`n = p^k\*(C'\fR for some integer \f(CW\*(C`p > 1, k > 1\*(C'\fR, and 0 otherwise. The k returned is the largest possible. This can be used in a boolean statement to determine if \f(CW\*(C`n\*(C'\fR is a perfect power. .PP If given two arguments \f(CW\*(C`n\*(C'\fR and \f(CW\*(C`k\*(C'\fR, returns 1 if \f(CW\*(C`n\*(C'\fR is a \f(CW\*(C`k\-th\*(C'\fR power, and 0 otherwise. For example, if \f(CW\*(C`k=2\*(C'\fR then this detects perfect squares. .PP This corresponds to Pari/GP's \f(CW\*(C`ispower\*(C'\fR function, with the limitations of only integer arguments and no third argument may be given to return the root. .SS "lucas_sequence" .IX Subsection "lucas_sequence" .Vb 1 \& my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k) .Ve .PP Computes \f(CW\*(C`U_k\*(C'\fR, \f(CW\*(C`V_k\*(C'\fR, and \f(CW\*(C`Q_k\*(C'\fR for the Lucas sequence defined by \&\f(CW\*(C`P\*(C'\fR,\f(CW\*(C`Q\*(C'\fR, modulo \f(CW\*(C`n\*(C'\fR. The modular Lucas sequence is used in a number of primality tests and proofs. The following conditions must hold: \&\f(CW\*(C` |P| < n\*(C'\fR ; \&\f(CW\*(C` |Q| < n\*(C'\fR ; \&\f(CW\*(C` k >= 0\*(C'\fR ; \&\f(CW\*(C` n >= 2\*(C'\fR. .SS "gcd" .IX Subsection "gcd" Given a list of integers, returns the greatest common divisor. This is often used to test for coprimality . .SS "lcm" .IX Subsection "lcm" Given a list of integers, returns the least common multiple. Note that we follow the semantics of Mathematica, Pari, and Perl 6, re: .PP .Vb 2 \& lcm(0, n) = 0 Any zero in list results in zero return \& lcm(n,\-m) = lcm(n, m) We use the absolute values .Ve .SS "gcdext" .IX Subsection "gcdext" Given two integers \f(CW\*(C`x\*(C'\fR and \f(CW\*(C`y\*(C'\fR, returns \f(CW\*(C`u,v,d\*(C'\fR such that \f(CW\*(C`d = gcd(x,y)\*(C'\fR and \f(CW\*(C`u*x + v*y = d\*(C'\fR. This uses the extended Euclidian algorithm to compute the values satisfying Be\*'zout's Identity. .PP This corresponds to Pari's \f(CW\*(C`gcdext\*(C'\fR function, which was renamed from \&\f(CW\*(C`bezout\*(C'\fR out Pari 2.6. The results will hence match \*(L"bezout\*(R" in Math::Pari. .SS "chinese" .IX Subsection "chinese" .Vb 1 \& say chinese( [14,643], [254,419], [87,733] ); # 87041638 .Ve .PP Solves a system of simultaneous congruences using the Chinese Remainder Theorem (with extension to non-coprime moduli). A list of \f(CW\*(C`[a,n]\*(C'\fR pairs are taken as input, each representing an equation \f(CW\*(C`x X a mod n\*(C'\fR. If no solution exists, \f(CW\*(C`undef\*(C'\fR is returned. If a solution is returned, the modulus is equal to the lcm of all the given moduli (see \*(L"lcm\*(R". In the standard case where all values of \f(CW\*(C`n\*(C'\fR are coprime, this is just the product. The \f(CW\*(C`n\*(C'\fR values must be positive integers, while the \f(CW\*(C`a\*(C'\fR values are integers. .PP Comparison to similar functions in other software: .PP .Vb 2 \& 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, ...}] .Ve .SS "vecsum" .IX Subsection "vecsum" .Vb 1 \& say "Totient sum 500,000: ", vecsum(euler_phi(0,500_000)); .Ve .PP Returns the sum of all arguments, each of which must be an integer. This is similar to List::Util's \*(L"sum0\*(R" 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. \f(CW\*(C`vecsum\*(C'\fR sums (signed) integers and returns the untruncated result. Processing is done on native integers while possible. .SS "vecprod" .IX Subsection "vecprod" .Vb 1 \& say "Totient product 5,000: ", vecprod(euler_phi(1,5_000)); .Ve .PP Returns the product of all arguments, each of which must be an integer. This is similar to List::Util's \*(L"product\*(R" in List::Util function, but keeps all results as integers and automatically switches to bigints if needed. .SS "vecmin" .IX Subsection "vecmin" .Vb 1 \& say "Smallest Totient 100k\-200k: ", vecmin(euler_phi(100_000,200_000)); .Ve .PP Returns the minimum of all arguments, each of which must be an integer. This is similar to List::Util's \*(L"min\*(R" 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. \f(CW\*(C`vecmin\*(C'\fR validates and compares all results as integers. The validation step will make it a little slower than \*(L"min\*(R" in List::Util but this prevents accidental and unintentional use of floats. .SS "vecmax" .IX Subsection "vecmax" .Vb 1 \& say "Largest Totient 100k\-200k: ", vecmax(euler_phi(100_000,200_000)); .Ve .PP Returns the maximum of all arguments, each of which must be an integer. This is similar to List::Util's \*(L"max\*(R" 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. \f(CW\*(C`vecmax\*(C'\fR validates and compares all results as integers. The validation step will make it a little slower than \*(L"max\*(R" in List::Util but this prevents accidental and unintentional use of floats. .SS "vecreduce" .IX Subsection "vecreduce" .Vb 2 \& say "Count of non\-zero elements: ", vecreduce { $a + !!$b } (0,@v); \& my $checksum = vecreduce { $a ^ $b } @{twin_primes(1000000)}; .Ve .PP Does a reduce operation via left fold. Takes a block and a list as arguments. The block uses the special local variables \f(CW\*(C`a\*(C'\fR and \f(CW\*(C`b\*(C'\fR 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 \f(CW\*(C`undef\*(C'\fR will be returned with an empty list. .PP The interface is exactly the same as \*(L"reduce\*(R" 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. .PP 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. .SS "invmod" .IX Subsection "invmod" .Vb 1 \& say "The inverse of 42 mod 2017 = ", invmod(42,2017); .Ve .PP Given two integers \f(CW\*(C`a\*(C'\fR and \f(CW\*(C`n\*(C'\fR, return the inverse of \f(CW\*(C`a\*(C'\fR modulo \f(CW\*(C`n\*(C'\fR. If not defined, undef is returned. If defined, then the return value multiplied by \f(CW\*(C`a\*(C'\fR equals \f(CW1\fR modulo \f(CW\*(C`n\*(C'\fR. .PP This results correspond to the Pari result of \f(CW\*(C`lift(Mod(1/a,n))\*(C'\fR. The semantics with respect to negative arguments match Pari. Notably, a negative \f(CW\*(C`n\*(C'\fR is negated, which is different from Math::BigInt, but in both cases the return value is still congruent to \f(CW1\fR modulo \f(CW\*(C`n\*(C'\fR as expected. .SS "valuation" .IX Subsection "valuation" .Vb 1 \& say "$n is divisible by 2 ", valuation($n,2), " times."; .Ve .PP Given integers \f(CW\*(C`n\*(C'\fR and \f(CW\*(C`k\*(C'\fR, returns the numbers of times \f(CW\*(C`n\*(C'\fR is divisible by \f(CW\*(C`k\*(C'\fR. This is a very limited version of the algebraic valuation meaning, just applied to integers. This corresponds to Pari's \f(CW\*(C`valuation\*(C'\fR function. \&\f(CW0\fR is returned if \f(CW\*(C`n\*(C'\fR or \f(CW\*(C`k\*(C'\fR is one of the values \f(CW\*(C`\-1\*(C'\fR, \f(CW0\fR, or \f(CW1\fR. .SS "hammingweight" .IX Subsection "hammingweight" Given an integer \f(CW\*(C`n\*(C'\fR, returns the binary Hamming weight of \f(CWabs(n)\fR. This is also called the population count, and is the number of 1s in the binary representation. This corresponds to Pari's \f(CW\*(C`hammingweight\*(C'\fR function for \&\f(CW\*(C`t_INT\*(C'\fR arguments. .SS "moebius" .IX Subsection "moebius" .Vb 3 \& 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)); .Ve .PP Returns X(n), the Mo\*:bius function (also known as the Moebius, Mobius, or MoebiusMu function) for an integer input. This function is 1 if \&\f(CW\*(C`n = 1\*(C'\fR, 0 if \f(CW\*(C`n\*(C'\fR is not square free (i.e. \f(CW\*(C`n\*(C'\fR has a repeated factor), and \f(CW\*(C`\-1^t\*(C'\fR if \f(CW\*(C`n\*(C'\fR is a product of \f(CW\*(C`t\*(C'\fR distinct primes. This is an important function in prime number theory. Like \s-1SAGE,\s0 we define \&\f(CW\*(C`moebius(0) = 0\*(C'\fR for convenience. .PP If called with two arguments, they define a range \f(CW\*(C`low\*(C'\fR to \f(CW\*(C`high\*(C'\fR, and the function returns an array with the value of the Mo\*: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 Dele\*'glise and Rivat (1996) algorithm 4.1, which is a segmented version of Lioen and van de Lune (1994) algorithm 3.2. .PP 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). .SS "mertens" .IX Subsection "mertens" .Vb 1 \& say "Mertens(10M) = ", mertens(10_000_000); # = 1037 .Ve .PP Returns M(n), the Mertens function for a non-negative integer input. This function is defined as \f(CW\*(C`sum(moebius(1..n))\*(C'\fR, but calculated more efficiently for large inputs. For example, computing Mertens(100M) takes: .PP .Vb 4 \& time approx mem \& 0.4s 0.1MB mertens(100_000_000) \& 5.6s 880MB vecsum(moebius(1,100_000_000)) \& 98s 0MB $sum += moebius($_) for 1..100_000_000 .Ve .PP 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 \f(CW\*(C`n^1/2\*(C'\fR version of Dele\*'glise and Rivat (1996), which involves calculating all moebius values to \f(CW\*(C`n^1/2\*(C'\fR, which in turn will require prime sieving to \f(CW\*(C`n^1/4\*(C'\fR. .PP Various algorithms exist for this, using differing quantities of X(n). The simplest way is to efficiently sum all \f(CW\*(C`n\*(C'\fR values. Benito and Varona (2008) show a clever and simple method that only requires \f(CW\*(C`n/3\*(C'\fR values. Dele\*'glise and Rivat (1996) describe a segmented method using only \f(CW\*(C`n^1/3\*(C'\fR values. The current implementation does a simple non-segmented \f(CW\*(C`n^1/2\*(C'\fR 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. .SS "euler_phi" .IX Subsection "euler_phi" .Vb 1 \& say "The Euler totient of $n is ", euler_phi($n); .Ve .PP Returns X(n), the Euler totient function (also called Euler's phi or phi function) for an integer value. This is an arithmetic function which counts the number of positive integers less than or equal to \f(CW\*(C`n\*(C'\fR that are relatively prime to \f(CW\*(C`n\*(C'\fR. Given the definition used, \f(CW\*(C`euler_phi\*(C'\fR will return 0 for all \&\f(CW\*(C`n < 1\*(C'\fR. This follows the logic used by \s-1SAGE. \s0 Mathematica and Pari return \f(CW\*(C`euler_phi(\-n)\*(C'\fR for \f(CW\*(C`n < 0\*(C'\fR. Mathematica returns 0 for \f(CW\*(C`n = 0\*(C'\fR, Pari pre\-2.6.2 raises and exception, and Pari 2.6.2 and newer returns 2. .PP If called with two arguments, they define a range \f(CW\*(C`low\*(C'\fR to \f(CW\*(C`high\*(C'\fR, and the function returns an array with the totient of every n from low to high inclusive. .SS "jordan_totient" .IX Subsection "jordan_totient" .Vb 1 \& say "Jordan\*(Aqs totient J_$k($n) is ", jordan_totient($k, $n); .Ve .PP Returns Jordan's totient function for a given integer value. Jordan's totient is a generalization of Euler's totient, where \f(CW\*(C`jordan_totient(1,$n) == euler_totient($n)\*(C'\fR This counts the number of k\-tuples less than or equal to n that form a coprime tuple with n. As with \f(CW\*(C`euler_phi\*(C'\fR, 0 is returned for all \f(CW\*(C`n < 1\*(C'\fR. This function can be used to generate some other useful functions, such as the Dedekind psi function, where \f(CW\*(C`psi(n) = J(2,n) / J(1,n)\*(C'\fR. .SS "exp_mangoldt" .IX Subsection "exp_mangoldt" .Vb 1 \& say "exp(lambda($_)) = ", exp_mangoldt($_) for 1 .. 100; .Ve .PP Returns \s-1EXP\s0(X(n)), the exponential of the Mangoldt function (also known as von Mangoldt's function) for an integer value. The Mangoldt function is equal to log p if n is prime or a power of a prime, and 0 otherwise. We return the exponential so all results are integers. Hence the return value for \f(CW\*(C`exp_mangoldt\*(C'\fR is: .PP .Vb 2 \& p if n = p^m for some prime p and integer m >= 1 \& 1 otherwise. .Ve .SS "liouville" .IX Subsection "liouville" Returns X(n), the Liouville function for a non-negative integer input. This is \-1 raised to X(n) (the total number of prime factors). .SS "chebyshev_theta" .IX Subsection "chebyshev_theta" .Vb 1 \& say chebyshev_theta(10000); .Ve .PP Returns X(n), the first Chebyshev function for a non-negative integer input. This is the sum of the logarithm of each prime where \f(CW\*(C`p <= n\*(C'\fR. This is effectively: .PP .Vb 1 \& my $s = 0; forprimes { $s += log($_) } $n; return $s; .Ve .SS "chebyshev_psi" .IX Subsection "chebyshev_psi" .Vb 1 \& say chebyshev_psi(10000); .Ve .PP Returns X(n), the second Chebyshev function for a non-negative integer input. This is the sum of the logarithm of each prime power where \f(CW\*(C`p^k <= n\*(C'\fR for an integer k. An alternate but slower computation is as the summatory Mangoldt function, such as: .PP .Vb 1 \& my $s = 0; for (1..$n) { $s += log(exp_mangoldt($_)) } return $s; .Ve .SS "divisor_sum" .IX Subsection "divisor_sum" .Vb 3 \& 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); .Ve .PP This function takes a positive integer as input and returns the sum of its divisors, including 1 and itself. An optional second argument \f(CW\*(C`k\*(C'\fR may be given, which will result in the sum of the \f(CW\*(C`k\-th\*(C'\fR powers of the divisors to be returned. .PP This is known as the sigma function (see Hardy and Wright section 16.7, or \s-1OEIS A000203\s0). The \s-1API\s0 is identical to Pari/GP's \f(CW\*(C`sigma\*(C'\fR function. This function is useful for calculating things like aliquot sums, abundant numbers, perfect numbers, etc. .PP 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 \f(CW\*(C`sumdiv\*(C'\fR function. .PP An example of the 5th Jordan totient (\s-1OEIS A059378\s0): .PP .Vb 1 \& divisor_sum( $n, sub { my $d=shift; $d**5 * moebius($n/$d); } ); .Ve .PP though we have a function \*(L"jordan_totient\*(R" which is more efficient. .PP 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. .SS "primorial" .IX Subsection "primorial" .Vb 1 \& $prim = primorial(11); # 11# = 2*3*5*7*11 = 2310 .Ve .PP Returns the primorial \f(CW\*(C`n#\*(C'\fR of the positive integer input, defined as the product of the prime numbers less than or equal to \f(CW\*(C`n\*(C'\fR. This is the \&\s-1OEIS\s0 series A034386 : primorial numbers second definition. .PP .Vb 2 \& primorial(0) == 1 \& primorial($n) == pn_primorial( prime_count($n) ) .Ve .PP The result will be a Math::BigInt object if it is larger than the native bit size. .PP Be careful about which version (\f(CW\*(C`primorial\*(C'\fR or \f(CW\*(C`pn_primorial\*(C'\fR) 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. \&\s-1OEIS,\s0 Wikipedia, and Mathworld are all consistent, and these functions should match that terminology. This function should return the same result as the \&\f(CW\*(C`mpz_primorial_ui\*(C'\fR function added in \s-1GMP 5.1.\s0 .SS "pn_primorial" .IX Subsection "pn_primorial" .Vb 1 \& $prim = pn_primorial(5); # p_5# = 2*3*5*7*11 = 2310 .Ve .PP Returns the primorial number \f(CW\*(C`p_n#\*(C'\fR of the positive integer input, defined as the product of the first \f(CW\*(C`n\*(C'\fR prime numbers (compare to the factorial, which is the product of the first \f(CW\*(C`n\*(C'\fR natural numbers). This is the \&\s-1OEIS\s0 series A002110 : primorial numbers first definition. .PP .Vb 2 \& pn_primorial(0) == 1 \& pn_primorial($n) == primorial( nth_prime($n) ) .Ve .PP The result will be a Math::BigInt object if it is larger than the native bit size. .SS "consecutive_integer_lcm" .IX Subsection "consecutive_integer_lcm" .Vb 1 \& $lcm = consecutive_integer_lcm($n); .Ve .PP Given an unsigned integer argument, returns the least common multiple of all integers from 1 to \f(CW\*(C`n\*(C'\fR. This can be done by manipulation of the primes up to \f(CW\*(C`n\*(C'\fR, resulting in much faster and memory-friendly results than using a factorial. .SS "partitions" .IX Subsection "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 \f(CW\*(C`numbpart\*(C'\fR function and Mathematica's \f(CW\*(C`PartitionsP\*(C'\fR function. The values produced in order are \s-1OEIS\s0 series A000041 . .PP This uses a combinatorial calculation, which means it will not be very fast compared to Pari, Mathematica, or \s-1FLINT\s0 which use the Rademacher formula using multi-precision floating point. In 10 seconds: .PP .Vb 6 \& 70 Integer::Partition \& 90 MPU forpart { $n++ } \& 10_000 MPU pure Perl partitions \& 250_000 MPU GMP partitions \& 35_000_000 Pari\*(Aqs numbpart \& 500_000_000 Jonathan Bober\*(Aqs partitions_c.cc v0.6 .Ve .PP If you want the enumerated partitions, see \*(L"forpart\*(R". .SS "carmichael_lambda" .IX Subsection "carmichael_lambda" Returns the Carmichael function (also called the reduced totient function, or Carmichael X(n)) of a positive integer argument. It is the smallest positive integer \f(CW\*(C`m\*(C'\fR such that \f(CW\*(C`a^m = 1 mod n\*(C'\fR for every integer \f(CW\*(C`a\*(C'\fR coprime to \f(CW\*(C`n\*(C'\fR. This is \s-1OEIS\s0 series A002322 . .SS "kronecker" .IX Subsection "kronecker" Returns the Kronecker symbol \f(CW\*(C`(a|n)\*(C'\fR for two integers. The possible return values with their meanings for odd prime \f(CW\*(C`n\*(C'\fR are: .PP .Vb 3 \& 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) .Ve .PP The Kronecker symbol is an extension of the Jacobi symbol to all integer values of \f(CW\*(C`n\*(C'\fR from the latter's domain of positive odd values of \f(CW\*(C`n\*(C'\fR. The Jacobi symbol is itself an extension of the Legendre symbol, which is only defined for odd prime values of \f(CW\*(C`n\*(C'\fR. This corresponds to Pari's \&\f(CW\*(C`kronecker(a,n)\*(C'\fR function, Mathematica's \f(CW\*(C`KroneckerSymbol[n,m]\*(C'\fR function, and \s-1GMP\s0's \f(CW\*(C`mpz_kronecker(a,n)\*(C'\fR, \f(CW\*(C`mpz_jacobi(a,n)\*(C'\fR, and \&\f(CW\*(C`mpz_legendre(a,n)\*(C'\fR functions. .SS "factorial" .IX Subsection "factorial" Given positive integer argument \f(CW\*(C`n\*(C'\fR, returns the factorial of \f(CW\*(C`n\*(C'\fR, defined as the product of the integers 1 to \f(CW\*(C`n\*(C'\fR with the special case of \f(CW\*(C`factorial(0) = 1\*(C'\fR. This corresponds to Pari's \f(CWfactorial(n)\fR and Mathematica's \f(CW\*(C`Factorial[n]\*(C'\fR functions. .SS "binomial" .IX Subsection "binomial" Given integer arguments \f(CW\*(C`n\*(C'\fR and \f(CW\*(C`k\*(C'\fR, returns the binomial coefficient \&\f(CW\*(C`n*(n\-1)*...*(n\-k+1)/k!\*(C'\fR, also known as the choose function. Negative arguments use the Kronenburg extensions . This corresponds to Pari's \f(CW\*(C`binomial(n,k)\*(C'\fR function, Mathematica's \&\f(CW\*(C`Binomial[n,k]\*(C'\fR function, and \s-1GMP\s0's \f(CW\*(C`mpz_bin_ui\*(C'\fR function. .PP For negative arguments, this matches Mathematica. Pari does not implement the \f(CW\*(C`n < 0, k <= n\*(C'\fR extension and instead returns \f(CW0\fR for this case. \s-1GMP\s0's \s-1API\s0 does not allow negative \f(CW\*(C`k\*(C'\fR but otherwise matches. Math::BigInt does not implement any extensions and the results for \&\f(CW\*(C`n < 0, k \*(C'\fR 0> are undefined. .SS "bernfrac" .IX Subsection "bernfrac" Returns the Bernoulli number \f(CW\*(C`B_n\*(C'\fR for an integer argument \f(CW\*(C`n\*(C'\fR, as a rational number represented by two Math::BigInt objects. B_1 = 1/2. This corresponds to Pari's \f(CWbernfrac(n)\fR and Mathematica's \f(CW\*(C`BernoulliB\*(C'\fR functions. .PP This currently uses the simple Brent-Harvey recurrence, so will not be nearly as fast as Pari or Mathematica which use high-precision values of Pi and Zeta. With Math::Prime::Util::GMP installed it is, however, faster than Math::Pari which uses an older algorithm. .SS "bernreal" .IX Subsection "bernreal" Returns the Bernoulli number \f(CW\*(C`B_n\*(C'\fR for an integer argument \f(CW\*(C`n\*(C'\fR, as a Math::BigFloat object using the default precision. An optional second argument may be given specifying the precision to be used. .SS "stirling" .IX Subsection "stirling" .Vb 2 \& say "s(14,2) = ", stirling(14, 2); \& say "S(14,2) = ", stirling(14, 2, 2); .Ve .PP Returns the Stirling numbers of either the first kind (default) or second kind (with a third argument of 2). It takes two non-negative integer arguments \f(CW\*(C`n\*(C'\fR and \f(CW\*(C`k\*(C'\fR. This corresponds to Pari's \f(CW\*(C`stirling(n,k,{type})\*(C'\fR function and Mathematica's \f(CW\*(C`StirlingS1\*(C'\fR / \f(CW\*(C`StirlingS2\*(C'\fR functions. .PP Stirling numbers of the first kind are \f(CW\*(C`\-1^(n\-k)\*(C'\fR times the number of permutations of \f(CW\*(C`n\*(C'\fR symbols with exactly \f(CW\*(C`k\*(C'\fR cycles. Stirling numbers of the second kind are the number of ways to partition a set of \f(CW\*(C`n\*(C'\fR elements into \f(CW\*(C`k\*(C'\fR non-empty subsets. .SS "znorder" .IX Subsection "znorder" .Vb 1 \& $order = znorder(2, next_prime(10**16)\-6); .Ve .PP Given two positive integers \f(CW\*(C`a\*(C'\fR and \f(CW\*(C`n\*(C'\fR, returns the multiplicative order of \f(CW\*(C`a\*(C'\fR modulo \f(CW\*(C`n\*(C'\fR. This is the smallest positive integer \f(CW\*(C`k\*(C'\fR such that \&\f(CW\*(C`a^k X 1 mod n\*(C'\fR. Returns 1 if \f(CW\*(C`a = 1\*(C'\fR. Returns undef if \f(CW\*(C`a = 0\*(C'\fR or if \&\f(CW\*(C`a\*(C'\fR and \f(CW\*(C`n\*(C'\fR are not coprime, since no value will result in 1 mod n. This corresponds to Pari's \f(CW\*(C`znorder(Mod(a,n))\*(C'\fR function and Mathematica's \&\f(CW\*(C`MultiplicativeOrder[a,n]\*(C'\fR function. .SS "znprimroot" .IX Subsection "znprimroot" Given a positive integer \f(CW\*(C`n\*(C'\fR, returns the smallest primitive root of \f(CW\*(C`(Z/nZ)^*\*(C'\fR, or \f(CW\*(C`undef\*(C'\fR if no root exists. A root exists when \&\f(CW\*(C`euler_phi($n) == carmichael_lambda($n)\*(C'\fR, which will be true for all prime \f(CW\*(C`n\*(C'\fR and some composites. .PP \&\s-1OEIS A033948\s0 is a sequence of integers where the primitive root exists, while \s-1OEIS A046145\s0 is a list of the smallest primitive roots, which is what this function produces. .SS "znlog" .IX Subsection "znlog" .Vb 1 \& $k = znlog($a, $g, $p) .Ve .PP Returns the integer \f(CW\*(C`k\*(C'\fR that solves the equation \f(CW\*(C`a = g^k mod p\*(C'\fR, or undef if no solution is found. This is the discrete logarithm problem. .PP 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 relatively fast Shanks \&\s-1BSGS,\s0 as well as trial and Pollard's \s-1DLP\s0 Rho. .PP The \s-1PP\s0 implementation is less sophisticated, with only a memory-heavy \s-1BSGS\s0 being used. .SS "legendre_phi" .IX Subsection "legendre_phi" .Vb 1 \& $phi = legendre_phi(1000000000, 41); .Ve .PP Given a non-negative integer \f(CW\*(C`n\*(C'\fR and a non-negative prime number \f(CW\*(C`a\*(C'\fR, returns the Legendre phi function (also called Legendre's sum). This is the count of positive integers <= \f(CW\*(C`n\*(C'\fR which are not divisible by any of the first \f(CW\*(C`a\*(C'\fR primes. .SH "RANDOM PRIMES" .IX Header "RANDOM PRIMES" .SS "random_prime" .IX Subsection "random_prime" .Vb 2 \& my $small_prime = random_prime(1000); # random prime <= limit \& my $rand_prime = random_prime(100, 10000); # random prime within a range .Ve .PP 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. .PP 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 \&\f(CW\*(C`PRIMEINC\*(C'\fR, 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 \f(CW\*(C`random_prime(1000)\*(C'\fR, the numbers 2, 3, 487, 631, and 997 all have the same probability of being returned. .PP The configuration option \f(CW\*(C`use_primeinc\*(C'\fR can be set to override this and use the \s-1PRIMEINC\s0 algorithm for non-trivial sizes. This applies to all random prime functions. Never use this for crypto or if uniformly random primes are desired, but if you really don't care and just want any old prime in the range, setting this may make this run 2\-4x faster. .PP For small numbers, a random index selection is done, which gives ideal uniformity and is very efficient with small inputs. For ranges larger than this ~16\-bit threshold but within the native bit size, a Monte Carlo method is used (multiple calls to \f(CW\*(C`irand\*(C'\fR will be made if necessary). This also gives ideal uniformity and can be very fast for reasonably sized ranges. For even larger numbers, we partition the range, choose a random partition, then select a random prime from the partition. This gives some loss of uniformity but results in many fewer bits of randomness being consumed as well as being much faster. .PP If an \f(CW\*(C`irand\*(C'\fR function has been set via \*(L"prime_set_config\*(R", it will be used to construct any ranged random numbers needed. The function should return a uniformly random 32\-bit integer, which is how the irand functions exported by Math::Random::Secure, Math::Random::MT, Math::Random::ISAAC, and most other modules behave. .PP If no \f(CW\*(C`irand\*(C'\fR function was set, then Bytes::Random::Secure is used with a non-blocking seed. This will create good quality random numbers, so there should be little reason to change unless one is generating long-term keys, where using the blocking random source may be preferred. .PP Examples of various ways to set your own irand function: .PP .Vb 2 \& # System rand. You probably don\*(Aqt want to do this. \& prime_set_config(irand => sub { int(rand(4294967296)) }); \& \& # Math::Random::MTwist. Fastest RNG by quite a bit. \& use Math::Random::MTwist; \& prime_set_config(irand => \e&Math::Random::MTwist::_irand32); \& \& # Math::Random::Secure. Uses ISAAC and strong seed methods. \& use Math::Random::Secure; \& prime_set_config(irand => \e&Math::Random::Secure::irand); \& \& # Bytes::Random::Secure (OO interface with full control of options): \& use Bytes::Random::Secure (); \& BEGIN { \& my $rng = Bytes::Random::Secure\->new( Bits => 512 ); \& sub irand { return $rng\->irand; } \& } \& prime_set_config(irand => \e&irand); \& \& # Crypt::Random. Uses Pari and /dev/random. *VERY* slow. \& use Crypt::Random qw/makerandom/; \& prime_set_config(irand => sub { makerandom(Size=>32, Uniform=>1); }); \& \& # Net::Random. You probably don\*(Aqt want to use this, but if you do: \& use Net::Random; \& { my $rng = Net::Random\->new(src=>"fourmilab.ch",max=>0xFFFFFFFF); \& sub nr_irand { return $rng\->get(1); } } \& prime_set_config(irand => \e&nr_irand); \& \& # Go back to MPU\*(Aqs default configuration \& prime_set_config(irand => undef); .Ve .SS "random_ndigit_prime" .IX Subsection "random_ndigit_prime" .Vb 1 \& say "My 4\-digit prime number is: ", random_ndigit_prime(4); .Ve .PP Selects a random n\-digit prime, where the input is an integer number of digits. One of the primes within that range (e.g. 1000 \- 9999 for 4\-digits) will be uniformly selected using the \f(CW\*(C`irand\*(C'\fR function as described above. .PP 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 \f(CW\*(C`nobigint\*(C'\fR 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. .SS "random_nbit_prime" .IX Subsection "random_nbit_prime" .Vb 1 \& my $bigprime = random_nbit_prime(512); .Ve .PP Selects a random n\-bit prime, where the input is an integer number of bits. A prime with the nth bit set will be uniformly selected, with randomness supplied via calls to the \f(CW\*(C`irand\*(C'\fR function as described above. .PP For bit sizes of 64 and lower, \*(L"random_prime\*(R" 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 \&\*(L"random_prime\*(R" case while running slightly faster (in contrast, for large bit sizes \*(L"random_prime\*(R" selects a random upper partition then loops on the values within the partition, which very slightly skews the results towards smaller numbers). .PP The \f(CW\*(C`irand\*(C'\fR function is used for randomness, so all the discussion in \&\*(L"random_prime\*(R" about that applies here. The result will be a BigInt if the number of bits is greater than the native bit size. For better performance with large bit sizes, install Math::Prime::Util::GMP. .SS "random_strong_prime" .IX Subsection "random_strong_prime" .Vb 1 \& my $bigprime = random_strong_prime(512); .Ve .PP Constructs an n\-bit strong prime using Gordon's algorithm. We consider a strong prime \fIp\fR to be one where .IP "\(bu" 4 \&\fIp\fR is large. This function requires at least 128 bits. .IP "\(bu" 4 \&\fIp\-1\fR has a large prime factor \fIr\fR. .IP "\(bu" 4 \&\fIp+1\fR has a large prime factor \fIs\fR .IP "\(bu" 4 \&\fIr\-1\fR has a large prime factor \fIt\fR .PP 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 \s-1ECM. \s0 Second, modern factoring methods like \s-1GNFS\s0 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. .PP Similar to \*(L"random_nbit_prime\*(R", 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. .SS "random_proven_prime" .IX Subsection "random_proven_prime" .Vb 1 \& my $bigprime = random_proven_prime(512); .Ve .PP Constructs an n\-bit random proven prime. Internally this may use \&\*(L"is_provable_prime\*(R"(\*(L"random_nbit_prime\*(R") or \&\*(L"random_maurer_prime\*(R" depending on the platform and bit size. .SS "random_proven_prime_with_cert" .IX Subsection "random_proven_prime_with_cert" .Vb 1 \& my($n, $cert) = random_proven_prime_with_cert(512) .Ve .PP Similar to \*(L"random_proven_prime\*(R", 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 \*(L"prime_certificate\*(R" or \&\*(L"is_provable_prime_with_cert\*(R", and can be parsed by \*(L"verify_prime\*(R" or any other software that understands \s-1MPU\s0 primality certificates. .SS "random_maurer_prime" .IX Subsection "random_maurer_prime" .Vb 1 \& my $bigprime = random_maurer_prime(512); .Ve .PP 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 \*(L"random_nbit_prime\*(R", 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 \&\*(L"random_shawe_taylor_prime\*(R". .PP The differences between this function and that in Crypt::Primes are described in the \*(L"\s-1SEE ALSO\*(R"\s0 section. .PP Internally this additionally runs the \s-1BPSW\s0 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. .PP An alternative to this function is to run \*(L"is_provable_prime\*(R" on the result of \*(L"random_nbit_prime\*(R", which will provide more diversity and will be faster up to 512 or so bits. Maurer's method should be much faster for large bit sizes (larger than 2048). If you don't need absolutely proven results, then using \*(L"random_nbit_prime\*(R" followed by additional tests (\*(L"is_strong_pseudoprime\*(R" and/or \*(L"is_frobenius_underwood_pseudoprime\*(R") should be much faster. .SS "random_maurer_prime_with_cert" .IX Subsection "random_maurer_prime_with_cert" .Vb 1 \& my($n, $cert) = random_maurer_prime_with_cert(512) .Ve .PP As with \*(L"random_maurer_prime\*(R", 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 \*(L"prime_certificate\*(R" or \&\*(L"is_provable_prime_with_cert\*(R", and can be parsed by \*(L"verify_prime\*(R" or any other software that understands \s-1MPU\s0 primality certificates. The proof construction consists of a single chain of \f(CW\*(C`BLS3\*(C'\fR types. .SS "random_shawe_taylor_prime" .IX Subsection "random_shawe_taylor_prime" .Vb 1 \& my $bigprime = random_shawe_taylor_prime(8192); .Ve .PP Construct an n\-bit provable prime, using the Shawe-Taylor algorithm in section C.6 of \s-1FIPS 186\-4. \s0 This uses 512 bits of randomness and \s-1SHA\-256\s0 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 \s-1FIPS 186\-4\s0 algorithm. .PP Similar to \*(L"random_nbit_prime\*(R", 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 \&\*(L"random_maurer_prime\*(R" and \*(L"random_proven_prime\*(R". .PP Internally this additionally runs the \s-1BPSW\s0 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. .SS "random_shawe_taylor_prime_with_cert" .IX Subsection "random_shawe_taylor_prime_with_cert" .Vb 1 \& my($n, $cert) = random_shawe_taylor_prime_with_cert(4096) .Ve .PP As with \*(L"random_shawe_taylor_prime\*(R", 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 \*(L"prime_certificate\*(R" or \&\*(L"is_provable_prime_with_cert\*(R", and can be parsed by \*(L"verify_prime\*(R" or any other software that understands \s-1MPU\s0 primality certificates. The proof construction consists of a single chain of \f(CW\*(C`Pocklington\*(C'\fR types. .SH "UTILITY FUNCTIONS" .IX Header "UTILITY FUNCTIONS" .SS "prime_precalc" .IX Subsection "prime_precalc" .Vb 1 \& prime_precalc( 1_000_000_000 ); .Ve .PP 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. .SS "prime_memfree" .IX Subsection "prime_memfree" .Vb 1 \& prime_memfree; .Ve .PP Frees any extra memory the module may have allocated. Like with \&\f(CW\*(C`prime_precalc\*(C'\fR, 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. .SS "Math::Prime::Util::MemFree\->new" .IX Subsection "Math::Prime::Util::MemFree->new" .Vb 2 \& my $mf = Math::Prime::Util::MemFree\->new; \& # perform operations. When $mf goes out of scope, memory will be recovered. .Ve .PP This is a more robust way of making sure any cached memory is freed, as it will be handled by the last \f(CW\*(C`MemFree\*(C'\fR 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. .SS "prime_get_config" .IX Subsection "prime_get_config" .Vb 1 \& my $cached_up_to = prime_get_config\->{\*(Aqprecalc_to\*(Aq}; .Ve .PP 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: .PP .Vb 11 \& 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) \& use_primeinc allow the PRIMEINC random prime algorithm .Ve .SS "prime_set_config" .IX Subsection "prime_set_config" .Vb 1 \& prime_set_config( assume_rh => 1 ); .Ve .PP Allows setting of some parameters. Currently the only parameters are: .PP .Vb 8 \& 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, \& which means using Pure Perl code for big numbers. Set \& to 0 to disable GMP, set to 1 to re\-enable. \& You probably will never want to do this. \& \& assume_rh Allows functions to assume the Riemann hypothesis is \& true if set to 1. This defaults to 0. Currently this \& setting only impacts prime count lower and upper \& bounds, but could later be applied to other areas such \& as primality testing. A later version may also have a \& way to indicate whether no RH, RH, GRH, or ERH is to \& be assumed. \& \& irand Takes a code ref to an irand function returning a \& uniform number between 0 and 2**32\-1. This will be \& used for all random number generation in the module. \& \& use_primeinc When generating random primes, allow the PRIMEINC algorithm \& to be used. This can be 2\-4x faster than the default \& methods, but gives bad uniformity. .Ve .SH "FACTORING FUNCTIONS" .IX Header "FACTORING FUNCTIONS" .SS "factor" .IX Subsection "factor" .Vb 2 \& my @factors = factor(3_369_738_766_071_892_021); \& # returns (204518747,16476429743) .Ve .PP Produces the prime factors of a positive number input, in numerical order. The product of the returned factors will be equal to the input. \f(CW\*(C`n = 1\*(C'\fR will return an empty list, and \f(CW\*(C`n = 0\*(C'\fR will return 0. This matches Pari. .PP In scalar context, returns X(n), the total number of prime factors (\s-1OEIS A001222\s0 ). This corresponds to Pari's \f(CWbigomega(n)\fR function and Mathematica's \&\f(CW\*(C`PrimeOmega[n]\*(C'\fR function. This is same result that we would get if we evaluated the resulting array in scalar context. .PP The current algorithm does a little trial division, a check for perfect powers, followed by combinations of Pollard's Rho, \s-1SQUFOF,\s0 and Pollard's p\-1. The combination is applied to each non-prime factor found. .PP 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 \fBvery\fR slow for \*(L"hard\*(R" numbers. Installing the Math::Prime::Util::GMP module will speed up bigint factoring a \fBlot\fR, 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 \s-1GMP\s0 or Pari version of bigint if possible (e.g. \f(CW\*(C`use bigint try => \*(AqGMP,Pari\*(Aq\*(C'\fR), which will run 2\-3x faster (though still 100x slower than the real \s-1GMP\s0 code). .SS "factor_exp" .IX Subsection "factor_exp" .Vb 4 \& 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) .Ve .PP 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 \f(CW\*(C`FactorInteger[n]\*(C'\fR and Pari/GP's \f(CW\*(C`factorint\*(C'\fR functions return. Note that Math::Pari transposes the Pari result matrix. .PP In scalar context, returns X(n), the number of unique prime factors (\s-1OEIS A001221\s0 ). This corresponds to Pari's \f(CWomega(n)\fR function and Mathematica's \&\f(CW\*(C`PrimeNu[n]\*(C'\fR function. This is same result that we would get if we evaluated the resulting array in scalar context. .PP The internals are identical to \*(L"factor\*(R", so all comments there apply. Just the way the factors are arranged is different. .SS "divisors" .IX Subsection "divisors" .Vb 1 \& my @divisors = divisors(30); # returns (1, 2, 3, 5, 6, 10, 15, 30) .Ve .PP 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 \f(CW\*(C`divisors\*(C'\fR and Mathematica's \f(CW\*(C`Divisors[n]\*(C'\fR functions. .PP In scalar context this returns the sigma0 function, the sigma function (see Hardy and Wright section 16.7, or \s-1OEIS A000203\s0). This is the same result as evaluating the array in scalar context. .PP Also see the \*(L"for_divisors\*(R" functions for looping over the divisors. .SS "trial_factor" .IX Subsection "trial_factor" .Vb 1 \& my @factors = trial_factor($n); .Ve .PP 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 \f(CW*_factor\fR routines, this is not exported unless explicitly requested. .SS "fermat_factor" .IX Subsection "fermat_factor" .Vb 1 \& my @factors = fermat_factor($n); .Ve .PP 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). .SS "holf_factor" .IX Subsection "holf_factor" .Vb 1 \& my @factors = holf_factor($n); .Ve .PP 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. .SS "squfof_factor" .IX Subsection "squfof_factor" .Vb 1 \& my @factors = squfof_factor($n); .Ve .PP 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. .SS "prho_factor" .IX Subsection "prho_factor" .SS "pbrent_factor" .IX Subsection "pbrent_factor" .Vb 2 \& my @factors = prho_factor($n); \& my @factors = pbrent_factor($n); \& \& # Use a very small number of rounds \& my @factors = prho_factor($n, 1000); .Ve .PP 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. .SS "pminus1_factor" .IX Subsection "pminus1_factor" .Vb 3 \& 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 .Ve .PP Produces factors, not necessarily prime, of the positive number input. This is Pollard's \f(CW\*(C`p\-1\*(C'\fR method, using two stages with default smoothness settings of 1_000_000 for B1, and \f(CW\*(C`10 * B1\*(C'\fR for B2. This method can rapidly find a factor \f(CW\*(C`p\*(C'\fR of \f(CW\*(C`n\*(C'\fR where \f(CW\*(C`p\-1\*(C'\fR is smooth (it has no large factors). .SS "pplus1_factor" .IX Subsection "pplus1_factor" .Vb 2 \& my @factors = pplus1_factor($n); \& my @factors = pplus1_factor($n, 1_000); # set B1 smoothness .Ve .PP Produces factors, not necessarily prime, of the positive number input. This is Williams' \f(CW\*(C`p+1\*(C'\fR method, using one stage and two predefined initial points. .SS "ecm_factor" .IX Subsection "ecm_factor" .Vb 2 \& my @factors = ecm_factor($n); \& my @factors = ecm_factor($n, 100, 400, 10); # B1, B2, # of curves .Ve .PP Produces factors, not necessarily prime, of the positive number input. This is the elliptic curve method using two stages. .SH "MATHEMATICAL FUNCTIONS" .IX Header "MATHEMATICAL FUNCTIONS" .SS "ExponentialIntegral" .IX Subsection "ExponentialIntegral" .Vb 1 \& my $Ei = ExponentialIntegral($x); .Ve .PP Given a non-zero floating point input \f(CW\*(C`x\*(C'\fR, this returns the real-valued exponential integral of \f(CW\*(C`x\*(C'\fR, defined as the integral of \f(CW\*(C`e^t/t dt\*(C'\fR from \f(CW\*(C`\-infinity\*(C'\fR to \f(CW\*(C`x\*(C'\fR. .PP If the bignum module has been loaded, all inputs will be treated as if they were Math::BigFloat objects. .PP For non\-BigInt/BigFloat objects, the result should be accurate to at least 14 digits. .PP For BigInt / BigFloat objects, we first check to see if Math::MPFR is available. If so, then it is used since it is very fast and has high accuracy. Accuracy when using \s-1MPFR\s0 will be equal to the \f(CW\*(C`accuracy()\*(C'\fR value of the input (or the default BigFloat accuracy, which is 40 by default). .PP \&\s-1MPFR\s0 is used for positive inputs only. If Math::MPFR is not available or the input is negative, then other methods are used: continued fractions (\f(CW\*(C`x < \-1\*(C'\fR), rational Chebyshev approximation (\f(CW\*(C` \-1 < x < 0\*(C'\fR), a convergent series (small positive \f(CW\*(C`x\*(C'\fR), or an asymptotic divergent series (large positive \f(CW\*(C`x\*(C'\fR). Accuracy should be at least 14 digits. .SS "LogarithmicIntegral" .IX Subsection "LogarithmicIntegral" .Vb 1 \& my $li = LogarithmicIntegral($x) .Ve .PP Given a positive floating point input, returns the floating point logarithmic integral of \f(CW\*(C`x\*(C'\fR, defined as the integral of \f(CW\*(C`dt/ln t\*(C'\fR from \f(CW0\fR to \f(CW\*(C`x\*(C'\fR. If given a negative input, the function will croak. The function returns 0 at \f(CW\*(C`x = 0\*(C'\fR, and \f(CW\*(C`\-infinity\*(C'\fR at \f(CW\*(C`x = 1\*(C'\fR. .PP This is often known as \f(CWli(x)\fR. A related function is the offset logarithmic integral, sometimes known as \f(CWLi(x)\fR which avoids the singularity at 1. It may be defined as \f(CW\*(C`Li(x) = li(x) \- li(2)\*(C'\fR. Crandall and Pomerance use the term \f(CW\*(C`li0\*(C'\fR for this function, and define \f(CW\*(C`li(x) = Li0(x) \- li0(2)\*(C'\fR. Due to this terminology confusion, it is important to check which exact definition is being used. .PP If the bignum module has been loaded, all inputs will be treated as if they were Math::BigFloat objects. .PP For non\-BigInt/BigFloat objects, the result should be accurate to at least 14 digits. .PP For BigInt / BigFloat objects, we first check to see if Math::MPFR is available. If so, then it is used, as it will return results much faster and can be more accurate. Accuracy when using \s-1MPFR\s0 will be equal to the \&\f(CW\*(C`accuracy()\*(C'\fR value of the input (or the default BigFloat accuracy, which is 40 by default). .PP \&\s-1MPFR\s0 is used for inputs greater than 1 only. If Math::MPFR is not installed or the input is less than 1, results will be calculated as \&\f(CW\*(C`Ei(ln x)\*(C'\fR. .SS "RiemannZeta" .IX Subsection "RiemannZeta" .Vb 1 \& my $z = RiemannZeta($s); .Ve .PP Given a floating point input \f(CW\*(C`s\*(C'\fR where \f(CW\*(C`s >= 0\*(C'\fR, returns the floating point value of X(s)\-1, where X(s) is the Riemann zeta function. One is subtracted to ensure maximum precision for large values of \f(CW\*(C`s\*(C'\fR. The zeta function is the sum from k=1 to infinity of \f(CW\*(C`1 / k^s\*(C'\fR. This function only uses real arguments, so is basically the Euler Zeta function. .PP If the bignum module has been loaded, all inputs will be treated as if they were Math::BigFloat objects. .PP For non\-BigInt/BigFloat objects, the result should be accurate to at least 14 digits. The \s-1XS\s0 code uses a rational Chebyshev approximation between 0.5 and 5, and a series for other values. The \s-1PP\s0 code uses an identical series for all values. .PP For BigInt / BigFloat objects, we first check to see if the Math::MPFR module is installed. If so, then it is used, as it will return results much faster and can be more accurate. Accuracy when using \s-1MPFR\s0 will be equal to the \&\f(CW\*(C`accuracy()\*(C'\fR value of the input (or the default BigFloat accuracy, which is 40 by default). .PP If Math::MPFR is not installed, then results are calculated using either Borwein (1991) algorithm 2, or the basic series. Full input accuracy is attempted, but Math::BigFloat \&\s-1RT 43692\s0 produces incorrect high-accuracy computations without the fix. It is also very slow. I highly recommend installing Math::MPFR for BigFloat computations. .SS "RiemannR" .IX Subsection "RiemannR" .Vb 1 \& my $r = RiemannR($x); .Ve .PP 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. .PP If the bignum module has been loaded, all inputs will be treated as if they were Math::BigFloat objects. .PP For non\-BigInt/BigFloat objects, the result should be accurate to at least 14 digits. .PP For BigInt / BigFloat objects, we first check to see if the Math::MPFR module is installed. If so, then it is used, as it will return results much faster and can be more accurate. Accuracy when using \s-1MPFR\s0 will be equal to the \&\f(CW\*(C`accuracy()\*(C'\fR value of the input (or the default BigFloat accuracy, which is 40 by default). Accuracy without \s-1MPFR\s0 should be 35 digits. .SS "LambertW" .IX Subsection "LambertW" Returns the principal branch of the Lambert W function of a real value. Given a value \f(CW\*(C`k\*(C'\fR this solves for \f(CW\*(C`W\*(C'\fR in the equation \f(CW\*(C`k = We^W\*(C'\fR. The input must not be less than \f(CW\*(C`\-1/e\*(C'\fR. This corresponds to Pari's \f(CW\*(C`lambertw\*(C'\fR function and Mathematica's \f(CW\*(C`LambertW\*(C'\fR function. .SS "Pi" .IX Subsection "Pi" .Vb 2 \& my $tau = 2 * Pi; # $tau = 6.28318530717959 \& my $tau = 2 * Pi(40); # $tau = 6.283185307179586476925286766559005768394 .Ve .PP With no arguments, returns the value of Pi as an \s-1NV. \s0 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 \s-1NV\s0 if the number of digits fits in an \s-1NV \s0(typically 15 or less), or a Math::BigFloat object otherwise. .PP For sizes over 10k digits, having one of Math::MPFR, Math::Prime::Util::GMP, or Math::BigInt::GMP installed will help performance. For sizes over 50k one of the first two are highly recommended. .SH "EXAMPLES" .IX Header "EXAMPLES" Print strong pseudoprimes to base 17 up to 10M: .PP .Vb 2 \& # Similar to A001262\*(Aqs isStrongPsp function, but much faster \& perl \-MMath::Prime::Util=:all \-E \*(Aqforcomposites { say if is_strong_pseudoprime($_,17) } 10000000;\*(Aq .Ve .PP Print some primes above 64\-bit range: .PP .Vb 1 \& perl \-MMath::Prime::Util=:all \-Mbigint \-E \*(Aqmy $start=100000000000000000000; say join "\en", @{primes($start,$start+1000)}\*(Aq \& \& # Another way \& perl \-MMath::Prime::Util=:all \-E \*(Aqforprimes { say } "100000000000000000039", "100000000000000000993"\*(Aq \& \& # Similar using Math::Pari: \& # perl \-MMath::Pari=:int,PARI,nextprime \-E \*(Aqmy $start = PARI "100000000000000000000"; my $end = $start+1000; my $p=nextprime($start); while ($p <= $end) { say $p; $p = nextprime($p+1); }\*(Aq .Ve .PP Generate Carmichael numbers (\s-1OEIS A002997\s0 ): .PP .Vb 1 \& perl \-MMath::Prime::Util=:all \-E \*(Aqforoddcomposites { say if $_ % carmichael_lambda($_) == 1 } 1e6;\*(Aq .Ve .PP Examining the X3(x) function of Planat and Sole\*' (2011): .PP .Vb 12 \& 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); .Ve .PP Construct and use a Sophie-Germain prime iterator: .PP .Vb 10 \& 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\->(), "\en" for 1 .. 10000; .Ve .PP Project Euler, problem 3 (Largest prime factor): .PP .Vb 3 \& use Math::Prime::Util qw/factor/; \& use bigint; # Only necessary for 32\-bit machines. \& say 0+(factor(600851475143))[\-1] .Ve .PP Project Euler, problem 7 (10001st prime): .PP .Vb 2 \& use Math::Prime::Util qw/nth_prime/; \& say nth_prime(10_001); .Ve .PP Project Euler, problem 10 (summation of primes): .PP .Vb 7 \& use Math::Prime::Util qw/forprimes/; \& my $sum = 0; \& forprimes { $sum += $_ } 2_000_000; \& say $sum; \& # Or (fine for this size, not good for very large values) \& use Math::Prime::Util qw/vecsum primes/; \& say vecsum( @{primes(2_000_000)} ); .Ve .PP Project Euler, problem 21 (Amicable numbers): .PP .Vb 10 \& 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 ); .Ve .PP Project Euler, problem 41 (Pandigital prime), brute force command line: .PP .Vb 1 \& perl \-MMath::Prime::Util=primes \-MList::Util=first \-E \*(Aqsay first { /1/&&/2/&&/3/&&/4/&&/5/&&/6/&&/7/} reverse @{primes(1000000,9999999)};\*(Aq .Ve .PP Project Euler, problem 47 (Distinct primes factors): .PP .Vb 5 \& 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; .Ve .PP Project Euler, problem 69, stupid brute force solution (about 1 second): .PP .Vb 7 \& 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"; .Ve .PP Here is the right way to do \s-1PE\s0 problem 69 (under 0.03s): .PP .Vb 4 \& use Math::Prime::Util qw/pn_primorial/; \& my $n = 0; \& $n++ while pn_primorial($n+1) < 1000000; \& say pn_primorial($n); .Ve .PP Project Euler, problem 187, stupid brute force solution, 2 to 3 minutes: .PP .Vb 4 \& use Math::Prime::Util qw/forcomposites factor/; \& my $nsemis = 0; \& forcomposites { $nsemis++ if scalar factor($_) == 2; } int(10**8)\-1; \& say $nsemis; .Ve .PP Here is one of the best ways for \s-1PE187: \s0 under 20 milliseconds from the command line. Much faster than Pari, and competitive with Mathematica. .PP .Vb 8 \& 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; .Ve .PP To get the result of \*(L"matches\*(R" in Math::Factor::XS: .PP .Vb 7 \& 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]X$_\->[1]" } matches($n)); .Ve .PP or its \f(CW\*(C`matches\*(C'\fR function with the \f(CW\*(C`skip_multiples\*(C'\fR option: .PP .Vb 6 \& 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; } \& } .Ve .PP Compute \s-1OEIS A054903\s0 just like CRG4s Pari example: .PP .Vb 4 \& use Math::Prime::Util qw/forcomposite divisor_sum/; \& forcomposites { \& say if divisor_sum($_)+6 == divisor_sum($_+6) \& } 9,1e7; .Ve .PP Construct the table shown in \s-1OEIS A046147\s0 : .PP .Vb 10 \& 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); \& } \& } .Ve .PP Find the 7\-digit palindromic primes in the first 20k digits of Pi: .PP .Vb 6 \& 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])(\ed)(\ed)\ed\e4\e3\e2)/g) { \& say "$1 at ",pos($pi)\-7 if is_prime($1) \& } .Ve .PP The Bell numbers B_n: .PP .Vb 2 \& sub B { my $n = shift; vecsum(map { stirling($n,$_,2) } 0..$n) } \& say "$_ ",B($_) for 1..50; .Ve .SH "PRIMALITY TESTING NOTES" .IX Header "PRIMALITY TESTING NOTES" Above \f(CW\*(C`2^64\*(C'\fR, \*(L"is_prob_prime\*(R" performs an extra-strong \&\s-1BPSW\s0 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, \s-1SAGE, FLINT,\s0 etc., then this function should be appropriate for you. \*(L"is_prime\*(R" will do the same \s-1BPSW\s0 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. \*(L"is_provable_prime\*(R" constructs a primality proof. If a certificate is requested, then either \s-1BLS75\s0 theorem 5 or \s-1ECPP\s0 is performed. Without a certificate, the method is implementation specific (currently it is identical, but later releases may use \s-1APRCL\s0). With Math::Prime::Util::GMP installed, this is quite fast through 300 or so digits. .PP Math systems 30 years ago typically used Miller-Rabin tests with \f(CW\*(C`k\*(C'\fR bases (usually fixed bases, sometimes random) for primality testing, but these have generally been replaced by some form of \s-1BPSW\s0 as used in this module. See Pinch's 1993 paper for examples of why using \f(CW\*(C`k\*(C'\fR M\-R tests leads to poor results. The three exceptions in common contemporary use I am aware of are: .IP "libtommath" 4 .IX Item "libtommath" Uses the first \f(CW\*(C`k\*(C'\fR 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 \s-1CPAN\s0 Crypto module. .IP "\s-1GMP/MPIR\s0" 4 .IX Item "GMP/MPIR" Uses a set of \f(CW\*(C`k\*(C'\fR static-random bases. The bases are randomly chosen using a \s-1PRNG\s0 that is seeded identically each call (the seed changes with each release). This offers a very slight advantage over using the first \f(CW\*(C`k\*(C'\fR prime bases, but not much. See, for example, Nicely's mpz_probab_prime_p pseudoprimes page. .IP "Math::Pari (not recent Pari/GP)" 4 .IX Item "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 \s-1PRNG\s0 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 \s-1BPSW / APRCL,\s0 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. .PP Basically the problem is that it is just too easy to get counterexamples from running \f(CW\*(C`k\*(C'\fR M\-R tests, forcing one to use a very large number of tests (at least 20) to avoid frequent false results. Using the \s-1BPSW\s0 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 \fIstill\fR be much faster. .PP Using \f(CW\*(C`k\*(C'\fR 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 \*(L"pick my number\*(R" 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. .SH "LIMITATIONS" .IX Header "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 \&\f(CW\*(C`~ 2^49\*(C'\fR). The best solution is updating to a more recent Perl. .PP 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 \f(CW\*(C`prime_precalc\*(C'\fR is called before using \f(CW\*(C`primes\*(C'\fR, \&\f(CW\*(C`prime_count\*(C'\fR, or \f(CW\*(C`nth_prime\*(C'\fR with large inputs. This is \fBonly\fR an issue if you use non-Cygwin Win32 \fBand\fR call these routines from within Perl threads. .PP Because the loop functions like \*(L"forprimes\*(R" use \f(CW\*(C`MULTICALL\*(C'\fR, there is some odd behavior with anonymous sub creation inside the block. This is shared with most \s-1XS\s0 modules that use \f(CW\*(C`MULTICALL\*(C'\fR, and is rarely seen because it is such an unusual use. An example is: .PP .Vb 2 \& forprimes { my $var = "p is $_"; push @subs, sub {say $var}; } 50; \& $_\->() for @subs; .Ve .PP This can be worked around by using double braces for the function, e.g. \&\f(CW\*(C`forprimes {{ ... }} 50\*(C'\fR. .SH "SEE ALSO" .IX Header "SEE ALSO" This section describes other \s-1CPAN\s0 modules available that have some feature overlap with this one. Also see the \*(L"\s-1REFERENCES\*(R"\s0 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. .PP I will use SoE to indicate the Sieve of Eratosthenes, and \s-1MPU\s0 to denote this module (Math::Prime::Util). Some quick alternatives I can recommend if you don't want to use \s-1MPU:\s0 .IP "\(bu" 4 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. .IP "\(bu" 4 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 \fIvery\fR slow. .IP "\(bu" 4 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). .PP Math::Prime::XS has \f(CW\*(C`is_prime\*(C'\fR and \f(CW\*(C`primes\*(C'\fR functionality. There is no bigint support. The \f(CW\*(C`is_prime\*(C'\fR function uses well-written trial division, meaning it is very fast for small numbers, but terribly slow for large 64\-bit numbers. \s-1MPU\s0 is similarly fast with small numbers, but becomes faster as the size increases. \&\s-1MPXS\s0's prime sieve is an unoptimized non-segmented SoE which returns an array. Sieve bases larger than \f(CW\*(C`10^7\*(C'\fR start taking inordinately long and using a lot of memory (gigabytes beyond \f(CW\*(C`10^10\*(C'\fR). E.g. \f(CW\*(C`primes(10**9, 10**9+1000)\*(C'\fR takes 36 seconds with \s-1MPXS,\s0 but only 0.0001 seconds with \s-1MPU.\s0 .PP Math::Prime::FastSieve supports \f(CW\*(C`primes\*(C'\fR, \f(CW\*(C`is_prime\*(C'\fR, \f(CW\*(C`next_prime\*(C'\fR, \&\f(CW\*(C`prev_prime\*(C'\fR, \f(CW\*(C`prime_count\*(C'\fR, and \f(CW\*(C`nth_prime\*(C'\fR. The caveat is that all functions only work within the sieved range, so are limited to about \f(CW\*(C`10^10\*(C'\fR. It uses a fast SoE to generate the main sieve. The sieve is 2\-3x slower than the base sieve for \s-1MPU,\s0 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 \s-1MPU\s0 using \&\f(CW\*(C`prime_precalc\*(C'\fR but is not required. .PP Bit::Vector supports the \f(CW\*(C`primes\*(C'\fR and \f(CW\*(C`prime_count\*(C'\fR functionality in a somewhat similar way to Math::Prime::FastSieve. It is the slowest of all the \s-1XS\s0 sieves, and has the most memory use. It is faster than pure Perl code. .PP Crypt::Primes supports \f(CW\*(C`random_maurer_prime\*(C'\fR functionality. \s-1MPU\s0 has more options for random primes (n\-digit, n\-bit, ranged, and strong) in addition to Maurer's algorithm. \s-1MPU\s0 does not have the critical bug \&\s-1RT81858\s0 . \s-1MPU\s0 has a more uniform distribution as well as return a larger subset of primes (\s-1RT81871\s0 ). \&\s-1MPU\s0 does not depend on Math::Pari though can run slow for bigints unless the Math::BigInt::GMP or Math::BigInt::Pari modules are installed. Having Math::Prime::Util::GMP installed also helps performance for \s-1MPU.\s0 Crypt::Primes is hardcoded to use Crypt::Random, while \s-1MPU\s0 uses Bytes::Random::Secure, and also allows plugging in a random function. This is more flexible, faster, has fewer dependencies, and uses a \s-1CSPRNG\s0 for security. \s-1MPU\s0 can return a primality certificate. What Crypt::Primes has that \s-1MPU\s0 does not is the ability to return a generator. .PP Math::Factor::XS calculates prime factors and factors, which correspond to the \*(L"factor\*(R" and \*(L"divisors\*(R" functions of \s-1MPU. \s0 These 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 \f(CW\*(C`count_prime_factors\*(C'\fR can be done in \s-1MPU\s0 using \f(CW\*(C`scalar factor($n)\*(C'\fR. See the \*(L"\s-1EXAMPLES\*(R"\s0 section for a 2\-line function replicating \f(CW\*(C`matches\*(C'\fR. .PP Math::Big version 1.12 includes \f(CW\*(C`primes\*(C'\fR functionality. The current code is only usable for very tiny inputs as it is incredibly slow and uses lots of memory. \s-1RT81986\s0 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 \s-1MPU.\s0 .PP 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 \s-1MPU\s0 or Math::Factor::XS. 19\-digit semiprimes will take \fIhours\fR versus \s-1MPU\s0's single milliseconds. .PP Math::Factoring is a placeholder module for bigint factoring. Version 0.02 only supports trial division (the Pollard-Rho method does not work). .PP Math::Prime::TiedArray allows random access to a tied primes array, almost identically to what \s-1MPU\s0 provides in Math::Prime::Util::PrimeArray. \s-1MPU\s0 has attempted to fix Math::Prime::TiedArray's shift bug (\s-1RT58151\s0 ). \s-1MPU\s0 is typically much faster and will use less memory, but there are some cases where \&\s-1MP:TA\s0 is faster (\s-1MP:TA\s0 stores all entries up to the largest request, while \&\s-1MPU:PA\s0 stores only a window around the last request). .PP List::Gen is very interesting and includes a built-in primes iterator as well as a \f(CW\*(C`is_prime\*(C'\fR filter for arbitrary sequences. Unfortunately both are very slow. .PP Math::Primality supports \f(CW\*(C`is_prime\*(C'\fR, \f(CW\*(C`is_pseudoprime\*(C'\fR, \&\f(CW\*(C`is_strong_pseudoprime\*(C'\fR, \f(CW\*(C`is_strong_lucas_pseudoprime\*(C'\fR, \f(CW\*(C`next_prime\*(C'\fR, \&\f(CW\*(C`prev_prime\*(C'\fR, \f(CW\*(C`prime_count\*(C'\fR, and \f(CW\*(C`is_aks_prime\*(C'\fR functionality. This is a great little module that implements primality functionality. It was the first \s-1CPAN\s0 module to support the \s-1BPSW\s0 test. All inputs are processed using \s-1GMP,\s0 so it of course supports bigints. In fact, Math::Primality was made originally with bigints in mind, while \s-1MPU\s0 was originally targeted to native integers, but both have added better support for the other. The main differences are extra functionality (\s-1MPU\s0 has more functions) and performance. With native integer inputs, \s-1MPU\s0 is generally much faster, especially with \*(L"prime_count\*(R". For bigints, \&\s-1MPU\s0 is slower unless the Math::Prime::Util::GMP module is installed, in which case \s-1MPU\s0 is ~2x faster. Math::Primality also installs a \f(CW\*(C`primes.pl\*(C'\fR program, but it has much less functionality than the one included with \s-1MPU.\s0 .PP Math::NumSeq does not have a one-to-one mapping between functions in \s-1MPU,\s0 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 \fIlot\fR of sequences. For those methods that overlap, \s-1MPU\s0 is usually much faster. Importantly, most of the sequences in Math::NumSeq are limited to 32\-bit indices. .PP \&\*(L"cr_combine\*(R" in Math::ModInt::ChineseRemainder is similar to \s-1MPU\s0's \*(L"chinese\*(R", and in fact they use the same algorithm. The former module uses caching of moduli to speed up further operations. \s-1MPU\s0 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). .PP For combinations and permutations there are many alternatives. One difference with nearly all of them is that \s-1MPU\s0's \*(L"forcomb\*(R" and \&\*(L"forperm\*(R" 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 \s-1MPU \s0(they use the same user-block structure but twiddle the user array each call). .PP Math::Pari supports a lot of features, with a great deal of overlap. In general, \s-1MPU\s0 will be faster for native 64\-bit integers, while it's differs for bigints (Pari will always be faster if Math::Prime::Util::GMP is not installed; with it, it varies by function). Note that Pari extends many of these functions to other spaces (Gaussian integers, complex numbers, vectors, matrices, polynomials, etc.) which are beyond the realm of this module. Some of the highlights: .ie n .IP """isprime""" 4 .el .IP "\f(CWisprime\fR" 4 .IX Item "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 \s-1GMP\s0's \f(CW\*(C`is_probab_prime\*(C'\fR). This has a greater chance of false positives compared to the \s-1BPSW\s0 test \*(-- some composites such as \&\f(CW9\fR, \f(CW88831\fR, \f(CW38503\fR, etc. (\s-1OEIS A141768\s0 ) have a surprisingly high chance of being indicated prime. Using \f(CW\*(C`isprime($n,1)\*(C'\fR will perform an \f(CW\*(C`n\-1\*(C'\fR proof, but this becomes unreasonably slow past 70 or so digits. .Sp If Math::Pari is built using Pari 2.3.5 (this requires manual configuration) then the primality tests are completely different. Using \&\f(CW\*(C`ispseudoprime\*(C'\fR will perform a \s-1BPSW\s0 test and is quite a bit faster than the older test. \f(CW\*(C`isprime\*(C'\fR now does an APR-CL proof (fast, but no certificate). .Sp Math::Primality uses a strong \s-1BPSW\s0 test, which is the standard \s-1BPSW\s0 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 \s-1BPSW\s0 test for its \&\f(CW\*(C`ispseudoprime\*(C'\fR 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 \s-1BPSW\s0 test, which has an even lower chance of counterexample. With Math::Prime::Util::GMP, \f(CW\*(C`is_prime\*(C'\fR adds 1 to 5 extra M\-R tests using random bases, which further reduces the probability of a composite being allowed to pass. .ie n .IP """primepi""" 4 .el .IP "\f(CWprimepi\fR" 4 .IX Item "primepi" Only available with version 2.3 of Pari. Similar to \s-1MPU\s0's \*(L"prime_count\*(R" function in \s-1API,\s0 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 \s-1MPU.\s0 .ie n .IP """primes""" 4 .el .IP "\f(CWprimes\fR" 4 .IX Item "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 \s-1MPU,\s0 but doesn't support segmenting. .ie n .IP """factorint""" 4 .el .IP "\f(CWfactorint\fR" 4 .IX Item "factorint" Similar to \s-1MPU\s0's \*(L"factor_exp\*(R" though with a slightly different return. \&\s-1MPU\s0 offers \*(L"factor\*(R" for a linear array of prime factors where n = p1 * p2 * p3 * ... as (p1,p2,p3,...) and \*(L"factor_exp\*(R" 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 \s-1MPU\s0 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 \s-1MPU.\s0 .ie n .IP """divisors""" 4 .el .IP "\f(CWdivisors\fR" 4 .IX Item "divisors" Similar to \s-1MPU\s0's \*(L"divisors\*(R". .ie n .IP """forprime"", ""forcomposite"", ""fordiv"", ""sumdiv""" 4 .el .IP "\f(CWforprime\fR, \f(CWforcomposite\fR, \f(CWfordiv\fR, \f(CWsumdiv\fR" 4 .IX Item "forprime, forcomposite, fordiv, sumdiv" Similar to \s-1MPU\s0's \*(L"forprimes\*(R", \*(L"forcomposites\*(R", \*(L"fordivisors\*(R", and \&\*(L"divisor_sum\*(R". .ie n .IP """eulerphi"", ""moebius""" 4 .el .IP "\f(CWeulerphi\fR, \f(CWmoebius\fR" 4 .IX Item "eulerphi, moebius" Similar to \s-1MPU\s0's \*(L"euler_phi\*(R" and \*(L"moebius\*(R". \s-1MPU\s0 is 2\-20x faster for native integers. \s-1MPU\s0 also supported range inputs, which can be much more efficient. Without Math::Prime::Util::GMP installed, \s-1MPU\s0 is very slow with bigints. With it installed, it is about 2x slower than Math::Pari. .ie n .IP """gcd"", ""lcm"", ""kronecker"", ""znorder"", ""znprimroot"", ""znlog""" 4 .el .IP "\f(CWgcd\fR, \f(CWlcm\fR, \f(CWkronecker\fR, \f(CWznorder\fR, \f(CWznprimroot\fR, \f(CWznlog\fR" 4 .IX Item "gcd, lcm, kronecker, znorder, znprimroot, znlog" Similar to \s-1MPU\s0's \*(L"gcd\*(R", \*(L"lcm\*(R", \*(L"kronecker\*(R", \*(L"znorder\*(R", \&\*(L"znprimroot\*(R", and \*(L"znlog\*(R". Pari's \f(CW\*(C`znprimroot\*(C'\fR 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). \s-1MPU\s0's \*(L"znprimroot\*(R" will always return the smallest root if it exists, and \f(CW\*(C`undef\*(C'\fR otherwise. Similarly, \s-1MPU\s0's \*(L"znlog\*(R" will return the smallest \f(CW\*(C`x\*(C'\fR and work with non-primitive-root \f(CW\*(C`g\*(C'\fR, which is similar to Pari/GP 2.6, but not the older versions in Math::Pari. The performance of \*(L"znlog\*(R" is fairly good compared to older Pari/GP, but much worse than 2.6's new methods. .ie n .IP """sigma""" 4 .el .IP "\f(CWsigma\fR" 4 .IX Item "sigma" Similar to \s-1MPU\s0's \*(L"divisor_sum\*(R". \s-1MPU\s0 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. .ie n .IP """numbpart"", ""forpart""" 4 .el .IP "\f(CWnumbpart\fR, \f(CWforpart\fR" 4 .IX Item "numbpart, forpart" Similar to \s-1MPU\s0's \*(L"partitions\*(R" and \*(L"forpart\*(R". These functions were introduced in Pari 2.3 and 2.6, hence are not in Math::Pari. \f(CW\*(C`numbpart\*(C'\fR produce identical results to \f(CW\*(C`partitions\*(C'\fR, but Pari is \fImuch\fR faster. forpart is very similar to Pari's function, but produces a different ordering (\s-1MPU\s0 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. .ie n .IP """eint1""" 4 .el .IP "\f(CWeint1\fR" 4 .IX Item "eint1" Similar to \s-1MPU\s0's \*(L"ExponentialIntegral\*(R". .ie n .IP """zeta""" 4 .el .IP "\f(CWzeta\fR" 4 .IX Item "zeta" \&\s-1MPU\s0 has \*(L"RiemannZeta\*(R" which takes non-negative real inputs, while Pari's function supports negative and complex inputs. .PP 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 \s-1MPU. \s0 For bigints, Math::Pari may be superior and it rarely has any performance surprises. Some of the unique features \s-1MPU\s0 offers include super fast prime counts, nth_prime, \&\s-1ECPP\s0 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. .SH "PERFORMANCE" .IX Header "PERFORMANCE" First, for those looking for the state of the art non-Perl solutions: .IP "Primality testing" 4 .IX Item "Primality testing" For general numbers smaller than 2000 or so digits, \s-1MPU\s0 is the fastest solution I am aware of (it is faster than Pari 2.7, \s-1PFGW,\s0 and \s-1FLINT\s0). For very large inputs, \&\s-1PFGW\s0 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 \s-1BPSW\s0 test however, and there are quite a few counterexamples for a given base of its \s-1PRP\s0 test, so it is commonly used for fast filtering of large candidates. A test such as the \s-1BPSW\s0 test in this module is then recommended. .IP "Primality proofs" 4 .IX Item "Primality proofs" Primo 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 , the \s-1APRCL\s0 from the modern Pari package, or the standalone \s-1ECPP\s0 from this module with large polynomial set. .IP "Factoring" 4 .IX Item "Factoring" yafu , msieve , and gmp-ecm are all good choices for large inputs. The factoring code in this module (and all other \s-1CPAN\s0 modules) is very limited compared to those. .IP "Primes" 4 .IX Item "Primes" primesieve and 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. .Sp Note that the Sieve of Atkin is \fInot\fR 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 \fIprimegen\fR, which runs on par with the SoE in this module. The SoE's in Pari, yafu, and primesieve are all faster. .IP "Prime Counts and Nth Prime" 4 .IX Item "Prime Counts and Nth Prime" Outside of private research implementations doing prime counts for \&\f(CW\*(C`n > 2^64\*(C'\fR, this module should be close to state of the art in performance, and supports results up to \f(CW\*(C`2^64\*(C'\fR. Further performance improvements are planned, as well as expansion to larger values. .Sp 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. .SS "\s-1PRIME COUNTS\s0" .IX Subsection "PRIME COUNTS" Counting the primes to \f(CW\*(C`800_000_000\*(C'\fR (800 million): .PP .Vb 10 \& Time (s) Module Version Notes \& \-\-\-\-\-\-\-\-\- \-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\- \-\-\-\-\-\-\- \-\-\-\-\-\-\-\-\-\-\- \& 0.001 Math::Prime::Util 0.37 using extended LMO \& 0.007 Math::Prime::Util 0.12 using Lehmer\*(Aqs method \& 0.27 Math::Prime::Util 0.17 segmented mod\-30 sieve \& 0.39 Math::Prime::Util::PP 0.24 Perl (Lehmer\*(Aqs 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 .Ve .PP Python's standard modules are very slow: \s-1MPMATH\s0 v0.17 \f(CW\*(C`primepi\*(C'\fR takes 169.5s and 25+ \s-1GB\s0 of \s-1RAM. \s0 SymPy 0.7.1 \f(CW\*(C`primepi\*(C'\fR 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 \s-1NUMPY\s0 in 2.8s. .SS "\s-1PRIMALITY TESTING\s0" .IX Subsection "PRIMALITY TESTING" .IP "Small inputs: is_prime from 1 to 20M" 4 .IX Item "Small inputs: is_prime from 1 to 20M" .Vb 7 \& 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) .Ve .IP "Large native inputs: is_prime from 10^16 to 10^16 + 20M" 4 .IX Item "Large native inputs: is_prime from 10^16 to 10^16 + 20M" .Vb 5 \& 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. .Ve .IP "bigints: is_prime from 10^100 to 10^100 + 0.2M" 4 .IX Item "bigints: is_prime from 10^100 to 10^100 + 0.2M" .Vb 7 \& 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) .Ve .IP "\(bu" 4 \&\s-1MPU\s0 is consistently the fastest solution, and performs the most stringent probable prime tests on bigints. .IP "\(bu" 4 Math::Primality has a lot of overhead that makes it quite slow for native size integers. With bigints we finally see it work well. .IP "\(bu" 4 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. \s-1MPU\s0 is still faster, but clearly the time for native integers is dominated by the calling overhead. .SS "\s-1FACTORING\s0" .IX Subsection "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 \s-1GMP\s0 module, performance is pretty similar from 20 through 70 digits, which the caveat that the current \s-1MPU\s0 factoring uses more memory for 60+ digit numbers. .PP This slide presentation has a lot of data on 64\-bit and \s-1GMP\s0 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 \s-1SQUFOF\s0 is typically the method of choice for 11\-18 digits (I've seen claims that a lightweight \s-1QS\s0 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 , msieve , gmp-ecm , \&\s-1GGNFS\s0 , and Pari . The latest yafu should cover most uses, with \s-1GGNFS\s0 likely only providing a benefit for numbers large enough to warrant distributed processing. .SS "\s-1PRIMALITY PROVING\s0" .IX Subsection "PRIMALITY PROVING" The \f(CW\*(C`n\-1\*(C'\fR 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 \s-1MPU::GMP\s0 contain an \&\s-1ECPP\s0 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 . .SS "\s-1RANDOM PRIME GENERATION\s0" .IX Subsection "RANDOM PRIME GENERATION" Seconds per prime for random prime generation on a circa\-2009 workstation, with Math::BigInt::GMP, Math::Prime::Util::GMP, and Math::Random::ISAAC::XS installed. .PP .Vb 10 \& bits random +testing rand_prov Maurer Shw\-Tylr CPMaurer \& \-\-\-\-\- \-\-\-\-\-\-\-\- \-\-\-\-\-\-\-\- \-\-\-\-\-\-\-\-\- \-\-\-\-\-\-\-\- \-\-\-\-\-\-\-\- \-\-\-\-\-\-\-\- \& 64 0.0001 +0.000008 0.0002 0.0001 0.010 0.022 \& 128 0.0020 +0.00023 0.011 0.063 0.028 0.057 \& 256 0.0034 +0.0004 0.058 0.13 0.042 0.16 \& 512 0.0097 +0.0012 0.28 0.28 0.085 0.41 \& 1024 0.060 +0.0060 0.65 0.65 0.24 2.19 \& 2048 0.57 +0.039 4.8 4.8 1.0 10.99 \& 4096 6.24 +0.25 31.9 31.9 8.2 79.71 \& 8192 58.6 +1.61 234.0 234.0 112.9 947.3 \& \& random = random_nbit_prime (results pass BPSW) \& random+ = additional time for 3 M\-R and a Frobenius test \& rand_prov = random_proven_prime \& maurer = random_maurer_prime \& Shw\-Tylr = random_shawe_taylor_prime \& CPMaurer = Crypt::Primes::maurer .Ve .PP \&\*(L"random_nbit_prime\*(R" is reasonably fast, and for most purposes should suffice. If good uniformity isn't important, the \f(CW\*(C`use_primeinc\*(C'\fR config option can be set and double the speed. 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 \s-1BPSW,\s0 three more M\-R tests, and a Frobenius test is \fIextraordinarily\fR small. .PP \&\*(L"random_proven_prime\*(R" provides a randomly selected prime with an optional certificate, without specifying the particular method. Below 512 bits, using \*(L"is_provable_prime\*(R"(\*(L"random_nbit_prime\*(R") is typically faster than Maurer's algorithm, but becomes quite slow as the bit size increases. This leaves the decision of the exact method of proving the result to the implementation. .PP \&\*(L"random_maurer_prime\*(R" 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. .PP \&\*(L"random_shawe_taylor_prime\*(R" similarly constructs a provable prime. It uses a simpler construction method. The implementation uses a single large random seed followed by \s-1SHA\-256\s0 as specified by \s-1FIPS 186\-4. \s0 As seen, it is a bit faster than the Maurer implementation. .PP \&\*(L"maurer\*(R" in Crypt::Primes times are included for comparison. It is pretty fast for small sizes but gets slow as the size increases. It does not perform any primality checks on the intermediate results or the final result (I highly recommended you run a primality test on the output). Additionally important for servers, \*(L"maurer\*(R" in Crypt::Primes uses excessive system entropy and can grind to a halt if \f(CW\*(C`/dev/random\*(C'\fR is exhausted (it can take \fBdays\fR to return). The times above are on a machine running \&\s-1HAVEGED\s0 so never waits for entropy. Without this, the times would be much higher. .SH "AUTHORS" .IX Header "AUTHORS" Dana Jacobsen .SH "ACKNOWLEDGEMENTS" .IX Header "ACKNOWLEDGEMENTS" Eratosthenes of Cyrene provided the elegant and simple algorithm for finding primes. .PP Terje Mathisen, A.R. Quesada, and B. Van Pelt all had useful ideas which I used in my wheel sieve. .PP The \s-1SQUFOF\s0 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. .PP The \s-1LMO\s0 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. .SH "REFERENCES" .IX Header "REFERENCES" .IP "\(bu" 4 Christian Bau, \*(L"The Extended Meissel-Lehmer Algorithm\*(R", 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 \s-1LMO. \s0 .IP "\(bu" 4 Manuel Benito and Juan L. Varona, \*(L"Recursive formulas related to the summation of the Mo\*:bius function\*(R", \fIThe Open Mathematics Journal\fR, v1, pp 25\-34, 2007. Among many other things, shows a simple formula for computing the Mertens functions with only n/3 Mo\*:bius values (not as fast as Dele\*'glise and Rivat, but really simple). .IP "\(bu" 4 John Brillhart, D. H. Lehmer, and J. L. Selfridge, \*(L"New Primality Criteria and Factorizations of 2^m +/\- 1\*(R", Mathematics of Computation, v29, n130, Apr 1975, pp 620\-647. .IP "\(bu" 4 W. J. Cody and Henry C. Thacher, Jr., \*(L"Rational Chebyshev Approximations for the Exponential Integral E_1(x)\*(R", \fIMathematics of Computation\fR, v22, pp 641\-649, 1968. .IP "\(bu" 4 W. J. Cody and Henry C. Thacher, Jr., \*(L"Chebyshev approximations for the exponential integral Ei(x)\*(R", \fIMathematics of Computation\fR, v23, pp 289\-303, 1969. .IP "\(bu" 4 W. J. Cody, K. E. Hillstrom, and Henry C. Thacher Jr., \*(L"Chebyshev Approximations for the Riemann Zeta Function\*(R", \*(L"Mathematics of Computation\*(R", v25, n115, pp 537\-547, July 1971. .IP "\(bu" 4 Henri Cohen, \*(L"A Course in Computational Algebraic Number Theory\*(R", Springer, 1996. Practical computational number theory from the team lead of Pari . Lots of explicit algorithms. .IP "\(bu" 4 Marc Dele\*'glise and Joo\*:l Rivat, \*(L"Computing the summation of the Mo\*:bius function\*(R", \fIExperimental Mathematics\fR, v5, n4, pp 291\-295, 1996. Enhances the Mo\*:bius computation in Lioen/van de Lune, and gives a very efficient way to compute the Mertens function. .IP "\(bu" 4 Pierre Dusart, \*(L"Autour de la fonction qui compte le nombre de nombres premiers\*(R", PhD thesis, 1998. In French. The mathematics is readable and highly recommended reading if you're interested in prime number bounds. .IP "\(bu" 4 Pierre Dusart, \*(L"Estimates of Some Functions Over Primes without R.H.\*(R", preprint, 2010. Updates to the best non-RH bounds for prime count and nth prime. .IP "\(bu" 4 Pierre-Alain Fouque and Mehdi Tibouchi, \*(L"Close to Uniform Prime Number Generation With Fewer Random Bits\*(R", 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. .IP "\(bu" 4 Walter M. Lioen and Jan van de Lune, \*(L"Systematic Computations on Mertens' Conjecture and Dirichlet's Divisor Problem by Vectorized Sieving\*(R", in \fIFrom Universal Morphisms to Megabytes\fR, Centrum voor Wiskunde en Informatica, pp. 421\-432, 1994. Describes a nice way to compute a range of Mo\*:bius values. .IP "\(bu" 4 Ueli M. Maurer, \*(L"Fast Generation of Prime Numbers and Secure Public-Key Cryptographic Parameters\*(R", 1995. Generating random provable primes by building up the prime. .IP "\(bu" 4 Gabriel Mincu, \*(L"An Asymptotic Expansion\*(R", \fIJournal of Inequalities in Pure and Applied Mathematics\fR, v4, n2, 2003. A very readable account of Cipolla's 1902 nth prime approximation. .IP "\(bu" 4 \&\s-1OEIS:\s0 Primorial .IP "\(bu" 4 Vincent Pegoraro and Philipp Slusallek, \*(L"On the Evaluation of the Complex-Valued Exponential Integral\*(R", \fIJournal of Graphics, \s-1GPU,\s0 and Game Tools\fR, v15, n3, pp 183\-198, 2011. .IP "\(bu" 4 William H. Press et al., \*(L"Numerical Recipes\*(R", 3rd edition. .IP "\(bu" 4 Hans Riesel, \*(L"Prime Numbers and Computer Methods for Factorization\*(R", Birkh?user, 2nd edition, 1994. Lots of information, some code, easy to follow. .IP "\(bu" 4 David M. Smith, \*(L"Multiple-Precision Exponential Integral and Related Functions\*(R", \fI\s-1ACM\s0 Transactions on Mathematical Software\fR, v37, n4, 2011. .IP "\(bu" 4 Douglas A. Stoll and Patrick Demichel , \*(L"The impact of X(s) complex zeros on X(x) for x < 10^{10^{13}}\*(R", \*(L"Mathematics of Computation\*(R", v80, n276, pp 2381\-2394, October 2011. .SH "COPYRIGHT" .IX Header "COPYRIGHT" Copyright 2011\-2014 by Dana Jacobsen .PP This program is free software; you can redistribute it and/or modify it under the same terms as Perl itself.