Seite 1 von 2

Primzahl Prüfung Miller-Rabin Test

Verfasst: 30.12.2017 18:27
von NicknameFJ
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: Alles auswählen

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: Alles auswählen

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

Re: Primzahl Prüfung

Verfasst: 30.12.2017 19:12
von STARGÅTE
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: Alles auswählen

Debug isprime(11037271757) ; 11037271757 ist eine Primzahl! Es wird aber 0 zurückgegeben.
Selbst bei mehr Loops bleibt es Falsch.

Re: Primzahl Prüfung

Verfasst: 30.12.2017 19:32
von NicknameFJ
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: Alles auswählen

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

Re: Primzahl Prüfung

Verfasst: 30.12.2017 19:41
von GPI
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.

Re: Primzahl Prüfung

Verfasst: 30.12.2017 19:44
von STARGÅTE
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: Alles auswählen

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.

Re: Primzahl Prüfung Miller-Rabin Test

Verfasst: 31.12.2017 02:04
von NicknameFJ
@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

Re: Primzahl Prüfung Miller-Rabin Test

Verfasst: 31.12.2017 12:12
von ccode_new
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: Alles auswählen

#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

Re: Primzahl Prüfung Miller-Rabin Test

Verfasst: 31.12.2017 15:24
von NicknameFJ
@ccode_new: Der Ablauf bei Dir wirkt eleganter :allright:

Ich habe mir etwas davon abgeschaut.

Danke.

Code oben geändert.


NicknameFJ

Re: Primzahl Prüfung Miller-Rabin Test

Verfasst: 03.01.2018 22:37
von Helle
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: Alles auswählen

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

Re: Primzahl Prüfung Miller-Rabin Test

Verfasst: 03.01.2018 23:43
von NicTheQuick
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=nu ... +1+and+1e8