Elementare Zahlentheorie: Blatt 5

Aufgabe 1

Sage bringt bereits eine solche Funktion power, weshalb ich der Funktion einen anderen Namen gegeben habe.

def my_power(a, n):
    # Triviale Fälle
    if n == 0 or a == 1:
        return 1
    if a == 0:
        return 0

    result = 1
    while n > 0:
        # Wenn n ungerade ist:
        if n % 2 == 1:
            result *= a
        # Setze n auf die größte ganze Zahl ≤ n/2
        n = n//2
        # und quadriere a
        a *= a
    return result

Zum Testen:

for i in range(0,100):
    assert(2**i == my_power(2,i))
    assert(3**i == my_power(3,i))
    assert(371**i == my_power(371,i))

Aufgabe 2

Die folgende Funktion liftet bei jedem Durchlauf der Schleife eine Lösung modulo \(p^n\) mithilfe des Newton-Verfahrens zu einer Lösung modulo \(p^{n+1}\):

def lift_mod(x1, a, p, n):
    """Solve x^2 ≡ a mod p^n for some prime p and a given solution x0 of
    x^2 ≡ a mod p.
    """
    xi = x1
    for n in range(1,100+1):
        # Wir müssen Integer(...) schreiben, weil Sage sonst beim nächsten
        # Schleifendurchlauf xi als Restklasse modulo p^(n-1) interpretiert, wir
        # also nichts neues mehr erhalten würden.  Dann hätten wir einfach
        # xi=x1 für alle n.
        xi = Integer(xi - (xi^2-2)*Mod(2*xi, p^n)^-1) % (p^n)
    return xi

Man sieht leicht, dass \(3^2\equiv 2\bmod 7\) ist. Damit ist natürlich auch \(4\) eine Lösung modulo \(7\). Wir erhalten also ein solches \(x\) mittels

a = lift_mod(3, 2, 7, 100); b = lift_mod(4, 2, 7, 100)
print a
print b

# Lösungen überprüfen
assert(a^2 % 7^100 == 2)
assert(b^2 % 7^100 == 2)

Aufgabe 3

Die drei Primzahlen finden wir mittels

can = [p for p in range(10**6+1, 10**6+100, 4) if is_prime(p)]; can

Man beachte, dass in Sage sowohl x**y, als auch x^y verwendet werden kann, um \(x^y\) zu berechnen.

Nun müssen jeweils ein Darstellung als Summe zweier Quadrate finden. Dazu probieren wir einfach durch:

for p in can:
    for a in range(1, int(sqrt(p))):
        if 2*a*a > p:
            break
        if Integer(p - a*a).is_square():
            print "{} = {}² + {}²".format(p, a, int(sqrt(p-a*a)))