sqrt() für Integer

[sqrt (Gast)]

Wie funktioniert diese sqrt() Implementierung für Integer?

https://www.mikrocontroller.net/attachment/431705/int_sqrt.c

[Jörg W. (dl8dtl) (Moderator)]

Vielleicht postest du ja auch mal die URL des dazu gehörigen Threads?

[sqrt (Gast)]

Jörg W. schrieb:
> Vielleicht postest du ja auch mal die URL des dazu gehörigen
> Threads?

Ich wünschte die hätte ich noch. Es ging in dem Thread allerdings 
ohnehin nicht um diese Funktionen. Die wurden nur am Rande gepostet.

[Theor (Gast)]

https://github.com/sensorium/Mozzi/blob/master/cogl_sqrti.h

[Theor (Gast)]

http://www.codecodex.com/wiki/Calculate_an_integer_square_root

Scheint aber keine eigentliche Erklärung des mathematischen 
Hintergrundes zu enthalten.

[sqrt (Gast)]

Theor schrieb:
> http://www.codecodex.com/wiki/Calculate_an_integer_square_root
>
> Scheint aber keine eigentliche Erklärung des mathematischen
> Hintergrundes zu enthalten.

Danke! Der Link zur Erklärung hinter diesem Link ist inzwischen zwar 
down, aber mit kleiner Zeitreise findet sich folgendes:
https://web.archive.org/web/20020125100441/http://www.embedded.com:80/98/9802fe2.htm

Scheint Raketenwissenschaft zu sein, die NASA benutzt ihn nämlich auch 
;-) Muss ich mir wachen Auges nochmal gründlich durchlesen

[Theor (Gast)]

Witzig. Die Maschine, die diesen Algo in den 60ern implementierte diente 
als Gegengewicht in einem MG bei Ralleys. :-)))))

[Karl K. (karl2go)]

Achwas Raketenwissenschaft. Geh den Code einfach mal durch.

Beispiel n = 25, 1 Byte

root = 0
rem = 25
place = 0x40 => 64, na dämmerts, genau 8²

while (place > remainder)
  place = place >> 2;
Es wird die Quadratzahl gesucht, die kleiner n ist und ein Quadrat von 
2^k, 64 >> 2 = 16 = 4²

place = 16

    while (place)
solange place nicht Null (in richtigen Programmiersprachen würde man 
schreiben while place > 0)

    if (remainder >= root + place)
25 >= 0 + 16, passt, also

            remainder = remainder - root - place;
25 - 0 - 16 = 9
            root = root + (place << 1);
0 + (16 << 1) = 32, klingt komisch ist aber erstmal so

und immer
        root = root >> 1;
32 >> 1 = 16, wieder zurück
        place = place >> 2;
place = 4 = 2², und die nächstkleinere Quadratzahl probieren

    if (remainder >= root + place)
9 >= 16 + 4, passt nicht, also weiter mit

        root = root >> 1;
16 >> 1 = 8
        place = place >> 2;
place = 1 = 1²

    if (remainder >= root + place)
9 >= 8 + 1, passt, also

            remainder = remainder - root - place;
9 - 8 - 1 = 0, nu gugge, es bleibt kein Rest, klar 25 = 5²
            root = root + (place << 1);
8 + (1 << 1) = 10

und wieder
        root = root >> 1;
10 >> 1 = 5, sieht schonmal gut aus
        place = place >> 2;
place = 0, damit simmer fertisch.


es bleibt
root = 5 = sqrt(25)
rem = 0

Und für große Zahlen geht das genauso, nur langweiliger.

[Egon D. (Gast)]

sqrt schrieb:

> Wie funktioniert diese sqrt() Implementierung
> für Integer?

1. Iterativ.
2. Schätzungsweise mit der binomischen Formel
   (a+2^k)^2 = a^2 + 2^(2k)*a + 2^(2k)

[Yalu X. (yalu) (Moderator)]

Diese Funktion versteht jeder:

uint32_t isqrt32a(uint32_t n) {
  uint32_t root = 0;

  for(uint32_t bitval=0x8000; bitval; bitval>>=1) {
    if ((root + bitval) * (root + bitval) <= n)
      root += bitval;
  }
  return root;
}

