math.h - wo finde ich den Quelltext zur Funktionsweise vom Wurzelziehen?

[Hansen (Gast)]

Hallo und schönen Ostermontag!

Ich bin auf der Suche nach der Funktionsweise vom Wurzelziehen in der 
math.h?

Hintergrund: Ich möchte das Wurzelziehen umgehen und meine Rechnung 
stattdessen mit einer Polynomanpassung, welche ich vorher in Excel 
berechne, lösen. Es geht um ein Termoelement - die Rechnung dazu richtet 
sich nach einem Wiki-Artikel, welches ein Ersatzpolynom erstellt.

Ich würde daher gerne den rechnerischen Aufwand für den Mikrocontroller 
abschätzen können. Kann ich irgendwo den Quelltext dazu sehen? Ich finde 
in der math.h nicht den verweis auf eine math.c oder ähnliches.

Kann mir evtl. einer sogar sagen, in welchem Rahmen sich der Aufwand für 
den Mikrocontroller beim Wurzelziehen bewegt?

Dazu vielleicht noch im Vergleich, ob der mathematische Aufwand von 
Wurzelziehen gegen Polynomen in Fließkommaarithemtik größer/kleiner ist?


Bin um jeden Tip dankbar!

[... (Gast)]

http://sourceware.org/newlib

[... (Gast)]

PS:

Könnte auch sein, dass die Funktion sqrt direkt im Quellcode vom gcc zu 
finden ist.

[Johann L. (gjlayde)]

Die math.h enthält nur Prototypen. Da gibt es keine "Funktionsweise".

Je nach verwendeter libc/libm-Implementierung (glibc, newlib, ...) 
können unterschiedliche Algorithmen zur Anwendung kommen.

Mögliche Kandidaten sind z.B. CORDIC, Fixpunktiteration, ...

Der Aufwand auf einem µC richtet sich nach der konkreten Implementierung 
und den Fähigkeiten, die der µC mitbringt (FPU, ...)

Bei einer reinen Software-Implementierung der Wurzel liegt der Aufwand 
in der Größenordnung einer Software-Implementierung einer Division (bis 
auf O-Konstanten).

[Johann L. (gjlayde)]

... schrieb:
> PS:
>
> Könnte auch sein, dass die Funktion sqrt direkt im Quellcode vom gcc zu
> finden ist.

Nein. Entweder in der libc-Implementierung oder wenn was zur Compilezeit 
benötigt wird (Wurzeln aus Konstanten) in gmp/mpfr oder für komplexe 
Wurzeln mpc.

[Hansen (Gast)]

Danke schonmal für die Antworten!

Habe mir jetzt mal Cordic und Fixpunktiteration angeguckt - da muss man 
sich ja wirklich erstmal reinlesen um das zu verstehen.

Daher muss ich nochmal anders fragen - kann mir jemand rein aus 
Erfahrung sagen, ob man das Wurzelziehen auf einem uC OHNE FPU vermeiden 
sollte?

[Johann L. (gjlayde)]

Hansen schrieb:

> Daher muss ich nochmal anders fragen - kann mir jemand rein aus
> Erfahrung sagen, ob man das Wurzelziehen auf einem uC OHNE FPU vermeiden
> sollte?

Welche Antwort erwartest du auf so eine Wischi-Waschi-Frage?

Das hängt ab von

- Dem µC
- Wieviel Resourcen der hat
- Von der zu lösenden Aufgabe
- Von den Nebenbedingungen (Wertebereich, Genauigkeit,
  gibt es Alternativen, sind die teurer, ...?)

[Michelle K. (Firma: electronica@tdnet) (michellekonzack)]

Also ich bin hier beim Aufbau eines MPPT Solar-Converters mit
einem *NXP Cortex M0 LPC11C14* und der ist schneller, als der
Zugriff auf eine Tabelle im I²C NVRAM/EEPROM.

Grüße
Michelle

[Tobias (Gast)]

Ein interessanter Artikel ist da auch "Magical Square Root 
Implementation In Quake III" [1]. Ich habe aber noch nie ausprobiert, 
wie genau die Werte sind.


[1] http://www.codemaestro.com/reviews/9

