wake-up-neo.net

Optimierungen für pow () mit const nicht ganzzahligem Exponenten?

Ich habe Hotspots in meinem Code, an denen ich pow() mache und die ungefähr 10-20% meiner Ausführungszeit beanspruchen.

Meine Eingabe für pow(x,y) ist sehr spezifisch, daher frage ich mich, ob es eine Möglichkeit gibt, zwei pow()-Approximationen (eine für jeden Exponenten) mit höherer Leistung zu würfeln:

  • Ich habe zwei konstante Exponenten: 2.4 und 1/2.4.
  • Wenn der Exponent 2,4 ist, liegt x im Bereich (0,090473935, 1,0].
  • Wenn der Exponent 1/2,4 ist, liegt x im Bereich (0,0031308, 1,0].
  • Ich verwende SSE/AVX float Vektoren. Wenn die Plattformspezifikationen genutzt werden können, klicken Sie auf!

Eine maximale Fehlerrate um 0,01% ist ideal, obwohl ich mich auch für Algorithmen mit voller Genauigkeit (für float) interessiere.

Ich verwende bereits eine schnelle pow()Approximation , aber diese Einschränkungen werden nicht berücksichtigt. Kann man es besser machen?

57
Cory Nelson

In der IEEE 754-Hacking-Ader gibt es eine andere Lösung, die schneller und weniger "magisch" ist. Sie erreicht eine Fehlerspanne von 0,08% in etwa einem Dutzend Taktzyklen (für den Fall von p = 2,4 bei einer Intel Merom-CPU).

Fließkommazahlen wurden ursprünglich als Annäherung an Logarithmen entwickelt. Sie können also den ganzzahligen Wert als Annäherung an log2 verwenden. Dies ist etwas portabel durch Anwenden des Befehls convert-from-integer auf einen Gleitkommawert möglich, um einen anderen Gleitkommawert zu erhalten.

Um die pow-Berechnung abzuschließen, können Sie mit einem konstanten Faktor multiplizieren und den Logarithmus mit der Anweisung convert-to-integer zurück konvertieren. Bei SSE lauten die entsprechenden Anweisungen cvtdq2ps und cvtps2dq.

Es ist jedoch nicht ganz so einfach. Das Exponentenfeld in IEEE 754 ist vorzeichenbehaftet, wobei der Vorspannungswert 127 einen Exponenten von Null darstellt. Diese Verzerrung muss entfernt werden, bevor Sie den Logarithmus multiplizieren, und vor der Potenzierung erneut hinzugefügt. Darüber hinaus funktioniert die Bias-Anpassung durch Subtraktion nicht bei Null. Glücklicherweise können beide Einstellungen durch vorheriges Multiplizieren mit einem konstanten Faktor erreicht werden.

