C78E 3718 1375 E408 A2C3
CCEA 7D67 D672 420A 14E6

last update:

Da das vorhin in der Übung etwas durcheinander gegangen ist, habe ich mal die durch $y^2=x^3+3x+5$ gegebene elliptische Kurve $E$ mit den Punkten $P=(-1,1)$ und $Q=(4,9)$ geplottet, inklusive der Konstruktion von $P+Q=-\left(\frac{11}{25},\frac{237}{125}\right)$. Das geht wie folgt in Sage:

# x ist bereits eine Variable, aber y kennt Sage nicht.  Wir legen also als
# erstes eine neue Variable an.
y = var('y')

# Die ellipitsche Kurve E
E = implicit_plot(y^2 == x^3+3*x+5, (x,-2, 4.5), (y, -10, 10), color="black")
# und die beiden gegebenen Punkte
P = point((-1,1), color="red", size=20) + text("P", (-1.125,1.5), color="red")
Q = point((4,9), color="red", size=20) + text("Q", (3.875, 9.3), color="red")

# Gerade durch P und Q
g = line([(-3.5, -3), (4+5/8,10)], color="green") + text("g", (1.5, 5.7), color="green")
# sowie der weitere Schnittpunkt von g und E
R = point((-11/25, 237/125), color="blue", size=20) + text("R", (-0.3, 1.5))

PQ = point((-11/25, -237/125), color="blue", size=20) + text("P+Q", (-0.2, -1.5))
s = line([(-11/25, 237/125), (-11/25, -237/125)], color="grey", linestyle=":")
Plot = E+P+Q+g+R+PQ+s
# Ordentliches Format wählen
Plot.set_aspect_ratio(0.25)
# und anzeigen
Plot

Damit erhält man Diagramm zur Addition auf E

Falls ihr das mal braucht: Solche Plots kann man mit (in diesem Beispiel) Plot.save("dateiname.endung") abspeichern, wobei man als Endung unter anderem png, pdf und svg benutzen kann und Sage dann Dateien in dem entsprechenden Format erstellt.

Aufgabe 1

Der LLT ist schnell implementiert:

def is_prime_lucas_lehmer(p):
    """Check whether 2^p-1 is a Mersenne prime using the Lucas-Lehmer test."""
    s = Mod(4, 2^p - 1)
    for i in range(3, p+1): s = s^2 - 2
    return s == 0

Um eine Übersicht über die Laufzeite zu erhalten, lassen wir den Test mit Primzahlen verschiedener größe laufen. Zum Stoppen der Laufzeit benutzen wir die time()-Funktion.

from time import time
print "| $p$  | Prime (Yes/No) | Laufzeit des LLT (s) |"
print "|------|----------------|----------------------|"
for n in [1..10] + [100, 500, 1000, 1500, 2000, 3000, 4000, 5000]:
    p = nth_prime(n)
    print "| {} |".format(str(p).rjust(5, " ")),
    start = time()
    b = is_prime_lucas_lehmer(p)
    duration = time() - start
    if b:
        print "      Yes     ",
    else:
        print "      No      ",
    print "| {} |".format(str(duration).ljust(len("Laufzeit des LLT (s)"), " "))

Als Ergebnis erhalten wir die folgende Tabelle (die zweite Spalte $M_p$ habe ich mal weggelassen, da die Zahl doch recht groß werden…):

$p$ Prime (Yes/No) Laufzeit des LLT (s)
2 No 0.00597095489502
3 Yes 0.000696897506714
5 Yes 0.000622987747192
7 Yes 0.000731945037842
11 No 0.000629901885986
13 Yes 0.000751972198486
17 Yes 0.00091290473938
19 Yes 0.000904083251953
23 No 0.000869989395142
29 No 0.000170946121216
541 No 0.00384306907654
3571 No 0.133982896805
7919 No 0.841007947922
12553 No 1.47259402275
17389 No 3.35499310493
27449 No 10.6544179916
37813 No 23.6767170429
48611 No 42.8488218784

Aufgabe 2

