Bitte Kritik zu meinem Code

neuernutzer · May 28th 2023

Hi,

ich habe leider jetzt erst gesehen, dass du ein paar Punkte von Dennis89 aufgegriffen hast. Ich lasse sie oben trotzdem mal stehen.

Quote from __blackjack__

Bei `intervall_pruefen()` stimmt die Beschreibung nicht. Da steht für alle ungeraden `h`, aber wenn `ausgabe` den Wert `True` hat, werden auch die geraden `h` geprüft. Was stimmt denn nun? Der Code oder die Dokumentation?

`ausgabe` hat Dennis89 ja bereits angesprochen. Im Code ist dann auch noch mal ein Kommentar der nötig ist, weil das so ein bisschend verwirrend ist. Kann es sein das in beiden Fällen nur die ungeraden `h` bearbeitet werden sollen? Dann steht da fast der gleiche Code doppelt. Der sollte da nur einmal stehen und nur das `print()` sollte bedingt ausgeführt werden.

Ja, da ist mir was durcheinander gegangen. Es sollten erstmal nur die ungeraden h geprüft werden, unabhängig von ausgabe.

Quote from __blackjack__

Es sieht auch so ein bisschen danach aus, als wenn hier keine selbst gebastelte Lösung, sondern vielleicht Logging verwendet werden sollte.

Was ist das?

Quote from __blackjack__

`exp` in `bosma()` wird auch berechnet wenn das gar nicht verwendet wird. Das sollte in den entsprechenden ``if``-Zweig wandern.

Verstehe ich.

Auch die anderen Punkte habe ich größtenteils nachvollzogen.

Eine Frage habe ich zu raise ValueError. Ich verstehe den Hintergrund, eine Funktion sollte immer die gleiche Art von Werten zurückgeben.

Wenn ich nun aber sagen wir 10.000 Koeffizienten h überprüfe, dann kann es sein, dass nur ein einziger nicht mit der vorgelegten Liste von Primzahlen überprüft werden kann, alle anderen aber schon. Wenn dann ein ValueError auftritt, wird die komplette Berechnung doch sofort abgebrochen. Ich würde aber gerne auch nachhalten, wieviele Kandidaten ich nicht berechnen konnte.

Danke sehr! Ich gehe nun weiter den Code durch.

__blackjack__ · May 28th 2023

neuernutzer Logging ist das Protokollieren von Ereignissen, unter anderem üblicherweise mit verschiedenen Leveln über die man steuern kann was ausgegeben wird und was nicht. Da gibt's Bibliotheken für. In der Python-Standardbibliothek das `logging`-Modul oder auch das externe `loguru`. Letzteres nehme ich ganz gerne für Programme und Dienste.

Was wäre denn so schlimm daran das ``if`` in der Schleife zu haben? Was passieren kann wenn man der Code stattdessen zweimal schreibt oder kopiert hast Du ja gesehen: Wenn man dabei Fehler macht, dann machen die beiden ”Kopien” nicht mehr das gleiche. Redundanter Code macht mehr Arbeit beim Schreiben und beim Verändern, und ist fehleranfälliger.

Ob bei einer Ausnahme die komplette Berechnung abgebrochen wird, hängt davon ab was der Code der das aufruft mit der Ausnahme macht. Bei -1 funktioniert der Indexzugriff nicht, da wäre ein `IndexError` die Folge, also nichts gewonnen. Und bei dem `None` ignoriert der Code Argumente mit denen `jacobisymbol()` gar nicht hätte aufgerufen werden dürfen, was im Grunde ja ein Programmierfehler ist, der dann einfach ignoriert wird. Und wenn man diese Fälle auch irgendwie aufzeichnen will, müsste man Code schreiben der explizit auf die -1 oder das `None` prüft. Wenn man das macht, könnte man an der Stelle genau so gut Code schreiben der statt einen besonderen, magischen Rückgabewert zu prüfen, die Ausnahme behandelt. Ausnahmen wurden erfunden damit man keine magischen Fehlerwerte braucht, und die immer überall prüfen muss, oder falls man das nicht tut, mit falschen Werten weitergearbeitet wird, als wäre nichts passiert.

neuernutzer · May 28th 2023

Quote from __blackjack__

neuernutzer

Was wäre denn so schlimm daran das ``if`` in der Schleife zu haben? Was passieren kann wenn man der Code stattdessen zweimal schreibt oder kopiert hast Du ja gesehen: Wenn man dabei Fehler macht, dann machen die beiden ”Kopien” nicht mehr das gleiche. Redundanter Code macht mehr Arbeit beim Schreiben und beim Verändern, und ist fehleranfälliger.