x^p
= exp2( p * log2( x ) )
= exp2( p * ( log2( x ) + 127 - 127 ) - 127 + 127 )
= cvtps2dq( p * ( log2( x ) + 127 - 127 - 127 / p ) )
= cvtps2dq( p * ( log2( x ) + 127 - log2( exp2( 127 - 127 / p ) ) )
= cvtps2dq( p * ( log2( x * exp2( 127 / p - 127 ) ) + 127 ) )
= cvtps2dq( p * ( cvtdq2ps( x * exp2( 127 / p - 127 ) ) ) )

exp2( 127 / p - 127 ) ist der konstante Faktor. Diese Funktion ist eher spezialisiert: Sie funktioniert nicht mit kleinen fraktionalen Exponenten, da der konstante Faktor exponentiell mit dem Inversen des Exponenten wächst und überläuft. Bei negativen Exponenten funktioniert es nicht. Große Exponenten führen zu hohen Fehlern, da die Mantissenbits durch die Multiplikation mit den Exponentenbits vermischt werden.

Aber es sind nur 4 schnelle Anweisungen. Vor dem Multiplizieren von "Ganzzahl" (in den Logarithmus) konvertieren, Power-Multiplizieren, in "Ganzzahl" (vom Logarithmus) konvertieren. Konvertierungen sind bei dieser Implementierung von SSE sehr schnell. Wir können auch einen zusätzlichen konstanten Koeffizienten in die erste Multiplikation einquetschen.

template< unsigned expnum, unsigned expden, unsigned coeffnum, unsigned coeffden >
__m128 fastpow( __m128 arg ) {
        __m128 ret = arg;
//      std::printf( "arg = %,vg\n", ret );
        // Apply a constant pre-correction factor.
        ret = _mm_mul_ps( ret, _mm_set1_ps( exp2( 127. * expden / expnum - 127. )
                * pow( 1. * coeffnum / coeffden, 1. * expden / expnum ) ) );
//      std::printf( "scaled = %,vg\n", ret );
        // Reinterpret arg as integer to obtain logarithm.
        asm ( "cvtdq2ps %1, %0" : "=x" (ret) : "x" (ret) );
//      std::printf( "log = %,vg\n", ret );
        // Multiply logarithm by power.
        ret = _mm_mul_ps( ret, _mm_set1_ps( 1. * expnum / expden ) );
//      std::printf( "powered = %,vg\n", ret );
        // Convert back to "integer" to exponentiate.
        asm ( "cvtps2dq %1, %0" : "=x" (ret) : "x" (ret) );
//      std::printf( "result = %,vg\n", ret );
        return ret;
}

Einige Versuche mit Exponent = 2,4 zeigen, dass dies ständig um etwa 5% überschätzt wird. (Die Routine wird garantiert immer überschätzt.) Sie können einfach mit 0,95 multiplizieren, aber mit ein paar weiteren Anweisungen erhalten Sie eine Genauigkeit von 4 Dezimalstellen, was für Grafiken ausreichend sein sollte.

Der Schlüssel ist, das Überschätzte mit einem Unterschätzten abzugleichen und den Durchschnitt zu ermitteln.

  • Berechne x ^ 0,8: vier Anweisungen, Fehler ~ + 3%.
  • Berechne x ^ -0.4: eine rsqrtps. (Dies ist ziemlich genau genug, opfert aber die Fähigkeit, mit Null zu arbeiten.)
  • Berechne x ^ 0,4: eine mulps.
  • Berechne x ^ -0.2: eine rsqrtps.
  • Berechne x ^ 2: one mulps.
  • Berechne x ^ 3: one mulps.
  • x ^ 2,4 = x ^ 2 * x ^ 0,4: eine mulps. Dies ist das Überschätzen.
  • x ^ 2,4 = x ^ 3 * x ^ -0,4 * x ^ -0,2: zwei mulps. Dies ist das unterschätzte.
  • Den Durchschnitt der obigen Werte berechnen: eine addps, eine mulps.

Anweisungszähler: vierzehn, einschließlich zwei Konvertierungen mit Latenz = 5 und zwei reziproken Quadratwurzelschätzungen mit Durchsatz = 4.

Um den Durchschnitt richtig zu ermitteln, wollen wir die Schätzungen mit ihren erwarteten Fehlern gewichten. Die Unterschätzung erhöht den Fehler auf eine Potenz von 0,6 gegen 0,4, so dass wir erwarten, dass er 1,5-fach als fehlerhaft gilt. Die Gewichtung fügt keine Anweisungen hinzu. es kann im Vorfaktor erfolgen. Man nenne den Koeffizienten a: a ^ 0,5 = 1,5 a ^ -0,75 und a = 1,38316186.

Der endgültige Fehler ist ungefähr 0,015% oder 2 Größenordnungen besser als das erste fastpow-Ergebnis. Die Laufzeit beträgt etwa ein Dutzend Zyklen für eine ausgelastete Schleife mit volatile-Quell- und Zielvariablen. Obwohl sie die Iterationen überlappt, wird bei der Verwendung in der realen Welt auch Parallelität auf Befehlsebene zu sehen sein. Bei SIMD ist dies ein Durchsatz von einem Skalarergebnis pro 3 Zyklen!

int main() {
        __m128 const x0 = _mm_set_ps( 0.01, 1, 5, 1234.567 );
        std::printf( "Input: %,vg\n", x0 );

        // Approx 5% accuracy from one call. Always an overestimate.
        __m128 x1 = fastpow< 24, 10, 1, 1 >( x0 );
        std::printf( "Direct x^2.4: %,vg\n", x1 );

        // Lower exponents provide lower initial error, but too low causes overflow.
        __m128 xf = fastpow< 8, 10, int( 1.38316186 * 1e9 ), int( 1e9 ) >( x0 );
        std::printf( "1.38 x^0.8: %,vg\n", xf );

        // Imprecise 4-cycle sqrt is still far better than fastpow, good enough.
        __m128 xfm4 = _mm_rsqrt_ps( xf );
        __m128 xf4 = _mm_mul_ps( xf, xfm4 );

        // Precisely calculate x^2 and x^3
        __m128 x2 = _mm_mul_ps( x0, x0 );
        __m128 x3 = _mm_mul_ps( x2, x0 );

        // Overestimate of x^2 * x^0.4
        x2 = _mm_mul_ps( x2, xf4 );

        // Get x^-0.2 from x^0.4. Combine with x^-0.4 into x^-0.6 and x^2.4.
        __m128 xfm2 = _mm_rsqrt_ps( xf4 );
        x3 = _mm_mul_ps( x3, xfm4 );
        x3 = _mm_mul_ps( x3, xfm2 );

        std::printf( "x^2 * x^0.4: %,vg\n", x2 );
        std::printf( "x^3 / x^0.6: %,vg\n", x3 );
        x2 = _mm_mul_ps( _mm_add_ps( x2, x3 ), _mm_set1_ps( 1/ 1.960131704207789 ) );
        // Final accuracy about 0.015%, 200x better than x^0.8 calculation.
        std::printf( "average = %,vg\n", x2 );
}

Tja… Entschuldigung, dass ich das früher nicht posten konnte. Und die Erweiterung auf x ^ 1/2.4 bleibt als Übung übrig; v).


Update mit Statistiken

Ich habe ein kleines Testgeschirr und zwei x implementiert(512) Fälle, die den oben genannten entsprechen.

#include <cstdio>
#include <xmmintrin.h>
#include <cmath>
#include <cfloat>
#include <algorithm>
using namespace std;

template< unsigned expnum, unsigned expden, unsigned coeffnum, unsigned coeffden >
__m128 fastpow( __m128 arg ) {
    __m128 ret = arg;
//  std::printf( "arg = %,vg\n", ret );
    // Apply a constant pre-correction factor.
    ret = _mm_mul_ps( ret, _mm_set1_ps( exp2( 127. * expden / expnum - 127. )
        * pow( 1. * coeffnum / coeffden, 1. * expden / expnum ) ) );
//  std::printf( "scaled = %,vg\n", ret );
    // Reinterpret arg as integer to obtain logarithm.
    asm ( "cvtdq2ps %1, %0" : "=x" (ret) : "x" (ret) );
//  std::printf( "log = %,vg\n", ret );
    // Multiply logarithm by power.
    ret = _mm_mul_ps( ret, _mm_set1_ps( 1. * expnum / expden ) );
//  std::printf( "powered = %,vg\n", ret );
    // Convert back to "integer" to exponentiate.
    asm ( "cvtps2dq %1, %0" : "=x" (ret) : "x" (ret) );
//  std::printf( "result = %,vg\n", ret );
    return ret;
}

