Algorithmensammlung: Zahlentheorie: Quadratwurzel mit Rest

Algorithmensammlung: Zahlentheorie

• Euklidischer Algorithmus
• Sieb des Eratosthenes
• Primfaktorisierung
• Fibonacci-Folge
• Quadratwurzel mit Rest
• Goldener Schnitt
• Signum

In der Zahlentheorie ist es oft erforderlich, ganzzahlig Wurzeln großer Zahlen zu ziehen, z. B. um festzustellen, ob eine Zahl eine Quadratzahl ist. Mit Gleitkommazahlen können zwar große Zahlen dargestellt werden, aber sobald eine Zahl eine größere Bitlänge als die Länge der Mantisse hat, reicht die Genauigkeit nicht mehr aus. Bei einer Mantissenlänge von 53 Bit, wie z. B. bei Python, kommt es ab 2⁵⁴ zu Abweichungen. Wenn ein Integer-Datentyp mit beliebiger Genauigkeit, wie in GNU Multiple Precision Arithmetic Library oder in einer Programmiersprache wie Python, zur Verfügung steht, kann die Wurzel aus n durch Iteration ausgehend von einem Startwert $x_{0}$ mit Rest berechnet werden:

$x_{k+1}={\frac {1}{2}}(x_{k}+{\frac {n}{x_{k}}})$ wird wiederholt, bis $0<=n-x^{2}<=x\cdot 2$ ist.

Da $(x+1)^{2}=x^{2}+2\cdot x+1$ ist, kann der Rest $n-x^{2}$ maximal $2\cdot x$ sein.

Pseudocode:

k = bitlength(n)   // entspricht 1 + int(log2(n))
if k mod 2 == 1    // wenn ungerade
    k = k - 1
// Startwert: 
x = (n + 2^k) div (2 * 2^(k div 2))
d = n - x*x
// Iteration bis der genaue Wert erreicht ist:
while d < 0 or d > 2 * x
    x = (x + n div x) div 2
    d = n - x*x
return x, d    // Ergebnis: x = ganzzahliger Anteil der Wurzel, d = Rest

Python

# Berechnung Quadratwurzel von n mit Rest durch Iteration, Parameter: (integer) n >= 0, 
# Rueckgabewert: Tupel (int(Quadratwurzel), Rest)
def isqrtrem (n):
    if not isinstance(n, int): raise TypeError()
    if n < 2:
        if n < 0: raise ValueError("Negative argument for isqrt: " + str(n))
        return n, 0
    # bestimme Bitlänge von n -> k
    # berechne Startwert als Mittelwert aus n und der nächsten 4er-Potenz (4**(k/2)) 
    # dividiert durch die Wurzel derselben Potenz (2**(k/2)):
    k = n.bit_length()
    k ^= (k & 1)                             # wenn ungerade, invertiere das letzte Bit
    x = (n + (1 << k)) >> ((k >> 1) + 1)    # x = (n + 2**k) // (2**(k/2) * 2)
     # Iteration bis der genaue Wert erreicht ist:
    while True:
        x = (x + n//x) >> 1
        d = n - x*x
        if 0 <= d <= (x << 1):    # 0 <= d and d <= (x * 2)
            return x, d           # ganzzahliger Anteil der Quadratwurzel x und Rest d = n - x*x

# ohne Rueckgabe des Rests:
def isqrt(n):
    """Integer square root"""
    return isqrtrem(n)[0]

# Pruefen, ob n eine Quadratzahl ist 
# Rückgabewert: (True, Wurzel), wenn n eine Quadratzahl ist, sonst False
def is_square(n):
    if not isinstance(n, int): raise TypeError()    
    if n < 2:
        if n < 0: return False
        return True, n
    k = 0
    while (n & 3) == 0:
        n >>= 2
        k += 1
    # pruefen, ob Divisionsreste durch 8, 9, 5, 7 in Frage kommen:
    if (n & 7) != 1: return False
    if not ((n % 9) in [0, 1, 4, 7]): return False
    if (n % 5) in [2, 3]: return False 
    if (n % 7) in [3, 5, 6]: return False 
    root, rem = isqrtrem(n)
    if rem: return False 
    return True, root << k    # Ergebnis kann auch mit if is_square(n) ausgewertet werden