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 $_), $_->[0], split /,\s*/, $_->[1]]
   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[0] 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!

Updates

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,

http://probableprime.org/download/example-primality.c

shows a 2-base hash solution for 2^49 and 3-base hash for 64-bit. I wrote that almost two years ago, and maybe I should publish it. I've mentioned it on a few forums, but honestly I think BPSW is a better solution for 64-bit.