Fastest code to find the next prime

Question

The problem is as follows.

Input: An integer n

Output: The smallest prime bigger than n.

The challenge is to give the fastest code possible to do this. I will test the code on values starting at size roughly 10^8 10^200 and doubling in size until it takes more than one minute 10 seconds on my computer.

The winning code will find the next prime for the largest input size.

By way of comparison, a simple sieve written in python is able to find the next prime bigger than 10^8 in about 20 seconds.

The requirement that I can test it on my 4GB RAM ubuntu computer is strict. All code must be free (in both senses) and if it uses libraries they must also be free and easily installable. Any false primes reported will immediately disqualify the submission.

I will award separate commendations for the winners in each programming language too if the code is entirely written in that language without the use of external libraries. I will also keep a running table of the fastest times as the competition goes on so people can see how they are doing.

Table so far

• Python. An astonishing 357 digit prime 343239883006530485749095039954069660863471765007165270469723172959277159169882802606127982033072727748864815569574042901856099399985832190628701414555752857600000000000000000000000000000000000000002872284792758930912601189043411951050852357613658978971208596097634095500808832510259693761982135208603287199546795000697807728609476163156438356035166156820611 was the final number under 10 seconds using the code supplied by primo. Will anyone beat this first entry?

Answer 1

Python ~451 digits

This is part of the library I wrote for a semiprime factorization problem, with unnecessary functions removed. It uses the Baillie-PSW primality test, which is technically a probabilistic test, but to date, there are no known pseudoprimes – and there's even a cash reward if you're able to find one (or for supplying a proof that none exist).

Edit: I hadn't realized that Python has built-in modular exponentiation. Replacing my own for the built-in results in a performance boost of about 33%.

my_math.py

# legendre symbol (a|m)
# note: returns m-1 if a is a non-residue, instead of -1
def legendre(a, m):
  return pow(a, (m-1) >> 1, m)

# strong probable prime
def is_sprp(n, b=2):
  d = n-1
  s = 0
  while d&1 == 0:
    s += 1
    d >>= 1

  x = pow(b, d, n)
  if x == 1 or x == n-1:
    return True

  for r in range(1, s):
    x = (x * x)%n
    if x == 1:
      return False
    elif x == n-1:
      return True

  return False

# lucas probable prime
# assumes D = 1 (mod 4), (D|n) = -1
def is_lucas_prp(n, D):
  P = 1
  Q = (1-D) >> 2

  # n+1 = 2**r*s where s is odd
  s = n+1
  r = 0
  while s&1 == 0:
    r += 1
    s >>= 1

  # calculate the bit reversal of (odd) s
  # e.g. 19 (10011) <=> 25 (11001)
  t = 0
  while s > 0:
    if s&1:
      t += 1
      s -= 1
    else:
      t <<= 1
      s >>= 1

  # use the same bit reversal process to calculate the sth Lucas number
  # keep track of q = Q**n as we go
  U = 0
  V = 2
  q = 1
  # mod_inv(2, n)
  inv_2 = (n+1) >> 1
  while t > 0:
    if t&1 == 1:
      # U, V of n+1
      U, V = ((U + V) * inv_2)%n, ((D*U + V) * inv_2)%n
      q = (q * Q)%n
      t -= 1
    else:
      # U, V of n*2
      U, V = (U * V)%n, (V * V - 2 * q)%n
      q = (q * q)%n
      t >>= 1

  # double s until we have the 2**r*sth Lucas number
  while r > 0:
      U, V = (U * V)%n, (V * V - 2 * q)%n
      q = (q * q)%n
      r -= 1

  # primality check
  # if n is prime, n divides the n+1st Lucas number, given the assumptions
  return U == 0

# primes less than 212
small_primes = set([
    2,  3,  5,  7, 11, 13, 17, 19, 23, 29,
   31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
   73, 79, 83, 89, 97,101,103,107,109,113,
  127,131,137,139,149,151,157,163,167,173,
  179,181,191,193,197,199,211])

# pre-calced sieve of eratosthenes for n = 2, 3, 5, 7
indices = [
    1, 11, 13, 17, 19, 23, 29, 31, 37, 41,
   43, 47, 53, 59, 61, 67, 71, 73, 79, 83,
   89, 97,101,103,107,109,113,121,127,131,
  137,139,143,149,151,157,163,167,169,173,
  179,181,187,191,193,197,199,209]

# distances between sieve values
offsets = [
  10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6,
   6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4,
   2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6,
   4, 2, 4, 6, 2, 6, 4, 2, 4, 2,10, 2]

max_int = 2147483647

# an 'almost certain' primality check
def is_prime(n):
  if n < 212:
    return n in small_primes

  for p in small_primes:
    if n%p == 0:
      return False

  # if n is a 32-bit integer, perform full trial division
  if n <= max_int:
    i = 211
    while i*i < n:
      for o in offsets:
        i += o
        if n%i == 0:
          return False
    return True

  # Baillie-PSW
  # this is technically a probabalistic test, but there are no known pseudoprimes
  if not is_sprp(n): return False
  a = 5
  s = 2
  while legendre(a, n) != n-1:
    s = -s
    a = s-a
  return is_lucas_prp(n, a)

# next prime strictly larger than n
def next_prime(n):
  if n < 2:
    return 2
  # first odd larger than n
  n = (n + 1) | 1
  if n < 212:
    while True:
      if n in small_primes:
        return n
      n += 2

  # find our position in the sieve rotation via binary search
  x = int(n%210)
  s = 0
  e = 47
  m = 24
  while m != e:
    if indices[m] < x:
      s = m
      m = (s + e + 1) >> 1
    else:
      e = m
      m = (s + e) >> 1

  i = int(n + (indices[m] - x))
  # adjust offsets
  offs = offsets[m:]+offsets[:m]
  while True:
    for o in offs:
      if is_prime(i):
        return i
      i += o

