Mersenne number primality test program

[Sternenlicht]

Ein Programm zum testen von Zahlen auf Primalität mit Hilfe des Lucas Lehmer Tests der Form: 2^(P) - 1.

Marin Mersenne stellte ca: 1626 die mathematische Theorie auf, daß obige Formel, vorausgesetzt P ist eine Primzahl, wieder eine Primzahl ergibt. Für P = 3, 5, 7 stimmt dies auch, doch für 11 passt die Formel beweisbar schon nicht mehr.
Die Wahrscheinlichkeit nach dieser Formel dennoch eine Mersenne Primzahl zu finden, ist so hoch, daß sie heute noch angewendet, oder besser gesagt getestet wird.

Beispiel: 2^(5)-1 sei prime: 2*2*2*2*2 = 32: 32 - 1 = 31: is prime.
Beispiel: 2^(11)-1 sei prime: 2*2*2*2*2*2*2*2*2*2*2 = 2048: 2048 - 1 = 2047: is not prime: 2047 = 23 * 89

Die getesteten Zahlen befinden sich mittlerweile im millionenstelligem Bereich. Die größte bisher gefundene Form einer Mersenne Primzahl lautet am Samstag, den 16.02.2018: 2^(82,589,933)-1

May Be Prime!

[Pedaltreter]

Hallo Sternenlicht, ich habe auf Deiner Homepage gesehen, dass Du Dich für die Programmiersprache Pascal interessierst. Vielleicht kannst Du mir ja helfen. Vor einiger Zeit haben wir eine Seminarsaufgabe erhalten. Wir sollten ein Programm schreiben, welches nach dem Sieb von Eratosthenes Primzahlen ermittelt. Die Aufgabe zu lösen, war nicht schwer.

Nun möchte ich das Programm gern weiterentwickeln, stoße dabei jedoch sehr schnell an meine Kenntnisgrenzen. Mit Pascal, ob mit Delphi oder Lazarus, kann ich nur zuverlässig Zahlen bis zum 15 bzw. 18-stelligem Bereich berechnen. Int64 kann nur bis zu einem positiven Wert von 9223372036854775807 rechnen. Extended zwar bis 1.1E4932, allerdings bleibt die Genauigkeit ab der 18. oder 19. Stelle auf der Strecke.

Meine Frage an Dich, gibt es in Pascal eine Möglichkeit oder Lösung, mit viel größeren Zahlen zu rechnen? Ich möchte zum Beispiel 100 - stellige Zahlen auf Primalität untersuchen. Hast Du eine Idee?

Liebe Grüße, Pedaltreter

[Sternenlicht]

Hallo lieber Pedaltreter,

entschuldige, wenn meine Antwort so lange gedauert hat. Ich war reichlich beschäftigt.

Deine Anfrage habe ich vielleicht nicht so ganz verstanden. Primzahlen nach dem Aussiebungsverfahren von Eratosthenes zu ermitteln, ist die effektivste, bekannteste und zugleich schnellste Methode. Sie verlangt allerdings auch sehr viel Speicherplatz.

Du hast zum Beispiel die Zahl 9223372036854775807, der maximale Wert für eine int64 Variable, angeführt. Deine Anfrage hat mich dazu bewegt, einmal nach dem Verfahren von Eratosthenes alle Primzahlen in diesem Wertebereich (2 bis 9223372036854775807) zu ermitteln. Es war mir zunächst nicht möglich, alle Zahlen zu speichern, weil ein enormer Speicher bereitgestellt werden musste. Im ersten Schritt sollten alle Zahlen zwischen 2 und 9223372036854775807 aufgelistet und festgehalten werden. Das erfordert ein Array der Länge 9223372036854775807 und alles int64 Zahlen. Das sind mehrere Gigabyte. Also habe ich den Datentyp verändert, indem ich einfach mit boolschen Werten gearbeitet habe. Doch auch dies erforderte einen enorm hohen Speicherbedarf. Das Aussieben selbst, erfordert eine Rechnerleistung mit einer guten CPU nicht mehr als eine Minute.

