PureBoard
http://forums.purebasic.com/german/

Primzahl Prüfung Miller-Rabin Test
http://forums.purebasic.com/german/viewtopic.php?f=8&t=30523
Seite 1 von 2

Autor:  NicknameFJ [ 30.12.2017 18:27 ]
Betreff des Beitrags:  Primzahl Prüfung Miller-Rabin Test

Hallo zusammen,

für ein Projekt habe ich eine schnelle (die Geschwindigkeit war zumindest für meine Zwecke mehr als ausreichend) Prüfung gebraucht ob eine Zahl (bis ca. 2^62) also positive Quad eine Primzahl ist oder nicht.

Ich habe mich mit dem Miller-Rabin Test vertraut gemacht und die Procedure IsPrime() geschrieben.

Da der Miller-Rabin Test nicht sicher sagt ob eine Zahl eine Primzahl ist, jedoch sicher sagt wenn es keine Primzahl ist, wurde der Parameter "Loop" mit eingebaut mit der an der Wahrscheinlichkeitsschraube gedreht werden kann. Damit kann die Wahrscheinlichkeit einer Fehlaussage bis zur gewünschten kleinen Wahrscheinlichkeit beeinflusst werden. Ohne Angabe ist im schlechtesten Fall eine maximale Fehlerwahrscheinlichkeit von ca. 0,0001 % gegeben. Je höher der Loop gesetzt wird desto unwahrscheinlicher ist eine Fehleinschätzung.

Die Procedure gibt #False (wenn keine Primzahl) oder #True (wenn Primzahl mit der beschriebenen Fehlereinschätzung) zurück.

Wer es gebrauchen kann, kann es nutzen (Ich wollte auch mal etwas an die Community zurückgeben).

Für Fehlerhinweise wäre ich dankbar.


NicknameFJ


Sollen nur Zahlen bis QUAD geprüft werden ist die untenstehende Variante von Helle deutlich schneller da unser Assembler-Spezialist Helle dort den BigInt Teil durch 64/128-Bit Routinen ersetzt hat.

@Helle: Danke dir dafür.

Code:
EnableExplicit

; Pfad zur LIB anpassen

IncludeFile "O:\Includes\BigInt\BigInt.pbi"


; Nutzung der BIGInt Lib von wilbert aus dem englischen Forum
; http://www.purebasic.fr/english/viewtopic.php?f=12&t=61309



Procedure isPrime(Zahl.q, Loop = 10)
 
  ; Ist Z keine Primzahl gibt die Proc.                   0 zurück
 
  ; Ist Z eine Primzahl  gibt die Proc.                   1 zurück
  ;    Rückgabe 1 kann auch bedingt durch den Algoritmus
  ;    Miller-Rabin eine Zahl anzeigen die nicht
  ;    prim ist. Die Fehlerwahrscheinlichkeit kann durch
  ;    duch Erhöhung des Parameters Loop weiter gesenkt
  ;    werden.
 
  Protected d.q, r.q, a.q, x.q, lp.i, CheckBZE.i
  Protected.BigInt::BigInt BX, BA, BD, BZahl, BEins, BZahl_1, n0, BZero
 
  If Zahl = 2
    ProcedureReturn #True
   
  ElseIf Zahl < 2 Or Zahl % 2 = 0
    ProcedureReturn #False
   
  EndIf           
 
  If Loop < 1 : Loop = 1 : EndIf
 
 
  BigInt::SetValue(BEins,1)   
  BigInt::SetValue(BZero,0)
  BigInt::SetValue(BZahl_1,Zahl - 1)
 
  d = Zahl - 1
  r = 0
 
  While d % 2 = 0   
    d = d / 2
    r  + 1   
  Wend 
 
 
  While Loop > 0
   
    Loop -1
   
    a = Random(Zahl-1,2) 
   
   
    BigInt::SetValue(BA,a)
    BiGInt::SetValue(BD,d)
    BiGInt::SetValue(BZahl,Zahl)

    BIGInt::ModPow(BX,BA,BD,BZahl)
   
    If BigInt::Compare(BX,BEins) = 0 Or BigInt::Compare(BX,BZahl_1) = 0
      Continue
    EndIf
   
    lp = r
   
    While lp > 0
     
     
     BigInt::ModMul(BX,BX,BZahl)
     
      If BigInt::Compare(BX,BEins) = 0
        ProcedureReturn #False
      EndIf
     
      CheckBZE = BigInt::Compare(BX,BZahl_1)
     
      If  CheckBZE = 0
        Break               
      EndIf
     
      lp - 1       
     
    Wend
   
    If  CheckBZE <> 0
      ProcedureReturn #False
     
    EndIf
   
   
  Wend
 
  ProcedureReturn #True
 