[Michel (Gast)]

Wenn du wissen möchtest, wieviel Aufwand das Wurzelziehen ist, bestimme 
doch im Simulator die Zeit, die ein typische (tm) Wurzeloperation 
benötigt.

[Jobst M. (jobstens-de)]

Y(n+1) = ((X/Yn)+Yn)/2

Das wendet man entweder so lange an, bis sich keine Änderung mehr ergibt 
oder lässt es einfach eine bestimmte Dauer laufen.


Auf einem 8051 (4MIPS) schaffe ich es, geschätzt ca. 10 64-Bit 
Festkommazahl-Wurzelziehungen pro Sekunde zu machen.


Gruß

Jobst

[W.S. (Gast)]

Besorge dir mal das steinalte Buch "Algorithmen der Mikrorechentechnik" 
(oder so ähnlich) vom Autorenteam Lampe, Jorke, Wengel.

Die haben dort sowohl ne komplette Gleitkommaarithmetik (single, also 32 
Bit) als auch das Wurzelziehen durchexerziert. Als Basis diente damals 
der Z80, aber das hat ja nix mit den Prinzipien zu tun. Als Aufwand 
für's Wurzelziehen kann man so etwa den für eine Gleitkommadivision 
annehmen. Es geht also recht schnell.

W.S

[Yalu X. (yalu) (Moderator)]

Hansen schrieb:
> Hintergrund: Ich möchte das Wurzelziehen umgehen und meine Rechnung
> stattdessen mit einer Polynomanpassung, welche ich vorher in Excel
> berechne, lösen.

Wenn du die Wurzelberechnung mittels einer von Excel erzeugten Polynom-
anpassung umgehst, wozu brauchst du dann noch den Wurzelalgorithmus?

[Joachim (Gast)]

Die schöne Michelle schrieb:

>Also ich bin hier beim Aufbau eines MPPT Solar-Converters mit
>einem NXP Cortex M0 LPC11C14 und der ist schneller, als der
>Zugriff auf eine Tabelle im I²C NVRAM/EEPROM.

>Grüße
>Michelle

Hast aber vergessen zu sagen, daß das 'ne 32-bit Familie ist. Für mich 
hörte sich die Frage eher danach an ob es auf einem 8-bit System Sinn 
macht.

[MoTheG (Gast)]

Zum Wurzelziehen von Ganzzahlen/Integers:
http://www-user.tu-chemnitz.de/~heha/Mikrocontroller/Wurzel.htm
Hinzu kommen noch zwei Verfahren: Heron

$$Y(n+1)=\frac{\frac{X}{Y(n)}+Y(n)}{2}$$

und Ausprobieren mit Teile-und-Herrsche.
Bei beiden Verfahren kann man Schätzwerte, Obergrenze, Untergrenze durch 
Halbierung der Stellenanzahl erzeugen. (Kostet aber auch Zeit O(n)=n)

[Moritz G. (mbydq2)]

im Beitrag #2171280:

> http://www-user.tu-chemnitz.de/~heha/Mikrocontroller/Wurzel.htm

Von den dort dargestellten Verfahren, halte ich nur das 
Doppelwurzel-Verfahren (keine Ahnung wie es wirklich heißt) für 
sinnvoll.
Da µC heute in Hardware multiplizieren sollten (die paar Transistoren 
kosten doch keinen Unterschied), dürfte systematisches Ausprobieren 
trotz des Mangels an Eleganz nicht schlecht dar stehen.
Tabellen sind viel zu teuer für µC, ich halte sie nur für sinnvoll, wenn 
die Funktion dem Programmierer überhaupt nicht bekannt ist.

[Hagen R. (hagen)]

>Tabellen sind viel zu teuer für µC, ich halte sie nur für sinnvoll, wenn
>die Funktion dem Programmierer überhaupt nicht bekannt ist.

Aha, hast du deinen Satz mal verstanden ?

Du hälst es nur für sinnvoll das der Programmierer Lookuptabellen
benutzt wenn der Programmierer selber, der sie jetzt benutzen möchte,
garnicht weiß wie und was er da zu benutzen hat.