Wenn ich nun bedenke, dass 9223372036854775807 eine 19-stellige Zahl ist, wirst Du nicht mit 100-stelligen Zahlen nach dem Sieb von Eratosthenes arbeiten wollen. Dies stößt unweigerlich an seine Grenzen. Um nun 100-stellige Zahlen auf Primalität untersuchen zu wollen, wird ein anderes Verfahren nötig sein. Vielleicht nach dem Prinzip der Probedivision, oder das Anwenden eines Lucas-Lehmer-Tests. Allerdings wird eine hohe Zahl auch entsprechend viel Zeit in Anspruch nehmen. Das Testen von millionenstelligen Zahlen dauert Tage bis mehrere Wochen.

Ich gebe Dir allerdings Recht, dass 100-stellige Zahlen unweigerlich ein Problem darstellen. Wie Du schon erwähnt hast, hört die Genauigkeit auch bei Extended - Variablen bei ca. 18 Stellen auf. Insofern ist der int64 - Wertebereich mit seinen 19 Stellen vielleicht sogar besser. Trotzdem mit 100-stelligen Zahlen zu arbeiten, bedeuten Rundungsfehler, auf der bei der Überprüfung einer Zahl auf Primalität keinen Verlass mehr wäre.

Etwas erfolgversprechender scheinen stattdessen String- Variablen zu sein, oder in Pascal die Nutzung von TList (Pro Stelle einen Listeneintrag). Ich persönlich würde der Einfachheit halber zunächst String-Variablen ausprobieren. Allerdings wird es dabei notwendig werden, für die einzelnen Berechnungsschritte eigene Routinen zu schreiben. Zum Beispiel stellt sich die Frage, wie mit String- Variablen dividiert wird.

Vielleicht konnte ich Dir mit meinen Ausführungen einige Ideen mitgeben, wie Du weiter verfahren möchtest?

Herzliche Grüße, das Sternenlicht...

[Pedaltreter]

Hallo Sternenlicht, vielen Dank! Deine Ausführungen waren sehr ideenreich. Wie ich weiter verfahren möchte? Tatsächlich möchte ich mit mehr Stellen im Dezimalsystem arbeiten, als der Wertebereich von int64 oder Extended mit einer Genauigkeit von 18 Stellen maximal zuläßt.

Deine Idee, TList oder String Variablen zu benutzen, erscheint mir als eine gute Idee. Allerdings weiss ich noch nicht, wie ich mit diesen Datentypen umgehen soll, wenn ich mathematische Berechnungen durchführen möchte. Arbeite ich mit Integer - Variablen so kann ich einfach Minuent - Subtrahend angeben, bei der Division Zahl1 / Zahl2. Wie soll ich jedoch mit String - Variablen umgehen, wenn diese zwar Zahlen enthalten, jedoch in Pascal als String ausgewertet werden. Hast Du da eine Idee?

Liebe Grüße, Pedaltreter

[Sternenlicht]

Lieber Pedaltreter,

selbstverständlich ist es nicht ganz so einfach, mit String - Variablen zu rechnen. Um mit solch einem Datentyp rechnen zu können, musst Du eine eigene Routine schreiben. Es ist ganz klar, denn der Compiler kann von sich aus nicht wissen, was Du machen möchtest. Eine String - Variable kann generell auch andere Zeichen enthalten, als nur Zahlen. Stell Dir vor, man würde dem Compiler anweisen, er solle 16mai2021nachApril mit 15 multiplizieren. Man sieht also sehr schnell, dass String - Datentypen auch als solche behandelt werden müssen.

Es besteht jedoch generell die Möglichkeit, auf Teile eines Strings zuzugreifen. Beispiel: Ich habe zwei Variablen:

Code: Alles auswählen

Var test : String;
      Zahl : Integer;
Begin
  test := '544';

So kann ich mit der Anweisung:

Code: Alles auswählen

Zahl := strToInt(test[2]);

der Variable Zahl den Wert 4 zuweisen. Mit der Anweisung:

Code: Alles auswählen

Zahl := strToInt(test);