def sum_of_proper_divisors(n):
    if n == 0:
        return 0

    # Zahlen vom Typ `int` (welche Python standardmäßig für ganze Zahlen
    # verwendet) haben keine Methode `divisors()`.  Daher müssen wir den
    # Sage-Typ `Integer` verwenden.
    # Summiere über alle Teiler bis auf den letzten (n selbst).
    return sum(Integer(n).divisors()[:-1])

for n in [1..1000]:
    f = sum_of_proper_divisors
    lst = [n]
    # Berechne nur die ersten 50 Folgenglieder
    for k in [1..50]:
        # Wenn die Folge ab jetzt nur noch Nullen enthält, oder eine
        # vollkommene Zahl, die sich dementsprechend unendlich oft wiederholt,
        # bringt es nichts, diesen Teil der Folge auszugeben.  Daher beenden
        # wir die Schleife in diesem Fall vorzeitig.
        if lst[-1] == 0 or len(lst) >= 2 and lst[-1] == lst[-2]:
            break

        # Ansonsten fügen wir das nächste Folgenglied hinzu.
        lst.append(f(lst[-1]))

    # Falls die Folge periodisch mit einer Periode der Länge 1 ist, dann geben
    # wir das an.
    if lst[-1] == lst[-2]:
        print lst, "und weiter mit Periode", lst[-1]
    else:
        print lst

Wenn eine Folge eine Periode größerer Länge hat, kann man dies natürlich auch leicht feststellen. Man überprüft einfach ob f(lst[-1]) schon in lst auftaucht.

Für diejenigen unter euch, die gerne ein bisschen tüfteln, gibt bei Project Euler unter anderem eine gut zu diesem Thema passende Aufgabe.

Etwas verspätet ist hier meine Lösung zu Aufgabe 3:

# Laut Skript ist $n$ genau dann eine Kongruenzzahl, wenn die elliptische Kurve
# $y^2=x(x^2-n^2)$ eine rationale Lösung $(x,y)$ mit $y≠0$ besitzt.
for n in [37,41,53]:
    # Erzeuge die richtige elliptische Kurve y^2=x(x^2-n^2)
    E = EllipticCurve(QQ, [0,0,0,-n^2,0])
    # und wähle einen (möglichst allgemeinen) Punkt auf dieser.
    (x,y,z) = E.an_element()
    # Überprüfen, ob der Punkt wirklich auf der Kurve liegt und nicht der
    # unendlich ferne Punkt ist.
    assert(z == 1 and y^2 == x*(x^2-n^2))

    # Berechne die Seitenlängen des Dreiecks
    a = abs(y / x)
    c = abs(2*x/a - a)
    b = sqrt(c^2 - a^2)
    # und seine Fläche
    A = a*b/2
    # Nochmal überprüfen, dass wirklich nichts schiefgegangen ist
    assert(A == n)

    # und das Ergebnis schön formattiert ausgeben lassen.
    print """Die elliptische Kurve ${}$ enhält den nicht-trivialen rationalen
Punkt ${}$. Also hat das Dreick mit (rationalen) Seiten
$a={}$, $b ={}$ und $c={}$
den Flächeninhalt $\\frac{{ab}}2={}$.
""".format(latex(E), latex((x,y)), latex(a), latex(b), latex(c), A)

was die folgende Ausgabe liefert:

Die elliptische Kurve $y^2=x^{3}-1369x$ enhält den nicht-trivialen rationalen Punkt $\left(-\frac{1764}{21025}, -\frac{32672766}{3048625}\right)$. Also hat das Dreick mit (rationalen) Seiten $a=\frac{777923}{6090}$, $b =\frac{450660}{777923}$ und $c=\frac{605170417321}{4737551070}$ den Flächeninhalt $\frac{ab}2=37$.

Die elliptische Kurve $y^2=x^{3}-1681x$ enhält den nicht-trivialen rationalen Punkt $\left(-9, 120\right)$. Also hat das Dreick mit (rationalen) Seiten $a=\frac{40}{3}$, $b =\frac{123}{20}$ und $c=\frac{881}{60}$ den Flächeninhalt $\frac{ab}2=41$.

