Math::...e::Util(3pm)

💖 NAME 💖

Math::Prime::Util — utilities related to prime numbers, including fast sieves and factoring. 🧚✨

This is the 3pm manual for that whole vibe, bestie. We're in the group chat now.

---

🌸 VERSION 🌸

Version 0.74

---

✨ SYNOPSIS ✨

okay so nothing is exported by default, periodt. You either list the functions you want or you grab :all. The code below is the entire synopsis, kept exactly as-is because code is sacred and I do not play with it. 💅

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  # Free any memory used by the module.
  prime_memfree;

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

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

---

🌈 DESCRIPTION 🌈

A module for number theory in Perl. 💖 This includes prime sieving, primality tests, primality proofs, integer factoring, counts / bounds / approximations for primes, nth primes, and twin primes, random prime generation, and much more. It's giving "math girlie who actually knows what she's doing." 👑

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

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

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

Two scripts are also included and installed by default:

• primes.pl — displays primes between start and end values or expressions, with many options for filtering (e.g. twin, safe, circular, good, lucky, etc.). Use --help to see all the options. ✨
• factor.pl — operates similarly to the GNU factor program. It supports bigint and expression inputs. 🔥

---

🦋 ENVIRONMENT VARIABLES 🦋

There are two environment variables that affect operation. 👀 These are typically used for validation of the different methods or to simulate systems that have different support. (Yes bestie, keep reading, the vibes continue below.)

🛡️ MPU_NO_XS

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

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

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

🌊 MPU_NO_GMP

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

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

🎀 MPU_DEVNAMES

If set to 1, the PP package will be loaded on startup rather than on demand, and the package aliases "MPU", "PP", "GMP" will be used for the main, Perl, and GMP packages respectively. Normally you wouldn't want this both for aggressive namespace pollution and for performance (there is often no need to load the huge PP module). But it is convenient if one wants to call the different paths explicitly.

Regarding performance, on a 2020 Macbook M1, normal startup time is about 10 milliseconds. With this option set it becomes 45 milliseconds. This is the reason the PP code is only loaded if needed. For many purposes this amount of time is trivial, but slower computers or more time critical short applications will care.

---

✨ BIGNUM SUPPORT ✨

Okay bestie, so by default all functions support bigints. For performance, you should install Math::Prime::Util::GMP, which will be automatically used as a backend. It's giving speed. 💅

The default bigint class is Math::BigInt, which is not particularly speedy but is available by default in all Perl distributions, and is well tested. If you want to try something different, you can install and use Math::GMPz or Math::GMP which will be much faster. You can have this module use and return them using, for example:

prime_set_config(bigint => Math::GMPz);
my $n = next_prime(~0);
say "$n ",ref($n);
# 18446744073709551629 Math::GMPz

If you use Math::BigInt, I highly recommend also installing one of Math::BigInt::GMPz, Math::BigInt::GMP, or Math::BigInt::LTM. 🔥

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

