Hat man, wie im letzten Post beschrieben, Mehl oder eine Backmischung beschafft, kann man sich ans Backen begeben.

Grundsätzlich ist es nicht schwer ein Brot zu machen:

• 500g Mehl, 360g Wasser, 10g Salz und 2g Trockenhefe vermischen.
• Maschinell 10 Minuten oder manuell 20 Minuten durchkneten.
• Gehen lassen, bis sich das Volumen etwa verdoppelt hat (ca. 5 Stunden).
• Oberfläche mit Mehl bestäuben oder den Finger anfeuchten.
• Mit dem Finger eine ca. 1cm tiefe Delle in den Teig drücken. Wenn die Delle bleibt ist das Brot lang genug gegangen.
• Ab in den auf 240℃ vorgeheizten Ofen tun, Temperatur auf 220℃ runterdrehen.
• Bis zur gewünschten Bräune backen. Anschließend auf einem Rost abkühlen lassen.

Möchte man ein richtig gutes Brot backen, so kann man eine ganze Reihe weiterer Dinge beachten. Statt jetzt selbst zu beschrieben, was man beachten sollte, verweise ich hier lieber auf ein paar exzellente Quellen, die mir sehr geholfen haben.

Plötzblog und Brotbackbuch

Einen Haufen toller Rezepte und eine Unmenge von Hintergrundwissen gibt es im Plötzblog. Darüber habe ich meinen Einstieg gefunden. Der Betreiber, Lutz Geissler, hat zudem mehrere Bücher zum Thema geschrieben.

Das Brotbackbuch Nr. 1 habe ich selbst und kann es klar empfehlen. Auch hier gibt es eine Reihe toller Rezepte, die alle überaus ausführlich und auch für Anfänger leicht zu verstehen sind. Zudem werden für ein tieferes Verständnis nötigen Grundlagen erklärt und es gibt eine ganze Reihe von Tipps, z.B. woran es liegen kann, wenn ein Brot bröckelig geworden ist.

Sowohl beim Plötzblog, als auch beim Brotbackbuch stört mich jedoch, dass jedes Rezept eine andere Menge Brot ergibt.

Flour, Water, Salt, and Yeast

Deshalb habe ich in den letzten Monaten meist auf Flour, Water, Salt, and Yeast von Ken Forkish zurückgegriffen. Das Buch ist exzellent in dem was es beschreibt, deckt aber weniger ab als das Brotbackbuch. Die Rezepte sind standardisiert und basieren immer auf 1kg Mehl, werden alle idealerweise in einer gusseisernen Auflaufform bei 245℃ für 30 Minuten mit und 20 Minuten ohne Deckel gebacken.

Gerade wenn man gerne saftige Brote isst, sind die Rezepte aufgrund des meist recht langen Gehens und des hohen Wasseranteils ideal.

Neben den Rezepten, die wie im Brotbackbuch, jeweils alle wichtigen Details enthalten, wird am Anfang noch erklärt, wie man gute Brote ohne Kneten machen kann. Statt den Teig zu kneten wird er nämlich bei aller Rezepten in FWSY mehrfach gefaltet. Dazu braucht man keine Knetmaschine und es ist erheblich weniger aufwendig, als den Teig traditionell per Hand zu kneten. Meiner Erfahrung nach ist das Falten für ein Brot ein Arbeitsaufwand von unter 10 Minuten, also mit der Verwendung einer Knetmaschine vergleichbar, nur ist die Arbeit auf drei bis fünf kurze Abschnitte über ca. zwei Stunden verteilt.

Die Rezepte in FWSY wirken vielleicht etwas langweilig, weil sie alle nur kleine Variationen von zwei Grundlegenden Rezepten sind. Man kann sie aber leicht selbst variieren, indem man beispielsweise eingeweichte Getreidekörner, Sonnenblumenkerne, Leinsamen, Sesam, Mohn, Haferflocken oder was einem sonst so einfällt in den Teig oder auf den fertigen Laib tut.

Wenn man eine große Menge in den Teig tut sollte man allerdings darauf achten, dass der Teig nicht zu flüssig wird. Sonst kann es passieren, dass das Brot beim Schneiden immer am Brotmesser klebt. Das tut zwar dem Geschmack keinen Abbruch, aber es ist doch sehr unpraktisch alle paar Scheiben das Brotmesser reinigen zu müssen.