Die elliptische Kurve $y^2=x^{3}-2809x$ enhält den nicht-trivialen rationalen Punkt $\left(-\frac{1158313156}{35343025}, -\frac{50101876246422}{210114283625}\right)$. Also hat das Dreick mit (rationalen) Seiten $a=\frac{1472112483}{202332130}$, $b =\frac{21447205780}{1472112483}$ und $c=\frac{4850493897329785961}{297855654284978790}$ den Flächeninhalt $\frac{ab}2=53$.

In Anbetracht dieser Seitenlängen verwundert es nicht, dass die Suche beim letzten Blatt zwar $41$, aber weder $37$, noch $53$ als Kongruenzzahl gefunden hat.

Aufgabe 1

def find_solution(a,b,c, N=10, all=False):
    """Löse ax^2+by^2+cz^2=0 mit 1≤x,y,z≤N, nicht alle 0, wobei standardmäßig
    N=10 ist."""
    results = []

    # Iteriere einfach über alle möglichen (x,y,z) im zu relevanten Bereich.
    for x in range(1, N):
        for y in range(1, N):
            for z in range(1, N):
                if x == y == z == 0:
                    continue
                if a*x^2 + b*y^2 + c*z^2 == 0:
                    if all:
                        results.append((x,y,z))
                    else:
                        return (x,y,z)

    return results

# Suche nach allen Lösungen mit Koordinaten ≤ 10
for lst in [(3, -10, 13), (5, 7, -13), (-7, 11, 13)]:
    print find_solution(*lst, all=True)

Das liefert die Ausgabe
[(1, 8, 7), (3, 2, 1), (6, 4, 2), (7, 4, 1), (9, 6, 3)]
[(1, 4, 3), (2, 8, 6), (3, 1, 2), (5, 7, 6), (6, 2, 4), (9, 3, 6)]
[(3, 1, 2), (4, 3, 1), (6, 2, 4), (8, 6, 2), (9, 3, 6)]

Aufgabe 2

Alternativ zur Funktion aus dem Skript geht es auch auch etwas kürzer:

def fundamental_solution(n):
    "Finde die Fundamentallösung von x^2-nx^2=1, wobei n kein Quadrat ist."
    # Definiere Polynomring über ℚ
    R.<x> = QQ[]
    # und erweitere ℚ um √n
    K.<a> = NumberField(x^2-n)
    # Berechne die multiplikative Gruppe von K
    G = K.unit_group()

    # G.1 ist der Erzeuger der zu ℤ isomorphen Untergruppe von K^*
    return list(K(G.1))

x, y = fundamental_solution(571)
assert(x^2-571*y^2 == 1)
(x,y)

Mit den Mitteln aus der Vorlesung kann man das Problem aber mit dem Programm aus dem Skript lösen:

def Pell_fs(n) :
    """Returns the fundamental solution of xˆ2−ny ˆ2=1."""
    U = lambda a : matrix (ZZ, 2, [a, 1, 1, 0])
    phi = lambda x : 1 / (x - floor(x))
    w = QQbar(sqrt(n))
    s = phi(w); lst = [w, s]
    k = 2; t = phi(s)
    while is_even(k) or t != s:
        lst.append(t)
        k += 1; t = phi(t)
    cf = map(lambda t : floor(t) , lst)
    M = prod(U(a) for a in cf)*U(cf[0])^-1
    return M[0,0], M[1,0], len(cf)-1

(x,y,n) = Pell_fs(571); (x,y,n)

Das Ergebnis kann man schön ausgeben lassen mittels

print 'x = {:27,}\ny = {:27,}'.format(int(x), int(y)).replace(",", " ")

Das liefert als Ausgabe x und y als 27-stellige Zahl mit führenden Leerzeichen und Zifferngruppierung.

Aufgabe 3

Finde mithilfe der Parametrisierung pythagoräische Tripel und sammle die dabei gefundenen Kongruenzzahlen. Dabei werden zunächst quadratfreie Kongruenzzahlen gesucht und diese anschließend mit allen in frage kommenden Quadraten multipliziert werden.

