Streckenlänge ohne Wurzel berechnen

uatu schrieb am 07.01.2015:Möglicherweise reicht Dir dieses Verfahren: http://netstorage.iar.com/SuppDB/Public/SUPPORT/000419/AN-G-002.pdf .

Das ist mir tatsächlich erst jetzt aufgefallen. Ich muß leider meine Begeisterung über diese schnelle Wurzelfunktion zurücknehmen, die liefert ziemlich schräge Werte.

Beispielsweise die Wurzel aus 9. Die Funktion rechnet daraus (Kommentare dahinter von mir):

unsigned int root(unsigned int x){ // x = 9
unsigned int a,b;
b = x; // b = 9
a = x = 0x3f; // a, x = 0x3f
x = b/x; // x = 9 / 0x3f = 0
a = x = (x+a)>>1; // a, x = 0x3f >> 1 = 0x1f
x = b/x; // x = 9 / 0x1f = 0
a = x = (x+a)>>1; // a, x = 0x1f >> 1 = 0x0f
x = b/x; // x = 9 / 0x0f = 0
x = (x+a)>>1; // x = 0x0f >> 1 = 0x07
return(x); // x = 7

Die Wurzel aus 9 ist also 7, aha... ok, der Fehler liegt nur bei 4, absolut gesehen also recht gering.

Naja, wollte ich nur noch kurz berichten.

Z.

Wie wäre ein etwas verbessertes Verfahren (abgeleitet von der Manhattandistanz):
1
Strecken vergleichen
|x2-x1| >= |y2-y1|
2
wenn ja, dann ist eine Näherung
s= |x2-x1| + |y2-y1| / 2
wenn nein, dann
s= |y2-y1|+ |x2-x1| / 2

Fehler sollte <10% (1,1) bleiben

@erbus

Ja, könnte was dran sein, das werde ich mir mal genauer anschauen.

Z.

@zaeld:

zaeld schrieb:... die liefert ziemlich schräge Werte ...

Der Startwert der Funktion (0x3f, d.h. 63) ist so gewählt, dass er über das gesamte vorgesehene Intervall (0..65535) eine möglichst geringe Abweichung ergibt. Dadurch lässt die Genauigkeit in den Randbereichen des Intervalls nach (mit 9 befindest Du Dich in der Gegend des ersten Zehntausendstels des Intervalls).

Falls das Intervall, in dem sich Deine Werte üblicherweise bewegen, merklich kleiner ist, kannst Du die Genauigkeit bei kleinen Eingangswerten deutlich verbessern (was allerdings mit einem Genauigkeitsverlust bei grossen Werten verbunden ist -- die Du jedoch u.U. gar nicht brauchst)., indem Du einen kleineren Startwert verwendest. Z.B. mit dem Startwert 15 ergibt sich für 9 das korrekte Ergebnis, mit dem Startwert 31 ergibt sich 4. Für grafische Anwendungen, wo vermutlich mit Werten in der Gegend bis 4000 zu rechnen ist, dürfte 15 ein günstiger Startwert sein.

Ich wette meine Tastatur gegen einen Guitar-Hero-Controller, dass Pythagoras besser performt als alle hier angebotenen Alternativen. Von den Lookup-Tabellen mal abgesehen, da ist Pythagoras aber genauer.

@Dahergelaufen: Es geht darum -- die Gründe hatte zaeld ja erläutert -- die Berechnung ohne Nutzung der normalen Quadratwurzel-Funktion durchzuführen. Ich kenne ähnliche Problemstellungen im Zusammenhang mit einfachen (dafür aber auch billigen) CPUs, die keine Hardware-Quadratwurzel-Funktion haben. Ohne Quadratwurzel kein Pythagoras.

@uatu

Ach so, da gibt es eine Hardware-Einschränkung. Teufel, das hatte ich überlesen. Ich dachte, es geht rein um Performance, und da dürfte Pythagoras wirklich schwer zu schlagen sein(vorausgesetzt, natürlich, dass man eine Wurzelfunktion zur Verfügung hat. :) )

@Dahergelaufen: Bei zaeld gibt es keine Hardware-Einschränkung, aber er möchte -- was ich gut verstehen kann -- nicht die ganze entsprechende (Floatingpoint-) Bibliothek wegen der Nutzung einer einzigen Funktion einbinden.

Spricht eigentlich was dagegen die Länge mit Winkelfunktionen zu berechnen?
Es gibt Näherungsformeln arctan und sinus mit sehr kleinen Fehlern.

Also man könnte natürlich auch einfach auf die SSE-Routinen des Prozessors zurückgreifen. Vermutlich ist das schneller, als irgendwelche handgestrickten Lösungen.

Man könnte aber natürlich auch einfach die inverse-square-routine aus QUAKE etwas abändern, um die Wurzel zu berechnen, ohne auf irgendwelche Bibliotheken zurückgreifen zu müssen. Der Code ist fast 1:1 aus Wiki übernommen.

#include <stdio.h> float sqr_approx(float n); int main(void) { int i; for (i=1;i<512;i++) printf("\nsqr(%d): %3.3f",i, sqr_approx((float) i)); printf("\n"); return 0; } float sqr_approx(float n) { int i; float x, y; x=n*0.5f; y=n; i= *(int *) &y; i=0x5f3759df- (i >> 1); y= *(float *) &i; y=y*(1.5f -(x*y*y)); y=y*(1.5f -(x*y*y)); return 1/y; }

Falls die Funktion zu genau sein sollte, bzw. man auch mit weniger Genauigkeit auskommt, dann kann man die Funktion auf eine Iteration beschränken.

/* .... */ y=y*(1.5f -(x*y*y)); /* möglicherweise schon genau genug */ /* y=y*(1.5f -(x*y*y)); */ return 1/y; }

