Forum: Mikrocontroller und Digitale Elektronik Quadratwurzelberechnung optimieren


von Wooog (Gast)


Lesenswert?

Hi,

mangels FPU mache ich die Berechnung von Quadratwurzeln selbst:

- die Zahl "m", aus der die Wurzel gezogen werden soll, wird um 4 Bit 
rechts geshiftet

- das Ergebnis wird in einer Schleife mit sich selbst multipilizert und 
dann (per sukzessiver Approximation) abhängig davon, ob es größer oder 
kleiner als "m" ist immer so weit verändert, dass es sich "m" weiter 
annähhert.

Nach etwa 10 Schleifendurchläufen breche ich dann ab und verwende das 
Ergebnis.

Jetzt ließe sich die Genauigkeit noch deutlich erhöhhen (bzw. die Anzahl 
der maximalen Schleifendurchläufe verringern) wenn ich mit der initialen 
Ermittlung des Startwertes näher am Ergebnis dran wäre als meine 
Division durch 16.

Wie könnte ich hier vorgehen um mit meinem Startwert genauer zu sein?

von Daniel (Gast)


Lesenswert?

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?

von Wooog (Gast)


Lesenswert?

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

von max (Gast)


Lesenswert?

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)

von Uwe Bonnes (Gast)


Lesenswert?

CLZ: _C_ount _l_eading _z_eros

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

von max (Gast)


Lesenswert?

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

von Daniel (Gast)


Lesenswert?

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

von Jens (Gast)


Lesenswert?

Was ist mit dem Heron-Verfahren?

von max (Gast)


Lesenswert?

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.

von DerAlbi (Gast)


Lesenswert?

ööö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.

von DerAlbi (Gast)


Angehängte Dateien:

Lesenswert?

...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.

von Klaus W. (mfgkw)


Lesenswert?

DerAlbi schrieb:
> Außerdem wird man nur bekloppt vollgemeckter hier... blabla. kein bock
> auf Codekritik.

Dann lass es doch.

Oder geh in ein Forum, wo man sagt was man tolles hat, es aber nicht 
zeigt.

von DerAlbi (Gast)


Angehängte Dateien:

Lesenswert?

-.-

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^^

von max (Gast)


Lesenswert?

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.

von DerAlbi (Gast)


Lesenswert?

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

von max (Gast)


Lesenswert?

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

von DerAlbi (Gast)


Lesenswert?

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^^

von DerAlbi (Gast)


Lesenswert?

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 :-)

von DerAlbi (Gast)


Lesenswert?

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

Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.