Problem jest następujący.

Dane wejściowe: liczba całkowitan

Wyjście: najmniejsza liczba pierwsza większa niż n.

Wyzwanie polega na podaniu najszybszego możliwego kodu. Przetestuję kod na wartościach zaczynających się od rozmiaru z grubsza10^8 10^200 i podwajających rozmiar, aż zajmie to więcej niż minutę i 10 sekund na moim komputerze.

Zwycięski kod znajdzie kolejną liczbę pierwszą dla największego rozmiaru wejściowego.

Dla porównania, proste sito napisane w pythonie jest w stanie znaleźć kolejną liczbę pierwszą większą niż 10^8za około 20sekundy.

Wymóg, że mogę to przetestować na moim 4-GB RAM Ubuntu jest ścisły. Cały kod musi być wolny (w obu aspektach), a jeśli korzysta z bibliotek, muszą być również bezpłatne i łatwe do zainstalowania. Wszelkie zgłoszone fałszywe liczby pierwsze natychmiast zdyskwalifikują zgłoszenie.

Przyznam osobne wyróżnienia dla zwycięzców w każdym języku programowania, jeśli kod jest w całości napisany w tym języku bez użycia zewnętrznych bibliotek. Będę również utrzymywać bieżący stół w najszybszym czasie w trakcie trwania zawodów, aby ludzie mogli zobaczyć, jak sobie radzą.

Tabela do tej pory

• Pyton. Zdumiewającą 357liczbą pierwszą 343239883006530485749095039954069660863471765007165270469723172959277159169882802606127982033072727748864815569574042901856099399985832190628701414555752857600000000000000000000000000000000000000002872284792758930912601189043411951050852357613658978971208596097634095500808832510259693761982135208603287199546795000697807728609476163156438356035166156820611była liczba końcowa poniżej 10 sekund z użyciem kodu dostarczonego przez primo. Czy ktoś pokona ten pierwszy wpis?

Python ~ 451 cyfr

Jest to część biblioteki, którą napisałem dla problemu podziału na czynniki pierwsze , z usuniętymi niepotrzebnymi funkcjami. Wykorzystuje test pierwotności Baillie-PSW , który jest technicznie testem probabilistycznym, ale do tej pory nie są znane pseudopierwsze - a nawet nagroda pieniężna, jeśli możesz ją znaleźć (lub za dostarczenie dowodu, że nie istnieje) .

Edycja : Nie zdawałem sobie sprawy, że Python ma wbudowane modułowe potęgowanie. Zamiana własnego na wbudowane powoduje wzrost wydajności o około 33%.

moja_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

Przykładowy skrypt testowy:

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)

Wybrano współczynnik 317, ponieważ jest to w przybliżeniu pierwiastek kwadratowy z 10000, dodając około 2,5 cyfry na iterację (i ponieważ podwojenie było zbyt wolne, aby usiąść). Dane wyjściowe pokazują bieżącą liczbę cyfr i czas potrzebny.

Przykładowe wyniki:

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

Cały kod jest teraz kompatybilny z Python 3.

C ++ z GMP: 567 cyfr

Wykorzystuje implementację Millera-Rabina w GMP. Może zwrócić fałszywie dodatni, ale powodzenia trafienie w jedno z prawdopodobieństwem 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;
  }
}

Znajduje 10^200 * 2^1216 + 361liczbę pierwszą (567 cyfr), zanim z czasem uruchomię ją na moim wolnym laptopie.

Perl z modułem GMP, 1300 cyfr

Korzystanie z mojego modułu Math :: Prime :: Util i jego zaplecza GMP . Dokładny punkt podziału zależy od komputera i tego, czy masz najnowszą bibliotekę GMP. Cały kod jest bezpłatny (moduły znajdują się na github i CPAN, a GMP jest swobodnie dostępny). Uruchomiłem je na Ubuntu AWS, a także na Ubuntu na pulpicie (i Fedorze, AIX i NetBSD itp.).

Kod podstawowy znajduje się w C i C + GMP. next_prime z MPU widzi, że liczba jest za duża i przesyła ją do zaplecza GMP (lub czystego kodu Perla, jeśli zaplecze nie jest zainstalowane). To przekreślije i konwertuje na mpz, i przekształci wynik z powrotem w typ obiektu wejściowego lub Math :: BigInt. sam next_prime:

• koło mod 30
• śledzi resztę mod 23 #, dzięki czemu może tworzyć natywne moduły dla liczb pierwszych do 23
• prawdopodobny test główny na rzeczy, które je przechodzą.

Prawdopodobny test główny to:

• sprawdź małe dzielniki za pomocą mpz_gcd_ui (w 64-bitowych dwóch z nich sprawdź do 101)
• sprawdź małe dzielniki, używając pojedynczo obliczonych dużych pierwiastków. Jest to liczba pierwsza do 10 000 lub 40 000 w zależności od rozmiaru wejściowego.
• dla wartości większych niż 2 ^ 1600 wykonuje dalszy podział próbny za pomocą sita drzewnego. Można to zrobić bardziej efektywnie.
• wreszcie ES BPSW jest zakończony (test Millera-Rabina z bazą 2, a następnie wyjątkowo silny test Lucasa ).

Wszystko przed ES BPSW to tylko optymalizacja, której oczywiście chcemy na next_prime. next_prime jest również zaimplementowany w Perlu przy użyciu modułu Math :: BigInt (w rdzeniu z opcjonalnymi zapleczami Pari i GMP). To robi AES BPSW (jak Pari), ale nie jest tak zoptymalizowany.

Myślałem o zaletach wersji opartej na sicie częściowym, wykorzystując zakres, na przykład, 2 zalet. Po prostu nie jestem pewien, czy to naprawdę byłoby lepsze, ponieważ przez większość czasu robiliśmy niepotrzebne przesiewanie, ponieważ przerwa była niewielka, a czasami dla dużej przerwy musielibyśmy to powtarzać wiele razy.

Biblioteka implementuje ECPP (w tym certyfikaty), abyśmy mogli przetestować wynik, ale 1200 cyfr jest naprawdę zbyt duże dla niewielkiego domyślnego zestawu dołączonych wielomianów (istnieje metoda pobierania większych zestawów - proofy zajęłyby nieco mniej 15 minut, co jest nieco szybsze niż APR-CL Pari, ale nieco wolniejsze niż mpz_aprcl WraithX). Jedną wadą ECPP w porównaniu z APR-CL jest to, że ma większą zmienność czasową, więc istnieje duża szansa, że ​​przekroczy on 10 sekund w stosunku do pewnej liczby, zanim upłynie średni czas. Z dowodem uważam, że jesteśmy ograniczeni do czegoś w zakresie 400 cyfr, chyba że zezwolimy na oprogramowanie wielowątkowe.

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

Postanowiłem spróbować z tą samą sekwencją, której używał primo. Doszedł do 1191 cyfr, ponieważ w tym miejscu osiągnęliśmy lukę 18138. Testowałem również kod primo przy użyciu najnowszego my_math.py. Dochodzi do 630 cyfr z sekwencją 10 ^ e i 641 z jego sekwencją. Bardzo imponujące, jak na kompaktowy, całkowicie Pythonowy kod bez wielu testów wstępnych.