Emodul

@emodul:

emodul schrieb:float x, y;

Es ging gerade darum, die Floatingpoint-Bibliothek nicht einzubinden.

uatu schrieb:Es ging gerade darum, die Floatingpoint-Bibliothek nicht einzubinden.

Schau dir den Code nochmals an. Es wird zwar mit Floating-Point-Zahlen gerechnet, aber es werden eben keine (prozessorintensiven) Funktionen (z.B. aus math.h) eingebunden.

Emodul

@emodul: Ich kenne den Quake-Code. McNeal hatte den Vorschlag hier bereits Anfang Januar gepostet.

uatu schrieb:Ich kenne den Quake-Code. McNeal hatte den Vorschlag hier bereits Anfang Januar

Also eigentlich ist das ja gar nicht von Cormack, er hat es nur in QUAKE verwendet und bekannt gemacht, aber alle berufen sich halt darauf. Und wie McNeal schon schrieb, ist es vermutlich sinnvoller (wie ich auch schon schrieb) sich mal die zur Verfügung stehenden Prozessorfunktionen anzuschauen.

Besser als der von dir gepostete Vorschlag ist der "QUAKE"-Code aber allemal.

Emodul

@emodul: Ich weiss nicht, was so schwer daran zu begreifen ist, dass zaeld nach einer Integer-Implementierung sucht. Es geht nicht um maximale Geschwindigkeit (da kämen ganz andere Überlegungen ins Spiel), sondern nur um eine Lösung, die für den gegebenen Zweck schnell und genau genug ist.

Sofern man sowieso Floatingpoint-Berechnungen einsetzen will, ist der Quake-Algorithmus übrigens in Bezug auf die Performance grundsätzlich nicht sinnvoll, falls -- was bei der hier im Thread behandelten Längenberechnung der Fall ist -- die Quadratwurzel, und nicht der Kehrwehrt der Quadratwurzel gesucht wird. Die Floatingpoint-Division am Ende allein (!) ist (mit gewissen Unterschieden je nach Hardware-Implementierung) etwa genauso rechenintensiv wie die Floatingpoint-Quadratwurzel.

Ich weiss nicht, was so schwer daran zu begreifen ist, dass zaeld nach einer Integer-Implementierung sucht. Es geht nicht um maximale Geschwindigkeit (da kämen ganz andere Überlegungen ins Spiel), sondern nur um eine Lösung, die für den gegebenen Zweck schnell und genau genug ist.

Also wenn ich die Beiträge von zaeld richtig verstehe, dann möchte er primär auf die Wurzelfunktion aus der Bibliothek verzichten und nach Möglichkeit auch nicht mit Floatingpoint-Zahlen rechen, nur um diese dann wieder konvertieren zu müssen.

Dass man keine Funktionen aus math.h (oder was auch immer) nutzen möchte, das kann ich nachvollziehen, aber die Berechnung nur mit Integer-Werten durchzuführen, das halte ich für einen Fehler.

Die nachfolgende Funktion rechnet mit Integerwerten. Ist aber langsamer als die "QUAKE"-Version.

#include <stdio.h> unsigned int int_sqrt(unsigned long n); int main(void) { unsigned long i; for (i=1;i<256;i++) { printf("\ninteger_square_root(%u): %u", i, int_sqrt(i)); } printf("\n"); return 0; } /* taken from: www.codecodex.com/wiki/Calculate_an_integer_square_root */ unsigned int int_sqrt (unsigned long n) { unsigned int c = 0x8000; unsigned int g = 0x8000; for(;;) { if(g*g > n) g ^= c; c >>= 1; if(c == 0) return g; g |= c; } }

Emodul

@emodul: Das ist die Ausgangsversion des Shift-and-Add-Algorithmus, den ich hier erwähnt hatte. Der Algorithmus lässt sich noch erheblich optimieren (u.a. ohne Multiplikation), und wird in dieser Form oft verwendet, um die Quadratwurzel-Funktion auf CPUs ohne Hardware-Division zu implementieren.

Bei 32-Bit-Eingangswerten ist er bei CPUs mit Hardware-Division in der Tat auch in der optimierten Version performancemässig einigen anderen Algorithmen unterlegen. Da der Algorithmus Max-Bit-Zahl-des-Eingangswerts/2 Iterationen (mit nur wenigen, extrem einfachen Operationen) benötigt, reduziert sich die Iterationszahl bei 16-Bit-Eingangswerten (was für den hier im Thread behandelten Fall vermutlich ausreicht) auf 8. Dann dürfte die Performance gegenüber den Alternativen (ausser der Hardware-Quadratwurzel-Funktion selbst) durchaus brauchbar sein.

@zaeld: Das Problem des Integer-Fast-Square-Root-Algorithmus mit kleinen Eingangswerten lässt sich ziemlich gut (besser als mit einer grundsätzlichen Herabsetzung des Startwertes, die ich weiter oben vorgeschlagen hatte) mit einem Startwert von 0x0f für Eingangswerte < 1024, und 0x3f, wie bisher, für alle grösseren Eingangswerte lösen. Die Zeile a = x = 0x3f; müsste also nur mit einer entsprechenden Bedingung ergänzt werden. Das entspricht quasi einer 2-Eintrags-Lookup-Table. Alles andere bleibt unverändert.

Ergänzung zu meinem vorhergehenden Beitrag: Es geht sogar noch besser: Verwendet man für Eingangswerte < 1024 den Startwert 0x0f und für grössere Eingangswerte den Startwert 0x7f (statt 0x3f ) bleibt die absolute Abweichung zum exakten Wert über den gesamten Wertebereich <= 1. Die Verbesserung ist so erheblich, dass man sie evtl. sogar den Autoren des Papers mitteilen könnte.