Math::P...il::GMP(3) User Contributed Perl Documentation Math::P...il::GMP(3) NAME Math::Prime::Util::GMP - Utilities related to prime numbers and factoring, using GMP VERSION Version 0.53 SYNOPSIS use Math::Prime::Util::GMP ':all'; my $n = "115792089237316195423570985008687907853269984665640564039457584007913129639937"; # This doesn't impact the operation of the module at all, but does let you # enter big number arguments directly as well as enter (e.g.): 2**2048 + 1. use bigint; # These return 0 for composite, 2 for prime, and 1 for probably prime # Numbers under 2^64 will return 0 or 2. # is_prob_prime does a BPSW primality test for numbers > 2^64 # is_prime adds some MR tests and a quick test to try to prove the result # is_provable_prime will spend a lot of effort on proving primality say "$n is probably prime" if is_prob_prime($n); say "$n is ", qw(composite prob_prime def_prime)[is_prime($n)]; say "$n is definitely prime" if is_provable_prime($n) == 2; # Miller-Rabin and strong Lucas-Selfridge pseudoprime tests say "$n is a prime or spsp-2/7/61" if is_strong_pseudoprime($n, 2, 7, 61); 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); # Return array reference to primes in a range. my $aref = primes( 10 ** 200, 10 ** 200 + 10000 ); $next = next_prime($n); # next prime > n $prev = prev_prime($n); # previous prime < n # Primorials and lcm say "23# is ", primorial(23); say "The product of the first 47 primes is ", pn_primorial(47); say "lcm(1..1000) is ", consecutive_integer_lcm(1000); # Find prime factors of big numbers @factors = factor(5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497); # Finer control over factoring. # These stop after finding one factor or exceeding their limit. # # optional arguments o1, o2, ... @factors = trial_factor($n); # test up to o1 @factors = prho_factor($n); # no more than o1 rounds @factors = pbrent_factor($n); # no more than o1 rounds @factors = holf_factor($n); # no more than o1 rounds @factors = squfof_factor($n); # no more than o1 rounds @factors = pminus1_factor($n); # o1 = smoothness limit, o2 = stage 2 limit @factors = ecm_factor($n); # o1 = B1, o2 = # of curves @factors = qs_factor($n); # (no arguments) DESCRIPTION A module for number theory in Perl using GMP. This includes primality tests, getting primes in a range, factoring, and more. While it certainly can be used directly, the main purpose of this module is for Math::Prime::Util. That module will automatically load this one if it is installed, greatly speeding up many of its operations on big numbers. Inputs and outputs for big numbers are via strings, so you do not need to use a bigint package in your program. However if you do use bigints, in- puts will be converted internally so there is no need to convert before a call. Output results are returned as either Perl scalars (for native-size) or strings (for bigints). Math::Prime::Util tries to reconvert all strings back into the callers bigint type if possible, which makes it more conve- nient for calculations. The various "is_*_pseudoprime" tests are more appropriately called "is_*_probable_prime" or "is_*_prp". They return 1 if the input is a prob- able prime based on their test. The naming convention is historical and follows Pari, Math::Primality, and some other math packages. The modern definition of pseudoprime is a composite that passes the test, rather than any number. FUNCTIONS is_prob_prime my $prob_prime = is_prob_prime($n); # Returns 0 (composite), 2 (prime), or 1 (probably prime) Takes a positive number as input and returns back either 0 (composite), 2 (definitely prime), or 1 (probably prime). For inputs below "2^64" the test is deterministic, so the possible return values are 0 (composite) or 2 (definitely prime). For inputs above "2^64", a probabilistic test is performed. Only 0 (com- posite) and 1 (probably prime) are returned. The current implementation uses the Baillie-PSW (BPSW) test. There is a possibility that composites may be returned marked prime, but since the test was published in 1980, not a single BPSW pseudoprime has been found, so it is extremely likely to be prime. While we believe (Pomerance 1984) that an infinite number of coun- terexamples exist, there is a weak conjecture (Martin) that none exist un- der 10000 digits. In more detail, we are using the extra-strong Lucas test (Grantham 2000) using the Baillie parameter selection method (see OEIS A217719). Previous versions of this module used the strong Lucas test with Selfridge parame- ters, but the extra-strong version produces fewer pseudoprimes while run- ning 1.2 - 1.5x faster. It is slightly stronger than the test used in Pari <http://pari.math.u-bordeaux.fr/faq.html#primetest>. is_prime say "$n is prime!" if is_prime($n); Takes a positive number as input and returns back either 0 (composite), 2 (definitely prime), or 1 (probably prime). Composites will act exactly like "is_prob_prime", as will numbers less than "2^64". For numbers larger than "2^64", some additional tests are performed on probable primes to see if they can be proven by another means. This call walks the line between the performance of "is_prob_prime" and the certainty of "is_provable_prime". Those calls may be more appropriate in some cases. What this function does is give most of the performance of the former, while adding more certainty. For finer tuning of this tradeoff, you may instead use "is_prob_prime" followed by additional probable prime tests such as "miller_rabin_random" and/or "is_frobenius_underwood_pseudo- prime". As with "is_prob_prime", a BPSW test is first performed. This is determin- istic for all 64-bit numbers. Next, if the number is a Proth or LLR form, then a proof is constructed. If the result is still "probably prime" and the input is smaller than the Sorenson/Webster (2015) deterministic Miller-Rabin limit (approximately 82 bits) then the 11 or 12 Miller-Rabin tests are performed and the result is confirmed. For larger inputs that are still "probably prime" but under 200 bits, a quick BLS75 "n-1" primal- ity proof is attempted. This is tuned to give up if the result cannot be quickly determined, and results in success rates of ~80% at 80 bits, ~30% at 128 bits, and ~13% at 160 bits. Lastly, for results still "probably prime", an additional random-base Miller-Rabin test is performed. The result is that many numbers will return 2 (definitely prime), and the numbers that return 1 (probably prime) have gone through more tests than "is_prob_prime" while not taking too long. For cryptographic key generation, you may want even more testing for proba- ble primes (NIST recommends a few more additional M-R tests than we per- form). The function "miller_rabin_random" is made for this. Alternately, a different test such as "is_frobenius_underwood_pseudoprime" can be used. Even better, use "is_provable_prime" which should be reasonably fast for sizes under 2048 bits. Typically for key generation one wants random primes, and there are many functions for that. is_provable_prime say "$n is definitely prime!" if is_provable_prime($n) == 2; Takes a positive number as input and returns back either 0 (composite), 2 (definitely prime), or 1 (probably prime). A great deal of effort is taken to return either 0 or 2 for all numbers. The current method first uses BPSW to find composites and provide a deter- ministic answer for tiny numbers (under "2^64"). If no certificate is re- quired, LLR and Proth tests can be run, and small numbers (under approxi- mately "2^82") can be satisfied with a deterministic Miller-Rabin test. If the result is still not determined, a quick BLS75 "n-1" test is attempted, followed by ECPP. The time required for primes of different input sizes on a circa-2009 work- station averages about "3ms" for 30-digits, "5ms" for 40-digit, "20ms" for 60-digit, "50ms" for 80-digit, "100ms" for 100-digit, "2s" for 200-digit, and 400-digit inputs about a minute. Expect a lot of time variation for larger inputs. You can see progress indication if verbose is turned on (some at level 1, and a lot at level 2). A certificate can be obtained along with the result using the "is_prov- able_prime_with_cert" method. There is no appreciable extra performance cost for returning a certificate. is_provable_prime_with_cert Takes a positive number as input and returns back an array with two ele- ments. The result will be one of: (0, '') The input is composite. (1, '') The input is probably prime but we could not prove it. This is a failure in our ability to factor some necessary element in a reasonable time, not a significant proof failure (in other words, it remains a probable prime). (2, '...') The input is prime, and the certificate contains all the information necessary to verify this. The certificate is a text representation containing all the necessary in- formation to verify the primality of the input in a reasonable time. The result can be used with "verify_prime" in Math::Prime::Util for verifica- tion. Proof types used include: ECPP BLS3 BLS15 BLS5 Small is_pseudoprime Takes an integer "n" and zero or more positive bases as input. Returns 1 if the input is a probable prime to each base, 0 if not. This is the sim- ple Fermat primality test, checking that "b^(n-1) = 1 mod n". If no bases are given, base 2 is used. All bases must be 2 or greater. Removing primes, given base 2 this produces the sequence OEIS A001567 <http://oeis.org/A001567>. "n=2" is always returned as a pseudoprime. No other restrictions are used, so for instance, 4 is a Fermat pseudoprime base 5. is_euler_pseudoprime Takes an integer "n" and zero or more non-zero positive bases as input. Returns 1 if the input is an Euler probable prime to each base, 0 if not. This is the Euler test, sometimes called the Euler-Jacobi test, checking that "kronecker(b,n) = b^((n-1)/2) mod n" If no bases are given, base 2 is used. All bases must be 2 or greater. Removing primes, given base 2 this produces the sequence OEIS A047713 <http://oeis.org/A047713>. "n=2" is always returned as an Euler-Jacobi pseudoprime. All other even "n" return 0. All bases not coprime to "n" will return 0. is_strong_pseudoprime my $maybe_prime = is_strong_pseudoprime($n); my $probably_prime = is_strong_pseudoprime($n, 2, 3, 5, 7, 11, 13, 17); Takes an integer "n" and zero or more non-zero positive bases as input. Returns 1 if the input is a strong probable prime to each base, 0 if not. This is often called the Miller-Rabin test. If no bases are given, base 2 is used. All bases must be 2 or greater. 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 strong that the number is actually prime. Both the input number and the bases may be big integers. This is usually used in combination with other tests to make either stronger tests (e.g. the strong BPSW test) or deterministic results for numbers less than some verified limit (e.g. Jaeschke showed in 1993 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 multi- ple bases, there are some math packages that just use multiple MR tests for primality testing, though in the early 1990s almost all serious software switched to the BPSW test. "n=2" is always returned as strong pseudoprime. All other even "n" return 0 (most sources define this test for odd "n"). Returning 0 for all even composite "n" makes the function match most other implementations including Math::Primality's "is_strong_pseudoprime" function. If base modulo n <= 1 or base modulo n = n-1, then the result will be 1. The test is not well defined if "base = 0 mod n" (e.g. Crandall and Pomer- ance page 135), so we define it to return 1, which allows bases to be larger than n if desired while still returning meaningful results for all n. For example, is_strong_pseudoprime(367, 1101) would return 0 if this change was not made. This is important for using deterministic base sets with large bases. Note that this only matters if the base used is greater than or equal to "n". A 0 result will not be re- turned for a prime "n", regardless of the base. miller_rabin_random my $maybe_prime = miller_rabin_random($n, 10); # 10 random bases Takes an integer "n" as input and a positive number "k" indicating the num- ber of bases to use. Performs "k" Miller-Rabin tests using uniform random bases between 2 and "n-2". This is the correct way to perform "k" Miller-Rabin tests, rather than the common but broken method of using the first "k" primes. An optional third argument may be given, which is a seed to use. The seed should be a number either in decimal, binary with a leading "0b", hex with a leading "0x", or octal with a leading 0. It will be converted to a GMP integer, so may be large. Typically this is not necessary, but crypto- graphic applications may prefer the ability to use this, and it allows re- peatable test results. There is no check for duplicate bases. Input sizes below 65 bits make lit- tle sense for this function since is_prob_prime is deterministic at that size. For numbers of 65+ bits, the chance of duplicate bases is quite small. The exponentiation approximation for the birthday problem gives a probability of less than 2e-16 for 100 random bases to have a duplicate with a 65-bit input, and less than 2e-35 with a 128-bit input. is_lucas_pseudoprime is_strong_lucas_pseudoprime Given an integer "n", returns 1 if the input is a standard or strong Lucas probable prime. The Selfridge method of choosing D, P, and Q are used (some sources call this a Lucas-Selfridge test). This is one half of the BPSW primality test (the Miller-Rabin strong probable prime test with base 2 being the other half). The canonical BPSW test (page 1401 of Baillie and Wagstaff (1980)) uses the strong Lucas test with Selfridge parameters, but in practice a variety of Lucas tests with different parameters are used by tests calling themselves BPSW. The standard Lucas test implemented here corresponds to the Lucas test de- scribed in FIPS 186-4 section C.3.3, though uses a slightly more efficient calculation. Since the standard Lucas-Selfridge test is a subset of the strong Lucas-Selfridge test, I recommend using the strong test rather than the standard test for cryptographic purposes. It is often slightly faster, has over 4x fewer pseudoprimes, and is the method recommended by Baillie and Wagstaff in their 1980 paper. is_extra_strong_lucas_pseudoprime Given an integer "n", returns 1 if the input is an extra-strong Lucas prob- able prime. This is defined in Grantham (2000), and is a slightly more stringent test than the strong Lucas test, though because different parame- ters are used the pseudoprimes are not a subset. As expected by the extra conditions, the number of pseudoprimes is less than 2/3 that of the strong Lucas-Selfridge test. Runtime performance is 1.2 to 1.5x faster than the strong Lucas test. The parameters are selected using the Baillie-OEIS method: P = 3; Q = 1; while ( jacobi( P*P-4, n ) != -1 ) P += 1; is_almost_extra_strong_lucas_pseudoprime Given an integer "n", returns 1 if the input is an "almost" extra-strong Lucas probable prime. This is the classic extra-strong Lucas test but without calculating the U sequence. This makes it very fast, although as the input increases in size the time converges to the conventional ex- tra-strong implementation: at 30 digits this routine is about 15% faster, at 300 digits it is only 2% faster. With the current implementations, there is little reason to prefer this un- less trying to reproduce specific results. The extra-strong implementation has been optimized to use similar features, removing most of the perfor- mance advantage. An optional second argument (must be between 1 and 256) indicates the in- crement amount for P parameter selection. The default value of one yields the method described in "is_extra_strong_lucas_pseudoprime". A value of 2 yields the method used in Pari <http://pari.math.u-bor- deaux.fr/faq.html#primetest>. Because the "U = 0" condition is ignored, this produces about 5% more pseudoprimes than the extra-strong Lucas test. However this is still only 66% of the number produced by the strong Lucas-Selfridge test. No BPSW counterexamples have been found with any of the Lucas tests described. is_euler_plumb_pseudoprime Given an integer "n", returns 1 if the input passes Colin Plumb's Euler Criterion primality test. Pseudoprimes to this test are a subset of the base 2 Fermat and Euler tests, but a superset of the base 2 strong pseudo- prime (Miller-Rabin) test. The main reason for this test is that is very slightly more efficient than other probable prime tests. is_perrin_pseudoprime Given an integer "n", returns 1 if "n" divides P(n) where P(n) is the Per- rin number of "n". The Perrin sequence is defined by "P(n) = P(n-2) + P(n-3)" with "P(0) = 3, P(1) = 0, P(2) = 2". This is not a commonly used test, as it runs slower than most of the other probable prime tests and offers little benefit, especially over combined tests like "is_bpsw_prime", "is_frobenius_underwood_pseudoprime", and "is_frobenius_khashin_pseudoprime". An optional second argument "r" indicates whether to run additional tests. With "r=1", "P(-n) = -1 mod n" is also verified, creating the "minimal re- stricted" test. With "r=2", the full signature is also tested using the Adams and Shanks (1982) rules (without the quadratic form test). With "r=3", the full signature is tested using the Grantham (2000) test, which additionally does not allow pseudoprimes to be divisible by 2 or 23. The minimal restricted pseudoprime sequence is OEIS A018187 <http://oeis.org/A018187>. is_frobenius_pseudoprime Given an integer "n" and two optional integer parameters "a" and "b", re- turns 1 if "n" is positive and a Frobenius probable prime with respect to the polynomial "x^2 - ax + b". Without the parameters, "b = 2" and "a" is the least positive odd number such that "(a^2-4b|n) = -1". This selection has no pseudoprimes below "2^64" and none known. In any case, the discrim- inant "a^2-4b" must not be a perfect square. is_frobenius_underwood_pseudoprime Given an integer "n", returns 1 if "n" is positive and passes the efficient Frobenius test of Paul Underwood, and returns 0 otherwise. This selects a parameter "a" as the least non-negative integer such that "(a^2-4|n)=-1", then verifies that "(x+2)^(n+1) = 2a + 5 mod (x^2-ax+1,n)". This combines a Fermat and Lucas test at a computational cost of about 2.5x a strong pseudoprime test. This makes it similar to, but faster than, a standard Frobenius test. This test is deterministic (no randomness is used). There are no known pseudoprimes to this test. This test also has no overlap with the BPSW test, making it a very effective method for adding additional certainty. is_frobenius_khashin_pseudoprime Given an integer "n", returns 1 if "n" is positive and passes the Frobenius test of Sergey Khashin, and returns 0 otherwise. The test verifies "n" is not a perfect square, selects the parameter "c" as the smallest odd prime such that "(c|n)=-1", then verifies that "(1+D)^n = (1-D) mod n" where "D = sqrt(c) mod n". This test is deterministic (no randomness is used). There are no known pseudoprimes to this test. is_bpsw_prime Given an integer input, returns 0 (composite), 2 (definitely prime), or 1 (probably prime), using the BPSW primality test (extra-strong variant). This function does the extra-strong BPSW test and nothing more. That is, it will skip all pretests and any extra work that the "is_prob_prime" test may add. This saves some time if the input has no small factors, such as testing results that have been sieved. is_trial_prime Takes a positive number as input, and returns 1 if the input is not divisi- ble by any prime less than or equal to the square root of the input. This is simple trial division by primes. There is little practical use as we have much more efficient primality and compositeness tests, but it can be a useful comparison. is_aks_prime say "$n is definitely prime" if is_aks_prime($n); Takes a positive number as input, and returns 1 if the input passes the Agrawal-Kayal-Saxena (AKS) primality test. This is a deterministic uncon- ditional primality test which runs in polynomial time for general input. The particular method used is theorem 4.1 from Bernstein (2003). This is substantially faster than the original AKS publication, the later version with improvements by Lenstra (sometimes called the V6 paper), or the later improvements of Voloch and Bornemann. It is, by a large order, faster than any other known implementation as of early 2017. For theoretical analysis of the primality task, AKS is extremely important. In practice, it is essentially useless. Estimated run time for a 150 digit input is over 2 days, making the case that while the algorithmic complexity growth is polynomial, the constants are extremely high. It will take years for for numbers that ECPP or APR-CL can prove in seconds. With the "verbose" option set to 1, the chosen "r" and "s" values are printed before the test starts. With "verbose" set to 2 or higher, each of the "s" tests results in a "." output as the test runs, allowing progress to be monitored. Typically you should use "is_provable_prime" and let it decide the method. is_mersenne_prime say "2^607-1 (M607) is a Mersenne prime" if is_mersenne_prime(607); Takes a positive number "p" as input and returns 1 if "2^p-1" is prime. After some pre-testing, the Lucas-Lehmer test is performed. This is a de- terministic unconditional test that runs very fast compared to other pri- mality methods for numbers of comparable size, and vastly faster than any known general-form primality proof methods. is_llr_prime Takes a positive number "n" as input and returns one of: 0 (definitely com- posite), 2 (definitely prime), or -1 (test does not indicate anything). This implements the Lucas-Lehmer-Riesel test for fast deterministic primal- ity testing on numbers of the form "k * 2^n - 1". If the input is not of this form or if "k >= 2^n" then -1 will be returned as the test does not apply. If "k = 1" then this is a Mersenne number and the Lucas-Lehmer test is used. Otherwise, the LLR test is performed. While not as fast as the Lucas-Lehmer test for Mersenne numbers, it is almost as fast as a single strong pseudoprime test (i.e. Miller-Rabin test) while giving a certain answer. is_proth_prime Takes a positive number "n" as input and returns one of: 0 (definitely com- posite), 2 (definitely prime), or -1 (test does not indicate anything). This applies Proth's theorem for fast Las Vegas primality testing on num- bers of the form "k * 2^n + 1". If the input is not of this form or if "k >= 2^n" then -1 will be returned as the test does not apply. Otherwise, a search is performed to find a quadratic nonresidue modulo "n". If none can be found after a brief search, -1 is returned as no conclusion can be reached. Otherwise, Proth's theorem is checked which conclusively indi- cates primality. While not as fast as the Lucas-Lehmer test for Mersenne numbers, it is almost as fast as a single strong pseudoprime test (i.e. Miller-Rabin test) while giving a certain answer. is_gaussian_prime Given two integers "a" and "b" as input, returns 0 (definitely composite), 1 (probably prime), or 2 (definitely prime) to indicate whether the complex number "a+bi" is a Gaussian prime. The "is_prime" function is used inter- nally to make the determination. is_miller_prime say "$n is definitely prime" if is_miller_prime($n); say "$n is definitely prime assuming the GRH" if is_miller_prime($n, 1); Takes a positive number as input, and returns 1 if the input passes the de- terministic Miller test. An optional second argument indicates whether the Generalized Riemann Hypothesis should be assumed, and defaults to 0. Set- ting the verbose flag to 2 or higher will show how many bases are used. The unconditional test is exponential time, while the conditional test (as- suming the GRH) is polynomial time. This is a very slow method in practice, and generally should not be used. The asymptotic complexity of the GRH version is good in theory, matching ECPP, but in practice it is much slower. The number of bases used by the unconditional test grows quite rapidly, impractically many past about 160 bits, and overflows a 64-bit integer at 456 bits -- sizes that are trivial for the unconditional APR-CL and ECPP tests. is_nminus1_prime say "$n is definitely prime" if is_nminus1_prime($n); Takes a positive number as input, and returns 1 if the input passes either theorem 5 or theorem 7 of the Brillhart-Lehmer-Selfridge primality test. This is a deterministic unconditional primality test which requires factor- ing "n-1" to a linear factor less than the cube root of the input. For small inputs (under 40 digits) this is typically very easy, and some num- bers will naturally lead to this being very fast. As the input grows, this method slows down rapidly. This method is most appropriate for numbers of the form "k+1" where "k" can be easily factored. This is a test specifically for this proof method. The return values are: 1 We constructed a primality proof, hence C<n> is definitely prime. 0 We were unable to construct a proof, hence no conclusive result. An answer of 0 might be because "n" is composite, but it could also be a failure to factor or find the right parameters. Typically you should use "is_provable_prime" and let it decide the method and returns definitely prime, probably prime, or definitely composite. is_nplus1_prime Takes a positive number as input, and returns 1 if the input passes either theorem 17 or theorem 19 of the Brillhart-Lehmer-Selfridge primality test. This is a deterministic unconditional primality test which requires factor- ing "n+1" to a linear factor less than the cube root of the input. For small inputs (under 40 digits) this is typically very easy, and some num- bers will naturally lead to this being very fast. As the input grows, this method slows down rapidly. Disregarding factoring, this is slightly slower than the "n-1" methods. It is most appropriate for numbers of the form "k-1" where "k" can be easily factored. This is a test specifically for this proof method. The return values are: 1 We constructed a primality proof, hence C<n> is definitely prime. 0 We were unable to construct a proof, hence no conclusive result. An answer of 0 might be because "n" is composite, but it could also be a failure to factor or find the right parameters. Typically you should use "is_provable_prime" and let it decide the method and returns definitely prime, probably prime, or definitely composite. is_bls75_prime Takes a positive number as input, and returns 1 if the input passes one of the tests from the Brillhart-Lehmer-Selfridge (1975) paper. These use par- tial factoring of "n-1" and "n+1". Currently the implementation will use one of: N-1 Corollary 1, Theorem 5, Theorem 7 N+1 Corollary 8, Theorem 17, Theorem 19 Comb Theorem 20 This is appropriate for cases where either "n-1" or "n+1" can be easily factored, or when both of them have many small factors. This is a test specifically for this proof method. The return values are: 1 We constructed a primality proof, hence C<n> is definitely prime. 0 We were unable to construct a proof, hence no conclusive result. An answer of 0 might be because "n" is composite, but it could also be a failure to factor or find the right parameters. Typically you should use "is_provable_prime" and let it decide the method and returns definitely prime, probably prime, or definitely composite. is_ecpp_prime say "$n is definitely prime" if is_ecpp_prime($n); Takes a positive number as input, and returns 1 if the input passes the ECPP primality test. This is the Atkin-Morain Elliptic Curve Primality Proving algorithm. It is the fastest primality proving method in Math::Prime::Util. This implementation uses a "factor all strategy" (FAS) with backtracking. A limited set of about 500 precalculated discriminants are used, which works well for inputs up to 300 digits, and for many inputs up to one thou- sand digits. Having a larger set will help with large numbers (a set of 2650 is available on github in the "xt/" directory). A future implementa- tion may include code to generate class polynomials as needed. This is a test specifically for this proof method. The return values are: 1 We constructed a primality proof, hence C<n> is definitely prime. 0 We were unable to construct a proof, hence no conclusive result. An answer of 0 might be because "n" is composite, but it could also be a failure to factor or find the right parameters. Typically you should use "is_provable_prime" and let it decide the method and returns definitely prime, probably prime, or definitely composite. primes my $aref1 = primes( 1_000_000 ); my $aref2 = primes( 2 ** 448, 2 ** 448 + 10000 ); say join ",", @{primes( 2**2048, 2**2048 + 10000 )}; Returns all the primes between the lower and upper limits (inclusive), with a lower limit of 2 if none is given. An array reference is returned, matching the signature of the function of the same name in Math::Prime::Util. Values above 64-bit are extra-strong BPSW probable primes. prime_count Returns the number of primes up to "n" (single argument) or between "lo" and "hi" (two arguments). The values are inclusive. The method is simple sieving followed by primality testing. This is appro- priate for small ranges and is useful for very large arguments. The Math::Prime::Util module has much more sophisticated methods for 64-bit ar- guments. prime_count_lower prime_count_upper Returns lower or upper bounds for the prime count of the input "n". Bounds use Dusart 2010, BA1/4the 2014, BA1/4the 2015, and Axler 2017. sieve_primes my @primes = sieve_primes(2**100, 2**100 + 10000); my @candidates = sieve_primes(2**1000, 2**1000 + 10000, 40000); Given two arguments "low" and "high", this returns the primes in the inter- val (inclusive) as a list. It operates similar to primes, though must al- ways have an lower and upper bound and returns a list. With three arguments "low", "high", and "limit", this does a partial sieve over the inclusive range and returns the list that pass the sieve. If "limit" is less than 2 then it is identical to the two-argument version, in that a primality test will be performed after sieving. Otherwise, sieving is performed up to and including "limit". The two-argument version is typically only used internally and adds little functionality. The three-argument version is quite useful for applications that want to apply their own primality or other tests, and wish to have a list of values in the range with no small factors. This is quite common for applications involving prime gaps. Also see "sieve_range". sieve_range my @candidates = sieve_range(2**1000, 10000, 40000); Given a start value "n", and native unsigned integers "width" and "depth", a sieve of maximum depth "depth" is done for the "width" consecutive num- bers beginning with "n". An array of offsets from the start is returned. The returned list contains those offsets in the range "n" to "n+width-1" where "n + offset" has no prime factors less than or equal to "depth". This function is very similar to the three argument form of "sieve_primes". The differences are using "(n,width)" instead of "(low,high)", and most im- portantly returning small offsets from the start value rather than the val- ues themselves. This can substantially reduce overhead for multi-thousand digit numbers. sieve_twin_primes my @primes = sieve_twin_primes(2**1000, 2**1000 + 500000); Given two arguments "low" and "high", this returns each lower twin prime in the interval (inclusive). The result is a list, not a reference. This does a partial sieve of the range, removes any non-twin candidates, then checks that each pair are both BPSW probable primes. This is substan- tially more efficient than sieving for all primes followed by removing those that are not twin primes. sieve_prime_cluster # Find some prime septuplets my @s = sieve_prime_cluster(2**100, 2**100+1e12, 2,6,8,12,18,20); Efficiently finds prime clusters between the first two arguments "low" and "high" (inclusive). The remaining arguments describe the cluster. The cluster values must be even, less than 31 bits, and strictly increasing. Given a cluster set "C", the returned values are all primes in the range where "p+c" is prime for all "c" in the cluster set "C". The cluster is described as offsets from 0, with the implicit prime at 0. Hence an empty list is asking for all primes (the cluster "p+0"). A list with the single value 2 will find all twin primes (the cluster where "p+0" and "p+2" are prime). The list "2,6,8" will find prime quadruplets. Note that there is no requirement that the list denote a constellation (a clus- ter with minimal distance) -- the list "42,92,606" is just fine. For long clusters, e.g. OEIS series A213601 <http://oeis.org/A213601> prime 12-tuplets, this will be immensely more efficient than filtering out the cluster from a list of primes. For that example, a range of "10^13" takes less than a second to search -- thousands of times faster than filtering results from primes or twin primes. Shorter clusters are not quite this efficient, and the overhead for returning large arrays should not be ig- nored. next_prime $n = next_prime($n); Returns the prime following the input number (the smallest prime number that is greater than the input number). The function "is_prob_prime" is used to determine when a prime is found, hence the result is a probable prime (using BPSW). For large inputs this function is quite a bit faster than GMP's "mpz_nextprime" or Pari's "nextprime". prev_prime $n = prev_prime($n); Returns the prime preceding the input number (the largest prime number that is less than the input number). 0 is returned if the input is 2 or lower. The function "is_prob_prime" is used to determine when a prime is found, hence the result is a probable prime (using BPSW). surround_primes ($dprev, $dnext) = surround_primes($n); Returns the distances to the previous and next primes of the input "n". This is slightly more efficient than calling both "prev_prime" and "next_prime", and returning the distances as native integers is more effi- cient with large inputs. If an optional second argument "d" is given, and the input "n" is larger than "2^64", then if a SPSP-2 is found in the range "n-d" to "n+d" (inclu- sive) then it will be returned with the other argument set to 0. Other- wise, the first SPSP-2 values found are returned. This feature is espe- cially useful for prime gap searches as well as finding the nearest prime to a value. Note that with a non-zero second argument, the values returned have not un- dergone a full BPSW test; just sieving and a SPSP-2 test. next_twin_prime Returns the start of the next twin prime after the input "n". The returned value will always be greater than the input. For a return value of "t", both "t" and "t+2" will be a probable prime (using BPSW). random_nbit_prime say "random 512-bit prime: ", random_nbit_prime(512); Returns a randomly selected prime of exactly "n" bits. "undef" is returned if "n" is less than 2. The returned prime has passed the "is_prob_prime" (extra strong BPSW) test. random_safe_prime Returns a randomly selected safe prime of exactly "n" bits. "n" must be at least 3. The returned value "p" will be of the form "2q+1" where "q" has passed the extra strong BPSW test. "p" is guaranteed to be a prime if "q" is prime. Setting verbose level to 3 or higher will produce progress output not un- like "openssl". A "." for each candidate that passes pretests, a "+" for those where one is likely prime, and "*" when both are likely prime with only confirmation tests remaining. This generates safe primes about 4-10x faster than "openssl"'s "dhparam". random_strong_prime say "random 512-bit strong prime: ", random_strong_prime(512); Returns a randomly selected strong prime of exactly "n" bits. "n" must be at least 256. The returned prime has passed the "is_prob_prime" (extra strong BPSW) test. Given the returned prime "p", "p+1", "q=p-1", and "q-1" will all have a large factor. This makes using factoring methods such as p-1 and p+1 much harder. Gordon's algorithm is used. The value of using strong primes is questionable over proper random primes when the number of bits is at least 1024. random_ndigit_prime say "random 200-digit prime: ", random_ndigit_prime(200); Returns a randomly selected prime of exactly "n" digits. "undef" is re- turned if "n" is less than 1. The returned prime has passed the "is_prob_prime" (extra strong BPSW) test. random_prime say random_prime(1000, 2000); # prime between 1000 and 2000 inclusive Returns a random prime in the interval "[a,b]" or "undef" if no prime is in the range. The returned prime has passed the "is_prob_prime" (extra strong BPSW) test. The random prime functions use the internal CSPRNG for randomness. This is currently ISAAC-32 but will likely change to ChaCha20 in a later release. This corresponds to Mathematica's "RandomPrime[{min,max}]" function. This is a superset of Pari's randomprime(n) function, where our interval API is more convenient for cryptographic functions. random_maurer_prime say "random 512-bit proven prime: ", random_maurer_prime(512); Returns an n-bit proven prime using Ueli Maurer's FastPrime algorithm (1995). This results in uniform random selection of a proven prime, though not every n-bit prime can be generated with this algorithm. "undef" is returned if "n" is less than 2. As a safety check, internally the extra strong BPSW test is additionally run on each intermediate and the final result. random_shawe_taylor_prime say "random 512-bit proven prime: ", random_shawe_taylor_prime(512); Returns an n-bit proven prime using the Shawe-Taylor algorithm (1986) from section C.6 of FIPS 186-4, although using our CSPRNG rather than SHA-256. This is a slightly simpler and older (1986) method than Maurer's algorithm. It is a bit faster than Maurer's method but has a smaller subset of re- turned primes. "undef" is returned if "n" is less than 2. As a safety check, internally the extra strong BPSW test is additionally run on each intermediate and the final result. random_maurer_prime_with_cert Like "random_maurer_prime" but also returns a string certificate. random_shawe_taylor_prime_with_cert Like "random_shawe_taylor_prime" but also returns a string certificate. lucasu say "Fibonacci($_) = ", lucasu(1,-1,$_) for 0..100; Given integers "P", "Q", and the non-negative integer "k", computes "U_k" for the Lucas sequence defined by "P","Q". These include the Fibonacci numbers ("1,-1"), the Pell numbers ("2,-1"), the Jacobsthal numbers ("1,-2"), the Mersenne numbers ("3,2"), and more. This corresponds to OpenPFGW's "lucasU" function and gmpy2's "lucasu" func- tion. lucasv say "Lucas($_) = ", lucasv(1,-1,$_) for 0..100; Given integers "P", "Q", and the non-negative integer "k", computes "V_k" for the Lucas sequence defined by "P","Q". These include the Lucas numbers ("1,-1"). This corresponds to OpenPFGW's "lucasV" function and gmpy2's "lucasv" func- tion. lucasuv Given integers "P", "Q", and the non-negative integer "k", return ("U_k" modulo "n", "V_k" modulo "n") for the Lucas sequence defined by "P","Q". lucasumod Given integers "P", "Q", the non-negative integer "k", and the integer "n", computes "U_k" modulo "|n|" for the Lucas sequence defined by "P","Q". lucasvmod Given integers "P", "Q", the non-negative integer "k", and the integer "n", computes "V_k" modulo "|n|" for the Lucas sequence defined by "P","Q". This can be particularly efficient, especially if "Q=1". lucasuvmod Given integers "P", "Q", the non-negative integer "k", and the integer "n", returns ("U_k" modulo "|n|", "V_k" modulo "|n|") for the Lucas sequence de- fined by "P","Q". lucas_sequence my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k) Computes "U_k", "V_k", and "Q_k" for the Lucas sequence defined by "P","Q", modulo "n". The modular Lucas sequence is used in a number of primality tests and proofs. "k" must be non-negative, and "n" must be greater than zero. "P" and "Q" are restricted to native signed integers. The newer function "lucasuvmod" accepts bigint values. primorial $p = primorial($n); Given an unsigned integer argument, returns the product of the prime num- bers which are less than or equal to "n". This definition of "n#" follows OEIS series A034386 <http://oeis.org/A034386> and Wikipedia: Primorial def- inition for natural numbers <http://en.wikipedia.org/wiki/Primorial#Defini- tion_for_natural_numbers>. pn_primorial $p = pn_primorial($n) Given an unsigned integer argument, returns the product of the first "n" prime numbers. This definition of "p_n#" follows OEIS series A002110 <http://oeis.org/A002110> and Wikipedia: Primorial definition for prime numbers <http://en.wikipedia.org/wiki/Primorial#Definition_for_prime_num- bers>. The two are related with the relationships: pn_primorial($n) == primorial( nth_prime($n) ) primorial($n) == pn_primorial( prime_count($n) ) factorial Given positive integer argument "n", returns the factorial of "n", defined as the product of the integers 1 to "n" with the special case of "factor- ial(0) = 1". This corresponds to Pari's factorial(n) and Mathematica's "Factorial[n]" functions. multifactorial Given two positive integer arguments "n" and "m", returns "n!^(m)", the multifactorial. "m=1" is the standard "factorial" while "m=2" is the dou- ble factorial. While the factorial is the product of all integers "n" and below, the multifactorial skips those without the same parity as "n mod m". Hence multifactorial(n,2) = n * (n-2) * (n-4) * ... The multifactorials are the OEIS series (m=1) OEIS series A000142 <http://oeis.org/A000142>, (m=2) OEIS series A006882 <http://oeis.org/A006882>, (m=3) OEIS series A007661 <http://oeis.org/A007661>, (m=4) OEIS series A007662 <http://oeis.org/A007662>, ... Also see Peter Luschny's excellent writeup <http://oeis.org/wiki/User:Pe- ter_Luschny/Multifactorials>. factorialmod Given positive integer "n" and integer "m", returns "n! mod |m|". This is much faster than computing the large factorial(n) followed by a mod opera- tion. The implementation is a modular product of prime powers, including grouping of like powers to reduce the number of exponentiations. While this is quite efficient, the fast multipoint evaluation method (see Johannson (2012), used in FLINT) is significantly better. factorial_sum Given positive integer argument "n", returns the factorial sum of "n". This is defined as the sum of factorial(k) for "0 .. n-1". These are equivalent, though this function is faster: factorial_sum($n) == vecsum(map{ factorial($_) } 0..$n-1) This is sometimes called the left factorial, confusingly also used for "subfactorial". This is OEIS series A003422 <http://oeis.org/A003422>. subfactorial Given positive integer argument "n", returns the subfactorial of "n". This is also called the derangement number, and occasionally the left factorial. This is OEIS series A000166 <http://oeis.org/A000166>. This corresponds to Mathematica's "Subfactorial[n]" function. falling_factorial Given two integers "x" and "n", with "n" non-negative, returns the falling factorial of "n". falling_factorial(x,n) = x * (x-1) * (x-2) * ... * (x-(n-1)) This corresponds to Mathematica's "FactorialPower[x,n]" function. rising_factorial Given two integers "x" and "n", with "n" non-negative, returns the rising factorial of "n". rising_factorial(x,n) = x * (x+1) * (x+2) * ... * (x+(n-1)) This corresponds to Mathematica's "Pochhammer[x,n]" function. gcd Given a list of integers, returns the greatest common divisor. This is of- ten used to test for coprimality <https://oeis.org/wiki/Coprimality>. An empty list will return 0. lcm Given a list of integers, returns the least common multiple. An empty list will return 1. gcdext Given two integers "x" and "y", returns "u,v,d" such that "d = gcd(x,y)" and "u*x + v*y = d". This uses the extended Euclidian algorithm to compute the values satisfying BA(C)zout's Identity. This corresponds to Pari's "gcdext" function, which was renamed from "be- zout" in Pari 2.6. The results will hence match "bezout" in Math::Pari. chinese say chinese( [14,643], [254,419], [87,733] ); # 87041638 Solves a system of simultaneous congruences using the Chinese Remainder Theorem (with extension to non-coprime moduli). A list of "[a,n]" pairs are taken as input, each representing an equation "x a! a mod |n|". If no solution exists, "undef" is returned. If a solution is returned, the modu- lus is equal to the lcm of all the given moduli (see "lcm"). In the stan- dard case where all values of "n" are coprime, this is just the product. The "n" values must be non-zero integers, while the "a" values are inte- gers. chinese2 Functions like "chinese" but returns two items: the remainder and the modu- lus. If a solution exists, the second value (the final modulus) is equal to the lcm of the absolute values of all the given moduli. If no solution exists, both return values will be "undef". vecsum Returns the sum of all arguments, each of which must be an integer. vecprod Returns the product of all arguments, each of which must be an integer. kronecker Returns the Kronecker symbol "(a|n)" for two integers. The possible return values with their meanings for odd prime "n" are: 0 a = 0 mod n 1 a is a quadratic residue modulo n (a = x^2 mod n for some x) -1 a is a quadratic non-residue modulo n The Kronecker symbol is an extension of the Jacobi symbol to all integer values of "n" from the latter's domain of positive odd values of "n". The Jacobi symbol is itself an extension of the Legendre symbol, which is only defined for odd prime values of "n". This corresponds to Pari's "kro- necker(a,n)" function and Mathematica's "KroneckerSymbol[n,m]" function. If "n" is not an odd prime, then the result does not necessarily indicate whether "a" is a quadratic residue mod "n". Using the function "is_qr" will return correct results for any "n", though could be slower. binomial Given integer arguments "n" and "k", returns the binomial coefficient "n*(n-1)*...*(n-k+1)/k!", also known as the choose function. Negative ar- guments use the Kronenburg extensions <http://arxiv.org/abs/1105.3689/>. This corresponds to Mathematica's "Binomial[n,k]" function, Pari's "bino- mial(n,k)" function, and GMP's "mpz_bin_ui" function. The second argument "k" must fit in an unsigned long. For negative arguments, this matches Mathematica. Pari does not implement the "n < 0, k <= n" extension and instead returns 0 for this case. GMP's API does not allow negative "k" but otherwise matches. "Math::BigInt" ver- sion 1.999816 and later supports negative arguments with similar semantics. Prior to this, "n < 0, k > 0" was undefined. binomialmod Given integer arguments "n", "k", and "m", returns "binomial(n,k) mod m" efficiently. All requirements for binomial must be met including the size of k and the treatment of negative "n" and "k". The current implementation is not optimized, but can greatly reduce the overhead associated with the very large intermediate. addreal =head2 subreal =head2 mulreal =head2 divreal Returns the corresponding basic math operation applied to the two inputs. An optional third argument indicates the number of significant digits (de- fault 40) with the result rounded. logreal Returns the natural logarithm of the input "n". An optional second argu- ment indicates the number of significant digits (default 40) with the re- sult rounded. For logreal(2) we use Formula 25 from Gourdon and Sebah (2010). For other values we use AGM (Sasaki and Kanada theta method). Performance is 100-1000x faster than Math::BigFloat's GMP backend. It is 10x slower than Pari/GP 2.10 and MPFR. Negative inputs are returned as "-log(-n)", which matches bignum. Pari/GP and Mathematica return "log(-n) + Pi*i". expreal Returns "e^n" for the input "n". An optional second argument indicates the number of significant digits (default 40) with the result rounded. The implementation computes sinh(n), then "e^x" from that. powreal Returns "x^n" for the inputs "x" and "n". An optional third argument indi- cates the number of significant digits (default 40) with the result rounded. Like "logreal" and "expreal", this is a basic math function that is not available from the GMP library but implemented in MPFR. Since the latter is not always available, this can be useful to have. A negative "x" with non-integer "n" will return undef as the result is com- plex. "x=0" with negative "n" will also return undef as this is dividing by zero. rootreal Returns the "n"-th root of "x" for the inputs "n" and "x". An optional third argument indicates the number of significant digits (default 40) with the result rounded. This is just powreal(x^(1/n), but allows "n" to be easily specified in full precision. Asking for the 0-th root will return undef, as will negative "x" with roots other than 1 or -1. agmreal Returns the Arithmetic-Geometric mean (AGM) of "a" and "b". An optional third argument indicates the number of significant digits (default 40) with the result rounded. Examples of use include elementary constants (e.g. "Pi" and "e"), logs, ex- ponentials, trigonometric functions, elliptic integrals, computing pendulum periods, and more. This corresponds to Pari's "agm(x,y)" function, limited to positive reals (Pari also handles negative, complex, p-adic, and power series arguments). bernfrac Returns the Bernoulli number "B_n" for an integer argument "n", as a ratio- nal number. Two values are returned, the numerator and denominator. B_1 = 1/2. This corresponds to Pari's bernfrac(n) and Mathematica's "BernoulliB" functions. bernreal Returns the Bernoulli number "B_n" for an integer argument "n", as a string floating point. An optional second argument indicates the number of sig- nificant digits to be used, with the result rounded. The default is 40 digits. This corresponds to Pari's "bernreal" function. bernvec say join ", ", map { join "/" ,@$_; } bernvec(8); # 1/1, 1/6, -1/30, 1/42, -1/30, 5/66, -691/2730, 7/6, -3617/510 Returns an array of Bernoulli numbers "B_0, B_2, ..., B_2n" for an integer argument "n". Each array element is an 2-element array reference contain- ing the numerator and denominator. This will be faster than calling bernfrac for each value. The computed values are stored in a cache (not incremental), so later calls to to "bern- frac", "bernreal", "bernvec", and "powersum" with an "n" in the range will be substantially faster. When called in void context this just calculates and caches the values. The cache is not thread aware. This corresponds to Pari's bernvec(n) function, albeit noticeably slower. harmfrac Returns the Harmonic number "H_n" for an integer argument "n", as a ratio- nal number. Two values are returned, the numerator and denominator. num- bers are the sum of reciprocals of the first "n" natural numbers: "1 + 1/2 + 1/3 + ... + 1/n". This corresponds to Mathematica's "HarmonicNumber" function. harmreal Returns the Harmonic number "H_n" for an integer argument "n", as a string floating point. An optional second argument indicates the number of digits to be preserved past the decimal place, with a default of 40. stirling say "s(14,2) = ", stirling(14, 2); say "S(14,2) = ", stirling(14, 2, 2); Returns the Stirling numbers of either the first kind (default), the second kind, or the third kind (the unsigned Lah numbers), with the kind selected as an optional third argument. It takes two non-negative integer arguments "n" and "k" plus the optional "type". This corresponds to Pari's "stir- ling(n,k,{type})" function and Mathematica's "StirlingS1" / "StirlingS2" functions. Stirling numbers of the first kind are "-1^(n-k)" times the number of per- mutations of "n" symbols with exactly "k" cycles. Stirling numbers of the second kind are the number of ways to partition a set of "n" elements into "k" non-empty subsets. The Lah numbers are the number of ways to split a set of "n" elements into "k" non-empty lists. zeta Given a positive integer or float "n", returns the real Riemann Zeta value as a string floating point. An optional second argument indicates the num- ber of digits past the decimal point (default 40). The implementation is algorithm 2 of Borwein (1991). Performance with in- teger inputs is good, but floating point arguments with high precision will be slower than methods using MPFR. Math::Prime::Util will try to use Math::MPFR if possible. li Given a positive integer or float "n", returns the real Logarithmic Inte- gral as a string floating point. An optional second argument indicates the number of significant digits (default 40) with the result rounded. The implementation uses Ramanujan's series. This corresponds to Mathemat- ica's "Li" function. ei Given a positive integer or float "n", returns the real Exponential Inte- gral as a string floating point. An optional second argument indicates the number of significant digits (default 40) with the result rounded. The implementation is simply li(exp(x)). riemannr Given a positive integer or float "n", returns the real Riemann R function as a string floating point. An optional second argument indicates the num- ber of significant digits (default 40) with the result rounded. The implementation is the standard Gram series. This corresponds to Mathe- matica's "RiemannR" function. lambertw Given a float "x", returns the principal branch of the Lambert W function. This solves for "W" in the equation "x = W*exp(W)". The input must not be less than "-1/e". This corresponds to Pari's "lambertw" function and Math- ematica's "ProductLog" / "LambertW" function. znorder $order = znorder(17, "100000000000000000000000065"); Given two positive integers "a" and "n", returns the multiplicative order of "a" modulo "n". This is the smallest positive integer "k" such that "a^k a! 1 mod n". Returns 1 if "a = 1". Returns undef if "a = 0" or if "a" and "n" are not coprime, since no value will result in 1 mod n. This corresponds to Pari's "znorder(Mod(a,n))" function and Mathematica's "Mul- tiplicativeOrder[a,n]" function. znprimroot Given a positive integer "n", returns the smallest primitive root of "(Z/nZ)^*", or "undef" if no root exists. A root exists when "totient($n) == carmichael_lambda($n)", which is true only if "n one of {2, 4, p^k, 2p^k}" for odd prime p. OEIS A033948 <http://oeis.org/A033948> is a sequence of integers where the primitive root exists, while OEIS A046145 <http://oeis.org/A046145> is a list of the smallest primitive roots, which is what this function produces. This is similar to Pari/GP's "znprimroot", but our function always returns the smallest root. Pari/GP does not if "n = 2p^k". is_primitive_root Given two non-negative numbers "a" and "n", returns 1 if "a" is a primitive root modulo "|n|", and 0 if not. If "a" is a primitive root, then to- tient(n) is the smallest "e" for which "a^e = 1 mod n". is_semiprime Given a positive integer "n", returns 1 if "n" is a semiprime, 0 otherwise. A semiprime is the product of exactly two primes. The boolean result is the same as "scalar(factor(n)) == 2", but this func- tion performs shortcuts that can greatly speed up the operation. is_almost_prime say is_almost_prime(6,2169229601); # True if n has exactly 6 factors Given positive integers "k" and "n", returns 1 if "n" has exactly "k" prime factors, and 0 otherwise. With "k=1", this is a standard primality test. With "k=2", this is the same as "is_semiprime". is_smooth Given two non-negative integer inputs "n" and "k", returns 1 if "|n|" is "k"-smooth, and 0 otherwise. This uses the OEIS definition: Returns true if no prime factors of "n" are larger than "k". is_rough Given two non-negative integer inputs "n" and "k", returns 1 if "|n|" is "k"-rough, and 0 otherwise. This uses the OEIS definition: Returns true if no prime factors of "n" are smaller than "k". is_powerful Given an integer input "n" and an optional non-negative integer input "k", returns 1 if "n" is "k"-powerful, and 0 otherwise. If "k" is omitted, "k=2" is used. A k-powerful number is a positive integer where all prime factors appear at least "k" times. All positive integers are therefore 0- and 1-powerful. "n=1" is powerful for all "k". 0 is returned for all negative or zero val- ues of "n". is_practical Given a non-negative integer "n", returns 1 if "n" is a practical number, and 0 otherwise. A practical number is a positive integer "n" such that all smaller positive integers can be represented as sums of distinct divi- sors of "n". This is OEIS series A005153 <http://oeis.org/A005153>. is_carmichael Given a positive integer "n", returns 1 if "n" is a Carmichael number, 0 otherwise. These are composites that satisfy "b^(n-1) a! 1 mod n" for all "1 < b < n" relatively prime to "n". Alternately Korselt's theorem says these are composites such that "n" is square-free and "p-1" divides "n-1" for all prime divisors "p" of "n". Inputs greater than 50 digits use a probabilistic test to avoid fully fac- toring the input. is_fundamental Given a positive integer "n", returns 1 if "n" is a fundamental discrimi- nant, 0 otherwise. is_totient Given an integer "n", returns 1 if there exists an integer "x" where "to- tient(x) == n". is_polygonal Given integers "x" and "s", return 1 if x is an s-gonal number, 0 other- wise. "s" must be greater than 2. polygonal_nth Given integers "x" and "s", return N if "x" is the "N-th" s-gonal number, 0 otherwise. is_square_free Given integer "n", returns 1 if "|n|" has no repeated factor. is_powerfree Given an integer "n" and an optional non-negative integer "k", returns 1 if "|n|" has no divisor "d^k", and returns 0 otherwise. This determines if "|n|" has any "k"-th (or higher) powers in the prime factorization. "k" defaults to 2. next_powerfree Given an integer "n" and an optional non-negative integer "k", returns the next k-powerfree integer. The result will be the smallest k-powerfree in- teger greater than "n". "k" defaults to 2. prev_powerfree Given an integer "n" and an optional non-negative integer "k", returns the previous k-powerfree integer. The result will be the largest k-powerfree integer smaller than "n", unless "n <= 1", in which case undef is returned. "k" defaults to 2. powerfree_count Given an integer "n" and an optional non-negative integer "k", returns the number of k-powerfree positive integers less than or equal to "n". "k" de- faults to 2. nth_powerfree Given a non-negative integer "n" and an optional non-negative integer "k", returns the "n"-th k-powerfree number. If "k" is omitted, "k=2" is used. Returns undef if "k" is less than 2 or "n=0". Returns 1 for "n=1". powerful_count my $npowerful3 = powerful_count(2**32, 3); Given two non-negative integer inputs "n" and "k", returns the total number of "k"-powerful numbers from 1 to "n" inclusive. If "k" is omitted or zero, "k=2" is used. is_perfect_power Given an integer "n", returns 1 if "n" is a perfect power, and 0 otherwise. That is, if "n = c^k" for some integers "c" and "k" with "k" greater than 1. The results match the mpz_perfect_power_p(n) function of GMP 4.0+. Follow- ing GMP, SAGE, and FLINT, we treat -1, 0, and 1 as perfect powers. For positive integers, this is OEIS series A001597 <http://oeis.org/A001597>. next_perfect_power Given an integer "n", returns the smallest perfect power greater than "n". Similar in API to "next_prime", but returns the next perfect power with ex- ponent greater than 1. Starting from 0 this gives the sequence "1,4,8,9,16,25,...". Negative inputs are supported, with the result being the nearest value greater than "n" where "is_perfect_power" returns true. prev_perfect_power Given an integer "n", returns the greatest perfect power less than "n". Similar in API to "prev_prime", but returns the previous perfect power with exponent greater than 1. Negative inputs are supported, with the result being the nearest value less than "n" where "is_perfect_power" returns true. perfect_power_count Returns the number of perfect powers up to "n" (single argument) or between "lo" and "hi" (two arguments). The values are inclusive, and must be non-negative. By convention, numbers less than 1 are not counted. This is OEIS series A069623 <http://oeis.org/A069623>. nth_perfect_power Given a non-negative integer "n", returns the "n"-th perfect power. Since the perfect power count can be calculated quite quickly, we use suc- cessive approximations with corrections provided by the count to calculate the "n"-th perfect power quite rapidly. Similar to "perfect_power_count", the convention is to exclude all integers less than 1. Hence "n=0" returns undef and "n=1" returns 1. nth_perfect_power_approx Given a non-negative integer "n", quickly returns a good estimate of the "n"-th perfect power. prime_power_count Returns the number of positive prime powers up to "n" (single argument) or between "lo" and "hi" (two arguments). The values are inclusive. This is OEIS series A025528 <http://oeis.org/A025528>. powersum faulhaber_sum say powersum(111,18); # 41588295196841906092077874022002239896 Given "n" and "p", returns the sum of the "p"-th powers of the first "n" positive integers. That is, the sum of "k^p" for "k = 1 .. n". Bernoulli's formula using the Bernoulli numbers is used. This uses the Bernoulli number cache. Hence if calling multiple times, the largest "p" value should be used first, or "bernvec(p >> 1)" should be called in void context first. Additionally, because the cache is used, this is not thread safe in gen- eral. "bernvec" should be called once on a single thread, then future calls to "powersum" and the other functions can be safely done on multiple threads if and only if within that range. "powersum" and "faulhaber_sum" are aliases for the same function. sigma say "Sum of divisors of $n:", sigma( $n ); say "sigma_2($n) = ", sigma($n, 2); say "Number of divisors: sigma_0($n) = ", sigma($n, 0); This function takes a positive integer as input and returns the sum of its divisors, including 1 and itself. An optional second argument "k" may be given, which will result in the sum of the "k-th" powers of the divisors to be returned. This is known as the sigma function (see Hardy and Wright section 16.7, or OEIS A000203). The API is identical to Pari/GP's "sigma" function. This function is useful for calculating things like aliquot sums, abundant num- bers, perfect numbers, etc. ramanujan_tau Takes a positive integer as input and returns the value of Ramanujan's tau function. The result is a signed integer. This corresponds to Mathemat- ica's "RamanujanTau" function and Pari's "ramanujantau" function. valuation say "$n is divisible by 2 ", valuation($n,2), " times."; Given integers "n" and "k", returns the numbers of times "|n|" is divisible by "k". This is a very limited version of the algebraic valuation -- here it is just applied to integers. "k" must be greater than 1. "|n| = 0" returns undef, and "|n| = 1" returns zero. This corresponds to Pari and SAGE's "valuation" function. is_qr Given two integers "a" and "n", returns 1 if "a" is a quadratic residue modulo "|n|", and 0 otherwise. A return value of 1 indicates there exists an "x" where "a = x^2 mod |n|". For odd primes, this is similar to checking "a==0 || kronecker(a,n) == 1". For all values, this will be equal to "sqrtmod(a,n) != undef", with possi- bly better performance. If "n = 0" the function returns undef. hammingweight Given an integer "n", returns the binary Hamming weight of abs(n). This is also called the population count, and is the number of 1s in the binary representation. This corresponds to Pari's "hammingweight" function for "t_INT" arguments. moebius say "$n is square free" if moebius($n) != 0; $sum += moebius($_) for (1..200); say "Mertens(200) = $sum"; say "Mertens(2000) = ", vecsum(moebius(0,2000)); Returns I1/4(n), the MA<paragraph>bius function (also known as the Moebius, Mobius, or MoebiusMu function) for an integer input. This function is 1 if "n = 1", 0 if "n" is not square free (i.e. "n" has a repeated factor), and "-1^t" if "n" is a product of "t" distinct primes. This is an important function in prime number theory. Like SAGE, we define "moebius(0) = 0" for convenience. If called with two arguments, they define a range "low" to "high", and the function returns an array with the value of the MA<paragraph>bius function for every n from low to high inclusive. invmod say "The inverse of 42 mod 2017 = ", invmod(42,2017); Given two integers "a" and "n", return the inverse of "a" modulo "n". If not defined, undef is returned. If defined, then the return value multi- plied by "a" equals 1 modulo "n". sqrtmod Given two integers "a" and "n", return the square root of "a" mod "|n|". If no square root exists, undef is returned. If defined, the return value "r" will always satisfy "r^2 = a mod |n|". If the modulus is prime, the function will always return "r", the smaller of the two square roots (the other being "-r mod |n|". If the modulus is composite, one of possibly many square roots will be returned, and it will not necessarily be the smallest. negmod Given two integers "a" and "n", return "-a mod |n|". This is similar to "submod(0,$a,$n)". addmod Given three integers "a", "b", and "n", return "(a+b) mod |n|". This is particularly useful when dealing with numbers that are larger than a half-word but still native size. No Perl bigint package is needed and this can be 10-200x faster than using one. Like other mod functions, undef is returned if "n=0". submod Given three integers "a", "b", and "n", return "(a-b) mod |n|". This is particularly useful when dealing with numbers that are larger than a half-word but still native size. No bigint package is needed and this can be 10-200x faster than using one. mulmod Given three integers "a", "b", and "n", return "(a*b) mod |n|". This is particularly useful when "n" fits in a native integer. No bigint package is needed and this can be 10-200x faster than using one. powmod Given three integers "a", "b", and "n", return "(a ** b) mod |n|". Typi- cally binary exponentiation is used, so the process is very efficient. With native size inputs, no bigint library is needed. divmod Given three integers "a", "b", and "n", return "(a/b) mod |n|". This is done as "(a * (1/b mod |n|)) mod |n|". If "n=0" or no inverse of "b mod |n|" exists then undef if returned. muladdmod Given four integers "a", "b", "c", and "n", return "(a*b+c) mod |n|". mulsubmod Given four integers "a", "b", "c", and "n", return "(a*b-c) mod |n|". consecutive_integer_lcm $lcm = consecutive_integer_lcm($n); Given an unsigned integer argument, returns the least common multiple of all integers from 1 to "n". This can be done by manipulation of the primes up to "n", resulting in much faster and memory-friendly results than using factorials. partitions Calculates the partition function p(n) for a non-negative integer input. This is the number of ways of writing the integer n as a sum of positive integers, without restrictions. This corresponds to Pari's "numbpart" function and Mathematica's "PartitionsP" function. The values produced in order are OEIS series A000041 <http://oeis.org/A000041>. This uses a combinatorial calculation, which means it will not be very fast compared to Pari, Mathematica, or FLINT which use the Rademacher formula using multi-precision floating point. In 10 seconds, the pure Perl version can produce partitions(10_000) while with Math::Prime::Util::GMP it can do partitions(220_000). In contrast, in about 10 seconds Pari can solve numb- part(22_000_000). If you want the enumerated partitions, see "forpart" in Math::Prime::Util or Integer::Partition. These are fast and memory efficient iterators, but not practical for producing the partition number for values over 100 or so. numtoperm @p = numtoperm(10,654321); # @p=(1,8,2,7,6,5,3,4,9,0) Given a non-negative integer "n" and integer "k", return the rank "k" lexi- cographic permutation of "n" elements. "k" will be interpreted as mod "n!". permtonum $k = permtonum([1,8,2,7,6,5,3,4,9,0]); # $k = 654321 Given an array reference containing integers from 0 to "n", returns the lexicographic permutation rank of the set. This is the inverse of the "numtoperm" function. All integers up to "n" must be present. The result will be between 0 and "n!-1". Pi Takes a positive integer argument "n" and returns the constant Pi with that many digits (including the leading 3). Rounding is performed. The implementation uses either AGM or Ramanujan/Chudnovsky with binary splitting, depending on the number of digits. It is a little over 2x faster than MPFR, similar in speed to Pari/GP, but about 1.5x slower than Xue's Chudnovsky demo from the GMP web site. Specialized programs such as "y-cruncher" are even faster. Note there is a non-trivial amount of overhead in turning the result into a string, as well as even more if using the ntheory module which further con- verts the result into a Math::BigFloat object. Called in void context, this just calculates and caches the result. Euler Takes a positive integer argument "n" and returns Euler's constant with that many digits. Rounding is performed. The implementation is Brent-McMillan algorithm B, just like Pari/GP. Per- formance is about 3x faster than Pari/GP, but 2-10x slower than MPFR which uses binary splitting. Called in void context, this just calculates and caches the result. exp_mangoldt say "exp(lambda($_)) = ", exp_mangoldt($_) for 1 .. 100; Returns EXP(I(n)), the exponential of the Mangoldt function (also known as von Mangoldt's function) for an integer value. The Mangoldt function is equal to log p if n is prime or a power of a prime, and 0 otherwise. We return the exponential so all results are integers. Hence the return value for "exp_mangoldt" is: p if n = p^m for some prime p and integer m >= 1 1 otherwise. totient say "The Euler totient of $n is ", totient($n); Returns I(n), the Euler totient function (also called Euler's phi or phi function) for an integer value. This is an arithmetic function which counts the number of positive integers less than or equal to "n" that are relatively prime to "n". Given the definition used, "totient" will return 0 for all "n < 1". This follows the logic used by SAGE. Mathematica and Pari return totient(-n) for "n < 0". Mathematica returns 0 for "n = 0", Pari pre-2.6.2 raises and exception, and Pari 2.6.2 and newer returns 2. jordan_totient say "Jordan's totient J_$k($n) is ", jordan_totient($k, $n); Returns Jordan's totient function for a given integer value. Jordan's to- tient is a generalization of Euler's totient, where "jordan_totient(1,$n) == euler_totient($n)" This counts the number of k-tuples less than or equal to n that form a coprime tuple with n. As with "totient", 0 is returned for all "n < 1". This function can be used to generate some other useful functions, such as the Dedekind psi function, where "psi(n) = J(2,n) / J(1,n)". carmichael_lambda Returns the Carmichael function (also called the reduced totient function, or Carmichael I>>(n)) of a positive integer argument. It is the smallest positive integer "m" such that "a^m = 1 mod n" for every integer "a" co- prime to "n". This is OEIS series A002322 <http://oeis.org/A002322>. liouville Returns I>>(n), the Liouville function for a non-negative integer input. This is -1 raised to I(C)(n) (the total number of prime factors). prime_bigomega say "$n has ", prime_bigomega($n), " total factors"; Returns I(C)(|n|), the prime omega function returning the total number of prime factors of n including multiplicities. The result is identical to "scalar(factor($n))". This corresponds to Pari's "bigomega" function and Mathematica's "PrimeOmega[n]" function. prime_omega say "$n has ", prime_omega($n), " distinct factors"; Returns I(|n|), the prime omega function returning the number of distinct prime factors of n. The result is identical to "scalar(factor_exp($n))". This corresponds to Pari's "omega" function and Mathematica's "PrimeNu[n]" function. is_divisible Given integers "n" and "d", returns 1 if "n" is exactly divisible by "d", and 0 otherwise. This corresponds to the GMP function "mpz_divisible_p". This includes its semantics with "d=0" which returns 0 unless "n=0". More than one divisor can be given, e.g. "is_divisible(1001,2,3,5,7)", al- lowing one to test multiple divisors with one call. The result is 1 if "n" is exactly divisible by any of the "d" values, and 0 if it is divisible by none of them. is_congruent Given integers "n", "c", and "d", returns 1 if "n" is congruent to "c" mod- ulo "d", and 0 otherwise. This corresponds to the GMP function "mpz_congruent_p". This includes its semantics with "d=0" which returns 0 unless "n=c". is_power say "$n is a perfect square" if is_power($n, 2); say "$n is a perfect cube" if is_power($n, 3); say "$n is a ", is_power($n), "-th power"; Given a single positive integer input "n", returns k if "n = p^k" for some integer "p > 1, k > 1", and 0 otherwise. The k returned is the largest possible. This can be used in a boolean statement to determine if "n" is a perfect power. If given two arguments "n" and "k", returns 1 if "n" is a "k-th" power, and 0 otherwise. For example, if "k=2" then this detects perfect squares. This corresponds to Pari/GP's "ispower" function, with the limitations of only integer arguments and no third argument may be given to return the root. is_square Given a positive integer "n", returns 1 if "n" is a perfect square, 0 oth- erwise. This is identical to "is_power(n,2)". This corresponds to Pari/GP's "issquare" function. is_prime_power Given an integer input "n", returns "k" if "n = p^k" for some prime p, and zero otherwise. This corresponds to Pari/GP's "isprimepower" function. sqrtint Returns the truncated integer part of the square root of "n". This corresponds to Pari/GP's "sqrtint" function. rootint Given "n" and "k", returns the truncated integer part of the "k-th" root of "n". This corresponds to Pari/GP's "sqrtnint" function. logint Given "n" and "b", returns the integer base-"b" logarithm of "n". This is the largest integer "e" such that "b^e <= n". This corresponds to Pari/GP's "logint" function. powint Given integers "a" and "b", returns "a^b". For <0^0> we return 1. The exponent "b" is converted into an unsigned long. mulint Given integers "a" and "b", returns "a * b". addint Given integers "a" and "b", returns "a + b". subint Given integers "a" and "b", returns "a - b". divint Given integers "a" and "b", returns the quotient "a / b". Floor division is used, so q is rounded towards -inf and the remainder has the same sign as the divisor "b". This is the same as modern "bdiv" in Math::BigInt and the GMP "fdiv" functions, but not the same as Pari/GP's "\\" operator. cdivint Given integers "a" and "b", returns the quotient "a / b" when using ceiling division. The quotient is rounded up towards +infinity and the remainder will have the opposite sign as the divisor "b". modint Given integers "a" and "b", returns the modulo "a % b". C<r = a - b * floor(a / b)> Floor division is used, so q is rounded towards -inf and r has the same sign as the divisor "b". This is the same as modern "bmod" in Math::BigInt and the GMP "fdiv" functions, but not the same as Pari/GP's "%" operator. divrem my($quo, $rem) = divrem($a, $b); Given integers "a" and "b", returns a list of two items: the Euclidean quo- tient and the Euclidean remainder. This corresponds to Pari/GP's "divrem" function. There is no explicit function in Math::BigInt that gives this division method for signed inputs. tdivrem Given integers "a" and "b", returns a list of two items: the truncated quo- tient and the truncated remainder. The resulting pair will match "btdiv" in Math::BigInt and "btmod" in Math::BigInt. fdivrem Given integers "a" and "b", returns a list of two items: the floor quotient and the floor remainder. cdivrem Given integers "a" and "b", returns a list of two items: the ceiling quo- tient and the ceiling remainder. lshiftint Given an integer "n" and an optional signed integer number of bits "k", perform a left shift of "n" by "k" bits. If the second argument is not provided, it is assumed to be 1. This is equivalent to multiplying by "2^k". This corresponds to Pari/GP's "shift" function with a positive number of bits, and Mathematica's "BitShiftLeft" function. rshiftint Given an integer "n" and an optional signed integer number of bits "k", perform a right shift of "n" by "k" bits. If the second argument is not provided, it is assumed to be 1. This is equivalent to truncated division by "2^k". With a negative "n", the result is equal to "-rshiftint(-n,k)". This corresponds to Pari/GP's "shift" function with a negative number of bits, and Mathematica's "BitShiftRight" function. The result is equal to dividing by the power of 2 using "tdivrem" or GMP's "mpz_tdiv_q_2exp". rashiftint Given an integer "n" and an optional signed integer number of bits "k", perform a signed arithmetic right shift of "n" by "k" bits. If the second argument is not provided, it is assumed to be 1. This is equivalent to floor division by "2^k". For non-negative "n", this is always equal to "rshiftint". With negative arguments it is similar to Math::BigInt#brsft, Python, and Java's BigInteger, which use floor division by "2^k". The result is equal to dividing by the power of 2 using "divint" or GMP's "mpz_fdiv_q_2exp". absint Given integer "n", return "|n|", i.e. the absolute value of "n". negint Given integer "n", return "-n". signint Given an integer "n", returns the sign of "n". Returns -1, 0, or 1 if "n" is negative, zero, or positive respectively. cmpint Given integers "a" and "b", returns -1, 0, or 1 if "a" is respectively less than, equal to, or greater than "b". cmpabsint As "cmpint" but using the absolute values of "a" and "b". add1int Given integer "n", return "n+1". sub1int Given integer "n", return "n-1". setbit Given integer "n" and a positive integer bit position "k", return "n" with bit "k" set. clrbit Given integer "n" and a positive integer bit position "k", return "n" with bit "k" cleared. notbit Given integer "n" and a positive integer bit position "k", return "n" with bit "k" complemented. tstbit Given integer "n" and a positive integer bit position "k", return the value of bit "k". bitand Given integers "a" and "b", return bitwise "a AND b". bitor Given integers "a" and "b", return bitwise "a OR b". bitxor Given integers "a" and "b", return bitwise "a XOR b". bitnot Given integer "n", return bitwise "NOT n". factor @factors = factor(640552686568398413516426919223357728279912327120302109778516984973296910867431808451611740398561987580967216226094312377767778241368426651540749005659); # Returns an array of 11 factors Returns a list of prime factors of a positive number, in numerical order. The special cases of "n = 0" and "n = 1" will return "n". Like most advanced factoring programs, a mix of methods is used. This in- cludes trial division for small factors, perfect power detection, Pollard's Rho, Pollard's P-1 with various smoothness and stage settings, Hart's OLF (a Fermat variant), ECM (elliptic curve method), and QS (quadratic sieve). Certainly improvements could be designed for this algorithm (suggestions are welcome). In practice, this factors 26-digit semiprimes in under "100ms", 36-digit semiprimes in under one second. Arbitrary integers are factored faster. It is many orders of magnitude faster than any other factoring module on CPAN circa 2013. It is comparable in speed to Math::Pari's "factorint" for most inputs. If you want better factoring in general, I recommend looking at the stand- alone programs yafu <http://sourceforge.net/projects/yafu/>, msieve <http://sourceforge.net/projects/msieve/>, gmp-ecm <http://ecm.gforge.in- ria.fr/>, and GGNFS <http://sourceforge.net/projects/ggnfs/>. divisors my @divisors = divisors(30); # returns (1, 2, 3, 5, 6, 10, 15, 30) Produces all the divisors of a positive number input, including 1 and the input number. The divisors are a power set of multiplications of the prime factors, returned as a sorted list with no duplications. The result is identical to that of Pari's "divisors" and Mathematica's "Divisors[n]" functions. In scalar context this returns the sigma0 function (OEIS A000005), and is the same result as evaluating the returned array in scalar context (but much more efficient). The result then corresponds to Pari's "numdiv" and Mathematica's "DivisorSigma[0,n]" functions. An optional second positive integer argument "k" indicates that the results should not include any value larger than "k". This is especially useful when the number has thousands of divisors and we may only be interested in the small ones. trial_factor my @factors = trial_factor($n); my @factors = trial_factor($n, 1000); Given a positive number input, tries to discover a factor using trial divi- sion. The resulting array will contain either two factors (it succeeded) or the original number (no factor was found). In either case, multiplying @factors yields the original input. An optional divisor limit may be given as the second parameter. Factoring will stop when the input is a prime, one factor is found, or the input has been tested for divisibility with all primes less than or equal to the limit. If no limit is given, then "2**31-1" will be used. This is a good and fast initial test, and will be very fast for small num- bers (e.g. under 1 million). For larger numbers, faster methods for com- plete factoring have been known since the 17th century. For inputs larger than about 1000 digits, a dynamic product/remainder tree is used, which is faster than GMP's native methods. This helps when prun- ing composites or looking for very small factors. prho_factor my @factors = prho_factor($n); my @factors = prho_factor($n, 100_000_000); Given a positive number input, tries to discover a factor using Pollard's Rho method. The resulting array will contain either two factors (it suc- ceeded) or the original number (no factor was found). In either case, mul- tiplying @factors yields the original input. An optional number of rounds may be given as the second parameter. Factoring will stop when the input is a prime, one factor has been found, or the number of rounds has been ex- ceeded. This is the Pollard Rho method with "f = x^2 + 3" and default rounds 64M. It is very good at finding small factors. Typically "pbrent_factor" will be preferred as it behaves similarly but runs quite a bit faster. They use different parameters however, so are not completely identical. pbrent_factor my @factors = pbrent_factor($n); my @factors = pbrent_factor($n, 100_000_000); Given a positive number input, tries to discover a factor using Pollard's Rho method with Brent's algorithm. The resulting array will contain either two factors (it succeeded) or the original number (no factor was found). In either case, multiplying @factors yields the original input. An op- tional number of rounds may be given as the second parameter. Factoring will stop when the input is a prime, one factor has been found, or the num- ber of rounds has been exceeded. This is the Pollard Rho method using Brent's modified cycle detection, de- layed "gcd" computations, and backtracking. It is essentially Algorithm P''2 from Brent (1980). Parameters used are "f = x^2 + 3" and default rounds 64M. It is very good at finding small factors. pminus1_factor my @factors = pminus1_factor($n); # Set B1 smoothness to 10M, second stage automatically set. my @factors = pminus1_factor($n, 10_000_000); # Run p-1 with B1 = 10M, B2 = 100M. my @factors = pminus1_factor($n, 10_000_000, 100_000_000); Given a positive number input, tries to discover a factor using Pollard's "p-1" method. The resulting array will contain either two factors (it suc- ceeded) or the original number (no factor was found). In either case, mul- tiplying @factors yields the original input. An optional first stage smoothness factor (B1) may be given as the second parameter. This will be the smoothness limit B1 for the first stage, and will use "10*B1" for the second stage limit B2. If a third parameter is given, it will be used as the second stage limit B2. Factoring will stop when the input is a prime, one factor has been found, or the algorithm fails to find a factor with the given smoothness. This is Pollard's "p-1" method using a default smoothness of 5M and a sec- ond stage of "B2 = 10 * B1". It can quickly find a factor "p" of the input "n" if the number "p-1" factors into small primes. For example "n = 22095311209999409685885162322219" has the factor "p = 3916587618943361", where "p-1 = 2^7 * 5 * 47 * 59 * 3137 * 703499", so this method will find a factor in the first stage if "B1 >= 703499" or in the second stage if "B1 >= 3137" and "B2 >= 703499". The implementation is written from scratch using the basic algorithm in- cluding a second stage as described in Montgomery 1987. It is faster than most simple implementations I have seen (many of which are written assuming native precision inputs), but slower than Ben Buhrow's code used in earlier versions of yafu <http://sourceforge.net/projects/yafu/>, and nowhere close to the speed of the version included with modern GMP-ECM with large B val- ues (it is actually quite a bit faster than GMP-ECM with small smoothness values). pplus1_factor my @factors = pplus1_factor($n); Given a positive number input, tries to discover a factor using Williams' "p+1" method. The resulting array will contain either two factors (it suc- ceeded) or the original number (no factor was found). In either case, mul- tiplying @factors yields the original input. An optional first stage smoothness factor (B1) may be given as the second parameter. This will be the smoothness limit B1 for the first stage. Factoring will stop when the input is a prime, one factor has been found, or the algorithm fails to find a factor with the given smoothness. cheb_factor my @factors = cheb_factor($n); my @factors = cheb_factor($n, 1_000); # set B1 smoothness Given a positive number input, tries to discover a factor using properties of Chebyshev polynomials (particularly that "T_mn(x) = T_m(T_n(x))") and their relationship with the Lucas sequence. This works when "p-1" or "p+1" is smooth. holf_factor my @factors = holf_factor($n); my @factors = holf_factor($n, 100_000_000); Given a positive number input, tries to discover a factor using Hart's OLF method. The resulting array will contain either two factors (it succeeded) or the original number (no factor was found). In either case, multiplying @factors yields the original input. An optional number of rounds may be given as the second parameter. Factoring will stop when the input is a prime, one factor has been found, or the number of rounds has been ex- ceeded. This is Hart's One Line Factorization method, which is a variant of Fer- mat's algorithm. A premultiplier of 480 is used. It is very good at fac- toring numbers that are close to perfect squares, or small numbers. Very naive methods of picking RSA parameters sometimes yield numbers in this form, so it can be useful to run a few rounds to check. For example, the number: 18548676741817250104151622545580576823736636896432849057 \ 10984160646722888555430591384041316374473729421512365598 \ 29709849969346650897776687202384767704706338162219624578 \ 777915220190863619885201763980069247978050169295918863 was proposed by someone as an RSA key. It is indeed composed of two dis- tinct prime numbers of similar bit length. Most factoring methods will take a very long time to break this. However one factor is almost exactly 5x larger than the other, allowing HOLF to factor this 222-digit semiprime in only a few milliseconds. squfof_factor my @factors = squfof_factor($n); my @factors = squfof_factor($n, 100_000_000); Given a positive number input, tries to discover a factor using Shanks' square forms factorization method (usually known as SQUFOF). The resulting array will contain either two factors (it succeeded) or the original number (no factor was found). In either case, multiplying @factors yields the original input. An optional number of rounds may be given as the second parameter. Factoring will stop when the input is a prime, one factor has been found, or the number of rounds has been exceeded. This is Daniel Shanks' SQUFOF (square forms factorization) algorithm. The particular implementation is a non-racing multiple-multiplier version, based on code ideas of Ben Buhrow and Jason Papadopoulos as well as many others. SQUFOF is often the preferred method for small numbers, and Math::Prime::Util as well as many other packages use it was the default method for native size (e.g. 32-bit or 64-bit) numbers after trial divi- sion. The GMP version used in this module will work for larger values, but my testing indicates it is generally slower than the "prho" and "pbrent" implementations. ecm_factor my @factors = ecm_factor($n); my @factors = ecm_factor($n, 12500); # B1 = 12500 my @factors = ecm_factor($n, 12500, 10); # B1 = 12500, curves = 10 Given a positive number input, tries to discover a factor using ECM. The resulting array will contain either two factors (it succeeded) or the orig- inal number (no factor was found). In either case, multiplying @factors yields the original input. An optional maximum smoothness may be given as the second parameter, which relates to the size of factor to search for. An optional third parameter indicates the number of random curves to use at each smoothness value being searched. This is an implementation of Hendrik Lenstra's elliptic curve factoring method, usually referred to as ECM. The implementation is reasonable, us- ing projective coordinates, Montgomery's PRAC heuristic for EC multiplica- tion, and two stages. It is much slower than the latest GMP-ECM, but still quite useful for factoring reasonably sized inputs. qs_factor my @factors = qs_factor($n); Given a positive number input, tries to discover factors using QS (the qua- dratic sieve). The resulting array will contain one or more numbers such that multiplying @factors yields the original input. Typically multiple factors will be produced, unlike the other "..._factor" routines. The current implementation is a modified version of SIMPQS, a predecessor to the QS in FLINT, and was written by William Hart in 2006. It will not operate on input less than 30 digits. The memory use for large inputs is more than desired, so other methods such as "pbrent_factor", "pminus1_fac- tor", and "ecm_factor" are recommended to begin with to filter out small factors. However, it is substantially faster than the other methods on large inputs having large factors, and is the method of choice for 35+ digit semiprimes. todigits Given an integer "n", return an array of digits of "|n|". An optional sec- ond integer argument specifies a base (default 10). For example, given a base of 2, this returns an array of binary digits of "n". An optional third argument specifies a length for the returned array. The result will be either have upper digits truncated or have leading zeros added. todigits(0) returns an empty array. The base must be at least 2, and is limited to an int. Length must be at least zero and is limited to an int. fromdigits Takes an array reference of base "b" digits in "MSB -" LSB> order, and an optional unsigned integer base "b" (default 10). Returns the equivalent base-10 integer. Alternately, takes a string of digits in base "b" instead of the the array reference, with both lower and upper case denoting digits 11 through 36. Each input digit must be an unsigned int. The base must be at least 2. seed_csprng Takes a non-negative integer "nbytes" and a string "data" as input. These are used to seed the internal CSPRNG used for random functions, including the random prime functions. Ideally this is 16-256 bytes of good entropy. Currently the CSPRNG is ISAAC-32, and the maximum number of seed bytes used is 1024. The CSPRNG will likely change to ChaCha20 in a later release. is_csprng_well_seeded Returns true if the CSPRNG has been seeded with 16 or more bytes (128 bits). There is no measurement of how "good" the input was. On startup the module will attempt to seed the CSPRNG from "/dev/urandom", so this function will return true if that was successful, but false other- wise. urandomb $n32 = urandomb(32); # Classic irand32, returns a UV $n = urandomb(1024); # Random integer less than 2^1024 Given a number of bits "b", returns a random unsigned integer less than "2^b". The result will be uniformly distributed between 0 and "2^b-1" in- clusive. This is similar to the GMP function "mpz_urandomb". urandomm $n = urandomm(100); # random integer in [0,99] $n = urandomm(1024); # random integer in [0,1023] Given a positive integer "n", returns a random unsigned integer less than "n". The results will be uniformly distributed between 0 and "n-1" inclu- sive. This is similar to the GMP function "mpz_urandomm". urandomr $n = urandomr(100, 110); # Random number [100,110] $nb = urandomr(2**24,2**25-1); # Random 25-bit number $nd = urandomr(10**24,10**25-1); # Random 25-digit number Given values "low" and "high", returns a uniform random unsigned integer in the range "[low,high]". Both inputs must be non-negative. If "low > high" then function will return "undef". Note that the range is inclusive, so "low", "high", and each integer between them have an equal probability of appearing. irand $n32 = irand; # random 32-bit integer Returns a random 32-bit integer using the CSPRNG. Performance is similar to "rand" in Math::Random::MTwist and Math::Ran- dom::Xorshift. It is somewhat faster than casting system "rand" to a 32-bit int. It is noticeably faster than Math::Random::ISAAC, Math::Ran- dom::ISAAC::XS, Math::Random::MT, Math::Random::MT::Auto, and Crypt::PRNG. irand64 $n64 = irand64; # random 64-bit integer Returns a random 64-bit integer using the CSPRNG (on 64-bit Perl). drand $f = drand; # random floating point value in [0,1) $r = drand(25); # random floating point value in [0,25) Returns a random NV (Perl's native floating point) using the CSPRNG. The number of bits returned is equal to the mantissa bits of the NV type used for the Perl build, with a max of 64. By default Perl uses doubles and the returned values have 53 bits. If Perl is built with long double support and the long doubles have a larger mantissa, then more bits are used. This gives substantially better quality random numbers than the default Perl "rand" function. Among other things, on modern Perl's, "rand" uses drand48, which gives 32 bits of decent random values and 16 more bits of known patterns (e.g. the 48th bit alternates, the 47th has a period of 4, etc.). There are much better choices for standard random number genera- tors, such as the Mersenne Twister from Math::Random::MTwist. Performance is similar to "rand" in Math::Random::MTwist and Math::Ran- dom::Xorshift. It is 1.5 - 2x slower than core "rand" (as are the other modules). random_bytes $str = random_bytes(32); # 32 random bytes Given an unsigned number "n" of bytes, returns a binary string filled with random data from the CSPRNG. Performance for getting 256 byte strings: Module/Method Rate Type ------------- --------- ---------------------- Data::Entropy::Algorithms 2027/s CSPRNG - AES Counter Crypt::Random 6649/s CSPRNG - /dev/urandom Bytes::Random 9217/s drand48 Bytes::Random::Secure 23043/s CSPRNG - ISAAC Math::Random::ISAAC::XS 58377/s CSPRNG - ISAAC rand+pack 82977/s drand48 Crypt::PRNG 298567/s CSPRNG - Fortuna Bytes::Random::XS 383354/s drand48 ntheory 770364/s CSPRNG - ChaCha20 Math::Random::MTwist 1890151/s Mersenne Twister Math::Prime::Util::GMP 2045715/s CSPRNG - ISAAC Each of the CSPRNG modules should be high quality. There are no known flaws in any of ISAAC, AES CTR, ChaCha20, or Fortuna. The industry seems to be standardizing on ChaCha20 (e.g. BSD, Linux, TLS). SEE ALSO Math::Prime::Util Has many more functions, lots of fast code for dealing with native-precision arguments (including much faster primes using sieves), and will use this module when needed for big numbers. Using Math::Prime::Util rather than this module directly is recommended. Math::Primality (version 0.08) A Perl module with support for the strong Miller-Rabin test, strong Lucas-Selfridge test, the BPSW probable prime test, next_prime / prev_prime, the AKS primality test, and prime_count. It uses Math::GMPz to do all the calculations, so is faster than pure Perl bignums, but a lot slower than XS+GMP. The prime_count function is only usable for very small inputs, but the other functions are quite good for big numbers. Make sure to use version 0.05 or newer. Math::Pari Supports quite a bit of the same functionality (and much more). See "SEE ALSO" in Math::Prime::Util for more detailed information on how the modules compare. yafu <http://sourceforge.net/projects/yafu/>, msieve <http://source- forge.net/projects/msieve/>, gmp-ecm <http://ecm.gforge.inria.fr/>, GGNFS <http://sourceforge.net/projects/ggnfs/> Good general purpose factoring utilities. These will be faster than this module, and much better as the factor increases in size. Primo <http://www.ellipsa.eu/public/primo/primo.html> is the state of the art in freely available (though not open source!) primality proving pro- grams. If you have 1000+ digit numbers to prove, you want to use this. mpz_aprcl <http://sourceforge.net/projects/mpzaprcl/> Open source APR-CL primality proof implementation. Fast primality proving, though without cer- tificates. REFERENCES Robert Baillie and Samuel S. Wagstaff, Jr., "Lucas Pseudoprimes", Mathemat- ics of Computation, v35 n152, October 1980, pp 1391-1417. <http://mpqs.free.fr/LucasPseudoprimes.pdf> Daniel J. Bernstein, "Proving Primality After Agrawal-Kayal-Saxena", preprint, Jan 2003. <http://cr.yp.to/papers/aks.pdf> Jon Grantham, "Frobenius Pseudoprimes", Mathematics of Computation, v70 n234, March 2000, pp 873-891. <http://www.ams.org/jour- nals/mcom/2001-70-234/S0025-5718-00-01197-2/> John Brillhart, D. H. Lehmer, and J. L. Selfridge, "New Primality Criteria and Factorizations of 2^m +/- 1", Mathematics of Computation, v29, n130, Apr 1975, pp 620-647. <http://www.ams.org/jour- nals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf> Richard P. Brent, "An improved Monte Carlo factorization algorithm", BIT 20, 1980, pp. 176-184. <http://www.cs.ox.ac.uk/peo- ple/richard.brent/pd/rpb051i.pdf> Peter L. Montgomery, "Speeding the Pollard and Elliptic Curve Methods of Factorization", Mathematics of Computation, v48, n177, Jan 1987, pp 243-264. <http://www.ams.org/jour- nals/mcom/1987-48-177/S0025-5718-1987-0866113-7/> Richard P. Brent, "Parallel Algorithms for Integer Factorisation", in Num- ber Theory and Cryptography, Cambridge University Press, 1990, pp 26-37. <http://www.cs.ox.ac.uk/people/richard.brent/pd/rpb115.pdf> Richard P. Brent, "Some Parallel Algorithms for Integer Factorisation", in Proc. Third Australian Supercomputer Conference, 1999. (Note: there are multiple versions of this paper) <http://www.cs.ox.ac.uk/peo- ple/richard.brent/pd/rpb193.pdf> William B. Hart, "A One Line Factoring Algorithm", preprint. <http://wstein.org/home/wstein/www/home/wbhart/onelinefactor.pdf> Daniel Shanks, "SQUFOF notes", unpublished notes, transcribed by Stephen McMath. <http://www.usna.edu/Users/math/wdj/mcmath/shanks_squfof.pdf> Jason E. Gower and Samuel S. Wagstaff, Jr, "Square Form Factorization", Mathematics of Computation, v77, 2008, pages 551-588. <http://homes.ce- rias.purdue.edu/~ssw/squfof.pdf> A.O.L. Atkin and F. Morain, "Elliptic Curves and primality proving", Mathe- matics of Computation, v61, 1993, pages 29-68. <http://www.ams.org/jour- nals/mcom/1993-61-203/S0025-5718-1993-1199989-X/> R.G.E. Pinch, "Some Primality Testing Algorithms", June 1993. Describes the primality testing methods used by many CAS systems and how most were compromised. Gives recommendations for primality testing APIs. <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.33.4409> AUTHORS Dana Jacobsen <dana@acm.org> Jason Papadopoulos wrote the "tinyqs" code which is basically unchanged. William Hart wrote SIMPQS which is the basis for the QS code. ACKNOWLEDGEMENTS Obviously none of this would be possible without the mathematicians who created and published their work. Eratosthenes, Gauss, Euler, Riemann, Fermat, Lucas, Baillie, Pollard, Brent, Montgomery, Shanks, Hart, Wagstaff, Dixon, Pomerance, A.K. Lenstra, H. W. Lenstra Jr., Atkin, Knuth, etc. The GNU GMP team, whose product allows me to concentrate on coding high-level algorithms and not worry about any of the details of how modular exponentiation and the like happen, and still get decent performance for my purposes. Ben Buhrow and Jason Papadopoulos deserve special mention for their open source factoring tools, which are both readable and fast. In particular I am leveraging their work on SQUFOF in the current implementation. They are a huge resource to the community. Jonathan Leto and Bob Kuo, who wrote and distributed the Math::Primality module on CPAN. Their implementation of BPSW provided the motivation I needed to do it in this module and Math::Prime::Util. I also used their module quite a bit for testing against. Paul Zimmermann's papers and GMP-ECM code were of great value for my pro- jective ECM implementation, as well as the many papers by Brent and Mont- gomery. COPYRIGHT Copyright 2011-2026 by Dana Jacobsen <dana@acm.org> This program is free software; you can redistribute it and/or modify it un- der the same terms as Perl itself. SIMPQS Copyright 2006, William Hart. SIMPQS is distributed under GPL v2+. perl v5.42.3 2026-03-08 Math::P...il::GMP(3)
NAME | VERSION | SYNOPSIS | DESCRIPTION | FUNCTIONS | SEE ALSO | REFERENCES | AUTHORS | ACKNOWLEDGEMENTS | COPYRIGHT
Want to link to this manual page? Use this URL:
<https://man.FreeBSD.org/cgi/man.cgi?query=Math::Prime::Util::GMP&sektion=3&manpath=FreeBSD+15.1-RELEASE+and+Ports.quarterly>