Berechnung der ganzzahligen Quadratwurzel

Es musste überprüft werden, ob die Ganzzahl ein Quadrat ist, und wenn ja, die Wurzel berechnen. Und ich möchte dies in ganzzahliger Arithmetik tun. Es ist klar, dass die Newton-Methode in ganzen Zahlen implementiert werden kann, aber sie erfordert bei jedem Schritt eine Division. Aber kann es auch anders sein? Finden Sie das Quadratwurzel-Modulo mit der Potenz von zwei und prüfen Sie, ob es sich um eine gewöhnliche Quadratwurzel handelt.

Sie können sich auf ungerade Zahlen beschränken: Wenn für eine gerade Zahl die Anzahl der niedrigstwertigen Nullstellen ungerade ist, gibt es keine Wurzel. Wenn sie jedoch gerade ist, können Sie die Zahl nach rechts verschieben, die Wurzel der ungeraden Zahl berechnen und die ursprüngliche Anzahl der Nullbits zur linken Hälfte zurückschieben.

Für ungerade N und 2 ^k , k> 3, wenn N ≡ 1 mod 8 ist, dann gibt es 4 verschiedene Wurzeln modulo 2 ^k , andernfalls gibt es keine Wurzeln. Wir brauchen die kleinste der vier Wurzeln von x. Die anderen drei Wurzeln sind 2 ^k - x, 2 ^{k - 1} + x und 2 ^k - 2 ^{k - 1} - x

Ich möchte etwas Ähnliches wie die Berechnung des inversen Modulos 2 ^k - Verdoppelung der Anzahl gültiger Bits pro Iteration.

Angenommen, wir haben bereits eine Wurzel x ₀ von N modulo 2 ^k : N - x ₀ ² = 2 ^k a
Und wir wollen x ₁ = x ₀ + 2 ^k-1 y so finden, dass es in N - x ₁ ² mehr niedrigere Nullbits gibt.
N - (x ₀ + 2 ^{k - 1} y) ² = 2 ^k a - 2 ^k x ₀ * y - 2 2 ^k - ² y ²
Teilen Sie durch 2 ^k : a - x ₀ * y - 2 ^{k - 2} y ²
Und gleich 0 modulo 2 ^k-2 : y = a * x ₀ ^-1 mod 2 ^k-2
Erhalten x ₁ = x ₀ + 2 ^k-1 a * (x ₀ ^-1 mod 2 ^k-2 )
Und schließlich x ₁ = x ₀ + (N - x ₀ ² ) / 2 * (x ₀ ^-1 mod 2 ^k-2 )

Von den k Bits in der nächsten Iteration erhalten Sie 2 (k-1) Bits. Gleichzeitig betrachten wir bei jeder Iteration das Gegenteil zur Wurzel.

Testcode:

uint8_t sqr16(uint16_t n) { if (n % 8 != 1) return 0; uint16_t sqr = (n + 1) / 2; //4 bit uint16_t inv = 2 - sqr; sqr += inv * (n-sqr*sqr)/2; //6 bit inv *= 2 - sqr * inv; sqr += inv * (n-sqr*sqr)/2; //10 bit //inv *= 2 - sqr * inv; if (sqr & 256) sqr = 0u - sqr; sqr = (uint8_t)sqr; // lowest root if (n == sqr*sqr) return sqr; return 0; } 

Wenn wir ein paar Iterationen hinzufügen, erhalten wir die Wurzel von uint_64