hat die Variable Zahl vom Typ Integer, den Wert 544. usw...

Um mit String - Datentypen rechnen zu können, musst Du also, wie ich Eingangs erwähnte, eigene Routinen schreiben. Du möchtest Dir zur Aufgabe machen, so habe ich es jedenfalls verstanden, mit 100-stelligen Zahlen arbeiten, um diese auf Primalität zu untersuchen. An dieser Stelle möchte ich die einzelnen Rechenschritte nicht erwähnen. Sicher ist jedoch, dass man generell mit den vier Grundrechenarten (Addieren, Subtrahieren, Multiplizieren und Dividieren) auskommen sollte, um diese Aufgaben zu lösen. Im ersten Schritt würde ich, um die einfachste Rechneart dieser vier Grundrechenarten heraus zu suchen, zunächst das Addieren ausprobieren. Nimm dazu einfach einen Schreibblock, einen Stift und addiere schriftlich 54321 mit 12345 miteinander. So wirst Du schnell herausfinden, welcher Algorithmus nun in Pascal, bzw. generell in einer Programmiersprache anzuwenden wäre, um eine Additionsaufgabe zu lösen. Übertrage einfach die schriftliche Addition in einen Algorithmus in Pascal. Im zweiten Schritt würde ich genauso verfahren und das Subtrahieren ausprobieren. Generell können so alle vier Grundrechenarten nach dem Algorithmus der schriftlichen Addition, Subtraktion, Multiplikation und Division gelöst werden.

Vielleicht konnte ich Dir mit dieser Idee weiterhelfen?

Herzliche Grüße, das Sternenlicht...

[Pedaltreter]

Hallo Sternenlicht,

ja klar, ich habe verstanden, wie Du das meinst. Es ist, so wie ich es inzwischen überschlagen habe, recht aufwändig. Mit der schriftlichen Addition und Subtraktion bin ich noch klar gekommen und habe es mehr recht als schlecht in das Standard - Pascal von Lazarus übertragen können.
Mit der Multiplikation bin ich gar nicht klar gekommen, geschweige denn die Division. Vielleicht bin ich auch zu dumm.

Hast Du einen Tipp für mich?
Liebe Grüße, Pedaltreter

[Sternenlicht]

Lieber Pedaltreter,

na klar hätte ich einen Tipp für Dich. Ich werde mich einfach selbst mal daran probieren und Dir in den kommenden Tagen ein Beispielscript zukommen lassen.

Bis die Tage, das Sternenlicht...

[Sternenlicht]

Lieber Pedaltreter,

sorry, es hat einige Tage gedauert, weil ich mir zu den Skripten einige Gedanken machen musste. Außerdem habe ich mit Pascal schon längere Zeit nicht mehr gearbeitet. Dennoch hat es recht gut geklappt.

Weiter unten habe ich Dir die Unit "bigInt" angehängt. Diese kann einfach in ein anderes Projekt eingebunden werden. Die Unit "bigInt" wird einfach in einer Uses - Klausel definiert und mit der Variable "intGlobal" angesprochen. Beispiel: wurzel := intGlobal.WurzelInteger('144');

Einzig die Routine "divisionReal" habe ich aus zeitlichen Gründen nicht mehr fertig gestellt, weil ich mir außerdem dachte, dass Du auch gerne etwas probieren möchtest. Eigentlich kann die Routine genauso arbeiten wie die "divisionInteger" - Routine, nur dass die Kommastellen berücksichtigt werden müssen. Dies sollte jedoch kein größeres Problem mehr sein.

Außerdem muss bedacht werden, dass die Routine "wurzelInteger" immer zur nächsten Ganzzahl aufrundet. Zum Beispiel kann die Wurzel aus 7 keine Ganzzahl sein, daher ist das Ergebnis 3; denn 2*2 = 4 und 3*3 = 9.

Du kannst außerdem die Routinen in einer ausführbaren Datei ausprobieren. Das gesamte Lazarus - Projekt habe ich Dir in ein Zip - Archiv gepackt und hier zum Download aus der Cloud bereitgestellt.