EndProcedure



CompilerIf #PB_Compiler_IsMainFile
 
Define Zahl.q

For Zahl = 2 To 1000
  If isprime(Zahl)
    Debug  Str(Zahl) + " ist prim"
  EndIf
 
Next

CompilerEndIf




[EDIT]
Nachfolgend noch eine Variante die komplett mit BigInt arbeitet. Der Größe der zu prüfenden Zahl ist damit theoretisch nach oben hin keine Grenze gesetzt.

In der BigInt - Lib ist die Konstante #BigIntBit auf 2048 gesetzt. Daher keine Zahlen größer als 2^1024 -1 prüfen da ansonsten innerhalb dere Berechnung ein Überlauf auftreten kann und das Ergebnis dann nicht mehr stimmt. Ggfs. muss die Konstante innerhalb der Lib auf einen höheren Wert gesetzt werden.

Code:
EnableExplicit

; Pfad zur LIB anpassen

IncludeFile "O:\Includes\BigInt\BigInt.pbi"


; Nutzung der BIGInt Lib von wilbert aus dem englischen Forum
; http://www.purebasic.fr/english/viewtopic.php?f=12&t=61309



Procedure isPrime(*pZahl.BigInt::BigInt, Loop.i = 10)
 
  ; wird 0 zurückgegeben ist die zu testende Zahl keine Primzahl
 
  ; wird 1 zurückgegeben ist die zu testende Zahl eine Primzahl
  ;    Rückgabe 1 kann auch bedingt durch den Algoritmus
  ;    Miller-Rabin eine Zahl anzeigen die nicht
  ;    prim ist. Die Fehlerwahrscheinlichkeit kann durch
  ;    die Erhöhung des Parameters Loop weiter gesenkt
  ;    werden.
 
  Protected r.q, lp.i, CheckBZE.i, length, BitMask = $8000000000000000
  Protected.BigInt::BigInt BX, BA, BD, BZahl, BEins, BZahl_1, BZero, BZwei, BReminder
 
 
  BigInt::Assign(BZahl,*pZahl)   
 
  BigInt::SetValue(BZwei,2)
  BigInt::SetValue(BEins,1)   
  BigInt::SetValue(BZero,0)
 
 
  If BigInt::Compare(*pZahl,BZwei) = 0
    ProcedureReturn #True
  EndIf         
 
  length = BigInt::#BigIntBits >> 6
 
  If length       
    If *pZahl\q[length-1] & BitMask
      ProcedureReturn #False
    EndIf
  EndIf
 
  If *pZahl\q[0] & 1 = 0
    ProcedureReturn #False
  EndIf
 
 
  If Loop < 1 : Loop = 1 : EndIf
 
  BigInt::Assign(BZahl_1,BZahl)
  BigInt::Subtract(BZahl_1,BEins)
 
  BigInt::Assign(BD,BZahl_1)
 
  r = 0
 
  While BD\q[0] & 1 = 0   
    BigInt::Divide(BD,BD,BZwei,BReminder)
    r + 1
  Wend
 
 
  While Loop > 0
   
    Loop - 1
   
    length = 0
   
    While length = 0 Or (length = 1 And BA\q[0] < 2)
     
      RandomData(@BA,BigInt::#BigIntBits >> 6)   
     
      If BA\q[(BigInt::#BigIntBits >> 6)-1] & BitMask
        BigInt::NEG(BA)
      EndIf
     
      BigInt::Divide(BA,BA,BZahl_1,BReminder)
      BigInt::Assign(BA,BReminder)
      length = BigInt::NumberSize(BA)
     
    Wend
   
   
    BigInt::ModPow(BX,BA,BD,BZahl)
   
    If BigInt::Compare(BX,BEins) = 0 Or BigInt::Compare(BX,BZahl_1) = 0
      Continue
    EndIf
   
    lp = r       
   
    While lp > 0
     
      BigInt::ModMul(BX,BX,BZahl)
     
      If BigInt::Compare(BX,BEins) = 0
        ProcedureReturn #False
      EndIf
     
      CheckBZE = BigInt::Compare(BX,BZahl_1)
     
      If  CheckBZE = 0
        Break               
      EndIf
     
      lp -1
     
    Wend
   
    If CheckBZE <> 0
      ProcedureReturn #False     
    EndIf
   
  Wend
 
  ProcedureReturn #True
 
EndProcedure


CompilerIf #PB_Compiler_IsMainFile
 
 
  Procedure.s GetDec(*Original.BigInt::BigInt)
   
    Protected Result.s, BitMask.q, Sign, length
    Protected.BigInt::BigInt Teiler, r, n
   
    BitMask = $8000000000000000   
    length = BigInt::#BigIntBits >> 6
   
    BigInt::Assign(n,*Original)       
    BigInt::SetValue(Teiler,10)
       
    If n\q[length-1] & BitMask
      BigInt::Neg(n)
      Sign = #True
    EndIf
   
    Repeat
      BigInt::Divide(n, n, Teiler, r)
      Result = Str(r\q[0]) + Result           
    Until BigInt::isZero(n) <> 0
   
    If Sign
      Result= "-" + Result 
    EndIf
   
    ProcedureReturn Result
   
  EndProcedure
 
 
 
 
  Define.BigInt::BigInt n0
 
  BigInt::SetHexValue(n0,"7ffffffffffffffffffffffffff")
 
 
  ; #BigIntBits ist in der BigInt-Lib auf 2048 gesetzt
 
  ; daher maximal  2^1024-1     = 1797693134862315907729305190789024733617976978
  ; übergeben                     9423065727343008115773267580550096313270847732
  ;                               2407536021120113879871393357658789768814416622
  ;                               4928474306394741243777678934248654852763022196
  ;                               0124609411945308295208500576883815068234246288
  ;                               1473913110540827237163350510684586298239947245
  ;                               938479716304835356329624224137215
 
  ; für größere Werte muss die #BigIntBits-Konstante erhöht werden
 
 
  Debug "zu prüfende Zahl: $" + BigInt::GetHex(n0) + "   = " + GetDec(n0)                                                     
  If  isPrime(n0)
    Debug "ist eine Primzahl"
  Else
    Debug "ist keine Primzahl"
  EndIf
 
CompilerEndIf

Autor:  STARGÅTE [ 30.12.2017 19:12 ]
Betreff des Beitrags:  Re: Primzahl Prüfung

Hallo NicknameFJ,

es wäre gut, wenn du Sachen wie Protected und Typendefinitionen verwendest, sonst gibt es bei andere ggf. große Probleme.
Vorallem wenn man es in x86 kompiliert, hast du überall Longs und nicht Quads, wie du wolltest.

Zum Code selbst:
Ich kenne mich mit dem Miller-Rabin Test nicht aus, aber ein Test mit Random() halte ich für riskannt, da es nun mal nie 100% Sicherheit gibt.
Ich habe auch mal mit ein paar Nummern getestet und dabei auf Probleme gestoßen:
Code:
Debug isprime(11037271757) ; 11037271757 ist eine Primzahl! Es wird aber 0 zurückgegeben.

Selbst bei mehr Loops bleibt es Falsch.

Autor:  NicknameFJ [ 30.12.2017 19:32 ]
Betreff des Beitrags:  Re: Primzahl Prüfung

Hallo Stargate,

Typdeklarationen usw. baue ich noch ein. Da hast du Recht, es ist wichtig.

Ob 1103727757 eine Primzahl ist weiß ich nicht, aber:

sofern der Algo richtig funktioniert und er #False zurückgibt nutzt es nichts die Loops zu erhöhen. Dann kann sich nichts mehr ändern. Nimm mal 221 (ist keine Primzahl) als Zahl und lasse es mit Loop 1 mehrere hundert male durchlaufen und dann mal mit Loop 10 und zähle die Fälle in denen die Procedure #true zurückgibt. Wenn #false zurückgegeben wird ist es keine Primzahl (so sollte zumindest der Miller-Rabin Algo sein), nur bei #true ist es unsicher, dafür aber relativ schnell.

Dem Algorithmus sollte eine zufällige Zahl eingespielt werden, dieses Prinzip liegt dem Algo zugrunde wenn ich es richtig verstanden habe. Es gibt lt. Wikipedia auch eine deterministische Variante wo bestimmte Primzahlen zur Ermittlung eingespielt werden. Aber da bin ich erst später dahinter gekommen.

Aus Wikipedia:

Code:
Oft wird a zufällig gewählt, der MRT zählt in dieser Form zur Klasse der Monte-Carlo-Algorithmen. Durch wiederholtes Testen mit verschiedenen a kann man die Wahrscheinlichkeit eines Irrtums beliebig klein machen. Es gibt auch deterministische Varianten des MRT, bei denen durch geeignete Wahl der a ein Irrtum ausgeschlossen wird.


Aber ich schaue mir dass alles noch mal an.

Danke für dein Feedback.


Grüße

NicknameFJ

Autor:  GPI [ 30.12.2017 19:41 ]
Betreff des Beitrags:  Re: Primzahl Prüfung

Meine Lösung wäre wohl, eine Datenbank mit allen Primzahlen anzulegen und dann nachschauen, ob sie drin ist.
Für was genau brauchst du diese Funktion? Oft findet man eine andere bessere Lösung, wenn man weiß, um was es geht.

Autor:  STARGÅTE [ 30.12.2017 19:44 ]
Betreff des Beitrags:  Re: Primzahl Prüfung

Zitat:
und er #False zurückgibt nutzt es nichts die Loops zu erhöhen

Ah ok, dann gibt es wohl ein Problem mit dem Algo. --> Is 11037271757 prime

Ich denke das Problem ist hier:
Code:
For lp = c To 62
    Result = Result * Result
    Result = Result % m

Hier wird es ein Überlauf geben, wenn Result (vor der Quadrierung) schon mehr als eine Long einnimmt, was ja durchaus passiert, weil du Random() theoretisch auch recht groß werden kann für große Zahlen, dann kommt natürlich Müll raus.

ggf muss du dein Test maximal 2^31 beschränken.

Autor:  NicknameFJ [ 31.12.2017 02:04 ]
Betreff des Beitrags:  Re: Primzahl Prüfung Miller-Rabin Test

@GPI: Die Procedure funktioniert ja bis auf den von Stargate aufgedeckten OverFlow Error, der jetzt bei Auftreten -1 zurückgibt, wunderbar. Auf die Nutzung einer Datenbank mit allen Primzahlen habe ich ganz bewusst verzichtet.

@Stargate: Danke Dir für das Aufzeigen des Overflow Error. An einen möglichen Überlauf an dieser Stelle habe ich leider nicht gedacht.
Nur das a vor Aufruf der BinExpMod Procedure auf 2^31-1 zu begrenzen genügt leider nicht da auch der Exponent je nach zu prüfender Zahl ziemlich groß werden kann.

[EDIT:]

Ich verwende jetzt die BigInt-Lib von wilbert aus dem englischen Forum. Somit dürften die Overflow Probleme behoben sein.


Grüße

NicknameFJ

Autor:  ccode_new [ 31.12.2017 12:12 ]
Betreff des Beitrags:  Re: Primzahl Prüfung Miller-Rabin Test

Hallo NicknameFJ!

NicknameFJ hat geschrieben:
Evtl. müsste die Procedur auf eine BigInt - Lib. umgestellt werden um diese Überläufe zu vermeiden.


Ein erhöhen von AllocateMemory() bringt dir auch nur dann etwas wenn du es schaffst den kompletten Algorithmus zu ändern (Aufteilung sehr großer Zahlen, sonst: "Numerical overflow: too many digits.") und du musst dir auch eine eigene Pow()-Funktion (Return = Double = zu klein) implementieren.

Test:
Code:
#Verbindung = 0
#Wahrscheinlich_PRIMZAHL = 1
 
Procedure Miller_Rabin(n, k)
  Protected d=n-1, s, x, r
  If n=2
    ProcedureReturn #Wahrscheinlich_PRIMZAHL
  ElseIf n%2=0 Or n<2
    ProcedureReturn #Verbindung
  EndIf
  While d%2=0
    d/2
    s+1
  Wend
  While k>0
    k-1
    x=Int(Pow(2+Random(n-4),d))%n
    If x=1 Or x=n-1: Continue: EndIf
    For r=1 To s-1
      x=(x*x)%n
      If x=1
        ProcedureReturn #Verbindung
      EndIf
      If x=n-1
        Break
      EndIf
    Next
    If x<>n-1
      ProcedureReturn #Verbindung
    EndIf
  Wend
  ProcedureReturn #Wahrscheinlich_PRIMZAHL
EndProcedure


"Der Miller-Rabin-Test ist ein Monte-Carlo-Algorithmus zum Testen, ob Zahlen prim sind oder nicht."

Anbei:
Schönes neues Jahr an Alle!

Bild

Autor:  NicknameFJ [ 31.12.2017 15:24 ]
Betreff des Beitrags:  Re: Primzahl Prüfung Miller-Rabin Test

@ccode_new: Der Ablauf bei Dir wirkt eleganter :allright:

Ich habe mir etwas davon abgeschaut.

Danke.

Code oben geändert.


NicknameFJ

Autor:  Helle [ 03.01.2018 22:37 ]
Betreff des Beitrags:  Re: Primzahl Prüfung Miller-Rabin Test

Da der Eingabewert bei diesem Programm durch PB auf Quad begrenzt ist, kann der Multiplikations-Überlauf auch ohne externe Lib vermieden werden. Natürlich mit Assembler :mrgreen: ! Wenn die CPU zwei 64-Bit-Werte multipliziert, steht das Ergebnis 128-bittig zur Verfügung. Dies kann für Modulo (also Division) prima genutzt werden. Sieht dann mit NicknameFJs (neuem) Code so aus:
Code:
EnableExplicit

Procedure isPrime(Zahl.q, Loop = 10)   ;64-Bit
 
  ; Ist Z keine Primzahl gibt die Proc.                   0 zurück
 
  ; Ist Z eine Primzahl  gibt die Proc.                   1 zurück
  ;    Rückgabe 1 kann auch bedingt durch den Algoritmus
  ;    Miller-Rabin eine Zahl anzeigen die nicht
  ;    prim ist. Die Fehlerwahrscheinlichkeit kann durch
  ;    duch Erhöhung des Parameters Loop weiter gesenkt
  ;    werden.
 
  Protected.q d, r, a, lp, BX, BA, BD, BZahl, BEins, BZahl_1, BZero
  ;Protected.BigInt::BigInt BX, BA, BD, BZahl, BEins, BZahl_1, BZero
 
  If Zahl = 2
    ProcedureReturn #True
   
  ElseIf Zahl < 2 Or Zahl % 2 = 0
    ProcedureReturn #False
   
  EndIf           
 
  If Loop < 1 : Loop = 1 : EndIf
 
 
  ;BigInt::SetValue(BEins,1)   
  ;BigInt::SetValue(BZero,0)
  ;BigInt::SetValue(BZahl_1,Zahl - 1)
  BEins = 1
  BZero = 0
  BZahl_1 = Zahl - 1

  d = Zahl - 1
  r = 0
 
  While d % 2 = 0   
    d = d / 2      ;d=Faktor
    r + 1          ;r=Exponent
  Wend
 
 
  While Loop > 0
   
    Loop -1
   
    a = Random(Zahl-1,2)
   
   
    ;BigInt::SetValue(BA,a)
    ;BiGInt::SetValue(BD,d)
    ;BiGInt::SetValue(BZahl,Zahl)
    BA = a                   ;könnte man sich auch schenken
    BD = d
    BZahl = Zahl

    ;BIGInt::ModPow(BX,BA,BD,BZahl)
    !mov r8,1                ;Startwert BX
    !mov r9,[p.v_BD]
    !bsr rcx,r9              ;höchstes gesetztes Bit
    !add rcx,1               ;Schleifenzähler
    !mov r10,[p.v_BA]
    !mov r11,[p.v_BZahl]
   !.L1:
    !shr r9,1
   !jnc @f                   ;das aktuelle Bit0 war nicht gesetzt
    !mov rax,r8
    !mul r10                 ;rax*r10, Produkt in rdx:rax
    !div r11                 ;rdx:rax / r11
    !mov r8,rdx              ;Modulo=rdx (Remainder)
   !@@:
    !mov rax,r10
    !mul rax
    !div r11
    !mov r10,rdx
    !sub rcx,1
   !jnz .L1
    !mov [p.v_BX],r8

    ;If BigInt::Compare(BX,BEins) = 0 Or BigInt::Compare(BX,BZahl_1) = 0
    If BX = BEins Or BX = BZahl_1
      Continue
    EndIf
   
    lp = r
   
    While lp > 0
     
     
     ;BigInt::ModMul(BX,BX,BZahl)
      !mov rax,[p.v_BX]
      !mov rcx,[p.v_BZahl]
      !mul rax
      !div rcx
      !mov [p.v_BX],rdx

      ;If BigInt::Compare(BX,BEins) = 0
      If BX = BEins
        ProcedureReturn #False
      EndIf
     
      ;CheckBZE = BigInt::Compare(BX,BZahl_1)
     
      ;If  CheckBZE = 0
      If BX = BZahl_1
        Break               
      EndIf
     
      lp - 1       
     
    Wend
   
    ;If  CheckBZE <> 0
    If BX <> BZahl_1
      ProcedureReturn #False
     
    EndIf
   
   
  Wend
 
  ProcedureReturn #True
 
EndProcedure



CompilerIf #PB_Compiler_IsMainFile
 
  Define.q Zahl, ZahlMin, ZahlMax, TA, TE, Res1, Res2, Prim, i

  ;für Tests
  ;Zahl = 15                    ;3*5
  ;Zahl = 91                    ;7*13
  ;Zahl = 561                   ;1.Carmicheal-Zahl  3*11*17
  ;Zahl = 65537
  ;Zahl = 11037271757
  ;Zahl = 2305843009213693951 ;Mersenne M9=2305843009213693951
  Zahl = 9223372036854775783 ;größtstmögliche hier zu testende Primzahl

  Res1 = isPrime(Zahl)

  ZahlMin = 2
  ZahlMax = 1E6
  ;Test, wieviel Primzahlen im gewählten Bereich gefunden wurden
  ;hier geht es nicht primär um die Zeit (da gibt es viel schnellere Routinen), sondern ob die Anzahl stimmt
  TA = ElapsedMilliseconds()
  For i = ZahlMin To ZahlMax
    Res2 = isPrime(i)
    If Res2
      Prim + 1
    EndIf
  Next
  TE = ElapsedMilliseconds() - TA

  MessageRequester("isPrime-Test-ohne Lib", "Ergebnis für Zahl: " + Str(Res1) + #LFCR$ + "Zwischen ZahlMin und ZahlMax wurden " + Str(Prim) + " Primzahlen gefunden." + #LFCR$ + "Dafür benötigte Zeit: " + Str(TE) + "ms")

CompilerEndIf


Bringt auch schönen Speed-Gewinn (Faktor mitunter zweistellig). Die 2 Assembler-Parts können auch als Macro/Procedure ausgeführt werden.

Gruß
Helle

Autor:  NicTheQuick [ 03.01.2018 23:43 ]
Betreff des Beitrags:  Re: Primzahl Prüfung Miller-Rabin Test

Cool, hab es mal getestet mit ZahlMax = 1E8. Ich hab 5761455 Primzahlen in 39676 ms gefunden. Müsste man jetzt noch parallelisieren. Scheint auch zu stimmen: http://www.wolframalpha.com/input/?i=number+of+primes+between+1+and+1e8

Seite 1 von 2 Alle Zeiten sind UTC + 1 Stunde [ Sommerzeit ]
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group
http://www.phpbb.com/