Le problème est le suivant.

Entrée: un entiern

Sortie: Le plus petit nombre premier plus grand que n.

Le défi est de donner le code le plus rapide possible pour ce faire. Je vais tester le code sur des valeurs commençant à peu près10^8 à la taille 10^200et doublant de taille jusqu'à ce que cela prenne plus d' une minute 10 secondes sur mon ordinateur.

Le code gagnant trouvera le premier nombre premier pour la plus grande taille d'entrée.

À titre de comparaison, un simple tamis écrit en python est capable de trouver le premier nombre supérieur plus gros qu'en quelques secondes 10^8environ 20.

L'obligation de pouvoir le tester sur mon ordinateur ubuntu de 4 Go de RAM est stricte. Tout le code doit être gratuit (dans les deux sens) et s'il utilise des bibliothèques, il doit également être gratuit et facilement installable. Tout faux nombre premier signalé disqualifiera immédiatement la soumission.

J'attribuerai également des mentions élogieuses aux gagnants dans chaque langage de programmation si le code est entièrement écrit dans ce langage sans utiliser de bibliothèques externes. Je garderai également un tableau des meilleurs temps au fur et à mesure de la compétition afin que les gens puissent voir comment ils vont.

Table jusqu'ici

• Python. Un 357nombre premier étonnant 343239883006530485749095039954069660863471765007165270469723172959277159169882802606127982033072727748864815569574042901856099399985832190628701414555752857600000000000000000000000000000000000000002872284792758930912601189043411951050852357613658978971208596097634095500808832510259693761982135208603287199546795000697807728609476163156438356035166156820611était le nombre final inférieur à 10 secondes utilisant le code fourni par primo. Est-ce que quelqu'un battra cette première entrée?

Python ~ 451 chiffres

Cela fait partie de la bibliothèque que j'ai écrite pour un problème de factorisation semi - primaire , avec des fonctions inutiles supprimées. Il utilise le test de primalité Baillie-PSW , qui est techniquement un test probabiliste, mais à ce jour, il n'y a pas de pseudoprimes connus - et il y a même une récompense en espèces si vous en trouvez un (ou pour fournir une preuve qu'il n'en existe pas) .

Edit : je n'avais pas réalisé que Python avait une exponentiation modulaire intégrée. Remplacer le mien pour les résultats intégrés entraîne une augmentation des performances d'environ 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

Un exemple de script de test:

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)

Un facteur de 317 a été choisi, car il s'agit approximativement de la racine carrée de 10000, ajoutant environ 2,5 chiffres par itération (et parce que le doublement était trop lent pour s'asseoir). La sortie affiche le nombre actuel de chiffres et le temps nécessaire.

Exemples de résultats:

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

Tout le code est désormais compatible python 3.

C ++ avec GMP: 567 chiffres

Utilise l'implémentation Miller-Rabin dans GMP. Cela pourrait renvoyer un faux positif, mais bonne chance en frapper un avec une probabilité de 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;
  }
}

Trouve le premier 10^200 * 2^1216 + 361(567 chiffres) avant de courir dans le temps sur mon ordinateur portable lent.

Perl avec module GMP, 1300 chiffres

Utiliser mon module Math :: Prime :: Util et son back-end GMP . Le point de croisement exact dépendra de votre ordinateur et de votre dernière bibliothèque GMP. Tout le code est gratuit (les modules sont sur github et CPAN, et GMP est disponible gratuitement). Je les ai exécutés sur Ubuntu d'AWS ainsi que sur Ubuntu de bureau (et Fedora, et AIX, et NetBSD, etc.).

Le code principal est en C et C + GMP. next_prime de MPU voit que le nombre est trop grand et le transmet au backend GMP (ou du code Perl pur si le backend n'est pas installé). Cela chaîne et convertit en mpz, et convertira le résultat en arrière dans le type d'objet d'entrée ou Math :: BigInt. next_prime lui-même fait:

• une roue mod 30
• garde la trace du mod 23 # restant afin qu'il puisse faire des modules natifs pour les nombres premiers jusqu'à 23
• test primaire probable sur les choses qui réussissent.

Le test premier probable est:

• vérifier les minuscules diviseurs à l'aide de mpz_gcd_ui (dans deux de 64 bits, vérifier jusqu'à 101)
• vérifier les petits diviseurs à l'aide de grandes primitives calculées individuellement. Il s'agit soit d'amorces jusqu'à 10k ou 40k selon la taille d'entrée.
• pour des valeurs supérieures à 2 ^ 1600, effectue une nouvelle division d'essai à l'aide d'un tamis. Cela pourrait être fait plus efficacement.
• enfin l'ES BPSW est fait (test de Miller-Rabin avec base 2 suivi d'un test de Lucas extra fort ).

Tout ce qui précède l'ES BPSW n'est qu'une optimisation, ce que nous voulons bien sûr pour next_prime. next_prime est également implémenté en Perl en utilisant le module Math :: BigInt (en core avec des backends optionnels Pari et GMP). Cela fait AES BPSW (comme Pari) mais n'est pas aussi optimisé.

J'ai réfléchi aux avantages d'une version à tamis partiel, en utilisant par exemple une gamme de 2 avantages. Je ne suis tout simplement pas sûr que ce serait vraiment mieux, car la plupart du temps, nous effectuions un tamisage inutile car l'écart était petit, et parfois pour un grand écart, nous devions le répéter plusieurs fois.

La bibliothèque implémente ECPP (y compris les certificats) afin que nous puissions exécuter une épreuve sur le résultat, mais 1200 chiffres est vraiment trop grand pour le petit ensemble par défaut de polynômes inclus (il existe une méthode pour télécharger des ensembles plus grands - les épreuves prendraient un peu moins de 15 minutes, ce qui est un peu plus rapide que l'APR-CL de Pari mais un peu plus lent que le mpz_aprcl de WraithX). Un inconvénient d'ECPP par rapport à APR-CL est qu'il a plus de variance de temps, donc il y a de fortes chances qu'il dépasse 10 secondes sur un certain nombre avant que le temps moyen n'y arrive. Avec une preuve, je pense que nous sommes limités à quelque chose dans la plage de 400 chiffres, sauf si nous autorisons les logiciels multithreads.

#!/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;
}

J'ai décidé d'essayer avec la même séquence utilisée par primo. Il a atteint 1191 chiffres, car c'est là que nous avons atteint un écart de 18138. J'ai également testé le code de primo en utilisant le dernier my_math.py. Il atteint 630 chiffres avec la séquence 10 ^ e et 641 avec sa séquence. Très impressionnant pour du code compact tout-Python sans beaucoup de pré-tests.