Schnellster Code, um die nächste Primzahl zu finden

17

Das Problem ist wie folgt.

Eingabe: Eine ganze Zahln

Output: Die kleinste Primzahl größer als n.

Die Herausforderung besteht darin, den schnellstmöglichen Code dafür anzugeben. Ich werde den Code auf Werten testen, die ungefähr10^8 bei der Größe beginnen 10^200und sich verdoppeln, bis es auf meinem Computer mehr als eine Minute und 10 Sekunden dauert .

Der Gewinncode findet die nächste Primzahl für die größte Eingabegröße.

Im Vergleich dazu , ein einfaches Sieb in Python geschrieben ist in der Lage die nächste Primzahl größer zu finden , als 10^8in etwa 20Sekunden.

Die Anforderung, dass ich es auf meinem 4 GB RAM Ubuntu-Computer testen kann, ist streng. Der gesamte Code muss (in beiden Richtungen) frei sein, und wenn Bibliotheken verwendet werden, müssen sie auch frei und leicht installierbar sein. Alle gemeldeten falschen Primzahlen disqualifizieren die Einsendung sofort.

Ich werde auch für die Gewinner in jeder Programmiersprache separate Auszeichnungen vergeben, wenn der Code vollständig in dieser Sprache geschrieben ist und keine externen Bibliotheken verwendet werden. Ich werde auch einen Lauftisch mit den schnellsten Zeiten halten, während der Wettbewerb weitergeht, damit die Leute sehen können, wie es ihnen geht.

Tisch so weit

  • Python. Eine erstaunliche 357Ziffer Primzahl 343239883006530485749095039954069660863471765007165270469723172959277159169882802606127982033072727748864815569574042901856099399985832190628701414555752857600000000000000000000000000000000000000002872284792758930912601189043411951050852357613658978971208596097634095500808832510259693761982135208603287199546795000697807728609476163156438356035166156820611war die letzte Zahl unter 10 Sekunden mit dem Code von primo. Wird jemand diesen ersten Eintrag schlagen?
felipa
quelle
1
Fast ein genaues Duplikat von Code-Challenge: The Nearest Prime
Peter Taylor
@PeterTaylor Bei dieser Frage geht es meiner Meinung nach um Zeitkomplexität. Hier geht es um die praktische Geschwindigkeit in Sekunden. Ich denke, diese beiden Dinge können sehr unterschiedlich sein.
Felipa
Klar, wenn man sich an kleine Testfälle hält. Da sich jedoch niemand die Mühe gemacht hat, AKS für die andere Frage zu implementieren, erhalten Sie die gleichen Antworten.
Peter Taylor
3
@PeterTaylor erlaube mir nicht zuzustimmen. Schließlich sollten 90% des Datenverkehrs einer Website von Suchmaschinen stammen . Eine Google-Suche nach schneller Semiprime-Faktorisierung und Multiple Polynomial Quadratic Sieve gibt das ursprüngliche Problem zurück, bei dem ich meinen Code an Platz 2 bzw. 4 entnommen habe. Ich stelle mir irgendwann vor, dass dieses Problem auch ziemlich hoch sein wird fast next prime function.
Primo
1
Ich denke, das OP hat seine Tests der Antworten nicht aktualisiert ...
mbomb007

Antworten:

21

Python ~ 451 Stellen

Dies ist ein Teil der Bibliothek, die ich für ein Semiprime-Faktorisierungsproblem geschrieben habe , bei dem unnötige Funktionen entfernt wurden. Es wird der Baillie-PSW-Primalitätstest verwendet , der technisch gesehen ein probabilistischer Test ist. Bisher sind jedoch keine Pseudoprimes bekannt - und es gibt sogar eine Geldprämie, wenn Sie eine finden können (oder wenn Sie einen Beweis vorlegen, dass keine existieren). .

Edit : Ich hatte nicht bemerkt, dass Python modulare Potenzierung eingebaut hat. Das Ersetzen meiner eigenen durch die eingebauten führt zu einer Leistungssteigerung von ca. 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

Ein Beispiel für ein Testskript:

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)

