THE RANT / THE SCHPLOG
Schmorp's POD Blog a.k.a. THE RANT
a.k.a. the blog that cannot decide on a name

This document was first published 2016-03-07 09:21:06, and last modified 2016-03-13 03:53:59.

# Utility time: find the next prime larger than

Today it's utility time: I'm publishing here a little utility that sits on my disk for many many years now, but hasn't been published for the lack of a more suitable medium.

The sole purpose of it is to find the next primes equal to or larger than a given number. Not primes as in RSA primes (i.e. with a lot of digits), but primes as in hash table sizes or other applications, such as scheduling recurring tasks with a minimum of overlap.

It's not something that I need often, but occasionally I just need one or more prime larger than, say, 10000 or so, and searching in the Debian GNU/Linux packages didn't find me an obvious candidate. So I wrote one myself: primenext (well, for that reason, and also because I wanted to learn how `Math::BigInt` works. I don't like it very much).

## Usage

The usage is very simple:

```\$ primenext 65536 | head
65537
65539
65543
65551
65557
65563
65579
65581
65587
65599```

It works with numbers up to `2**64`, and although the main goal was for the program to be as simple and short as possible, it is still "fast enough" for finding magic program constants:

```\$ time primenext 10_000000_000000_000000 | head
10000000000000000051
10000000000000000087
10000000000000000091
10000000000000000097
10000000000000000099
10000000000000000147
10000000000000000169
10000000000000000273
10000000000000000297
10000000000000000307

real    0m0.054s
user    0m0.052s
sys     0m0.000s```

## Theory, errrrr, Practise!

I won't go into the maths behind this program - every child knows the Miller-Rabin primality test which often is used to generate really large random primes for cryptography applications. It's probabilistic, but if you can make the chances of a bad result almost infinitely less likely than getting hit by lighting in the next minute, it's usually fine. Still, I wanted something better, so the program is using a deterministic Miller-Rabin primality test.

For this, you need a list of strong pseudoprimes and the limit to which they can reliably detect primes. Only one set is needed, but for speed reasons, a number of these sets is implemented, with the script choosing the appropriate one:

```use Math::BigInt try => "GMP";

my @sprp = (
[              "341531", "9345883071009581737"],
[           "520924141", "15, 750068417525532"],
[        "109134866497", "2, 45650740, 3722628058"],
[      "47636622961201", "2, 2570940, 211991001, 3749873356"],
[    "3770579582154547", "2, 2570940, 880937, 610386380, 4130785767"],
["18446744073709551615", "2, 325, 9375, 28178, 450775, 9780504, 1795265022"],
#   ^^^^^^^^^^^^^^^^^^^^^^  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
#          limit .............. pseudo primes
);

\$_ = [map +(new Math::BigInt \$_), \$_->, split /,\s*/, \$_->]
for @sprp;```

Wikipedia has a list, but this page has a better one. As you can see, there is still room for improvement, although I would suggest not using `Math::BigInt` for record searches.

The main program simply starts at \$ARGV and iterates, switching to ever larger SPRP sets as required:

```my \$n = (new Math::BigInt shift) | 1;

my @bases;
my \$max = 0;

while () {
while (\$n > \$max) {
@sprp
or die "range exceeded\n";
(\$max, @bases) = @{ shift @sprp };
}
print "\$n\n" if isprime \$n;
\$n += 2;
}```

The actual maths is carried out by `isprime`:

```# requires @bases to be set properly before call
sub isprime(\$) {
my (\$n) = @_;

\$n % \$_ || \$n == \$_ || return
for 3, 5, 7; # 2 not tested

return 1 if \$n < 121;

my (\$d, \$s) = \$n - 1;
\$s++ until \$d & (1<<\$s);
\$d >>= \$s;

a:
for my \$a (@bases) {
my \$x = \$a->copy->bmodpow (\$d, \$n);
next if \$x == 1 || \$x == \$n - 1;
for (1 .. \$s - 1) {
\$x = \$x * \$x % \$n;
return if \$x == 1; # warn "\$n == 0 mod ", Math::BigInt::bgcd (\$a ** (\$d * 2 ** \$_ / 2) - 1, \$n);
next a if \$x == \$n - 1;
}
return;
}

1
}```

Anyways, not much else to be said, you can download the `primenext` source code and be happy with it or not.

Ah, one last thing, the illustrious Steve Worley who is responsible for many of the record entries also has a nifty way to do really fast primality tests, for those times where you need it, using a smallish lookup table and a single test.

## Hints welcome

Are there better SPRP lists? Probably not, but if you find one (or an altogether better method), tell me. Likewise, if you find a bug or some material improvement in this script, please tell me as well!

As b_jonas pointed out on IRC, the `bsdgames` package contains a primes program that is a slight superset of the primenext program. Also, it contains rot13, which is "why it's installed everywhere".

Also, Dana Jacobsen sent a mail with so much additional interesting info that I'll just quote it here in full (any editorial mistakes are mine):

ntheory (aka `Math::Prime::Util`) comes with a primes.pl program. It isn't quite the same, but you can do:

`primes.pl 65536 +100`

to get primes between 65536 and 65536+100. It also takes expressions for things like `primes.pl 2**9 "nth_prime(105)"`

You could also use ntheory for the next_prime loop (or `forprimes` or an iterator), but then you don't have as much fun. You do get support for bigints, results use ES BPSW, and you could check each with `is_provable_prime` if you really wanted, which shouldn't slow things down until 200 or so digits.

For 64-bit, Wojciech's page is the best reference. There are some results for larger than 64-bit now, which is on Wikipedia.

Regarding hashing, there are some better results than Worley's. For 32-bit ntheory uses a single-base hashed solution similar to Forisek and Jancina. For larger that 32-bit,