__m128 pow125_4( __m128 arg ) {
    // Lower exponents provide lower initial error, but too low causes overflow.
    __m128 xf = fastpow< 4, 5, int( 1.38316186 * 1e9 ), int( 1e9 ) >( arg );

    // Imprecise 4-cycle sqrt is still far better than fastpow, good enough.
    __m128 xfm4 = _mm_rsqrt_ps( xf );
    __m128 xf4 = _mm_mul_ps( xf, xfm4 );

    // Precisely calculate x^2 and x^3
    __m128 x2 = _mm_mul_ps( arg, arg );
    __m128 x3 = _mm_mul_ps( x2, arg );

    // Overestimate of x^2 * x^0.4
    x2 = _mm_mul_ps( x2, xf4 );

    // Get x^-0.2 from x^0.4, and square it for x^-0.4. Combine into x^-0.6.
    __m128 xfm2 = _mm_rsqrt_ps( xf4 );
    x3 = _mm_mul_ps( x3, xfm4 );
    x3 = _mm_mul_ps( x3, xfm2 );

    return _mm_mul_ps( _mm_add_ps( x2, x3 ), _mm_set1_ps( 1/ 1.960131704207789 * 0.9999 ) );
}

__m128 pow512_2( __m128 arg ) {
    // 5/12 is too small, so compute the sqrt of 10/12 instead.
    __m128 x = fastpow< 5, 6, int( 0.992245 * 1e9 ), int( 1e9 ) >( arg );
    return _mm_mul_ps( _mm_rsqrt_ps( x ), x );
}

__m128 pow512_4( __m128 arg ) {
    // 5/12 is too small, so compute the 4th root of 20/12 instead.
    // 20/12 = 5/3 = 1 + 2/3 = 2 - 1/3. 2/3 is a suitable argument for fastpow.
    // weighting coefficient: a^-1/2 = 2 a; a = 2^-2/3
    __m128 xf = fastpow< 2, 3, int( 0.629960524947437 * 1e9 ), int( 1e9 ) >( arg );
    __m128 xover = _mm_mul_ps( arg, xf );

    __m128 xfm1 = _mm_rsqrt_ps( xf );
    __m128 x2 = _mm_mul_ps( arg, arg );
    __m128 xunder = _mm_mul_ps( x2, xfm1 );

    // sqrt2 * over + 2 * sqrt2 * under
    __m128 xavg = _mm_mul_ps( _mm_set1_ps( 1/( 3 * 0.629960524947437 ) * 0.999852 ),
                                _mm_add_ps( xover, xunder ) );

    xavg = _mm_mul_ps( xavg, _mm_rsqrt_ps( xavg ) );
    xavg = _mm_mul_ps( xavg, _mm_rsqrt_ps( xavg ) );
    return xavg;
}

__m128 mm_succ_ps( __m128 arg ) {
    return (__m128) _mm_add_epi32( (__m128i) arg, _mm_set1_epi32( 4 ) );
}

void test_pow( double p, __m128 (*f)( __m128 ) ) {
    __m128 arg;

    for ( arg = _mm_set1_ps( FLT_MIN / FLT_EPSILON );
            ! isfinite( _mm_cvtss_f32( f( arg ) ) );
            arg = mm_succ_ps( arg ) ) ;

    for ( ; _mm_cvtss_f32( f( arg ) ) == 0;
            arg = mm_succ_ps( arg ) ) ;

    std::printf( "Domain from %g\n", _mm_cvtss_f32( arg ) );

    int n;
    int const bucket_size = 1 << 25;
    do {
        float max_error = 0;
        double total_error = 0, cum_error = 0;
        for ( n = 0; n != bucket_size; ++ n ) {
            float result = _mm_cvtss_f32( f( arg ) );

            if ( ! isfinite( result ) ) break;

            float actual = ::powf( _mm_cvtss_f32( arg ), p );

            float error = ( result - actual ) / actual;
            cum_error += error;
            error = std::abs( error );
            max_error = std::max( max_error, error );
            total_error += error;

            arg = mm_succ_ps( arg );
        }

        std::printf( "error max = %8g\t" "avg = %8g\t" "|avg| = %8g\t" "to %8g\n",
                    max_error, cum_error / n, total_error / n, _mm_cvtss_f32( arg ) );
    } while ( n == bucket_size );
}

int main() {
    std::printf( "4 insn x^12/5:\n" );
    test_pow( 12./5, & fastpow< 12, 5, 1059, 1000 > );
    std::printf( "14 insn x^12/5:\n" );
    test_pow( 12./5, & pow125_4 );
    std::printf( "6 insn x^5/12:\n" );
    test_pow( 5./12, & pow512_2 );
    std::printf( "14 insn x^5/12:\n" );
    test_pow( 5./12, & pow512_4 );
}

Ausgabe:

4 insn x^12/5:
Domain from 1.36909e-23
error max =      inf    avg =      inf  |avg| =      inf    to 8.97249e-19
error max =  2267.14    avg =  139.175  |avg| =  139.193    to 5.88021e-14
error max = 0.123606    avg = -0.000102963  |avg| = 0.0371122   to 3.85365e-09
error max = 0.123607    avg = -0.000108978  |avg| = 0.0368548   to 0.000252553
error max =  0.12361    avg = 7.28909e-05   |avg| = 0.037507    to  16.5513
error max = 0.123612    avg = -0.000258619  |avg| = 0.0365618   to 1.08471e+06
error max = 0.123611    avg = 8.70966e-05   |avg| = 0.0374369   to 7.10874e+10
error max =  0.12361    avg = -0.000103047  |avg| = 0.0371122   to 4.65878e+15
error max = 0.123609    avg =      nan  |avg| =      nan    to 1.16469e+16
14 insn x^12/5:
Domain from 1.42795e-19
error max =      inf    avg =      nan  |avg| =      nan    to 9.35823e-15
error max = 0.000936462 avg = 2.0202e-05    |avg| = 0.000133764 to 6.13301e-10
error max = 0.000792752 avg = 1.45717e-05   |avg| = 0.000129936 to 4.01933e-05
error max = 0.000791785 avg = 7.0132e-06    |avg| = 0.000129923 to  2.63411
error max = 0.000787589 avg = 1.20745e-05   |avg| = 0.000129347 to   172629
error max = 0.000786553 avg = 1.62351e-05   |avg| = 0.000132397 to 1.13134e+10
error max = 0.000785586 avg = 8.25205e-06   |avg| = 0.00013037  to 6.98147e+12
6 insn x^5/12:
Domain from 9.86076e-32
error max = 0.0284339   avg = 0.000441158   |avg| = 0.00967327  to 6.46235e-27
error max = 0.0284342   avg = -5.79938e-06  |avg| = 0.00897913  to 4.23516e-22
error max = 0.0284341   avg = -0.000140706  |avg| = 0.00897084  to 2.77556e-17
error max = 0.028434    avg = 0.000440504   |avg| = 0.00967325  to 1.81899e-12
error max = 0.0284339   avg = -6.11153e-06  |avg| = 0.00897915  to 1.19209e-07
error max = 0.0284298   avg = -0.000140597  |avg| = 0.00897084  to 0.0078125
error max = 0.0284371   avg = 0.000439748   |avg| = 0.00967319  to      512
error max = 0.028437    avg = -7.74294e-06  |avg| = 0.00897924  to 3.35544e+07
error max = 0.0284369   avg = -0.000142036  |avg| = 0.00897089  to 2.19902e+12
error max = 0.0284368   avg = 0.000439183   |avg| = 0.0096732   to 1.44115e+17
error max = 0.0284367   avg = -7.41244e-06  |avg| = 0.00897923  to 9.44473e+21
error max = 0.0284366   avg = -0.000141706  |avg| = 0.00897088  to 6.1897e+26
error max = 0.485129    avg = -0.0401671    |avg| = 0.048422    to 4.05648e+31
error max = 0.994932    avg = -0.891494 |avg| = 0.891494    to 2.65846e+36
error max = 0.999329    avg =      nan  |avg| =      nan    to       -0
14 insn x^5/12:
Domain from 2.64698e-23
error max =  0.13556    avg = 0.00125936    |avg| = 0.00354677  to 1.73472e-18
error max = 0.000564988 avg = 2.51458e-06   |avg| = 0.000113709 to 1.13687e-13
error max = 0.000565065 avg = -1.49258e-06  |avg| = 0.000112553 to 7.45058e-09
error max = 0.000565143 avg = 1.5293e-06    |avg| = 0.000112864 to 0.000488281
error max = 0.000565298 avg = 2.76457e-06   |avg| = 0.000113713 to       32
error max = 0.000565453 avg = -1.61276e-06  |avg| = 0.000112561 to 2.09715e+06
error max = 0.000565531 avg = 1.42628e-06   |avg| = 0.000112866 to 1.37439e+11
error max = 0.000565686 avg = 2.71505e-06   |avg| = 0.000113715 to 9.0072e+15
error max = 0.000565763 avg = -1.56586e-06  |avg| = 0.000112415 to 1.84467e+19

Ich vermute, die Genauigkeit des genaueren 5/12 wird durch die rsqrt-Operation eingeschränkt.

22
Potatoswatter

Eine andere Antwort, weil diese sich sehr von meiner vorherigen Antwort unterscheidet, und diese ist schnell lodernd. Der relative Fehler ist 3e-8. Wünschen Sie mehr Genauigkeit? Fügen Sie ein paar weitere Chebychev-Begriffe hinzu. Es ist am besten, die Reihenfolge ungerade zu halten, da dies zu einer kleinen Diskontinuität zwischen 2 ^ n-Epsilon und 2 ^ n + Epsilon führt.

#include <stdlib.h>
#include <math.h>

// Returns x^(5/12) for x in [1,2), to within 3e-8 (relative error).
// Want more precision? Add more Chebychev polynomial coefs.
double pow512norm (
   double x)
{
   static const int N = 8;

   // Chebychev polynomial terms.
   // Non-zero terms calculated via
   //   integrate (2/pi)*ChebyshevT[n,u]/sqrt(1-u^2)*((u+3)/2)^(5/12)
   //   from -1 to 1
   // Zeroth term is similar except it uses 1/pi rather than 2/pi.
   static const double Cn[N] = { 
       1.1758200232996901923,
       0.16665763094889061230,
      -0.0083154894939042125035,
       0.00075187976780420279038,
      // Wolfram alpha doesn't want to compute the remaining terms
      // to more precision (it times out).
      -0.0000832402,
       0.0000102292,
      -1.3401e-6,
       1.83334e-7};

   double Tn[N];

   double u = 2.0*x - 3.0;

   Tn[0] = 1.0;
   Tn[1] = u;
   for (int ii = 2; ii < N; ++ii) {
      Tn[ii] = 2*u*Tn[ii-1] - Tn[ii-2];
   }   

   double y = 0.0;
   for (int ii = N-1; ii >= 0; --ii) {
      y += Cn[ii]*Tn[ii];
   }   

   return y;
}


// Returns x^(5/12) to within 3e-8 (relative error).
double pow512 (
   double x)
{
   static const double pow2_512[12] = {
      1.0,
      pow(2.0, 5.0/12.0),
      pow(4.0, 5.0/12.0),
      pow(8.0, 5.0/12.0),
      pow(16.0, 5.0/12.0),
      pow(32.0, 5.0/12.0),
      pow(64.0, 5.0/12.0),
      pow(128.0, 5.0/12.0),
      pow(256.0, 5.0/12.0),
      pow(512.0, 5.0/12.0),
      pow(1024.0, 5.0/12.0),
      pow(2048.0, 5.0/12.0)
   };

   double s;
   int iexp;

   s = frexp (x, &iexp);
   s *= 2.0;
   iexp -= 1;

   div_t qr = div (iexp, 12);
   if (qr.rem < 0) {
      qr.quot -= 1;
      qr.rem += 12;
   }

   return ldexp (pow512norm(s)*pow2_512[qr.rem], 5*qr.quot);
}

Nachtrag: Was ist hier los?
Auf Anfrage wird im Folgenden erläutert, wie der obige Code funktioniert.