Viel Spaß beim Ausprobieren und liebe Grüße,

das Sternenlicht...

---

Code: Alles auswählen

unit bigint;

{$mode objfpc}{$H+}

interface

uses Classes, SysUtils, strUtils;

Type TIntegerGlobal = class(TObject)
  Public
    Function wurzelInteger(Const Number: String): String;
    Function additionInteger(Const Number1, Number2: String): String;
    Function subtraktionInteger(Const number1: String; Const number2: String): String;
    Function multiplikationInteger(Const number1: String; Const number2: String): String;
    Function divisionInteger(number1, number2: String): String;
    Function isGreaterThan(number1, number2: String): Boolean;
    Function isGreaterOrEqualThan(number1, number2: String): Boolean;
End;

Var intGlobal : TIntegerGlobal;

implementation

Function TIntegerGlobal.isGreaterThan(number1, number2: String): Boolean;
Var Index    : Integer;
    greater : Boolean;
Begin
  greater := false;
  While (number1[1] = '0') And (Length(number1) > 1) Do Delete(number1, 1, 1);
  While (number2[1] = '0') And (Length(number2) > 1) Do Delete(number2, 1, 1);
  If Length(number1) < Length(number2)Then greater := true;
  If Length(number1) = Length(number2) Then
  Begin
    Index := 1;
    While (number1[Index] = number2[Index]) And (Index < Length(number1)) Do inc(Index);
    If number2[Index] > number1[Index] Then greater := true;
  end;
  isGreaterThan := greater;
end;

Function TIntegerGlobal.isGreaterOrEqualThan(number1, number2: String): Boolean;
Var Index          : Integer;
    greaterOrEqual : Boolean;
Begin
  Index := 1;
  greaterOrEqual := false;
  While (number1[1] = '0') And (Length(number1) > 1) Do Delete(number1, 1, 1);
  While (number2[1] = '0') And (Length(number2) > 1) Do Delete(number2, 1, 1);
  If Length(number1) < Length(number2)Then greaterOrEqual := true;
  If Length(number1) = Length(number2) Then
  Begin
    While (number1[Index] = number2[Index]) And (Index < Length(number1)) Do inc(Index);
    If number2[Index] >= number1[Index] Then greaterOrEqual := true;
  end;
  isGreaterOrEqualThan := greaterOrEqual;
end;

Function TIntegerGlobal.additionInteger(Const number1: String; Const number2: String) : String;
Var ergebnis, erg  : String;
    p1, p2, p3, p4 : int64;
    Index, Temp    : int64;
Begin
  ergebnis := '';
  temp := 0;
  p1 := Length(number1) + 1;
  p3 := Length(number2) + 1;
  While (p1 > 1) And (p3 > 1) Do
  Begin
    p1 := p1 - 1; p2 := p1;
    p3 := p3 - 1; p4 := p3;
    Index := 0;
    While (p1 > 1) And (p3 > 1) And (Index < 17) Do Begin Inc(Index); Dec(p1); Dec(p3); End;
    erg := intToStr(strToInt64(Copy(number1, p1, p2 - p1 + 1)) + strToInt64(Copy(number2, p3, p4 - p3 + 1)) + Temp);
    While(Length(erg) < p2 - p1 + 1)Do erg := '0' + erg;
    If(Length(erg) = p2 - p1 + 1) Then
    Begin
      ergebnis := erg + ergebnis;
      Temp := 0;
    End
    Else Begin
      ergebnis := Copy(erg, 2, Length(erg) - 1) + ergebnis;
      Temp := 1;
    End;
  End;
  While(Temp = 1) Do
  Begin
    If(p1 > 1)Then
    Begin
      Dec(p1);
      Temp := Temp + strToInt64(Copy(number1, p1, 1));
      If(Temp > 9) Then
      Begin
        Temp := Temp - 10;
        ergebnis := intToStr(Temp) + ergebnis;
        Temp := 1;
      End
      Else Begin
        ergebnis := intToStr(Temp) + ergebnis;
        Temp := 0;
      End;
    end;
    If(p3 > 1) Then
    Begin
      Dec(p3);
      Temp := Temp + strToInt64(Copy(number2, p3, 1));
      If(Temp > 9) Then
      Begin
        Temp := Temp - 10;
        ergebnis := intToStr(Temp) + ergebnis;
        Temp := 1;
      End
      Else Begin
        ergebnis := intToStr(Temp) + ergebnis;
        Temp := 0;
      End;
    End;
    If(p1 = 1) And (p3 = 1) And (Temp = 1) Then
    Begin
      Temp := 0;
      ergebnis := '1' + ergebnis;
    End;
  End;
  If(p1 > 1) Then
  Begin
    p2 := p1 - 1;
    p1 := 1;
    ergebnis := Copy(number1, 1, p2 - p1 + 1) + ergebnis;
  End;
  If(p3 > 1) Then
  Begin
    p4 := p3 - 1;
    p3 := 1;
    ergebnis := Copy(number2, 1, p4 - p3 + 1) + ergebnis;
  End;
  additionInteger := ergebnis;