A sample test script:

from time import clock
from my_math import *

n = i = 317**79
while True:
  i *= 317
  time1 = clock()
  n, o = next_prime(i), n
  span = clock()-time1
  if span > 10:
    break
  print(len(str(n)), span)
print(o)

A factor of 317 was chosen, because it's approximately the square root of 10000, adding roughly 2.5 digits per iteration (and because doubling was too slow to sit through). Output shows the current number of digits, and the time taken.

Sample results:

201 0.13121248650317288
203 0.059535499623555505
206 0.9157767258129175
208 0.2583420518529589
211 0.15367400046653978
213 0.32343915218274955
216 1.3962866788935466
218 0.5986165839513125
221 0.973842206202185
223 2.346910291671148
...
428 0.932809896229827
431 4.345940056627313
433 9.511724255457068
436 6.089835998709333
438 1.3793498894412721
441 4.290633027381972
443 3.5102506044762833
446 3.1629148397352083
448 3.364759208223404
451 7.34668009481652
1551197868099891386459896063244381932060770425565921999885096817830297496627504652115239001983985153119775350914638552307445919773021758654815641382344720913548160379485681746575245251059529720935264144339378936233043585239478807971817857394193701584822359805681429741446927344534491412763713568490429195862973508863067230162660278070962484418979417980291904500349345162151774412157280412235743457342694749679453616265540134456421369622519723266737913

All code is now python 3 compatible.

Answer 2

C++ with GMP: 567 digits

Uses the Miller-Rabin implementation in GMP. It might return a false positive, but good luck actually hitting one with probability 2^-200.

#include <gmp.h>
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>

double time() {
  struct timeval t;
  gettimeofday(&t, NULL);
  return t.tv_usec  * 1e-6 + t.tv_sec;
}

int main(int argc, char *argv[]) {
  mpz_t n, m;
  mpz_init_set_ui(n, 10);
  mpz_pow_ui(n, n, 200);
  mpz_init(m);
  for (int i = 0; true; i++, mpz_mul_ui(n, n, 2)) {
    double start = time();
    for (mpz_add_ui(m, n, 1); !mpz_millerrabin(m, 100); mpz_add_ui(m, m, 2)) ;
    double t = time() - start;
    gmp_printf("%d %Zd %f\n", i, m, t);
    if (t > 10.0) break;
  }
}

Finds the prime 10^200 * 2^1216 + 361 (567 digits) before running over time on my slow laptop.

Answer 3

Perl with GMP module, 1300 digits

Using my module Math::Prime::Util and its GMP back end. The exact crossover point will depend on your computer and whether you have the latest GMP library. All code is free (the modules are on github and CPAN, and GMP is freely available). I've run them on AWS's Ubuntu as well as a desktop Ubuntu (and Fedora, and AIX, and NetBSD, etc.).

The core code is in C and C+GMP. next_prime from MPU sees the number is too large and forwards it to the GMP back end (or pure Perl code if the back end isn't installed). That stringiifies and converts to a mpz, and will convert the result back into the input object type or Math::BigInt. next_prime itself does:

• a mod 30 wheel
• keeps track of the remainder mod 23# so it can do native modulos for primes up to 23
• probable prime test on things that pass these.

The probable prime test is:

• check for tiny divisors using mpz_gcd_ui (in 64-bit two of these check up to 101)
• check for small divisors using singly-calculated large primorials. This is either primes up to 10k or 40k depending on the input size.
• for values larger than 2^1600, performs further trial division using a treesieve. This could be done more efficiently.
• finally, ES BPSW is done (Miller-Rabin test with base 2 followed by an extra strong Lucas test).

Everything prior to the ES BPSW is just optimization, which of course we want for next_prime. next_prime is also implemented in Perl using the Math::BigInt module (in core with optional Pari and GMP back ends). That does AES BPSW (like Pari) but isn't as optimized.

I have thought about the merits of a partial-sieve-based version, using a range of, for instance, 2 merits. I'm just not sure if this would really be better, as most of the time we'd be doing unnecessary sieving as the gap was small, and sometimes for a large gap we'd have to repeat it multiple times.

The library implements ECPP (including certificates) so we could run a proof on the result, but 1200 digits is really too large for the tiny default set of included polynomials (there is a method for downloading larger sets -- proofs would take a bit under 15 minutes which is a bit faster than Pari's APR-CL but a bit slower than WraithX's mpz_aprcl). One downside of ECPP vs. APR-CL is that it has more time variance so there is a good chance of it exceeding 10 seconds on some number before the average time gets there. With a proof I think we're limited to something in the 400 digit range unless we allow multi-threaded software.

#!/usr/bin/env perl
use warnings;
use strict;
use Math::Prime::Util ":all";
use Math::Prime::Util::GMP;  # Barf if the backend isn't installed
use Time::HiRes qw(gettimeofday tv_interval);
use Math::GMP;

my $n = Math::GMP->new(10) ** 200;
while (1) {
  my $start = [gettimeofday];
  my $np = next_prime($n);
  my $sec = tv_interval($start);
  my $len = length($n);
  die "next_prime $len = +",$np-$n," in $sec seconds\n" if $sec > 10;
  warn "  next_prime $len = +",$np-$n," in $sec seconds\n";
  $n *= 10;
}

I decided to try with the same sequence used by primo. It got to 1191 digits, as that is where we hit a gap of 18138. I also tested primo's code using the latest my_math.py. It gets to 630 digits with the 10^e sequence and 641 with his sequence. Very impressive for compact all-Python code without lots of pretests.