sin() und arccos() : Assembler für die FPU

[helgew]

Hi,

besser wäre der Titel gewesen: "Optimierung von Funktionen: gekämpft und verloren." :wall:
Nunja, ihr kennt das, der Frust muss raus, und deshalb will ich euch das Werk des Wochenendes nicht unterschlagen, das hier für sin() etwa <10%, für arccos() immerhin ~30% Geschwindigkeitsgewinn bringt. Wenn sich jemand mit x86 Assembler und FPU-Befehlen vertraut machen will, könnte der code ja sogar noch zu etwas nütze sein. Für arccos ist der Fehler über das volle Intervall < 4ppm, für sin() ist er < 25ppm. Enjoy.

Delphi-Quellcode:

unit fmath;

interface

  function arccos_s6(x: single):single;
  function sin_6(x: single):single;
  function RndSingleParity(s: single): Cardinal;



implementation



var
  xs : single;

const
  a = +1.493415037;
  b = -0.159437656;
  c = +0.48560268e-1;
  d = -0.199230397e-1;
  e = +0.94214809e-2;
  f = -0.48413448e-2;
  g = +0.2627085e-2;

function arccos_s6(x: single):single; 
begin 
  if x > 0 then
    begin 
        if x > 1 then x := 1;
        xs :=   x - 0.42;
        result :=   sqrt(1-x) * ((((((g*xs+f)*xs+e)*xs+d)*xs+c)*xs+b)*xs + a);
    end 
  else
    begin 
        if x < -1 then x := -1;
        xs := - x - 0.42;
        result := - sqrt(1+x) * ((((((g*xs+f)*xs+e)*xs+d)*xs+c)*xs+b)*xs + a) + pi;
    end;
end;



var
  RNDSP_EXP : Cardinal;
  RNDSP_LOOP: integer;

function RndSingleParity(s: single): Cardinal; // zero: even, non-zero: odd
var
  c : Cardinal absolute s;
begin
  c := ((c + (1 shl 23)) and $3FFFFFFF) shl 1;
  RNDSP_EXP := c shr 24;
  c := (c  shl 7) ;{or (1 shl 31);} // naträglich auskommentiert
  
  if RNDSP_EXP < 25
    then for RNDSP_LOOP := RNDSP_EXP downto 1 do c := c shl 1
    else c := 0;
  
  result := c;
end;



var
  fcw16  : Word;
  fbuf32 : single;
  ibuf32 : cardinal absolute fbuf32;

const
  sin_6_ipi: single  = 1/pi;
  sin_6_a:single = 1;// sin(0.5*Pi);
  sin_6_c:single = -4.93429; // fitted values
  sin_6_e:single =  4.04506;
  sin_6_g:single = -1.23285;
  sin_6_p:single = 0.5; // argument offset for the polynomial
  sin_6_shift:single = 1E01; // must be even, higher orders of magnitude result in larger errors

function sin_6(x: single):single;
begin 
  asm
    // set FRNDINT mode to truncate
    fstcw fcw16;
    mov ax, fcw16 ;
    and ax, $f0ff ;             // Clears bits 8-11.
    or ax, $0c00 ;              // Rounding control=%11, Precision = %00.
    mov fcw16 , ax;
    fldcw fcw16;

    // perform
    fld [sin_6_p]               // keep 0.5 for later
    fld x;                      // |x|||
    fld [sin_6_ipi];            // |ipi|x|
    fmul;                       // |ipi*x||

    fld [sin_6_shift]
    fadd                        // add 10^2 to avoid the crack at the zero point. don't forget this anomaly,
                                //    limit input range to <10 pi , this shall do for most applications within -4pi..4pi
    fld st(0);                  // |ipi*x|ipi*x|
    frndint                     // |int(ipi*x)|ipi*x|
    fst fbuf32                  // store for subsequent parity analysis
    fsub                        // |ipi*x- int(ipi*x)||  = |modulo(ipi*x)||
    fabs                        // |abs(modulo(ipi*x)|| 0..1
    fsub                        // |abs(modulo(ipi*x)-0.5|| = |x||||||
    fld st(0)
    fmul                        // |x^2|||

    fld st(0)                   // |x^2|x^2||
    fmul st(1), st(0)           // |x^2|x^4||
    fld st(0)                   // |x^2|x^2|x^4||
    fmul st(0), st(2)           // |x^6|x^2|x^4||

    fld [sin_6_g]               // |g|x^6|x^2|x^4||
    fmulp st(1), st(0)          // |g*x^6|x^2|x^4||
    fld [sin_6_e]               // |e|g*x^6|x^2|x^4||
    fmulp st(3), st(0)          // |g*x^6|x^2|e*x^4||
    fld [sin_6_c]               // |c|g*x^6|x^2|e*x^4||
    fmulp st(2), st(0)          // |g*x^6|c*x^2|e*x^4||
    fld [sin_6_a]               // |a|g*x^6|c*x^2|e*x^4||
    fadd                        // |a + g*x^6|c*x^2|e*x^4||
    fadd                        // |a + g*x^6 + c*x^2|e*x^4||
    fadd                        // |a + g*x^6 + c*x^2 + e*x^4||

    fstp  [ESP]                 // return st(0) = a + c*x^2 + e*x^4 + g*x^6
  end;

  //      find int-part parity and correct result sign (32bit IEEE float only)
  //  
  ibuf32    := ((ibuf32 + (1 shl 23)) and $3FFFFFFF) shl 1;
  RNDSP_EXP := ibuf32 shr 24;
  ibuf32    := (ibuf32  shl 7) ;{or (1 shl 31);} // naträglich auskommentiert

  if RNDSP_EXP < 25
    then 
      begin 
          for RNDSP_LOOP := RNDSP_EXP downto 1 do 
            ibuf32 := ibuf32 shl 1;
      end
    else 
      begin
          ibuf32 := 0;
      end;

  //      now ibuf is 0 for even integer-parts of the normalized argument
  //      zero for positive half-waves
  //
  if LongBool(ibuf32) then result := -result;