Ein Faktor von 317 wurde gewählt, weil er ungefähr die Quadratwurzel von 10000ist und ungefähr 2,5 Stellen pro Iteration hinzufügt (und weil das Verdoppeln zu langsam war, um durchzuhalten). Die Ausgabe zeigt die aktuelle Anzahl der Stellen und die benötigte Zeit.

Beispielergebnisse:

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

Der gesamte Code ist jetzt Python 3-kompatibel.

primo
quelle
Das geht erstaunlich schnell! Ich werde es in ein paar Tagen korrekt ausführen und die Größe verdoppeln (und einen deterministischen Primalitätstest durchführen) und die größte Zahl in die Tabelle aufnehmen. Ich vermute jedoch, dass Sie bereits der Gewinner sind.
Felipa
1
FWIW, in Sage, next_prime((2^520)*(10^200))ungefähr 15 Sekunden auf meiner Maschine, also auf den ersten Blick ist das ziemlich beeindruckend. Es next_prime((2^520)*(10^200),proof=False)dauert jedoch ... 0,4 Sekunden, da nur auf Pseudoprimalität geprüft wird. Ihr Anspruch : „Es gibt keine bekannte Pseudoprimzahlen“ zu überzeugen , ist verschwindend als die Anzahl der Bits geht über 64. 357 Ziffern, ich bin nicht einmal entfernt durch einen Mangel an Gegenbeispielen überzeugt.
Stand
@boothby es ist erwähnenswert, dass dies die gleiche Methode ist, die von Maple verwendet wird . Dass die Methode vor 33 Jahren veröffentlicht wurde und noch keine Pseudoprimes bekannt sind, spricht für ihre Genauigkeit.
Primo
1
Deshalb benutze ich Salbei. "Nicht bekannt, dass es versagt" ist eigentlich nicht dasselbe wie "bekannt, dass es funktioniert". Angenommen, es gab eine falsche Pseudoprime mit weniger als 400 Stellen. Es würde Billionen von Jahren dauern, um es zu finden - aber es wäre immer noch da und würde jeden Versuch vereiteln, 'Pseudoprime = Prime' zu beweisen. Ich werde immer "Lösungen" ablehnen, die probabalistische Methoden ohne Garantie anwenden. Monte Carlo? Sichere Sache. "Es ist primär, weil ein Zauberer mir sagte, dass es wahrscheinlich war"? Nee.
Stand
1
@boothby Sie müssen eine Antwort hinzufügen, damit wir sie kommentieren können :)
felipa
6

C ++ mit GMP: 567 Stellen

Verwendet die Miller-Rabin-Implementierung in GMP. Es könnte ein falsch positives Ergebnis liefern, aber viel Glück beim Treffen eines mit einer Wahrscheinlichkeit von 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;
  }
}

Findet die Primzahl 10^200 * 2^1216 + 361(567 Stellen), bevor sie mit der Zeit auf meinem langsamen Laptop ausgeführt wird.

Keith Randall
quelle
3

Perl mit GMP-Modul, 1300 Stellen

Mit meinem Modul Math :: Prime :: Util und seinem GMP- Backend . Der genaue Übergangspunkt hängt von Ihrem Computer ab und davon, ob Sie über die neueste GMP-Bibliothek verfügen. Der gesamte Code ist kostenlos (die Module sind auf Github und CPAN und GMP ist frei verfügbar). Ich habe sie auf AWSs Ubuntu sowie auf einem Desktop-Ubuntu (und Fedora, AIX und NetBSD usw.) ausgeführt.

Der Kerncode ist in C und C + GMP. next_prime von MPU erkennt, dass die Nummer zu groß ist und leitet sie an das GMP-Back-End weiter (oder reinen Perl-Code, wenn das Back-End nicht installiert ist). Dadurch wird das Ergebnis in einen mpz-Wert umgewandelt und wieder in den Eingabeobjekttyp oder Math :: BigInt konvertiert. next_prime selbst macht:

  • ein mod 30 rad
  • Verfolgt den Rest von Mod 23 #, damit native Module für Primzahlen bis 23 erstellt werden können
  • wahrscheinlicher Primetest auf Dingen, die diese bestehen.

