Quadratwurzelberechnung optimieren

Was für Datentypen verwendest du?
Was ist der (übliche) Wertebereich für die Eingabe?
Wie führst du die sukzessive Approximation durch?
Auf welchem Prozessor soll das schnell laufen?

Daniel schrieb:
> Was für Datentypen verwendest du?

32 Bit signed

> Was ist der (übliche) Wertebereich für die Eingabe?

10240..1024000000, eine Häufung für einen bestimmten, kleineren Bereich 
ist nicht vorhersagbar.

> Wie führst du die sukzessive Approximation durch?

Grob vereinfacht so:

r=m>>4; // -> initialen Wert ermitteln

for (i=0; i<10; i++)
{
   inter=r*r;
   if (inter>m) r_max=r;
   else if (inter<m) r_min=r;
   else return r<<5;
   r=(r_max+r_min)>>1;
}
return r<<5;

> Auf welchem Prozessor soll das schnell laufen?

ARM Cortex M3 - am liebsten aber Portabel in reinem C, ohne Assembler

Beim Cortex M3 hast du natürlich einen Vorteil durch CLZ (sicher auch 
als Intrinsic in deinem C/C++ Compiler verfügbar)

Dann kannst du den Startwert besser annähern als mit /16 in dem du dir 
klar machst, dass bei der Multiplikation zweier Zahlen sich deren 
effektive Länge ca verdoppelt.

Die effektive Länge von m ist 32-CLZ(m) die Länge von Wurzel aus m ist 
daher nicht kleiner als LowerGaussian((32-CLZ(m)) / 2) und nicht größer 
als LowerGaussian((32-CLZ(m)) / 2)+1

Du kannst also deinen Wertebereich auf 1 << (((32-CLZ(m)) / 2)-1) bis 1 
<< ((32-CLZ(m)) / 2) einschränken. Starte in der Mitte, dann hast du 
durchschnittlich den den "kürzeszten" Weg. 3 << (((32-CLZ(m)) / 2)-2)

(hoffe ich hab keine allzugroßen Typos im Schnelltipp)

CLZ: _C_ount _l_eading _z_eros

Aber was genau ist LowerGaussian() und aus welcher Bibliothek kommt die 
Funktion?

Lower Gaussian ist die mathematische Bezeichnung für die Rundung auf die 
nächstgelegene Ganzzahl mit kleinerem Absolutwert.

Man könnte auch einfach m>>((32-CLZ(m))>>1) als erste Näherung nehmen.

Was ist mit dem Heron-Verfahren?

Ist bescheiden.
zwar beliebig wählbar genau (inklusive einer garantierten Genauigkeit 
nach x Schritten) und nciht auf die Quadratwurzel beschränkt, aber hier 
zu aufwendig. Riesiger Minuspunkt: die Division.

öööhm... jemand itneressiert an einer exakten Fixkomma-Sqrt in 176 
Takten auf nem AVR32..? Ist leider TemplateMagie in C++ -> aber um das 
Prinzip zu verstehen.. sollte in C konvertierbar sein. Der Algo benötigt 
genau 1 Division.

...hab mir damals ein PDF darüber gemacht, weil ich es für wirklich sehr 
schnell hielt^^
Code geb ich ungern raus... zumal er den Anforderungen nicht entspricht. 
(->kein C)
Außerdem wird man nur bekloppt vollgemeckter hier... blabla. kein bock 
auf Codekritik.

-.-

Typ - Datentyp (un)signed int/short
Fracbits: Anzahl der Kommastellen-bits
Multiplikation mit floats -> bitte händisch durch fixkommarechung 
ersetzten (macht sonst meine lib)
Fixkommat-Typ-Rechungen und Initialisierungen bitte händisch umändern.

mac -> multiply accumulate  Mac(a,b,c) = a*b+c
div64u... nuja.. damit muss man selbst klarkommen. ich hab eine 
Divisionsroutine geschrieben die mit einem Hardeware-32Bit-Div eine gute 
näherung des Ergebnisses bringt. Was der Compiler bei einem echten 64Div 
einsetzt... dauert so lange wie die wurzel selbst^^

DerAlbi schrieb:
> ...hab mir damals ein PDF darüber gemacht

Sorry, du hast in entscheidenden Punkt einen Fehler:

#1  Heronformel 2mal ineinander eingesetzt.
# bn+1 = 1/2 * (bn + x/bn)  =  (x^2 + bn) / (2bn) //Gleichnamig

beim Umformen sind dir x und bn durcheinander geraten. Es müsste heissen
= (bn^2 + x) / 2bn
Multiplikation von bn mit bn/bn(=1), und Addition beider Brüche.

Das verhindet deine Vereinfachung auf einen Bruch ohne Abhängigkeit vom 
vorhergehenden Iterationsschritt.

Das is echt lachhaft :-D Irgendwie haste recht.... :-D
Kann jemand erklären warums trotzdem funktioniert?
:-)

Das tut es nicht.
Beispiel bei dem das veränderte Verfahren fehlverhält:

x:=16, Startwert: bn:=8
-> bn+1 = 16.5
-> bn+2 = 8.25
-> bn+3 = 16.02
-> bn+4 = 8.49
...
Es konvergiert nicht gegen 4

Ich sterb gleich. Was bitte ist denn da jetzt los.
Ich hab damit nen 3D Renderer gebaut - und das Bild stimmt^^ DAS KANN 
doch gar nicht falschn sein...
HÄÄÄÄÄÄÄÄÄÄÄÄÄÄ

..ich arbeite dran^^

Fehler gefunden :-)
die Formel mit dem gemeinsamen Nenner ist falsch. Also die letzte Zeile.
das was rauskommen soll lautet:

( x^2 + 6*bn^2*x + bn^4 ) / (4*bn*x + 4*bn^3)


BEISPIEL: Wurzel aus 14

ld(14) = 3.   ==> bn = 2^1.5

(14^2 + 48*14 + 64) / (4 * 2^1.5 * 14 + 4*2^4.5) = 3.744

Realer wert: 3.742


Ich weiß echt nicht, warum es funktioniert :-)

Also ich kanns immernoch nicht nachvollziehen :-/ Wie bin ich denn da 
drauf gekommen... die Heron-Formel ist ursprünglich korrekt, und wie 
sich herausgestellt hat...
http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Bakhshali_approximation
ist meine Endformel auch korrekt. Wie kommt man da auf mathematisch 
nachvollziehbaren Weg hin? :-D