def find_congruent_numbers(N=100, k=100):
    # Menge der gefundenen Kongruenzzahlen
    congruent_numbers = set()

    for p in range(1, N):
        for q in range(1, p):
            if gcd(p, q) != 1:
                continue

            # a = (p^2-q^2), b = 2*p*q ⇝ n = ab/2 = (p^2-q^2)pq
            n = (p^2-q^2)*p*q

            if n.is_integral():
                n = n.squarefree_part()
                if n > k:
                    continue
                congruent_numbers.add(n)

    for c in list(congruent_numbers):
        for square in [x^2 for x in ([2..9] + [2*3, 2*4, 3*3, 2*5])]:
            if c*square <= k:
                congruent_numbers.add(c*square)

    return congruent_numbers

# Berechne nun für 0<p,q<1000 die auftretenden Kongruenzzahlen.
S = find_congruent_numbers(1000); print S
print len(S)

So findet man 33 der insgesamt 50 Kongruenzzahlen ≤100: 5, 6, 7, 13, 14, 15, 20, 21, 22, 24, 28, 30, 34, 39, 41, 45, 46, 52, 54, 55, 56, 60, 63, 65, 69, 70, 78, 80, 84, 85, 86, 88, 96.

Aufgabe 1

Standardmäßig gibt es in Sage nur eine symbolische Variable. Wir legen also zunächst zwei neue Variablen x1 und x2 an (die von Sage als x_1 bzw. x_2 ausgegeben und per LaTeX in $x_1$ bzw. $x_2$ übersetzt wird, wenn man die Option Typeset bei Sage aktiviert).

# Setze x = (x1, x2)
x = x1, x2 = var("x_1 x_2")
# und entsprechend u = (u1, u2).
u = u1, u2 = (sin(pi/2*x1)/cos(pi/2*x2), sin(pi/2*x2)/cos(pi/2*x1))
# Die Jacobideterminante ist also
D = jacobian(u, x).det(); D

was Sage als \[ \frac{1}{4} \, \pi^{2} - \frac{\pi^{2} \sin\left(\frac{1}{2} \, \pi x_{1}\right)^{2} \sin\left(\frac{1}{2} \, \pi x_{2}\right)^{2}}{4 \, \cos\left(\frac{1}{2} \, \pi x_{1}\right)^{2} \cos\left(\frac{1}{2} \, \pi x_{2}\right)^{2}} \] ausgibt. Andererseits ist \[ 1-u_1^2u_2^2=1-\frac{\sin\left(\frac12\,πx_1\right)^2\sin\left(\frac12\,πx_2\right)^2}{\cos\left(\frac12\,πx_1\right)^2\cos\left(\frac12\,πx_2\right)^2} \] Also ist

(D / (1-u1^2*u2^2)).full_simplify()

das gesuchte Ergebnis $\frac14\pi^2$.

Aufgabe 4

Zur Veranschaulichung und zum Testen wählen wir zunächst einen primitiven Vektor

b = vector([132, 162, 43, 72, 25, 341])
# Sicherstellen, das der Vektor primitiv ist
assert(gcd(b) == 1)

Nun der Code, der zu einem solchen Vektor eine in $\operatorname{GL}(n,ℤ)$ invertierbare Matrix konstruiert:

