Jacobisymbol: effiziente Implementierung

Neue Frage »

Malcang Auf diesen Beitrag antworten »
Jacobisymbol: effiziente Implementierung
Hallo,

ich grüble aktuell an einer effizienten Implementierung des Jacobisymbols. Ich habe dazu zwar bereits viel gesucht, aber noch nicht das Richtige gefunden. Mir fällt das Thema leider auch immer wieder schwer, daher frage ich um euren Rat. Ich habe bisher folgendes gemacht (jeweils programmiert in phyton):

Eine Funktion prime_decomp(n), die mir für die Liste
code:
1:
 (p_1, alpha_1) ..., (p_n, alpha_n)

ausgibt. (Hier achte ich aktuell nicht auf Effizienz, daher habe ich den Code hier nicht geschrieben).

Eine Funktion isprime(n), welche mir zwei Werte zurückgibt. True bzw. False, falls die Zahl prim ist bzw. zusammengesetzt, und die Liste der Teiler aus prime_decomp(n).

So, dann habe ich eine Funktion zur Berechnung des Legendresymbols. Dies funktioniert nur wenn eine Primzahl ist. Berechnet wird dann das Euler-Kriterium:
code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
def legendre(a, p):
    if a % p == 0:   return 0
   
    if (( p == 2) or (isprime(p)[0] == False) ):
        print(p, "ist keine ungerade Primzahl.")
        return
    res = pow(a, (p-1) // 2, p)
    if res == p-1:   return -1
    return res


Diese Codes funktionieren auch.
Was mir Sorgen macht, ist der folgende Code:
code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
def jacobi(a, n):
    # Nur für ungerade n
    if n % 2 == 0:   return

    
    erg, teiler = isprime(n)

    # Falls n eine Primzahl ist, wird einfach das legendresymbol zurückgegeben.
    if erg == True:
        return legendre(a, n)
    
    #Nun werden die Ergänzungssätze des quadratischen Reziprozitätsgesetztes überprüft
    if a == 2:
        if n % 8 == 1 or n % 8 == n-1:   return 1
        return -1
    if a == -1:
        return pow(-1, (n-1) // 2)

    # Der Spezialfall a = 1 liefert 1 zurück
    if a == 1:   return 1

    # Der "Zähler" darf "modulo Nenner" berechnet werden
    if a > n:   a = a % n

    # Hier kommt nun die Anwendung des QRP
    if a % 4 == 1 or n % 4 == 1:   return jacobi(n, a)
    return -jacobi(n, a)


Das funktioniert aber noch nicht. Die Eingabe
code:
1:
jacobi( 57, 171)
sollte mir wegen ja 0 ausgeben. Aber er gibt nichts zurück. Außerdem fehlt mir ja noch die "Multiplikativität", aber ich sehe ehrlich gesagt nicht, wie ich das hier einbringe verwirrt

Kann mir jemand helfen?
IfindU Auf diesen Beitrag antworten »

Der Code ist in einer Endlosrekursion gefangen. Am Ende der Funktion ruft er die Funktion noch einmal auf: mit den gleichen Werten verdreht. Am Ende setzt du wegen und rufst dann wieder die Funktion auf, diesmal mit 0 als zweites Argument. Dafür ist Jacobi nicht definiert.

Ich würde ganz oben eine Prüfung einbauen, ob ist. Wenn ja, dann die Jacobi einfach 0. Ansonsten sieht das nach einer guten Erklärung aus: https://planetmath.org/calculatingthejacobisymbol.

Damit du besseres Gefühl hast, was deine Funktion macht: nutze einen Debugger oder gib dir aus, wo er ist: Ganz am Anfang print("Ich starte Jacobi-Berechnung mit den Werten a = .. und n = ... ").

Dann siehst du was er macht und erkennst schneller wenn er was seltsames macht.
Malcang Auf diesen Beitrag antworten »

Hallo IfindU,

vielen Dank für die Hinweise.
Ich habe den Code mit einigen Sätzen versehen:
code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
def jacobi(a, n):
    print("--------------------------------")
    print("Berechne J({},{})".format(a, n) )
    if n == 1 or n % 2 == 0:
        print("Nenner muss ungerade und ungleich 1 sein.")
        return
    if a == 2:
        print("Zähler ist 2")
        if n % 8 == 1 or n % 8 == n-1:
            print("1, denn Nenner = {} mod 8".format(n % 8))
            return 1
        print("-1, denn Nenner = {} mod 8".format(n % 8))
        return -1
    if a == -1:
        print("Zähler ist -1, also 1. Ergänzungssatz")
        return pow(-1, (n-1)//2)
    if a == 1:
        print("Zähler ist 1")
        return 1
    if a > n:
        print("Zähler ist größer Nenner, also umdrehen")
        return jacobi(a % n, n)
    faktor = pow(-1, ((n-1) // 2) * (a-1) // 2)
    print("Kein Spezialfall erreicht, also QRP mit Faktor", faktor)
    return faktor*jacobi(n, a)


Nun muss ich aber noch den Fall unterbringen, dass der Zähler eine gerade Zahl ist. Da bastle ich gerade noch dran.

Edit:
Ich habe nun folgenden Code, der auch zu funktionieren scheint. Jedenfalls habe ich es mit einigen Werten ausprobiert und mit denen von wolframalpha verglichen.
code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
def jacobi(a, n):
    print("--------------------------------")
    print("Berechne J({},{})".format(a, n) )
    if n == 1 or n % 2 == 0:
        print("Nenner muss ungerade und ungleich 1 sein.")
        return
    if a == 2:
        print("Zähler ist 2")
        if n % 8 == 1 or n % 8 == 7:
            print("1, denn Nenner = {} mod 8".format(n % 8))
            return 1
        else:
            print("-1, denn Nenner = {} mod 8".format(n % 8))
            return -1
    if a == 0:
        return 0
    if a % 2 == 0:
        print("Zähler ist gerade")
        k = 0
        while (a % 2 == 0):
            k += 1
            a = a // 2
        print("Berechne nun J(2, {})^{} * J({}, {})".format(n, k, a, n))
        return jacobi(2, n) ** k * jacobi(a, n)
    if a == -1:
        print("Zähler ist -1, also 1. Ergänzungssatz")
        return pow(-1, (n-1)//2)
    if a == 1:
        print("Zähler ist 1")
        return 1
    if a > n:
        print("Zähler ist größer Nenner, also umdrehen")
        return jacobi(a % n, n)
    faktor = pow(-1, ((n-1) // 2) * (a-1) // 2)
    print("Kein Spezialfall erreicht, also QRP mit Faktor", faktor)
    return faktor*jacobi(n, a)
IfindU Auf diesen Beitrag antworten »

Du kannst ja mal die Funktion hier nehmen: Link. Und einfach mal 1000 random-Werte für deine und die andere Funktion einsetzen und schauen, ob immer das gleiche rauskommt.

Edit: Oder einfach die Library nutzen: SymPy
Malcang Auf diesen Beitrag antworten »

Zitat:
Original von IfindU
Du kannst ja mal die Funktion hier nehmen: Link.

Die hatte ich schonmal ausprobiert, machte mir aber auch Probleme, wenn ich auch gerade nicht mehr weiß, was genau Big Laugh

Zitat:


Edit: Oder einfach die Library nutzen: SymPy


Das ich da nicht drauf gekommen bin geschockt
Die nehme ich auch ab jetzt.

Danke für die geduldige Hilfe!
IfindU Auf diesen Beitrag antworten »

Schön, dass es geklappt hat Freude Ich hoffe nur, du hast es implementiert, weil du es üben/verstehen wolltest und nicht, weil du es einfach gebraucht hast und einfach die Bibliothek mit allen diesen Funktionen übersehen hast smile
 
 
Neue Frage »
Antworten »



Verwandte Themen

Die Beliebtesten »
Die Größten »
Die Neuesten »