Na super, da habe ich jahrelang wohl was falsch gemacht, hätte ich
vorher gewusst das der Shit den ich verzapfe garnicht funktionieren muß.

Eine MCU die durchschnittlich 100 Takte pro arithmetischen Befehl
benötigt im Vergleich dazu aber nur 1 Takt um eine 32Bit Variable aus
dem Speicher zu laden und von diesem Speicher noch 256kB hat, wird bei
fast allen arithmetischen komplexeren Formeln mit Lookuptabellen die
performanteste Lösung erzeugen. Warum ? Weil es nicht vom Programmierer
direkt abhängt sondern objektiv vom Aufwand/Nutzen Verhältnis, das uns
die Architektur der MCU aufzwingt. Je mehr Erfahrung und Kenntnisse der
Programmierrer nun hat je besser wird er auf dem Instrument spielen
können.


1.) deine Herleitung ist falsch
2.) deine Annahme ein Entwickler mit mögl. wenig Wissen wäre das beste
ist ebenfalls falsch

Gruß Hagen

[Moritz G. (mbydq2)]

Hagen Re schrieb:

> Aha, hast du deinen Satz mal verstanden ?
 Sie sollten sich zunächst selbst fragen, ob Sie meinen Satz verstanden 
haben. Fragen Sie das Nächste-mal doch einfach nach bevor Sie OT gehen.

> Eine MCU die durchschnittlich 100 Takte pro arithmetischen Befehl
> benötigt im Vergleich dazu aber nur 1 Takt um eine 32Bit Variable aus
> dem Speicher zu laden und von diesem Speicher noch 256kB hat, wird bei
> fast allen arithmetischen komplexeren Formeln mit Lookuptabellen die
> performanteste Lösung erzeugen.
 Motto: Ich konstruiere mir ein bizarres Gegenbeispiel.
Ich dachte eher an einen 8-bit µC der <256KB Speicher hat, nicht an 
einen 32-bit Mikroprozessor mit, in Ihrem Fall, offenbar sehr schnellem 
Speicher.

[Moritz G. (mbydq2)]

Für

$$a \in (1;1{,}5)$$

geht

$$b = \frac{a-1}{2}+1$$

und dann das übliche

$$r = \frac{b + \frac{a}{b}}{2} \approx \sqrt{a}$$

nur ein mal.
Das klingt erst ein mal nach einem unnützen Spezialfall, ist aber sehr 
nützlich, wenn man die

$$\sqrt[2^n]{x}$$

auf viele Nachkommastellen genau berechnen will, was wiederum nützlich 
ist, wenn man einen Logarithmus berechnen will.

[Moritz G. (mbydq2)]

Moritz G. schrieb:

> Das klingt erst ein mal nach einem unnützen Spezialfall, ist aber sehr
> nützlich

..., da alle binären zahlen leicht auf <2 zu schieben sind, und viele 
auf zwischen 1 und 1,5_d = 1,1_b zu schieben sind, nämlich alle mit 
10...

Ich habe es nicht getestet aber man könnte es vermutlich auch so:

//2.Wurzelziehen
int sqrt(int r){
//führende nullen zählen; nur für 16bit=lengthof(r)
   char n=0;
   int t=r;
   if(r > 0) //oberste Stelle ist frei
    for( ; t > 0; ++n) t <<= 1;
    n=16-n;
  t=r;
    if(n&&1){
   r-=1<<--n;//ginge auch durch XOR
   r>>=1;
   r+=1<<--n;//ginge auch durch XOR
   r>>=n>>1;
   }
  else{
   r>>=3;//(2 hoch 3)~~(3 hoch 2)
   n-=3;
   r-=1<<--n;//ginge auch durch XOR
   r>>=1;
   r+=1<<--n;//ginge auch durch XOR
   r>>=n>>1;
   r*=3;
   }
  r+=r+t/r>>1;//einzige Division; sollte unter 1% Abweichung haben.
  return(r);

machen.
Eine Division ist immernoch nötig.

[Moritz G. (mbydq2)]

Moritz G. schrieb:

Schon Fehler gefunden: (ich brauche dringend Linux + GCC damit ich mal 
auf die Schnelle ein Test-konsolen-programm schreiben kann)

>   if(r > 0) //oberste Stelle ist frei
>   for( ; t > 0; ++n) t <<= 1;
 Kein Abbruch bei ungültigen Werten.

>   else{
>    r>>=3;//(2 hoch 3)~~(3 hoch 2)
>    n-=3;
>    r-=1<<--n;//ginge auch durch XOR
>    r>>=1;
>    r+=1<<--n;//ginge auch durch XOR
>    r>>=n>>1;
>    r*=3;
>    }
 "--n" geht hier nicht, weil n dann zu klein für folgendes n>>1.

int sqrt(int r){
//führende nullen zählen; nur für 16bit=lengthof(r)
   char n=0;
   int t=r;
   if(r > 0) for( ; t > 0; ++n) t <<= 1;//oberste Stelle ist frei
   else return(0);//Fehler 1/2
   n=16-n;
   t=r;
   if(n&&1){
    r-=1<<--n;//ginge auch durch XOR
    r>>=1;
    r+=1<<--n;//ginge auch durch XOR
    r>>=n>>1;
   }
   else{
    r>>=3;//(2 hoch 3)~~(3 hoch 2)
    n-=3;
    r-=1<<n-1;//ginge auch durch XOR - Fehler 2/2
    r>>=1;
    r+=1<<n-1;//ginge auch durch XOR - Fehler 2/2
    r>>=n>>1;
    r*=3;
   }
   r+=r+t/r>>1;//einzige Division; sollte unter 1% Abweichung haben.
   return(r);

Es geht mir auch mehr um das Prinzip. Ich wollte keine C-Funktion 
veröffentlichen, der Algorithmus ist so gut zu beschreiben.
Für große Werte geht auch eine Gerade als erste Näherung g=mx+b 
m=a/4096.

[Arne (Gast)]

http://supp.iar.com/Support/?note=18180

[Moritz G. (mbydq2)]

Das Verfahren aus den Beiträgen von:
06.05.2011 21:50
07.05.2011 02:48
hat einen fiesen Haken:
Es funktioniert nicht für Werte wie:
0000.001x_b, 0000.1xxx_b, 001x.xxxx_b
erst bei 1xxx.xxxx_b kommt etwas sinnvolles heraus, aber das ist leider 
bereits eine negative zahl und damit dann doch wieder falsch.
=> Das Verfahren ist für Werte < 2^15 ungeeignet. Werte > 2^15 lassen 
sich aber gut durch wenige Geraden annähern. Damit ist das Verfahren nur 
ein interessantes Kuriosum, da es nur bei unganzzahlig Stelligen 
Binärwerten funktioniert. Wenn es um Fließkommazahlen geht, ist die 
Sache wieder eine andere.

Ich komme zu dem Schluss, dass auf den folgenden Intervallen folgende 
Verfahren die beste erste Näherung bieten:
[0,30]: gezielt aus 2+2^2=6 Möglichkeiten Auswählen
(30,65535]: Geraden, intervallhalbierendes Ausprobieren
Exemplarisch: (25,128)

$$w = \frac{\frac{235 \cdot r}{256} + 69}{16}$$

Die durchschnittliche Abweichung beträgt (absolut) 0,1554.

[Moritz G. (mbydq2)]

Arne schrieb:
> http://supp.iar.com/Support/?note=18180

Das ist das Heron-Verfahren aus Post:
Autor: Jobst
Datum: 25.04.2011 19:59
Y(n+1) = ((X/Yn)+Yn)/2

Bei geeigneter Wahl des ersten Schätzwertes, braucht man es nur ein mal 
anwenden um auf <<0,5% Abweichung zu kommen.

[joker (Gast)]

Moritz G. schrieb:
> Bei geeigneter Wahl des ersten Schätzwertes, braucht man es nur ein mal
> anwenden um auf <<0,5% Abweichung zu kommen.

Und wie lautet der erste Schätzwert?

[Moritz G. (mbydq2)]

joker schrieb:

> Und wie lautet der erste Schätzwert?

Ja, das ist natürlich das Problem.
Damit befasse ich mich ja in meinen letzten Beiträgen. Wie kann man ohne 
Division einen guten ersten Schätzwert finden? Da gibt es einige 
Möglichkeiten, ich habe alle mir bekannten genannt.
Die oben genannte Gerade kommt auf einem Intervall (a,b) b=5*a dem 
Ergebnis äußerst nahe.

Ich habe eine Verbesserung für die als C-Code veröffentlichte Methode 
gefunden:

   else{
     r>>=1;
     --n;
     r-=1<<n-1;//ginge auch durch XOR
     r>>=1;
     r+=1<<n-1;//ginge auch durch XOR
     r>>=n>>1;
     r*=3;//ließe sich in ASM schnell machen
     r>>=1
   }

Das nutzt jetzt:

$$\sqrt{2} \approx \frac{3}{2}$$

und geht daher auch bei kleineren Zahlen.

[Moritz G. (mbydq2)]

Hier noch weitere hilfreiche Links zum Thema:
Beitrag "Wurzel ziehen - VHDL"
Beitrag "[ASM] (schnelle) Integer Wurzel 32bit"
http://de.wikipedia.org/wiki/Schriftliches_Wurzelziehen
Man kommt also ganz ohne Näherungsverfahren aus.(fast schade, da es so 
viele Möglichkeiten gibt)
Schneller als O(n)=n geht es wohl nicht, wenn man es recht genau haben 
will.
Heron kommt wegen der Division eigentlich nicht in Betracht.

[joker (Gast)]

Hier eine fertige Routine

http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29

[Jobst M. (jobstens-de)]

Du könntest auch ein Verfahren mit sukzessiver Approximation anwenden.


Gruß

Jobst

[Hagen R. (hagen)]

Jobst M. schrieb:
> Du könntest auch ein Verfahren mit sukzessiver Approximation anwenden.
>
>
> Gruß
>
> Jobst

auch Newton Iteration genannt falls man auf Suche ist

[Jobst M. (jobstens-de)]

Hagen Re schrieb:
> Newton Iteration

Hmmm ... nö, so meinte ich das eigentlich nicht.


1. - Ergebnisregister auf 0 setzen
2. - erstes Bit im Ergebnisregister setzen.
     (Bzw. dort anfangen, wo es Sinn macht)
3. - Ergebnisregister quadrieren
4. - Ergebnis der Quadrierung mit dem Radikand vergleichen
5. - zu groß? -> Bit im Ergebnisregister wieder löschen
6. - genau gleich? -> Ergebnis gefunden
7. - nächstes Bit im Ergebnisregister setzen
8. - Solange nicht Ende erreicht bei 3 weiter machen.



Gruß

Jobst

[Moritz G. (mbydq2)]

Hallo Jobst,
das wäre was ich in 05.05.2011 08:19 meinte.
Auch das Heron oder Newton Verfahren funktionieren so ähnlich. Das Heron 
Verfahren kann so hergeleitet werden.

In Anlehnung an das zuvor von mir beschriebene Verfahren zum Erzeugen 
von Schätzwerten \ Startwerten hier eine Verbesserung:
Genutzt wird wie zuvor eine Zerlegung der Zahl unter der Wurzel in ein 
Produkt aus einer kleinen Zahl [0,35;1,4] und 4^n. Wegen sqrt(4^n)=2^n.
Danach müssen wir noch die Wurzel aus einer Zahl zwischen 0,35 und 1,4 
ziehen.
Das machen wir mit zwei Näherungen. Die erste Näherung ist die von Oben.
Die zweite Näherung für das Intervall 0,35 bis 0,98 ist:
15*(x/8-(x^2)/16)
Alternativ kann man Näherungen von 0,7 bis 0,98 und 1 bis 1,4 nehmen, 
wenn man zusätzlich eine 2 abspaltet. Dann muss man aber sqrt(2) ~= 3/2 
~= 181/128 = (53+128)/128 = 1+(53/128) benutzen und:
(13x-5x^2)/8 oder einfach eine Gerade.
Insgesamt kann man so einen sehr guten Startwert finden.

Das

if(n&&1){
    r-=1<<--n;//ginge auch durch XOR
    r>>=1;
    r+=1<<--n;//ginge auch durch XOR

müsste

if(n&1){
    r-=1<<(n-1);//ginge auch durch XOR
    r>>=1;
    r+=1<<(n-1);//ginge auch durch XOR

sein.
Das ist die Halbierung der Nachkommastellen, also die Näherung für 
(1;1,4].

[amateur (Gast)]

Wenn Du den originalen Quelltext zum Wurzelschlagen haben willst, musst 
Du Dir die Sourcen zur Mathematikbibliothek herunterladen. Ist ja 
offener Quellcode.
Im stillen Kämmerlein gäbe es noch die Möglichkeit irgendwie einen 
Funktionsaufruf in Dein Programm einzubauen und in der Ausgabedatei 
(.lst, .asm oder so) und dem Funktionsaufruf zu folgen.

[Anonymous U. (gastt)]

Hallo,

Ich hab mir einen code überlegt wie man Wurzeln aus Ganzzahlen ziehen 
kann. Was haltet ihr davon?

//Squareroot input=return; 4byte=2byte
unsigned short isqurt(unsigned long input){
  unsigned long rtemp;
  unsigned short lbound;
  unsigned short ubound;
  char maxgexp;
  unsigned short output=0;
  //estimation of bounds with max even exponent: squrt(2^2*n)=2^n
  rtemp = input;
  maxgexp = 0;
  while (rtemp){  //find index of highest exponent
    ++maxgexp;
    rtemp >>= 1;
  }
  --maxgexp;  //exponent = index-1
  if (maxgexp & 1)  //max even exponent
    --maxgexp;
  lbound = (1 << (maxgexp >> 1));  //lbound = 2^(maxgexp/2)
  ubound = (1 << ((maxgexp >> 1)+2));  //lbound = 2^((maxgexp/2)+2)
  //aproaching: calculate with value between boundaries and compare
  while (ubound - lbound){
    output = ((ubound - lbound) >> 1) + lbound;
    rtemp = output*output;
    if (rtemp <= input){
      if (lbound == output)
        break;
      lbound = output;
    }
    else
      ubound = output;
  }
  return output;
}

[c-hater (Gast)]

Hansen schrieb:

> Daher muss ich nochmal anders fragen - kann mir jemand rein aus
> Erfahrung sagen, ob man das Wurzelziehen auf einem uC OHNE FPU vermeiden
> sollte?

Was für eine schwachsinnige Frage.

Wenn man mathematisch die Möglichkeit hat, das Wurzelziehen zu 
vermeiden, dann wird man es tun. Wenn nicht, dann stellt sich die Frage 
garnicht, denn dann MUSS man es einfach tun.

Man könnte dann allenfalls noch über den Algorithmus diskutieren, mit 
dem man es tut. Das hängt aber von sehr konkreten Umständen ab, über die 
du in typischer Troll-Manier natürlich kein Wort verloren hast...

Nur als Hilfe zur Formulierung: hauptsächlich kommt es auf den 
Wertebereich und die geforderte Auflösung (bzw. den zulässigen Fehler) 
innerhalb des Wertebereichs an. Also genau so wie bei jeder anderen 
mathematischen Funktion auch, es sollte also einen echten Programmierer 
nicht wirklich überraschen, dass das wichtig sein könnte...

[Rufus Τ. F. (rufus)]

c-hater schrieb:
> Was für eine schwachsinnige Frage.

Um das festzustellen, hast Du jetzt geschlagene vier Jahre gebraucht?

[c-hater (Gast)]

Rufus Τ. Firefly schrieb:
> c-hater schrieb:
>> Was für eine schwachsinnige Frage.
>
> Um das festzustellen, hast Du jetzt geschlagene vier Jahre gebraucht?

Nein, das wußte ich sofort, nachdem ich die Frage gelesen hatte.

Dass die Frage vier Jahre alt war, war mir nicht bewußt. Ändert aber an 
der Korrektheit meiner Aussage eigentlich auch rein garnichts. Die ist 
heute so gültig wie vor vier Jahren, vor vierzig Jahren oder auch in 
hundert Jahren...

Das ist das Schöne an korrekten Aussagen, sie veralten einfach nicht...