Übersicht
Der obige Code definiert zwei Funktionen, double pow512norm (double x) und double pow512 (double x). Letzteres ist der Einstiegspunkt in die Suite. Dies ist die Funktion, die der Benutzercode aufrufen sollte, um x ^ (5/12) zu berechnen. Die Funktion pow512norm(x) verwendet Chebyshev-Polynome, um sich an x ​​^ (5/12) zu orientieren, jedoch nur für x im Bereich [1,2]. (Verwenden Sie pow512norm(x) für Werte von x außerhalb dieses Bereichs, und das Ergebnis ist Müll.)

Die Funktion pow512(x) teilt den eingehenden xin ein Paar (double s, int n) auf, so dass x = s * 2^n und 1 ≤ sname __ <2 ist. Eine weitere Partitionierung von nin (int q, unsigned int r), so dass n = 12*q + r und rkleiner als 12 ist, lässt mich das Problem des Auffindens von x ^ (5/12) in Teile aufteilen:

  1. x^(5/12)=(s^(5/12))*((2^n)^(5/12)) via (u v) ^ a = (u ^ a) (v ^ a) für positive u, v und reelle a.
  2. s^(5/12) wird über pow512norm(s) berechnet.
  3. (2^n)^(5/12)=(2^(12*q+r))^(5/12) durch Vertretung.
  4. 2^(12*q+r)=(2^(12*q))*(2^r) über u^(a+b)=(u^a)*(u^b) für positive u, real a, b.
  5. (2^(12*q+r))^(5/12)=(2^(5*q))*((2^r)^(5/12)) über einige weitere Manipulationen.
  6. (2^r)^(5/12) wird von der Nachschlagetabelle pow2_512 berechnet.
  7. Berechnen Sie pow512norm(s)*pow2_512[qr.rem] und wir sind fast da. Hier ist qr.rem der in Schritt 3 oben berechnete rWert. Sie müssen dies nur mit 2^(5*q) multiplizieren, um das gewünschte Ergebnis zu erzielen.
  8. Genau das macht die mathematische Bibliotheksfunktion ldexpname__.

Funktionsannäherung
Ziel ist es, eine leicht berechenbare Näherung von f (x) = x ^ (5/12) zu finden, die für das anstehende Problem "gut genug" ist. Unsere Annäherung sollte in gewissem Sinne nahe an f(x) liegen. Rhetorische Frage: Was bedeutet "nahe"? Zwei konkurrierende Interpretationen minimieren den mittleren quadratischen Fehler gegenüber dem minimalen absoluten Fehler.

Ich werde eine Aktienmarktanalogie verwenden, um den Unterschied zwischen diesen beiden zu beschreiben. Angenommen, Sie möchten für Ihren eventuellen Ruhestand sparen. Wenn Sie in Ihren Zwanzigern sind, sollten Sie am besten in Aktien oder Börsenfonds investieren. Dies ist darauf zurückzuführen, dass der Aktienmarkt über einen ausreichend langen Zeitraum im Durchschnitt alle anderen Anlageformen übertrifft. Wir haben jedoch alle Zeiten erlebt, in denen es sehr schlecht ist, Geld in Aktien zu investieren. Wenn Sie in den Fünfzigern oder Sechzigern sind (oder in den Vierzigern, wenn Sie jung in Rente gehen möchten), müssen Sie etwas konservativer investieren. Diese Abwärtsbewegungen können sich auf Ihr Altersportfolio auswirken.

Zurück zur Annäherung an die Funktion: Als Verbraucher einer Annäherung machen Sie sich normalerweise Sorgen um den Worst-Case-Fehler und nicht um die Leistung "im Durchschnitt". Verwenden Sie eine Näherung, die so konstruiert ist, dass sie im Durchschnitt die beste Leistung liefert (z. B. kleinste Quadrate), und das Murphy-Gesetz schreibt vor, dass Ihr Programm eine ganze Menge Zeit damit verbringen wird, die Näherung genau dort zu verwenden, wo die Leistung weit unter dem Durchschnitt liegt. Was Sie wollen, ist eine minimax-Approximation, die den maximalen absoluten Fehler in einer Domäne minimiert. Für eine gute Mathematikbibliothek wird ein Minimax-Ansatz und nicht ein Verfahren der kleinsten Quadrate verwendet, da die Autoren der Mathematikbibliothek dadurch eine garantierte Leistung ihrer Bibliothek erhalten.

Mathematische Bibliotheken verwenden normalerweise ein Polynom oder ein rationales Polynom, um eine Funktion f(x) über eine Domäne a ≤ x ≤ b anzunähern. Angenommen, die Funktion f(x) ist für diese Domäne analytisch, und Sie möchten die Funktion durch ein Polynom p(x) vom Grad N approximieren. Für einen bestimmten Grad N gibt es einige magisches, einzigartiges Polynom p(x) so, dass p(x)-f(x) N + 2 Extrema über [a, b] hat und die absoluten Werte dieser Werte sind N + 2-Extreme sind alle gleich. Das Finden dieses magischen Polynoms p(x) ist der heilige Gral der Funktionsapproximatoren.

Ich habe diesen heiligen Gral nicht für dich gefunden. Ich habe stattdessen eine Chebyshev-Annäherung verwendet. Die Chebyshev-Polynome der ersten Art sind ein orthogonaler (aber nicht orthonormaler) Satz von Polynomen mit einigen sehr schönen Merkmalen, wenn es um die Funktionsannäherung geht. Die Chebyshev-Annäherung liegt oft sehr nahe an diesem magischen Polynom p (x). (In der Tat beginnt der Remez-Austauschalgorithmus, der feststellt, dass das Polynom des Heiligen Grals normalerweise mit einer Chebyshev-Näherung beginnt.)