def unimodular_matrix_with_column(column):
    # Erweitere die Einheitsmatrix über ℤ um die gegebene Spalte.
    #
    # Spalten kann man in Sage nicht zu einer Matrix hinzufügen, daher fügen
    # wir den gegebenen Vektor als Zeile ein und transponieren das Ergebnis.
    M = matrix.identity(ZZ, len(column)).insert_row(0, column).transpose()

    # Solange die erste Spalte nicht `n` nur einen von Null verschiedenen
    # Eintrag hat:
    while sum(M.column(0)[1:]) != 0:
        # Sortiere die Matrix nach der 1. Spalte (ändert höchstens das
        # Vorzeichen der Determinante)
        M = matrix(sorted(M.rows(), reverse=True))

        i = m = -1
        for j in range(len(M.rows())-1, -1, -1):
            r = M.row(j)

            # Finde die Zeile mit dem kleinsten positiven Eintrag
            if r[0] == 0:
                continue
            if i == -1:
                i = j
                m = r[0]

            # Reduziere nun den 1. Eintrag der darüber liegenden Zeilen modulo m
            if j != i and M.row(j)[0] >= m:
                M.add_multiple_of_row(j, i, -(M.row(j)[0]//m))
                # Dies ändert die Determinante nicht, M (ohne die 1. Spalte)
                # ist also in jedem Schritt in GL(n,ℤ).

    # Das Ergebnis ist die Inverse der konstruierten Matrix ohne die 1. Spalte.
    return M.matrix_from_columns([1..len(column)]).inverse()

Testen wir das ganze:

A = gcd_matrix(b)
assert(abs(det(A)) == 1)
assert(A.column(0) == b)
A

Damit erhalten wir in diesem Fall \[ A = \begin{pmatrix} 132 & 0 & 5 & 0 & 58 \\
43 & 1 & 1 & 0 & 19 \\
75 & 0 & 3 & 1 & 33 \\
25 & 0 & 1 & 0 & 11 \\
341 & 0 & 13 & 0 & 150 \end{pmatrix} \]

Aufgabe 5

Zunächst einmal ist klar, dass die gesuchte Zahl der Koeffizient von $x^{10000·100}$ in der Potenzreihenentwicklung von \[ \frac1{(1-x^1)(1-x^2)(1-x^5)(1-x^{10})(1-x^{20})(1-x^{50})(1-x^{100})(1-x^{200})} \] um $0$ ist. Diesen berechnet man geschickterweise über das (Cauchy-)Produkt der Reihen \[ \sum_{i=0}^∞ x^i,\quad \sum_{i=0}^∞ x^{2i},\quad \sum_{i=0}^∞ x^{5i},\quad …,\quad \sum_{i=0}^∞ x^{200i}. \] Das ist zu viel Aufwand, um es per Hand zu machen. Per Computer kann man das aber leicht erledigen:

denominations = [1, 2, 5, 10, 20, 50, 100, 200]

def count_combinations(value, denom=denominations):
    """Berechne die Anzahl der Möglichkeiten, `value` als
    ℤ_{≥0}-Linearkombination der Werte aus `denom` zu schreiben."""

    c = [0]*(value+1)
    c[0] = 1

    for k in denom:
        for i in [0..value-k]:
            c[k+i] += c[i]

    return c[value]

Das gesuchte Ergebnis $99\,341\,140\,660\,285\,639\,188\,927\,260\,001$ erhalten wir nun durch

count_combinations(10000*100)

Das Problem, eine Zahl als $ℤ_{≥0}$-Linearkombination von gegebenen Zahlen mit möglichst kleiner Koeffizientensumme zu schreiben heißt auch change-making problem.

Hier die Musterlösung:

def find_squares( p):
    """
    Return the solution (x,y) of 0<x<=y of p=x^2+y^2 if it exists,
    otherwise throw an exception. Input must be prime p = +1 mod 4.
    """
    assert p%4 == 1, 'Error: %d must be a prime = 1 mod 4' % p
    # Find a solution of x^2=-1 mod p
    # using Wilson's theorem
    x = 1
    for j in range(2,(p-1)/2):
        x *= j
        x = x%p

    # Modify x if necessary so that p > x > p/2
    if 2*x < p:
        x = p - x

    # compute r_l as in the preceding theorem
    a=p; b=x
    while b*b > p:
        r=a%b; a=b; b=r
          
    a = int( math.sqrt(p-b*b))
    return (a,b) if a < b else (b,a)

can = [p for p in primes(10^6, 10^6+100) if p%4 == 1]
can
# Lösungen berechnen
sol = [ (p,find_squares( p)) for p in can]

# Lösungen ausgeben und überprüfen
for p,(a,b) in sol:
    print p, a, b, p == a^2+b^2

Aufgabe 2

Zunächst einmal lösen wir die Kongruenz modulo $7$. Da können wir sofort zwei Lösungen erraten: $3$ und $4$. Diese können wir nun sukzessive mithilfe des Newton-Verfahrens liften. Wir schreiben $i_n(a)$ für einen (beliebigen) Repräsentanten von $a^{-1}\bmod 7^n$.

  • Sei $x_1 = 3$.
  • $x_1 - (x_1^2-2)·i_1(2x_1) ≡ 3 - (9-2)·i_1(6) ≡ 3-7·(-1) ≡ 10\bmod 7^2$, setze also $x_1:=10$.
  • $x_2 - (x_2^2-2)·i_2(2x_2) ≡ … ≡ 108\bmod 7^3$, setze $x_3:=108$
  • $x_3 - (x_3^2-2)·i_3(2x_3) ≡ … ≡ 2166\bmod 7^3$, setze $x_4:=2166$

und so weiter. Entsprechend könnten wir natürlich auch mit $x_1=4$ anfangen. Auf diese Weise könnten wir die andere Lösung modulo $7^{100}$ bestimmen.

Die Berechnung ist zu aufwendig, um sie per Hand durchzuführen. Siehe auch Sage-Programm.

Eine kurze Übersicht über das Wurzelziehen modulo $m$. Allgemeine quadratische Kongruenzen kann man dann natürlich mithilfe der $p$-$q$-Formel lösen. Für das Wurzelziehen modulo $p$ mit einer Primzahl $p≡1\bmod 4$ kann man den Tonelli-Shanks-Algorithmus verwenden, aber den werdet ihr im Rahmen der Vorlesung „Elementare Zahlentheorie“ sicher nicht brauchen.

Bei den Tests aus dem SoSe 2008 könnt ihr euch ansehen, wie Klausuraufgaben in der Zahlentheorie aussehen könnten.

Aufgaben

  1. Bestimme die Lösungsmenge von $x≡1 \bmod 9,\quad x≡8\bmod 21,\quad x≡4\bmod 7$.
  2. Bestimme eine Lösung von $x^2+2x-4≡0\bmod 55$.
  3. Löse $x^2≡22\bmod 27$.
  4. Finde alle Primitivwurzeln modulo $17$. Berechne dazu geeignete Potenzen um zu überprüfen, ob ein Element Ordnung $φ(17)$ hat. Ist eine Primitivwurzel $g$ gefunden, kann man die restlichen leicht hinschreiben. (Für welche $x$ ist $g^x$ wieder eine Primitivwurzel?)

Beweise

Blatt 4, Aufgabe 1

$\newcommand{ord}{\operatorname{ord}}$ Sei $p≠2$ eine Primzahl und $w$ eine Primitivwurzel modulo $p^2$, also $\ord_{p^2} w=φ(p^2)=p(p-1)$. Zudem gebe es für jedes $α≥1$ eine zu $p$ teilerfremde Zahl $n_α$, so dass \[ w^{φ(p^{α-1})} = 1+n_α p^{α-1} \] gilt.

Der Induktionsanfang ist also bereits erledigt. Nehmen wir nun an, dass $w$ eine Primitivwurzel modulo $p^{α-1}$ ist, also $\ord_{p^{α-1}} w=φ(p^{α-1})=p^{α-2}(p-1)$. Das bedeutet, dass $w^{p^{α-2}(p-1)}≡1\bmod p^{α-1}$ (das gilt nach Euler immer) und $w^{p^{α-2}(p-1)/q}\not ≡1\bmod p^{α-1}$ ist für alle echten Teiler $q$ von $p^{α-2}(p-1)$.

Wir müssen nun zeigen, dass $n:=\ord_{p^α}(w)≥φ(p^α)$ ist, dass also $w^n\not ≡ 1\bmod p^α$ falls $n=\frac{φ(p^α)}q$ und $q$ ein echter Teiler von $φ(p^α)=p^{α-1}(p-1)$ ist. Da aus $w^n≡1\bmod p$ folgt, dass $p-1$ ein Teiler von $n$ ist (warum?), genügt es schon zu zeigen, dass \[ w^{p^{α-1}(p-1)/p}=w^{p^{α-2}(p-1)}\not ≡ 1\bmod p^α \] gilt. Das folgt aber sofort aus der Existenz von $n_α$, denn $1+n_αp^{α-1}\not ≡1\bmod p^α$, da $n_α$ teilerfremd zu $p$ ist.

Aufgabe 2

Im Beweis der Existenz von $n_α$ wird benutzt, dass $p\mid {p \choose 2}$ ist, was aber für $p=2$ falsch ist, denn ${2\choose 2} = 1$.