Take the 2-minute tour ×
Programming Puzzles & Code Golf Stack Exchange is a question and answer site for programming puzzle enthusiasts and code golfers. It's 100% free, no registration required.

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?
share|improve this question
1  
Almost an exact duplicate of Code-Challenge:The Nearest Prime –  Peter Taylor Feb 16 '13 at 11:00
    
@PeterTaylor That question is about time complexity I think. This is about practical speed in seconds. I think those two things can be quite different. –  felipa Feb 16 '13 at 11:03
    
Sure, if you stick to small test cases. But since no-one bothered to implement AKS for the other question, you're going to get the same answers. –  Peter Taylor Feb 16 '13 at 11:12
    
@PeterTaylor There is already an amazing new answer here. I will have to test with huge test cases now! –  felipa Feb 16 '13 at 11:25
2  
@PeterTaylor allow me to disagree. Eventually, 90% of a site's traffic should come from search engines. A google search for fast semiprime factorization and Multiple Polynomial Quadratic Sieve return the original problem I've taken my code from at place #2 and #4 respectively. I imagine at some point, this problem will rank fairly high for fast next prime function as well. –  primo Feb 18 '13 at 15:04
show 4 more comments

2 Answers

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.

share|improve this answer
    
That is astonishingly fast! I will run it properly with doubling size in a few days (and a deterministic primality test) and put the biggest number in the table. I suspect you may already be the winner though. –  felipa Feb 16 '13 at 11:08
    
FWIW, in Sage, next_prime((2^520)*(10^200)) about 15 seconds on my machine, so at first blush this is quite impressive. However... next_prime((2^520)*(10^200),proof=False) takes 0.4 seconds because it only checks for pseudoprimality. Your claim "there are no known pseudoprimes" is vanishingly convincing as the number of bits goes over 64. For 357 digits, I'm not even remotely convinced by a lack of counterexamples. –  boothby Feb 19 '13 at 3:36
    
@boothby it's worth noting that this is the very same method used by Maple. That the method was published 33 years ago, and there still no known pseudoprimes speaks for its degree of accuracy. –  primo Feb 19 '13 at 13:09
    
This is why I use Sage. "Not known to fail" is really actually not the same as "known to work". Suppose there was one false pseudoprime under 400 digits. It would take trillions of years to find it -- but it'd still be there, foiling any attempt to prove 'pseudoprime = prime'. I will always downvote "solutions" that use probabalistic methods with zero guarantee. Monte Carlo? Sure thing. "It's prime 'cause a wizard told me it probably was"? Nope. –  boothby Feb 19 '13 at 17:10
    
@boothby Is sage running AKS to check? –  felipa Feb 19 '13 at 20:50
show 5 more comments

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.