Math::...l::GMP(3pm)

🦋 Math::Prime::Util::GMP (3pm) — User Contributed Perl Documentation 🦋

💖 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 🫶

okay so bestie, this is a module for number theory in Perl using GMP. That means primality tests, grabbing primes in a range, factoring, and a whole bunch more. 💅

Like, you can use it directly, but its main character purpose is to serve Math::Prime::Util. That module will automatically load this one if it's installed, and it speeds up a ton of its big-number operations. No cap. 🔥

Inputs and outputs for big numbers? They're strings, so you don't need to drag a bigint package into your program. But if you do use bigints, inputs get converted internally anyway — no pre-conversion needed before a call. Output results come back as either Perl scalars (for native-size) or strings (for bigints). Math::Prime::Util tries to reconvert all strings back into the caller's bigint type if possible, which makes calculations way more convenient. ✨

Yeah, the is_*_pseudoprime tests are honestly better named is_*_probable_prime or is_*_prp. They hand you back 1 if the input is a probable prime according to their test. The naming is historical — it follows Pari, Math::Primality, and some other math packages. And FYI, the modern definition of pseudoprime is a composite that passes the test, not just any number. 💀

👑 FUNCTIONS 👑

🧚 is_prob_prime 🧚

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

okay so you give it a positive number and it hands back 0 (composite), 2 (definitely prime), or 1 (probably prime).

For inputs below 2^64, the test is deterministic — so your only possible returns are 0 (composite) or 2 (definitely prime). no maybes, bestie. 💅

For inputs above 2^64, it switches to a probabilistic test. You'll only see 0 (composite) and 1 (probably prime). The current implementation uses the Baillie-PSW (BPSW) test. Could a composite get marked prime? Theoretically yes, but since the test was published in 1980, not a single BPSW pseudoprime has been found — so it's extremely likely to be prime. We believe (Pomerance 1984) that an infinite number of counterexamples exist, but there's a weak conjecture (Martin) that none exist under 10000 digits. 😱

In more detail: we're using the extra-strong Lucas test (Grantham 2000) with the Baillie parameter selection method (see OEIS A217719). Previous versions of this module used the strong Lucas test with Selfridge parameters, but the extra-strong version produces fewer pseudoprimes while running 1.2 - 1.5x faster. It's 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);

You feed it a positive number and it returns 0 (composite), 2 (definitely prime), or 1 (probably prime). Composites act exactly like is_prob_prime, and so do numbers less than 2^64. For numbers bigger than 2^64, it runs some extra tests on probable primes to see if they can be proven by another means. 🔍

This call lives in the sweet spot between the speed of is_prob_prime and the certainty of is_provable_prime — those might be a better fit in some cases. What it does is give you most of the performance of the former while adding more certainty. Want to fine-tune that tradeoff? Use is_prob_prime followed by additional probable prime tests like miller_rabin_random and/or is_frobenius_underwood_pseudoprime.

Just like is_prob_prime, a BPSW test goes first. That's deterministic for all 64-bit numbers. Next, if the number is a Proth or LLR form, a proof gets 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 primality proof is attempted. This is tuned to give up if the result can't be quickly determined, and it 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.

So many numbers will return 2 (definitely prime), and the ones that return 1 (probably prime) have gone through more tests than is_prob_prime while not taking too long.

For cryptographic key generation, you might want even more testing for probable primes (NIST recommends a few more additional M-R tests than we perform). The function miller_rabin_random is made exactly for this. Alternately, a different test like is_frobenius_underwood_pseudoprime can be used. Even better? Use is_provable_prime — it should be reasonably fast for sizes under 2048 bits. Typically for key generation you want random primes, and there are many functions for that. 🔐

🌟 is_provable_prime 🌟

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

Give it a positive number and it returns 0 (composite), 2 (definitely prime), or 1 (probably prime). It takes a great deal of effort to return either 0 or 2 for all numbers. periodt. 💅

The current method first uses BPSW to find composites and give a deterministic answer for tiny numbers (under 2^64). If no certificate is required, LLR and Proth tests can be run, and small numbers (under approximately 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.

On a circa-2009 workstation, the time required for primes of different input sizes 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 take 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). ⏳

You can also get a certificate along with the result using the is_provable_prime_with_cert method. There's no appreciable extra performance cost for returning a certificate. 🎀

🦋 is_provable_prime_with_cert 🦋

Give it a positive number and it returns an array with two elements. 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 that contains all the necessary information to verify the primality of the input in a reasonable time. You can use the result with verify_prime in Math::Prime::Util for verification. Proof types used include:

ECPP
BLS3
BLS15
BLS5
Small

🍓 is_pseudoprime 🍓

Give it an integer n and zero or more positive bases. It returns 1 if the input is a probable prime to each base, 0 if not. This is the simple Fermat primality test — it checks that b^(n-1) = 1 mod n.

No bases given? Base 2 is used. All bases must be 2 or greater.

If you remove primes, with 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 🌈

Give it an integer n and zero or more non-zero positive bases. It 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 — it checks that kronecker(b,n) = b^((n-1)/2) mod n

No bases given? Base 2 is used. All bases must be 2 or greater.

If you remove primes, with 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, and 0 if it's not. Bestie this is literally the Miller-Rabin test, the classic, the icon 💅.

If you don't hand it any bases, it just uses base 2. And every base you do give has to be 2 or greater, no exceptions.

So if it returns 0? Yeah that number is 100% composite, no takebacks. But if it returns 1 😭 it's either an actual prime or a strong pseudoprime to all the bases you gave it. Give it enough distinct bases and the odds get extremely strong that it's genuinely prime. Slay.

Both the input number and the bases are allowed to be big integers. Go big or go home ✨.

This one is usually used in combo with other tests, either to build stronger tests (like the strong BPSW test) or to get deterministic results for numbers under some verified limit (for example, Jaeschke showed back in 1993 that no more than three selected bases are needed for correct primality results on any 32-bit number 👑). And because the chance of sailing through multiple bases is so tiny, some math packages literally just run multiple MR tests for primality. Though ngl, by the early 1990s almost all serious software had switched over to the BPSW test. It's giving glow-up.

n=2 always comes back as a strong pseudoprime, always. Every other even n returns 0 (most sources define this test for odd n only). Returning 0 for all even composite n is what makes this function match most other implementations, including Math::Primality's is_strong_pseudoprime function. Consistency queen 🫶.

If base modulo n <= 1 or base modulo n = n-1, then the result is gonna be 1. And listen, the test is not well defined if base = 0 mod n (like, Crandall and Pomerance page 135 vibes), so we just define it to return 1 — which lets bases be bigger than n while still giving meaningful results for all n. For example:

  is_strong_pseudoprime(367, 1101)

would return 0 if that change wasn't made. This matters big time for using deterministic base sets with large bases. Note that it only matters if the base you use is greater than or equal to n. And a 0 result will never come back for a prime n, no matter what base you throw at it. Periodt 💖.

---

🎲 miller_rabin_random 🎲

  my $maybe_prime = miller_rabin_random($n, 10); # 10 random bases

Takes an integer n as input plus a positive number k telling it how many bases to use. Then it runs k Miller-Rabin tests with uniform random bases between 2 and n-2. This is the correct way to do k Miller-Rabin tests 😤, not the super common but honestly broken method of just using the first k primes. Like, why would you do that.

