Dieses Kapitel stellt keine eigenständige Einführung in die elementare Zahlentheorie dar, sondern eine computermathematische Begleitung. Die Vorgehensweise wird meist sein, Definitionen algorithmisch direkt umzusetzen. Sind die Algorithmen ineffizient, weil sie beispielsweise eine törichte erschöpfende Suche beinhalten, werden daraufhin effizientere Algorithmen entwickelt.
In Python ist die Ganzzahldivision als floor-Division implementiert. Sie unterscheidet sich streng genommen von der Ganzzahldivision in der Zahlentheorie, die als euklidische Division definiert ist. Man kann sie mithilfe der floor-Division wie folgt erhalten:
def sgn(x):
return 0 if x == 0 else -1 if x < 0 else 1
def div_euc(x, y):
return sgn(y)*(x//abs(y))
def mod_euc(x, y):
return x%abs(y)
Wie man sieht, unterscheiden sich floor-Division und euklidische Division bei einem negativen Divisor. Solange keine negativen Divisoren auftreten, wollen wir auf die Benutzung der euklidischen Division verzichten.
Die Teilerrelation »x teilt n«
ist genau dann erfüllt, wenn die Ganzzahldivision
n//x keinen Rest lässt. Man definiert demnach:
def divides(x, n):
return n%x == 0
Man kann die Korrektheit dieser Berechnung anhand der Definition
x | n :⇔ ∃k ∈ ℤ: n = kx
prüfen. Beschränken wir uns zunächst auf positive Zahlen. Eine Quantifizierung über unendlich viele ganze Zahlen ist nicht machbar. Allerdings ist die Existenzquantifizierung für k > n per se nicht erfüllt.
def divides(x, n):
return any(n == k*x for k in range(0, n + 1))
Eine geringfügige Modifikation erlaubt negative Zahlen:
def divides(x, n):
return any(n == k*x for k in range(-abs(n), abs(n) + 1))
Zur Ermittlung der Teilerliste von n filtert man einfach alle Teiler aus dem Bereich von 1 bis einschließlich n aus:
def divisors(n):
return [x for x in range(1, abs(n) + 1) if n%x == 0]
Schließt man die trivialen Teiler 1 und n aus der Teilerliste von n aus, spricht man von der Liste der echten Teiler. Wird von jedem der echten Teiler abermals die Liste der echten Teiler bestimmt, entsteht eine Baumstruktur.
def proper_divisors(n):
return [x for x in range(2, abs(n)) if n%x == 0]
def divisor_tree(n, indent = [], node = ""):
L = proper_divisors(n)
print("".join(indent[:-1]) + node + str(n))
if len(L) != 0:
last = len(L) - 1
for (i, d) in enumerate(L):
if i == last:
divisor_tree(d, indent + [" "], "└─ ")
else:
divisor_tree(d, indent + ["│ "], "├─ ")
Der Wert der Teileranzahlfunktion ist schlicht die Länge der Teilerliste:
def divisor_count(n):
return len(divisors(n))
Entsprechend erhält man die Teilersumme als Summe der Elemente:
def divisor_sum(n):
return sum(divisors(n))
Die beiden Funktionen sind Spezialfälle der Teilerfunktion:
def sigma(k, n):
return sum(d**k for d in divisors(n))
Mit der Teileranzahlfunktion lassen sich die hochzusammengesetzten Zahlen erzeugen.
def highly_composite_numbers(order_max):
order = 1; n = 1; count = 0
while order <= order_max:
c = divisor_count(n)
if c > count:
yield (order, n, c)
count = c; order += 1
n += 1
for (order, n, count) in highly_composite_numbers(10):
print("{:2} {:4} {:4}".format(order, n, count))
Wird zur Berechnung des ggT in direkter Weise dessen Definition herangezogen, ermöglichen Mengenoperationen eine sehr knappe Formulierung. Die Teilerlisten werden dabei als Teilermengen betrachtet, deren Schnittmenge gebildet wird, von der das Maximum zu ermitteln ist.
def gcd(a, b):
return max(set(divisors(a)) & set(divisors(b)))
Klassischer euklidischer Algorithmus:
def gcd(a, b):
if b == 0:
return a
elif a == 0:
return b
elif b < a:
return gcd(a - b, b)
else:
return gcd(a, b - a)
Weil die Funktion endrekursiv ist, kann sie direkt zu einem iterativen Algorithmus transformiert werden:
def gcd(a, b):
if a == 0:
return b
else:
while b != 0:
if b < a:
a = a - b
else:
b = b - a
return a
Moderner euklidischer Algorithmus:
def gcd(a, b):
return a if b == 0 else gcd(b, a%b)
Diese Funktion ist ebenfalls endrekursiv. Durch Transformation zu einem iterativen Algorithmus findet sich:
def gcd(a, b):
while b != 0:
h = a%b; a = b; b = h;
return a
Die Hilfsvariable h kann man in Python entfallen
lassen, indem die betreffende Zeile gegen
a, b = b, a%b
ersetzt wird.
In Pythons Standardbibliothek ist gcd im
Modul math zu finden.
Das kgV zweier Zahlen wird mittels der Formel
ggT(a, b) kgV(a, b) = |ab|
berechnet.
def lcm(a, b):
return abs(a*b)//gcd(a, b)
Das kgV von mehr als zwei Zahlen ist gemäß
from functools import reduce
def lcm(*t):
return reduce(lambda a, b: abs(a*b)//gcd(a, b), t)
zu berechnen.
In Pythons Standardbibliothek ist lcm im
Modul math zu finden.
Eine Zahl ist eine Primzahl, wenn sie genau zwei Teiler besitzt. Demzufolge lässt sich der Primzahltest bereits mittels der Teileranzahlfunktion definieren:
def is_prime(n):
return divisor_count(n) == 2
print([n for n in range(0, 20) if is_prime(n)])
Die Ausgabe:
[2, 3, 5, 7, 11, 13, 17, 19]
Zur Ermittlung der ersten m Primzahlen kann man
Werkzeuge aus dem Modul itertools einsetzen:
from itertools import count, islice
def primes(m):
return list(islice(filter(is_prime, count(0)), m))
Keine Zahl n besitzt einen echten Teiler k, bei dem k2 größer als n ist. Nutzt man diesen Umstand aus, erhält man eine wesentliche Verbesserung des Primzahltests:
def is_prime(n):
if n < 2: return False
k = 2
while k*k <= n:
if n%k == 0: return False
k += 1
return True
Eine weitere Verbesserung des Verfahrens erhält man wie folgt. Ist die Zahl n ungerade, kommt keine gerade Zahl k mehr als Teiler infrage, weil n sonst erst recht gerade wäre. Somit darf ab k=3 um 2 inkrementiert werden.
Ist die Zahl n weder durch zwei noch durch drei teilbar, kommt keine durch zwei oder drei teilbare Zahl k mehr als Teiler infrage, weil n sonst erst recht durch zwei oder drei teilbar wäre. In den Restklassen modulo 6 verbleiben damit nur noch 1 und 5. Somit darf ab k=5 wechselweise um 2 und 4 inkrementiert werden. Dieser Wechsel lässt sich elegant dadurch erzeugen, dass das Inkrement von 6 subtrahiert wird.
def is_prime(n):
if n < 4: return n > 1
if n%2 == 0 or n%3 == 0: return False
k = 5; i = 2
while k*k <= n:
if n%k == 0: return False
k += i;
i = 6 - i
return True
Ferner lässt sich noch die Quadrierung von k aus der
Schleife entfernen, dergestalt dass mit der Funktion isqrt
einmalig die abgerundete Wurzel von n berechnet wird.
Genauer gesagt ist die Bedingung k*k <= n
äquivalent zu k <= isqrt(n).
Mit der allgemeinen Äquivalenz von 0 ≤ x und
0 ≤ ⌊x⌋
findet sich nämlich die äquivalente Umformung
k2 ≤ n ⇔ k ≤ √n ⇔ 0 ≤ √n − k ⇔ 0 ≤ ⌊√n − k⌋ = ⌊√n⌋ − k ⇔ k ≤ ⌊√n⌋.
Die Probedivision verhält sich lediglich für kleine Zahlen vorteilhaft. Für größere Zahlen beschrieb Gary Miller im Jahr 1976 einen Test, den Michael Rabin im Jahr 1980 zu einem probalistischen Test modifizierte.
Die ungerade ganze Zahl n > 2 wird zuänchst in die Form n = 2sd + 1 gebracht. Sei außerdem 1 ≤ a < n. Man nennt n eine starke Pseudoprimzahl zur Basis a, falls mindestens eine der folgenden beiden Bedingungen erfüllt ist:
Andernfalls ist a ein Zeuge dafür, dass n definitiv zusammengesetzt sein muss. Für den Test wird a nun zufällig aus dem Bereich 2 ≤ a ≤ n−2 ausgewählt. Wird dieser Test k mal wiederholt, sinkt die Wahrscheinlichkeit dafür, dass eine zusammengesetzte Zahl fälschlich als prim gesehen wird, auf 4−k ab. Wir wählen k = 128 als Vorgabewert.
from random import randint
def is_witness(a, n, d, s):
if pow(a, d, n) == 1:
return False
else:
for r in range(0, s):
if pow(a, d*2**r, n) == n - 1: return False
return True
def is_probably_prime(n, count = 128):
if n == 2 or n == 3:
return True
elif n%2 == 0 or n < 2:
return False
s = 0; d = n - 1
while d%2 == 0:
s += 1; d = d//2
for i in range(count):
a = randint(2, n - 2)
if is_witness(a, n, d, s): return False
return True
Der Test-Algorithmus lässt sich noch ein wenig optimieren:
def is_witness(a, n, d, s):
x = pow(a, d, n)
if x == 1 or x == n - 1:
return False
else:
for i in range(s - 1):
x = (x*x)%n
if x == n - 1: return False
return True
Bei der Faktorisierung per Probedivision wird der kleinste Primteiler der Zahl ermittelt und so oft abgespalten wie dieser enthalten ist. Der Vorgang wird anschließend abermals auf die verbleibende Zahl angewendet. Das geht solange weiter, bis keine Primteiler mehr enthalten sind bzw. bis die verbleibende Zahl gleich eins geworden ist.
def factor(n):
assert n > 0
acc = []
p = 2
while n != 1:
e = 0
while n%p == 0:
n = n//p; e += 1
if e != 0: acc.append((p, e))
p += 1
return acc
Analog zum Primzahltest erhält man eine erhebliche Verbesserung, wenn man die Bedingung n ≠ 1 bzw. p ≤ n gegen p2 ≤ n ersetzt. Ist diese Bedingung nicht erfüllt, muss es sich bei n im Fall n ≠ 1 um einen Primfaktor handeln.
def factor(n):
acc = []
p = 2
while p*p <= n:
e = 0
while n%p == 0: n //= p; e += 1
if e != 0: acc.append((p, e))
p += 1
if n != 1: acc.append((n, 1))
return acc
Auch hier erhält man durch Beschränkung auf die Restklassen 1 und 5 modulo 6 nochmals eine geringfügige Beschleunigung.
def trial_div(acc, n, p):
e = 0
while n%p == 0: n //= p; e += 1
if e != 0: acc.append((p, e))
return n
def factor(n):
acc = []
n = trial_div(acc, n, 2)
n = trial_div(acc, n, 3)
p = 5; i = 2
while p*p <= n:
if n%p == 0:
n = trial_div(acc, n, p)
p += i; i = 6 - i
if n != 1: acc.append((n, 1))
return acc
Es bezeichnet σk(n) die Teilerfunktion. Bedeutsam sind ebenso ihre beiden Spezialfälle,
Zur Berechnung lässt sich ohne Weiteres die Definition der Teilerfunktion einsetzen:
def sigma(n, k = 1):
return sum(d**k for d in range(1, n + 1) if n%d == 0)
Man kann im Fortgang die Primfaktorzerlegung von n einsetzen. Daraufhin lässt sich die Summe faktorisieren, wobei die Faktoren Partialsummen von geometrischen Folgen darstellen. Setzt man noch jeweils die geometrische Summenformel ein, gelangt man zu:
def sigma(n, k = 1):
if k == 0:
return prod((e + 1) for (p, e) in factor(n))
else:
return prod((p**(k*e + k) - 1)//(p**k - 1) for (p, e) in factor(n))
Im Fall k = 0 kommt die geometrische Summenformel nicht zur Anwendung, bzw. es wäre der Grenzwert für k → 0 zu betrachten.
Der Wert φ(n) ist definiert als die Ordnung der Einheitengruppe des Restklassenrings ℤ/nℤ. Mit anderen Worten handelt es sich um die Anzahl der multiplikativ invertierbaren Restklassen modulo n. Demzufolge ist der Wert wie folgt zu ermitteln:
def is_invertible(x, n):
e = 1 if n != 1 else 0
return any((x*y)%n == e for y in range(0, n))
def units(n):
return [k for k in range(0, n) if is_invertible(k, n)]
def phi(n):
return len(units(n))
Bei n = 1 handelt es sich um einen Sonderfall, bei dem das additiv neutrale Element mit dem multiplikativ neutralen Element zusammenfällt.
Eine elementare Gesetzmäßigkeit ist, dass eine Zahl k genau dann multiplikativ invertierbar modulo n ist, wenn k koprim zu n ist, wenn also ggT(k, n) = 1 gilt. Daher gilt:
def phi(n):
return len([k for k in range(0, n) if gcd(n, k) == 1])
Der Wert lässt sich mit der Primfaktorzerlegung auch über eine Produktformel ermitteln:
from math import prod
def phi(n):
return prod(p**(e-1)*(p - 1) for (p, e) in factor(n))
Setzt man die Probedivision zur Primfaktorzerlegung ein, findet sich:
def phi(n):
assert n > 0
y = 1; p = 2
while p*p <= n:
e = 0
while n%p == 0: n //= p; e += 1
if e != 0: y *= p**(e-1)*(p - 1)
p += 1
if n != 1: y *= (n - 1)
return y
Der Wert λ(n) ist definiert als der Gruppenexponent der Einheitengruppe des Restklassenrings ℤ/nℤ. Die algorithmische Umsetzung dieser Definition:
def units(n):
return [x for x in range(0, n) if gcd(n, x) == 1]
def carm(n):
e = 1 if n != 1 else 0
m = 1; G = units(n)
while True:
if all(pow(g, m, n) == e for g in G): return m
m += 1
Die Carmichael-Funktion wurde carm statt
lambda genannt, weil lambda in Python
ein Schlüsselwort ist.
Der Wert der Möbiusfunktion ist gemäß ihrer Definition folgendermaßen zu ermitteln:
def mu(n):
assert n > 0
a = factor(n)
if any(e != 1 for (p, e) in a):
return 0
elif len(a)%2 == 0:
return 1
else:
return -1
Zwei Zahlen a, b heißen kongruent modulo m, kurz a ≡ b (mod m), wenn die Division durch m bei beiden Zahlen den gleichen Rest lässt. Weil dieser Rest das Ergebnis der Modulo-Operation ist, kann die Kongruenzrelation als
def is_cong(a, b, m):
return a%m == b%m
implementiert werden. Zwei Zahlen sind genau dann kongruent modulo m, wenn ihre Differenz durch m teilbar ist, womit man die Kongruenzrelation alternativ als
def is_cong(a, b, m):
return (a - b)%m == 0
implementieren kann. Diese Variante besitzt den Vorteil, eine Division weniger ausführen zu müssen, deren Berechnungen aufwendig sind.
Aufgabe sei es zunächst, sämtliche Lösungen der Kongruenz
x ≡ a (mod m)
zu bestimmen. Erfreulich ist hierbei der Umstand, dass der
Restklassenring ℤ/mℤ nur endlich viele Elemente besitzt,
deren kanonische Repräsentanten den Bereich
range(0, m) bilden. Ist m hinreichend
klein, kann die Lösungsmenge also durch erschöpfende Suche ermittelt
werden. Ohne Frage ist das töricht, weil die eindeutige Lösung
sinnfällig direkt als a%m vorliegt.
Nun betrachten wir aber die Aufgabe, dass die beiden Kongruenzen
x ≡ a (mod m),
x ≡ b (mod n)
gleichzeitig erfüllt sein sollen. Ein solches Problem bezeichnet man als System von Kongruenzen oder als simultane Kongruenz. Ähnlich wie bei einem Gleichungssystem handelt es sich schlicht um eine logische UND-Verknüpfung von Kongruenzen. Man kann es sogar als Gleichungssystem
[x]m = [a]m,
[x]n = [b]n
betrachten. Die kanonische Projektion
[x]m wird als x%m kodiert.
Es ist nun so, dass diese Projektion periodisch mit Periodenlänge
m ist. Infolge ist die Abbildung
x ↦ ([x]m, [x]n)
periodisch mit Periodenlänge kgV(m, n). Die erschöpfende Suche darf also als
def solve(a, m, b, n):
return [x for x in range(0, lcm(m, n))
if is_cong(x, a, m) and is_cong(x, b, n)]
implementiert werden. Die allgemeine Formulierung
def solve(s):
M = lcm(*(m for a, m in s))
return [x for x in range(0, M)
if all(is_cong(x, a, m) for a, m in s)]
löst ein System von beliebig vielen Kongruenzen.
Wie bereits im Abschnitt zur phi-Funktion ausgeführt, sind die Einheiten genau die zum Modul koprimen Repräsentanten. Die Berechnung der Einheitengruppe des Restklassenrings ℤ/nℤ kann somit als
def units(n):
return [x for x in range(0, n) if gcd(n, x) == 1]
erfolgen. Die Inverse eines Repräsentanten wird durch
def inv(x, n):
e = 1 if n != 1 else 0
for y in range(0, n):
if (x*y)%n == e: return y
berechnet. Das folgende Programm erstellt die Gruppentafel der Einheitengruppe:
def print_group_table(n, space = 2):
G = units(n)
fmt = lambda x: "{:>{}}".format(x, space)
print(fmt("*") + " | " + " ".join(fmt(y) for y in G))
print("─"*(space + 2 + (space + 1)*len(G)))
for x in G:
print(fmt(x) + " | " + " ".join(fmt((x*y)%n) for y in G))
Für manche n ist die Einheitengruppe von ℤ/nℤ eine zyklische Gruppe. Das Erzeugnis ⟨g⟩ wird durch
def generate(g, n):
acc = set()
x = 1 if n != 1 else 0
while not x in acc:
acc.add(x)
x = (x*g)%n
return acc
berechnet. Existiert ein Erzeuger, das heißt, ein g mit ⟨g⟩ = G, dann ist G zyklisch:
def is_cyclic(n):
G = units(n); setG = set(G)
return any(generate(g, n) == setG for g in G)
Die Erzeuger einer zyklischen Einheitengruppe bezeichnet man als Primitivwurzeln. Wir können sie mit
def generators(n):
G = units(n); setG = set(G)
return [g for g in G if generate(g, n) == setG]
berechnen.
Für manche Aufgaben benötigt man die Inverse einer gegebenen Funktion f: ℤ/mℤ → ℤ/mℤ. Weil f dabei nicht als Injektion vorausgesetzt werden soll, muss die Inverse immer das Urbild ihres Arguments berechnen. Für kleine m lässt sich dieses Problem in kunstloser Weise lösen, indem einmalig die gesamte invertierte Wertetabelle berechnet wird.
def inv(f, m):
d = {}
for x in range(m):
y = f(x, m)
if y in d:
d[y].append(x)
else:
d[y] = [x]
def finv(y):
return d[y] if y in d else []
return finv, sorted(d.keys())
def tabulate_inv(f, m):
finv, img = inv(f, m)
space = len(str(m - 1))
for n in img:
print("{:{}}: {}".format(n, space, finv(n)))
Beispielsweise berechnet
def sq(x, m):
return (x*x) % m
tabulate_inv(sq, 24)
die Wertetabelle der Quadratwurzel modulo 24. Das Programm erzeugt die Ausgabe:
0: [0, 12] 1: [1, 5, 7, 11, 13, 17, 19, 23] 4: [2, 10, 14, 22] 9: [3, 9, 15, 21] 12: [6, 18] 16: [4, 8, 16, 20]
Anhand der invertierten Wertetabelle lassen sich insofern die quadratischen Reste bestimmen. Man bezeichnet ein a als quadratischen Rest modulo m, wenn die Kongruenz x2 ≡ a (mod m) eine nichtleere Lösungsmenge besitzt. Das Programm
def quadratic_residues(m):
_, img = inv(sq, m)
coprime_t = [a for a in img if gcd(a, m) == 1]
coprime_f = [a for a in img if gcd(a, m) != 1]
return coprime_t, coprime_f
for m in range(1, 20):
coprime_t, coprime_f = quadratic_residues(m)
print("mod {:2}: {}, {}".format(m, coprime_t, coprime_f))
berechnet die Tafel der quadratischen Reste, wobei zum Modul teilerfremde quadratische Reste von den übrigen getrennt als Erstes aufgelistet werden.
Die Bisektion ist ein gutmütiges, vergleichsweise leicht einsichtiges numerisches Verfahren. Wir können sie durch Anpassung auch zur Berechnung der abgerundeten Wurzel einsetzen. Im Unterschied zur Bisektion bei reellen Funktionen liefert sie hier das exakte Ergebnis in einer endlichen Anzahl von Schritten.
def isqrt(n):
assert isinstance(n, int) and n >= 0
a = 0; b = n + 1
while a != b - 1:
m = (a + b)//2
if m*m <= n:
a = m
else:
b = m
return a
Beim Eintritt in die Funktion soll n ≥ 0 als Vorbedingung gelten. Die zielführende Schleifeninvariante ist a2 ≤ n < b2. Mit dem Austritt aus der Schleife gilt nun zusätzlich b = a + 1, und somit
a2 ≤ n < (a + 1)2 ⇔ a ≤ √n < a + 1 ⇔ 0 ≤ √n − a < 1 ⇔ a = ⌊√n⌋.
Im letzten Schritt wurde genutzt
y = ⌊x⌋ :⇔ y ∈ ℤ und 0 ≤ x − y < 1.
Wie gewünscht gilt a = ⌊√n⌋
beim Verlassen der Funktion. Es wurde isqrt nach einigen
Jahren dem Modul math hinzugefügt. Man erhält demgemäß
mit
assert all(isqrt(n) == math.isqrt(n) for n in range(0, 10000))
noch kurzerhand einen flüchtigen Realitätsabgleich.
Obgleich die Bisektion eine logarithmische Zeitkomplexität besitzt, kommt eine Anpassung des Heron-Verfahrens zumindest für große Zahlen noch wesentlich schneller zum Ziel. Eine schnell zu berechnende, möglichst gute Startnäherung verleiht dem Verfahren den Feinschliff.
def isqrt(n):
assert isinstance(n, int) and n >= 0
if n == 0: return 0
x = 1 << (n.bit_length() + 1 >> 1)
while True:
y = (x + n//x) >> 1
if y >= x:
return x
x = y
Die zielführende Schleifeninvariante ist √n < x + 1. Dass sie von der Zuweisung x := y unberührt bleibt, bestätigt die aus 0 ≤ (x − √n)2 gefolgerte Abschätzung
√n ≤ (x + n/x)/2 < ⌊(x + n/x)/2⌋ + 1 = ⌊(x + ⌊n/x⌋)/2⌋ + 1.
Beim Austritt aus der Funktion gilt außerdem
x ≤ y = ⌊(x + n/x)/2⌋ ≤ (x + n/x)/2,
was man für x > 0 zu x2 ≤ n umformen kann. Insgesamt gilt also beim Austritt wie gewünscht
x ≤ √n < x + 1 ⇔ 0 ≤ √n − x < 1 ⇔ x = ⌊√n⌋.
Schließlich ist noch zu prüfen, dass die Schleifeninvariante auch beim Eintritt in die Schleife erfüllt ist. Hierzu wird die Startnäherung in die Invariante eingesetzt. Weiterhin besteht die Beziehung
bit_length(n) = ⌊log2(n)⌋ + 1.
Zu prüfen ist also, dass jedes n ≥ 1 der Ungleichung
√n < 2⌊(⌊log2(n)⌋ + 2)/2⌋ + 1 = 2⌊log2(n)/2⌋ + 1 + 1
genügt. Sie lässt sich in die logarithmierte äquivalente Form
log2(√n − 1) − 1 < ⌊log2(n)/2⌋
bringen. Aufgrund der Identität x − 1 < ⌊x⌋ folgt die Ungleichung aus
log2(√n − 1) ≤ log2(n)/2 ⇔ (√n − 1)2 ≤ n ⇔ n ≥ 1/4.