end;

Function TintegerGlobal.subtraktionInteger(Const number1: String; Const number2: String): String;
Var ergebnis, erg  : String;
    p1, p2, p3, p4 : int64;
    Index, Temp    : int64;
Begin
  ergebnis := '';
  temp := 0;
  p1 := Length(number1) + 1;
  p3 := Length(number2) + 1;
  While (p1 > 1) And (p3 > 1) Do
  Begin
    p1 := p1 - 1; p2 := p1;
    p3 := p3 - 1; p4 := p3;
    Index := 0;
    While (p1 > 1) And (p3 > 1) And (Index < 17) Do Begin Inc(Index); Dec(p1); Dec(p3); End;
    erg := intToStr(strToInt64('1' + Copy(number1, p1, p2 - p1 + 1)) - strToInt64(Copy(number2, p3, p4 - p3 + 1)) - Temp);
    While(Length(erg) < p2 - p1 + 1)Do erg := '0' + erg;
    If(Length(erg) = p2 - p1 + 1) Then
    Begin
      ergebnis := erg + ergebnis;
      Temp := 1;
    End
    Else Begin
      ergebnis := Copy(erg, 2, Length(erg) - 1) + ergebnis;
      Temp := 0;
    End;
  End;
  While(Temp = 1) Do
  Begin
    Dec(p1);
    Temp := strToInt64('1' + Copy(number1, p1, 1)) - Temp;
    If(Temp > 9) Then
    Begin
      Temp := Temp - 10;
      ergebnis := intToStr(Temp) + ergebnis;
      Temp := 0;
    End
    Else Begin
      ergebnis := intToStr(Temp) + ergebnis;
      Temp := 1;
    End;
  End;
  If(p1 > 1) Then
  Begin
    p2 := p1 - 1;
    p1 := 1;
    ergebnis := Copy(number1, 1, p2 - p1 + 1) + ergebnis;
  End;
  subtraktionInteger := ergebnis;
end;

Function TIntegerGlobal.multiplikationInteger(Const number1: String; Const number2: String): String;
Var ergebnis, erg  : String;
    p1, p2, p3, p4 : int64;
    Index          : int64;
Begin
  p1 := Length(number1) + 1;
  ergebnis := '0';
  While (p1 > 1) Do
  Begin
    p1 := p1 - 1; p2 := p1;
    p3 := Length(number2) + 1;
    Index := 0;
    While (p1 > 1) And (Index < 8) Do Begin Inc(Index); Dec(p1); End;
    While (p3 > 1) Do
    Begin
      p3 := p3 - 1; p4 := p3;
      Index := 0;
      While (p3 > 1) And (Index < 8) Do Begin Inc(Index); Dec(p3); End;
      erg := intToStr(strToInt64(Copy(number1, p1, p2 - p1 + 1)) * strToInt64(Copy(number2, p3, p4 - p3 + 1)));
      For Index := p2 + 1 To Length(number1) Do erg := erg + '0';
      For Index := p4 + 1 To Length(number2) Do erg := erg + '0';
      ergebnis := intGlobal.additionInteger(ergebnis, erg);
    End;
  End;
  multiplikationInteger := ergebnis;
