Thema: 3te wurzel berechnen

[TomCat]

Hallo,

für die Quadratwurzel gibt es ja den FPU -Befehl: fsqrt.

Wie kann man aber z.b. die 3. Wurzel (Kubikwurzel) berechnen?


THX
TomCat

[Svenska]

Es gilt:
"nte Wurzel aus x" = "x hoch (1/n)" = "e hoch ( ln(x) / n )".

Du kannst also jede Wurzel (wie auch jede Potenz) auf die Exponential- und Logarithmusfunktion zurückführen.

Quelle: Wikipedia.

Nachtrag: Du kannst natürlich auch Näherungsverfahren anwenden.

Vermutlich ist es aber real sinnvoller, diese dann für Exponential- und Logarithmusfunktion zu benutzen, was im Übrigen vor allem sinnvoll ist, wenn die Hardware selbst keine FPU enthält oder die Genauigkeit davon nicht ausreicht (die meisten Mikrocontroller haben keine FPU, und die meisten mit FPU können nur float aber kein double) - oder aber die FPU nicht benutzt werden kann (z.B. im Kernel).

[TomCat]

Danke erstmal für die umfangreiche Antwort.

ich habe vergessen zu erwähnen, dass es um x86 ASM Programmierung geht.

mir geht's hier nicht Geschwindigkeit, sondern darum dies mit der FPU zu lösen.
bzw. wie macht dies ein Compiler?

 

[Svenska]

Hallo,

weder die 8087-FPU noch die C-Programmiersprache kennen einen Kubikwurzel-Befehl. Um die Kubikwurzel zu berechnen, musst du also die Funktion irgendwie nachbauen. Dazu brauchst du, wenn du die von mir genannte Identität benutzen möchtest, eine e-Funktion und einen natürlichen Logarithmus. Beide Funktionen sind im x87 nicht direkt vorhanden, lassen sich aber mit den gegebenen Methoden nachbauen.

Um rauszukriegen, wie ein Compiler das macht, kannst du ihn das einfach machen lassen:

Code: [Auswählen]

$ cat test.c
#include <stdio.h>
#include <math.h>

double kw(double x)
{
        return exp(log(x)/3.0);
}

int main()
{
        printf("Kubikwurzel: %.2f\n", kw(8));
}

Dann einfach durch den Compiler jagen und disassemblieren:

Code: [Auswählen]

$ i686-w64-mingw32-gcc -O3 -march=i486 -o test test.c
$ i686-w64-mingw32-objdump -d test.exe
[...]
00401510 <_kw>:
  401510:       83 ec 1c                sub    $0x1c,%esp
  401513:       dd 44 24 20             fldl   0x20(%esp)
  401517:       dd 1c 24                fstpl  (%esp)
  40151a:       e8 81 11 00 00          call   4026a0 <_log>
  40151f:       d8 35 14 40 40 00       fdivs  0x404014
  401525:       dd 5c 24 20             fstpl  0x20(%esp)
  401529:       83 c4 1c                add    $0x1c,%esp
  40152c:       e9 9f 0f 00 00          jmp    4024d0 <_exp>
  401531:       90                      nop