You can also pass an optional third argument, which is a seed. The seed should be a number in decimal, or binary with a leading 0b, or hex with a leading 0x, or octal with a leading 0. It gets converted to a GMP integer, so it can be huge. Usually you don't need it, but cryptographic applications might prefer having it, and it lets you get repeatable test results. Fancy 🔥.

There is no check for duplicate bases. And inputs below 65 bits make little sense here since is_prob_prime is already deterministic at that size. For numbers 65+ bits the chance of duplicate bases is pretty 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. Math is wild 🌸.

---

🧚 is_lucas_pseudoprime 🧚

🧚 is_strong_lucas_pseudoprime 🧚

Given an integer n, returns 1 if the input is a standard or strong Lucas probable prime. It uses the Selfridge method of choosing D, P, and Q (some sources call this a Lucas-Selfridge test). This is one half of the BPSW primality test — the other half being the Miller-Rabin strong probable prime test with base 2. The canonical BPSW test (page 1401 of Baillie and Wagstaff (1980)) uses the strong Lucas test with Selfridge parameters, but in practice a whole variety of Lucas tests with different parameters get used by things calling themselves BPSW 👀.

The standard Lucas test implemented here corresponds to the Lucas test described in FIPS 186-4 section C.3.3, though it uses a slightly more efficient calculation. And 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 one for cryptographic purposes. It's often slightly faster, has over 4x fewer pseudoprimes, and it's the method Baillie and Wagstaff recommended in their 1980 paper. The strong one ate, fr fr 💅.

---

✨ is_extra_strong_lucas_pseudoprime ✨

Given an integer n, returns 1 if the input is an extra-strong Lucas probable prime. This is defined in Grantham (2000), and it's a slightly more stringent test than the strong Lucas test — though because different parameters are used, the pseudoprimes are not a subset. As the extra conditions would lead you to expect, 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. Fast and picky, main character energy 😌.

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. That makes it very fast — although as the input grows, the time converges to the conventional extra-strong implementation. At 30 digits this routine is about 15% faster, and at 300 digits it's only 2% faster. Diminishing returns, we've all been there 🙈.

With the current implementations, there's little reason to prefer this unless you're trying to reproduce specific results. The extra-strong implementation has been optimized to use similar features, which removed most of the performance advantage.

There's an optional second argument (it must be between 1 and 256) that sets the increment amount for P parameter selection. The default value of one gives you the method described in is_extra_strong_lucas_pseudoprime. A value of 2 gives you the method used in Pari (http://pari.math.u-bordeaux.fr/faq.html#primetest).

Because the U = 0 condition is ignored, this produces about 5% more pseudoprimes than the extra-strong Lucas test. But that's still only 66% of the number produced by the strong Lucas-Selfridge test. And no BPSW counterexamples have ever been found with any of the Lucas tests described here. Not a single one 😭.

---

💜 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 pseudoprime (Miller-Rabin) test.

The main reason this test exists is that it's very slightly more efficient than the other probable prime tests. Efficiency icon 💅.

---

🍓 is_perrin_pseudoprime 🍓

Given an integer n, returns 1 if n divides P(n), where P(n) is the Perrin number of n. The Perrin sequence is defined by P(n) = P(n-2) + P(n-3) with P(0) = 3, P(1) = 0, P(2) = 2. Cute recursion honestly.

This is not a commonly used test, since it runs slower than most of the other probable prime tests and offers little benefit — especially compared to combined tests like is_bpsw_prime, is_frobenius_underwood_pseudoprime, and is_frobenius_khashin_pseudoprime.

There's an optional second argument r that says whether to run additional tests. With r=1, P(-n) = -1 mod n is also verified, creating the "minimal restricted" test. With r=2, the full signature is also tested using the Adams and Shanks (1982) rules (without the quadratic form test). With r=3, the full signature is 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, returns 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 discriminant a^2-4b must not be a perfect square. Rules are rules 💅.

---

🌈 is_frobenius_underwood_pseudoprime 🌈

Given an integer n, returns 1 if n is positive and passes the efficient Frobenius test of Paul Underwood, and 0 otherwise. It 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. Which makes it similar to, but faster than, a standard Frobenius test. Understood the assignment 🔥.

This test is deterministic (no randomness is used). There are no known pseudoprimes to this test. And it has no overlap with the BPSW test either, which makes it a very effective method for adding extra certainty. Bestie duo energy, honestly 🫶.

---

👑 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. Not one, periodt ✨.

---

💖 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). Three whole answers, we love options.

This function does the extra-strong BPSW test and nothing more. Like, it skips all pretests and any extra work that the is_prob_prime test might add. That saves time if the input has no small factors — for example, when testing results that have already been sieved. Efficient queen 🌸.

---

🎀 is_trial_prime 🎀

Takes a positive number as input, and returns 1 if the input isn't divisible by any prime less than or equal to the square root of the input. This is simple trial division by primes. There's little practical use for it since we have way more efficient primality and compositeness tests, but it can be a useful comparison. It's giving "learner test", and honestly respect 🙈.

---

⭐ 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 unconditional primality test that runs in polynomial time for general input. Polynomial time, we love to see it 💅.

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. And it is, by a large order, faster than any other known implementation as of early 2017. Ate and left no crumbs 🔥👑.

For real though, bestie — for theoretical analysis of the primality thing, AKS is extremely important. In practice? It's lowkey useless, no cap. 😭 Estimated runtime for a 150 digit input is over 2 days, which is the whole case: yes the algorithmic complexity growth is polynomial, but the constants are extremely high. It'll take years for numbers that ECPP or APR-CL can prove in seconds. It's giving "technically correct but emotionally devastating" 💀

With the verbose option set to 1, the chosen r and s values get printed before the test even starts. With verbose set to 2 or higher, each of the s tests spits out a . as it runs, so you can literally watch the progress. 🫶

Typically you should just use is_provable_prime and let it decide the method. Let the algorithm pick, you don't need that stress. ✨

---

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 deterministic unconditional test that runs very fast compared to other primality methods for numbers of comparable size, and vastly faster than any known general-form primality proof methods. She's fast, she's sure of herself, she slays. 💅

is_llr_prime 🎀

Takes a positive number n as input and returns one of:

• 0 — definitely composite
• 2 — definitely prime
• -1 — test does not indicate anything

This implements the Lucas-Lehmer-Riesel test for fast deterministic primality 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 gets returned because the test just doesn't apply — not its business, not its problem. 🙈 If k = 1 then this is a Mersenne number and the Lucas-Lehmer test is used instead. Otherwise, the LLR test is performed.

While not as fast as the Lucas-Lehmer test for Mersenne numbers, it's almost as fast as a single strong pseudoprime test (i.e. Miller-Rabin test) while giving a certain answer. Certainty AND speed? Main character energy fr fr. 👑

is_proth_prime 🌸

Takes a positive number n as input and returns one of:

• 0 — definitely composite
• 2 — definitely prime
• -1 — test does not indicate anything

This applies Proth's theorem for fast Las Vegas primality 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 since the test does not apply. Otherwise, a search happens to find a quadratic nonresidue modulo n. If none can be found after a brief search, -1 is returned because no conclusion can be reached. Otherwise, Proth's theorem is checked, which conclusively indicates 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
• 2 — definitely prime

