Hallo,
ich möchte euch bitten, mir als Programmierneuling Rückmeldung zu meinem Code zu geben. Der Hintergrund ist folgender:
Ich beschäftige mich in meinem Studium mit Zahlen der Form n = h * 2^k + 1, wobei h ungerade und kleiner 2^k ist. Nun kann man prüfen, ob die Zahl n eine Primzahl ist:
1) Gehe die Primzahlen sukzessive durch und bilde für jede Primzahl D das Jacobisymbol (n, D).
2) Ist das Jacobisymbol 1, gehe eine Primzahl weiter
3) Ist es 0, ist n keine Primzahl
4) Ist es -1, berechne die Kongruenz D^((n-1)/2) modulo n. Kommt dort -1 raus, ist n prim. Sonst nicht
Ich würde meinen Code mal päckchenweise hier posten in der Reihenfolge wie ich auch programmiert habe. Hoffe das ist so ok.
Hier aber auch als Spoiler der gesamte Code:
Ich habe mit einem anderen Skript die Primzahlen bis 10^7 gesiebt und in der Datei abgelegt. Dazu habe ich mir die beiden Funktionen listsave und listload geschrieben. Falls jemand den Code testen will, hänge ich die Datei gerne an.
import time
import matplotlib.pyplot as plt
import pickle
def listsave(liste, datei):
with open(datei, 'wb') as file:
pickle.dump(liste, file)
def listload(datei):
with open(datei, 'rb') as file:
return pickle.load(file)
primzahlen = listload('primes_up_to_10_7')Dann berechne ich das Jacobisymbol. Ist jetzt vielleicht zu mathematisch, das hier zu erklären. Aber die Werte die dort herauskommen habe ich auch mit denen von sympy.ntheory.residue_ntheory.jacobi_symbol(m, n) verglichen. Ob letztere schneller ist als meine Version ist mir nicht so wichtig, da ich das hier erstmal verstehen muss.
def jacobisymbol(a, n):
if n == 1 or n % 2 == 0:
return
if a == 2:
if n % 8 == 1 or n % 8 == 7:
return 1
else:
return -1
if a == -1:
return pow(-1, (n-1)//2)
if a == 1:
return 1
if a >= n:
return jacobisymbol(a % n, n)
if a == 0:
return 0
if a % 2 == 0:
k = 0
while (a % 2 == 0):
k += 1
a = a // 2
return jacobisymbol(2, n) ** k * jacobisymbol(a, n)
faktor = pow(-1, ((n-1) // 2) * (a-1) // 2)
return faktor * jacobisymbol(n, a)Der folgende Test ist Punkt 4 aus der Aufzählung oben. Hier ist mir eine Sache wichtig:
Anstatt das Jacobisymbol (h * 2^k + 1 , D) zu berechnen, kann man reduzieren und zwar:
h % D statt h
k % (D-1) statt k
Für statistische Zwecke lasse ich mir nicht nur das Ergebnis ausgeben, ob n eine Primzahl ist. Sondern auch die Primzahl D, mit der das Jacobisymbol -1 ergibt.
Sollten alle Primzahl in der vorgelegten List 1 ergeben, so ist die Liste einfach nicht lang genug. Dann gebe -1 zurück.
def bosma(h, k, primes):
for D in primes:
jacobi = jacobisymbol( (h % D)*2**(k % (D-1)) + 1, D)
exp = h*2**(k-1)
if jacobi == 0:
return [False, D]
if jacobi == -1:
n = 2*exp+1
if pow(int(D), exp, n) == n-1:
return [True, D]
return [False, D]
return -1Das funktioniert auch soweit schonmal. Dann habe ich mich an die Auswertung gemacht. Ich möchte einen Exponenten k vorgeben und für verschiedene h prüfen, ob h*2^k+1 prim ist oder nicht.
Dazu prüfe ich mit dem oben beschrieben Verfahren die Tupel (h, k) und lege das Ergebnis als Tupel (h, k, True/False) in einer Liste ab.
Dann akkumuliere ich die Werte in der Liste und plotte die Verteilung der Primzahlen.
Da ich schon Kommentare im Code habe, gebe ich diesen einfach mal komplett wieder.
Bitte einfach schimpfen wenn das nicht korrekt ist 
# Es wird ein Vektor angelegt, der laenge mal den Eintrag "False" enthält
def vektor_anlegen(laenge, belegung = False):
vektor = laenge * [False]
return vektor
# Diese Funktion bekommt eine Liste, die boole'sche Variablen enthält.
# Dann wird der Vektor mit akkumulierten Werten zurückgeliefert.
# Dies geht wohl auch mit itertools.accumulate()
def anzahl_ermitteln(vektor):
ywerte = []
anzahl = 0
for eintrag in vektor:
if eintrag == True:
anzahl += 1
ywerte.append(anzahl)
return ywerte
# Diese Funktion erhält eine Liste mit Tupeln der Form (h, k, True/False)
# aus der Funktion intervall_pruefen.
# Dann werden die einzelnen h in die Liste xwerte gesteckt,
# und der jeweils zugehörige boole'sche Wert in die Liste y_rohdaten.
# Aus dieser Liste werden dann über anzahl_ermitteln die akkumulierten y-Werte erzeugt
# und zurückgeliefert.
def liste_vorbereiten(vektor):
y_rohdaten = []
xwerte = []
for eintrag in vektor:
xwerte.append(eintrag[0])
y_rohdaten.append(eintrag[2])
ywerte = anzahl_ermitteln(y_rohdaten)
return [xwerte, ywerte]
# Diese Funktion prüft die Zahlen der Form h*2^k + 1 mit dem Bosma-Test.
# Dabei werden alle ungeraden h zwischen unter- und obergrenze durchlaufen.
# Die einzelnen Ergebnisse werden als Tupel (h, k, True/False) in der Liste ergebnis abgelegt.
def intervall_pruefen(untergrenze, obergrenze, exponent, liste, ausgabe=False):
if untergrenze % 2 == 0:
untergrenze += 1
print("Untergrenze auf {} angepasst".format(untergrenze))
if obergrenze == 0:
obergrenze = 2**exponent - 1
if obergrenze >= 2**exponent:
obergrenze = 2**exponent -1
print("Obergrenze auf {} angepasst".format(obergrenze))
ergebnis = []
anzahl = (obergrenze - untergrenze)//2
if ausgabe:
for h in range(untergrenze, obergrenze+1):
print((h-1)//2 / anzahl, end = '\r')
B = bosma(h, exponent, liste)
ergebnis.append((h, exponent, B[0]))
return ergebnis
# Diese Stelle wird nur erreicht, wenn ausgabe=False gesetzt wurde
for h in range(untergrenze, obergrenze+1, 2):
B = bosma(h, exponent, liste)
ergebnis.append((h, exponent, B[0]))
return ergebnis
def plotten(untergrenze, obergrenze, exponenten):
for k in exponenten:
print("Aktueller Exponent:", k)
vektor = intervall_pruefen(untergrenze, 0, k, primzahlen, True)
plotvektor = liste_vorbereiten(vektor)
plt.plot(plotvektor[0], plotvektor[1])
plt.show()
h_unter = 1
h_ober = -1
exponenten_unter = 2
exponenten_ober = 12
exponenten = list(range(exponenten_unter,exponenten_ober+1, 1))
plotten(h_unter, h_ober+1, exponenten)Es funktioniert auch wunderbar und ich bekomme schöne plots dabei raus. Aber Verbesserungen sind immer wichtig
Vielen Dank für eure Hilfe!
Hier einmal das Skript als Datei: statistik_in_einer_datei.py
Außerdem: die Primzahlliste. Gezippt, da sonst "ungültige Dateiendung" kommt.
