Schnelle Exp-Berechnung: Ist es möglich, die Genauigkeit zu verbessern, ohne zu viel Leistung zu verlieren?

Question

Ich probiere gerade die schnelle Funktion Exp(x) aus, die zuvor in diese Antwort auf eine SO-Frage zur Verbesserung der Rechengeschwindigkeit in C#:

public static double Exp(double x)
{
  var tmp = (long)(1512775 * x + 1072632447);
  return BitConverter.Int64BitsToDouble(tmp << 32);
}

Der Ausdruck verwendet einige IEEE-Gleitkomma-"Tricks" und ist in erster Linie für die Verwendung in neuronalen Gruppen gedacht. Die Funktion ist etwa 5 Mal schneller als die reguläre Math.Exp(x) Funktion.

Leider liegt die numerische Genauigkeit nur bei -4% bis +2% im Vergleich zur normalen Math.Exp(x) Funktion möchte ich im Idealfall eine Genauigkeit im unteren Prozentbereich erreichen.

Ich habe den Quotienten zwischen der approximativen und der regulären Exp-Funktion aufgetragen, und wie in der Grafik zu sehen ist, scheint sich der relative Unterschied mit praktisch konstanter Häufigkeit zu wiederholen.

Ist es möglich, diese Regelmäßigkeit zu nutzen, um die Genauigkeit der Funktion "fast exp" weiter zu verbessern, ohne die Berechnungsgeschwindigkeit wesentlich zu verringern, oder würde der Rechenaufwand einer Genauigkeitsverbesserung den rechnerischen Gewinn des ursprünglichen Ausdrucks überwiegen?