Seit ein paar Jahren backe ich häufig mein eigenes Brot. In diesem Rahmen habe ich einige Erfahrungen mit der Beschaffung von Mehl und anderen Zutaten gemacht, die ich hier einmal festhalten möchte.

Bestellen

Bestellt habe ich bisher bei mehlkaufen.de und der Blattert-Mühle.

mehlkaufen.de ist relativ teuer, hat dafür aber eine stattliche Auswahl von Backmischungen. Wenn man sich keine großen Gedanken machen möchte ist das eine schöne Sache. Wasser und Hefe dazu, vermischen, gehen lassen, backen. Fertig.

Von Backmischungen mit Körnern muss ich abraten. Wenn man die einfach so macht wie es drauf steht bekommt man ein Brot mit zum Teil sehr harten Körnern. Alternativ kann man die Körner gesondert einweichen (etwa über Nacht), dafür muss man die aber erst aussieben. Da ist es einfacher die Körner separat zu kaufen.

Bei der Blattert-Mühle gibt es Weizen- und Roggenmehl, gerade in größeren Packungen, sagen wir ab 5kg, ziemlich günstig mit (letzter mir bekannter Stand: knapp 0.90€/kg). Insbesondere kosten hier dunklere Mehle – also die mineralstoffreicheren, die mit einer größeren Typennummer – oder Vollkornmehle nicht mehr als die hellen.

Auch Dinkelmehl ist nicht übermäßig teuer. Dazu sollte ich aber sagen, dass ich bisher keinen geschmacklichen Unterschied zwischen Dinkel- und Weizenbrot feststellen konnte. Vielleicht fehlt mir da der direkte Vergleich, vielleicht müsste ich da dunklere Mehle probieren.

Im Laden kaufen

Angefangen habe ich mit Backmischungen von Aldi. Die gibt es nicht immer, sind aber sehr günstig. Da bezahlt man für 1kg Backmischung, also mehr als 1.5kg fertiges Brot weniger als anderswo 1kg Mehl kostet.

Bei Aldi gibt es auch teilweise relativ günstiges Mehl, das 405er ist aber was für Kuchen oder Teilchen, nicht für Brot/Brötchen/Pizza.

550er Weizenmehl gibt es aber zum Beispiel bei DM für 0.85€/kg. Backmischungen sind aber hier sündhaft teuer, teilweise um die 8€/kg.

Ansonsten kostet Mehl im Laden fast immer über 1.50€/kg. Bei kleinen Mengen ist das aber aufgrund der Versandkosten zum Teil günstiger oder zumindest nicht wesentlich teurer als im Versandhandel.

Packungsgrößen

Zum Brot backen braucht man relativ viel Mehl, verglichen mit dem, was sonst so haushaltsüblich ist. Ich habe in den letzten Jahren deutlich über 70kg Mehl verbraucht. Da lohnt es sich natürlich zu größeren Packungen zu greifen, nicht zuletzt, weil diese oft günstiger zu erwerben sind.

Dabei sollte man aber das eine oder anderen beachten. Zum einen kann es immer mal vorkommen, dass man zusammen mit dem Mehl auch Mehlmotten bekommt. Bei dem Mehl aus großen Mühlen, etwa bei dem aus dem Supermarkt, ist das seltener als bei kleinen Mühlen, die auch Endkunden beliefern. Zudem ist es ein Unterschied, ob die 1kg-Packung wegen Mottenbefall in den Müll wandert, oder eine 25kg-Packung.