"Schlimm" ist es nicht. Aber bei sovielen Werten dauert es doch wirklich länger, jedesmal die if-Abfrage durchzuführen. Dabei sind 10.000 viel zu wenig, sagen wir mal 2^20 Werte.
Jedenfalls war das mein Gedanke dahinter.

__blackjack__ · May 28th 2023

Naja, das fällt IMHO in die Kategorie „premature optimization“. Wie gesagt sieht man ja was passieren kann wenn man wegen so etwas Code zweimal schreibt — der muss dann nicht das gleiche machen wenn man nicht aufpasst. Und wenn man hier irgendwie Probleme mit der Geschwindigkeit hat, dann ist nicht ein einzelnes ``if`` das wirkliche Problem, sondern dass das in Python geschrieben ist. Solche Mikrooptimierungen sollte man nur machen wenn man da auch Messungen macht ob und was das am Ende tatsächlich bringt. Insbesondere wenn dadurch der Code fehleranfälliger wird.

neuernutzer · May 28th 2023

Quote from __blackjack__

Naja, das fällt IMHO in die Kategorie „premature optimization“. Wie gesagt sieht man ja was passieren kann wenn man wegen so etwas Code zweimal schreibt — der muss dann nicht das gleiche machen wenn man nicht aufpasst. Und wenn man hier irgendwie Probleme mit der Geschwindigkeit hat, dann ist nicht ein einzelnes ``if`` das wirkliche Problem, sondern dass das in Python geschrieben ist. Solche Mikrooptimierungen sollte man nur machen wenn man da auch Messungen macht ob und was das am Ende tatsächlich bringt. Insbesondere wenn dadurch der Code fehleranfälliger wird.

Das sehe ich ein. Python ist für den Zweck der Statistik sicher nicht die richtige Wahl. Für meine Zwecke ist es allerdings ausreichend. Trortdem möchte ich python verstehen und das beste rausholen.

Als nächstes stellt sich mir die Frage, ob ich das nun paralleliseren kann.

Das plotten kann erstmal in den Hintergrund rücken. Ich habe nun über den Nachmittag die Exponenten bis einschließlich 26 überprüfen können und die Ergebnisse pro Exponent in einer Datei abgelegt. Diese kann ich nun laden und plotten.

Jetzt frage ich mich, ob ich nicht bei meinen 4 Kernen jeweils 4 Exponenten gleichzeitig überprüfen kann. Das war der Hintergrund für den ersten Thread.

Lohnt es sich, das hier zu machen? Wie erkenne ich das?
Mit "Lohnen" meine ich nicht die allgemeine Fragestellung, ob python oder mein Rechner überhaupt schnell genug ist dafür, sondern prinzipiell.

__blackjack__ · May 29th 2023

neuernutzer Habe meinen Code jetzt tatsächlich mal ausprobiert. So viele Primzahlen mussten ja gar nicht generiert werden um das mal zu testen.

Und es läuft in die Ausnahme bei der Jacobisymbol-Funktion, weil die erste Primzahl, die 2, die Bedingung nicht erfüllt nicht gerade zu sein:

Python

...
  File "forum4.py", line 69, in bosma
    jacobi = jacobisymbol((h % prime) * 2 ** (k % (prime - 1)) + 1, prime)
             │             │   │              │    │                └ 2
             │             │   │              │    └ 2
             │             │   │              └ 2
             │             │   └ 2
             │             └ 1
             └ <function jacobisymbol at 0x7fa6137ae940>

  File "forum4.py", line 39, in jacobisymbol
    raise ValueError(f"expexted odd n ≠ 1, got {n!r}")

ValueError: expected odd n ≠ 1, got 2

Display More

Nach dem laden der Primzahlen könnte man die 2 einfach entfernen, dann läuft man nicht mehr in dieses Problem.

Wenn man sich anschaut wie man das am besten parallelisieren kann, fällt auf, das die Aufteilung des Codes auf Funktionen ungünstig ist. Eine Funktion sollte nur eine in sich geschlossene Sache machen. `plotten()` sollte Ergebnisse plotten, und nicht vorher auch noch die Berechnung machen. Ich habe das mal innerhalb der Funktion getrennt, das erst die Ergebnisse berechnet, und dann geplottet wird. Dann kann man statt der Schleife oder „list comprehension“ welche die Ergebnisse berechnet, `map()` beziehungsweise `starmap()` von `multiprocessing.Pool` verwenden und mehrere Exponenten parallel verarbeiten lassen.

Python

#!/usr/bin/env python3
import pickle
from itertools import accumulate, repeat
from multiprocessing import Pool
from pathlib import Path

import matplotlib.pyplot as plt

SELF_PATH = Path(__file__).parent
PRIMES_FILE_PATH = SELF_PATH / "primes_up_to_10_7"


def save(object_, file_path):
    with open(file_path, "wb") as file:
        pickle.dump(object_, file)