Wegen der Multiplikation in jedem Schleifendurchlauf ist diese Methode
für CPUs ohne Multiplizierer nicht so gut geeignet.

Die folgende Funktion tut genau dasselbe:

uint32_t isqrt32b(uint32_t n) {
  uint32_t root = 0, remainder = n;

  for(uint32_t bitval=0x8000; bitval; bitval>>=1) {
    uint32_t deltaremainder = 2 * root * bitval + bitval * bitval;
    if (remainder >= deltaremainder) {
      remainder -= deltaremainder;
      root += bitval;
    }
  }
  return root;
}

Im Gegensatz zu isqrt32a wird hier nur mit Zweierpotenzen multipliziert.

Erklärung:

Statt in jedem Schleifendurchlauf das Quadrat des aktuellen Kandidaten
für root zu berechnen (was eine Multiplikation erfordert), wird der Rest
(remainder), d.h. die Differenz zwischen dem Radikanden und dem Quadrat
des Kandidaten berechnet. Tut man dies rekursiv (d.h. unter Verwendung
des jeweils letzten Werts von remainder), ist bei jeder verbleibenden
Multiplikation mindestens einer der beiden Faktoren eine Zweierpotenz.

Die (nicht sehr komplizierte) Mathematik dahinter:

  remainder     = n - root**2
  root'         = root + bitval
  remainder'    = n - root'**2
                = n - (root + bitval)**2
                = n - root**2 - 2*root*bitval - bitval**2
  deltareminder = remainder - remainder'
                = 2*root*bitval + bitval**2
  remainder'    = remainder - deltaremainder

remainder' darf dabei nicht negativ, d.h deltareminder nicht größer als
remainder werden, sonst war root' zu groß. Das wird in der If-Abfrage
abgeprüft, ggf. wird der alte Wert von root beibehalten.

Man kann nun isqrt32b noch etwas optimieren, so dass

- alle Multiplikationen durch Shifts ersetzt werden,

- alle Shifts eine konstante Shift-Weite ≤ 2 haben (gut für CPUs ohne
  Barrel-Shifter) und

- die ersten paar Schleifendurchläufe, in denen die If-Bedingung
  unabhängig von root falsch ist, in eine eigene, etwas vereinfachte
  Schleife ausgelagert werden.

Heraus kommt die vom TE gepostete Variante.

[Christoph db1uq K. (christoph_kessler)]

https://www.oldcalculatormuseum.com/fridenstw.html
ein paar Fotos des beschriebenen elektromechanischen Rechners, 
Leistenbruch-hervorrufendes Gewicht, die Lampen im Raum werden dunkler 
dazu das Geräusch
"punch-punch-cachunk, punch-punch-cachunk-DING, 
punch-punch-DING-clang-clang"
sehr schön formuliert. Erinnert mich an meinen LO-15 Fernschreiber.

[Signalverarbeiter (Gast)]

Yalu X. schrieb:
> Die folgende Funktion tut genau dasselbe:

"Schriftliches Wurzelziehen", hatten wir das in der Schule immer 
genannt.
Wird das eigentlich noch gelehrt? Oder ist es dem G8 zum Opfer gefallen?

[Jörg W. (dl8dtl) (Moderator)]

Signalverarbeiter schrieb:
> Wird das eigentlich noch gelehrt? Oder ist es dem G8 zum Opfer gefallen?

Wenn überhaupt, wird es nur noch als Kuriosität gelehrt.

Wer braucht das schon noch? Früher hatte man Tafelwerke und 
Rechenschieber (schon da musste man das nicht mehr schriftlich können), 
heute gibt's Taschenrechner.

Einen Sinus oder Tangens kann man auch nicht „schriftlich“ ausrechnen.

[x^2 (Gast)]

Vermutlich lässt sich die Funktion nochmals effizienter abbilden falls 
eine Count Leading Zeros Instruktion (CLZ o.ä.) verfügbar ist.

[sqrt (Gast)]

Yalu X. schrieb:
> Heraus kommt die vom TE gepostete Variante.

Danke, sehr gut erklärt!