end;

Function TIntegerGlobal.divisionInteger(number1, number2: String): String;
Const maxInteger                 = '9223372036854775800';
Var pointer, Index, Zaehler      : Integer;
    Ergebnis, temp, Rest         : String;
Begin
  ergebnis := '0';
  Rest := '';
  While (number1[1] = '0') And (Length(number1) > 1) Do delete(number1, 1, 1);
  While (number2[1] = '0') And (Length(number2) > 1) Do delete(number2, 1, 1);
  For Index := 1 To Length(number2) - 1 Do Rest := Rest + number1[Index];
  For Index := Length(number2) To Length(number1) Do
  Begin
    Rest := Rest + number1[Index];
    temp := Rest;
    Zaehler := 0;
    While isGreaterOrEqualThan(number2, temp) Do
    Begin
      If isGreaterOrEqualThan(Rest, maxInteger) And isGreaterOrEqualThan(temp, maxInteger) Then
        temp := intToStr(strToInt64(Rest) - strToInt64(number2))
      Else
        temp := subtraktionInteger(Rest, number2);
      If isGreaterOrEqualThan('0', temp) Then inc(Zaehler);
      Rest := temp;
    end;
    If Zaehler > 0 Then
    Begin
      temp := intToStr(Zaehler);
      For pointer := Index + 1 To Length(number1) Do temp := temp + '0';
      ergebnis := intGlobal.additionInteger(ergebnis, temp);
    end;
  end;
  divisionInteger := Ergebnis;
end;

Function TIntegerGlobal.wurzelInteger(Const Number: String): String;
Var pos1, pos2                                    : Integer;
    quadrat, ergebnis, divident, rest, temp, zErg : String;
    found                                         : Boolean;
Begin
  quadrat := TrimLeftSet(Number, ['0']);
  If(Length(quadrat) Mod 2 = 1) Then quadrat := '0' + quadrat;
  pos1 := 1; Pos2 := 2;
  divident := quadrat[pos1] + quadrat[pos2];
  temp := '9';
  ergebnis := '';
  While isGreaterThan(divident, multiplikationInteger(temp, temp)) Do temp := subtraktionInteger(temp, '1');
  ergebnis := temp;
  rest := subtraktionInteger(Divident, multiplikationInteger(temp, temp));
  While (Pos2 < Length(quadrat)) Do
  Begin
    Inc(pos1, 2); Inc(pos2, 2);
    divident := rest + quadrat[pos1] + quadrat[pos2];
    temp := divisionInteger(divident, multiplikationInteger('20', ergebnis));
    found := false;
    Repeat
      zErg := multiplikationInteger(additionInteger(multiplikationInteger('20', ergebnis), temp), temp);
      If isGreaterOrEqualThan(zErg, divident) Then
        found := true
      Else
        temp := subtraktionInteger(temp, '1');
    Until found;
    temp := TrimLeftSet(temp, ['0']);
    If temp = '' Then temp := '0';
    ergebnis := ergebnis + temp;
    rest := subtraktionInteger(divident, zErg);
  end;
  wurzelInteger := ergebnis;
end;

end.

[Pedaltreter]

Hallo Sternenlicht, danke für den Code und Deine Hilfe. Die Skripte probiere ich direkt mal aus und gebe Dir danach eine Rückmeldung!
Danke!

[Pedaltreter]

Hallo Sternenlicht,

noch einmal herzlichen Dank für Deine durchaus kompetente Hilfe. Zu Deinen Skripten gibt es eine positive und negative Nachricht. Die positive Nachricht ist, dass das Skript richtig gut und zuverläßig zu gebrauchen ist. Die negative Nachricht hingegen, dass die Ausführungsgeschwindigkeit arg mit diesen Routinen nachläßt. Aber das ist nunmal so mit großen Zahlen. Je größer sie sind, desto länger dauern die Berechnungen.

Oder hättest Du auch da eine Idee?

Liebe Grüße, Pedaltreter