end;




end.

ps. der richtige Weg wäre, den SSE-Befehlssatz auszuschöpfen und für Anwendungen wie Bildtransformationen (sphärische Projektionen), an denen ich gerade auch arbeite, parallel zu rechnen, indem man blockweise transformiert. Doch dazu vielleicht später einmal ein Nachtrag, wenn wieder frische Motivation gewachsen ist.

Greeds
helgew

Edit: etwas aufgeräumt...

[himitsu]

Auch wenn es jetzt nichts mit dem Thema zu tun hat,
aber immerhin geht es hier ja grob um "Optimierungen":

Deine Codeformatierung und vorallem die Einrückung ist grauenhaft.

Zitat:

Delphi-Quellcode:

begin
if x > 0 then
  begin if x > 1 then x := 1;
        xs :=   x - 0.42;
        result :=   sqrt(1-x) * ((((((g*xs+f)*xs+e)*xs+d)*xs+c)*xs+b)*xs + a);
  end else
  begin if x < -1 then x := -1;
        xs := - x - 0.42;
        result := - sqrt(1+x) * ((((((g*xs+f)*xs+e)*xs+d)*xs+c)*xs+b)*xs + a) + pi;
  end;
end;

Hier wollte ich dich erst darauf hinweisen, daß Result bei x <= 0 undefiniert wäre
(vorallem wegen des "übergeordneten" if x > 0 then )
und wollte dich gleichzeitig fragen, ob da überhaupt richtige Ergebnisse rauskommen könnten.

Denn rein optisch hatte ich erst diesen Aufbau im Kopf, denn genau darauf läßt deine Formatierung auf den ersten Blick schließen.

Delphi-Quellcode:

begin
  if x > 0 then
  begin 
    if x > 1 then
    begin
      x := 1;
      xs :=   x - 0.42;
      result :=   sqrt(1-x) * ((((((g*xs+f)*xs+e)*xs+d)*xs+c)*xs+b)*xs + a);
    end 
    else
      if x < -1 then
      begin
        x := -1;
        xs := - x - 0.42;
        result := - sqrt(1+x) * ((((((g*xs+f)*xs+e)*xs+d)*xs+c)*xs+b)*xs + a) + pi;
      end;
  end;
end;

[jfheins]

Mich würde das Verfahren dahinter interessieren.

Ich hätte eine Taylor-Serie vermutet, aber das sieht irgendwie komplizierter aus, als als die Auswertung eines Polynoms an einer Stelle.

[helgew]

Zitat:

Zitat von himitsu

Deine Codeformatierung und vorallem die Einrückung ist grauenhaft.

Da hast du leider recht und von dir höre ich mir das auch gerne an. Als ich das programmiert habe, war ich zudem mit Maple und gnuplot beschäftigt. Da sollte das Produkt nicht drunter leiden, tut es aber, zumindest bei mir und wenn man es aus Spaß macht. Aufräumen kommt dann später :-)

Da ich mir aber zumindest einen kleinen Lernwert bei den vorgestellten Prozeduren verspreche, sollte ich mein posting auch leichter erschließbar und nach allgemein anerkannten Gesichtspunkten hin akzeptabel formatieren. Ich habe die gröbsten Unklarheiten nun umformatiert. für die globalen Variablen, die nur in einer Prozedur verwendet werden und eigentlich unter das implementation - Schüsselwort gehören würden, fällt mir aber nichts saubereres ein.

Zitat:

Zitat von jfheins

Mich würde das Verfahren dahinter interessieren.

Taylorreihen vom Rand des Intervalls ab sind a-priori eine der schlechtesten Entwicklungsansätze. Ganz gleich, aus welchem Bedürfnis heraus man eine Funktion in eine Reihe entwickelt (sei es stetige Differenzierbarkeit der Approximation, Invertierbarkeit oder Kürzbarkeit mit anderen Termen) ist es immer nützlich, das Gültigkeitsintervall vorher zu kennen. Ein kleines Beispiel, wo es unablässlich ist:

In der komplexen Ebene wird eine Approximation innerhalb eines Quadrates 1 auf i benötigt, die Reihentwicklung der betrachteten Funktion ist aber nur gültig innerhalb eines Kreises mit Radius 1 (man spricht von Konvergenzradius). Legt man den Kreismittepunkt auf eine Seitenmitte, hat man schon ungültige Bereiche, legt man sie auf eine Ecke, sind nurnoch ~78% des benötigten Intervalls abgedeckt.

Ist der Fehler der Approximation in der Zahlenebene etwa ein Paraboloid, würde man den Kreismittelpunkt in die Mitte des Quadrats legen. Ist er asymmetrisch, so schiebt man den steileren Teil weiter zum Rand hin.


Genaugenommen habe den Fehler des Polynoms der arccos-Funktion gegen arccos(x) / sqrt(1-x) (der Hinweis kam von hier) über das Intervall [0;1] integriert und dieses Funktional minimiert. So kommt die 0.42 als Entwicklungspunkt zustande, wobei der Fehler am unteren Ende kritischer ist. Da auch die einseitige Ableitung der arccos-funktion an den Intervallgrenzen divergiert, nimmt man zunächst eine sqrt()-Funktion. Der Quotient ist dann ein sehr langsam veränderliches Polynom, fast eine Gerade.

Eine ebenfalls beliebte Vorgehensweise ist das fitting der Parameter. Für sin_6() habe ich mir 1000 Sinuswerte von 0..pi berechnet und diese um das Maximum bei pi/2 herum polynomial approximiert. Dort fängt man an mit etwas wie "1-x^2". Durch das gerade Verhalten fallen alle ungeraden Potenzen raus, das berücksichtigt man schon im Ansatz.

Eine dritte Methode, die hier keine Anwendung findet, die ich aber auch schon praktiziert habe, ist das Optimieren nach einzelnen Parametern hin, auf das Fehlerfunktional oder den dem Betrage nach größten, auftretenden Fehler (globales Extremum). So kann es schon den Fehler um eine Größenordnung verringern, aus einem Koeffizienten von 1 eine 0.94 zu machen und bei der Berechnung ist das natürlich vollkommen nebensächlich, welchen Wert man dann einsetzt. An dem Punkt ist man auch frei, beliebige Funktionen zur Approximation des Ausgangsproblems heranzuziehen, zum Beispiel Gausskurven für bestimmte, gebrochenrationale Funktionen oder eben diese für das Integral einer Gausskurve, welches man so nicht bestimmen kann.

Delphi-Quellcode:

function RndSingleParity(...)

Wenn ich schon beim Erklären bin, noch ein paar Worte zu der wohl kryptischsten Funktion hier:
32bit floats werden IEEE -konform dargestellt als

• 1 Signumsbit (MSB)
• 8 Exponentbits
• 23 bit Mantisse (bis LSB)

00: 0 00000000 00000000000000000000000
01: 0 01111111 00000000000000000000000 ( wenn ich mir das so ansehe, könnte man sich das einfügen der eins sogar sparen... )
02: 0 10000000 00000000000000000000000
03: 0 10000000 10000000000000000000000
04: 0 10000001 00000000000000000000000
05: 0 10000001 01000000000000000000000
06: 0 10000001 10000000000000000000000
07: 0 10000001 11000000000000000000000
08: 0 10000010 00000000000000000000000
09: 0 10000010 00100000000000000000000


Je nach Architektur/Betriebssystem (wer weiß dazu näheres?) kann es u.U. noch Unterschiede geben, ob little / big endian gespeichert wird.
Glücklicherweise ist nach FRNDINT der Single-Wert nach vorzeichenbefreitem Integer übersetzbar und die Parität (gerade Zahlen haben binär eine 0 in der "Einser"-Stelle). Wo man diese Null oder Eins finden kann, sagt der Exponent, nur ist dieser bei 0 : $00 und bei 1 : $7F, bei 2 dann $80 und danach aufsteigend. Das hat rührt von der Interpretation der Kommastelle her und dass das MSB des Oktetts als Exponentenvorzeichen interpretiert wird. Da es maximal 24 Stellen gibt, an denen das Paritätsbit stehen kann, habe ich den Exponenten auf 5 bit getrimmt und wegen der Problematik mit 0, +/-1 um eins erhöht. Da für eins an einer Stelle gelesen wird, die vorher nicht definiert war, muss nach Extrahieren des Exponenten eine 1 gesetzt werden. Das shifting ist zwar nicht gerade elegant und die Schleife wäre besser mit einer schnellen Sprungtabelle zu ersetzen, doch ich weiß noch nicht, wie.
shl / shr sieht weiters nur Konstanten vor, da müsste man schon zur Laufzeit im Maschinencode herumschreiben und das lasse ich an der Stelle.

[helgew]

Ich habe eben noch interessanten code gefunden, den man übersetzt auch verwenden könnte. (sin, cos, log, ...)