pow512norm (x)
Diese Funktion verwendet die Chebyshev-Näherung, um ein Polynom p * (x) zu finden, das x ^ (5/12) approximiert. Hier verwende ich p * (x), um diese Chebyshev-Approximation von dem oben beschriebenen magischen Polynom p(x) zu unterscheiden. Die Chebyshev-Näherung p * (x) ist leicht zu finden; p(x) zu finden ist ein Bär. Die Chebyshev-Approximation p * (x) ist sum_i Cn [i] * Tn (i, x), wobei Cn [i] die Chebyshev-Koeffizienten und Tn (i, x) die bei x bewerteten Chebyshev-Polynome sind.

Ich habe Wolfram alpha verwendet, um die Chebyshev-Koeffizienten Cnfür mich zu finden. Zum Beispiel berechnet Cn[1] . Das erste Feld nach dem Eingabefeld hat die gewünschte Antwort, in diesem Fall 0,166658. Das sind nicht so viele Ziffern, wie ich möchte. Klicken Sie auf "mehr Ziffern" und voila, Sie erhalten eine ganze Reihe weiterer Ziffern. Wolfram alpha ist kostenlos; Es gibt eine Begrenzung für die Berechnung. Es trifft diese Grenze bei Bedingungen höherer Ordnung. (Wenn Sie mathematica kaufen oder darauf zugreifen können, können Sie diese Koeffizienten höherer Ordnung mit hoher Genauigkeit berechnen.)