Um das Mottenrisiko zu reduzieren bestellen manche Leute bevorzugt im Winter, wenn die Motten weniger aktiv sind und damit auch die Wahrscheinlichkeit sinkt, dass ein bestimmter Mehlsack verseucht ist. Zum Teil liest man auch, man solle das Mehl nach dem Kauf für ein paar Tage in den Gefrierschrank tun. Ob das wirklich hilft kann ich nicht sagen, die Motten müssen ja grundsätzlich immer irgendwie über den Winter kommen…

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 elliptische 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

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…):

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 formatiert ausgeben lassen.
    print """Die elliptische Kurve \\({}\\) enthält den nicht-trivialen rationalen
Punkt \\({}\\). Also hat das Dreieck 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\) enthält den nicht-trivialen rationalen Punkt \(\left(-\frac{1764}{21025}, -\frac{32672766}{3048625}\right)\). Also hat das Dreieck 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\) enthält den nicht-trivialen rationalen Punkt \(\left(-9, 120\right)\). Also hat das Dreieck 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\) enthält den nicht-trivialen rationalen Punkt \(\left(-\frac{1158313156}{35343025}, -\frac{50101876246422}{210114283625}\right)\). Also hat das Dreieck 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 \(\mathrm{GL}(n,\mathbb{Z})\) 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\cdot 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}^\infty x^i,\quad \sum_{i=0}^\infty x^{2i},\quad \sum_{i=0}^\infty x^{5i},\quad\ldots,\quad \sum_{i=0}^\infty 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 \(\mathbb{Z}_{\ge 0}\)-Linearkombination von gegebenen Zahlen mit möglichst kleiner Koeffizientensumme zu schreiben heißt auch change-making problem.

Lösungen

Aufgabe 2

low = 10**12
high = low+200
ps = []
for p in range(low+1, high, 2):
    if is_prime(p) and kronecker(5, p) == 1:
        ps.append(p)

print ps

Zusätzliche Übungen

Einige Aufgaben zum quadratischen Reziprozitätsgesetz

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)\cdot i_1(2x_1) \equiv 3 - (9-2)\cdot i_1(6) \equiv 3-7\cdot (-1) \equiv 10\bmod 7^2\), setze also \(x_1:=10\).
• \(x_2 - (x_2^2-2)\cdot i_2(2x_2) \equiv\ldots\equiv 108\bmod 7^3\), setze \(x_3:=108\)
• \(x_3 - (x_3^2-2)\cdot i_3(2x_3) \equiv\ldots\equiv 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]({{< ref "post/sose16-zt-05.md#aufgabe-2" >}}).

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\equiv 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\equiv 1 \bmod 9,\quad x\equiv 8\bmod 21,\quad x\equiv 4\bmod 7\).
2. Bestimme eine Lösung von \(x^2+2x-4\equiv 0\bmod 55\).
3. Löse \(x^2\equiv 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

Sei \(p\ne 2\) eine Primzahl und \(w\) eine Primitivwurzel modulo \(p^2\), also \(\mathrm{ord}_{p^2} w=φ(p^2)=p(p-1)\). Zudem gebe es für jedes \(α\ge 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 \(\mathrm{ord}_{p^{α-1}} w=φ(p^{α-1})=p^{α-2}(p-1)\). Das bedeutet, dass \(w^{p^{α-2}(p-1)}\equiv 1\bmod p^{α-1}\) (das gilt nach Euler immer) und \(w^{p^{α-2}(p-1)/q}\not \equiv 1\bmod p^{α-1}\) ist für alle echten Teiler \(q\) von \(p^{α-2}(p-1)\).

Wir müssen nun zeigen, dass \(n:=\mathrm{ord}_{p^α}(w)\ge φ(p^α)\) ist, dass also \(w^n\not \equiv 1\bmod p^α\) falls \(n=\frac{φ(p^α)}q\) und \(q\) ein echter Teiler von \(φ(p^α)=p^{α-1}(p-1)\) ist. Da aus \(w^n\equiv 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 \equiv 1\bmod p^α \] gilt. Das folgt aber sofort aus der Existenz von \(n_α\), denn \(1+n_αp^{α-1} \not \equiv 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\).

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)))

Aufgabe 5

def ps(n):
    """
    Compute how many percent of the primes below `n` have 2 as a primitive root.
    """
    prime_count = pr = 0
    for p in prime_range(3, n):
        prime_count += 1
        if Mod(2,p).is_primitive_root():
            pr += 1
    return N((pr / prime_count) * 100)

for i in [1..6]:
    print ps(10**i)

Aufgabe 1

crt([13,16], [51,49])

liefert das gesuchte Ergebnis. Alternativ kann man auch den chinesischen Restsatz manuell anwenden. In Sage geht das so:

(13*49*Integer(Mod(49, 51)^-1) + 16*51*Integer(Mod(51, 49)^-1)) % (49*51)

Aufgabe 2

Für Teil 2 erhalten wir aus den Lösungen modulo 71 bzw. 73 die Lösungen modulo 71·73=5183 mithilfe des chinesischen Restsatzes:

xs = [32,37]
ys = [36,41]
for a in xs:
    for b in ys:
        print crt([a, b], [71, 73])

Aufgabe 3

pr = [2,3,5]
counters = [0,0,0]
low = 50000
high = 1000000
for p in prime_range(low, high):
    for i in range(len(pr)):
        if Mod(pr[i], p).multiplicative_order() == p-1:
            counters[i] += 1

print [N(c/(prime_pi(high) - prime_pi(low))) for c in counters]

Aufgabe 5

Mit @interact könnt ihr euch die Primitivwurzeln module einer ausgewählten Primzahl zum Beispiel so anzeigen lassen:

def list_of_primitive_roots(p):
    result = []
    for a in [2..p-1]:
        if Mod(a, p).multiplicative_order() + 1 == p:
            result.append((p,a))
    return result

@interact
def pw_dist(p=slider([a for a in [3..100] if a.is_prime()])):
    show(scatter_plot([(a,1) for (b,a) in list_of_primitive_roots(p)],
                      markersize=10, xmin=0, xmax=p))

Alternativ könnt ihr euch auch die Plots für verschiedene Primzahlen gleichzeitig anzeigen lassen:

scatter_plot(sum([list_of_primitive_roots(p) for p in prime_range(5, 100)], []))

Aufgabe 1

u = Mod(2222222222222298763, 111111111111111111111112)
print u^(-1) * 1111111111111111111

print (16 + 1000) % 24, "Uhr"
print ["Mo", "Di", "Mit", "Do", "Fr", "Sa", "So"][(2 + 1041) % 7]

Aufgabe 2

for a in range(1,23):
    for b in range(1,23):
        # %2d teilt Sage mit, dass wir eine ganze Zahl ausgeben möchten, die
        # genau 2 Zeichen lang ist.  Bei einstelligen Zahlen wird mit
        # Leerzeichen aufgefüllt.  Das darauffolgende Prozentzeichen sagt, dass
        # der Parameter rechts für %2d eingesetzt werden soll.  Das Komma am
        # Ende bewirkt, dass Sage den sonst üblichen Zeilenumbruch unterdrückt.
        print "%2c" % (a*b % 23),
    # Damit das ganze eine Tabelle wird brauchen wir aber doch noch ein paar
    # Zeilenumbrüche.
    print ""

Aufgabe 3

Diese Aufgabe kann natürlich mit dem Chinesischen Restsatz gelöst werden. In Sage geht dies folgendermaßen:

x = crt([1,1,1,1,1,0], [2..7]); print x
# Überprüfen wir noch kurz das Ergebnis:
assert(x%2==1 and x%3==1 and x%4==1 and x%5==1 and x%6==1 and x%7==0)

Wie versprochen gibt es hier meine Sage-Programme zum 1. Übungsblatt. Das Skript findet ihr unter http://data.countnumber.de/ENTH/.

Aufgabe 1

Eine einfache Implementierung des Sieb des Eratosthenes ist die folgende Funktion sieve:

def sieve(n):
    """Berechne die Liste aller Primzahlen 2 <= p < n.
    Die gleiche Funktionalität steht in Sage auch als Funktion primes_first_n
    zur Verfügung.
    """
    candidates = [True] * n
    candidates[0] = candidates[1] = False
    # Da wir nur Vielfache von i ab i² streichen, können wir aufhören, sobald
    # i²>n ist.  Die Funktion range erwartet ganze Zahlen als Argument, so dass
    # wir die Gleitkommazahl, die sqrt(n) zurückgibt runden müssen.  Zudem
    # liefert range(a,b) die Liste
    #   a, a+1, a+2, …, b-2, b-1,
    # weshalb wir noch 1 addieren müssen um auch für Quadratzahlen n das
    # richtige Ergebnis zu bekommen.
    for i in range(2, int(sqrt(n))+1):
        if not candidates[i]:
            # Vielfache von zusammengesetzten Zahlen werden schon vorher als
            # Vielfache der Primteiler gestrichen, mache also mit dem nächsten
            # i weiter.
            continue
        # Streiche alle Vielfachen ab i², die durch i teilbaren Zahlen zwischen
        # i und i² haben einen Primfaktor p < i und wurden somit bereits
        # gestrichen.
        for j in range(i*i, n, i):
            candidates[j] = False
    # 
    return [p for (p, is_prime) in enumerate(candidates) if is_prime]

Nun liefert diese Funktion alle Primzahlen unterhalb einer gegebenen Schranke, gefragt war aber eine Liste der ersten 10 000 Primzahlen.

Das Problem können wir lösen, indem wir sieve mit einem großzügig gewählten n aufrufen, etwa n = 1000000, durch ausprobieren (einfach n größer machen, bis die Liste lang genug ist).

# Berechne die Liste der Primzahlen <10⁶ und wähle mit [:10000] die ersten 10000
# aus dieser Liste aus.
first_10000_primes = sieve(1000000)[:10000]

Wir können auch ein bisschen mogeln und uns von Sage die 10 000. Primzahl p angeben lassen und dann sieve(p) aufrufen:

p = nth_prime(10000)
first_10000_primes = sieve(p)

Aufgabe 2

Für diese Aufgabe war eine Funktion zu implementieren, die die Anzahl der Primzahlen unterhalb einer gegebenen Schranke zählt. Dazu können wir unsere Funktion sieve wiederverwenden.

# Die in Aufgabe 1 erzeugte Liste first_10000_primes enthält alle Informationen,
# die wir nun brauchen.  Wir müssen lediglich Zählen, wie viele der Primzahlen
# unter einer gegebenen Schranke liegen.

# Die vielleicht einfachste Möglichkeit ist die folgende:
def num_primes_below(n):
    """Ermittle, wie viele Primzahlen ≤ n es gibt."""
    if n > first_10000_primes[-1]:
        # Wenn die Grenze größer ist, als die 10000. Primzahl, so können wir
        # nicht sicher sein, ob wir wirklich alle Primzahlen < n kennen.  In
        # diesem Fall beenden wir die Funktion einfach mit einer Fehlermeldung.
        print "Grenze zu groß: %d" % n
        return
    counter = 0
    for p in first_10000_primes:
        if p > n:
            # Wir haben alle Primzahlen < n gezählt, sind also fertig.
            return counter
        counter += 1
    return counter

Aufgabe 4.1

def is_valid_isbn10(n):
    # Zunächst wandeln wir die Zahl n in eine Liste ihrer Ziffern um.  Wie
    # Freitag erwähnt geht das, indem wir zunächst n in eine Zeichenkette
    # umwandeln und anschließend die einzelnen Zeichen wieder in Zahlen:
    #   digits = map(int, str(n))
    # Sage beinhaltet jedoch bereits eine passende Funktion, die jedoch nur
    # für den Sage-Ganzzahltyp Integer, nicht den Python-Ganzzahltyp int zur
    # Verfügung steht:
    digits = Integer(n).digits()
    # Diese Funktion liefert jedoch die Ziffern in in umgekehrter Reihenfolge,
    # weshalb wir anschließend die Liste umdrehen, oder das restliche Programm
    # anpassen müssen.  Der Einfachheit halber drehen wir die Liste um:
    digits.reverse()
    # Reduzieren modulo m kann man in Sage mittels des % Operators:
    return sum([(i+1)*d for (i, d) in enumerate(digits)]) % 11 == 0

def find_16_divisor_isbns(low, high):
    """Bestimme die Liste aller gültigen ISBN-Nummern n mit genau 16 Teilern,
    wobei low ≤ n ≤ high ist.
    """
    result = []
    for n in range(low, high+1):
        if is_valid_isbn10(n) and len(divisors(n)) == 16:
            # An dieser Stelle haben wir eine Lösung gefunden und fügen sie
            # hinten an die Liste an.
            result.append(n)
    return result

lst = find_16_divisor_isbns(1441926535, 1441928509)
# Eigentlich sollte lst[11] das richtige Ergebnis sein (die Indizes starten bei
# 0), aber aus irgendeinem Grund ist bei mir die richtige ISBN-10 die 14.
lst[13]

Lösungen für Aufgaben 20, 21 und 23

Lösungen für Blatt 1

Hier gibt es Dateien zur Vorlesung „Graphentheorie“ im Wintersemester 2014/2015.

• Eine ausführliche Lösung von Aufgabe 25 von Übung 4
• Lösung von Aufgabe 51 von Übung 9
• Lösungen für Aufgaben 53, 54, 56 und 57 von Übung 10
• Lösungen für Aufgaben 60 bis 64 von Übung 11
• Lösungen für Aufgaben 65, 66, 69 und 70 von Übung 12