def load(file_path):
    with open(file_path, "rb") as file:
        return pickle.load(file)


def jacobisymbol(a, n):
    if n == 1 or n % 2 == 0:
        raise ValueError(f"expected odd n ≠ 1, got {n!r}")

    if a == 2:
        return 1 if n % 8 in [1, 7] else -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)

    return pow(-1, ((n - 1) // 2) * (a - 1) // 2) * jacobisymbol(n, a)


def bosma(h, k, primes):
    # print(f"Bosmatest für {h}*2^{k}+1")
    for prime in primes:
        jacobi = jacobisymbol((h % prime) * 2 ** (k % (prime - 1)) + 1, prime)
        if jacobi == 0:
            return False

        if jacobi == -1:
            exp = h * 2 ** (k - 1)
            n = 2 * exp + 1
            return pow(prime, exp, n) == n - 1

    raise ValueError(
        f"nicht genug Primzahlen für einen Test von {h=} und {k=}"
    )


def intervall_pruefen(untergrenze, exponent, primzahlen, ausgabe=False):
    """
    Diese Funktion prüft die Zahlen der Form :math:`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 einer
    Liste zurückgegeben.
    """
    print("Aktueller Exponent:", exponent)

    if untergrenze % 2 == 0:
        untergrenze += 1
        print(f"Untergrenze auf {untergrenze} angepasst")

    obergrenze = 2**exponent - 1
    anzahl = (obergrenze - untergrenze) // 2
    ergebnis = []
    for h in range(untergrenze, obergrenze + 1, 2):
        if ausgabe:
            print((h - 1) // 2 / anzahl, end="\r")
        ergebnis.append((h, exponent, bosma(h, exponent, primzahlen)))

    return ergebnis


def plotten(untergrenze, exponenten, primzahlen):
    # ergebnisse = [
    #     intervall_pruefen(untergrenze, exponent, primzahlen, False)
    #     for exponent in exponenten
    # ]
    with Pool() as pool:
        ergebnisse = pool.starmap(
            intervall_pruefen,
            zip(repeat(untergrenze), exponenten, repeat(primzahlen)),
        )

    for ergebnis in ergebnisse:
        plt.plot(
            [h for h, _, _ in ergebnis],
            list(
                accumulate(
                    int(bosma_result) for _, _, bosma_result in ergebnis
                )
            ),
        )
    plt.show()


def main():
    primzahlen = load(PRIMES_FILE_PATH)
    assert primzahlen[0] == 2
    #
    # `jacobisymbol()` erwartet ungerade `n`, darum kann/muss man die gerade
    # Primzahl 2 hier entfernen.
    #
    del primzahlen[0]
    
    plotten(1, range(2, 25), primzahlen)


if __name__ == "__main__":
    main()

Display More

neuernutzer · May 29th 2023

Quote from __blackjack__

neuernutzer Habe meinen Code jetzt tatsächlich mal ausprobiert. So viele Primzahlen mussten ja gar nicht generiert werden um das mal zu testen.

Das ist richtig, statistisch braucht man nur sehr wenige und auch nur sehr kleine. Leider ist das noch nicht bewiesen.

Zu deiner weiteren Anmerkung: Mit map habe ich mich noch nicht beschäftigt, das schaue ich mir aber an.

Ich hatte gestern abend den Code auch nochmal leicht geändert, sodass "plotten" jetzt wirklich nur den Vektor plottet, den es bekommt. Hier erstmal der ganze Code, Änderungen sind ab Zeile 113:

Python

import time
import matplotlib.pyplot as plt
import pickle
from pathlib import Path

PRIMES_FILE = Path("/home/neuernutzer/primes_up_to_10_7")
H_UNTER = 1
H_OBER = 0
EXPONENTEN_UNTER = 27
EXPONENTEN_OBER = 27
STEP = 1

def save_to_file(liste, datei):
    with open(datei, 'wb') as file:
        pickle.dump(liste, file)
        
def load_from_file(datei):
    with open(datei, 'rb') as file:
        return pickle.load(file)
    
def jacobisymbol(a, n):
    if n == 1 or n % 2 == 0:
        return
    if a == 2:
        return 1 if n % 8 in [1,7] else -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)

def bosma(h, k, primes):
    for D in primes:
        jacobi = jacobisymbol( (h % D)*2**(k % (D-1)) + 1, D)
        if jacobi == 0:
            return [False, D]
        if jacobi == -1:
            exp = h*2**(k-1)
            n = 2*exp+1
            return [True, D] if pow(D, exp, n) == n-1 else [False, D]
    return -1

# 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, primzahlen, ausgabe=False):
    if untergrenze % 2 == 0:
        untergrenze += 1
        print(f"Untergrenze auf {untergrenze} angepasst")
    if obergrenze == 0 or obergrenze >= 2**exponent:
        obergrenze = 2**exponent - 1
        print(f"Obergrenze auf {obergrenze} angepasst")
    ergebnis = []
    anzahl = (obergrenze - untergrenze)//2 + 1
    if ausgabe:
        for h in range(untergrenze, obergrenze+1, 2):
            print("Stand:", (h+1)//2 / anzahl, end='\r')
            B = bosma(h, exponent, primzahlen)
            ergebnis.append((h, exponent, B[0]))
        return ergebnis
    else:
        for h in range(untergrenze, obergrenze+1, 2):
            B = bosma(h, exponent, primzahlen)
            ergebnis.append((h, exponent, B[0]))
        return ergebnis

def vektor_erstellen(untergrenze, obergrenze, exponent, primzahlen):
    vektor = intervall_pruefen(untergrenze, obergrenze, exponent, primzahlen, False)
    plotvektor = liste_vorbereiten(vektor)
    return plotvektor

def plotten(plotvektor):
    plt.plot(plotvektor[0], plotvektor[1])
    plt.show()
    
    
def main():
    primzahlen = load_from_file(PRIMES_FILE)
    print("Primzahlen geladen")
    exponenten = range(EXPONENTEN_UNTER, EXPONENTEN_OBER+1, STEP)
    print("Exponenten geladen")
    for exponent in exponenten:
        start = time.time()
        vektor = vektor_erstellen(H_UNTER, 2**exponent-1, exponent, primzahlen)
        SAVEPATH = Path("/home/neuernutzer/Auswertungen/"+str(exponent))
        save_to_file(vektor, SAVEPATH)
        print(f"Exponent {exponent} gespeichert in {SAVEPATH}")
        print(f"Das dauerte {time.time() - start} Sekunden")
    
if __name__ == "__main__":
    main()

# SAVEPATH = Path("/home/neuernutzer/Auswertungen/"+"27")
# vektor = load_from_file(SAVEPATH)
# plotten(vektor)

Display More

Kommentiere ich die letzten drei Zeilen ein und die Zeilen 118 bis 134 aus, kann ich die einzelnen Dateien plotten lassen. So hatte ich mir das gewünscht.

Möglicherweise ist meine Vorstellung von Parallelisierung falsch, aber kann ich nun nicht sagen: Kern 1 prüft die Exponenten 1,5,9..., Kern 2 prüft 2,6,10... und so weiter?

__blackjack__ · May 29th 2023

neuernutzer Diese Aufteilung willst Du nicht selbst machen. Da kümmert sich der `Pool` drum, dass die Prozesse ausgelastet sind.

neuernutzer · May 29th 2023

Quote from __blackjack__

neuernutzer Diese Aufteilung willst Du nicht selbst machen. Da kümmert sich der `Pool` drum, dass die Prozesse ausgelastet sind.

Echt schade. Ich hatte wirklich gehofft auf
"Prüfe wieviele Kerne vorhanden sind, dann verteile entsprechend dieser".

Und dieser "Pool", ist das der, den du auch in deinem vorherigen Beitrag bereits im Code angebracht hast?

Edit: Ah, jetzt habe ich den Code bei mir mal laufen lassen. Das ist schon auf die Kerne verteilt.
Hab ich mich leider wieder alles viel zu einfach vorgestellt.

Wenn ich das richtig sehe (anhand des System Monitors), dann nimmt sich jetzt einfach ein Kern, der gerade frei wird, den nächsten offenen Exponenten.

Mein Gedanke war, dass ganze nun auch mal auf unsere Uni-internen Hochleistungsrechner (sagt man das so? :D) laufen zu lassen.

Ihr werdet jetzt sagen, dass python dafür die ganz falsche Wahl ist, das verstehe ich. In meiner Arbeit wird nunmal python benutzt und ich werde jetzt nicht extra dafür noch C++ lernen. Das ist aber auch nicht das wichtige, es sollen einfach Berechnungen und Statistiken angestellt werden, die über das handeslübliche Maß hinausgehen. Die Exponenten bis 50 zu untersuchen, ist mit einem Haushaltsrechner halt nicht möglich. Und auch 1000 verschiedene Werte h für den Exponenten k = 50.000 wird nicht mehr möglich sein.

neuernutzer · May 29th 2023

Noch was: Wenn ich mir einen festen Exponenten nehme, beispielsweise k = 50.000, und ich möchte die Koeffizienten h von 1001 bis 2001 überprüfen. Dann kann ich doch auch das parallelisieren? Wie würde das aussehen?

Bitte Kritik zu meinem Code

Similar Threads

Treiberproblem für LCD Anzeigen

Users Viewing This Thread