Die Chebyshev-Polynome Tn (x) werden im Array Tnberechnet. Abgesehen von etwas, das dem magischen Polynom p (x) sehr nahe kommt, besteht ein weiterer Grund für die Verwendung der Chebyshev-Näherung darin, dass die Werte dieser Chebyshev-Polynome leicht berechnet werden: Beginnen Sie mit Tn[0]=1 und Tn[1]=x und berechnen Sie Tn[i]=2*x*Tn[i-1] - Tn[i-2] anschließend iterativ. (Ich habe 'ii' als Indexvariable anstelle von 'i' in meinem Code verwendet. Ich verwende niemals 'i' als Variablennamen. Wie viele Wörter in der englischen Sprache haben ein 'i' im Wort? Wie viele haben zwei? aufeinander folgende "Ich?"

pow512 (x)
pow512 ist die Funktion, die der Benutzercode aufrufen soll. Die Grundlagen dieser Funktion wurden bereits oben beschrieben. Noch ein paar Details: Die Mathematik-Bibliotheksfunktion frexp(x) gibt den Bedeutungswert sund den Exponent iexpfür die Eingabe xzurück. (Kleinere Ausgabe: Ich möchte szwischen 1 und 2 zur Verwendung mit pow512norm, aber frexpgibt einen Wert zwischen 0,5 und 1 zurück.) Die mathematische Bibliotheksfunktion divgibt den Quotienten und den Rest für die Integer-Division in einem Quellvorgang zurück. Schließlich verwende ich die mathematische Bibliotheksfunktion ldexpname__, um die drei Teile zu einer endgültigen Antwort zusammenzusetzen.

32
David Hammen

Ian Stephenson schrieb diesen Code , von dem er behauptet, dass er pow() übertrifft. Er beschreibt die Idee wie folgt:

Pow wird im Wesentlichen mit Logs implementiert: pow(a,b)=x(logx(a)*b). wir brauchen also ein schnelles Log und einen schnellen Exponenten - es spielt keine Rolle, was x ist, also verwenden wir 2. Der Trick besteht darin, dass ein Gleitpunkt Nummer ist bereits im Protokollformat vorhanden:

a=M*2E

Wenn Sie das Protokoll von beiden Seiten aufnehmen, erhalten Sie:

log2(a)=log2(M)+E

oder einfacher:

log2(a)~=E

Mit anderen Worten, wenn wir die schwebende -Punktdarstellung einer Zahl nehmen und Den Exponent extrahieren, haben wir Etwas, das einen guten Ausgangspunkt Als sein Log. Es stellt sich heraus, dass, wenn wir Dies tun, indem sie die Bitmuster massieren, Die Mantisse letztendlich eine gute Annäherung an den Fehler gibt und sie Ziemlich gut wirkt Gut.

Dies sollte für einfache Beleuchtungsberechnungen gut genug sein, aber wenn Sie Etwas Besseres benötigen, können Sie Die Mantissa entnehmen und diese verwenden, um A zu berechnen quadratischer Korrekturfaktor , der ziemlich genau ist.

20
PengOne

Zunächst einmal wird die Verwendung von Schwimmern heutzutage auf den meisten Maschinen nicht viel kosten. In der Tat können Doppelte schneller sein. Ihre Leistung 1,0/2,4 ist 5/12 oder 1/3 * (1 + 1/4). Obwohl dies cbrt (einmalig) und sqrt (zweimal!) Aufruft, ist es dennoch doppelt so schnell wie mit pow (). (Optimierung: -O3, Compiler: i686-Apple-darwin10-g ++ - 4.2.1). 

#include <math.h> // cmath does not provide cbrt; C99 does.
double xpow512 (double x) {
  double cbrtx = cbrt(x);
  return cbrtx*sqrt(sqrt(cbrtx));
}
16
David Hammen

Dies kann Ihre Frage möglicherweise nicht beantworten.

Der 2.4f und 1/2.4f machen mich sehr misstrauisch, weil genau dies die Kräfte sind, die zum Konvertieren zwischen sRGB und einem linearen RGB-Farbraum verwendet werden. Sie versuchen also tatsächlich, das / speziell zu optimieren. Ich weiß nicht, weshalb Ihre Frage möglicherweise nicht beantwortet wird.

Wenn dies der Fall ist, versuchen Sie es mit einer Nachschlagetabelle. So etwas wie:

__attribute__((aligned(64))
static const unsigned short SRGB_TO_LINEAR[256] = { ... };
__attribute__((aligned(64))
static const unsigned short LINEAR_TO_SRGB[256] = { ... };

void apply_lut(const unsigned short lut[256], unsigned char *src, ...

Wenn Sie 16-Bit-Daten verwenden, ändern Sie sie entsprechend. Ich würde die Tabelle sowieso auf 16 Bit setzen, damit Sie das Ergebnis bei 8-Bit-Daten bei Bedarf rastern können. Dies funktioniert offensichtlich nicht sehr gut, wenn Ihre Daten von Anfang an Fließkomma sind - es ist jedoch nicht wirklich sinnvoll, sRGB-Daten in Fließkommazahlen zu speichern. Daher können Sie auch zuerst in 16-Bit/8-Bit konvertieren und führen Sie dann die Konvertierung von linear in sRGB durch.

(Der Grund, aus dem sRGB als Fließkomma nicht sinnvoll ist, besteht darin, dass HDR linear sein sollte und sRGB nur zum Speichern auf der Festplatte oder zum Anzeigen auf dem Bildschirm geeignet ist, für die Manipulation jedoch nicht.)

15
Dietrich Epp

Binomial series berücksichtigt einen konstanten Exponenten, aber Sie können ihn nur verwenden, wenn Sie alle Eingaben auf den Bereich [1,2] normalisieren können. (Beachten Sie, dass (1 + x) ^ a berechnet wird). Sie müssen einige Analysen durchführen, um zu entscheiden, wie viele Begriffe Sie für Ihre gewünschte Genauigkeit benötigen.

3
zvrba

Ich werde die Frage beantworten, die Sie wirklich stellen wollten, nämlich wie Sie eine schnelle sRGB-lineare RGB-Konvertierung durchführen können. Um dies präzise und effizient durchzuführen, können wir Polynomapproximationen verwenden. Die folgenden Polynomapproximierungen wurden mit Sollya erstellt und weisen im ungünstigsten Fall einen relativen Fehler von 0,0154% auf.

inline double poly7(double x, double a, double b, double c, double d,
                              double e, double f, double g, double h) {
    double ab, cd, ef, gh, abcd, efgh, x2, x4;
    x2 = x*x; x4 = x2*x2;
    ab = a*x + b; cd = c*x + d;
    ef = e*x + f; gh = g*x + h;
    abcd = ab*x2 + cd; efgh = ef*x2 + gh;
    return abcd*x4 + efgh;
}

inline double srgb_to_linear(double x) {
    if (x <= 0.04045) return x / 12.92;

    // Polynomial approximation of ((x+0.055)/1.055)^2.4.
    return poly7(x, 0.15237971711927983387,
                   -0.57235993072870072762,
                    0.92097986411523535821,
                   -0.90208229831912012386,
                    0.88348956209696805075,
                    0.48110797889132134175,
                    0.03563925285274562038,
                    0.00084585397227064120);
}

inline double linear_to_srgb(double x) {
    if (x <= 0.0031308) return x * 12.92;

    // Piecewise polynomial approximation (divided by x^3)
    // of 1.055 * x^(1/2.4) - 0.055.
    if (x <= 0.0523) return poly7(x, -6681.49576364495442248881,
                                      1224.97114922729451791383,
                                      -100.23413743425112443219,
                                         6.60361150127077944916,
                                         0.06114808961060447245,
                                        -0.00022244138470139442,
                                         0.00000041231840827815,
                                        -0.00000000035133685895) / (x*x*x);

    return poly7(x, -0.18730034115395793881,
                     0.64677431008037400417,
                    -0.99032868647877825286,
                     1.20939072663263713636,
                     0.33433459165487383613,
                    -0.01345095746411287783,
                     0.00044351684288719036,
                    -0.00000664263587520855) / (x*x*x);
}

Und die Sollya-Eingabe, aus der die Polynome generiert wurden:

suppressmessage(174);
f = ((x+0.055)/1.055)^2.4;
p0 = fpminimax(f, 7, [|D...|], [0.04045;1], relative);
p = fpminimax(f/(p0(1)+1e-18), 7, [|D...|], [0.04045;1], relative);
print("relative:", dirtyinfnorm((f-p)/f, [s;1]));
print("absolute:", dirtyinfnorm((f-p), [s;1]));
print(canonical(p));

s = 0.0523;
z = 3;
f = 1.055 * x^(1/2.4) - 0.055;

p = fpminimax(1.055 * (x^(z+1/2.4) - 0.055*x^z/1.055), 7, [|D...|], [0.0031308;s], relative)/x^z;
print("relative:", dirtyinfnorm((f-p)/f, [0.0031308;s]));
print("absolute:", dirtyinfnorm((f-p), [0.0031308;s]));
print(canonical(p));

p = fpminimax(1.055 * (x^(z+1/2.4) - 0.055*x^z/1.055), 7, [|D...|], [s;1], relative)/x^z;
print("relative:", dirtyinfnorm((f-p)/f, [s;1]));
print("absolute:", dirtyinfnorm((f-p), [s;1]));
print(canonical(p));
2
orlp

Das Folgende ist eine Idee, die Sie mit jeder der schnellen Berechnungsmethoden verwenden können. Ob es schneller geht, hängt davon ab, wie Ihre Daten ankommen. Sie können die Tatsache verwenden, dass Sie, wenn Sie x und pow(x, n) kennen, die Änderungsrate der Potenz verwenden können, um eine vernünftige Annäherung von pow(x + delta, n) für kleine delta mit einem einzigen Multiplikator zu berechnen und (mehr oder weniger) hinzuzufügen. Wenn aufeinanderfolgende Werte, die Sie Ihren Leistungsfunktionen zuführen, nahe genug sind, würde dies die vollen Kosten der genauen Berechnung über mehrere Funktionsaufrufe hinweg amortisieren. Beachten Sie, dass Sie keine Extra-Berechnung benötigen, um die Ableitung zu erhalten. Sie können dies erweitern, um die zweite Ableitung zu verwenden, sodass Sie eine quadratische verwenden können, wodurch die delta erhöht wird, die Sie verwenden könnten, und dennoch die gleiche Genauigkeit erhalten.

1
Permaquid

Für Exponenten von 2.4 können Sie entweder eine Nachschlagetabelle für alle Ihre 2.4-Werte und die Funktion lirp oder eine Funktion höherer Ordnung erstellen, um die Zwischenwerte einzugeben, wenn die Tabelle nicht genau genug war (im Grunde eine große Protokolltabelle).

Oder setzen Sie den Quadratwert * auf 2/5s, der den anfänglichen Quadratwert aus der ersten Hälfte der Funktion und dann den fünften Wurzelwert annehmen könnte. Für die fünfte Wurzel können Sie Newton verwenden oder einen anderen schnellen Näherungswert ausführen. Wenn Sie ehrlich gesagt an diesem Punkt angelangt sind, ist es wahrscheinlich besser, wenn Sie die exp- und log-Funktionen mit den entsprechenden abgekürzten Serienfunktionen selbst ausführen.

1
Michael Dorgan

Normalerweise wird powf(x, p) = x^p dadurch gelöst, dass x als x=2^(log2(x)) making powf(x,p) = 2^(p*log2(x)) umgeschrieben wird, wodurch das Problem in zwei Annäherungen exp2() & log2() umgewandelt wird. Dies hat den Vorteil, dass mit größeren Leistungen p gearbeitet wird. Der Nachteil besteht jedoch darin, dass dies nicht die optimale Lösung für eine konstante Leistung p und für eine angegebene Eingabebindung 0 ≤ x ≤ 1 ist.

Bei der Potenz p > 1 ist die Antwort ein triviales Minimax-Polynom über dem gebundenen 0 ≤ x ≤ 1, was für p = 12/5 = 2.4 der Fall ist, wie unten zu sehen ist:

float pow12_5(float x){
    float mp;
    // Minimax horner polynomials for x^(5/12), Note: choose the accurarcy required then implement with fma() [Fused Multiply Accumulates]
    // mp = 0x4.a84a38p-12 + x * (-0xd.e5648p-8 + x * (0xa.d82fep-4 + x * 0x6.062668p-4)); // 1.13705697e-3
    mp = 0x1.117542p-12 + x * (-0x5.91e6ap-8 + x * (0x8.0f50ep-4 + x * (0xa.aa231p-4 + x * (-0x2.62787p-4))));  // 2.6079002e-4
    // mp = 0x5.a522ap-16 + x * (-0x2.d997fcp-8 + x * (0x6.8f6d1p-4 + x * (0xf.21285p-4 + x * (-0x7.b5b248p-4 + x * 0x2.32b668p-4))));  // 8.61377e-5
    // mp = 0x2.4f5538p-16 + x * (-0x1.abcdecp-8 + x * (0x5.97464p-4 + x * (0x1.399edap0 + x * (-0x1.0d363ap0 + x * (0xa.a54a3p-4 + x * (-0x2.e8a77cp-4))))));  // 3.524655e-5
    return(mp);
}

Wenn p < 1 jedoch die Minimax-Näherung über den gebundenen 0 ≤ x ≤ 1 nicht angemessen auf die gewünschte Genauigkeit konvergiert. Eine Option [nicht wirklich] besteht darin, das Problem y=x^p=x^(p+m)/x^m neu zu schreiben, wobei m=1,2,3 eine positive ganze Zahl ist, wodurch die neue Potenzannäherung p > 1 erfolgt. Dies führt jedoch zu einer Division, die von Natur aus langsamer ist. 

Es gibt jedoch eine andere Möglichkeit, die Eingabe x als Gleitkommaexponent und Mantissenform zu zerlegen: 

x = mx* 2^(ex) where 1 ≤ mx < 2
y = x^(5/12) = mx^(5/12) * 2^((5/12)*ex), let ey = floor(5*ex/12), k = (5*ex) % 12
  = mx^(5/12) * 2^(k/12) * 2^(ey)

Die Minimax-Annäherung von mx^(5/12) über 1 ≤ mx < 2 konvergiert jetzt viel schneller als zuvor ohne Division, erfordert jedoch eine 12-Punkt-LUT für die 2^(k/12). Der Code ist unten:

float powk_12LUT[] = {0x1.0p0, 0x1.0f38fap0, 0x1.1f59acp0,  0x1.306fep0, 0x1.428a3p0, 0x1.55b81p0, 0x1.6a09e6p0, 0x1.7f910ep0, 0x1.965feap0, 0x1.ae89fap0, 0x1.c823ep0, 0x1.e3437ep0};
float pow5_12(float x){
    union{float f; uint32_t u;} v, e2;
    float poff, m, e, ei;
    int xe;

    v.f = x;
    xe = ((v.u >> 23) - 127);

    if(xe < -127) return(0.0f);

    // Calculate remainder k in 2^(k/12) to find LUT
    e = xe * (5.0f/12.0f);
    ei = floorf(e);
    poff = powk_12LUT[(int)(12.0f * (e - ei))];

    e2.u = ((int)ei + 127) << 23;   // Calculate the exponent
    v.u = (v.u & ~(0xFFuL << 23)) | (0x7FuL << 23); // Normalize exponent to zero

    // Approximate mx^(5/12) on [1,2), with appropriate degree minimax
    // m = 0x8.87592p-4 + v.f * (0x8.8f056p-4 + v.f * (-0x1.134044p-4));    // 7.6125e-4
    // m = 0x7.582138p-4 + v.f * (0xb.1666bp-4 + v.f * (-0x2.d21954p-4 + v.f * 0x6.3ea0cp-8));  // 8.4522726e-5
    m = 0x6.9465cp-4 + v.f * (0xd.43015p-4 + v.f * (-0x5.17b2a8p-4 + v.f * (0x1.6cb1f8p-4 + v.f * (-0x2.c5b76p-8))));   // 1.04091259e-5
    // m = 0x6.08242p-4 + v.f * (0xf.352bdp-4 + v.f * (-0x7.d0c1bp-4 + v.f * (0x3.4d153p-4 + v.f * (-0xc.f7a42p-8 + v.f * 0x1.5d840cp-8))));    // 1.367401e-6

    return(m * poff * e2.f);
}
0
nimig18