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

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

... 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 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, ...?)

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

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

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?

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.

>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

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.

Für

geht

und dann das übliche

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

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

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:

1

//2.Wurzelziehen

2

int sqrt(int r){

3

//führende nullen zählen; nur für 16bit=lengthof(r)

4

   char n=0;

5

   int t=r;

6

   if(r > 0) //oberste Stelle ist frei

7

    for( ; t > 0; ++n) t <<= 1;

8

    n=16-n;

9

  t=r;

10

    if(n&&1){

11

   r-=1<<--n;//ginge auch durch XOR

12

   r>>=1;

13

   r+=1<<--n;//ginge auch durch XOR

14

   r>>=n>>1;

15

   }

16

  else{

17

   r>>=3;//(2 hoch 3)~~(3 hoch 2)

18

   n-=3;

19

   r-=1<<--n;//ginge auch durch XOR

20

   r>>=1;

21

   r+=1<<--n;//ginge auch durch XOR

22

   r>>=n>>1;

23

   r*=3;

24

   }

25

  r+=r+t/r>>1;//einzige Division; sollte unter 1% Abweichung haben.

26

  return(r);

machen.
Eine Division ist immernoch nötig.

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.

1

int sqrt(int r){

2

//führende nullen zählen; nur für 16bit=lengthof(r)

3

   char n=0;

4

   int t=r;

5

   if(r > 0) for( ; t > 0; ++n) t <<= 1;//oberste Stelle ist frei

6

   else return(0);//Fehler 1/2

7

   n=16-n;

8

   t=r;

9

   if(n&&1){

10

    r-=1<<--n;//ginge auch durch XOR

11

    r>>=1;

12

    r+=1<<--n;//ginge auch durch XOR

13

    r>>=n>>1;

14

   }

15

   else{

16

    r>>=3;//(2 hoch 3)~~(3 hoch 2)

17

    n-=3;

18

    r-=1<<n-1;//ginge auch durch XOR - Fehler 2/2

19

    r>>=1;

20

    r+=1<<n-1;//ginge auch durch XOR - Fehler 2/2

21

    r>>=n>>1;

22

    r*=3;

23

   }

24

   r+=r+t/r>>1;//einzige Division; sollte unter 1% Abweichung haben.

25

   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.

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)

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

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

1

   else{

2

     r>>=1;

3

     --n;

4

     r-=1<<n-1;//ginge auch durch XOR

5

     r>>=1;

6

     r+=1<<n-1;//ginge auch durch XOR

7

     r>>=n>>1;

8

     r*=3;//ließe sich in ASM schnell machen

9

     r>>=1

10

   }

Das nutzt jetzt:

und geht daher auch bei kleineren Zahlen.

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.

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


Gruß

Jobst

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

auch Newton Iteration genannt falls man auf Suche ist

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

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

1

if(n&&1){

2

    r-=1<<--n;//ginge auch durch XOR

3

    r>>=1;

4

    r+=1<<--n;//ginge auch durch XOR

müsste

1

if(n&1){

2

    r-=1<<(n-1);//ginge auch durch XOR

3

    r>>=1;

4

    r+=1<<(n-1);//ginge auch durch XOR

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

Hallo,

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

1

//Squareroot input=return; 4byte=2byte

2

unsigned short isqurt(unsigned long input){

3

  unsigned long rtemp;

4

  unsigned short lbound;

5

  unsigned short ubound;

6

  char maxgexp;

7

  unsigned short output=0;

8

  //estimation of bounds with max even exponent: squrt(2^2*n)=2^n

9

  rtemp = input;

10

  maxgexp = 0;

11

  while (rtemp){  //find index of highest exponent

12

    ++maxgexp;

13

    rtemp >>= 1;

14

  }

15

  --maxgexp;  //exponent = index-1

16

  if (maxgexp & 1)  //max even exponent

17

    --maxgexp;

18

  lbound = (1 << (maxgexp >> 1));  //lbound = 2^(maxgexp/2)

19

  ubound = (1 << ((maxgexp >> 1)+2));  //lbound = 2^((maxgexp/2)+2)

20

  //aproaching: calculate with value between boundaries and compare

21

  while (ubound - lbound){

22

    output = ((ubound - lbound) >> 1) + lbound;

23

    rtemp = output*output;

24

    if (rtemp <= input){

25

      if (lbound == output)

26

        break;

27

      lbound = output;

28

    }

29

    else

30

      ubound = output;

31

  }

32

  return output;

33

}

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

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