Primzahlen und Faktorisierung

Posted by André Mühlnikel at March 14, 2009 12:26 p.m.

Dank projecteuler.net habe ich mal wieder angefangen, mich mit Perl auseinander zu setzen. Der größte Brocken bisher: Primzahl-Faktorisierung einer ziemlich großen Zahl. Nach ein paar ersten ziemlich erfolgversprechenden Versuchen, ein binäres Sieb_des_Eratosthenes zu bauen (also jede Zahl von 1 bis N wird durch ein Bit dargestellt - 0 = keine Primzahl, 1 = Primzahl), musste ich feststellen, dass das sehr schnell meinen Arbeitsspeicher sprengt (600 Mrd Bit = 75 GB), zumal das allermeiste davon nur 0en sind. Ziemlich unbefriedigend, also auf zur „Brute-Force“-Lösung des Siebes: man nehme alle bekannten Primzahlen (mindestens also „2“ als kleinste Primzahl) und probiere für jede neue Zahl aus, ob sie durch eine der Primzahlen teilbar ist. Das dauert aber UNENDLICH lange (jedenfalls fast). Also optimieren:

Erstmal ist der größte Primfaktor von N, den es geben kann genau N/2, wenn N keine Primzahl ist. Habe ich aber nach sqrt(N) noch keinen primen Teiler gefunden, ist N zwingend selbst eine Primzahl:

Sei N nicht Primzahl und sei M eine Primzahl, M > sqrt(N). Dann ist P = N / M immer kleiner als sqrt(N)! Entweder ist P selbst eine Primzahl oder in solche zerlegbar (die dann wiederum < P sind), jedoch hätte dann beim Suchen von nach Primfaktoren im Bereich 1 bis M bereits P oder seine Primfaktoren auftauchen müssen. Sind sie das, so muss ohnehin nicht bis M gesucht werden, ansonsten ist spätestens bei M Schluss.

Damit man nicht immer und immer wieder dieselben Primzahlen ermitteln muss – die ändern sich ja schließlich nicht – ist es sinnvoll, diese nach einmaliger Ermittlung in eine Datei zu schreiben und beim nächsten Lauf dort auszulesen, womit dann nur noch Rechenaufwand für bis dato unbekannte Primzahlen-Bereiche nötig wird. Daraus ergibt sich folgende Routine (@primes_in_file enthält die bereits aus vorhergehenden Läufen bekannten Primzahlen):

sub fetch_primes_upto($ ) {
  my $max = shift;
  my $countknown = scalar(@primes_in_file);
  my $countunknown = 0;
  my $maxknown = 1;
  if ($countknown > 0) { $maxknown = $primes_in_file[$countknown-1] };        # supposing @primes is already sorted ... (as it was built by previous runs)
  if ($max > $maxknown) {
    for (my $i = $maxknown + 1; $i<=$max; $i++) {
      if ( (($countknown eq 0) && ($countunknown eq 0)) || is_prime($i) ) {
        push(@primes_in_file, $i);
        $countunknown++;
        if (($countunknown % 1000) eq 0) {    #save to file from time to time
          save_known_primes(@primes_in_file);
        }
      }
    }
  }
  my @return = ();
  #filter only upto required amount
  foreach (@primes_in_file) {
    if ($_ <= $max) { push(@return, $_); }
  }
  return \@return;
}

Weiter optimieren kann man, indem man nur die jeweils nächste benötigte Primzahl ermittelt, das sieht dann so aus:

sub fetch_smallest_prime_between ($$) {
  my ($current,$max) = @_;
  my $countknown = scalar(@primes_in_file);
  my $maxknown = 1;
  if ($countknown > 0) { $maxknown = $primes_in_file[$countknown-1] };        # supposing @primes is already sorted ... (as it was built by previous runs)
  if ($current <= $maxknown) {
    # next prime already known -> return
    foreach (@primes_in_file) {
      if ($_ >= $current) { return $_; }
    }
    return $maxknown;
  } else {
    #next prime unknown ->calculate
    if ($max > $maxknown) {
      for (my $i = $maxknown+1; $i<=$max; $i++) {
        if ( ($countknown eq 0) || is_prime($i) ) {
          push(@primes_in_file,$i);
          if ((scalar(@primes_in_file) % 1000) eq 0) {                #save to file from time to time
            save_known_primes(@primes_in_file);
          }
          return $i;
        }
      }
    }
  }
  return 0;
}

Das eigentliche Brute-Force-Durchprobieren der bekannten Primzahlen macht die Routine is_prime($):

sub is_prime($){
  my $test = shift;
  if (is_known_prime($test)) {
    return 1;
  } else {
    my @testprimes = @{fetch_primes_upto( POSIX::ceil(sqrt($test)) )};
    foreach (@testprimes) {
      if (($test % $_) eq 0) { return 0; } # not a prime
    }
    return 1;
  }
}

Diese bedient sich einer Abkürzung über die Routine is_known_prime, um zu prüfen, ob eine suche nach (womöglich noch unbekannten) Primfaktoren für diese Zahl überhaupt lohnt:

sub is_known_prime ($) {
  my $test = shift;
  my $maxknown = $primes_in_file[$#primes_in_file];
  if ($test > $maxknown) { return 0; }
  foreach (@primes_in_file) {
    if ($_ eq $test) { return 1; }
    if ($_ > $test) { last; } #done, no chance to find it now any more
  }
  return 0;
}

Damit fällt dann das oben beschriebene (rekursive) Coding für die eigentliche Faktorisierung relativ schmal aus:

sub factorize($ ) {
  my $base = shift;
  my @factors = ();
  if (is_known_prime($base)) {
    push(@factors, $base);
  } else {
    my $max = $base;
    my $current = fetch_smallest_prime_between(2,$max);
    while ( $current > 0) {
      if (($base % $current) eq 0) {
        push(@factors, $current);
        my $rest = $base / $current;
        if  ($rest > 1) {
          push(@factors, @{factorize($rest)});
        }
        last;
      }  else {
        $current++;
        $current = fetch_smallest_prime_between($current,$base);
      }
    }
  }
  return \@factors;
}

Gemeinsam mit Funktionen zum Datei-IO und ausgestattet mit ein wenig Logging ergibt sich ein schönes Perl-Modul, das besagte Aufgabe in wenigen Sekunden löst. (Die Lösung müsst ihr aber schon noch selbst ermitteln!)

Und hier die bisher von mir ermittelten 134852 Primzahlen (immernoch das Ergebnis einer ganzen Nacht) – natürlich ohne jede Garantie auf Richtigkeit! (Alle Dateien können gerne nach Herzenslust verwendet werden – wenn man sich denn den Spaß am selber implementieren rauben will – natürlich alles nur als beta-Version zu verstehen!)

No comments | Defined tags for this entry: open source, Perl, programming, Project Euler

Comments

No comments