• Install a recent version of Math::Prime::Util::GMP, as that will vastly increase the speed of many of the functions. This does require the GMP (http://gmplib.org) library be installed on your system, but this increasingly comes pre-installed or easily available using the OS vendor package installation tool. 🦋
• Install and use Math::BigInt::GMP (or "GMPz" or "LTM"), then use use bigint try => 'GMP,GMPz,LTM,Pari' in your script, or on the command line e.g. -Mbigint=lib,GMP. Large modular exponentiation is much faster using the better backends, as are the math and approximation functions when called with very large inputs. 💖
• I have run these functions on many versions of Perl, and my experience is that if you're using anything older than Perl 5.14, I would recommend you upgrade if you are using bignums a lot. There are some brittle behaviors on 5.12.4 and earlier with bignums. For example, the default BigInt backend in older versions of Perl will sometimes convert small results to doubles, resulting in corrupted output. 😭

---

🌈 PRIMALITY TESTING 🌈

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

Following the semantics used by Pari/GP, all primality test functions allow a negative primary argument, but will return false. All inputs must be integers or an error is raised. No cap. 🫶

---

💖 FUNCTIONS 💖

🌟 is_prime

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

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

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

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

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

✨ primes

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

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

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

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

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

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

💫 next_prime

$n = next_prime($n);

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

🦋 prev_prime

$n = prev_prime($n);

Returns the prime preceding the input number (i.e. the largest prime that is strictly less than the input). undef is returned if the input is 2 or lower. 🙈

The behavior in various programs of the previous prime function is varied. Pari/GP and Math::Pari returns the input if it is prime, as does nearest_le in Math::Prime::FastSieve. When given an input such that the return value will be the first prime less than 2, Math::Prime::FastSieve, Math::Pari, Pari/GP, and older versions of MPU will return 0. Math::Primality and the current MPU will return undef. WolframAlpha returns -2. Maple gives a range error. Drama. 💀

🌟 forprimes

forprimes { say } 100,200;                  # print primes from 100 to 200

$sum=0;  forprimes { $sum += $_ } 100000;   # sum primes to 100k

forprimes { say if is_prime($_+2) } 10000;  # print twin primes to 10k

Given a block and either an end count or a start and end pair, calls the block for each prime in the range. Compared to getting a big array of primes and iterating through it, this is more memory efficient and perhaps more convenient. This will almost always be the fastest way to loop over a range of primes. Nesting and use in threads are allowed. 💅

Math::BigInt objects may be used for the range.

For some uses an iterator (prime_iterator, prime_iterator_object) or a tied array (Math::Prime::Util::PrimeArray) may be more convenient. Objects can be passed to functions, and allow early loop exits. 🎀

💜 forcomposites

forcomposites { say } 1000;
forcomposites { say } 2000,2020;

Given a block and either an end number or a start and end pair, calls the block for each composite in the inclusive range. The composites, OEIS A002808 (http://oeis.org/A002808), are the numbers greater than 1 which are not prime: 4, 6, 8, 9, 10, 12, 14, 15, .... 🍓

🌸 foroddcomposites

Similar to forcomposites, but skipping all even numbers. The odd composites, OEIS A071904 (http://oeis.org/A071904), are the numbers greater than 1 which are not prime and not divisible by two: 9, 15, 21, 25, 27, 33, 35, .... ✨

🎀 forsemiprimes

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

This is essentially equivalent to:

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

🦋 foralmostprimes

foralmostprimes  { say }  3,  1000,2000;  # 3-almost-primes in [1000,2000]

Similar to forprimes, forsemiprimes, etc. but takes an additional first argument k and loops through the inclusive range for only those numbers with exactly k factors. If k=1 these are the primes, if k=2 these are the semiprimes, if k=3 these are the integers in the range with exactly 3 prime factors, etc. 👑

This is functionally equivalent to:

for ($a .. $b) { if (is_almost_prime($k,$_)) { ... } }
# or
for ($a .. $b) { if (prime_bigomega($_) == $k) { ... } }

though significantly faster and avoids issues with large loop variables. 🔥

💫 forfactored

forfactored { say "$_: @_"; } 100;

Given a block and either an end number or start/end pair, calls the block for each number in the inclusive range. $_ is set to the number while @_ holds the factors. Especially for small inputs or large ranges, this can be faster than calling factor on each sequential value. 🌈

Similar to the arrays returned by similar functions such as forpart, the values in @_ are read-only. Any attempt to modify them will result in undefined behavior. ⚠️

This corresponds to the Pari/GP 2.10 forfactored function.

✨ forsquarefree

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

This corresponds to the Pari/GP 2.10 forsquarefree function.

💖 forsquarefreeint

Similar to forsquarefree, but only sieves for square-free integers in the range (in segments so very large ranges still use little memory). No factoring information is returned: the @_ variable is not set. In return it is 2 to 20 times faster. 💅

As with range functions such as foralmostprimes this can be much faster than calling is_square_free for each integer in a large range. 🎀

🌟 fordivisors

fordivisors { $prod *= $_ } $n;

Given a block and a non-negative number n, the block is called with $_ set to each divisor in sorted order. Also see divisor_sum. 💜

🦋 forpart

forpart { say "@_" } 25;           # unrestricted partitions
forpart { say "@_" } 25,{n=>5}     # ... with exactly 5 values
forpart { say "@_" } 25,{nmax=>5}  # ... with <=5 values

Given a non-negative number n, the block is called with @_ set to the array of additive integer partitions. The operation is very similar to the forpart function in Pari/GP 2.6.x, though the ordering is different. The ordering is lexicographic. Use partitions to get just the count of unrestricted partitions. 🌸

An optional hash reference may be given to produce restricted partitions. Each value must be a non-negative integer. The allowable keys are:

• n — lock it to exactly this many values, no more no less 💅
• amin — every single element has to be at least this value
• amax — every single element has to be at most this value
• nmin — the array is required to have at least this many values
• nmax — the array is required to have at most this many values
• prime — every element must be prime (non-zero) or non-prime (zero) ✨

Okay so just like forcomb and forperm, the partition return values are read-only. Try to modify them and you're in undefined behavior territory, bestie, and not the cute kind 😭

---

forcomp

It's giving forpart's sibling, except it iterates over integer compositions instead of partitions. Think of it as every ordering of partitions — or flip it around, partitions are basically an ordered subset of compositions. The ordering is lexicographic. And good news: all the options from forpart work here too 🎀

Lowkey a fun fact: the number of unrestricted compositions of n is 2^(n-1).

forcomb

Give it non-negative arguments n and k, and the block gets called with @_ set to the k element array of values from 0 to n-1, representing all the combinations in lexicographical order. The binomial function tells you the total count, but this function is the one that actually enumerates the choices for you, periodt.

Instead of handing it a data array, you feed it an integer for n. A super convenient way to map back to array elements is:

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

Here the block maps the combination array @_ to array values, the n argument is the array itself (it gets evaluated as a scalar, so it just coughs up the size, iykyk), and the k argument is how big you want your combinations to be 💖

And same as forpart and forperm: the index return values are read-only. Modifying them = undefined behavior, don't do it 🙈

If you don't supply that second argument k, then all k-subsets come back, starting with the smallest set k=0 and rolling all the way to k=n. Every k-subset is in lexicographical order. That, bestie, is the power set of n 👑

This corresponds to the Pari/GP 2.10 forsubset function.

forperm

Give it non-negative argument n, and the block is called with @_ set to the k element array of values from 0 to n-1, representing permutations in lexicographical order. Total number of calls? n!. Fr fr. 🌈

Instead of a data array, it takes an integer for n. A convenient way to map to array elements:

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

where the block maps the permutation array @_ to array values, and the argument for n is given the array, because it gets evaluated as a scalar and hence gives you the size ✨

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

forderange

Similar to forperm, except it iterates over derangements — that's the set of permutations skipping any that maps an element to its original position. No element gets to stay home, everyone's moving 🦋

formultiperm

# Show all anagrams of 'serpent':
formultiperm { say join("",@_) } [split(//,"serpent")];

Similar to forperm but it takes an array reference as an argument. That gets treated as a multiset, and the block will be called with each multiset permutation. So while the standard permutation iterator takes a scalar and hands you index permutations, this one takes the set itself. Iconic behavior honestly.

If all values are unique, the results will be the same as a standard permutation. Otherwise the results will be like a standard permutation but with duplicate entries removed. Yes, you could just generate all permutations and filter out the dupes, but bestie that is very slow for large sets. This iterator is way more efficient 🔥

There's no ordering requirement for the input array reference. The results will be in lexicographic order.

forsetproduct

forsetproduct { say "@_" } [1,2,3],[qw/a b c/],[qw/@ $ !/];

Takes zero or more array references as arguments and iterates over the set product (a.k.a. Cartesian product or cross product) of the lists. The given subroutine is repeatedly called with @_ set to the current list. Since no de-duplication is done, this is not literally a "set" product — it's giving set product energy tho 💅

Zero or one array references is valid, but ngl the result is not very interesting. And if any array reference is empty, the product is empty, so no subroutine calls happen at all. Nothing. Zero. 🫥

The subroutine gets an array whose values are aliased to the inputs, and they are not set to read-only. So modifying the array inside the subroutine will cause side-effects. You've been warned, gorgeous 💖

As with other iterators, the lastfor function will cause an early exit.

lastfor

forprimes { lastfor,return if $_ > 1000; $sum += $_; } 1e9;

Calling lastfor requests that the current for... loop stop after this call. Ideally it would act exactly like a last inside a loop, but technical reasons mean it does not exit the block early — hence you typically add a return if needed. Slight inconvenience, massive payoff, we move 😌

prime_iterator

my $it = prime_iterator;
$sum += $it->() for 1..100000;

Returns a closure-style iterator. The start value defaults to the first prime (2), but you can hand it an initial value as an argument, and then the first value returned will be the next prime greater than or equal to that argument. For example, this:

my $it = prime_iterator(200);  say $it->();  say $it->();

will return 211 followed by 223, because those are the next primes >= 200. On each call, the iterator returns the current value and increments to the next prime. Chef's kiss ✨

Other options include forprimes (more efficiency, less flexibility), Math::Prime::Util::PrimeIterator (an iterator with more functionality), or Math::Prime::Util::PrimeArray (a tied array). Pick your fighter 🎀

prime_iterator_object

my $it = prime_iterator_object;
while ($it->value < 100) { say $it->value; $it->next; }
$sum += $it->iterate for 1..100000;

Returns a Math::Prime::Util::PrimeIterator object. It's a shortcut that loads the package if needed, calls new, and returns the object. See the documentation for that package for the details. This object has more features than the simple one above (e.g. the iterator is bi-directional), and it also handles iterating across bigints. Main character object, truly 👑

prime_count

my $primepi = prime_count( 1_000 );
my $pirange = prime_count( 1_000, 10_000 );

Returns the Prime Count function Pi(n), also called primepi in some math packages. When you give it two arguments, it returns the inclusive count of primes between the ranges. E.g. (13,17) returns 2, (14,17) and (13,16) return 1, and (14,16) returns 0. Math is mathing 🌸

The current implementation decides based on the ranges whether to use a segmented sieve with a fast bit count, or the extended LMO algorithm. The former is preferred for small sizes as well as small ranges. The latter is much faster for large ranges.

The segmented sieve is very memory efficient and is quite fast even with large base values. Its complexity is approximately O(sqrt(a) + (b-a)), where the first term is typically negligible below ~ 10^11. Memory use is proportional only to sqrt(a), with total memory use under 1MB for any base under 10^14.

The extended LMO method has complexity approximately O(b^(2/3)) + O(a^(2/3)), and also uses low memory. A calculation of Pi(10^14) completes in a few seconds, Pi(10^15) in well under a minute, and Pi(10^16) in about one minute. In contrast, even parallel primesieve would take over a week on a similar machine to determine Pi(10^16). Yeah. Read that again. 💀

Also see the function prime_count_approx which gives a very good approximation to the prime count, and prime_count_lower and prime_count_upper which give tight bounds to the actual prime count. These functions return quickly for any input, including bigints 🔥

prime_count_upper

Returns a proven upper bound on the number of primes up to n. See prime_count_lower for details common to both functions.

prime_count_lower

my $lower_limit = prime_count_lower($n);
my $upper_limit = prime_count_upper($n);
#   $lower_limit  <=  prime_count(n)  <=  $upper_limit

Returns a proven lower bound on the number of primes up to n. These are analytical routines, so they take a fixed amount of time and no memory. The actual prime_count will always be equal to or between these numbers. Sandwich them, bestie 🥪

A common place these would be used is sizing an array to hold the first $n primes. It may be desirable to use a bit more memory than is necessary, to avoid calling prime_count.

These routines use verified tight limits below a range at least 2^35. For larger inputs various methods are used including Dusart (2010), Büthe (2014,2015), and Axler (2014). These bounds do not assume the Riemann Hypothesis. If the configuration option assume_rh has been set (it is off by default), then the Schoenfeld (1976) bounds can be used for very large values. Rizz of the mathematicians, honestly 💜

prime_count_approx

print "there are about ",
      prime_count_approx( 10 ** 18 ),
      " primes below one quintillion.\n";

Returns an approximation to the prime_count function, without having to generate any primes. For values under 10^36 this uses the Riemann R function, which is quite accurate: an error of less than 0.0005% is typical for input values over 2^32, and it decreases as the input gets larger. She ate that ✨

A slightly faster but much less accurate answer can be obtained by averaging the upper and lower bounds.

is_prime_power

Given an integer n, returns k if n = p^k for some prime p, and zero otherwise.

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

This corresponds to Pari/GP's isprimepower function. It is related to Mathematica's PrimePowerQ[n] function. These all return zero/false for n=1.

This is the OEIS series A246655 http://oeis.org/A246655 🦋

prime_powers

my $aref = prime_powers( 10**4 );

Given either two non-negative limits lo, hi, or one non-negative limit hi, returns an array reference with all prime powers between the limits (inclusive). With only one input, the lower limit is 2.

The array reference values will be all p^e where lo <= p^e <= hi with p prime and e >= 1. Hence this includes the primes as well as the higher powers of primes.

See also primes and prime_power_count.

next_prime_power

Given an integer n, returns the smallest prime power greater than |n|. Similar to next_prime, but also includes powers of primes. Glow up edition ✨

prev_prime_power

Given an integer n, returns the greatest prime power less than |n|. Similar to prev_prime, but also includes powers of primes. If given |n| less than 3, undef will be returned. No prime power for you 💀

prime_power_count

Given a single non-negative integer n, returns the count of prime powers less than or equal to n. If given two non-negative integers lo and hi, returns the count of prime powers between lo and hi inclusive. Counting era 🎀

These are prime powers where the exponent is greater than 0 — so like, prime powers but not counting 1, bestie. It's giving OEIS series A025528 <http://oeis.org/A025528>.

💫 prime_power_count_approx

You hand it a non-negative integer n and it quickly hands back a good estimate of how many prime powers are less than or equal to n. Speed queen behavior fr fr ✨

💫 prime_power_count_lower

Non-negative integer n goes in, a lower bound comes out — fast. The real count will always be greater than or equal to whatever it gives you. So it's a floor, not the whole truth 💅

💫 prime_power_count_upper

Same deal with the non-negative n, but now it's an upper bound, returned quickly. The real count will always be less than or equal to the result. Ceiling core 🎀

💫 nth_prime_power

Give it a non-negative integer n and you get the n-th prime power. That's it, that's the function, no notes 🦋

💫 nth_prime_power_approx

Non-negative n in, a good estimate of the n-th prime power out, and it's fast. Approximate but make it snappy 💖

💫 nth_prime_power_lower

Non-negative n → a lower bound on the n-th prime power, quick quick. The actual value will always be greater than or equal to the result.

💫 nth_prime_power_upper

Non-negative n → an upper bound on the n-th prime power, also quick. The actual value will always be less than or equal to the result. No cap 🌈

💫 twin_primes

Returns the lesser of the twin primes living between the lower and upper limits, inclusive — and if you don't give a lower limit, it lowkey defaults to 2. Iconic OEIS series A001359 <http://oeis.org/A001359> 💞

So picture a twin prime pair (p,q) where q = p + 2, p is prime and q is prime. This function represents the whole pair using just p. Which means: your bounds need to include p, and the list you get back will have p but not q. Sorry to q, she's not invited 🙈

It works exactly like the primes function, except only the first primes of twin prime pairs come back. And just like that function, an array reference is returned, not a plain array.

💫 twin_prime_count

Basically prime count's twin sister (pun intended 💀), except it counts twin primes — primes p where p+2 is also prime. You can call it with either:

• one single number → a count from 2 up to that argument, or
• two numbers → a range.

The primes being counted are the first value, so a range of (3,5) returns a count of two, because both 3 and 5 count as twin primes. A range of (12,13) returns a count of zero, because neither 12+2 nor 13+2 are prime. In contrast, primesieve requires all elements of a constellation to be inside the range to be counted, so it would return one for that first example (5 doesn't get counted because its pair 7 isn't in the range). Main character logic differences, iykyk ✨

There's no useful formula known for this, unlike prime counts. So we sieve for the answer, using some small table acceleration to make it less painful.

💫 twin_prime_count_approx

Returns an approximation to the twin prime count of n. It comes back fast and has a very small error for large values 🫶

The method used is conjecture B of Hardy and Littlewood 1922, as stated in Sebah and Gourdon 2002. For inputs under 10M, a correction factor is additionally applied to reduce the mean squared error. Details, we have them 💅

💫 semi_primes

Returns an array reference to semiprimes between the lower and upper limits, inclusive, with a lower limit of 4 if none is given. It's OEIS A001358 <http://oeis.org/A001358>. Semiprimes are composite integers that are products of exactly two primes — no more, no less 🎀

Works just like the primes function, and like that one, you get back an array reference.

💫 semiprime_count

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

There's a fast method that only requires computation up to the square root of the range end — unless the range is so small that just walking it is faster. Efficiency slay 🔥

💫 semiprime_count_approx

Returns an approximation to the semiprime count of n. Comes back quickly and is square root accurate for native size inputs.

The series of Crisan and Erban (2020) is used with a maximum of 19 terms. Truncation is performed at empirical good crossovers. Clamping is done as needed at crossovers to ensure monotonic results. Very "math girlies stay winning" of it 💜

💫 almost_primes

my $ref_to_3_almost_primes = almost_primes(3, 1000, 2000);

Takes a non-negative integer argument k, plus either one or two more non-negative integer arguments for the upper limit, or the lower and upper limits. Limits are inclusive. The k-almost-primes are the integers with exactly k prime factors.

Works just like the primes function, and again — array reference, babes.

The lore on k:

• k=1 → these are the primes (OEIS A000040 <http://oeis.org/A000040>).
• k=2 → these are the semiprimes (OEIS A001358 <http://oeis.org/A001358>).
• k=3 → these are the 3-almost-primes (OEIS A014612 <http://oeis.org/A014612>).
• k=4 → these are the 4-almost-primes (OEIS A014613 <http://oeis.org/A014613>).

OEIS sequences can be found through k=20. The pipeline is deep 🌸

💫 almost_prime_count

say almost_prime_count(3,10000); # number of 3-almost-primes <= 10000

Given non-negative integers k and n, returns the count of k-almost-prime numbers up to and including n. With k=1 this is the standard prime count. With k=2 this is the semiprime count. In general, it's the count of all integers through n that have exactly k prime factors.

Under the hood it uses nested prime count sums, and caching along with LMO prime counts to get quite reasonable speeds. Nested loops but make it cute ✨

💫 almost_prime_count_approx

A fast approximation of the k-almost-prime count of n.

Heads up bestie: the current implementation for n greater than 64-bit is not well tested. Proceed with caution 💀

💫 almost_prime_count_lower

Quickly returns a lower bound for the k-almost-prime count of n. The actual count will be greater than or equal to this result.

Same warning, no skipping: the current implementation for n greater than 64-bit is not well tested.

💫 almost_prime_count_upper

Quickly returns an upper bound for the k-almost-prime count of n. The actual count will be less than or equal to this result.

And again 😭 the current implementation for n greater than 64-bit is not well tested.

💫 omega_primes

Takes a non-negative integer argument k and either one or two additional non-negative integer arguments indicating the upper limit or lower and upper limits. Limits are inclusive. The k-omega-primes are positive integers with exactly k distinct prime factors, with possible multiplicity. So these numbers are divisible by exactly k different primes. Distinct, not total, ok? Read the label 🏷️💅

The k-omega-primes (not a common term, she's niche 🦋) are exactly those integers where prime_omega(n) == k. Compare that to k-almost-primes where prime_bigomega(n) == k.

• k=1 → these are the prime powers.
• k=2 → these are OEIS A007774 <http://oeis.org/A007774>.
• k=3 → these are OEIS A033992 <http://oeis.org/A033992>.

💫 omega_prime_count

Given non-negative integers k and n, returns the count of k-omega-prime numbers from 1 up to and including n. This is the count of all positive integers through n that are divisible by exactly k different primes.

The implementation uses nested loops over prime powers. Loops, but make them elegant ✨

Though we have defined prime_omega(0) = 1, it is not included. We see it but we don't count it 🙈

💫 ramanujan_primes

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

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

Real talk: generating Ramanujan primes takes some effort, including overhead to cover a range. This will be substantially slower than generating standard primes. Budget accordingly bestie 😭

💫 ramanujan_prime_count

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

While it's not nearly as efficient as prime_count, it does use a number of speedups that result in it being much more efficient than generating all the Ramanujan primes. Glow up arc 🔥

💫 ramanujan_prime_count_approx

A fast approximation of the count of Ramanujan primes under n. Quick, cute, approximate 💖

💫 ramanujan_prime_count_lower

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

💫 ramanujan_prime_count_upper

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

💫 sieve_range

my @candidates = sieve_range(2**1000, 10000, 40000);

Given a start value n, plus 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 smaller than itself and less than or equal to depth. So a depth of 2 will remove all even numbers (other than 2 itself if it's in the range). A depth of 3 will remove all numbers divisible by 2 or 3, other than those primes themselves. Filtered, filtered, slay 🧚

💫 sieve_prime_cluster

my @s = sieve_prime_cluster(1, 1e9, 2,6,8,12,18,20);

Efficiently finds prime clusters between the first two arguments low and high. The remaining arguments describe the cluster. Cluster values must be even, less than 31 bits, and strictly increasing — that's a hard requirement, not a suggestion 💅

Given a cluster set C, the returned values are all primes in the range where p+c is prime for each c in the cluster set C. For returned values under 2^64, all cluster values are definitely prime. Above this range, all cluster values are BPSW probable primes (no counterexamples known). So like, probably prime but we're not swearing on it up there 👑

This function returns an array rather than an array reference. Typically the number of returned values is much lower than for other primes functions, so this uses the more convenient array return. It has an identical signature to the function of the same name in Math::Prime::Util:GMP. Matching energy 👯

The cluster is described as offsets from 0, with the implicit prime at 0. Which means:

• an empty list is asking for all primes (the cluster p+0),
• a list with the single value 2 will find all twin primes (the cluster where p+0 and p+2 are prime),
• the list 2,6,8 will find prime quadruplets.

And note: there is no requirement that the list denote a constellation (a cluster with minimal distance) — the list 42,92,606 is just fine. Freeform, we love to see it 🎀

💫 sum_primes

Returns the summation of primes between the lower and upper limits, inclusive, with a lower limit of 2 if none is given. It's essentially similar to either of:

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

but it is much more efficient. Same answer, way less work, that's the vibe 💜

okay so real talk about the engine under the hood 💅 — the current implementation is a mix of small-table-enhanced sieve count for sums that fit in a UV, an efficient sieve count for small ranges, and a Legendre sum method, including XS support for 128-bit results.

And like, it's pretty efficient, no cap! BUT. The actual state of the art is Kim Walisch's primesum https://github.com/kimwalisch/primesum. That one is recommended for very large values, because it can be hundreds of times faster. Hundreds, bestie. 🫠

print_primes

print_primes(1_000_000);             # print the first 1 million primes
print_primes(1000, 2000);            # print primes in range
print_primes(2,1000,fileno(STDERR))  # print to a different descriptor

Give it one argument and it prints every prime from 2 up to n straight to standard out. Give it two arguments and it prints the primes between low and high to standard output. Give it three and it prints primes between low and high to the file descriptor you hand it. And if that file descriptor can't be written to? It deadass croaks with print_primes write error. 😭

The output it makes is literally identical to this:

forprimes { say } $low,$high;

The whole point of this function is pure efficiency. It's over 10x faster than using say, print, or printf, though much more limited in functionality. A later version may allow a file handle as the third argument. (manifesting 🙏)

nth_prime

say "The ten thousandth prime is ", nth_prime(10_000);

Returns the prime sitting at index n in the array of prime numbers. Or, put in the fancy way: it returns the smallest p such that Pi(p) >= n. ✨

Like basically every other program with similar functionality, this is one-based. So nth_prime(0) hands you back undef, and nth_prime(1) gives you 2.

For relatively small inputs (below 1 million or so), it does a sieve over a range containing the nth prime, then counts up to the number — that's pretty efficient in both time and memory. For larger values, it creates a low-biased estimate using the inverse logarithmic integral, uses a fast prime count, then sieves in the small difference.

And yeah, this method is thousands of times faster than generating primes and doesn't need any big tables of precomputed values, but it still can take a fair amount of time for large inputs, fr. Calculating the 10^12th prime takes about 1 second, the 10^13th prime takes under 10 seconds, and the 10^14th prime (3475385758524527) takes under 30 seconds. So think about whether a bound or an approximation would be acceptable, because those can be computed analytically. 🧠

If your result is larger than a native integer size (32-bit or 64-bit), it's going to take a very long time. A later version of Math::Prime::Util::GMP may include this functionality, which would actually help on 32-bit machines.

nth_prime_upper

Returns a proven upper bound on the Nth prime. Peep "nth_prime_lower" for the details that both functions share.

nth_prime_lower

my $lower_limit = nth_prime_lower($n);
my $upper_limit = nth_prime_upper($n);
# For all $n:   $lower_limit  <=  nth_prime($n)  <=  $upper_limit

Returns a proven lower bound on the Nth prime. No sieving happens here at all, so these are fast even for big ol' inputs. ⚡

For tiny values of n, exact answers are returned. For small inputs, an inverse of the opposite prime count bound is used. And for larger values, the Dusart (2010) and Axler (2013) bounds do the heavy lifting.

nth_prime_approx

say "The one trillionth prime is ~ ", nth_prime_approx(10**12);

Returns an approximation of the nth_prime function, without generating a single prime. Gorgeous. 🎀 For values where the nth prime is smaller than 2^64, the inverse Riemann R function is used. For larger values, the inverse logarithmic integral is used.

The value it returns will not necessarily be prime, bestie. And that goes for all the nth prime approximations below too — the returned value is close to the real value, but nobody's making any effort to coerce the result to the nearest set element.

nth_twin_prime

Returns the Nth twin prime. It's done via sieving and counting, so it's not very fast for large values. Lowkey a grind. 💀

nth_twin_prime_approx

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

nth_semiprime

Returns the Nth semiprime — basically where a forsemiprimes loop would end up after N iterations, except way more efficiently. 💅

nth_semiprime_approx

Returns an approximation to the Nth semiprime. This approximation is orders of magnitude better than the simple n log n / log log n approximation for large n. Like, for n=10^12 the simple estimate lands within 3.6%, but this function is within 0.000012%. Chew on that. 😳

nth_almost_prime

say "500th number with exactly 3 factors: ", nth_almost_prime(3,500);

A k-almost prime is a product of k prime numbers, counted with multiplicity. In other words: there are exactly k prime factors (and they do not have to be distinct from each other, we don't gatekeep like that).

Give it non-negative integers k and n, and it returns the n-th k-almost prime. With k=1 this is the nth prime. With k=2 this is the nth semiprime. See the pattern? 👑

The implementation does a binary search lookup with almost_prime_count, so it's reasonably efficient for large values.

undef is returned for n == 0 and for all k == 0 other than n == 1.

nth_almost_prime_approx

A fast approximation of the n-th k-almost prime. 🦋

nth_almost_prime_lower

Quickly returns a lower bound for the n-th k-almost prime. The actual nth k-almost-prime will be greater than or equal to this result.

nth_almost_prime_upper

Quickly returns an upper bound for the n-th k-almost prime. The actual nth k-almost-prime will be less than or equal to this result.

nth_omega_prime

Give it non-negative integers k and n, and it returns the n-th k-omega prime. That's the n-th integer divisible by exactly k different primes. (Different from the almost-primes, where repeats are allowed — clock that. ⏰)

The implementation does a binary search lookup with omega_prime_count, so it's reasonably efficient for large values.

undef is returned for n == 0 and for all k == 0 other than n == 1.

nth_ramanujan_prime

Returns the Nth Ramanujan prime. For reasonable size values of n, like under 10^8 or so, this is relatively efficient for single calls. But if you're making multiple calls, it'll be much more efficient to just get the list once. Save yourself the pain bestie. 🫶

nth_ramanujan_prime_approx

A fast approximation of the Nth Ramanujan prime.

nth_ramanujan_prime_lower

A fast lower limit on the Nth Ramanujan prime.

nth_ramanujan_prime_upper

A fast upper limit on the Nth Ramanujan prime.

is_pseudoprime

Give it an integer n and zero or more positive bases, and it returns 1 if n is positive and a probable prime to each base, and 0 otherwise. This is the simple Fermat primality test. Removing primes, given base 2 this produces the sequence OEIS A001567 http://oeis.org/A001567.

If you don't give it any bases, base 2 is used. All bases must be 2 or greater.

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

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

is_euler_pseudoprime

Give it an integer n and zero or more positive bases, and it returns 1 if n is positive and an Euler probable prime to each base, and 0 otherwise. This is the Euler test, sometimes called the Euler-Jacobi test. Removing primes, given base 2 this produces the sequence OEIS A047713 http://oeis.org/A047713.

If no bases are given, base 2 is used. All bases must be 2 or greater.

If 0 comes back, then the number really is a composite (for bases less than n). If 1 comes back, then it's either a prime or an Euler pseudoprime to all the given bases. Give it enough distinct bases and the chances become very high that the number is genuinely prime. ✨

This test forms the basis of the Solovay-Strassen test, which is a precursor to the Miller-Rabin test (which uses the strong probable prime test). And good news: there are no analogies to the Carmichael numbers for this test. For the Euler test, at most 1/2 of witnesses pass for a composite, while at most 1/4 pass for the strong pseudoprime test.

is_strong_pseudoprime

my $maybe_prime = is_strong_pseudoprime($n);
my $probably_prime = is_strong_pseudoprime($n, 2, 3, 5, 7, 11, 13, 17);

Give it an integer n and zero or more positive bases, and it returns 1 if n is positive and a strong probable prime to each base, and 0 otherwise.

If no bases are given, base 2 is used. All bases must be 2 or greater.

If 0 comes back, then the number really is a composite (for any base). If 1 comes back, then it's either a prime or a strong pseudoprime to all the given bases. Given enough distinct bases, the chances become very high that the number is actually prime. 🔥

This is usually used in combination with other tests to make either stronger tests (e.g. the strong BPSW test) or deterministic results for numbers less than some verified limit (e.g. it has long been known that no more than three selected bases are required to give correct primality test results for any 32-bit number). And given the small chances of passing multiple bases, there are some math packages that just use multiple MR tests for primality testing. Iconic shortcut tbh.

Even inputs other than 2 will always return 0 (composite). The algorithm does run with even input, but most sources define it only on odd input. Returning composite for all non-2 even input makes the function match most other implementations, including Math::Primality's is_strong_pseudoprime function.

Generally, the bases you actually care about are between 2 and n-2. Bases 1 and n-1 will return 1 for any odd composites. Most sources do not define the test for bases equal to 0 mod n, and many don't for any bases larger than n. We allow all bases, noting that the case base = 0 mod n is defined as 1. This allows primes to return 1 regardless of the base. 🎀

is_lucas_pseudoprime

Give it an integer n, and it returns 1 if n is positive and a standard Lucas probable prime using the Selfridge method of choosing D, P, and Q (some sources call this a Lucas-Selfridge pseudoprime). Removing primes, this produces the sequence OEIS A217120 http://oeis.org/A217120.

is_strong_lucas_pseudoprime

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

is_extra_strong_lucas_pseudoprime 💫

Hand it an integer n and it returns 1 when n is positive and counts as an extra strong Lucas probable prime, exactly as defined in Grantham 2000 (http://www.ams.org/mathscinet-getitem?mr=1680879). No cap, this test is stricter than the strong Lucas test — like, way more picky about who gets in — and it coughs up roughly 60% fewer pseudoprimes. Also? Performance is typically 20-30% faster than the strong Lucas test. Bestie it's faster AND stricter, that's iconic.

The parameters get picked via the Baillie-OEIS method (http://oeis.org/A217719): you increment P starting from 3 until jacobi(D,n) = -1. And if you strip out the primes, you're left staring at the sequence OEIS A217719 (http://oeis.org/A217719). Slay.

---

is_almost_extra_strong_lucas_pseudoprime 🌸

Basically the same vibe as is_extra_strong_lucas_pseudoprime, except it doesn't compute U, so it's a lil faster but also a lil weaker. Honestly? With the current implementations there's barely any reason to pick this one unless you're specifically trying to reproduce certain results. The extra-strong implementation got optimized to use similar features, which ate most of the performance advantage. So the gap is lowkey gone. 💀

There's an optional second argument (an integer between 1 and 256) that sets the increment amount for P parameter selection. The default value of 1 gives you the parameter selection described up in is_extra_strong_lucas_pseudoprime, which creates a pseudoprime sequence that's a superset of the latter's pseudoprime sequence OEIS A217719 (http://oeis.org/A217719). A value of 2 gives you the method Pari uses (http://pari.math.u-bordeaux.fr/faq.html#primetest).

Since the U = 0 condition gets ignored, this produces about 5% more pseudoprimes than the extra-strong Lucas test. But hold on — that's still only 66% of the number produced by the strong Lucas-Selfridge test. And for the record, no BPSW counterexamples have ever been found with any of the Lucas tests described here. Not one. 🫶

---

is_euler_plumb_pseudoprime 🎀

Give it an integer n, and it returns 1 if n is positive and passes Colin Plumb's Euler Criterion primality test, otherwise 0. 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. Set theory girlie behavior fr fr.

The whole reason this test exists? It's slightly more efficient than the other probable prime tests. That's it. That's the tea. ✨

---

is_perrin_pseudoprime 🔥

Given an integer n, returns 1 if n is positive and n divides P(n) where P(n) is the Perrin number of n, otherwise 0. 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 little recurrence, no notes.

Pseudoprimes are relatively rare here (the first two are 271441 and 904631), but infinitely many exist. And they overlap a lot with the base-2 pseudoprimes and strong pseudoprimes, which makes this test inferior to the Lucas or Frobenius tests when you're doing combined testing. The pseudoprime sequence is OEIS A013998 (http://oeis.org/A013998).

The implementation uses modular pre-filters, Montgomery math, and the Adams/Shanks doubling method. This is significantly more efficient than other known implementations. We ate. 👑

There's an optional second argument r that says whether to run extra 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_catalan_pseudoprime 🦋

Give it an integer n and it returns 1 if n is positive and (-1)^{(n-1)/2} * C_{(n-1)/2} is congruent to 2 mod n, where C_n is the nth Catalan number, otherwise 0. The nth Catalan number equals binomial(2n,n)/(n+1). All odd primes satisfy this condition, and only three known composites do. Three. 😭

The pseudoprime sequence is OEIS A163209 (http://oeis.org/A163209).

Real talk though: there's no known efficient method to perform the Catalan primality test, so this is a curiosity rather than a practical test. It's giving "fun fact" energy, not "use me in prod" energy. The implementation uses a method from Charles Greathouse IV (2015) and results from Aebi and Cairns (2008) to produce results many orders of magnitude faster than other known implementations — but it's still vastly slower than other compositeness tests. Like, majestic but slow. A sloth in a ballgown. 💅

---

is_frobenius_pseudoprime 💜

Give it an integer n plus two optional integer parameters a and b, and it returns 1 if n is positive and is a Frobenius probable prime with respect to the polynomial x^2 - ax + b, otherwise 0. Without the parameters, b = 2 and a is the least positive odd number such that (a^2-4b|n) = -1. That selection has no pseudoprimes below 2^64 and none known at all. In any case, the discriminant a^2-4b must not be a perfect square. Don't do it.

Some authors use the Fibonacci polynomial x^2-x-1 corresponding to (1,-1) as the default method for a Frobenius probable prime test. That creates a weaker test than most other parameter choices (e.g. over twenty times more pseudoprimes than (3,-5)), so it is not the default here. With the (1,-1) parameters, the pseudoprime sequence is OEIS A212424 (http://oeis.org/A212424).

The Frobenius test is stronger than the Lucas test. Any Frobenius (a,b) pseudoprime is also a Lucas (a,b) pseudoprime, but the converse is not true, as any Frobenius (a,b) pseudoprime is also a Fermat pseudoprime to the base |b|. With the default parameters you can see this lands similar to, but somewhat weaker than, the BPSW test this module uses (which uses the strong and extra-strong versions of the probable prime and Lucas tests respectively).

Also go peek at the more efficient is_frobenius_khashin_pseudoprime and is_frobenius_underwood_pseudoprime — those have no known counterexamples and run quite a bit faster. 👀

---

is_frobenius_underwood_pseudoprime ⭐

Given an integer n, returns 1 if n is positive and passes the efficient Frobenius test of Paul Underwood, otherwise 0. This picks 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 for a cost of only slightly more than 2 strong pseudoprime tests. Which makes it similar to, but faster than, a regular Frobenius test. Efficient queen. 🔥

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

---

is_frobenius_khashin_pseudoprime 🌈

Given an integer n, returns 1 if n is positive and passes the Frobenius test of Sergey Khashin, otherwise 0. 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. Clean, efficient, no notes.

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

---

miller_rabin_random 🙈

Give it an integer n and a positive integer k, and it returns 1 if n is positive and passes k Miller-Rabin tests using uniform random bases selected between 2 and n-2.

Please don't use this in place of is_prob_prime, is_prime, or is_provable_prime. Those functions will be faster and provide better results than running k Miller-Rabin tests. This one's here in case you want more assurances for non-proven primes — like cryptographic uses where the size is large enough that proven primes aren't desired.

---

is_prob_prime 🍓

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

Given an integer n, returns 0 (composite), 2 (definitely prime), or 1 (probably prime). Three whole answers, she's versatile.

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

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

---

is_bpsw_prime 💖

Given an integer n, returns 0 (composite), 2 (definitely prime), or 1 (probably prime), using the BPSW primality test (extra-strong variant). Normally one of the is_prime in Math::Prime::Util or is_prob_prime in Math::Prime::Util functions will suffice, but those functions do pre-tests to find easy composites. If you already know that's not necessary, then calling is_bpsw_prime may save you a small amount of time. Free speed, no cap.

---

is_provable_prime 👑

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

Given an integer n, returns 0 (composite), 2 (definitely prime), or 1 (probably prime). Same return values as is_prime and is_prob_prime, consistency besties. Note that numbers below 2^64 are considered proven by the deterministic set of Miller-Rabin bases or the BPSW test. Both of these have been tested against all the small (64-bit) composites and they do not return false positives.

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

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

Okay bestie, buckle up because this section is where Math::Prime::Util::GMP really shows off. ✨

Math::Prime::Util::GMP implements both the BLS75 theorem 5 test and ECPP (elliptic curve primality proving). 🎀 It'll typically try a quick "n-1" proof before it pulls out ECPP. You can get certificates with either method. The result? Proofs of 200-digit primes in under 1 second on average, and many hundreds of digits are genuinely possible. 💅 Which makes it significantly faster than Pari 2.1.7's "is_prime(n,1)", which is the default for Math::Pari. Not even a competition, fr fr.

---

💖 prime_certificate

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

Give it an integer "n" and it hands back a primality certificate as a multi-line string. 🦋 If we couldn't prove "n" prime, you get an empty string back (and heads up: "n" may or may not actually be composite — it's just unproven). You can eyeball that certificate yourself or hand it to "verify_prime" to check it. The format details live over in that function's description, so scroll down, bestie. 👇

---

💜 is_provable_prime_with_cert

Give it an integer "n" and it returns a two-element array: the result of "is_provable_prime" plus the certificate. The result values are:

0  definitely composite
1  probably prime
2  definitely prime

And the primality certificate comes along just like "prime_certificate" would give you. 🎀 If the first element is not 2, then the certificate is an empty string. Simple as that.

---

🌸 verify_prime

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

Hand it a primality certificate and it returns either 0 (nope, not verified) or 1 (verified, slay). 👑

Most of the computation is done in pure Perl with Math::BigInt, so you probably want to install and use Math::BigInt::GMP. And ECPP certificates will be faster with Math::Prime::Util::GMP doing its elliptic curve thing. 🧚

If the certificate is malformed, the routine will carp a warning in addition to returning 0. If the "verbose" option is set (peep "prime_set_config"), then when validation fails, the reason for the failure gets printed too, plus the 0. And if "verbose" is set to 2 or higher, you also get a message saying it succeeded and what type the certificate is. 💖

Now here's the fun part — the actual format. A certificate may have arbitrary text before the beginning (the primality routines from this module won't add any extra text, but this way verbose output from the prover can be safely stored inside a certificate). The certificate begins with this line:

[MPU - Primality Certificate]

Every line beginning with "#" is treated as a comment and ignored, and so are blank lines. A version number may follow, like:

Version 1.0

For all inputs, base 10 is the default, but at any point you can change it with a line like:

Base 16

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

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

Proof for:
N 100003

The text "Proof for:" is the signal that we're about to read an "N" value. Skipping comments and blank lines, the next line should be "N " followed by the number.

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

If Q is prime, then N is prime.

Some blocks have more than one Q value hanging off them, but most only have one. Each block has its own set of conditions that must be verified, and each one can be checked completely self-contained. That is — every block is independent of the others and they can be processed in any order. 🔥 But for the whole thing to count as a complete proof, every block has to verify successfully. The block types and their conditions are below, bestie. 🫶

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

Some other certificate formats (like Primo) use an ordered chain, where the first block must be for the initial "N", a single "Q" is given as the implied "N" for the next block, and so on. That simplifies validation a bit and strips some redundant info from the certificate — but it has no obvious way to add proof types like Lucas or the various BLS75 theorems that use multiple factors. So the choice here was: the most general solution is to have the certificate contain the set in any order, and let the verifier do the work of constructing the tree. 💅

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

💖 "Small"

Type Small
N 5791

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

💖 "BLS3"

Type BLS3
N  2297612322987260054928384863
Q  16501461106821092981
A  5

A simple n-1 style proof using BLS75 theorem 3. This block verifies if:

a  Q is odd
b  Q > 2
c  Q divides N-1
.  Let M = (N-1)/Q
d  MQ+1 = N
e  M > 0
f  2Q+1 > sqrt(N)
g  A^((N-1)/2) mod N = N-1
h  A^(M/2) mod N != N-1

💖 "Pocklington"

Type Pocklington
N  2297612322987260054928384863
Q  16501461106821092981
A  5

A simple n-1 style proof using generalized Pocklington. This one is more restrictive than BLS3 and way more restrictive than BLS5. It's Primo's type 1, and this module does not currently generate these blocks. This block verifies if:

a  Q divides N-1
.  Let M = (N-1)/Q
b  M > 0
c  M < Q
d  MQ+1 = N
e  A > 1
f  A^(N-1) mod N = 1
g  gcd(A^M - 1, N) = 1

💖 "BLS15"

Type BLS15
N  8087094497428743437627091507362881
Q  175806402118016161687545467551367
LP 1
LQ 22

A simple n+1 style proof using BLS75 theorem 15. This block verifies if:

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

💖 "BLS5"

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

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

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

💖 "ECPP"

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

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

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

---

🌈 is_aks_prime

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

Give it an integer "n" and you get 1 if "n" is positive and passes the Agrawal-Kayal-Saxena (AKS) primality test, 0 otherwise. This is a deterministic unconditional primality test that runs in polynomial time for general input. 👑

Now listen — while this is an important theoretical algorithm and a genuinely interesting example, it's hard to overstate just how impractically slow it is in practice. 💀 It is not used for any purpose in non-theoretical work, because it is literally millions of times slower than other algorithms. Straight from R.P. Brent, 2010:

"AKS is not a practical algorithm. ECPP is much faster."

This module also has ECPP, and indeed it is much faster. 😌

This implementation uses theorem 4.1 from Bernstein (2003). It runs substantially faster than the original, than the v6 revised paper with Lenstra improvements, or than the late 2002 improvements of Voloch and Bornemann. The GMP implementation uses a binary segmentation method for modular polynomial multiplication (see Bernstein's 2007 Quartic paper), which reduces to a single scalar multiplication — and that's exactly the thing GMP excels at. Because of this, the GMP implementation is likely to be faster once the input is larger than "2^33". 🎀

---

⭐ is_mersenne_prime

say "2^607-1 (M607) is a Mersenne prime" if is_mersenne_prime(607);

Give it an integer "p" and you get 1 if "p" is positive and the Mersenne number "2^p-1" is prime, 0 otherwise. Because an enormous effort has gone into testing these, a list of known Mersenne primes is used to speed things up. 🌸 Beyond the highest sequential Mersenne prime (currently 37,156,667), this performs pretesting followed by the Lucas-Lehmer test. So iconic of it. 💅

The Lucas-Lehmer test is a deterministic, unconditional test — no vibes-based guessing, it fr fr knows. And it runs hella fast compared to other primality methods for numbers of comparable size, and vastly faster than any known general-form primality proof method. But like… plot twist: while the test itself is quick, the GMP implementation is not nearly as fast as specialized programs like prime95. Also, since we lean on that table for "small" numbers, calling this function only ever actually happens for numbers with over 9.8 million digits. At that size? Bestie, just use prime95. It's giving "why are you here". 😭

---

is_ramanujan_prime

Give it an integer n and it returns 1 if n is positive and a Ramanujan prime, and 0 otherwise. Which means: anything that ramanujan_primes and nth_ramanujan_prime can spit out will give you 1, and everything else gets a 0.

Heads up 🫶 there's no cute simple function for this predicate, so what actually happens is Ramanujan primes are generated through at least n, then a search runs for n. Lowkey not it for multiple calls. 💀

is_gaussian_prime

say is_gaussian_prime(3,0);  # "2"  :  3   => 3 mod 4 => prime
say is_gaussian_prime(1,1);  # "2"  :  1+i => norm 2  => prime
say is_gaussian_prime(5,0);  # "0"  :  5   => 1 mod 4 => (2+i)(2-i)

Given two integers a and b, this returns either 0, 1, or 2, telling you whether n = a+bi is respectively a Gaussian composite, a probable Gaussian prime, or a definite Gaussian prime. This is true if and only if one of these is true:

• a = 0 and |b| is a prime congruent to 3 modulo 4.
• b = 0 and |a| is a prime congruent to 3 modulo 4.
• a and b are nonzero and a^2 + b^2 is prime.

is_delicate_prime

Hand it an integer n and you get 1 if n is positive and is a digitally delicate prime, 0 otherwise. These are numbers that are prime, but the second you change any single base-10 digit you always land on a composite. Fragile queen behavior. 💅

There's an optional second argument, the base base, which must be at least 2. That's the base used when changing digits to check for compositeness.

People call these "weakly prime" or "digitally delicate prime" numbers. And ngl, note that the first digit can be changed to a zero.

Variations that are not considered here include restricting the first digit change to non-zero values (OEIS A158124) and allowing leading zero digits to be changed ("widely DDPs").

This is the OEIS series A050249 <http://oeis.org/A050249>. With different bases, this is OEIS series A186995 <http://oeis.org/A186995>.

is_odd

Given an integer n, returns 1 if n is odd and 0 otherwise. That's it. That's the whole function. Iconic minimalism. ✨

is_even

Given an integer n, returns 1 if n is even and 0 otherwise. Yes bestie, it's just the other one. 🎀

is_divisible

Given integers n and d, returns 1 if |n| is exactly divisible by |d|, and 0 otherwise.

This matches the GMP function mpz_divisible_p. That includes its semantics when d=0, which returns 0 unless n=0.

You can pass more than one divisor, like is_divisible(1001,2,3,5,7), so you can test a whole squad of divisors in one 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.

is_congruent

Given integers n, c, and d, returns 1 if n is congruent to c modulo |d|, and 0 otherwise.

This matches the GMP function mpz_congruent_p, semantics and all — including when d=0, which returns 0 unless n=c.

is_perfect_number

Given integer n, returns 1 if n is a positive integer that equals the sum of its divisors excluding the number itself — or, same energy, a number equal to its aliquot sum. 🫶

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

Give it a single integer input n and it returns k if n = r^k for some integer r > 1, k > 1, and 0 otherwise. The k it hands back is the largest possible one. So you can just drop it in a boolean statement to check if n is a perfect power. Chef's kiss. 👑

If you give it an integer n and a non-negative integer k, it returns 1 if n is a k-th power and 0 otherwise. Like, if k=2 then boom, it's detecting perfect squares. Setting k=0 behaves like the first case: the largest root gets found and its value is returned.

If a third argument shows up, it must be a scalar reference. And if n is a k-th power, it gets set to the k-th root of n. For example:

my $n = 222657534574035968;
if (my $pow = is_power($n, 0, \my $root)) { say "$n = $root^$pow" }
# prints:  222657534574035968 = 2948^5

This corresponds to Pari/GP's ispower function with integer arguments.

is_square

Given an integer n, returns 1 if n is a perfect square, 0 otherwise. Literally identical to is_power(n,2), no cap.

This corresponds to Pari/GP's issquare function.

is_sum_of_squares

Given an integer n and an optional positive integer number of squares k, returns 1 if |n| can be written as the sum of exactly k squares. k defaults to 2. Every positive integer can be represented by 4 or more squares, so only k == 2 and k == 3 are the interesting cases, fr fr.

With k == 2 this produces the sequence OEIS A001481 <http://oeis.org/A001481>. With k == 3 this produces the sequence OEIS A000378 <http://oeis.org/A000378>.

is_powerfree

Given an integer n and an optional non-negative integer k, returns 1 if |n| has no divisor d^k, and 0 otherwise. Basically it checks whether |n| has any k-th (or higher) powers hanging out in its prime factorization. k defaults to 2.

• With k == 2 this produces the sequence of square-free integers OEIS A005117 <http://oeis.org/A005117>.
• With k == 3 this produces the sequence of cube-free integers OEIS A004709 <http://oeis.org/A004709>.
• With k == 4 this produces the sequence of biquadrate-free integers OEIS A046100 <http://oeis.org/A046100>.

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.

With k == 2 this produces the sequence OEIS A013928 <http://oeis.org/A013928>. With k == 3 this produces the sequence OEIS A060431 <http://oeis.org/A060431>.

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 gets used. It returns undef if k is less than 2 or n=0. And it returns 1 for n=1.

With k == 2 this produces the sequence OEIS A005117 <http://oeis.org/A005117>. With k == 3 this produces the sequence OEIS A004709 <http://oeis.org/A004709>.

powerfree_sum

Given an integer n and an optional non-negative integer k, returns the sum of k-powerfree positive integers less than or equal to n. k defaults to 2.

With k == 2 this produces the sequence OEIS A066779 <http://oeis.org/A066779>.

powerfree_part

Given an integer n and an optional non-negative integer k, returns the k-powerfree part of n. How? By removing the "excess" powers — in the prime factorization of n, any exponent E in P^E gets reduced down to P^(E % k). Or, put another way, all k-th powers get divided out. k defaults to 2.

When k == 2, this is sometimes also called core(n). It's the unique square-free integer d such that n/d is a square. ✨

• With k == 2 this produces the sequence OEIS A007913 <http://oeis.org/A007913>.
• With k == 3 this produces the sequence OEIS A050985 <http://oeis.org/A050985>.

With k == 2 (the default), this corresponds to Pari/GP's core function and Sage's squarefree_part function.

powerfree_part_sum

Given an integer n and an optional non-negative integer k, returns the sum of k-powerfree parts of all positive integers <= n. This is exactly equivalent to

vecsum(map { powerfree_part($_,$k) } 1..$n)

…except substantially faster, bestie. 🏎️

With k == 2 this produces the sequence OEIS A069891 <http://oeis.org/A069891>.

squarefree_kernel

Given an integer n, returns the square-free kernel of n. Also known as the integer radical. It's the largest square-free divisor of n, which is also the product of the distinct primes dividing n.

We choose to accept negative inputs here, and the result matches the input's sign. 🦋

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

sqrtint

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

This corresponds to Pari/GP's sqrtint function.

rootint

Given a non-negative integer n and positive exponent k, return the integer k-th root of n. That's the largest integer r such that r^k <= n.

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

Technically, if n is negative and k is odd, the root totally exists and equals sign(n) * |rootint(abs(n),k)|. But it was decided to follow the behavior of Pari/GP and Math::BigInt and disallow negative n instead. Rules are rules. 💅

This corresponds to Pari/GP's sqrtnint function.

logint

say "decimal digits: ", 1+logint($n, 10);
say "digits in base 12: ", 1+logint($n, 12);
my $be; my $e = logint(1000, 2, \$be);
say "largest power of 2 less than or equal to 1000:  2^$e = $be";

Given a non-zero positive integer n and an integer base b greater than 1, returns the largest integer e such that b^e <= n.

If a third argument is present, it must be a scalar reference. It will be set to b^e.

This corresponds to Pari/GP's logint function.

lshiftint

Given an integer n and an optional integer number of bits k, perform a left shift of n by k bits. If the second argument isn't provided, it's assumed to be 1. This is equivalent to multiplying by 2^k.

With negative n, it behaves exactly as described above. This is similar to how Perl behaves with use integer or use bigint — but raw Perl coerces the argument into an unsigned before left shifting, which is like, unlikely to ever be what you want. 🙈

If k is negative, a right shift is performed by |k| bits.

This corresponds to Pari/GP's shift function with a positive number of bits, and Mathematica's BitShiftLeft function.

rshiftint

Given an integer n and an optional integer number of bits k, perform a right shift of n by k bits. If the second argument isn't provided, it's assumed to be 1. This is equivalent to truncated division by 2^k.

With a negative n, the result is equal to -rshiftint(-n,k). Which means it is not "arithmetic right shift" or "logical right shift" the way those are commonly used with fixed-width registers in a particular bit format — instead it's treated as sign and magnitude, where just the magnitude gets right shifted. Different vibe, same digits. 💜

If k is negative, a left shift is performed by |k| bits.

For an actually interesting discussion of arithmetic right shift, bestie, go read Guy Steele's 1977 article "Arithmetic Shift Considered Harmful". 🎀

This is giving Pari/GP's shift function when you hand it a negative number of bits, and Mathematica's BitShiftRight. The result is literally the same as dividing by the power of 2 using tdivrem or GMP's mpz_tdiv_q_2exp.

rashiftint 💅

Give it an integer n and an optional integer number of bits k, and it does a signed arithmetic right shift of n by k bits. Don't pass that second arg? It just assumes 1. This is equivalent to floor division by 2^k.

And if k is negative? Plot twist, it does a left shift by |k| bits instead. 🦋

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

signint ✨

Given an integer n, returns the sign of n. So you get -1, 0, or 1 depending on whether n is negative, zero, or positive respectively. Iconic, simple, no notes. 👑

This corresponds to Pari/GP's sign function, GMP's mpz_sgn function, Raku's sign method, and Math::BigInt's sign method. Heads up: some of those extend to non-integers.

cmpint 🫶

Given integers a and b, returns -1, 0, or 1 if a is respectively less than, equal to, or greater than b.

The whole point of this one is to make sure Perl never silently converts the values to floating point, which can give wrong results, and also to save you from manually converting everything to bigints. It's giving safety, fr fr. 💖

This corresponds to Pari/GP's cmp function, GMP's mpz_cmp function, Math::BigInt's bcmp method, and Perl's <=> operator. Fun history fact: prior to version 6.2, GMP could return negative or positive values other than -1 and 1. 😭

addint 🌈

Given integers a and b, returns a + b.

So real talk, these integer arithmetic functions (addint, subint, mulint, add1int, sub1int, absint, negint) exist to offer exact integer arithmetic without overflow or NV conversion, while returning native integers when they fit, and bigints only when needed. It's literally the best of both worlds. Other choices include:

• Perl native operations. This is fine with small numbers, but once large enough, values will be converted to floating point (NV). This means incorrect results. 💀 Values larger than 64-bit are completely unsupported. One might expect 2^53 to be the usual point for "large enough", but not only is the NV type platform dependent, but very old 64-bit Perl will aggressively convert values to NV starting at 2^49 even with NV being a IEEE-754 double.
• use integer. Gives exact integer math as if we were using IV types in C. We are still left with 32-bit versus 64-bit platform differences, being restricted to signed type, and no support for larger values.
• Math::BigInt, Math::GMPz, etc. If one knows large values will be used, this is a good idea. Use bigint objects for all values, and all operations are methods on the objects and give correct results. This is functionally a good solution, but it will be 10 to 500 times slower as well as more memory. 🐌

All these functions accept native integers (IV/UV), bigints, and string representations of integers. Results will be in native types if possible, and as objects of the chosen bigint class otherwise. Best performance will still be had by native operations within range, or by using fast classes like Math::GMPz if most operations need it. We give correct behavior while only paying the performance penalty when needed, although there is still some overhead since we are not built into the language like Raku or Python.

subint 🍓

Given integers a and b, returns a - b.

add1int 🧚

Given integer n, returns n + 1.

sub1int 💜

Given integer n, returns n - 1.

mulint ⭐

Given integers a and b, returns a * b.

powint 🔥

Given an integer a and a non-negative integer b, returns a^b. And yes, 0^0 will return 1. 🤷

The exponent b is converted into an unsigned long.

divint 💖

Given integers a and b, returns the quotient a / b.

Floor division is used, so q is rounded towards -inf and the remainder has the same sign as the divisor b. This is the same as modern bdiv in Math::BigInt, GMP fdiv functions, and Python's integer division.

For negative inputs, this will not be identical to native Perl division, which oddly uses a truncated quotient and floored remainder. More importantly, consistent and correct 64-bit integer division in Perl is problematic. 😭

Pari/GP's \ integer division operator uses Euclidean division, which matches their divrem function. Our divint and modint operators both use floor division, which matches Raku and Python. We also have Euclidean, truncated, and ceiling division available via divrem, tdivrem, and cdivrem respectively.

modint 🌸

Given integers a and b, returns the modulo a % b.

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

Floor division is used, so q is rounded towards -inf and r has the same sign as the divisor b. This is the same as modern bmod in Math::BigInt and the GMP fdiv functions.

Like with divint, we use floor division, while Pari/GP uses Euclidean for their % integer remainder operator.

cdivint 🦋

Given integers a and b, returns the quotient a / b.

Ceiling division is used, so q is rounded towards +inf and the remainder has the opposite sign as the divisor b.

divrem ✨

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

Given integers a and b, returns a list of two items: the Euclidean quotient and the Euclidean remainder. The remainder is always non-negative (0 <= r < |b|), and the quotient is chosen to satisfy a = b*q + r.

This corresponds to Pari/GP's divrem function. There is no explicit function in Math::BigInt that gives this division method for signed inputs.

tdivrem 💅

Given integers a and b, returns a list of two items: the truncated quotient and the truncated remainder.

The resulting pair will match btdiv in Math::BigInt and btmod in Math::BigInt. This matches C99 "truncation toward zero" semantics as well.

fdivrem 🫶

Given integers a and b, returns a list of two items: the floored quotient and the floored remainder. The results will match the individual divint and modint functions, since they also use floored division.

This corresponds to Python's builtin divmod function, and Raku's builtin div and mod functions. The resulting pair will match bdiv in Math::BigInt and bmod in Math::BigInt.

cdivrem 👑

Given integers a and b, returns a list of two items: the ceiling quotient (rounded towards +inf) and the ceiling remainder. The remainder has the opposite sign from the divisor b. This allows one to perform division with rounding up.

absint 🌈

Given integer n, return |n|, i.e. the absolute value of n.

negint 💜

Given integer n, return -n.

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. Obsessed actually. 🎀

Also see lucasumod for fast computation mod n.

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

Also see lucasvmod for fast computation mod n.

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

lucasuv 🔥

($U, $V) = lucasuv(1,-2,17); # 17-th Jacobsthal, Jacobsthal-Lucas.

Given integers P, Q, and the non-negative integer k, computes both U_k and V_k for the Lucas sequence defined by P,Q. Generating both values is typically not much more time than one. Efficiency ate, no crumbs. 💅

Also see lucasuvmod for fast computation mod n.

gcd 🧚

Given a list of integers, returns the greatest common divisor. This is often used to test for coprimality https://oeis.org/wiki/Coprimality.

Each input n is treated as |n|.

lcm 🌸

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

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

gcdext 💖

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

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

chinese 🦋

say chinese( [14,643], [254,419], [87,733] );  # 87041638

Solves a system of simultaneous congruences using the Chinese Remainder Theorem (with extension to non-coprime moduli). A list of [a,n] pairs are taken as input, each representing an equation x ≡ a mod |n|. If no solution exists, undef is returned. If a solution is returned, the modulus is equal to the lcm of all the given moduli (see lcm). In the standard case where all values of n are coprime, this is just the product. The a values must be integers, while the n values must be non-zero integers. Like other mod functions, we use abs(n).

Comparison to similar functions in other software:

Math::ModInt::ChineseRemainder:
  cr_combine( mod(a1,m1), mod(a2,m2), ... )

Pari/GP:
  chinese( [Mod(a1,m1), Mod(a2,m2), ...] )

Mathematica:
  ChineseRemainder[{a1, a2, ...}, {m1, m2, ...}]

SAGE:
  crt( [a1,m1], [a2,m2], ... )
  crt(a1,m1,a2,m2,...)
  CRT_list( [a1,a2,...], [m1,m2,...] )

chinese2 ✨

Functions like chinese but returns two items: the remainder and the modulus. If a solution exists, the second value (the final modulus) is equal to the lcm of the absolute values of all the given moduli.

If no solution exists, both return values will be undef. Not me getting zero results and being sad about it. 😭

frobenius_number 👑

Finds the Frobenius number of a set of positive integers. This is the largest positive integer that cannot be represented as a non-negative linear combination of the input set. Each set element must be positive (all elements greater than zero) and setwise coprime: gcd(a1,a2,...,an) = 1.

This is sometimes called the "coin problem". 🪙

This corresponds to Mathematica's FrobeniusNumber function. Matching their API, we return -1 if any set element is 1.

💖 vecsum

say "Totient sum 500,000: ", vecsum(euler_phi(0,500_000));

Returns the sum of all its arguments, and every single one of them has to be an integer, okay? 😌

This is giving majorly similar-to-sum0 in List::Util energy, BUT — and this is the important part, bestie — there's one huge difference. List::Util yeets all its inputs into doubles and hands you back a double, which means your results go absolutely feral and wrong with large integers. vecsum sums (signed) integers and gives you the untruncated result. No rounding betrayal. 💅

The processing happens on native integers whenever it can, including using a 128-bit running sum in the C code. Queen behavior.

---

✨ vecprod

say "Totient product 5,000: ", vecprod(euler_phi(1,5_000));

Returns the product of all its arguments — again, they all have to be integers. This one's like product in List::Util, except it keeps everything as integers and automatically switches to bigints if it needs to. So it literally cannot lose your precision. 🫶

---

🦋 vecmin

say "Smallest Totient 100k-200k: ", vecmin(euler_phi(100_000,200_000));

Returns the minimum of all its arguments, each of which must be an integer. Very min-from-List::Util coded, except — say it with me — List::Util turns all inputs into doubles and returns a double, which gives you incorrect results with large integers. 🙈

vecmin validates and compares everything as integers. Yes, that validation step makes it a little slower than min in List::Util, but it's there so you can't accidentally, unintentionally trip into float land. It's a safety feature, not a bug, periodt.

---

🌸 vecmax

say "Largest Totient 100k-200k: ", vecmax(euler_phi(100_000,200_000));

Returns the maximum of all its arguments, each of which must be an integer. Same max-in-List::Util vibe, same confession: List::Util doubles everything and returns a double, and that gives incorrect results with large integers. 😭

vecmax validates and compares all results as integers. It's a touch slower than max in List::Util because of the validation, but that's the whole point — it prevents accidental and unintentional float nonsense. We stan an honest function.

---

💜 vecreduce

say "Count of non-zero elements: ", vecreduce { $a + !!$b } (0,@v);
my $checksum = vecreduce { $a ^ $b } @{twin_primes(1000000)};

Does a reduce operation via a left fold. Takes a block and a list as arguments. The block uses the special local variables a and b, which are the accumulation and the next element respectively, and whatever the block returns becomes the new accumulation. No initial element is used, so you get undef back if you pass an empty list. 🕳️

The interface is exactly the same as reduce in List::Util — this was done on purpose to increase portability and minimize confusion. If you want the full discourse on reduce with empty or singular-element lists, see chapter 7 of Higher Order Perl (or honestly, many other references). Friendly tip from the source: it's often a good idea to give an identity element as the first list argument. 🎀

Look, you could do vecmin, vecmax, vecsum, vecprod, etc. with this function, and it'd be fairly easy — but it will not be as efficient. Meanwhile there are a wide variety of other functions you can build easily with reduce, which makes it a genuinely useful tool. iykyk. 🔥

---

⭐ vecany

Returns true if any element of a list satisfies a block. See vecfirst. ✨

---

💖 vecall

Returns true if all elements of a list satisfy a block. See vecfirst.

---

🌈 vecnone

Returns true if no element of a list satisfies a block. See vecfirst.

---

🎀 vecnotall

Returns true if not all elements of a list satisfy a block. See vecfirst.

---

👑 vecfirst

say "all values are Carmichael" if vecall { is_carmichael($_) } @n;

Short circuit evaluations of a block over a list. Takes a block and a list as arguments. The block gets called with $_ set to each list element, and evaluation keeps going until either all the list values have been evaluated or the result condition can already be determined. For instance, in the vecall example above, evaluation stops the literal instant any value returns false. Efficiency queen. 💅

The interface is exactly the same as the any, all, none, notall and first functions in List::Util. This was done to increase portability and minimize confusion. Unlike the other vector functions like vecmin, vecmax, vecsum, etc., there is no added value to using these versus the ones from List::Util. They're here for convenience, that's it. 🫶

You can fairly easily map these operations to scalar(grep {...} @n), but that does not short-circuit and it's less obvious. So like… why would you. 😌

---

🔥 vecfirstidx

say "first Carmichael is index ", vecfirstidx { is_carmichael($_) } @n;

Returns the index of the first element in a list that evaluates to true. Just like vecfirst, except it hands you the index instead of the value. Returns -1 if the item could not be found. 🕵️

This interface matches firstidx and first_index from List::MoreUtils.

---

🧚 vecextract

say "Power set: ", join(" ",vecextract(\@v,$_)) for 0..2**scalar(@v)-1;
@word = vecextract(["a".."z"], [15, 17, 8, 12, 4]);

Extracts elements from an array reference based on a mask, and returns the result as an array. The mask is either an unsigned integer — which gets treated as a bit mask — or an array reference holding integer indices. 🎯

If the second argument is an integer, every bit set in the mask makes the corresponding element from the array reference come back to you. Bits are read from the right, so a mask of 1 returns the first element, while 5 returns the first and third. The mask may be a bigint. No cap.

If the second argument is an array reference, then its elements get used as zero-based indices into the first array. Duplicate values are allowed and the ordering is preserved. Given that Perl already has fully functional array slices in the language, this exists for completeness with Pari/GP. These two are equivalent:

vecextract($aref, $iref);
@$aref[@$iref];

---

🍓 vecuniq

my @vec = vecuniq(1,2,3,2,-10,-100,1);  # returns (1,2,3,-10,-100)

Given an array of integers, returns an array with all the duplicate entries removed. The original ordering is preserved. All values must be defined. ✨

This is similar to List::Util::uniqint (the integer comparison version of List::Util::uniq). Unlike the more generic List::Util::uniq and List::MoreUtils::XS::uniq, all inputs must be integers. With native integers, our function is 2-10x faster. Slay.

---

💜 vecfreq

# Produce frequency hash:
my %h = vecfreq(1,2,2,2,3,1,4);   #  (1=>2, 2=>3, 3=>1, 4=>1)
# Print most common value:
say vecreduce { $h{$a} > $h{$b} ? $a : $b } keys %h;

Given an array of items, returns a hash where each key holds a unique item, and the associated value is how many times it showed up in the array. Counting queen. 👑

This is identical to List::MoreUtils::frequency. It's typically faster when you give it only native integers.

It's also very similar to the Pari/GP function matreduce for vectors, and to Python's Counter. So if you know those, you already know her.

---

🦋 vecsingleton

my @solo = vecsingleton(1,4,17,1,17,-8);  # (4,-8)
# Same but slower:
my %h = vecfreq(@n);
my @onlyuniqs = grep { $h{$_} == 1 } @n;

Given an array of items, returns an array with all entries removed that appear more than once in the list. The original ordering is preserved. She said "if you're not unique, you're not coming." 💅

This is identical to List::MoreUtils::singleton. When you feed it only native integers, it's typically 2 to 10x faster. Iconic.

---

🌈 vecsort

my @sorted = vecsort(1,2,3,2,-10,-100,1);   # returns (-100,-10,1,1,2,2,3)
my @sorted = vecsort([1,2,3,2,-10,-100,1]); # same

Numerically (ascending) sorts a list of integers. The input is either a list, or a single array reference that holds the list. Simple as that. 🌸

All values must be defined and integers. They may be any mix of native IV, native UV, strings, and bigints.

Perl's built-in numerical sort can sometimes give incorrect results for typical cases we encounter — not her fault, but still. Prior to version 5.26 (2017), large 64-bit integers were turned into NV (floating point) types. With all current versions of Perl, strings get turned into NV types even if they are the text of a 64-bit integer. 😭

In scalar context, vecsort returns the number of items without sorting (but after input validation). This is expected and honestly what we typically want. Say you only want the number of divisors — you call in scalar context and get the count without needing the actual sort. Having the same results from $x = vecsort(5,6,7) and @v = vecsort(5,6,7); $x=@v; is exactly what we want. This contrasts with Perl's built-in sort, which in scalar context has undefined behaviour (in all current versions of perl it returns undef). In particular, this forces all programs to use a workaround if they want to return the results of sorting an array. See Perl 5 issue 12803 (https://github.com/Perl/perl5/issues/12803) for some discussion with no resolution. 💀

Using an array reference as input is slightly faster. ✨

This is almost always faster than Perl's built-in numerical sort: @a = sort { $a <=> $b } @a. See the performance section for more information.

---

🔥 vecsorti

my @arr = map { irand } 1..100000;
vecsorti \@arr;

Given an array reference of integers, numerically (ascending) sorts the integers in-place. The array reference is also returned for convenience. 🫶

This is more efficient than vecsort. Perl's sort has this optimization built-in when doing straightforward sorting on non-references.

---

🧚 vecequal

my $is_equal = vecequal( [1,2,-3,[4,5,undef]], [1,2,-3,[4,5,undef]] );

Compare two arrays for equality, nested arrays included. The values inside the two input array references must be either an array reference, a scalar, or undef. Simple integers get tested with integer comparison, while other scalars use string comparison. 🎯

Important: this is a vector comparison, not a set comparison, so ordering matters, bestie. For the sake of wider applicability, non-integers are allowed. Types other than integers and strings (e.g. floating point values) are not guaranteed to have consistent results. 🙈

No circular reference detection is performed. Like, at all. Don't come for her if you hand her a loop.

Performance with XS is 3x to 100x faster than perl looping or modules like Array::Compare, Data::Cmp, match::smart, List::Compare and Algorithm::Diff. Those modules have additional functionality, so this is not a complete comparison — just so we're being fair. 💅

---

⭐ vecmex

my $minimum_excluded = vecmex(0,1,2,4,6);  # returns 3

Given a list of non-negative integers, returns the smallest non-negative integer that is not in the list. "mex" is short for "minimum excluded". You can think of the list as a set, and the return value is the minimum of the set complement. Repeated values are allowed in the list, no problem. ✨

vecmex() = 0. vecmex(0,1,2,...,w) = w+1. Cute little edge cases, love that for us.

---

💜 vecpmex

my $minimum_excluded = vecpmex(1,2,4,6);  # returns 3

Given a list of positive integers, returns the smallest positive integer that is not in the list. "mex" is short for "minimum excluded". Same deal as before: the list can be seen as a set, and the return value is the minimum of the set complement. Repeated values are allowed in the list. 🦋

vecpmex() = 1. vecpmex(1,2,...,w) = w+1.

vecslide 🛝

@pairsum = vecslide {$a+$b} 1..5;  # returns (1+2,2+3,3+4,4+5)

say for vecslide { "$a->[0] $b->[1]" }
    (["hello","world"], ["goodbye","friends"], ["love","hate"]);
# hello friends
# goodbye hate

Bestie, you hand it a code block plus a list, and it runs that block once for every pair in the list, dropping the two values of each pair into the local $a and $b. Dynamic scoping, old school, iconic.

And no, the list isn't picky about what's in it — peep that second example, it's literally shipping arrays. Whatever you want, fr fr.

This is identical to List::MoreUtils::slide. Same vibe, different module 💅.

toset 🎀

my $set = toset(52,-6,14,-6,0);  # $set = [-6,0,14,52]
say "number of elements in set: ",scalar(@$set);
say "smallest value: ",$set->[0];
say "largest value: ",$set->[-1];

You give it a list of integers, it hands you back an array reference that is the integer set. Numerically sorted, duplicates yeeted, glow-up complete. ✨

The input array must contain only integers — signed integers, bigints, objects that evaluate to integers, strings representing integers, all of those are fine. This "set form" is the optimal form for all the set operations that follow, so convert early and thank yourself later.

Once your set is in this form, its size is just the length. Minimum and maximum? Trivial, it's $set->[0] and $set->[-1]. Everything in the output gets typed as either a native integer (IV or UV) or a bigint. No in-between, no chaos 💀.

setinsert 📥

my $s=[-10..-1,1..10];
setinsert($s, 0);            # $s is now [-10..10]
setinsert($s, [5,10,15,20]); # $s is now [-10..10,15,20]

Give it two array references of integers in set form, and it inserts every element of the second set into the first set — then returns the number of elements that were actually inserted.

Or: give it one array reference in set form followed by zero or more integer scalars (unordered? duplicates? whatever, doesn't matter), and it inserts all those list values into the first set, again returning how many got inserted. Same thing as wrapping your list in toset, just convenient and possibly more efficient. So nice of it. 🫶

You can think of this as an in-place setunion.

The one or two sets (array references) must be in set form (numerically sorted with no duplicates) or the results are undefined. Not "a little wrong", bestie — undefined. 🙈

setremove 🗑️

my $s=[-10..10];
setremove($s, 0);            # $s is now [-10..-1,1..10]
setremove($s, [5,10,15,20]); # $s is now [-10..-1,1..4,6..9]

Two array references of integers in set form: it removes every element of the second set from the first set and returns the number of elements that were removed.

Or a set form array reference followed by zero or more integer scalars (possibly unordered, possibly duplicating), removing all of those list values from the first set and returning how many got removed. Basically wrapping the list in toset but convenient and possibly more efficient.

This is the in-place setminus, in case you're keeping score.

Same deal as always: the one or two sets must be in set form (numerically sorted, no duplicates) or the results are undefined. 💀

setinvert 🔄

my $s=[-10..10];
setinvert($s, 0);            # $s is now [-10..-1,1..10]
setinvert($s, [5,10,15,20]); # $s is now [-10..-1,1..4,6..9,15,20]

Two array references of integers in set form: it inverts the containment status in the first set for each element of the second set. Which is a fancy way of saying — for each element of the second set, if it's not in the first set, insert it; if it is in the first set, remove it. Flip the switch, bestie. 🔀

Or a set form array reference followed by zero or more integer scalars (possibly unordered, possibly duplicating), which does the same as if the list was wrapped in toset.

It returns an integer telling you how many values were inserted minus how many were deleted.

This is the in-place setdelta. Look at her go 🦋.

The one or two sets (array references) must be in set form (numerically sorted with no duplicates) or the results are undefined. I keep saying it because it keeps being true.

setcontains 🔍

my $has_element = setcontains( [-12,1..20], 15 );
my $is_subset   = setcontains( [-12,1..20], [-12,5,10,15] );

Two sets (array references of numerically sorted, de-duplicated integers) walk in. It returns 1 or 0 telling you whether the second argument is a subset of the first set — i.e. whether every single element of the second argument is a member of the first set. No stragglers allowed ✋.

Or a set plus zero or more integers in any form (unordered? duplicated? sure!), which does the same as if you'd wrapped that list in toset.

If the first array reference is not in set form (numerically sorted with no duplicates, and no string forms), the result is undefined. It is unlikely to give a correct answer, like truly unlikely. Use toset to convert an arbitrary integer list into set form first. Do it. 😭

setcontainsany 🎯

# True if there is any intersection between the two sets
my $intersects = setcontainsany($set1,$set2);
my $has_one_of = setcontainsany( [-12,1..20], -14,0,1,100 );  # true

Two sets (array references of numerically sorted de-duplicated integers): returns 1 or 0 indicating whether any element of the second set is an element of the first set. Just one. One is enough. One crumb of overlap 😌.

Alternately, a set followed by a list of unordered integers does the same thing, as if the list was wrapped in toset.

There's a lil bit of functional duplication round here, and honestly? Obsessed. Checking for disjoint sets can be done with any of these:

my $dj1 = set_is_disjoint($set1, $set2);
my $dj2 = scalar(@{setintersect($set1, $set2)}) == 0;
my $dj3 = !setcontainsany($set1, $set2);

This function requires the array reference inputs be in set form or the result is undefined. But in return it can be thousands of times faster for large sets. Thousands. Not a typo. 🔥

setbinop ✖️

my $sumset = setbinop { $a + $b } [1,2,3], [2,3,4];  # [3,4,5,6,7]
my $difset = setbinop { $a - $b } [1,2,3], [2,3,4];  # [-3,-2,-1,0,1]
my $setsum = setbinop { $a + $b } [1,2,3];           # [2,3,4,5,6]

A code block plus two array references of integers: it treats them as integer sets and builds you a brand new set by running the cross product through that block. Only one array reference? It'll use it against itself, self-love era 💖.

The result comes back in set form (numerically sorted, no duplicates). The input sets are not aliased inside the block — so go ahead, mess with $a and $b all you want, nothing leaks outside the block.

This corresponds to Pari's setbinop function. Our function uses much less memory, as of Pari 2.18.1. We ate that one, ngl 👑.

sumset ➕

Give it two array references of integers, it treats them as integer sets and returns the sumset as a set — a sorted, de-duplicated array reference, gorge.

Only one array reference given? It'll be used for both. Which means applying sumset to a single set is a super common move, and honestly iconic.

This is equivalent to:

my %r;  my @A=(2,4,6,8);  my @B=(3,5,7);
forsetproduct { $r{vecsum(@_)}=undef; } \@A,\@B;
my $sumset = [vecsort(keys %r)];

or

my $sumset1 = setbinop { addint($a,$b) } [1,2,3];
my $sumset2 = setbinop { addint($a,$b) } [1,2,3], [2,3,4];

In Mathematica one can use Total[Tuples[A,B],{2}]. In Pari/GP one can use setbinop((a,b)-"a+b,X,Y)>.

setunion 🤝

Given exactly two array references of integers, it treats them as sets and returns the union as a set. The returned set has every element that appears in either input set. Everyone's invited to this party babes 🎉.

It's more efficient if the input is already in set form (numerically sorted, no duplicates), and the result will be in set form regardless.

This corresponds to Pari's setunion function, Mathematica's Union function, and Sage's union function on Set objects. Funny how they all agree on this one.

setintersect 💞

my $commonset = setintersect($set1,$set2);
my $is_disjoint = 0 == @$commonset;  # scalar size of the intersection

Exactly two array references of integers: treats them as sets and returns the intersection as a set. The returned set has all the elements that appear in both input sets. Mutuals only 💅.

More efficient if the input is in set form (numerically sorted, no duplicates). Result comes back in set form.

This corresponds to Pari's setintersect function, Mathematica's Intersection function, and Sage's intersection function on Set objects.

setminus ➖

Exactly two array references of integers: treats them as sets and returns the difference as a set. The returned set has all the elements that appear in the first set but not in the second. It's giving "you're not on the list" energy 🚫.

More efficient if the input is in set form (numerically sorted, no duplicates). Result will be in set form.

This corresponds to Pari's setminus function, Mathematica's Complement function, and Sage's difference function on Set objects.

setdelta 🔀

Exactly two array references of integers: treats them as sets and returns the symmetric difference as a set. The returned set has all the elements that appear in only one of the two input sets. Both-or-neither, no thirds, periodt.

More efficient if the input is in set form (numerically sorted, no duplicates). Result will be in set form.

This corresponds to Pari's setdelta function, Mathematica's SymmetricDifference function, and Sage's symmetric_difference function on Set objects.

is_sidon_set 🌸

Given an array reference of integers, it treats it as a set and returns 1 if it's a Sidon set (sometimes called a Sidon sequence), and 0 otherwise.

To be a Sidon set, all elements must be non-negative and all pair-wise sums a_i + a_j (i >= j) have to be unique. No twin sums allowed in this house.

All finite Sidon sets are Golomb rulers, and all Golomb rulers are Sidon. Cute mutual slay 🧚.

is_sumfree_set 🍓

Given an array reference of integers, it treats it as a set and returns 1 if it's a sum-free set, and 0 otherwise.

A sum-free set is one where no sum of two elements from the set is equal to any element of the set. Which is just saying: the set and its sumset are disjoint. They don't know each other. Strangers, fr fr.

set_is_disjoint 🙅

Given two array references of integers, it treats them as sets and returns 1 if the sets have no elements in common, 0 otherwise.

This corresponds to Mathematica's DisjointQ function.

set_is_equal 🪞

Given two array references of integers in set form, returns 1 if the sets have all elements in common, 0 otherwise.

This function works even if the inputs aren't sorted — she's flexible like that. But if they are sorted (proper set form) then vecequal can be used instead, and it's typically much faster.

set_is_subset 📚

Given two array references of integers in set form, returns 1 if the first set also contains all the elements of the second set, 0 otherwise.

The setcontains function can be used equivalently, and it doesn't require the second list to be in set form. So pick your fighter.

This corresponds to Mathematica's SubsetQ function (is B a subset of A).

set_is_proper_subset 🪆

Given two array references of integers in set form, returns 1 if the first set also contains all the elements of the second set but they are not equal, 0 otherwise. Basically subset but it has to be strict about it — the size of the first set must be strictly larger than the second. No equality here, we're exclusive 💜.

set_is_superset ✨

okay so you hand it two array references of integers in set form, right? it returns 1 if the second set ALSO has every single element of the first set, and 0 if not. that's the whole vibe.

and bestie, the setcontains function does the exact same thing — just with the arguments reversed. equivalently. 💅

set_is_proper_superset 👑

two array refs of integers in set form again — returns 1 if the second set contains all the elements of the first set but they're not equal, 0 otherwise. and the size of the second set has to be strictly larger than the first. no ties, it's giving strict. ✨

set_is_proper_intersection 🦋

two array refs of integers in set form. returns 1 if the two sets have at least one element in common AND each of the two sets has at least one element that isn't present in the other one. returns 0 otherwise. so like… they overlap but neither one swallows the other. 💜

todigits 🔢

say "product of digits of n: ", vecprod(todigits($n));

give it an integer n and it hands back an array of the digits of |n|. an optional second integer argument specifies a base (default 10). base 2? you get an array of binary digits of n. cute.

an optional third argument specifies a length for the returned array — the result will either have upper digits truncated or have leading zeros added. most often used with base 2, 8, or 16.

real talk: the values returned may be read-only. also todigits(0) returns an empty array. the base must be at least 2, and is limited to an int. length must be at least zero and is limited to an int.

this corresponds to Pari's digits and binary functions, and Mathematica's IntegerDigits function. 🌈

todigitstring 🎀

# arguments are:  input integer, base (optional), truncate (optional)
say "decimal 456 in hex is ", todigitstring(456, 16);
say "last 4 bits of $n are: ", todigitstring($n, 2, 4);

basically todigits but it gives you back a string. for bases <= 10, this is equivalent to joining the array returned by todigits.

the first argument n is the input integer, and the sign is just ignored, whatever. if no other arguments are given, this just returns the string of n. an optional second argument is the base base, which must be between 2 and 36. no prefix such as 0x will be added, and all bases over 9 use lower case a to z.

an optional third argument k requires the result to be exactly k digits. this truncates to the last k digits if the result has k or fewer digits, or zero extends if the result has more digits.

this corresponds to Mathematica's IntegerString function. 💖

fromdigits 💫

say "hex 1c8 in decimal is ", fromdigits("1c8", 16);
say "Base 3 array to number is: ", fromdigits([0,1,2,2,2,1,0],3);

this takes either a string or an array reference, plus an optional base (default 10). with a string, each character will be interpreted as a digit in the given base, with both upper and lower case denoting values 11 through 36. with an array reference, the values indicate the entries in that location, and values larger than the base are allowed (results are carried). the result is a number — either a native integer or a bigint.

this corresponds to Pari's fromdigits function and Mathematica's FromDigits function. 🧚

tozeckendorf 🐚

say tozeckendorf(24);                     #  "1000100"
say fromdigits(tozeckendorf(24),2);       #  68

given a non-negative integer n, return the Zeckendorf representation as a binary string. this represents n as a sum of nonconsecutive Fibonacci numbers. each set bit indicates summing the corresponding Fibonacci number, e.g. 24 = 21+3 = F(8)+F(4). F(0)=0 and F(1)=1 are not used. this is sometimes also called Fibbinary or the Fibonacci base.

the restriction that consecutive values are not used ("11" cannot appear) is required to create a unique mapping to the positive integers. a simple greedy algorithm suffices to construct the encoding.

say reverse(tozeckendorf($_)).'1'  for  1..20

shows the first twenty Fibonacci C1 codes (Fraenkel and Klein, 1996). this is an example of a self-synchronizing variable length code.

this corresponds to Mathematica's ZeckendorfRepresentation[n] function. Also see Math::NumSeq::Fibbinary and Data::BitStream::Code::Fibonacci. ✨

fromzeckendorf 🔄

say fromzeckendorf("1000100");            #  24
say fromzeckendorf(todigitstring(68,2));  #  24

given a binary string in Zeckendorf representation, return the corresponding integer. the string may not contain anything other than the characters 0 and 1, and must not contain 11. the resulting number is the sum of the Fibonacci numbers in the position starting from the right (the Fibonacci index is offset by two, as F(0)=0 and F(1)=1 are not used). 🌸

sumdigits ➕

# Sum digits of primes to 1 million.
my $s=0; forprimes { $s += sumdigits($_); } 1e6; say $s;

given an input n, return the sum of the digits of n. any non-digit characters of n are ignored (yep, including negative signs and decimal points). this is similar to the command vecsum(split(//,$n)) but faster, allows non-positive-integer inputs, and can sum in other bases.

an optional second argument indicates the base of the input number. this defaults to 10, and must be between 2 and 36. any character that is outside the range 0 to base-1 will be ignored.

if no base is given and the input number n begins with 0x or 0b then it will be interpreted as a string in base 16 or 2 respectively.

regardless of the base, the output sum is a decimal number.

this is similar but not identical to Pari's sumdigits function from version 2.8 and later. the Pari/GP function always takes the input as a decimal number, uses the optional base as a base to first convert to, then sums the digits. this can be done with either vecsum(todigits($n, $base)) or sumdigits(todigitstring($n,$base)). Math::BigInt version 1.999818 has a similar digitsum function. 💜

valuation 🪜

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

given integer n and non-negative integer k, returns the number of times n is divisible by k. this is a very limited version of the algebraic valuation — here it is just applied to integers.

k must be greater than 1. |n| is used, |n| = 0 returns undef, and |n| = 1 returns zero.

this corresponds to Pari and SAGE's valuation function. ⭐

hammingweight ⚖️

given an integer n, returns the binary Hamming weight of abs(n). this is also called the population count, and is the number of 1s in the binary representation. this corresponds to Pari's hammingweight function for t_INT arguments.

is_square_free 🧼

say "$n has no repeating factors" if is_square_free($n);

given integer n, returns 1 if |n| has no repeated factor.

is_cyclic 🔁

given integer n, returns 1 if n is positive and cyclic in the number theory sense, and returns 0 otherwise. a cyclic number n has only one group of order n. n and φ(n) are relatively prime.

this is the OEIS series A003277 <http://oeis.org/A003277>.

is_carmichael 🎭

for (1..1e6) { say if is_carmichael($_) } # Carmichaels under 1,000,000

given an integer n, returns 1 if n is positive and a Carmichael number, and returns 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.

for inputs larger than 50 digits after removing very small factors, this uses a probabilistic test since factoring the number could take unreasonably long. the first 150 primes are used for testing. any that divide n are checked for square-free-ness and the Korselt condition, while those that do not divide n are used as the pseudoprime base. the chances of a non-Carmichael passing this test are less than 2^-150.

this is the OEIS series A002997 <http://oeis.org/A002997>.

is_quasi_carmichael 🌀

given an integer n, returns 0 if n is negative or not a quasi-Carmichael number, and returns the number of bases otherwise. these are square-free composites that satisfy p+b divides n+b for all prime factors p of n and for one or more non-zero integer b.

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

is_semiprime ✌️

given an integer n, returns 1 if n is positive and a semiprime, and returns 0 otherwise. a semiprime is the product of exactly two primes.

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

is_almost_prime 🧮

say is_almost_prime(6,2169229601);  # True if n has exactly 6 factors

given non-negative 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.

functionally identical but possibly faster than prime_bigomega(n) == k.

is_omega_prime 🔍

say is_omega_prime(6,2169229601);  # True if n has 6 distinct factors

given non-negative integers k and n, returns 1 if n has exactly k distinct prime factors (not counting multiplicity), and 0 otherwise. with k=1, this is the same as is_prime_power.

functionally identical but possibly faster than prime_omega(n) == k.

is_chen_prime 🌟

given non-negative integer n return 1 if n is a Chen prime. that is, if n is prime and n+2 is either a prime or semi-prime.

is_fundamental 🏛️

given an integer d, returns 1 if d is a fundamental discriminant, 0 otherwise. we consider 1 to be a fundamental discriminant.

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

this corresponds to Pari's isfundamental function.

is_totient 🪞

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

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

also see inverse_totient which gives the count or list of values that produce a given totient. this function is more efficient than getting the full count or list. ✨

is_pillai 💊

given a non-negative integer n, if there exists a v where v! % n == n-1 and n % v != 1, then the least v is returned. otherwise 0.

for n prime, non-zero return values give OEIS series A063980 <http://oeis.org/A063980>. the non-zero values returned produce OEIS series A063828 <http://oeis.org/A063828>.

is_polygonal 🔷

given an integer x and a positive integer s greater than 2, return 1 if x is an s-gonal number, and return 0 otherwise.

if a third argument is present, it must be a scalar reference. it will be set to n if x is the nth s-gonal number. if the function returns 0, then it will be unchanged.

this corresponds to Pari's ispolygonal function.

💖 is_congruent_number

Hand it a non-negative integer n and it hands you back 1 if n is literally the area of a rational right triangle, and 0 if it's not. That's it, that's the whole vibe.

It answers the congruent number problem using Tunnell's theorem — which is riding hard on the Birch Swinnerton-Dyer conjecture, bestie, we're in deep math territory here. It also runs an extensive filter for known non-congruent families, pulling from the works of Bastien (1915), Lagrange (1974), Monsky (1990), Serf (1991), Iskra (1996), Feng (1996), Reinholz et al. (2013), Cheng and Guo (2018 and 2019), Das and Saikia (2020), and Evink (2021). A whole century of mathematicians in one function, we love to see it. ✨

---

🌸 cornacchia

Give it non-negative integers d and n and it goes looking for solutions (x,y) to the equation x^2 + d y^2 = n. If there's no solution? It just returns undef. No drama.

When n is a prime, this is done using Cornacchia's algorithm.

For non-prime n, it uses a combination of Cornacchia-Smith on all roots, plus a loop to chase down the harder cases. Which means — bestie, listen — we will always return a solution. 💅

There will often be multiple solutions, but only one is returned. Don't be greedy.

---

🦋 contfrac

  my @CF = contfrac(415,93);
  # CF = (4,2,6,7)  =>  4+(1/(2+1/(6+1/7))) = 415/93
  #                     ^     ^    ^   ^

Give it an integer n and a positive integer d, and it gives you back a list with the simple continued fraction representation of the rational n / d.

This is a subset of Pari's contfrac function, Mathematica's ContinuedFraction[n/d] function, and Sage's continued_fraction function. Look at us, matching the icons. 👑

---

🎀 from_contfrac

  my($N,$D) = from_contfrac(4,2,6,7);  # N = 415, D = 93

Give it an array of integers representing the simple continued fraction, and it gives you back the rational n / d as two integers.

Heads up on the rules though: the first input value represents the whole part, and it may be zero or negative. But all successive input values must be non-negative and non-zero. No exceptions, periodt.

This is a subset of Pari's contfracpnqn function, Mathematica's FromContinuedFraction[list] function, and one value of Sage's convergent(n) method.

---

🧚 next_calkin_wilf

  ($n,$d) = next_calkin_wilf($n,$d);

Give it two positive coprime integers n and d representing the rational n / d, and it returns the next value in the breadth-first traversal of the Calkin-Wilf tree of rationals, as a two-element list. Tree traversal but make it cute. 🌈

The Calkin-Wilf tree has an entry for all positive rationals in lowest form, with each one appearing only once. It's not a binary search tree over the positive rationals like the Stern-Brocot tree is — but it is more efficient to traverse in both depth and breadth order. Trade-offs, bestie. 🤷

This corresponds to Julia's Nemo next_calkin_wilf function.

Also! You can efficiently iterate through OEIS series A002487 OEIS series A002487 <http://oeis.org/A002487> with this. Slay.

---

🫶 next_stern_brocot

  ($n,$d) = next_stern_brocot($n,$d);

Give it two positive coprime integers n and d representing the rational n / d, and it returns the next value in the breadth-first traversal of the Stern-Brocot tree of rationals, as a two-element list.

The Stern-Brocot tree also has an entry for all positive rationals in lowest form, with each one appearing only once. Read left-to-right on each row, the numbers appear in ascending order. It is a binary search tree over the positive rationals (this was exactly Brocot's motivation, king behavior honestly). It is not as efficient as next_calkin_wilf though. 💀

This produces OEIS series A007305 OEIS series A007305 <http://oeis.org/A007305> (numerators) and OEIS series A047679 OEIS series A047679 <http://oeis.org/A047679> (denominators).

---

💜 calkin_wilf_n

  my $idx = calkin_wilf_n($n,$d);

Give it two positive coprime integers n and d representing the rational n / d, and it returns the index in the breadth-first traversal of the Calkin-Wilf tree of rationals.

This corresponds to the xy_to_n method in Math::PlanePath::RationalsTree with tree_type = 'CW'>. ✨

---

🌟 stern_brocot_n

  my $idx = stern_brocot_n($n,$d);

Give it two positive coprime integers n and d representing the rational n / d, and it returns the index in the breadth-first traversal of the Stern-Brocot tree of rationals.

This corresponds to the xy_to_n method in Math::PlanePath::RationalsTree with tree_type = 'SB'>. 🎯

---

🍓 nth_calkin_wilf

  ($n,$d) = nth_calkin_wilf($idx);

Give it a positive integer i, and it returns the rational in the corresponding index in the breadth-first traversal of the Calkin-Wilf tree of rationals.

This corresponds to the n_to_xy method in Math::PlanePath::RationalsTree with tree_type = 'CW'>. 🌸

---

🦄 nth_stern_brocot

  ($n,$d) = nth_stern_brocot($idx);

Give it a positive integer i, and it returns the rational in the corresponding index in the breadth-first traversal of the Stern-Brocot tree of rationals.

This corresponds to the n_to_xy method in Math::PlanePath::RationalsTree with tree_type = 'SB'>. 💫

---

🌈 nth_stern_diatomic

  $n = nth_stern_diatomic($idx);

Give it a non-negative integer i, and it returns the i-th Stern diatomic number. This is sometimes called fusc(i) (Dijkstra), Stern's diatomic series, or the Stern-Brocot sequence. Quick note bestie — that last name can be easily confused with the Stern-Brocot tree, so watch out. 👀

This corresponds to Sidef's fusc function. See also next_calkin_wilf where the sequence of numerators generates this sequence.

This produces OEIS series A002487 OEIS series A002487 <http://oeis.org/A002487>. Iconic.

---

🧁 farey

  #    F[3] = 0/1  1/3  1/2  2/3  1/1
  #
  say scalar farey(3);   #  5
  my @F3 = farey(3);     #  ([0,1], [1,3], [1,2], [2,3], [1,1])
  my $F33 = farey(3,3);  #  [2/3] = $F3[3]
  # Print the list in readable form
  say join " ",map { join "/",@$_ } farey(3);

Give it a single positive integer n and it returns the Farey sequence of order n. In scalar context, it returns the length without computing terms. In array context, it returns a list with each rational as a 2-entry array reference. Context-sensitive king. 👑

Give it two values instead — a positive integer n and a non-negative integer k — and it returns the k-th entry of the order n Farey sequence. The index starts at zero, so it matches using the full list as an array. If k is larger than the number of entries, undef is returned. No crash, just vibes.

This corresponds to Mathematica's FareySequence function (their two-argument version is 1-based rather than 0-based, so mind the offset).

The lengths are OEIS series A005728 OEIS series A005728 <http://oeis.org/A005728>. The numerators are OEIS series A006842 OEIS series A006842 <http://oeis.org/A006842>. The denominators are OEIS series A006843 OEIS series A006843 <http://oeis.org/A006843>. ✨

---

🎀 next_farey

  my $next = next_farey(9,[5,9]);  # returns [4,7]

Give it a positive integer n and a 2-element array reference containing a non-negative integer p and a positive integer q, and it returns the next rational appearing after p/q in the order n Farey sequence. Returns undef if p/q is greater than or equal to one. That's the ceiling, bestie, nothing above it. 🙅

---

💎 farey_rank

  my $rank = farey_rank(9,[5,9]);  # returns 15

Give it a positive integer n and a 2-element array reference containing a non-negative integer p and a positive integer q, and it returns the number of rationals less than p/q in the order n Farey sequence. The given fraction does not need to be an entry in the sequence, nor does it need to be in reduced form. So chill, it's flexible. 💅

A unit fraction will return the totient sum of n. Any fraction greater than one will return the length of the order n sequence, as expected.

Many OEIS sequences can be produced from this, including OEIS series A005728 OEIS series A005728 <http://oeis.org/A005728> (<= 1), OEIS series A049806 OEIS series A049806 <http://oeis.org/A049806> (<= 1/2), OEIS series A049807 OEIS series A049807 <http://oeis.org/A049807> (<= 1/3), OEIS series A049808 OEIS series A049808 <http://oeis.org/A049808> (<= 1/4), ..., OEIS series A049805 OEIS series A049805 <http://oeis.org/A049805> (<= 1/k),

---

🧮 prime_bigomega

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

Give it a non-negative integer n, and it returns Ω(|n|), the prime Omega function. This is the total number of prime factors of n including multiplicities. The result is identical to scalar(factor($n)). The return value is a read-only constant.

This corresponds to Pari's bigomega function and Mathematica's PrimeOmega[n] function. 🎯

---

🌸 prime_omega

  say "$n has ", prime_omega($n), " distinct factors";

Give it a non-negative integer n, and it returns ω(|n|), the prime omega function. This is the number of distinct prime factors of n. The result is identical to scalar(factor_exp($n)). The return value is a read-only constant.

This corresponds to Pari's omega function and Mathematica's PrimeNu[n] function. Notice it's the lowercase sibling of the one above — big Omega counts multiplicities, little omega doesn't. Sisters, not twins. 💜

---

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

Give it a single integer n, and it returns μ(|n|), the Möbius function (also known as the Moebius, Mobius, or MoebiusMu function — so many spellings, one queen). This function is 1 if n = 1, 0 if n is not square-free (i.e. n has a repeated factor), and (-1)^t if n is a product of t distinct primes. This is an important function in prime number theory. Like SAGE, we define moebius(0) = 0 for convenience. 🫶

If given two integers low and high, they define a range, and the function returns an array with the value of the Möbius function for every |n| from low to high inclusive.

⚠️ Large values of high will result in a lot of memory use. Like, a lot. Be careful bestie, your RAM has feelings.

The algorithm used for ranges is Deléglise and Rivat (1996) algorithm 4.1, which is a segmented version of Lioen and van de Lune (1994) algorithm 3.2. Math heritage. 📜

Negative ranges are possible, e.g. moebius(-30,-20) will return moebius(|n|) for -30, -29, -28, ..., -20.

The return values are read-only constants. This should almost never come up, but it means trying to modify aliased return values will cause an exception (modifying the returned scalar or array is fine). So don't try to edit the constants, that's an instant vibe check fail. 💀

---

👑 mertens

  say "Mertens(10M) = ", mertens(10_000_000);   # = 1037

Give it a non-negative integer n, and it returns M(n), the Mertens function. This function is defined as sum(moebius(1..n)), but calculated more efficiently for large inputs. For example, computing Mertens(100M) takes:

   time    approx mem
     0.01s     0.1MB   mertens(100_000_000)
     1.3s    880MB     vecsum(moebius(1,100_000_000))
    16s        0MB     $sum += moebius($_) for 1..100_000_000

Look at those numbers. 0.01s and 0.1MB versus 16s. It's not even close. 🔥

The summation of individual terms via factoring is quite expensive in time, though it uses O(1) space. Using the range version of moebius is much faster, but returns a 100M element array which, even though they are shared constants, is not good for memory at this size. In comparison, this function will generate the equivalent output via a sieving method that is relatively memory frugal and very fast. The current method is a simple n^1/2 version of Deléglise and Rivat (1996), which involves calculating all moebius values to n^1/2, which in turn will require prime sieving to n^1/4. Efficiency is the moment, and this function ate. 💅✨

okay so bestie, there's genuinely a whole buffet of algorithms for this, and the only thing they disagree on is how many μ(n) values they need to snack on 💅

The most basic vibe is just efficiently summing all n values. Then Benito and Varona (2008) pulled up with a clever, simple method that only needs n/3 values — lowkey iconic. Deléglise and Rivat (1996) describe a segmented method using only n^1/3 values. Our current implementation does a simple, non-segmented n^1/2 version of their method. Kuznetsov (2011) gives an alternate method that he says is even faster. And Helfgott and Thompson (2020) give a fast method based on advanced prime count algorithms. So many choices, all valid, we love options 🫶

---

🌸 euler_phi 🌸

say "The Euler totient of $n is ", euler_phi($n);

Hand it a single integer n and it gives you back φ(n), the Euler totient function (also called Euler's phi function or just the phi function, no big deal, she has many names 👑). It's an arithmetic function that counts how many positive integers less than or equal to n are relatively prime to n.

Because of the definition being used here, euler_phi returns 0 for all n < 1. That's the same logic SAGE uses. Meanwhile Mathematica and Pari return euler_phi(-n) for n < 0. Mathematica returns 0 for n = 0, Pari pre-2.6.2 raises an exception (drama 💀), and Pari 2.6.2 and newer returns 2.

If you call it with two integer arguments, low and high, those define an inclusive range — and you get back a list with the totient of every n from low to high, inclusive. Bulk order, slay.

✨ inverse_totient ✨

In array context, give it a non-negative integer n and it returns the complete list of values x where euler_phi(x) = n. Careful though, bestie — this can be a memory intensive operation if there are a lot of values. It's giving main character energy but also it's giving RAM usage 😭

In scalar context, it returns just the count of values. That's faster and uses substantially less memory. The whole list-vs-scalar distinction is the same situation as factor and divisors.

This roughly corresponds to the Maple function InverseTotient, and the hidden Mathematica function EulerPhiInverse. The algorithm in use here is from Max Alekseyev (2016). Respect the source 🫡

💜 jordan_totient 💜

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

Give it non-negative integers k and n, and you get Jordan's k-th totient function for n. Jordan's totient is a generalization of Euler's totient, where jordan_totient(1,$n) == euler_phi($n). It counts the number of k-tuples less than or equal to n that form a coprime tuple with n. Exactly like euler_phi, 0 is returned for all n < 1.

You can use this to generate other useful functions, like the Dedekind psi function, where psi(n) = J(2,n) / J(1,n). Little building block moment 🧚

🎀 sumtotient 🎀

Give it a non-negative integer n and it returns the summatory Euler totient function. It's defined as sum(euler_phi(1..n)), except it's calculated way more efficiently. Same answer, better attitude.

A sub-linear time recursion is implemented, using O(n^2/3) memory. Memory use is restricted so that growth becomes approximately linear above 10^13.

This is OEIS series A002088 (http://oeis.org/A002088).

🌈 ramanujan_sum 🌈

Give it two non-negative integers k and n, and you get Ramanujan's sum. This is the sum of the nth powers of the primitive k-th roots of unity. Gorgeous. Mathematical.

Heads up though: this is not related to Ramanujan summation for divergent series. Different Ramanujan, different thing, do not mix them up 🙈

⭐ exp_mangoldt ⭐

say "exp(lambda($_)) = ", exp_mangoldt($_) for 1 .. 100;

Give it a non-negative integer n and it returns EXP(Λ(n)), the exponential of the Mangoldt function (also known as von Mangoldt's function, so many aliases, very famous). The Mangoldt function equals log p if n is prime or a power of a prime, and 0 otherwise. We hand back the exponential so all the results are integers — no messy floats in the group chat.

So the return value for exp_mangoldt is:

p   if n = p^m for some prime p and integer m >= 1
1   otherwise.

🦋 liouville 🦋

Give it a non-negative integer n and it returns λ(n), the Liouville function. That's -1 raised to Ω(n), where Ω(n) is the total number of prime factors.

This corresponds to Mathematica's LiouvilleLambda[n] function. You can also compute it in Pari/GP as (-1)^bigomega(n). Options, we love options.

🌸 sumliouville 🌸

Give it a non-negative integer n and it returns L(n), the summatory Liouville function. It's defined as sum(liouville(1..n)), but calculated much more efficiently. Efficiency is the whole personality here.

There are a number of relations to the mertens function — they're basically mutuals.

This is OEIS series A002819 (http://oeis.org/A002819).

💅 chebyshev_theta 💅

say chebyshev_theta(10000);

Give it a non-negative integer n and it returns θ(n), the first Chebyshev function. That's the sum of the logarithm of each prime where p <= n. Effectively:

my $s = 0;  forprimes { $s += log($_) } $n;  return $s;

...except we compute it more efficiently and more accurately. Best of both worlds, no cap.

🔥 chebyshev_psi 🔥

say chebyshev_psi(10000);

Give it a non-negative integer n and it returns ψ(n), the second Chebyshev function. That's the sum of the logarithm of each prime power where p^k <= n for an integer k. Effectively:

my $s = 0;  for (1..$n) { $s += log(exp_mangoldt($_)) }  return $s;

...but, again, computed more efficiently and accurately. We compute it as a Neumaier sum from k = 1 .. floor(log2(n)) of chebyshev_theta(n^(1/k)). Fancy summing technique, she ate 🔥

👑 divisor_sum 👑

say "Sum of divisors of $n:", divisor_sum( $n );
say "sigma_2($n) = ", divisor_sum($n, 2);
say "Number of divisors: sigma_0($n) = ", divisor_sum($n, 0);

Give it a single non-negative integer n and it returns the sum of the divisors of n, including 1 and itself (yes, n is its own divisor, self-love 💖). We return 0 for n=0.

You can optionally pass a second non-negative integer k, which says the sum should use the k-th powers of the divisors.

This is known as the sigma function (peep Hardy and Wright section 16.7). The API is identical to Pari/GP's sigma function, and not dissimilar to Mathematica's DivisorSigma[k,n] function. Super handy for calculating things like aliquot sums, abundant numbers, perfect numbers, etc.

With various k values, the results are the OEIS sequences: OEIS series A000005 (http://oeis.org/A000005) for k=0 (number of divisors), OEIS series A000203 (http://oeis.org/A000203) for k=1 (sum of divisors), OEIS series A001157 (http://oeis.org/A001157) for k=2 (sum of squares of divisors), OEIS series A001158 (http://oeis.org/A001158) for k=3 (sum of cubes of divisors), etc. Numbered list of legends honestly.

The second argument may also be a code reference 🫶 — it gets called for each divisor and the results are summed. That lets you compute other functions, but it'll be less efficient than passing the numeric second argument. This corresponds to Pari/GP's sumdiv function.

An example of the 5th Jordan totient (OEIS A059378):

divisor_sum( $n, sub { my $d=shift; $d**5 * moebius($n/$d); } );

...though we do have a function jordan_totient which is more efficient, so like, use that one if you can bestie.

For numeric second arguments (the sigma computations), the result will be a bigint if necessary. For the code reference case, you have to take care to return bigints if overflow will be a concern. That one's on you 😭

🎀 ramanujan_tau 🎀

Give it an integer n and you get the value of Ramanujan's tau function. The result is a signed integer. Zero is returned for negative n. This corresponds to Pari v2.8's tauramanujan function and Mathematica's RamanujanTau function.

Right now this uses a simple method based on divisor sums, which — real talk — does not have a good computational growth rate. Pari's implementation uses Hurwitz class numbers and is more efficient for large inputs. Honest about its own limits, we respect that 🌸

🌈 primorial 🌈

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

Give it a non-negative integer n and it returns the primorial n#, defined as the product of the prime numbers less than or equal to n. This is the OEIS series A034386 (http://oeis.org/A034386): primorial numbers, second definition.

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

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

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

⭐ pn_primorial ⭐

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

Give it a non-negative integer n and it returns the primorial number p_n#, defined as the product of the first n prime numbers. (Compare to the factorial, which is the product of the first n natural numbers — same energy, different set 💅) This is the OEIS series A002110 (http://oeis.org/A002110): primorial numbers, first definition.

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

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

🦋 consecutive_integer_lcm 🦋

$lcm = consecutive_integer_lcm($n);

Give it a non-negative integer n and it returns the least common multiple of all integers from 1 to n. You can do this by manipulating the primes up to n, which gives way faster and more memory-friendly results than going through a factorial. Glow up achieved.

This is OEIS series A003418 (http://oeis.org/A003418). Matching that series, we define consecutive_integer_lcm(0) = 1.

🌸 partitions 🌸

Give it an integer n and it returns the partition function p(n). If n is negative, 0 is returned. This is the number of ways of writing the integer n as a sum of positive integers, with no restrictions whatsoever. Freeform. Chaotic. We love it.

This corresponds to Pari's numbpart function and Mathematica's PartitionsP function. The values produced in order are OEIS series A000041 (http://oeis.org/A000041).

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

            70    Integer::Partition
            90    MPU forpart { $n++ }
        15_000    MPU pure Perl partitions
       280_000    MPU GMP partitions
    35_000_000    Pari 2.6 numbpart
   500_000_000    Jonathan Bober's partitions_c.cc v0.6
 1_400_000_000    Pari 2.8 numbpart

I know, I know, the numbers 😭 anyway.

If you want the enumerated partitions, see forpart.

💜 lucky_numbers 💜

Give it a single 32-/64-bit non-negative integer n and it returns an array reference of values up to the input n (inclusive) which remain after the lucky number sieve originally defined by Gardiner, Lazarus, Metropolis, and Ulam. This is OEIS series A000959 (http://oeis.org/A000959). Big name energy for a sieve.

If given two non-negative integers, lo and hi, it returns sieve results between the two ranges inclusive. That's identical to the above, except it doesn't include any numbers less than lo. Currently there is very little time savings, but it does use less memory. Small win, still a win 🫶

A surprising number of the asymptotic properties of the primes are shared with this sieve, though the resulting sets are quite different. Not me finding out the primes and the lucky numbers have the same lore 😭 cute though.

okay so real talk bestie — there is no current algorithm for efficiently sieving a segment of lucky numbers 💀. Nobody's cracked that one yet. But the method used here is orders of magnitude faster than the ones linked on OEIS as of early 2023, so like, we're eating at least 🍓. CPU time growth is similar to prime sieving, about n log n. Memory use is linear with size and uses about n/25 bytes for the internal sieve. it's giving efficient, no cap.

is_lucky 🍀

Give it an integer n and it returns 1 if n is in the set of lucky numbers, and 0 if it's not. The process is analogous to trial division using the lucky numbers less than n/log(n). For inputs that don't get quickly discarded, the performance is essentially the same as generating the nth lucky number nearest to the input. iykyk.

lucky_count 🍀

Hand it a single non-negative integer n and it returns the count of lucky numbers less than or equal to n. Give it two non-negative integers lo and hi instead, and it returns the count of lucky numbers between lo and hi inclusive. both modes, one function, we love to see it ✨

lucky_count_approx 🍀

Given a single non-negative integer n, this quickly returns a good estimate of the count of lucky numbers less than or equal to n. estimate — not exact, bestie, don't come for me.

lucky_count_lower 🍀

Given a single non-negative integer n, quickly returns a lower bound of the count of lucky numbers less than or equal to n. The actual count will always be greater than or equal to the result. guaranteed, fr fr.

lucky_count_upper 🍀

Given a single non-negative integer n, quickly returns an upper bound of the count of lucky numbers less than or equal to n. The actual count will always be less than or equal to the result. so cozy, so safe 🫶.

nth_lucky 🍀

Given a non-negative integer n, returns the n-th lucky number ✨. This is done by sieving lucky numbers to n then performing a reverse calculation to determine the value at the nth position. This is much more efficient than generating all the lucky numbers up to the nth position, but it's much slower than nth_prime. 💀

the way that nth_lucky is just out here being slow compared to nth_prime... iconic behavior, honestly.

nth_lucky_approx 🍀

Given a single non-negative integer n, quickly returns a good estimate of the n-th lucky number.

nth_lucky_lower 🍀

Given a single non-negative integer n, quickly returns a lower bound of the n-th lucky number. The actual value will always be greater than or equal to the result.

nth_lucky_upper 🍀

Given a single non-negative integer n, quickly returns an upper bound of the n-th lucky number. The actual value will always be less than or equal to the result.

---

minimal_goldbach_pair 💎

Given a single non-negative integer n, returns the smallest prime p such that p + q = n and both p and q are primes. Only the single value p is returned, with q = n-p and p <= q. Both p and q are prime. two primes, one sum, zero drama ✨

undef is returned if no such p exists. This will happen for values less than 4 and for all odd n where n != 2+q for a prime q. The Goldbach Conjecture famously states that a p exists for all even n greater than 2. 🌈

This function is reasonably fast even for larger values of n as it can terminate after the first pair is found 🔥. On Macbook M1, average time is under 1 microsecond for 32-bit even inputs, under 10 microseconds for 64-bit even inputs, and 1 millisecond for 105 bit even inputs. those numbers?? she's fast.

goldbach_pair_count 💎

Given a single non-negative integer n, returns the number of pairs of primes p and q where p <= q and p + q = n.

If no such pairs exist, 0 is returned.

goldbach_pairs 💎

Given a single non-negative integer n, returns a list containing each p for all prime pairs p and q where p <= q and p + q = n. The number of elements returned is the same as goldbach_pair_count. 👑

If no such pairs exist, an empty list is returned. empty list is still a list bestie.

---

is_happy 😊

Given a single non-negative integer n, returns the number of iterations required for the map of sum of squared base-10 digits to converge to 1, or 0 if it does not converge to the value 1.

This returns the height using the OEIS A090425 definition of height, which is zero for non-happy numbers, 1 for n=1, 2 for numbers that produce 1 after a single iteration, etc. This is one more than the definitions used in many papers (e.g. Cai and Zhou 2008) where n=1 is considered to have height 0. so if you're comparing to a paper, mind the off-by-one, okay?? don't get caught lacking.

An optional base and exponent may be given (default base 10 exponent 2). The base must be between 2 and 36, and the exponent between 0 and 10. The input n is read as a decimal number, so giving input such as 1001 will be treated as the decimal 1001 regardless of base. read that twice bestie, it's a real gotcha.

With base 10 and exponent 2, non-zero values produce OEIS series A007770 <http://oeis.org/A007770>. The values themselves produce OEIS series A090425 <http://oeis.org/A090425>. 🌸

---

is_smooth 🧈

  my $is_23_smooth = is_smooth($n, 23);

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.

The values for n=0 and n=1 use the definition along with noting that factor(0) returns 0 and factor(1) returns an empty list.

The result is identical to:

  sub is_smooth { my($n,$k)=@_; return 0+(vecnone { $_ > $k } factor($n)); }

but shortcuts are taken to avoid fully factoring if possible. lazy in the productive way, periodt.

This corresponds to Mathematica's SmoothIntegerQ[n] resource function.

is_rough 🪨

  my $is_23_rough = is_rough($n, 23);

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.

The values for n=0 and n=1 use the definition along with noting that factor(0) returns 0 and factor(1) returns an empty list.

The result is identical to:

  sub is_rough { my($n,$k)=@_; return 0+(vecnone { $_ < $k } factor($n)); }

but shortcuts are taken to avoid fully factoring if possible.

okay so smooth = no big factors, rough = no small factors. literally sisters 💅. don't mix them up.

---

is_powerful 💪

  my $all_factors_cubes_or_higher = is_powerful($n, 3);

Given an integer n and an optional non-negative integer 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.

With k=2 this corresponds to Pari's ispowerful function for positive values of n. Pari chooses to define 0 as powerful and uses abs(n). we do NOT do that here — note the difference bestie.

While we can easily code this as a one line function using vecall and factor_exp, this is significantly faster and doesn't need to fully factor the input. 🔥

powerful_numbers 💪

  my $arrayref_pn1 = powerful_numbers(20);      # 1,4,8,9,16
  my $arrayref_pn2 = powerful_numbers(20,40);   # 25,27,32,36
  my $arrayref_pn3 = powerful_numbers(1,70,3);  # 1,8,16,27,32,64

Given a single non-negative integer n, returns an array reference with all 2-powerful integers from 1 to n inclusive.

Given two non-negative integers lo and hi, returns an array reference with all 2-powerful integers from lo to hi inclusive.

Given three non-negative integers lo, hi, and k, returns an array reference with all k-powerful integers from lo to hi inclusive.

  # Alternate solutions for values 1-n

  # Simple, but very slow for high $n.
  for (1..$n) { say if is_powerful($_,$k); }

  # Not so bad, especially for high $k.
  for (1..powerful_count($n,$k)) { say nth_powerful($_,$k); }

  # Best by far.
  say for @{powerful_numbers(1,$n,$k)};

Note that n <= 0 are non-powerful. zero? not powerful. negative? also not powerful. bye 🙈

powerful_count 💪

  my $npowerful3 = powerful_count(2**32, 3);

Given an integer n and an optional non-negative integer k, returns the total number of k-powerful numbers from 1 to n inclusive. If k is omitted, k=2 is used.

sumpowerful 💪

Given an integer n and an optional non-negative integer k, returns the sum of positive integer k-powerful numbers less than or equal to n. That is, the sum for all x, 1 <= x <= n, where x is a k-powerful number. If k is omitted, k=2 is used.

nth_powerful 💪

Given a non-negative integer n and an optional non-negative integer k, returns the n-th k-powerful number. If k is omitted, k=2 is used. For all k, returns undef for n=0 and 1 for n=1. yes, n=0 is undef, that's just the lore.

---

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. backwards and forwards, she does both 🦋

perfect_power_count ⭐

Given a non-negative integer n, returns the number of integers not exceeding n which are perfect powers. If given two non-negative integers lo and hi, returns the count of perfect powers between lo and hi inclusive.

By convention, numbers less than 1 are not counted.

This can be calculated extremely quickly (less than 100ns per call for native size integers), so in most cases there is no need for the approximations or bounds. girlies the approxes are RIGHT THERE but you don't need them, that's the flex 🔥

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

perfect_power_count_approx ⭐

Given a non-negative integer n, quickly returns a good estimate of the count of perfect powers less than or equal to n.

perfect_power_count_lower ⭐

Given a non-negative integer n, quickly returns a lower bound of the count of perfect powers less than or equal to n. The actual count will always be greater than or equal to the result.

perfect_power_count_upper ⭐

Given a non-negative integer n, quickly returns an upper bound of the count of perfect powers less than or equal to n. The actual count will always be less than or equal to the result.

nth_perfect_power ⭐

Given a non-negative integer n, returns the n-th perfect power. 👑

Since the perfect power count can be calculated extremely quickly, using inverse interpolation can calculate the n-th perfect power quite rapidly. fast math go brrr.

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 ✨

Bestie, hand it a non-negative integer n and it quickly gives you a good estimate of the n-th perfect power. That's it. No cap. 💅

nth_perfect_power_lower 🌸

Same vibe: give it a non-negative integer n, and it quickly returns a lower bound of the n-th perfect power. The actual value will always be greater than or equal to the result. So it's giving "floor-ish but fast" energy. 🫶

nth_perfect_power_upper 🦋

Give it a non-negative integer n, and it quickly returns an upper bound of the n-th perfect power. The actual value will always be less than or equal to the result. Iconic ceiling behavior. ✨

next_chen_prime 👑

Given a non-negative integer n, return the smallest Chen prime strictly greater than n. This will be a prime p: p > n, where p+2 is either a prime or a semiprime. Chen primes said "we're lowkey famous." 🔥

smooth_count 🍓

Given non-negative integer inputs n and k, returns the number of integers between 1 and n inclusive, that have no prime factor larger than k.

For all n, smooth_count(n,0) = smooth_count(n,1) = 1. For all k, smooth_count(0,k) = 0 and smooth_count(1,k) = 1.

This is equivalent to, but much faster than, vecsum( map { is_smooth($_,$k) } 1..$n ). So like, same result, way less waiting. 💖

rough_count 🧚

Given non-negative integer inputs n and k, returns the number of integers between 1 and n inclusive, that have no prime factor less than k.

For all n, rough_count(n,0) = rough_count(n,1) = rough_count(n,2) = n. For all k, rough_count(0,k) = 0 and rough_count(1,k) = 1.

This is equivalent to, but much faster than, vecsum( map { is_rough($_,$k) } 1..$n ). Rough but efficient. 😌

is_practical 💜

Given an integer n, returns 1 if n is a practical number, and returns 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 at http://oeis.org/A005153. Practical? More like booked and busy. 👑

carmichael_lambda 🌈

Given a non-negative integer n, returns the Carmichael function (also called the reduced totient function, or Carmichael λ(n)). This 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 at http://oeis.org/A002322.

This corresponds to Mathematica's CarmichaelLambda[n] function. It can be computed in Pari/GP as lcm(znstar(n)[2]). Math girlies, this one ate. 💅

kronecker ⭐

Given two integers a and n, returns the Kronecker symbol (a|n). The possible return values with their meanings for odd prime n are:

 0   a = 0 mod n
 1   a is a quadratic residue mod n       (x^2 = a mod n for some x)
-1   a is a quadratic non-residue mod n   (no x where x^2 = a mod n)

The Kronecker symbol is an extension of the Jacobi symbol to all integer values of n from the latter's domain of positive odd values of n. The Jacobi symbol is itself an extension of the Legendre symbol, which is only defined for odd prime values of n. This corresponds to Pari's kronecker(a,n) function, Mathematica's KroneckerSymbol[n,m] function, and GMP's mpz_kronecker(a,n), mpz_jacobi(a,n), and mpz_legendre(a,n) functions.

If n is not an odd prime, then the result does not necessarily indicate whether a is a quadratic residue mod n. Using the function is_qr will return correct results for any n, though could be slower. So don't get caught lacking by non-prime n, bestie. 😭

factorial 🎀

Given a non-negative integer n, returns the factorial of n, defined as the product of the integers 1 to n with the special case of factorial(0) = 1. This corresponds to Pari's factorial(n) and Mathematica's Factorial[n] functions. Classic. Iconic. Never not useful. ✨

subfactorial 🦋

Given a non-negative integer n, returns the subfactorial of n, which is the number of derangements of n objects. This is the number of permutations of n items where each item is not allowed to stay in its starting position.

This is OEIS series A000166 at http://oeis.org/A000166. This corresponds to Mathematica's Subfactorial[n] function. Derangements are so "not me staying where I started." 🙈

binomial 💖

Given two integers n and k, returns the binomial coefficient n*(n-1)*...*(n-k+1)/k!, also known as the choose function. Negative arguments use the Kronenburg extensions at http://arxiv.org/abs/1105.3689/. This corresponds to Pari's binomial(n,k) function, Mathematica's Binomial[n,k] function, and GMP's mpz_bin_ui function.

For negative arguments, this matches Mathematica. Pari does not implement the n < 0, k <= n extension and instead returns 0 for this case. GMP's API does not allow negative k but otherwise matches. Math::BigInt version 1.999816 and later supports negative arguments with similar semantics. Prior to this, n < 0, k > 0 was undefined.

Bestie, that "undefined" era was giving lawless behavior. Now it's at least specified. 💀

falling_factorial 🌸

Given two integers x and n, with n non-negative, returns the falling factorial of n.

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

This corresponds to Mathematica's FactorialPower[x,n] function. Falling, but make it factorial. 🍂

rising_factorial 🧚

Given two integers x and n, with n non-negative, returns the rising factorial of n.

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

This corresponds to Mathematica's Pochhammer[x,n] function. Rising like my serotonin when the code compiles. ✨

powersum 👑

  say powersum(100,1);  #     5050 = vecsum(1..100)
  say powersum(100,2);  #   338350 = vecsum(map{powint($_,2)} 1..100)
  say powersum(100,3);  # 25502500 = vecsum(map{powint($_,3)} 1..100)

Given two non-negative integers n and k, returns the sum of k-th powers of the first n positive integers.

With k=2 this is (OEIS A000330 at http://oeis.org/A000330). With k=3 this is (OEIS A000537 at http://oeis.org/A000537). With k=4 this is (OEIS A000538 at http://oeis.org/A000538). OEIS sequences can be found through k=8.

This corresponds to the faulhaber_sum(n,k) function in Math::AnyNum and Pari's dirpowerssum(n,k) function using integer arguments. Summing with main character energy. 🔥

hclassno 💜

Given an integer n, returns 12 times the Hurwitz-Kronecker class number. This will always be an integer due to the pre-multiplication by 12. The result is 0 for negative n and all n congruent to 1 or 2 mod 4.

This is related to Pari's qfbhclassno(n) where hclassno(n) for positive n equals 12 * qfbhclassno(n) in Pari/GP. This is OEIS A259825 at http://oeis.org/A259825. Multiplying by 12 so it stays integer? That's the math girlie version of setting boundaries. 💅

bernfrac 🎀

my($num,$den) = bernfrac(12);  # returns (-691,2730)

Returns the Bernoulli number B_n for an integer argument n, as a rational number represented by two integers. B_1 is chosen as 1/2, which is the same as Pari's bernfrac(n) and Mathematica's BernoulliB functions.

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

This corresponds to Pari's bernfrac function and Mathematica's BernoulliB function. Bernoulli said "install the GMP backend, bestie." 😭

bernreal 🌈

Given a non-negative integer n, returns the Bernoulli number B_n as a Math::BigFloat object using the default precision. An optional second argument may be given specifying the precision to be used.

This corresponds to Pari's bernreal function. Floaty but make it Bernoulli. 🫶

stirling ⭐

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

Given two 32-/64-bit non-negative integers n and k, plus an optional third argument kind (1, 2, or 3, with the default being 1), returns the Stirling number of the given kind. The third kind are the unsigned Lah numbers. This corresponds to Pari's stirling(n,k,{type}) function and Mathematica's StirlingS1 / StirlingS2 functions.

Stirling numbers of the first kind are -1^(n-k) times the number of permutations of n symbols with exactly k cycles. Stirling numbers of the second kind are the number of ways to partition a set of n elements into k non-empty subsets. The Lah numbers are the number of ways to split a set of n elements into k non-empty lists. Kind 1, 2, and 3? Versatile icon. 👑

fubini 🦋

Given a non-negative integer n, returns the Fubini number of n, also called the ordered Bell numbers, or the number of ordered partitions of n. It is the count of rankings of n items allowing for ties.

This is the OEIS series A000670 at http://oeis.org/A000670. Ties allowed? Inclusive and cute. 💖

harmfrac 🌸

my($num,$den) = harmfrac(12);  # returns (86021,27720)

Given a non-negative integer n, returns the Harmonic number H_n as a rational number represented by two integers. The harmonic numbers are the sum of reciprocals of the first n natural numbers: 1 + 1/2 + 1/3 + ... + 1/n.

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

This corresponds to Mathematica's HarmonicNumber function. Harmonic but not the music kind. 🎶

harmreal 🧚

Given a non-negative integer n, returns the Harmonic number H_n as a Math::BigFloat object using the default precision. An optional second integer argument may be given specifying the precision to be used.

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

Lowkey a speedrun tech: lower precision, faster compute. But only because the manpage says an asymptotic formula may be used. Not me adding extra lore. 😌

legendre_phi 💜

$phi = legendre_phi(1000000000, 41);

Given two non-negative integers n and a, returns the Legendre phi function (also called the Legendre sum). This is the count of positive integers <= n which are not divisible by any of the first a primes.

This corresponds to the legendre_phi(n,a) function in SAGE, and the --phi n a feature of primecount. Phi-nomenal. 💅

inverse_li 🔥

$approx_prime_count = inverse_li(1000000000);

Given a non-negative integer n, returns the least integer value k such that Li(k) >= n. Since the logarithmic integral Li(n) is a good approximation to the number of primes less than n, this function is a good simple approximation to the nth prime.

Approximating prime counts? She's giving quick-and-dirty but make it mathematically valid. ✨

inverse_li_nv ⚡

$faster_approx_prime_count = inverse_li_nv(1000000000);

With input x and output both in NV (floating point), computes the inverse of the logarithmic integral. This should be very fast, as everything is done in native long double precision, no Perl bigints or bigfloats are involved, and the computed result is returned as an NV.

The inverse_li function uses this to start, then ensures the integer return value is the closest inverse of the integer result of the LogarithmicIntegral function. While this is a small amount of extra time for small inputs, once we have to go to Perl and use BigInt / BigFloat, the extra time can be significant.

Floating point speed demon behavior. No bigints, no bigfloats, no waiting. 🚀

numtoperm 🎀

@p = numtoperm(10,654321);  # @p=(1,8,2,7,6,5,3,4,9,0)

Given a non-negative integer n and integer k, return the rank k lexicographic permutation of n elements. k will be interpreted as mod n!.

This will match iteration number k (zero based) of forperm. Permutation girlies, this is your moment. 🦋

This corresponds to Pari's <code>numtoperm(n,k)</code> function (and bestie, Pari 2.6 and later use the same lexicographic ordering, so you're all synced up ✨).

---

🌸 permtonum 🌸

  $k = permtonum([1,8,2,7,6,5,3,4,9,0]);  # $k = 654321

okay so you hand it an array reference 💖 that contains every integer from 0 to n-1, just in some kind of ordering, and it gives you back the lexicographic permutation rank of that set. It's literally the inverse of the numtoperm function — they're the inverse couple, your honor 💅

This will match iteration number k (zero based 👀) of forperm, and the result is gonna land somewhere between 0 and n!-1.

This corresponds to Pari's permtonum(n) function (and yes, Pari 2.6 and later use the same lexicographic ordering 🫶).

🦋 randperm 🦋

  @p = randperm(100);   # returns shuffled 0..99
  @p = randperm(100,4)  # returns 4 elements from shuffled 0..99
  @s = @data[randperm(1+$#data)];    # shuffle an array
  @p = @data[randperm(1+$#data,2)];  # pick 2 from an array

Give it a single non-negative integer n, and it hands you a random permutation of the integers from 0 to n-1. iconic behavior honestly 🌈

Optionally you can slide in a second non-negative integer argument k, and then the returned list will only have k elements in it. and listen — this is way more efficient than just truncating the whole full shuffled list, no cap ✨

The randomness comes from our CSPRNG 🔥

The slicing technique shown in that last example is similar to vecsample.

💜 shuffle 💜

  @shuffled = shuffle(@data);

Takes a list as input and returns a random permutation of that list. Like randperm, the randomness comes from our CSPRNG.

This function is functionally identical to the shuffle function in List::Util. The only difference? The random source (Chacha20, with better randomness, a larger period, and a larger state 🎀). Which does make it slower, ngl.

If you want the entire shuffled array, this is faster than slicing with randperm like the example above shows. But if you only need fewer elements (a "pick" or a "sample" 🍓), then vecsample or slicing with randperm will be much more efficient.

🧚 vecsample 🧚

  $oneof = vecsample(1,@data);  # Select one random value
  @twoof = vecsample(2,@data);  # Select two random values

Takes a non-negative integer k plus a list, and returns k randomly selected values from that list. The randomness comes from our CSPRNG 💖

Here's the cute part ✨ if the input is exactly two elements (k and one other) and the second value is an array reference, then we'll use it as the input list:

  $oneof = vecsample(1, $arrayref);
  @twoof = vecsample(1, \@data);

This can be a big performance increase if your input list is large (like 2x at 1000 elements, and it can be 10x with more 👑). Yeah, there might be a little confusion when you're sampling a list with exactly one element, where that element is an array reference, but lowkey this is assumed to be a rare case.

This is similar to sample from List::Util, choose_multiple from Rust rand, and Raku's pick.

---

🌈 MODULAR ARITHMETIC 🌈

✨ OVERVIEW ✨

Functions for fast modular arithmetic are provided, and yes bestie, the whole lineup: add, subtract, multiply, divide, power, square root, nth root, inverse. Plus, fast modular calculation of factorial, binomial, and Lucas sequences are provided too 🎀 See "MODULAR FUNCTIONS" for more functions that operate mod n.

Semantics mostly follow Pari/GP, though in some cases they will indicate an error while we just return undef 😭

• We use the absolute value of the modulus.
• A modulus of zero returns undef.
• A modulus of 1 will return 0.
• If a modular result doesn't exist, we return undef.

💅 negmod 💅

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

This is similar to submod(0,$a,$n) or $n ? modint(-$a,absint($n)) : undef.

🌸 addmod 🌸

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

🦋 submod 🦋

Given three integers a, b, and n, return (a-b) mod |n|.

🌈 mulmod 🌈

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

🎀 muladdmod 🎀

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

💖 mulsubmod 💖

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

🍓 divmod 🍓

Given three integers a, b, and n, return (a/b) mod |n|. This is done as (a * (1/b mod |n|)) mod |n|. And if no inverse of b mod |n| exists? then undef is returned 🙈

👑 powmod 👑

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

powmod(a,-b,n) is calculated as powmod(invmod(a,n),b,n). If 1/a mod |n| does not exist, undef is returned.

✨ sqrtmod ✨

Given two integers a and n, return the square root of a mod |n|. If no square root exists, undef is returned. If defined, the return value r will always satisfy r^2 = a mod |n| 🫶

If the modulus is prime, the function will always return r, the smaller of the two square roots (the other one being -r mod |n| 💅). If the modulus is composite, one of possibly many square roots will be returned, and it will not necessarily be the smallest. just a heads up, no cap.

🔥 allsqrtmod 🔥

Given two integers a and n, returns a sorted list of all modular square roots of a mod |n|. If no square root exists, an empty list is returned.

Some inputs will return very many roots 😭 like for example, a = p^4, n = 24 * p^4 for prime p, has many roots, and sqrtmod(89**8, 24*89**8) has over 500 million. yeah. over five hundred million. bestie that's a lot.

💜 rootmod 💜

Given three integers a, k, and n, returns a k-th root of a modulo |n|, or undef if one does not exist. If defined, the return value r will satisfy r^k = a mod |n|. There is no guarantee that the smallest root will be returned though 🙈

For some composites with large prime powers this may not be efficient.

rootmod(a,-k,n) is calculated as rootmod(invmod(a,n),k,n). If 1/a mod |n| does not exist, undef is returned.

🌸 allrootmod 🌸

Given three integers a, k, and n, returns a sorted list of all modular k-th root of a modulo |n|. If no root exists, an empty list is returned.

Similar to allsqrtmod, some inputs have millions or billions of roots, so it might not be able to successfully return them all 💀

🧚 invmod 🧚

  say "The inverse of 42 mod 2017 = ", invmod(42,2017);

Given two integers a and n, return the inverse of a modulo |n|. If not defined, undef is returned. If defined, then the return value multiplied by a equals 1 modulo |n| ✨

The results correspond to the Pari result of lift(Mod(1/a,n)). The semantics with respect to negative arguments match Pari. Notably, a negative n is negated, which is different from Math::BigInt, but in both cases the return value is still congruent to 1 modulo n as expected 🫶

Mathematica uses Powermod[a, -1, n], where n must be positive.

🎀 factorialmod 🎀

Given a non-negative integer n and an integer m, returns n! mod |m|. This is much faster than computing the large factorial(n) followed by a mod operation 👑

While very efficient, this is not state of the art. Currently, Fredrik Johansson's fast multi-point polynomial evaluation method as used in FLINT is the fastest known implementation. This becomes noticeable for n > 10^8 or so, and the O(n^.5) versus O(n) complexity is very apparent with large n 😭

🌈 binomialmod 🌈

Given integer arguments n, k, and m, returns binomial(n,k) mod |m|. This is much faster than computing the large binomial(n,k) followed by a mod operation 🔥

|m| does not need to be prime. The result is extended to negative n. Negative k will return zero.

This corresponds to Mathematica's BinomialMod[n,m,p] function. It has similar functionality to Max Alekseyev's binomod.gp Pari routine.

💖 lucasumod 💖

Given integers P, Q, the non-negative integer k, and the integer n, efficiently compute lucasu(P,Q,k) mod |n|.

This corresponds to gmpy2's lucasu_mod function.

When (P,Q) = (1,-1) this returns the modular Fibonacci sequence. This corresponds to Sidef's fibmod function ✨

🦋 lucasvmod 🦋

Given integers P, Q, the non-negative integer k, and the integer n, efficiently compute lucasv(P,Q,k) mod |n|.

This corresponds to gmpy2's lucasv_mod function.

⭐ lucasuvmod ⭐

  # Compute the 5000-th Fibonacci and Lucas numbers, mod 1001
  ($U,$V) = lucasuvmod(1, -1, 5000, 1001);

Given integers P, Q, the non-negative integer k, and the integer n, efficiently compute the k-th value of U(P,Q) mod |n| and V(P,Q) mod |n| 💅

This is similar to the lucas_sequence function, but uses a more consistent argument order and does not return Q_k.

🍓 lucas_sequence 🍓

  my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k);

🚨 lucas_sequence() is deprecated. Use lucasuvmod() instead, bestie, we're moving on 💅

Computes U_k, V_k, and Q_k for the Lucas sequence defined by P,Q, modulo |n|. The modular Lucas sequence is used in a number of primality tests and proofs. k must be non-negative, and n must be non-zero.

🌸 pisano_period 🌸

Given a non-negative integer n, returns the period of the Fibonacci sequence modulo n. The modular Fibonacci numbers can be produced using lucasumod(1,-1,k,n). They are periodic for any integer n, and the Pisano period is the length of the repeating sequence ✨

This is the OEIS series A001175 (http://oeis.org/A001175) 🫶

💖 MODULAR FUNCTIONS 💖

OVERVIEW

okay so bestie, there are more functions here and they all work mod n 🌸. They keep the same vibe as the others when it comes to the modulus: the absolute value gets used, and if your modulus is 0 you just get undef back, sorry babe 😭.

⚠️ but heads up!! the behavior with n = 1 is not always the same. it's giving "read each one individually", not "assume they all match" 💅

---

🔓 znlog

$k = znlog($a, $g, $p)

Returns the integer k that solves the equation a = g^k mod |p|, or undef if there's no solution to be found 🕵️♀️. This is literally the discrete logarithm problem, the icon itself 👑

The implementation for native integers first applies Silver-Pohlig-Hellman on the group order to possibly reduce the problem to a set of smaller problems ✨. Then the solutions get performed using a mixture of trial, Shanks' BSGS, and Pollard's DLP Rho — a whole little squad fr fr 🫶

The PP implementation is less sophisticated, with only a memory-heavy BSGS being used. ngl it's giving "budget version" but it still shows up 💀

---

🔢 znorder

$order = znorder(2, next_prime(10**16)-6);

Given two positive integers a and n, this returns the multiplicative order of a modulo |n| 🌸. Translation: it's the smallest positive integer k such that a^k ≡ 1 mod |n|.

It returns undef if n = 0, if a = 0, or if a and n are not coprime — because no value can result in 1 mod n, so there's just nothing to give you bestie 😭. And it returns 1 if a = 1 or if n = 1 ✨.

🚨 note that the n = 1 case differs from other mod functions, because the return value is a positive integer, not an integer mod n. like, it did NOT get the memo 💅

This corresponds to Pari's znorder(Mod(a,n)) function and Mathematica's MultiplicativeOrder[a,n] function. very "same energy, different font" of them 🎀

---

🌱 znprimroot

Given an integer n, where n is treated as |n|, this returns the smallest primitive root of (Z/nZ)^*, or undef if no root exists 🙈.

A root exists when euler_phi($n) == carmichael_lambda($n), which will be true only if n is one of {2, 4, p^k, 2p^k} for odd prime p. very specific guest list, no plus-ones 👑

Like other modular functions, if n = 0 the function returns undef. consistent queen 💅

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

---

🌈 is_primitive_root

Given two integers a and n, returns 1 if a is a primitive root modulo |n|, and 0 if not. So it's giving yes/no energy 🌸

If a is a primitive root, then euler_phi(n) is the smallest e for which a^e = 1 mod n ✨

Like other modular functions, if n = 0 the function returns undef. yes bestie, same rule as always 🫶

---

🍓 qnr

Given an integer n, this returns the least quadratic non-residue modulo |n| 💜. That means the smallest integer a where there does not exist an integer b such that a = b^2 mod |n| 🙅♀️

Like other modular functions, if n = 0 the function returns undef. you already knew that though, main character consistency 👑

This is OEIS A020649 <http://oeis.org/A020649>. For primes it is OEIS A053760 <http://oeis.org/A053760> 🎀

---

✨ is_qr

Given two integers a and n, returns 1 if a is a quadratic residue</quadratic> modulo |n|, and 0 otherwise.

A return value of 1 indicates there exists an x where a = x^2 mod |n| 🔥 so like, the square root situation is solvable, slay 💅

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

For all values, this will be equal to sqrtmod(a,n) != undef, with possibly better performance ⚡ so if you're just checking existence and not actually needing the root, this is your girl fr fr

Like other modular functions, if n = 0 the function returns undef. 🌸 literally every single one of these, i fear 💀

OVERVIEW

okay so bestie, storytime, because the RNG lore here is genuinely important 🎀

before version 5.20, Perl's rand was just calling your system's rand function. which meant it acted different on every machine, and was, like, almost always a bad choice. for 5.20, Perl standardized on drand48 and shipped the source with it so there'd be zero system dependencies. improvement? sure. good PRNG? not even a little. it's only got 32 bits of actual random values and it flops a bunch of statistical tests. see http://www.pcg-random.org/statistical-tests.html if you wanna see it get dragged. 💀

there are so many better choices for a standard random number generator — Mersenne Twister, PCG, Xoroshiro128+ — and someday maybe Perl will get one of those to replace drand48. until then: Math::Random::MTwist has a ton of features and excellent performance, or, hi, this module 💅

and since we're often out here dealing with random primes for cryptographic purposes, we have ✨extra requirements✨. this module uses a CSPRNG for its random stream. specifically ChaCha20, the same algorithm BSD's arc4random uses, and what /dev/urandom uses on BSD and Linux 4.8+. seeding happens at startup using the Win32 Crypto API (on Windows), /dev/urandom, /dev/random, or Crypt::PRNG — whichever one it finds first, bestie. 🌸

we use the original ChaCha definition, not RFC7539. which means a 64-bit counter, so a period of 2^72 bytes, or 2^68 calls to drand or irand64. compare that to the cute little 2^48 period of Perl's drand48. ✨ we've got a 512-bit state, which is significantly bigger than the 48-bit drand48 state. when seeding, 320 bits (40 bytes) get used. and among other things, this means all 52! permutations of a shuffled card deck are actually possible — which is not true of shuffle in List::Util. 👑

and no, do not come for me about performance, because one might think a CSPRNG would be slow, but benchmarking says no. the speed of irand, irand64, and drand is within 20% of the fastest existing modules using non-CSPRNG methods, and 2 to 20 times faster than most. and honestly? while a faster underlying RNG is useful, the Perl call interface overhead is the majority of the time for these calls anyway. carefully tuning that interface is critical for any module. periodt.

for performance on large amounts of data, peep the tables down in random_bytes.

each thread gets its own context, so seeding in one thread has literally no impact on the other threads. besides being more secure, this is better for performance than one single context with locks. if you need explicit control of multiple independent streams, then bestie, use a more specific module — I believe Crypt::PRNG (part of CryptX) and Bytes::Random::Secure are good alternatives. 🫶

and using the :rand export option will define rand and srand as similar but improved versions of the system functions of the same name, plus irand and irand64. iconic.

---

irand

$n32 = irand;     # random 32-bit integer

returns a random 32-bit integer, using the CSPRNG. slay. 🎲

irand64

$n64 = irand64;   # random 64-bit integer

returns a random 64-bit integer using the CSPRNG (on 64-bit Perl). on a 32-bit Perl, it gives you the maximum UV bits, which will only be 32. 🙈

drand

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

returns a random NV (Perl's native floating point) using the CSPRNG. the API is similar to Perl's rand but gives better results, obviously. ✨

the number of actual random bits will be equal to the number of mantissa bits in the NV type. for IEEE-754 doubles that means 53 bits, and it can go to 64 or 113 bits with long double / quadmath support. the _nvmantbits function lets you see how many bits are used. 🔍

this gives substantially better quality random numbers than the default Perl rand. among other things, on modern Perls, rand uses drand48, which has 32 bits of not-very-good randomness plus 16 more bits of obvious patterns (like, the 48th bit alternates, the 47th has a period of 4, etc. 💀 ). output from rand fails at least 5 tests from the TestU01 SmallCrush suite, while our function easily passes. ate and left no crumbs.

with the :rand tag, this function is additionally exported as rand.

random_bytes

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

given an unsigned number n of bytes, returns a string filled with random data from the CSPRNG. performance for large quantities, bestie:

Module/Method                  Rate   Type
-------------             ---------   ----------------------

Math::Prime::Util::GMP    1067 MB/s   CSPRNG - ISAAC
ntheory random_bytes       384 MB/s   CSPRNG - ChaCha20
Crypt::PRNG                140 MB/s   CSPRNG - Fortuna
Crypt::OpenSSL::Random      32 MB/s   CSPRNG - SHA1 counter
Math::Random::ISAAC::XS     15 MB/s   CSPRNG - ISAAC
ntheory entropy_bytes       13 MB/s   CSPRNG - /dev/urandom
Crypt::Random               12 MB/s   CSPRNG - /dev/urandom
Crypt::Urandom              12 MB/s   CSPRNG - /dev/urandom
Bytes::Random::Secure        6 MB/s   CSPRNG - ISAAC
ntheory pure perl ISAAC      5 MB/s   CSPRNG - ISAAC (no XS)
Math::Random::ISAAC::PP      2.5 MB/s CSPRNG - ISAAC (no XS)
ntheory pure perl ChaCha     1.0 MB/s CSPRNG - ChaCha20 (no XS)
Data::Entropy::Algorithms    0.5 MB/s CSPRNG - AES-CTR

Math::Random::MTwist       927 MB/s   PRNG - Mersenne Twister
Bytes::Random::XS          109 MB/s   PRNG - drand48
pack CORE::rand             25 MB/s   PRNG - drand48 (no XS)
Bytes::Random                2.6 MB/s PRNG - drand48 (no XS)

entropy_bytes

similar to random_bytes, but it goes directly to the entropy source. not normally recommended, because it can consume shared system resources and it's typically slow — on the computer that produced that random_bytes chart up there, using dd generated the same 13 MB/s performance as our entropy_bytes function.

the actual performance will be highly system dependent. 🌈

urandomb

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

give it a number of bits b, and it returns a random unsigned integer less than 2^b. the result will be uniformly distributed between 0 and 2^b-1 inclusive. no bias, we don't do that here. 💅

urandomm

$n = urandomm(100);    # random integer in [0,99]
$n = urandomm(1024);   # random integer in [0,1023]

given a positive integer n, returns a random unsigned integer less than n. the results will be uniformly distributed between 0 and n-1 inclusive. care is taken to prevent modulo bias. 🎀

csrand

takes a binary string data as input and seeds the internal CSPRNG. this is not normally needed, because system entropy is used as a seed on startup. for best security this should be 16-128 bytes of good entropy. no more than 1024 bytes will be used.

with no argument, it reseeds using system entropy, which is preferred. ✨

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

srand

takes a single UV argument and seeds the CSPRNG with it, and also returns it. if no argument is given, a new UV seed is constructed. heads up though bestie: this creates a very weak seed from a cryptographic standpoint, so it's useful for testing or simulations, but csrand is recommended, or just keep using the system entropy default seed. 🫶

the API is nearly identical to the system function srand. it uses a UV which can be 64-bit rather than always 32-bit. the behaviour for undef, empty string, empty list, etc. is slightly different (we treat these as 0).

this function is not exported with the :all tag, but is with :rand.

if the secure configuration has been set, this function will croak. manual seeding using srand is not compatible with cryptographic security. no cap.

rand

an alias for drand, not exported unless the :rand tag is used.

random_factored_integer

my($n, $factors) = random_factored_integer(1000000);

given a positive non-zero input n, returns a uniform random integer in the range 1 to n, along with an array reference containing the factors. 🌸

this uses Kalai's algorithm for generating random integers along with their factorization, and it is much faster than the naive method of generating random integers followed by a factorization. a later implementation may use Bach's more efficient algorithm. iykyk. 💜

✨ RANDOM PRIMES ✨

💖 random_prime

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

okay so this one hands you back a pseudo-randomly picked prime that's greater than or equal to the lower limit and less than or equal to the upper limit. no lower limit given? then 2 is implied, bestie. and if there are genuinely no primes living in that range, you get undef back. sad but honest. 😭

the whole goal here is a uniform distribution of the primes in the range — meaning every single prime in that range has the exact same chance of showing up. this is why it throws algorithms like "PRIMEINC" out of consideration entirely: sure, they're efficient, but their output is very non-random and we are simply not doing that. it also means the numbers themselves won't be evenly distributed, because primes are literally not evenly distributed, that's just their vibe. put differently: these random prime functions return a uniformly selected prime out of the set of primes in the range. so with random_prime(1000), the numbers 2, 3, 487, 631, and 997 all have the same probability of being the one. iconic behavior honestly 🎀

how it actually rolls depends on size:

• small numbers — a random index selection is done, which gives ideal uniformity and is super efficient with small inputs.
• ranges bigger than this ~16-bit threshold but still within the native bit size — a Monte Carlo method gets used. also ideal uniformity, and it can be really fast for reasonably sized ranges. ✨
• even larger numbers — we partition the range, pick a random partition, then select a random prime from inside that partition. this gives some loss of uniformity but consumes many fewer bits of randomness and is much faster.

---

💅 random_ndigit_prime

say "My 4-digit prime number is: ", random_ndigit_prime(4);

selects a random n-digit prime, where the input is an integer number of digits. one of the primes in that range (like 1000 - 9999 for 4-digits) gets uniformly selected. periodt.

if the number of digits is greater than or equal to the maximum native type, then the result comes back as a BigInt. however, if the "nobigint" configuration option is on, then output gets restricted to native size numbers, and asking for more digits than natively supported results in an error. heads up 💀

for better performance with large bit sizes, install Math::Prime::Util::GMP. go on, do it.

---

🦋 random_nbit_prime

my $bigprime = random_nbit_prime(512);

selects a random n-bit prime, where the input is an integer number of bits. a prime with the nth bit set will be uniformly selected.

for bit sizes of 64 and lower, random_prime is used, which gives completely uniform results in this range. for sizes larger than 64, Algorithm 1 of Fouque and Tibouchi (2011) is used — we select a random odd number for the lower bits, then loop selecting random upper bits until the result is prime. this allows a more uniform distribution than the general random_prime case while running slightly faster. (in contrast, for large bit sizes random_prime selects a random upper partition then loops on the values within the partition, which very slightly skews the results towards smaller numbers. just so you know. 👀)

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

---

🧚 random_safe_prime

my $bigprime = random_safe_prime(512);

produces an n-bit safe prime. that's a prime p where p = 2q+1 and q is also prime. the way that's such a cute little relationship 😭

these types of primes are sometimes useful for discrete logarithm based cryptography, and can be generated more efficiently using simultaneous sieving.

---

🔥 random_strong_prime

my $bigprime = random_strong_prime(512);

constructs an n-bit strong prime using Gordon's algorithm. we consider a strong prime p to be one where:

• p is large. this function requires at least 128 bits.
• p-1 has a large prime factor r.
• p+1 has a large prime factor s
• r-1 has a large prime factor t

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

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

---

👑 random_proven_prime

my $bigprime = random_proven_prime(512);

constructs an n-bit random proven prime. internally this may use is_provable_prime(random_nbit_prime) or random_maurer_prime depending on the platform and bit size. it's giving adaptive, we love to see it.

---

💜 random_proven_prime_with_cert

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

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

---

⭐ random_maurer_prime

my $bigprime = random_maurer_prime(512);

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

the performance with Math::Prime::Util::GMP installed is hundreds of times faster, so it is highly recommended. like, actually please install it.

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

internally this additionally runs the BPSW probable prime test on every partial result, and constructs a primality certificate for the final result, which is verified. these provide additional checks that the resulting value has been properly constructed. double-checking the homework 💅

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

---

🌸 random_maurer_prime_with_cert

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

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

---

🌈 random_shawe_taylor_prime

my $bigprime = random_shawe_taylor_prime(8192);

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

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

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

---

🎀 random_shawe_taylor_prime_with_cert

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

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

---

🍓 random_semiprime

takes a positive integer number of bits bits, returns a random semiprime of exactly bits bits. the result has exactly two prime factors (hence semiprime). no cap.

the factors will be approximately equal size, which is typical for cryptographic use. for example, a 64-bit semiprime of this type is the product of two 32-bit primes. bits must be 4 or greater.

some effort is taken to select uniformly from the universe of bits-bit semiprimes. this takes slightly longer than some methods that do not select uniformly. worth it though 💖

---

💀 random_unrestricted_semiprime

takes a positive integer number of bits bits, returns a random semiprime of exactly bits bits. the result has exactly two prime factors (hence semiprime).

the factors are uniformly selected from the universe of all bits-bit semiprimes. this means semiprimes with one factor equal to 2 will be most common, 3 next most common, etc. bits must be 3 or greater. (yeah, 3 instead of 4 — that's not a typo, that's just the ✨unrestricted✨ life.)

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

✨ UTILITY FUNCTIONS ✨

okay bestie, these are the lil' housekeeping functions — the ones that tell the module how to organise its brain and when to clean up after itself. not the flashiest part of the page but honestly? iconic. 💅

prime_precalc

prime_precalc( 1_000_000_000 );

This lets the module prep itself for fast operation up to a specific number. Like, you're giving it a heads-up: "hey, we're gonna be working with numbers this big, get ready." 🏃♀️

Not necessary to call — genuinely optional — but it gives you more control over when memory gets allocated, and in some cases it gives faster results when you're making multiple calls. In the current implementation this will calculate a sieve for all numbers up to the number you specify.

prime_memfree

prime_memfree;

Frees any extra memory the module may have allocated. Just like with prime_precalc, you don't have to call it, but if you're done making calls, or you just want things tidied up, you can use this. 🧹✨

The object method might be a better choice for complicated uses — keep reading, bestie.

Math::Prime::Util::MemFree->new

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

This is the more robust way of making sure any cached memory gets freed, because it's handled by the last MemFree object leaving scope. so like — the moment $mf goes poof, the memory gets recovered. 🦋

Why does that matter? Because if your routines were chillin' inside an eval that died, things will still get cleaned up. no leaks, no drama. And if you call another function that uses a MemFree object, the cache will stay in place because you still have an object. main character memory management fr fr. 👑

prime_get_config

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

# Print all configuration
my $r=prime_get_config();
say "$_ $r->{$_}" for sort (keys %$r);

Returns a reference to a hash of the current settings. Heads up: the hash is a copy of the configuration, so changing it has no effect. it's giving "look but don't touch". 🙈

The settings include:

• verbose — verbose level. 1 or more will result in extra output.
• bigintclass — the bigint type name (default Math::BigInt).
• precalc_to — primes up to this number are calculated.
• maxbits — the maximum number of bits for native operations.
• xs — 0 or 1, indicating the XS code is available.
• gmp — 0 or 1, indicating GMP code is available.
• maxparam — the largest value for most functions, without bigint.
• maxdigits — the max digits in a number, without bigint.
• maxprime — the largest representable prime, without bigint.
• maxprimeidx — the index of maxprime, without bigint.
• assume_rh — whether to assume the Riemann hypothesis (default 0).
• secure — disable ability to manually seed the CSPRNG.

prime_set_config

prime_set_config( assume_rh => 1 );

prime_set_config(bigint=>Math::GMPz);

Lets you set some of those parameters. Currently the options are:

• verbose — The default setting of 0 will generate no extra output. Setting to 1 or higher results in extra output. For example, at setting 1 the AKS algorithm will indicate the chosen r and s values. At setting 2 it will output a sequence of dots indicating progress. Similarly, for random_maurer_prime, setting 3 shows real time progress. Factoring large numbers is another place where verbose settings can give progress indications.
• bigint — You can give either a single object (e.g. a value of the class you want), or a comma separated list of class names. The first class we can load will be used for all operations that use a bigint. A warning will be produced if one was not found.
• trybigint — Exactly the same behavior as bigint but no warning will be output if we couldn't load anything from the list.
• xs — Allows turning off the XS code, forcing the Pure Perl code to be used. Set to 0 to disable XS, set to 1 to re-enable. You probably will never want to do this.
• gmp — Allows turning off the use of Math::Prime::Util::GMP, which means using Pure Perl code for big numbers. Set to 0 to disable GMP, set to 1 to re-enable. You probably will never want to do this.
• assume_rh — Allows functions to assume the Riemann hypothesis is true if set to 1. This defaults to 0. Currently this setting only impacts prime count lower and upper bounds, but could later be applied to other areas such as primality testing. A later version may also have a way to indicate whether no RH, RH, GRH, or ERH is to be assumed.
• secure — The CSPRNG may no longer be manually seeded. Once set, this option cannot be disabled. srand will croak if called, and csrand will croak if called with any arguments. csrand with no arguments is still allowed, as that will use system entropy without giving anything to the caller. The point of this option is to ensure that any called functions do not try to control the RNG.

🚨 bestie listen — that secure one is a one-way door. once you flip it you literally cannot unflip it, so be sure before you commit. the whole vibe of it is locking the RNG down so nothing downstream can meddle with it. periodt. 🔒

🔥 FACTORING FUNCTIONS 🔥

The moment you've all been waiting for: breaking big numbers into their lil prime pieces bestie 💅✨

---

factor</p>

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

Okay so this bad boy produces the prime factors of a positive number input, and it gives them to you in numerical order — the product of everything it returns will equal your input, always. Periodt.

• n = 1 → returns an empty list 🫶
• n = 0 → returns 0

(this matches Pari, iykyk)

In scalar context, it returns Ω(n), which is the total number of prime factors (that's OEIS A001222). This lines up with Pari's bigomega(n) function and Mathematica's PrimeOmega[n] function, and it's literally the same result you'd get by evaluating the resulting array in scalar context.

The algorithm currently does a little trial division, checks for perfect powers, then runs combinations of Pollard's Rho, SQUFOF, and Pollard's p-1 — and that combo gets applied to every non-prime factor it finds.

Factoring bigints works with pure Perl, which is honestly super handy on 32-bit machines for numbers just over the 32-bit limit 🌸. But heads up: it can be very slow for "hard" numbers 😭. Installing the Math::Prime::Util::GMP module speeds up bigint factoring a lot, and all future effort on large number factoring will live in that module. If you don't have it for some reason, use the GMP or Pari version of bigint if you can (like use bigint try => 'GMP,Pari') — that runs 2-3x faster (though still 100x slower than the real GMP code 💀).

---

factor_exp</p>

my @factor_exponent_pairs = factor_exp(29513484000);
# returns ([2,5], [3,4], [5,3], [7,2], [11,1], [13,2])
# factor(29513484000)
# returns (2,2,2,2,2,3,3,3,3,5,5,5,7,7,11,13,13)

This one gives you pairs of prime factors and their exponents, in numerical factor order. Some algorithms find this form way more convenient, and honestly fair. It's the same shape that Mathematica's FactorInteger[n] and Pari/GP's factorint functions return. Small note: Math::Pari transposes the Pari result matrix, so watch out for that 🙈.

In scalar context, it returns ω(n), the number of unique prime factors (OEIS A001221). That corresponds to Pari's omega(n) and Mathematica's PrimeNu[n], and yep — same result as evaluating the array in scalar context.

The internals are identical to factor, so all those comments up there apply. Only the arrangement of the factors is different. Simple as that ✨

---

divisors</p>

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 💖. The divisors are a power set of multiplications of the prime factors, returned as a uniqued sorted list. The result is identical to Pari's divisors and Mathematica's Divisors[n] functions.

In scalar context this returns the sigma0 function (see Hardy and Wright section 16.7). That's OEIS A000005. Same result as evaluating the array in scalar context, but more efficient. This corresponds to Pari's numdiv and Mathematica's DivisorSigma[0,n].

Also peep the fordivisors function if you wanna loop over the divisors instead 👑

When n=0 we return the empty set (which is zero in scalar context).

An optional second positive integer argument k means the results will not include any value larger than k. This is especially useful when your number has thousands of divisors and you only care about the small ones. So clutch 🫶

---

trial_factor</p>

my @factors = trial_factor($n);

Produces the prime factors of a positive number input using trial division. The factors come in numerical order. For large inputs this will be very slow 😭 (it's giving brute force, bestie).

An optional second argument sets an upper limit for factors. Factors 2, 3, and 5 are always pulled out. Factors larger than the second argument will not be found — which means the last value in the list might be composite 👀.

Like all the specific-algorithm *_factor routines, this is not exported unless you explicitly request it.

---

fermat_factor</p>

my @factors = fermat_factor($n);

Produces factors — not necessarily prime — of the positive number input. The algorithm here is Knuth's algorithm C. For small inputs this is very fast ⚡, but it slows down quite rapidly as the number of digits goes up. It absolutely eats for inputs with a factor close to the midpoint (think a semiprime p*q where p and q have the same number of digits). Iconic behavior honestly.

---

holf_factor</p>

my @factors = holf_factor($n);

Produces factors — not necessarily prime — of the positive number input. You can optionally pass a number of rounds as the second parameter. It's possible the function just... won't find a factor 🙈, in which case a single element (the input) is returned.

This uses Hart's One Line Factorization with no premultiplier. It's an interesting alternative to Fermat's algorithm, and there are some inputs it can rapidly factor 🌈. Overall, though, it has the same advantages and disadvantages as Fermat's method.

---

lehman_factor</p>

my @factors = lehman_factor($n);

Produces factors — not necessarily prime — of the positive number input. An optional argument, defaulting to 0 (false), says whether to run trial division. Without trial division, it's possible the function will be unable to find a factor — in which case a single element, the input, is returned.

This is Warren D. Smith's Lehman core with minor modifications. It's limited to 42-bit inputs: n < 8796393022208 — so don't come at it with anything bigger, it said what it said 💅

---

squfof_factor</p>

my @factors = squfof_factor($n);

Produces factors — not necessarily prime — of the positive number input. Optional number of rounds goes in as the second parameter. Again, it's possible it won't find a factor 🙈, in which case a single element, the input, is returned. This function typically runs very fast ⚡🔥.

---

prho_factor</p>

Pollard's rho factoring algorithm. See pbrent_factor right below for the shared description of both functions — they're basically twins 👯

---

pbrent_factor</p>

my @factors = prho_factor($n);
my @factors = pbrent_factor($n);

# Use a very small number of rounds
my @factors = prho_factor($n, 1000);

Produces factors — not necessarily prime — of the positive number input. An optional number of rounds goes in as the second parameter. These two attempt to find a single factor using Pollard's Rho algorithm — either the original version or Brent's modified version.

These are more specialized, usually used for pre-factoring very large inputs, because they are very fast at snagging small factors. Main character energy at finding the little guys 🔍✨

---

pminus1_factor</p>

my @factors = pminus1_factor($n);
my @factors = pminus1_factor($n, 1_000);          # set B1 smoothness
my @factors = pminus1_factor($n, 1_000, 50_000);  # set B1 and B2

Produces factors — not necessarily prime — of the positive number input. This is Pollard's p-1 method, using two stages with default smoothness settings of 1_000_000 for B1 and 10 * B1 for B2.

This method can rapidly find a factor p of n where p-1 is smooth (meaning it has no large factors). So when the vibes are smooth, this one delivers 💖

---

pplus1_factor</p>

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

Produces factors — not necessarily prime — of the positive number input. This is Williams' p+1 method, using one stage and two predefined initial points. Simple, focused, does its thing ✨

---

cheb_factor</p>

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

Produces factors — not necessarily prime — of the positive number input. This one uses the properties of Chebyshev polynomials — particularly the fact that T_mn(x) = T_m(T_n(x)) — and their relationship with the Lucas sequence, to find factors if p-1 or p+1 is smooth.

It generally works better than our pplus1_factor, but is slower than our pminus1_factor. Trade-offs, bestie 🤷

---

ecm_factor</p>

my @factors = ecm_factor($n);
my @factors = ecm_factor($n, 100, 400, 10);      # B1, B2, # of curves

Produces factors — not necessarily prime — of the positive number input. This is the elliptic curve method using two stages 🌀💜

✨ MATHEMATICAL FUNCTIONS ✨

ExponentialIntegral

my $Ei = ExponentialIntegral($x);

Okay so bestie, hand it a non-zero floating point x and it gives you back the real-valued exponential integral of x — which is defined as the integral of e^t/t dt going from -infinity up to x. That's the whole vibe. 💜

For non-BigFloat inputs, the result should be accurate to at least 14 digits. No cap.

For BigFloat inputs, you only get full accuracy and full performance if Math::Prime::Util::GMP is installed. If that module isn't around, other methods step in and still give you at least 14 digits of accuracy:

• continued fractions (when x < -1)
• rational Chebyshev approximation (when -1 < x < 0)
• a convergent series (for small positive x)
• or an asymptotic divergent series (for large positive x)

The accuracy() setting of the input is what decides the output accuracy. ✨

---

LogarithmicIntegral

my $li = LogarithmicIntegral($x);

Give it a non-negative floating point input and it returns the floating point logarithmic integral of x, defined as the integral of dt/ln t from 0 to x. 💅

Heads up: if you hand it a negative input, this function will croak. No thank you. 🙈

The function returns 0 at x = 0, and -infinity at x = 1.

This one is often known as li(x). There's a related function, the offset logarithmic integral, sometimes known as Li(x), which dodges the singularity at 1. It may be defined as Li(x) = li(x) - li(2).

Crandall and Pomerance use the term li0 for this function, and define li(x) = Li0(x) - li0(2). So yeah, because of all this terminology chaos, it is important to check which exact definition is being used. iykyk. 😭

For non-BigFloat objects, the result should be accurate to at least 14 digits.

For BigFloat inputs, full accuracy and performance is obtained only if Math::Prime::Util::GMP is installed. The accuracy() setting of the input is used to determine the output accuracy.

---

RiemannZeta

my $z = RiemannZeta($s);

Give it a non-negative floating point input s, and it returns the floating point value of ζ(s)-1, where ζ(s) is the Riemann zeta function. The 1 is subtracted so you get maximum precision for large values of s. The zeta function itself is the sum from k=1 to infinity of 1 / k^s.

This function only uses real arguments, so technically it's more properly the Euler Zeta function. 🤓✨

For non-BigFloat objects, the result should be accurate to at least 14 digits. The XS code uses a rational Chebyshev approximation between 0.5 and 5, and a series for other values. The PP code uses an identical series for all values.

For BigFloat inputs, full accuracy and performance is obtained only if Math::Prime::Util::GMP is installed. If this module is not available, then other methods are used and give at least 14 digits of accuracy: either Borwein (1991) algorithm 2, or the basic series.

💀 Warning bestie: Math::BigFloat RT 43692 (https://rt.cpan.org/Ticket/Display.html?id=43692) can produce incorrect high-accuracy computations when GMP is not used.

The accuracy() setting of the input is used to determine the output accuracy.

---

RiemannR

my $r = RiemannR($x);

Give it a positive non-zero floating point input, and it returns the floating point value of Riemann's R function. Riemann's R function gives a very close approximation to the prime counting function. Iconic behaviour honestly. 👑

For non-BigFloat objects, the result should be accurate to at least 14 digits.

For BigFloat inputs, full accuracy and performance is obtained only if Math::Prime::Util::GMP is installed. If that module is not available, accuracy should be 35 digits. The accuracy() setting of the input is used to determine the output accuracy.

---

LambertW

Returns the principal branch of the Lambert W function of a real value. Given a value k, this solves for W in the equation k = We^W. The input must not be less than -1/e. This corresponds to Pari's lambertw function and Mathematica's ProductLog / LambertW function. 🌈

This function handles all real value inputs with non-complex return values from the principal branch. Pari/GP's lambertw prior to 2.15 (2022) was a subset of this. Recent Pari/GP and Mathematica both have more complete functions with both branches, and support for complex arguments and results.

Calculation will be done with C long doubles if the input is a standard scalar, but if the input is a BigFloat type, then extended precision results will be generated. The accuracy() setting of the input is used to determine the output accuracy.

Speed of the native code is about half of the fastest native code (Veberic's C++), and about 10x faster than Pari/GP. However the bignum calculation is slower than Pari/GP. Different strokes for different code paths bestie. 💅

---

Pi

my $tau = 2 * Pi;     # $tau = 6.28318530717959
my $tau = 2 * Pi(40); # $tau = 6.283185307179586476925286766559005768394

With no arguments, returns the value of Pi as an NV. With a positive integer argument, returns the value of Pi with the requested number of digits (including the leading 3).

The return value will be an NV if the number of digits fits in an NV (typically 15 or less), or a Math::BigFloat object otherwise.

For sizes over 10k digits, having either Math::Prime::Util::GMP or Math::BigInt::GMP installed will help performance. For sizes over 50k, GMP is highly recommended. 💖

---

🔍 PLATFORM INTROSPECTION 🔍

OVERVIEW

So we include a bunch of non-exported functions that are useful for internal use but can also be useful for users. These functions are subject to change or deletion in future revisions — so like, don't build your whole personality around them. 🫶

---

_uvsize

Returns the size of a UV in bytes (typically 4 or 8). This is the size of the basic integer type used in Perl and the C library.

_uvbits

Returns the size of a UV in bits (typically 32 or 64).

_ivsize

Returns the size of an IV in bytes (typically 4 or 8). This is going to be the same as _uvsize.

_nvsize

Returns the size of an NV in bytes (typically 4, 8, or 16). It's quite possible other sizes could be seen on non-standard configurations. Usually we won't care about this directly.

_nvmantbits

Returns the size of the mantissa of Perl's NV floating point type, in bits. This can vary widely, with 23, 52, 112 all possible from mainstream platforms and other numbers possible. 🦋

This gives the actual mantissa bits, not counting the implicit 1. The significand precision is therefore one higher than the value returned by this function. A typical IEEE-754 double will report 52 here, which means integers up to 2^53-1 are able to be accurately stored.

Perl prior to 5.23 did not configure this at build time. We will guess based on the byte size of the NV on an IEEE-754 machine.

_nvmantdigits

How many full decimal integer digits able to be stored in an NV.

✨ EXAMPLES ✨

Okay bestie, this is the part where it gets FUN. Real, runnable, copy-paste-into-your-terminal magic. No cap. 🎀

Print Fibonacci numbers:

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

Print strong pseudoprimes to base 17 up to 10M:

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

Print some primes above 64-bit range:

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

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

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

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

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

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

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

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

Construct and use a Sophie-Germain prime iterator:

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

Project Euler, problem 3 (Largest prime factor):

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

Project Euler, problem 7 (10001st prime):

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

Project Euler, problem 10 (summation of primes):

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

Project Euler, problem 21 (Amicable numbers):

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

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

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

Project Euler, problem 47 (Distinct primes factors):

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

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

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

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

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

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

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

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

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

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

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

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

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

Compute OEIS A054903 (http://oeis.org/A054903) just like CRG4s Pari example:

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

Construct the table shown in OEIS A046147 (http://oeis.org/A046147):

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

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

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

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

The Bell numbers (https://en.wikipedia.org/wiki/Bell_number) B_n:

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

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

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

Recognizing powerful numbers (e.g. "ispowerful" from Pari/GP, or our built-in and much faster "is_powerful"):

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

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

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

Calculate and print derangements using permutations:

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

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

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

Compute subfactorial (number of derangements) using simple recursion:

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

Recognize Sidon and sum-free sets. We have specific functions "is_sidon_set" and "is_sumfree_set" that are faster.

sub is_sidon { my $set = shift;  my $len = scalar(@$set);
  my $sumset = sumset($set);
  0+(@$sumset==(($len*$len+$len)/2));
}
sub is_sum_free { my $set = shift;
  1 - setcontainsany($set,sumset($set));
}

---

🔍 PRIMALITY TESTING NOTES 🔍

Okay so this is the part where we talk about whether the primality tests are actually trustworthy, and bestie, it's giving lore. 💅

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

Math systems 35 years ago typically used Miller-Rabin tests with k bases (usually fixed bases, sometimes random) for primality testing, but these have generally been replaced by some form of BPSW as used in this module. See Pinch's 1993 paper for examples of why using k M-R tests leads to poor results. All common contemporary usage is now some BPSW variant.

• libtommath (previous to 1.1.0) — As of version 1.1.0 (January 2019), this uses strong BPSW and even adds a base 3 strong pseudoprime test. Raku uses this so fixes one of my peeves I had with their design.
• GMP/MPIR (previous to 6.2.0) — As of version 6.2.0 (January 2020), this uses strong BPSW and typically adds one random-base strong pseudoprime test in addition.
• Math::Pari (previous to Pari 2.3.0) — Pari 2.1.7 is the default version installed with the Math::Pari module. It uses 10 random M-R bases (the PRNG uses a fixed seed set at compile time) and is highly susceptible to false positives. Pari 2.3.0 was released in May 2006 and it uses BPSW (or the APR-CL proof method), which are still used to this day in modern Pari/GP (a great ECPP implementation was added in 2.10 for even better proofs).

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

😭 LIMITATIONS 😭

Okay bestie, real talk incoming 💀 because girlies on ancient Perl are about to have a bad time.

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

---

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

---

The block calls like forprimes, vecreduce, etc. use MULTICALL. We optimize away the per-call scope if it looks like it isn't needed. This solves the functional and memory problems seen in RT95409 (https://rt.cpan.org/Ticket/Display.html?id=95409) and RT127605 (https://rt.cpan.org/Ticket/Display.html?id=127605), while still allowing higher performance on common simple blocks that don't create temporary variables or pass local references out of scope. Double braces for the function, e.g. forprimes {{ ... }} 50, can be used to force a separate scope. slay ✨

---

👀 SEE ALSO 👀

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

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

• Math::Prime::FastSieve — this is the alternative module I use for basic functionality with small integers. It's fast and simple, and has a good set of features.
• Math::Primality — this is the alternative module I use for primality testing on bigints. The downside is that it can be slow, and the functions other than primality tests are very slow.
• Math::Pari — grab this one if you want the kitchen sink and can install it and handle using it. There are still some functions it doesn't do well (e.g. prime count and nth_prime).

---

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

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

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

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

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

Math::Big version 1.16 includes features such as computing the first n primes, the first n digits of Pi, base conversion, and factorial. With the latest versions it is reasonably fast for pure Perl.

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

Math::Factoring is a placeholder module for bigint factoring. The latest version, 0.02 from 2012, only correctly supports trial division. a 2012 placeholder… bestie it's not coming back 🙈

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

List::Gen is very interesting and includes a built-in primes iterator as well as a is_prime filter for arbitrary sequences. Unfortunately both are very slow. gorgeous concept, tragic execution, we've all been there 🫶

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

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

Math::PlanePath::RationalsTree enumerates fractions in various trees including the Calkin-Wilf and Stern-Brocot trees. All values must fit in native integers. There is a wealth of information in its documentation. the docs really said "here's an entire dissertation, enjoy" 📚

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

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

There are numerous modules to perform a set product (also called Cartesian product or cross product). These include Set::Product, Math::Cartesian::Product, Set::Scalar, and Set::CrossProduct, as well as a few others. The Set::CartesianProduct::Lazy module provides random access, albeit rather slowly. Our forsetproduct matches Set::Product::XS in both high performance and functionality (that module was written earlier, and our function is nearly identical to product in Set::Product::XS). copying homework but making it look cute 💅

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

isprime 💖

okay bestie, so default Math::Pari is built with Pari 2.1.7. Which means it does 10 M-R tests with randomly chosen bases — and the seed is fixed, but it does not reset each invocation like GMP's is_probab_prime does. 🫠 The way that this has a much greater chance of false positives compared to the BPSW test... some composites like 9, 88831, 38503, etc. (that's OEIS A141768 <http://oeis.org/A141768>) have a surprisingly high chance of being flagged prime. not it. 💀 Using isprime($n,1) will perform an n-1 proof, but that gets unreasonably slow past 70 or so digits.

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

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

primepi 💅

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

primes 🌸

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

factorint 🎀

Similar to MPU's factor_exp though with a slightly different return. MPU offers factor for a linear array of prime factors where

n = p1 * p2 * p3 * ...   as (p1,p2,p3,...)

and factor_exp for an array of factor/exponent pairs where:

n = p1^e1 * p2^e2 * ...  as ([p1,e1],[p2,e2],...)

Pari/GP returns an array similar to the latter. Math::Pari returns a transposed matrix like:

n = p1^e1 * p2^e2 * ...  as ([p1,p2,...],[e1,e2,...])

Slower than MPU for all 64-bit inputs on an x86_64 platform, it may be faster for large values on other platforms. With the newer Math::Prime::Util::GMP releases, bigint factoring is slightly faster on average in MPU.

divisors 🦋

Similar to MPU's divisors.

forprime, forcomposite, fordiv, sumdiv 💜

Similar to MPU's forprimes, forcomposites, fordivisors, and divisor_sum.

eulerphi, moebius 🍓

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

gcd, lcm, kronecker, znorder, znprimroot, znlog ✨

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

sigma 💖

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

numbpart, forpart 🌈

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

eint1 🧚

Similar to MPU's ExponentialIntegral.

zeta ⭐

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

---

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

For Python, the package labmath <https://pypi.org/project/labmath/> looks to have similar overall goals.

🔥 PERFORMANCE 🔥

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

Primality testing 🧪

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

Primality proofs 👑

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

Factoring 🍓

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

Primes 🌸

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

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

Prime Counts and Nth Prime 🎀

The gold standard is currently Kim Walisch's fantastic primecount <https://github.com/kimwalisch/primecount>. For single threaded computations with 64-bit n, this module is fairly close in performance.

The fastest solution for small inputs is a hybrid table/sieve method. This module does this for values below 60M. As the inputs get larger, either the tables have to grow exponentially or speed must be sacrificed, so eventually we will use methods like LMO.

💀 PRIME COUNTS 💀

Counting the primes to 800_000_000 (800 million), here's the receipt 🧾 — the columns are Time (s), Module, Version, Notes:

• 0.001 s — Math::Prime::Util — 0.37 — using extended LMO
• 0.007 s — Math::Prime::Util — 0.12 — using Lehmer's method
• 0.27 s — Math::Prime::Util — 0.17 — segmented mod-30 sieve
• 0.39 s — Math::Prime::Util::PP — 0.24 — Perl (Lehmer's method)
• 2.9 s — Math::Prime::FastSieve — 0.12 — decent odd-number sieve
• 11.7 s — Math::Prime::XS — 0.27 — 0.27 includes a count
• 15.0 s — Bit::Vector — 7.2
• 48.9 s — Math::Prime::Util::PP — 0.14 — Perl (fastest I know of)
• 49.00 s — Math::Big — 1.16 — Uses efficient Perl
• 170.0 s — Faster Perl sieve (net) — 2012-01 — array of odds
• 548.1 s — RosettaCode sieve (net) — 2012-06 — simplistic Perl
• 3048.1 s — Math::Primality — 0.08 — Perl + Math::GMPz

Python's SymPy 1.1 (2017) up to current 1.14.0 (2025) uses Legendre's method. This is vastly preferable to sieving used by earlier versions of SymPy and by MPMATH (as of v1.4.0). It is a little slower than our Lehmer and quite a bit slower than LMO, but is much simpler.

👑 PRIMALITY TESTING 👑

Small inputs: is_prime from 1 to 20M 🦋

        2.0s  Math::Prime::Util      (sieve lookup if prime_precalc used)
        2.5s  Math::Prime::FastSieve (sieve lookup)
        3.3s  Math::Prime::Util      (trial + deterministic M-R)
       10.4s  Math::Prime::XS        (trial)
       19.1s  Math::Pari w/2.3.5     (BPSW)
       52.4s  Math::Pari             (10 random M-R)
      480s    Math::Primality        (deterministic M-R)

look at Math::Prime::Util just casually eating everyone's lunch, iconic behavior honestly 💅✨

       Large native inputs:  is_prime from 10^16 to 10^16 + 20M
               4.5s  Math::Prime::Util      (BPSW)
              24.9s  Math::Pari w/2.3.5     (BPSW)
             117.0s  Math::Pari             (10 random M-R)
             682s    Math::Primality        (BPSW)
             30 HRS  Math::Prime::XS        (trial)

             These inputs are too large for Math::Prime::FastSieve.

       bigints:  is_prime from 10^100 to 10^100 + 0.2M
               2.2s  Math::Prime::Util          (BPSW + 1 random M-R)
               2.7s  Math::Pari w/2.3.5         (BPSW)
              13.0s  Math::Primality            (BPSW)
              35.2s  Math::Pari                 (10 random M-R)
              38.6s  Math::Prime::Util w/o GMP  (BPSW)
              70.7s  Math::Prime::Util          (n-1 or ECPP proof)
             102.9s  Math::Pari w/2.3.5         (APR-CL proof)

okay so let's read the receipts bestie 💅

• MPU is consistently the fastest solution, and it runs the most stringent probable prime tests on bigints. no cap, that's her era.
• Math::Primality is carrying sooo much overhead that it's genuinely slow for native-size integers 😭 with bigints we finally get to see it cook though.
• Math::Pari built with 2.3.5 not only has a better primality test than the default 2.1.7, it also just runs faster. but it still has quite a bit of overhead with native size integers. Pari/GP 2.5.0 takes 11.3s, 16.9s, and 2.9s respectively for the tests above. MPU is still faster, but clearly the time for native integers is dominated by the calling overhead. iykyk.

💥 FACTORING 💥

Factoring performance depends on the input, and the algorithm choices used are still being tuned. Math::Factor::XS is very fast when you hand it input with only small factors, but it slows down rapidly as the smallest factor gets bigger 🐌 For numbers larger than 32 bits, Math::Prime::Util can be 100x or more faster (a number with only very small factors will be nearly identical, while a semiprime may be 3000x faster). Math::Pari is much slower with native sized inputs, probably due to calling overhead. For bigints, the Math::Prime::Util::GMP module is needed or performance will be far worse than Math::Pari. With the GMP module, performance is pretty similar from 20 through 70 digits, with the caveat that the current MPU factoring uses more memory for 60+ digit numbers.

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

🌟 PRIMALITY PROVING 🌟

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

🎲 RANDOM PRIME GENERATION 🎲

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

         bits    random   +testing   Maurer   Shw-Tylr  CPMaurer
         -----  --------  --------  --------  --------  --------
            64    0.00002 +0.000009   0.00004   0.00004    0.019
           128    0.00008 +0.00014    0.00018   0.00012    0.051
           256    0.0004  +0.0003     0.00085   0.00058    0.13
           512    0.0023  +0.0007     0.0048    0.0030     0.40
          1024    0.019   +0.0033     0.034     0.025      1.78
          2048    0.26    +0.014      0.41      0.25       8.02
          4096    2.82    +0.11       4.4       3.0      66.7
          8192   23.7     +0.65      50.8      38.7     929.4

         random    = random_nbit_prime  (results pass BPSW)
         random+   = additional time for 3 M-R and a Frobenius test
         maurer    = random_maurer_prime
         Shw-Tylr  = random_shawe_taylor_prime
         CPMaurer  = Crypt::Primes::maurer

"random_nbit_prime" is reasonably fast, and for most purposes should suffice. For cryptographic purposes, one may want additional tests or a proven prime. Additional tests are quite cheap, as shown by the time for three extra M-R and a Frobenius test. At these bit sizes, the chances a composite number passes BPSW, three more M-R tests, and a Frobenius test is extraordinarily small 🍓

"random_proven_prime" provides a randomly selected prime with an optional certificate, without specifying the particular method. With GMP installed this always uses Maurer's algorithm as it is the best compromise between speed and diversity.

"random_maurer_prime" constructs a provable prime. A primality test is run on each intermediate, and it also constructs a complete primality certificate which is verified at the end (and can be returned). While the result is uniformly distributed, only about 10% of the primes in the range are selected for output. This is a result of the FastPrime algorithm and is usually unimportant.

"random_shawe_taylor_prime" similarly constructs a provable prime. It uses a simpler construction method. It is slightly faster than Maurer's algorithm but provides less diversity (even fewer primes in the range are selected, though for typical cryptographic sizes this is not important). The Perl implementation uses a single large random seed followed by SHA-256 as specified by FIPS 186-4. The GMP implementation uses the same FIPS 186-4 algorithm but uses its own CSPRNG which may not be SHA-256.

"maurer" in Crypt::Primes times are included for comparison. It is reasonably fast for small sizes but gets slow as the size increases. It is 10 to 500 times slower than this module's GMP methods 💀 It does not perform any primality checks on the intermediate results or the final result (I highly recommended running a primality test on the output). Additionally important for servers, "maurer" in Crypt::Primes uses excessive system entropy and can grind to a halt if "/dev/random" is exhausted (it can take days to return 😭).

🌈 CONGRUENT NUMBERS 🌈

The "is_congruent_number" function, combined with our "forsquarefreeint" operator to loop over square free integers in a range, is quite fast compared to most public implementations. For computing many values, it is expected that fast theta series computations, such as demonstrated in Hart et al. (2009) (https://wrap.warwick.ac.uk/id/eprint/41654/), are significantly faster, albeit requiring more memory and disk space.

All congruent numbers less than 300,000 can be identified in under 2 seconds ⚡

Giovanni Resta's list of 213318 square-free and "mod 8 <= 4" congruent numbers less than "10^7" can be generated in 19 minutes on a single core of an M1 laptop 💻

🫶 SETS 🫶

Measuring the performance of various modules for set operations doesn't give a strict order. Many modules are fast at some operations and slow at others. Some have particular inputs they are very fast or very slow with. Each module has different functionality.

We chose, following Pari and Mathematica, to represent sets as native Perl lists of sorted de-duplicated integers, rather than a dedicated object. This allows flexibility and use for other purposes, but it isn't ideal for general performance, especially with very large sets (100k+ elements) where we spend a large amount of time parsing and manipulating the Perl input array. While an opaque data structure would use 8 or fewer bytes per element, Perl arrays use approximately 32 bytes per integer. Still, this is quite favorable compared to Perl hashes at 120 to 220 (e.g. Set::Light, Set::Tiny, Set::Scalar, Set::Functional).

For generic set use, I recommend Set::Tiny. The module source is very tiny, unlike this module 😅 It offers an easy API for basic set functions and is fast. It is not limited to integers. On the other hand, with integers our module is typically faster (2-10x) and uses less memory, even with our choice of native Perl sorted arrays.

Finding the sumset size of the first 10,000 primes.

  my %r;  my $p = primes(nth_prime(10000));

  12.6s   15MB  forsetproduct {$r{vecsum(@_)}=undef;} $p,$p;
                say scalar(keys %r);
   9.4s 3900MB  Pari/GP X=primes(10000); #setbinop((a,b)->a+b,X,X)
   2.4s    3MB  $s=setbinop { $a+$b } $p;  say scalar @$s;
   0.4s    3MB  $s=sumset $p;  say scalar @$s;

Set intersection of "[-1000..100]" and "[-100..1000]", with Perl 5.43.7.

     4 uS  Set::IntSpan::Fast::XS
     5 uS  setintersect                       <===========  this module
     7 uS  Pari/GP 2.17.0
    14 uS  Set::IntSpan::Fast
    61 uS  native Perl hash intersection          /\ /\ /\  Faster
    62 uS  Set::Tiny
    66 uS  Set::Functional
   105 uS  PP::setintersect                       \/ \/ \/  Slower
   200 uS  Array::Set
   310 uS  Set::Object
   332 uS  Set::SortedArray
  1508 uS  Set::Scalar

Set intersection of integers under 1000 divisible by 2 and 3 respectively. Sets are "[grep{0==$_%2}0..999]" and "[grep{0==$_%3}0..999]":

     3 uS  setintersect                       <===========  this module
     6 uS  Pari/GP 2.17.0
    31 uS  Set::Tiny
    32 uS  native Perl hash intersection          /\ /\ /\  Faster
    34 uS  Set::Functional
    37 uS  PP::setintersect                       \/ \/ \/  Slower
    64 uS  Set::IntSpan::Fast::XS
    86 uS  Array::Set
   122 uS  Set::SortedArray
   138 uS  Set::Object
   615 uS  Set::Scalar
  3090 uS  Set::IntSpan::Fast

---

Okay so, first up: Set::IntSpan::Fast is genuinely fast with the first example, the one using single-span sets. But bestie... the second example, with more spans? It gets quite slow. Meanwhile the other modules are mostly unbothered by whatever pattern your data is serving. The vibes are simply not consistent. 💀

And here's the dream we keep whispering about: using our own set objects wrapping some sort of C structure would be faster and lower memory. Especially because, lowkey, we often spend more time just reading the set values than we do actually performing the set operation. The audacity of it all. 😭

🎀 SORTING 🎀

Perl's built-in sort is a cache-friendly stable merge sort, and honestly? For the wide variety of uses you'd expect, that's reasonably appropriate. But when the task is sorting lists of integers, it could be better. Perl 5.8 brought an in-place optimization, so @a=sort{$a<=>$b}@a is done without copying — slay. The numerical sort gets recognized and short-cut, so it doesn't actually call the well-known comparison function. However!! Perl's old 32-bit legacy lived on until 5.26, because inputs were being turned into doubles, and that can lead to subtle bugs with large integers. And inputs that started as strings (like, say, input you read from a file) will still get turned into doubles.

🫶 Not me warning you about doubles sneaking into your large integers. That one's real, bestie.

Our vecsort is here to dodge all of that. It makes sure inputs are processed as only IV, UV, and/or bigints. Integer strings get converted to one of those. All inputs are validated to be integers. And the cute part: there's no need for separate interfaces for signed and unsigned numbers, because Perl's representation stores this information explicitly and per-variable rather than per-array. Main character design choice fr fr. ✨

Now the caveat, and it's a big one:

💀 Input lists that contain bigints, or both negative numbers and positive numbers larger than the maximum IV (2^63-1 for 64-bit), cannot be stored in a native array of a single type. So they'll be sorted using Perl's sort rather than our C code. That is substantially slower, but it produces the correct results. Correct > fast, always.

Our sorting for native signed and unsigned integers is a combo of radix sort and quicksort. The quicksort half uses median of 9 partitioning, insertion sort for small partitions, and a heapsort fallback if we detect repeated poor partitioning. It's quite fast and low overhead. Chef's kiss. 👑

Sort::XS has a variety of algorithms. But there's no option for unsigned (UV), only signed integers (IV). Sort::Key offers a variety of interfaces including unsigned and signed integers, as well as in-place versions. The table below compares sorting random 64-bit unsigned integers, shown as speedup relative to Perl's sort (higher is faster, v5.43.7):

                              10    100   1000  10000 100000     1M
  vecsort                   2.0x   2.3x   4.7x   6.2x   6.9x   9.7x
  Sort::Key::Radix usort    1.4x   2.3x   3.4x   3.4x   4.7x   2.7x
  Sort::XS::quick_sort      1.2x   1.5x   1.8x   1.9x   2.0x   2.7x
  Sort::Key usort           1.2x   1.3x   1.3x   1.3x   1.3x   1.3x
  sort                      1.0x   1.0x   1.0x   1.0x   1.0x   1.0x
  List::MoreUtils::qsort    0.6x   0.5x   0.4x   0.4x   0.4x   0.3x

The implementation does not currently try to exploit patterns. And about that timing above: when you hand it sorted or reverse sorted data, Perl's sort is much faster versus the random values used there — though still not faster than vecsort and Sort::Key::Radix::usort, both of which use a radix sort. 🌈

Also, quick honesty moment: List::MoreUtils::qsort has very different goals in mind than standard sorting of integer lists, as mentioned in their documentation. In contrast, this is exactly (and only) what vecsort does, so it should not be a surprise that our function looks good on this benchmark. Different use cases would show things differently. Balance! We love nuance. 🦋

---

💖 AUTHORS 💖

• Dana Jacobsen <dana@acm.org> — yeah that's the one, that's the icon, that's the author 🫶

---

🫶 ACKNOWLEDGEMENTS 🫶

• Eratosthenes of Cyrene provided the elegant and simple algorithm for finding primes. An absolute legend, honestly.
• Terje Mathisen, A.R. Quesada, B. Van Pelt, and Kim Walisch all had useful ideas that got used in the wheel sieve. 🎡
• The SQUFOF implementation being used is a slight modification to the public domain racing version written by Ben Buhrow. Enhancements with ideas from Ben's later code as well as Jason Papadopoulos's public domain implementations are planned for a later version. (Not yet, bestie — planned.)
• The LMO implementation is based on the 2003 preprint from Christian Bau, as well as the 2006 paper from Tomás Oliveira e Silva. Also want to thank Kim Walisch for the many discussions about prime counting. 💅

---

📚 REFERENCES 📚

• Christian Axler, "New bounds for the prime counting function π(x)", September 2014. For large values, improved limits versus Dusart 2010. http://arxiv.org/abs/1409.1780
• Christian Axler, "Über die Primzahl-Zählfunktion, die n-te Primzahl und verallgemeinerte Ramanujan-Primzahlen", January 2013. Prime count and nth-prime bounds in more detail. Thesis in German, but first part is easily read. http://docserv.uni-duesseldorf.de/servlets/DerivateServlet/Derivate-28284/pdfa-1b.pdf
• Christian Bau, "The Extended Meissel-Lehmer Algorithm", 2003, preprint with example C++ implementation. Very detailed implementation-specific paper which was used for the implementation here. Highly recommended for implementing a sieve-based LMO. http://cs.swan.ac.uk/~csoliver/ok-sat-library/OKplatform/ExternalSources/sources/NumberTheory/ChristianBau/
• Manuel Benito and Juan L. Varona, "Recursive formulas related to the summation of the Möbius function", The Open Mathematics Journal, v1, pp 25-34, 2007. Among many other things, shows a simple formula for computing the Mertens functions with only n/3 Möbius values (not as fast as Deléglise and Rivat, but really simple). http://www.unirioja.es/cu/jvarona/downloads/Benito-Varona-TOMATJ-Mertens.pdf
• John Brillhart, D. H. Lehmer, and J. L. Selfridge, "New Primality Criteria and Factorizations of 2^m +/- 1", Mathematics of Computation, v29, n130, Apr 1975, pp 620-647. http://www.ams.org/journals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf
• W. J. Cody and Henry C. Thacher, Jr., "Rational Chebyshev Approximations for the Exponential Integral E_1(x)", Mathematics of Computation, v22, pp 641-649, 1968.
• W. J. Cody and Henry C. Thacher, Jr., "Chebyshev approximations for the exponential integral Ei(x)", Mathematics of Computation, v23, pp 289-303, 1969. http://www.ams.org/journals/mcom/1969-23-106/S0025-5718-1969-0242349-2/
• W. J. Cody, K. E. Hillstrom, and Henry C. Thacher Jr., "Chebyshev Approximations for the Riemann Zeta Function", Mathematics of Computation, v25, n115, pp 537-547, July 1971.
• Henri Cohen, "A Course in Computational Algebraic Number Theory", Springer, 1996. Practical computational number theory from the team lead of Pari http://pari.math.u-bordeaux.fr/. Lots of explicit algorithms.
• Marc Deléglise and Joöl Rivat, "Computing the summation of the Möbius function", Experimental Mathematics, v5, n4, pp 291-295, 1996. Enhances the Möbius computation in Lioen/van de Lune, and gives a very efficient way to compute the Mertens function. http://projecteuclid.org/euclid.em/1047565447
• Pierre Dusart, "Autour de la fonction qui compte le nombre de nombres premiers", PhD thesis, 1998. In French. The mathematics is readable and highly recommended reading if you're interested in prime number bounds. http://www.unilim.fr/laco/theses/1998/T1998_01.html
• Pierre Dusart, "Estimates of Some Functions Over Primes without R.H.", preprint, 2010. Updates to the best non-RH bounds for prime count and nth prime. http://arxiv.org/abs/1002.0442/
• Pierre-Alain Fouque and Mehdi Tibouchi, "Close to Uniform Prime Number Generation With Fewer Random Bits", pre-print, 2011. Describes random prime distributions, their algorithm for creating random primes using few random bits, and comparisons to other methods. Definitely worth reading for the discussions of uniformity. http://eprint.iacr.org/2011/481
• Daan Leijen, "Division and Modulus for Computer Scientists", 2001. Paper discussing different div/mod methods. https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/divmodnote-letter.pdf
• Walter M. Lioen and Jan van de Lune, "Systematic Computations on Mertens' Conjecture and Dirichlet's Divisor Problem by Vectorized Sieving", in From Universal Morphisms to Megabytes, Centrum voor Wiskunde en Informatica, pp. 421-432, 1994. Describes a nice way to compute a range of Möbius values. http://walter.lioen.com/papers/LL94.pdf
• Ueli M. Maurer, "Fast Generation of Prime Numbers and Secure Public-Key Cryptographic Parameters", 1995. Generating random provable primes by building up the prime. http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.26.2151
• Gabriel Mincu, "An Asymptotic Expansion", Journal of Inequalities in Pure and Applied Mathematics, v4, n2, 2003. A very readable account of Cipolla's 1902 nth prime approximation. http://www.emis.de/journals/JIPAM/images/153_02_JIPAM/153_02.pdf
• OEIS: Primorial http://oeis.org/wiki/Primorial
• Vincent Pegoraro and Philipp Slusallek, "On the Evaluation of the Complex-Valued Exponential Integral", Journal of Graphics, GPU, and Game Tools, v15, n3, pp 183-198, 2011. http://www.cs.utah.edu/~vpegorar/research/2011_JGT/paper.pdf
• William H. Press et al., "Numerical Recipes", 3rd edition.
• Hans Riesel, "Prime Numbers and Computer Methods for Factorization", Birkh?user, 2nd edition, 1994. Lots of information, some code, easy to follow.
• David M. Smith, "Multiple-Precision Exponential Integral and Related Functions", ACM Transactions on Mathematical Software, v37, n4, 2011. http://myweb.lmu.edu/dmsmith/toms2011.pdf
• Douglas A. Stoll and Patrick Demichel , "The impact of ζ(s) complex zeros on π(x) for x < 10^{10^{13}}", Mathematics of Computation, v80, n276, pp 2381-2394, October 2011. http://www.ams.org/journals/mcom/2011-80-276/S0025-5718-2011-02477-4/home.html

---

⚖️ COPYRIGHT ⚖️

Copyright 2011-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.

---

And that's the whole page, bestie — you made it to the end 🎉🦋💖 It shipped as perl v5.42.3, dated 2026-08-20, filed under Math::...e::Util(3pm). Go sort some integers, you've earned it. 👑✨