(Nebenbei bemerkt, habe ich auch versucht eine der Alternativen Ansätze, die in der gleichen SO-Frage vorgeschlagen wurden, aber dieser Ansatz scheint in C# nicht rechnerisch effizient zu sein, zumindest nicht für den allgemeinen Fall).

UPDATE 14. MAI

Auf Anfrage von @Adriano habe ich nun einen sehr einfachen Benchmark durchgeführt. Ich habe 10 Millionen Berechnungen mit jeder der folgenden Alternativen durchgeführt exp Funktionen für Fließkommazahlen im Bereich [-100, 100]. Da der Wertebereich, an dem ich interessiert bin, von -20 bis 0 reicht, habe ich auch den Funktionswert bei x = -5 explizit aufgeführt. Hier sind die Ergebnisse:

      Math.Exp: 62.525 ms, exp(-5) = 0.00673794699908547
Empty function: 13.769 ms
     ExpNeural: 14.867 ms, exp(-5) = 0.00675211846828461
    ExpSeries8: 15.121 ms, exp(-5) = 0.00641270968867667
   ExpSeries16: 32.046 ms, exp(-5) = 0.00673666189488182
          exp1: 15.062 ms, exp(-5) = -12.3333325982094
          exp2: 15.090 ms, exp(-5) = 13.708332516253
          exp3: 16.251 ms, exp(-5) = -12.3333325982094
          exp4: 17.924 ms, exp(-5) = 728.368055056781
          exp5: 20.972 ms, exp(-5) = -6.13293614238501
          exp6: 24.212 ms, exp(-5) = 3.55518353166184
          exp7: 29.092 ms, exp(-5) = -1.8271053775984
      exp7 +/-: 38.482 ms, exp(-5) = 0.00695945286970704

ExpNeural ist gleichbedeutend mit dem Exp Funktion, die am Anfang dieses Textes beschrieben wird. ExpSerie8 es el Formulierung die ich ursprünglich als nicht sehr effizient unter .NET bezeichnet hatte; bei der Implementierung genau wie Neil war sie tatsächlich sehr schnell. ExpSerie16 ist die analoge Formel, aber mit 16 Multiplikationen anstelle von 8. exp1 über exp7 sind die verschiedenen Funktionen aus Adriano's Antwort unten. Die letzte Variante von exp7 ist eine Variante, bei der das Vorzeichen von x wird geprüft; bei negativem Ergebnis gibt die Funktion 1/exp(-x) stattdessen.

Leider ist keiner der beiden expN Funktionen, die von Adriano sind in dem von mir in Betracht gezogenen breiteren negativen Wertebereich ausreichend. Der Ansatz der Reihenentwicklung von Neil Coffey scheint in "meinem" Wertebereich besser geeignet zu sein, obwohl er bei größeren negativen Werten zu schnell abweicht x , insbesondere wenn "nur" 8 Multiplikationen verwendet werden.

Answer 1

Ich habe für meine Zwecke die folgende Funktion entwickelt, die schnell und genau den natürlichen Exponenten mit einfacher Genauigkeit berechnet. Die Funktion funktioniert im gesamten Bereich der Float-Werte. Der Code ist unter Visual Studio (x86) geschrieben.

_declspec(naked) float _vectorcall fexp(float x)
{
  static const float ct[8] =       // Constants table
  {
    1.44269502f,                   // lb(e)
    1.92596299E-8f,                // Correction to the value lb(e)
    -9.21120925E-4f,               // 16*b2
    0.115524396f,                  // 4*b1
    2.88539004f,                   // b0
    4.29496730E9f,                 // 2^32
    0.5f,                          // 0.5
    2.32830644E-10f                // 2^-32
  };
  _asm
  {
    mov ecx,offset ct              // ecx contains the address of constants tables
    vmulss xmm1,xmm0,[ecx]         // xmm1 = x*lb(e)
    vcvtss2si eax,xmm1             // eax = round(x*lb(e)) = k
    vcvtsi2ss xmm1,xmm1,eax        // xmm1 = k
    dec eax                        // of=1, if eax=80000000h, i.e. overflow
    jno exp_cont                   // Jump to form 2^k, if k is normal
    vucomiss xmm0,xmm0             // Compare xmm0 with itself to identify NaN
    jp exp_end                     // Complete with NaN result, if x=NaN
    vmovd eax,xmm0                 // eax contains the bits of x
    test eax,eax                   // sf=1, if x<0; of=0 always
      exp_break:                   // Get 0 for x<0 or Inf for x>0
    vxorps xmm0,xmm0,xmm0          // xmm0 = 0
    jle exp_end                    // Ready, if x<<0
    vrcpss xmm0,xmm0,xmm0          // xmm0 = Inf at case x>>0
      exp_end:                     // The result at xmm0 is ready
    ret                            // Return
      exp_cont:                    // Continue if |x| is not too big 
    vfmsub231ss xmm1,xmm0,[ecx]    // xmm1 = x*lb(e)-k = t/2 in the range from -0,5 to 0,5
    cdq                            // edx=-1, if x<0, othervice edx=0
    vfmadd132ss xmm0,xmm1,[ecx+4]  // xmm0 = t/2 (corrected value)
    and edx,8                      // edx=8, if x<0, othervice edx=0
    vmulss xmm2,xmm0,xmm0          // xmm2 = t^2/4 - the argument of polynomial
    vmovss xmm1,[ecx+8]            // Initialize the sum with highest coefficient 16*b2
    lea eax,[eax+8*edx+97]         // The exponent of 2^(k-31), if x>0, othervice 2^(k+33)
    vfmadd213ss xmm1,xmm2,[ecx+12] // xmm1 = 4*b1+4*b2*t^2
    test eax,eax                   // Test the sign of x
    jle exp_break                  // Result is 0 if the exponent is too small
    vfmsub213ss xmm1,xmm2,xmm0     // xmm1 = -t/2+b1*t^2+b2*t^4
    cmp eax,254                    // Check that the exponent is not too large
    ja exp_break                   // Jump to set Inf if overflow
    vaddss xmm1,xmm1,[ecx+16]      // xmm1 = b0-t/2+b1*t^2+b2*t^4 = f(t)-t/2
    shl eax,23                     // eax contains the bits of 2^(k-31) or 2^(k+33)
    vdivss xmm0,xmm0,xmm1          // xmm0 = t/(2*f(t)-t)
    vmovd xmm2,eax                 // xmm2 = 2^(k-31), if x>0; otherwice 2^(k+33)
    vaddss xmm0,xmm0,[ecx+24]      // xmm0 = t/(2*f(t)-t)+0,5
    vmulss xmm0,xmm0,xmm2          // xmm0 = e^x with shifted exponent of +-32
    vmulss xmm0,xmm0,[ecx+edx+20]  // xmm0 = e^x with corrected exponent
    ret                            // Return
  }
}