Der wahrscheinliche Primetest ist:

  • Mit mpz_gcd_ui nach winzigen Teilern suchen (in 64-Bit-Zwei von diesen bis zu 101 prüfen)
  • Mit einfach berechneten großen Primorials nach kleinen Teilern suchen. Dies sind entweder Primzahlen bis zu 10k oder 40k, abhängig von der Eingabegröße.
  • Für Werte größer als 2 ^ 1600 wird eine weitere Versuchsteilung mit einem Baumsieb durchgeführt. Dies könnte effizienter erfolgen.
  • Schließlich wird ES BPSW durchgeführt (Miller-Rabin-Test mit Basis 2, gefolgt von einem besonders starken Lucas-Test ).

Alles vor dem ES BPSW ist nur Optimierung, was wir natürlich für next_prime wollen. next_prime wird auch in Perl mithilfe des Math :: BigInt-Moduls implementiert (im Kern mit optionalen Pari- und GMP-Backends). Das macht AES BPSW (wie Pari), ist aber nicht so optimiert.

Ich habe über die Vorzüge einer auf einem Teilsieb basierenden Version nachgedacht und dabei eine Reihe von beispielsweise zwei Vorzügen verwendet. Ich bin mir nur nicht sicher, ob dies wirklich besser wäre, da wir die meiste Zeit unnötig gesiebt haben, da die Lücke klein war, und manchmal mussten wir sie für eine große Lücke mehrmals wiederholen.

Die Bibliothek implementiert ECPP (einschließlich Zertifikate), so dass wir einen Proof für das Ergebnis ausführen können, aber 1200 Stellen sind wirklich zu groß für den winzigen Standardsatz der enthaltenen Polynome (es gibt eine Methode zum Herunterladen größerer Sätze - Proofs würden etwas länger dauern) 15 Minuten, etwas schneller als der APR-CL von Paris, aber etwas langsamer als der mpz_aprcl von WraithX). Ein Nachteil von ECPP im Vergleich zu APR-CL ist, dass es mehr Zeitabweichungen gibt, so dass die Wahrscheinlichkeit hoch ist, dass es bei einer bestimmten Zahl 10 Sekunden überschreitet, bevor die durchschnittliche Zeit dort ankommt. Mit einem Beweis, denke ich, sind wir auf etwas im 400-stelligen Bereich beschränkt, sofern wir keine Multithread-Software zulassen.

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

Ich habe mich entschlossen, es mit der gleichen Sequenz zu versuchen, die auch primo verwendet. Es wurden 1191 Stellen erreicht, da wir hier eine Lücke von 18138 erreichten. Ich habe auch Primos Code mit der neuesten Datei my_math.py getestet. Es werden 630 Stellen mit der Sequenz 10 ^ e und 641 mit seiner Sequenz erreicht. Sehr beeindruckend für kompakten Python-Code ohne viele Vortests.

DanaJ
quelle
Ich komme immer noch nicht darüber hinweg, wie schnell dieses Modul ist. Es hat mein Interesse an Perl als Werkzeug zum Knacken von Zahlen wiederbelebt. Ich schreibe gerade Math::GMPauf eine Weise um, die mit der Erstellung / Zerstörung von mpz-Referenzen nicht ganz so verschwenderisch ist.
Primo
Die eigentliche Arbeit ist alles in C + GMP, also könnte alles auch in Python funktionieren. Python hat einige gravierende Vorteile gegenüber Perl 5 für große Zahlen, von denen ich mir eine Lösung wünschte. Math :: GMPz ist übrigens schneller als Math :: GMP und hat im Grunde die gesamte mpz-API verfügbar gemacht, obwohl sie spröder und manchmal etwas seltsam zu nennen ist. Einige Dinge in Math :: GMP zu korrigieren, steht auf meiner ToDo-Liste hinter zu vielen anderen Dingen. Was MPU betrifft, habe ich darüber nachgedacht, die Entwicklung zu invertieren und in zwei C-Bibliotheken umzuwandeln. Dann muss das Perl-Modul das nur verwenden. Es würde helfen, es woanders einzusetzen.
DanaJ
Ich mache gute Fortschritte. Die folgende Schleife läuft über 10 - mal so schnell , allein durch besseres Referenz - Management: $x = new Math::GMP(0); $x += 3 for 1..1000000. Ich werde auf cpan posten, wenn ich fertig bin. Du wirst einer der ersten sein, die es wissen;)
Primo