...to indicate whether the complex number a+bi is a Gaussian prime. The is_prime function is used internally to make the determination. iykyk 💜

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 deterministic Miller test. An optional second argument indicates whether the Generalized Riemann Hypothesis should be assumed, and it defaults to 0. Setting the verbose flag to 2 or higher will show how many bases are used. The unconditional test is exponential time, while the conditional test (assuming 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 it overflows a 64-bit integer at 456 bits — sizes that are trivial for the unconditional APR-CL and ECPP tests. So like... cute theory, rough reality. 💀

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 factoring 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 numbers 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 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 return definitely prime, probably prime, or definitely composite. Let the function do the thinking bestie, you deserve a break. 🫶

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 factoring 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 numbers 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 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 return 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 partial 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 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 return 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. She's the moment. She's the icon. 🌟

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 thousand digits. Having a larger set will help with large numbers (a set of 2650 is available on github in the xt/ directory). A future implementation 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 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 return 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. Same vibe, same shape. ✨

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 appropriate for small ranges and is useful for very large arguments. The Math::Prime::Util module has much more sophisticated methods for 64-bit arguments. Ngl, if you're living in 64-bit land, go check that one out. 💜

prime_count_lower

prime_count_upper 🌸

Returns lower or upper bounds for the prime count of the input n.

Bounds use Dusart 2010, Büthe 2014, Büthe 2015, and Axler 2017. Academic girlies really came through fr fr. 📚

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 interval (inclusive) as a list. It operates similar to primes, though it must always have a lower and upper bound and it 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 numbers beginning with n. An array of offsets from the start is returned.

The returned list contains those offsets in the range n to n+width-1 where n + offset has no prime factors less than or equal to depth. Clean list, no small factors, immaculate vibes. 🫧

okay so this one is literally the three-argument sieve_primes but make it ✨different✨ — it takes (n,width) instead of (low,high), and the bestie-exclusive plot twist: it hands you back small offsets from the start value instead of the actual values themselves. sounds minor, it is NOT — for multi-thousand digit numbers this cuts the overhead way down. 💅

🦋 sieve_twin_primes

my @primes = sieve_twin_primes(2**1000, 2**1000 + 500000);

hand it two arguments, low and high, and it gives you back every lower twin prime in that interval (inclusive, bestie, both ends). heads up: the result is a list, not a reference. no dereferencing, we're chilling. 🫶

how's it doing that? it does a partial sieve of the range, yeets out any non-twin candidates, then checks that each pair are both BPSW probable primes. the way this absolutely demolishes the naive approach of sieving for all primes and then filtering out the non-twins… iconic. substantially more efficient. periodt. 🔥

👑 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 sitting between the first two arguments, low and high (inclusive again, we stay consistent). the remaining arguments are the description of the cluster. house rules for those cluster values: they must be even, less than 31 bits, and strictly increasing. no exceptions, no notes. 😤

so given a cluster set C, the returned values are all the primes in the range where p+c is prime for all c in the cluster set C.

the cluster gets described as offsets from 0, with the implicit prime sitting at 0. so:

• an empty list is you asking for all primes — that's the cluster p+0
• a list with just the single value 2 finds all twin primes — that's the cluster where p+0 and p+2 are prime
• the list 2,6,8 finds prime quadruplets

and here's the freeing part: there is no requirement that your list denote a constellation (a cluster with minimal distance). the list 42,92,606? just fine. valid. we respect it. 🙈

for long clusters — think OEIS series A213601 (http://oeis.org/A213601), the prime 12-tuplets — this is immensely more efficient than filtering the cluster out of a list of primes. for that example, a range of 10^13 takes less than a second to search. that's thousands of times faster than filtering results from primes or twin primes. not me being dramatic, that's just the numbers. 💀

lowkey real talk tho: shorter clusters are not quite this efficient, and the overhead for returning large arrays should not be ignored. ship responsibly bestie.

---

🌸 next_prime

$n = next_prime($n);

returns the prime following the input number — specifically the smallest prime number that is greater than the input number. the function is_prob_prime is what decides when a prime has been found, so the result is a probable prime (using BPSW).

for large inputs this is quite a bit faster than GMP's mpz_nextprime or Pari's nextprime. slay. ✨

🌈 prev_prime

$n = prev_prime($n);

returns the prime preceding the input number — the largest prime number that is less than the input number. and if your input is 2 or lower? you get 0. that's the deal. 🙂‍↕️

same as above, is_prob_prime is what decides when a prime is found, so the result is a probable prime (using BPSW).

🧚 surround_primes

($dprev, $dnext) = surround_primes($n);

gives you the distances to the previous and next primes of the input n. it's slightly more efficient than calling both prev_prime and next_prime, and returning the distances as native integers is more efficient with large inputs. multitasking, but make it math. 💅

now — if you pass an optional second argument d, and the input n is larger than 2^64, then the whole vibe changes:

• if a SPSP-2 is found in the range n-d to n+d (inclusive), it will be returned with the other argument set to 0
• otherwise, the first SPSP-2 values found are returned

this feature is especially useful for prime gap searches as well as finding the nearest prime to a value. very main character of it honestly.

⚠️ not me warning you but: with a non-zero second argument, the values returned have not undergone a full BPSW test — just sieving and a SPSP-2 test. know that before you build anything on top of it, bestie. 😭

💖 next_twin_prime

returns the start of the next twin prime after the input n. and the returned value will always be greater than the input, no funny business.

so for a return value of t, both t and t+2 will be a probable prime (using BPSW). matching set, matching energy. 🔥

---

⭐ random_nbit_prime

say "random 512-bit prime: ", random_nbit_prime(512);

returns a randomly selected prime of exactly n bits. undef gets returned if n is less than 2 — so don't come at it with nonsense. the returned prime has passed the is_prob_prime (extra strong BPSW) test. certified-ish. 👑

🎀 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. and p is guaranteed to be a prime if q is prime. conditional slay, but slay nonetheless.

setting verbose level to 3 or higher produces progress output not unlike openssl:

• a . for each candidate that passes pretests
• a + for those where one is likely prime
• a * when both are likely prime with only confirmation tests remaining

it generates safe primes about 4-10x faster than openssl's dhparam. not a competition but we won. 💅

🦋 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, then p+1, q=p-1, and q-1 will all have a large factor. why do you care? because that makes using factoring methods such as p-1 and p+1 much harder. Gordon's algorithm is used.

lowkey honesty hour: the value of using strong primes is questionable over proper random primes when the number of bits is at least 1024. just so you know. 🙈

🌸 random_ndigit_prime

say "random 200-digit prime: ", random_ndigit_prime(200);

returns a randomly selected prime of exactly n digits. undef is returned if n is less than 1. returned prime has passed the is_prob_prime (extra strong BPSW) test. digits not bits, stay woke bestie. 💖

✨ 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 there is no prime 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 — so don't hardcode your assumptions about the RNG, okay?

this corresponds to Mathematica's RandomPrime[{min,max}] function. and it's a superset of Pari's randomprime(n) function, where our interval API is more convenient for cryptographic functions. glow up. 🔥

👑 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 — note this — not every n-bit prime can be generated with this algorithm.

undef is returned if n is less than 2. and as a safety check, internally the extra strong BPSW test is additionally run on each intermediate and the final result. double-checked, bestie. 🫶

🌈 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. and yeah it's 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 returned primes. tradeoffs, they exist.

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 it also returns a string certificate. two for one. periodt. 🎀

💜 random_shawe_taylor_prime_with_cert

like random_shawe_taylor_prime but it also returns a string certificate. matching energy with the one above. ✨

---

🍓 lucasu

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

given integers P, Q, and the non-negative integer k, it 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. a whole friend group of sequences honestly. 💖

this corresponds to OpenPFGW's lucasU function and gmpy2's lucasu function.

🌸 lucasv

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

given integers P, Q, and the non-negative integer k, computes V_k for the Lucas sequence defined by P,Q. these include the Lucas numbers (1,-1). iconic lineage. 🦋

this corresponds to OpenPFGW's lucasV function and gmpy2's lucasv function.

🌈 lucasuv

given integers P, Q, and the non-negative integer k, returns (U_k modulo n, V_k modulo n) for the Lucas sequence defined by P,Q. both of them, together, at once. main character duo. 💅

⭐ 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. note the |n| — absolute value, so negative n isn't a problem, bestie.

🎀 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 one can be particularly efficient, especially if Q=1. remember that. write it down. 🙈

🧚 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 defined by P,Q. the full combo meal. 👑

💜 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 — so this is a whole backbone function, lowkey. no cap.

house rules: k must be non-negative, and n must be greater than zero. P and Q are restricted to native signed integers. and the newer function lucasuvmod accepts bigint values — so if you need big, go there. 🫶

---

🔥 primorial

$p = primorial($n);

given an unsigned integer argument, returns the product of the prime numbers which are less than or equal to n. this definition of n# follows OEIS series A034386 (http://oeis.org/A034386) and Wikipedia: Primorial definition for natural numbers (http://en.wikipedia.org/wiki/Primorial#Definition_for_natural_numbers). receipts included. ✨

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

the two are related with the relationships:

pn_primorial($n)  ==   primorial( nth_prime($n) )
primorial($n)     ==   pn_primorial( prime_count($n) )

look at them. a matched set. they ate that. 💅

💜 factorial

given positive integer argument n, returns the factorial of n, defined as the product of the integers 1 to n — with the special case of factorial(0) = 1. yes, zero gets the vip treatment. 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 double 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. so:

multifactorial(n,2) = n * (n-2) * (n-4) * ...

see how it just … skips every other one? obsessed with the efficiency of it honestly. 🦋

The multifactorials are the OEIS series — for m=1, that's OEIS series A000142 (<http://oeis.org/A000142>), m=2 is A006882 (<http://oeis.org/A006882>), m=3 is A007661 (<http://oeis.org/A007661>), m=4 is A007662 (<http://oeis.org/A007662>), and it just keeps going, bestie.

Also pls go read Peter Luschny's absolutely iconic writeup (<http://oeis.org/wiki/User:Peter_Luschny/Multifactorials>) — it's giving scholarship. 🎀

---

factorialmod

Gimme a positive integer n and an integer m, and you get back n! mod |m| 💅. This is way faster than computing the whole giant factorial(n) and then doing a mod operation afterward, like some kind of amateur.

Under the hood it's a modular product of prime powers, with like powers grouped together so we do fewer exponentiations. Which is pretty efficient ngl — but the fast multipoint evaluation method (see Johannson (2012), the one FLINT uses) is significantly better. We're not that. We're still cute though. 🦋

---

factorial_sum

Hand it a positive integer n and it returns the factorial sum of n. That's the sum of factorial(k) for 0 .. n-1. These two are literally equivalent, except this function is faster:

factorial_sum($n) == vecsum(map{ factorial($_) } 0..$n-1)

Sometimes people call this the left factorial, which is confusing because that name is ALSO used for subfactorial. Chaos. Anyway, it's OEIS series A003422 (<http://oeis.org/A003422>). ✨

---

subfactorial

Give it a positive integer n and you get the subfactorial of n back 🫶. It's also known as the derangement number, and occasionally as the left factorial (yes, again, ik).

It's OEIS series A000166 (<http://oeis.org/A000166>), and it lines up with Mathematica's Subfactorial[n] function. No cap. 💜

---

falling_factorial

Two integers, x and n, where n is non-negative — and you get the falling factorial of n:

falling_factorial(x,n) = x * (x-1) * (x-2) * ... * (x-(n-1))

Matches Mathematica's FactorialPower[x,n] function. 🌸

---

rising_factorial

Same deal, two integers x and n, n non-negative — this one gives you the rising factorial of n:

rising_factorial(x,n) = x * (x+1) * (x+2) * ... * (x+(n-1))

Matches Mathematica's Pochhammer[x,n] function. Iconic behavior honestly. ⭐

---

gcd

Give it a list of integers, get the greatest common divisor back. Super commonly used to test for coprimality (<https://oeis.org/wiki/Coprimality>). And heads up — an empty list will return 0. 👑

---

lcm

List of integers in, least common multiple out. An empty list will return 1. Periodt. 💖

---

gcdext

Two integers x and y go in, and you get u,v,d out such that d = gcd(x,y) and u*x + v*y = d. It uses the extended Euclidian algorithm to find the values that satisfy Bézout's Identity. Math girlies stay winning. 🧚

This corresponds to Pari's gcdext function, which got renamed from bezout back in Pari 2.6. So the results will 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, bc we're thorough 💅). You feed it a list of [a,n] pairs, and each one represents an equation x ≡ a mod |n|.

• If no solution exists → undef comes back 😭
• If a solution DOES exist → the modulus equals the lcm of all the given moduli (see lcm)
• In the standard case where every n is coprime, that's just the product

The n values must be non-zero integers, while the a values are integers. Keep it together, bestie. 🌈

---

chinese2

Works just like chinese, but returns two items instead: the remainder and the modulus. If a solution exists, that second value (the final modulus) equals the lcm of the absolute values of all the given moduli.

And if no solution exists, both return values will be undef. Both of them. Gone. 🙈

---

vecsum

Returns the sum of all arguments — each one has to be an integer. 💜

---

vecprod

Returns the product of all arguments — each one has to be an integer. 💜

---

kronecker

Returns the Kronecker symbol (a|n) for two integers. Here's what the possible return values mean when n is an odd prime:

• 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 extends the Jacobi symbol to all integer values of n, whereas the Jacobi symbol only lived over positive odd values of n. And the Jacobi symbol is itself an extension of the Legendre symbol, which is only defined for odd prime values of n. It's extensions all the way down, fr fr. 🦋

This corresponds to Pari's kronecker(a,n) function and Mathematica's KroneckerSymbol[n,m] function.

⚠️ Real talk: if n is not an odd prime, the result does not necessarily tell you whether a is a quadratic residue mod n. If you need correct results for any n, use is_qr instead — though it could be slower. 💀

---

binomial

Integer arguments n and k in, binomial coefficient out: n*(n-1)*...*(n-k+1)/k!, aka the choose function. Negative arguments use the Kronenburg extensions (<http://arxiv.org/abs/1105.3689/>). This corresponds to Mathematica's Binomial[n,k] function, Pari's binomial(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. A few notes so you're not caught off guard:

• Pari does not implement the n < 0, k <= n extension — it just returns 0 for that case.
• GMP's API doesn't allow negative k, but otherwise matches.
• Math::BigInt version 1.999816 and later supports negative arguments with similar semantics. Before that, n < 0, k > 0 was undefined. 💀

---

binomialmod

Give it integer arguments n, k, and m, and it efficiently returns binomial(n,k) mod m ✨.

All the requirements for binomial must be met, including the size of k and how negative n and k are treated.

The current implementation is not optimized, but it can greatly reduce the overhead that comes with the very large intermediate. Which is a big W. 🫶

---

addreal, subreal, mulreal, divreal

Each of these returns the corresponding basic math operation applied to the two inputs. There's an optional third argument that sets the number of significant digits (default 40), with the result rounded. 💅

---

logreal

Returns the natural logarithm of the input n. An optional second argument sets the number of significant digits (default 40), with the result 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-wise it's 100-1000x faster than Math::BigFloat's GMP backend — but it's 10x slower than Pari/GP 2.10 and MPFR. We're fast, just not that fast.

Negative inputs come back as -log(-n), which matches bignum. Meanwhile Pari/GP and Mathematica return log(-n) + Pi*i. Different vibes, different answers. 🦋

---

expreal

Returns e^n for the input n. An optional second argument sets the number of significant digits (default 40), with the result rounded. ⭐

The implementation computes sinh(n), then pulls e^x out of that.

---

powreal

Returns x^n for the inputs x and n. An optional third argument sets the number of significant digits (default 40), with the result rounded. 💜

Just like logreal and expreal, this is a basic math function that the GMP library doesn't provide but MPFR implements. Since MPFR isn't always around, it's genuinely nice to have this.

⚠️ Two ways to get undef back: a negative x with a non-integer n (the result is complex), or x=0 with a negative n (that's dividing by zero, bestie, absolutely not). 💀

---

rootreal

Returns the n-th root of x for the inputs n and x. An optional third argument sets the number of significant digits (default 40), with the result rounded. 🎀

This is basically powreal(x^(1/n), except it lets you specify n in full precision easily.

⚠️ Asking for the 0-th root returns undef. So does a negative x with roots other than 1 or -1. 💀

---

agmreal

Returns the Arithmetic-Geometric mean (AGM) of a and b. An optional third argument sets the number of significant digits (default 40), with the result rounded. 🌈

Use cases are lowkey everywhere: elementary constants (like Pi and e), logs, exponentials, trigonometric functions, elliptic integrals, computing pendulum periods, and more. She does it all. ✨

This corresponds to Pari's agm(x,y) function, limited to positive reals (Pari also handles negative, complex, p-adic, and power series arguments — we do not).

---

bernfrac

Returns the Bernoulli number B_n for an integer argument n, as a rational number. Two values come back: the numerator and the denominator. And importantly, 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 sets the number of significant 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. Every element of that array is a 2-element array reference holding the numerator and the denominator. 🎀

This will be faster than calling bernfrac for each value individually. The computed values get stashed in a cache (and it's not incremental), so later calls to bernfrac, bernreal, bernvec, and powersum with an n in that range will be substantially faster. Free speed, we love to see it.

When you call it in void context, it just calculates and caches the values. Chill.

⚠️ Heads up: 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 rational number. Two values come back: the numerator and the denominator. The numbers are the sum of the 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 sets the number of digits preserved past the decimal place, and the default is 40. ✨

---

stirling

say "s(14,2) = ", stirling(14, 2);
say "S(14,2) = ", stirling(14, 2, 2);

Returns Stirling numbers of the first kind (the default), the second kind, or the third kind (the unsigned Lah numbers) — you pick the kind with an optional third argument. It takes two non-negative integer arguments, n and k, plus that optional type. This corresponds to Pari's stirling(n,k,{type}) function and Mathematica's StirlingS1 / StirlingS2 functions. 👑

So what are they actually counting? 💅

• Stirling numbers of the first kind are -1^(n-k) times the number of permutations of n symbols with exactly k cycles.
• Stirling numbers of the second kind are the number of ways to partition a set of n elements into k non-empty subsets.
• The Lah numbers are the number of ways to split a set of n elements into k non-empty lists.

Combinatorics girlies, this one's for you. 🌸🦋

zeta ✨

Hand it a positive integer or float n and it hands you back the real Riemann Zeta value as a string floating point. The optional second argument is how many digits past the decimal point you want — default is 40.

The implementation is algorithm 2 of Borwein (1991). Performance with integer 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 🦋

Give it a positive integer or float n, get the real Logarithmic Integral back as a string floating point. Optional second argument = number of significant digits (default 40) with the result rounded.

The implementation uses Ramanujan's series. This corresponds to Mathematica's Li function. iykyk 💅

ei 🌸

Positive integer or float n goes in, real Exponential Integral comes out as a string floating point. Optional second argument = number of significant digits (default 40) with the result rounded.

The implementation is simply li(exp(x)). Not even joking, that's it, it's giving one-liner 💀

riemannr ⭐

Positive integer or float n → the real Riemann R function as a string floating point. Optional second argument = number of significant digits (default 40) with the result rounded.

The implementation is the standard Gram series. This corresponds to Mathematica's RiemannR function.

lambertw 🧚

Give it a float x and you get 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 Mathematica's ProductLog / LambertW function.

znorder 👑

  $order = znorder(17, "100000000000000000000000065");

Given two positive integers a and n, returns the multiplicative order of a modulo n. That's the smallest positive integer k such that a^k ≡ 1 mod n. Returns 1 if a = 1. Returns undef if a = 0 or if a and n are not coprime, because bestie, no value will ever result in 1 mod n in that case. This corresponds to Pari's znorder(Mod(a,n)) function and Mathematica's MultiplicativeOrder[a,n] function.

znprimroot 🌈

Given a positive integer n, returns the smallest primitive root of (Z/nZ)^*, or undef if no root exists. A root exists when 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 totient(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. tiny, iconic, niche 🫶

The boolean result is the same as scalar(factor(n)) == 2, but this function 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 values 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 divisors 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) ≡ 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.

Heads up bestie: Inputs greater than 50 digits use a probabilistic test to avoid fully factoring the input. Not a vibe for exactness girlies.

is_fundamental 🌟

Given a positive integer n, returns 1 if n is a fundamental discriminant, 0 otherwise.

is_totient 🌸

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

is_polygonal ⭐

Given integers x and s, return 1 if x is an s-gonal number, 0 otherwise. 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 integer 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 defaults 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+. Following 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 exponent 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 successive 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.

okay listen up, this one has ✨consequences✨: 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 general. 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 numbers, perfect numbers, etc. 💅

🎀 ramanujan_tau 🎀

Give it a positive integer and it hands you back Ramanujan's tau function value — a signed integer, bestie, could be negative, don't @ me. This is the same energy as Mathematica's RamanujanTau and Pari's ramanujantau.

🦋 valuation 🦋

say "$n is divisible by 2 ", valuation($n,2), " times.";

Hand it integers n and k and it tells you how many times |n| is divisible by k. Heads up though — this is a very limited version of the algebraic valuation, like it's only applied to integers here, no fancy generalization.

k has to be greater than 1, periodt. |n| = 0 returns undef, and |n| = 1 returns zero.

This is giving Pari and SAGE's valuation function, iykyk.

✨ is_qr ✨

Give it two integers a and n and you get 1 if a is a quadratic residue modulo |n|, and 0 if it's not. A return of 1 literally means there exists some x where a = x^2 mod |n|. 💅

For odd primes, this is basically checking a==0 || kronecker(a,n) == 1.

For every value, this equals sqrtmod(a,n) != undef, just possibly with better performance. Slay.

If n = 0, the function returns undef. 🙈

🌸 hammingweight 🌸

Give it an integer n, get the binary Hamming weight of abs(n). It's also called the population count, and it's just the number of 1s in the binary representation. This matches 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 μ(n), the Möbius function — also known as the Moebius, Mobius, or MoebiusMu function, it has main character energy with all those names — for an integer input.

The rules: it's 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 a genuinely important function in prime number theory, no cap. Like SAGE, we define moebius(0) = 0 for convenience.

If you call it with two arguments, they define a range low to high, and it returns an array with the value of the Mö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, returns the inverse of a modulo n. If it's not defined, undef comes back. If it is defined, then the return value multiplied by a equals 1 modulo n. Math ate, honestly.

🧚 sqrtmod 🧚

Given two integers a and n, returns the square root of a mod |n|. If no square root exists, you get undef. If it is 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 one being -r mod |n|). If the modulus is composite, you get one of possibly many square roots, and it will not necessarily be the smallest. So don't come for it. 💀

🌈 negmod 🌈

Given two integers a and n, returns -a mod |n|.

This is similar to submod(0,$a,$n).

⭐ addmod ⭐

Given three integers a, b, and n, returns (a+b) mod |n|. This is highkey useful when you're dealing with numbers bigger than a half-word but still native size. No Perl bigint package needed and this can be 10-200x faster than using one. 🔥

Like the other mod functions, undef is returned if n=0.

💖 submod 💖

Given three integers a, b, and n, returns (a-b) mod |n|. Again, super useful for numbers larger than a half-word but still native size. No bigint package needed and this can be 10-200x faster than using one.

🎀 mulmod 🎀

Given three integers a, b, and n, returns (a*b) mod |n|. This one especially shines when n fits in a native integer. No bigint package needed and this can be 10-200x faster than using one.

🦋 powmod 🦋

Given three integers a, b, and n, returns (a ** b) mod |n|. Typically binary exponentiation is used, so the process is very efficient — obsessed. With native size inputs, no bigint library is needed at all.

✨ divmod ✨

Given three integers a, b, and n, returns (a/b) mod |n|. This gets done as (a * (1/b mod |n|)) mod |n|. If n=0 or no inverse of b mod |n| exists, then undef is returned. 🙈

🌸 muladdmod 🌸

Given four integers a, b, c, and n, returns (a*b+c) mod |n|.

💜 mulsubmod 💜

Given four integers a, b, c, and n, returns (a*b-c) mod |n|.

🍓 consecutive_integer_lcm 🍓

$lcm = consecutive_integer_lcm($n);

Give it an unsigned integer argument and it returns the least common multiple of all integers from 1 to n. It pulls this off by manipulation of the primes up to n, which means much faster and more memory-friendly results than using factorials. It's giving efficiency queen. 👑

🧚 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. It corresponds to Pari's numbpart function and Mathematica's PartitionsP function. The values produced in order are OEIS series A000041.

Now listen — 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 numbpart(22_000_000). So like, know what you're signing up for. 😭

If you want the enumerated partitions, see forpart in Math::Prime::Util or Integer::Partition. Those 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 lexicographic permutation of n elements. k gets interpreted as mod n!. Simple as that, bestie.

⭐ 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 — they're a matched set, periodt. 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's 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's a non-trivial amount of overhead in turning the result into a string, and even more if you're using the ntheory module which further converts the result into a Math::BigFloat object. It's giving extra steps, honestly.

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. Performance 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(Λ(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 φ(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 an exception, and Pari 2.6.2 and newer returns 2. Everyone's doing their own thing, it's giving chaos. 💀

🌸 jordan_totient 🌸

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

Returns Jordan's totient function for a given integer value. Jordan's totient is a generalization of Euler's totient, where jordan_totient(1,$n) == euler_totient($n). It 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 λ(n)) of a positive integer argument. It is the smallest positive integer m such that a^m = 1 mod n for every integer a coprime to n. This is OEIS series A002322.

🍓 liouville 🍓

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

🧚 prime_bigomega 🧚

say "$n has ", prime_bigomega($n), " total factors";

Returns Ω(|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 ω(|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), which lets you test multiple divisors with one single call. The result is 1 if n is exactly divisible by any of the d values, and 0 if it's divisible by none of them. One call, many divisors, we love efficiency. 💅

💖 is_congruent 💖

Given integers n, c, and d, returns 1 if n is congruent to c modulo 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. And bestie, the k it hands back is the largest one possible, no contest. So you can drop this straight into a boolean statement to figure out whether n is a perfect power. ✨

If you give it two arguments — n and k — it returns 1 if n is a k-th power and 0 if it's not. Like, if k=2 it's literally just detecting perfect squares. 💅

This is the everyday equivalent of Pari/GP's ispower function, but with some limits: only integer arguments allowed, and you can't pass a third argument to get the root back. Rip. 😭

is_square

Hand it a positive integer n, get 1 back if n is a perfect square, 0 otherwise. Yeah it's literally the same as is_power(n,2), we said what we said. 🎀

This matches up with Pari/GP's issquare function.

is_prime_power

Give it an integer n, and it returns k if n = p^k for some prime p, and zero if that's not the vibe.

This is the twin of Pari/GP's isprimepower function.

sqrtint

Returns the truncated integer part of the square root of n. Simple, clean, iconic. 🌸

Corresponds to Pari/GP's sqrtint function.

rootint

Give it n and k, and it returns the truncated integer part of the k-th root of n.

Corresponds to Pari/GP's sqrtnint function.

logint

Give it n and b, and you get the integer base-b logarithm of n. Translation: it's the largest integer e such that b^e <= n. 🧮

Corresponds to Pari/GP's logint function.

powint

Give it integers a and b, get a^b back. And for 0^0? We return 1. Periodt.

Fun little detail: the exponent b gets converted into an unsigned long. 👑

mulint

Give it integers a and b, get a * b back. Math but make it a function call.

addint

Give it integers a and b, get a + b back. 💖

subint

Give it integers a and b, get a - b back. No notes.

divint

Give it integers a and b, get the quotient a / b back.

Floor division is what's happening here, so q gets rounded towards -inf and the remainder keeps the same sign as the divisor b. This is exactly what modern bdiv in Math::BigInt does, and what GMP's fdiv functions do — but it is not the same as Pari/GP's \ operator. Don't get caught slipping. 🙈

cdivint

Give it integers a and b, and you get the quotient a / b under ceiling division. So the quotient rounds up towards +infinity, and the remainder ends up with the opposite sign of the divisor b.

modint

Give it integers a and b, get the modulo a % b back.

r = a - b * floor(a / b)

Floor division is the move here — q rounds towards -inf and r shares the same sign as the divisor b. That's identical to modern bmod in Math::BigInt and the GMP fdiv functions, but not the same as Pari/GP's % operator. Ik, ik, different vibes entirely.

divrem

my($quo, $rem) = divrem($a, $b);

Give it integers a and b, and it returns a list of two things: the Euclidean quotient and the Euclidean remainder. Two-for-one deal, slay. 🫶

This matches Pari/GP's divrem function. Note that Math::BigInt has no explicit function giving this division method for signed inputs — so this is a you-can't-get-this-elsewhere moment.

tdivrem

Give it integers a and b, get a list of two back: the truncated quotient and the truncated remainder.

That resulting pair will line up with btdiv in Math::BigInt and btmod in Math::BigInt.

fdivrem

Give it integers a and b, get a list of two: the floor quotient and the floor remainder.

cdivrem

Give it integers a and b, get a list of two: the ceiling quotient and the ceiling remainder. 🌈

lshiftint

Give it an integer n and an optional signed integer number of bits k, and it'll do a left shift of n by k bits. Not passing that second argument? It just assumes 1. This is equivalent to multiplying by 2^k. 🦋

This corresponds to Pari/GP's shift function when the number of bits is positive, and Mathematica's BitShiftLeft function.

rshiftint

Give it an integer n and an optional signed integer number of bits k, and it does a right shift of n by k bits. Skip that second arg and it assumes 1, same as before. This is equivalent to truncated division by 2^k. And with a negative n, the result equals -rshiftint(-n,k). 🍓

This corresponds to Pari/GP's shift function when the number of bits is negative, plus Mathematica's BitShiftRight function. The result is the same as dividing by that power of 2 using tdivrem or GMP's mpz_tdiv_q_2exp.

rashiftint

Give it an integer n and an optional signed integer number of bits k, and it performs a signed arithmetic right shift of n by k bits. Second argument missing? Assumed to be 1, as always. This is equivalent to floor division by 2^k. And for non-negative n, it's always equal to rshiftint.

With negative arguments it's similar to Math::BigInt#brsft, Python, and Java's BigInteger, all of which use floor division by 2^k. The result equals dividing by the power of 2 using divint or GMP's mpz_fdiv_q_2exp. 👍

absint

Give it integer n, get |n| back — the absolute value of n. Bye bye negativity. 💜

negint

Give it integer n, get -n back.

signint

Give it an integer n, get the sign of n. So it returns -1, 0, or 1 depending on whether n is negative, zero, or positive respectively. Vibe check in integer form.

cmpint

Give it integers a and b, get -1, 0, or 1 back depending on whether a is less than, equal to, or greater than b. 🎯

cmpabsint

Same as cmpint, but it uses the absolute values of a and b. The math is comparative, the insecurity is gone.

add1int

Give it integer n, get n+1 back. Tiny little +1, we love her.

sub1int

Give it integer n, get n-1 back. Tiny little -1, equally iconic.

setbit

Give it integer n and a positive integer bit position k, and it returns n with bit k set. 🔥

clrbit

Give it integer n and a positive integer bit position k, and it returns n with bit k cleared.

notbit

Give it integer n and a positive integer bit position k, and it returns n with bit k complemented.

tstbit

Give it integer n and a positive integer bit position k, and it returns the value of bit k. Just peeking, no touching. 👀

bitand

Give it integers a and b, get bitwise a AND b. ✨

bitor

Give it integers a and b, get bitwise a OR b.

bitxor

Give it integers a and b, get bitwise a XOR b. Sparks fly, bits flip.

bitnot

Give it integer n, get bitwise NOT n. Full send on the flipping. 💅

factor

@factors = factor(640552686568398413516426919223357728279912327120302109778516984973296910867431808451611740398561987580967216226094312377767778241368426651540749005659);
# Returns an array of 11 factors

Returns a list of prime factors of a positive number, sorted in numerical order. And the special cases? n = 0 and n = 1 will just return n. 🤷

Like basically every advanced factoring program out there, it mixes methods: trial division for small factors, perfect power detection, Pollard's Rho, Pollard's P-1 with various smoothness and stage settings, Hart's OLF (which is a Fermat variant), ECM (the elliptic curve method), and QS (quadratic sieve). And yes, there are certainly improvements that could be designed for this algorithm — suggestions are welcome, bestie, the inbox is open. 💌

In practice: 26-digit semiprimes get factored in under 100ms, and 36-digit semiprimes in under one second. Arbitrary integers go even faster. It's many orders of magnitude faster than any other factoring module on CPAN circa 2013, and it's comparable in speed to Math::Pari's factorint for most inputs. She ate, ngl.

Want better factoring in general? The recommendation on the table is to go look at the standalone programs yafu <http://sourceforge.net/projects/yafu/>, msieve <http://sourceforge.net/projects/msieve/>, gmp-ecm <http://ecm.gforge.inria.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 itself. Basically it's the power set of multiplications of the prime factors, returned as a sorted list with zero duplications. The result is identical to Pari's divisors and Mathematica's Divisors[n] functions. 🎀

Oh and in scalar context this returns the sigma0 function (OEIS A000005), and it's the same result as evaluating the returned array in scalar context — except way more efficient. That result then corresponds to Pari's numdiv and Mathematica's DivisorSigma[0,n] functions.

There's an optional second positive integer argument k that says the results shouldn't include any value larger than k. Super useful when your number has like a thousand divisors and you truly only care about the small ones. 💖

trial_factor

my @factors = trial_factor($n);
my @factors = trial_factor($n, 1000);

Given a positive number input, it tries to discover a factor using trial division. The resulting array will contain either two factors (it succeeded, we love that) or the original number (no factor was found, rip). Either way, multiplying @factors gets you back the original input. 🫶

You can hand in an optional divisor limit as the second parameter. Factoring stops when the input is a prime, when one factor is found, or when the input has been tested for divisibility against all primes less than or equal to the limit. No limit given? Then 2**31-1 gets used.

This is a great and fast initial test, and it's very fast for small numbers (think under 1 million). For larger numbers, well, faster methods for complete factoring have been known since the 17th century, so. 👑

For inputs larger than about 1000 digits, a dynamic product/remainder tree is used, which is faster than GMP's native methods. That really helps when you're pruning composites or hunting for very small factors.

prho_factor

my @factors = prho_factor($n);
my @factors = prho_factor($n, 100_000_000);

Given a positive number input, it tries to discover a factor using Pollard's Rho method. The resulting array contains either two factors (it succeeded) or the original number (no factor was found). Either way, multiplying @factors yields the original input. You can pass an optional number of rounds as the second parameter. Factoring stops when the input is a prime, when one factor has been found, or when the number of rounds has been exceeded. 🎯

This is the Pollard Rho method with f = x^2 + 3 and a default of 64M rounds. It is very good at finding small factors. Typically pbrent_factor will be the preferred pick since it behaves similarly but runs quite a bit faster — though the two use different parameters, so they are not completely identical. 🌸

pbrent_factor

my @factors = pbrent_factor($n);
my @factors = pbrent_factor($n, 100_000_000);

Given a positive number input, it 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). Either way, 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. 🦋

okay so this one is the Pollard Rho method but make it ✨fancy✨ — it's using Brent's modified cycle detection, delayed gcd computations, and backtracking. it is literally Algorithm P''2 from Brent (1980), periodt 💅. the parameters it rolls with are f = x^2 + 3 and default rounds 64M. it is very good at finding small factors, no cap 🎀

---

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

hand it a positive number and it tries to discover a factor using Pollard's p-1 method. the resulting array will contain either two factors (bestie, it succeeded 🎉) or the original number (no factor was found 😭). either way, multiplying @factors yields the original input — the math always mathing 💯.

an optional first stage smoothness factor (B1) may be given as the second parameter. that becomes the smoothness limit B1 for the first stage, and it'll use 10*B1 for the second stage limit B2. give it a third parameter and that gets used as the second stage limit B2. factoring stops when the input is a prime, one factor has been found, or the algorithm fails to find a factor with the given smoothness 🫶

so this is Pollard's p-1 method using a default smoothness of 5M and a second 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. cheating? no, just efficient 💅

the implementation is written from scratch using the basic algorithm including a second stage as described in Montgomery 1987. it's faster than most simple implementations i've seen (many of which are written assuming native precision inputs — not us, we do it properly 💖), 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 values (it is actually quite a bit faster than GMP-ECM with small smoothness values though, slay).

---

pplus1_factor

my @factors = pplus1_factor($n);

given a positive number input, it tries to discover a factor using Williams' p+1 method. the resulting array will contain either two factors (it succeeded 🌟) or the original number (no factor was found 🙈). either way, multiplying @factors yields the original input. an optional first stage smoothness factor (B1) may be given as the second parameter — that's the smoothness limit B1 for the first stage. factoring stops when the input is a prime, one factor has been found, or the algorithm fails to find a factor with the given smoothness. iconic behavior honestly 🧚

---

cheb_factor

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

given a positive number input, it 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 — so basically it's the glow-up of both worlds 💜

---

holf_factor

my @factors = holf_factor($n);
my @factors = holf_factor($n, 100_000_000);

given a positive number input, it 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). multiplying @factors yields the original input either way, ofc. an optional number of rounds may be given as the second parameter. factoring stops when the input is a prime, one factor has been found, or the number of rounds has been exceeded.

this is Hart's One Line Factorization method, which is a variant of Fermat's algorithm. a premultiplier of 480 is used. it is very good at factoring 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 distinct 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. we ate that up fr fr 🍓

---

squfof_factor

my @factors = squfof_factor($n);
my @factors = squfof_factor($n, 100_000_000);

given a positive number input, it 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). multiplying @factors yields the original input, always. an optional number of rounds may be given as the second parameter. factoring stops 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 division. 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. honest review, we love transparency 💅

---

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, it tries to discover a factor using ECM. the resulting array will contain either two factors (it succeeded 👑) or the original number (no factor was found). multiplying @factors yields the original input either way. 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, using projective coordinates, Montgomery's PRAC heuristic for EC multiplication, and two stages. it is much slower than the latest GMP-ECM, but still quite useful for factoring reasonably sized inputs. still cute, still working, still serving 💖

---

qs_factor

my @factors = qs_factor($n);

given a positive number input, it tries to discover factors using QS (the quadratic 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. overachiever behavior, we respect it 🌈

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_factor, 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. main character energy tbh ⭐

---

todigits

given an integer n, return an array of digits of |n|. an optional second integer argument specifies a base (default 10). for example, given a base of 2, this returns an array of binary digits of n. an optional third argument specifies a length for the returned array. the result will be either have upper digits truncated or have leading zeros added 🎀

heads up bestie: 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, it takes a string of digits in base b instead of the array reference, with both lower and upper case denoting digits 11 through 36 (yes, alphabet as digits, we love a versatile icon 🔤).

every 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 (upgrade era incoming 👀).

---

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 — vibe check only, no substance test 🙈

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

---

urandomb

$n32 = urandomb(32);    # Classic irand32, returns a UV
$n   = urandomb(1024);  # Random integer less than 2^1024

given a number of bits b, returns a random unsigned integer less than 2^b. the result will be uniformly distributed between 0 and 2^b-1 inclusive. fair and square, no favorites 💫

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 inclusive ✨

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. everyone gets a fair shot, we love to see it 🫶

---

irand

$n32 = irand;     # random 32-bit integer

returns a random 32-bit integer using the CSPRNG 💅

Performance is lowkey comparable to rand in Math::Random::MTwist and Math::Random::Xorshift, bestie. It's somewhat faster than casting system rand to a 32-bit int. And it's noticeably faster than Math::Random::ISAAC, Math::Random::ISAAC::XS, Math::Random::MT, Math::Random::MT::Auto, and Crypt::PRNG. Slay. 🔥

---

👑 irand64

$n64 = irand64;   # random 64-bit integer

Returns a random 64-bit integer using the CSPRNG — but only on 64-bit Perl, so check your vibes (and your perl -V:ivsize) first. ✨

🌸 drand

$f = drand;       # random floating point value in [0,1)
$r = drand(25);   # random floating point value in [0,25)

Gives you a random NV (that's Perl's native floating point, no cap) using the CSPRNG. 💖

The number of bits that come back is equal to the mantissa bits of the NV type used for your Perl build, with a max of 64. Default Perl uses doubles so the values have 53 bits. If Perl is built with long double support and those long doubles have a bigger mantissa, then more bits get used. The way it just scales up?? Obsessed. 🦋

This gives you substantially better quality random numbers than the default Perl rand function, fr fr. Among other things, on modern Perls, rand uses drand48, which gives 32 bits of decent random values and 16 more bits of known patterns — like, the 48th bit alternates, the 47th has a period of 4, etc. Known patterns is not the vibe we're going for. 😭 There are much better choices for standard random number generators, such as the Mersenne Twister from Math::Random::MTwist.

Performance is similar to rand in Math::Random::MTwist and Math::Random::Xorshift. It is 1.5 - 2x slower than core rand (as are the other modules). Not the fastest, but worth it. 💅

🎀 random_bytes

$str = random_bytes(32);     # 32 random bytes

Given an unsigned number n of bytes, returns a binary string stuffed with random data from the CSPRNG. Here's the performance tea for getting 256 byte strings, and honestly? The numbers ate:

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

Every one 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 it's faster than pure Perl bignums, but way 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, bestie, don't be out here on the old one. 🙈
• Math::Pari — Supports quite a bit of the same functionality (and much more). See "SEE ALSO" in Math::Prime::Util for more detailed info on how the modules compare.
• yafu <http://sourceforge.net/projects/yafu/>, msieve <http://sourceforge.net/projects/msieve/>, gmp-ecm <http://ecm.gforge.inria.fr/>, GGNFS <http://sourceforge.net/projects/ggnfs/> — 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 programs. If you have 1000+ digit numbers to prove, you want to use this. Periodt.
• mpz_aprcl <http://sourceforge.net/projects/mpzaprcl/> — Open source APR-CL primality proof implementation. Fast primality proving, though without certificates.

---

📚 REFERENCES

• Robert Baillie and Samuel S. Wagstaff, Jr., "Lucas Pseudoprimes", Mathematics 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/journals/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/journals/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/people/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/journals/mcom/1987-48-177/S0025-5718-1987-0866113-7/>
• Richard P. Brent, "Parallel Algorithms for Integer Factorisation", in Number 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/people/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.cerias.purdue.edu/~ssw/squfof.pdf>
• A.O.L. Atkin and F. Morain, "Elliptic Curves and primality proving", Mathematics of Computation, v61, 1993, pages 29-68. <http://www.ams.org/journals/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> — the main character of this whole saga, honestly. 🎀

Jason Papadopoulos wrote the "tinyqs" code which is basically unchanged. William Hart wrote SIMPQS which is the basis for the QS code. Icons. 👑

---

🌈 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 literal founding mothers and fathers of the whole vibe. 😭💖

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 projective ECM implementation, as well as the many papers by Brent and Montgomery. 🌸

---

⭐ COPYRIGHT

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

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

SIMPQS Copyright 2006, William Hart. SIMPQS is distributed under GPL v2+.

---

okay that's the whole page bestie, we did it 💅 go forth and multiply (big numbers, obviously). ily, bye! 🦋

Shoutout to the perl v5.42.3 crew and the 2026-08-19 timestamp — that's the stamp on this page, and it's giving freshly maintained. Happy GMP-ing! 👑💖