[...]
004026a0 <_log>:
  4026a0:       83 ec 4c                sub    $0x4c,%esp
  4026a3:       dd 44 24 50             fldl   0x50(%esp)
  4026a7:       d9 e5                   fxam   
  4026a9:       9b df e0                fstsw  %ax
  4026ac:       66 25 00 45             and    $0x4500,%ax
  4026b0:       66 3d 00 40             cmp    $0x4000,%ax
  4026b4:       0f 84 a6 00 00 00       je     402760 <_log+0xc0>
  4026ba:       89 c2                   mov    %eax,%edx
  4026bc:       d9 e5                   fxam   
  4026be:       9b df e0                fstsw  %ax
  4026c1:       f6 c4 02                test   $0x2,%ah
  4026c4:       75 2a                   jne    4026f0 <_log+0x50>
  4026c6:       66 81 fa 00 05          cmp    $0x500,%dx
  4026cb:       74 73                   je     402740 <_log+0xa0>
  4026cd:       66 81 fa 00 01          cmp    $0x100,%dx
  4026d2:       74 7c                   je     402750 <_log+0xb0>
  4026d4:       db 3c 24                fstpt  (%esp)
  4026d7:       e8 e4 00 00 00          call   4027c0 <___logl_internal>
  4026dc:       dd 5c 24 38             fstpl  0x38(%esp)
  4026e0:       dd 44 24 38             fldl   0x38(%esp)
  4026e4:       83 c4 4c                add    $0x4c,%esp
  4026e7:       c3                      ret   
  4026e8:       8d b4 26 00 00 00 00    lea    0x0(%esi,%eiz,1),%esi
  4026ef:       90                      nop
  4026f0:       dd 5c 24 30             fstpl  0x30(%esp)
  4026f4:       e8 6b 01 00 00          call   402864 <__errno>
  4026f9:       d9 05 90 42 40 00       flds   0x404290
  4026ff:       c7 00 21 00 00 00       movl   $0x21,(%eax)
  402705:       c7 44 24 04 88 42 40    movl   $0x404288,0x4(%esp)
  40270c:       00
  40270d:       dd 54 24 18             fstl   0x18(%esp)
  402711:       c7 04 24 01 00 00 00    movl   $0x1,(%esp)
  402718:       d9 5c 24 2c             fstps  0x2c(%esp)
  40271c:       d9 ee                   fldz   
  40271e:       dd 5c 24 10             fstpl  0x10(%esp)
  402722:       dd 44 24 30             fldl   0x30(%esp)
  402726:       dd 5c 24 08             fstpl  0x8(%esp)
  40272a:       e8 11 f1 ff ff          call   401840 <___mingw_raise_matherr>
  40272f:       d9 44 24 2c             flds   0x2c(%esp)
  402733:       83 c4 4c                add    $0x4c,%esp
  402736:       c3                      ret   
  402737:       8d b4 26 00 00 00 00    lea    0x0(%esi,%eiz,1),%esi
  40273e:       66 90                   xchg   %ax,%ax
  402740:       dd d8                   fstp   %st(0)
  402742:       d9 05 94 42 40 00       flds   0x404294
  402748:       83 c4 4c                add    $0x4c,%esp
  40274b:       c3                      ret   
  40274c:       8d 74 26 00             lea    0x0(%esi,%eiz,1),%esi
  402750:       dd d8                   fstp   %st(0)
  402752:       d9 05 90 42 40 00       flds   0x404290
  402758:       83 c4 4c                add    $0x4c,%esp
  40275b:       c3                      ret   
  40275c:       8d 74 26 00             lea    0x0(%esi,%eiz,1),%esi
  402760:       dd 5c 24 30             fstpl  0x30(%esp)
  402764:       e8 fb 00 00 00          call   402864 <__errno>
  402769:       d9 05 8c 42 40 00       flds   0x40428c
  40276f:       c7 00 22 00 00 00       movl   $0x22,(%eax)
  402775:       c7 44 24 04 88 42 40    movl   $0x404288,0x4(%esp)
  40277c:       00
  40277d:       dd 54 24 18             fstl   0x18(%esp)
  402781:       c7 04 24 03 00 00 00    movl   $0x3,(%esp)
  402788:       d9 5c 24 2c             fstps  0x2c(%esp)
  40278c:       d9 ee                   fldz   
  40278e:       dd 5c 24 10             fstpl  0x10(%esp)
  402792:       dd 44 24 30             fldl   0x30(%esp)
  402796:       dd 5c 24 08             fstpl  0x8(%esp)
  40279a:       e8 a1 f0 ff ff          call   401840 <___mingw_raise_matherr>
  40279f:       d9 44 24 2c             flds   0x2c(%esp)
  4027a3:       83 c4 4c                add    $0x4c,%esp
  4027a6:       c3                      ret   
  4027a7:       90                      nop
[...]

Ganz allgemein gilt: Wenn eine Funktion brauchst, die es in Hardware oder in der Programmiersprache nicht gibt, dann musst du dir einen Algorithmus (d.h. einen Satz an nacheinander auszuführenden Regeln) ausdenken, der diese Funktion mit den gegebenen Möglichkeiten nachbaut. Im Falle deiner Kubikwurzel sind das also Exponential- und Logarithmusfunktion. Das habe ich so hingeschrieben, und das hat der Compiler so gemacht. Dazu hat er entsprechende Bibliotheksfunktionen (die in libc/libm implementiert sind) benutzt.

Wie du an meinem Beispiel siehst, ist das nicht unbedingt trivial. Unter bestimmten Randbedingungen lässt sich das sicherlich noch deutlich weiter optimieren. Die FPU kann einen Logarithmus berechnen (FY2LX), die E-Funktion muss man irgendwie nachbauen (Stackoverflow findet eine Sequenz aus 9 Befehlen). Da ich mich mit der FPU nicht besonders auskenne, kann ich nur darauf verweisen. :-)

Gag am Rande: Wenn ich die Funktion "static" mache, so dass die Funktion nur innerhalb der Datei existiert, dann optimiert der Compiler die gesamte Berechnung raus, weil er kw(8) in dem Beispiel schon vorher ausrechnen kann - das Ergebnis ist konstant.

Gruß,
Svenska