Оптимизация для pow () с константной нецелой экспонентой?

В духе взлома IEEE 754, вот еще одно решение, более быстрое и менее «волшебное». Он обеспечивает погрешность 0,08% примерно за десяток тактовых циклов (для случая p = 2,4 на процессоре Intel Merom).

Числа с плавающей запятой изначально были изобретены как приближение к логарифмам, поэтому вы можете использовать целочисленное значение как приближение log2. Это в некоторой степени легко достижимо, применяя инструкцию преобразования из целого числа к значению с плавающей запятой, чтобы получить другое значение с плавающей запятой.

Чтобы завершить вычисление pow, вы можете умножить его на постоянный коэффициент и преобразовать логарифм обратно с помощью инструкции преобразования в целое число. В SSE соответствующие инструкции cvtdq2ps и cvtps2dq.

Однако все не так просто. Поле экспоненты в IEEE 754 подписано со значением смещения 127, представляющим экспоненту, равную нулю. Это смещение необходимо удалить перед умножением логарифма и снова добавить перед возведением в степень. Кроме того, корректировка смещения путем вычитания не будет работать на нуле. К счастью, обе корректировки могут быть достигнуты путем предварительного умножения на постоянный коэффициент.

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 ) - постоянный коэффициент. Эта функция довольно специализированная: она не будет работать с малыми дробными показателями степени, потому что постоянный множитель экспоненциально растет с обратным показателем экспоненты и будет переполняться. Это не сработает с отрицательными показателями. Большие экспоненты приводят к высокой ошибке, потому что биты мантиссы смешиваются с битами экспоненты при умножении.

Но это всего лишь 4 быстрые инструкции. Предварительное умножение, преобразование из «целого числа» (в логарифм), умножение в степени, преобразование в «целое число» (из логарифма). В этой реализации SSE преобразования происходят очень быстро. Мы также можем втиснуть дополнительный постоянный коэффициент в первое умножение.

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;
}

Несколько испытаний с показателем экспоненты = 2,4 показывают, что это постоянно завышение примерно на 5%. (Процедура всегда гарантированно переоценивает.) Вы можете просто умножить на 0,95, но еще несколько инструкций дадут нам точность около 4 десятичных знаков, чего должно хватить для графики.

Ключ состоит в том, чтобы сопоставить завышенную оценку с заниженной и взять среднее значение.

• Вычислить x ^ 0,8: четыре инструкции, ошибка ~ + 3%.
• Вычислите x ^ -0,4: один rsqrtps. (Это достаточно точно, но приносит в жертву способность работать с нулем.)
• Вычислить x ^ 0,4: один mulps.
• Вычислите x ^ -0,2: один rsqrtps.
• Вычислить x ^ 2: one mulps.
• Вычислить x ^ 3: one mulps.
• x ^ 2,4 = x ^ 2 * x ^ 0,4: один mulps. Это завышенная оценка.
• x ^ 2,4 = x ^ 3 * x ^ -0,4 * x ^ -0,2: два mulps. Это недооценка.
• Среднее значение выше: один addps, один mulps.

Подсчет инструкций: четырнадцать, включая два преобразования с задержкой = 5 и две оценки обратного квадратного корня с пропускной способностью = 4.

Чтобы правильно вычислить среднее значение, мы хотим взвесить оценки по их ожидаемым ошибкам. Недооценка увеличивает ошибку до степени 0,6 против 0,4, поэтому мы ожидаем, что она будет в 1,5 раза ошибочной. Взвешивание не добавляет никаких инструкций; это можно сделать в предварительном факторе. Обозначая коэффициент a: a ^ 0,5 = 1,5 a ^ -0,75, а a = 1,38316186.

Окончательная ошибка составляет около 0,015%, что на 2 порядка лучше, чем первоначальный результат fastpow. Время выполнения составляет около дюжины циклов для цикла занятости с volatile исходными и целевыми переменными… хотя итерации перекрываются, при реальном использовании также будет наблюдаться параллелизм на уровне инструкций. Учитывая SIMD, это пропускная способность одного скалярного результата за 3 цикла!

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 );
}

Что ж ... извините, я не смог опубликовать это раньше. И расширение его до x ^ 1 / 2.4 остается в качестве упражнения; v).

Обновить статистику

Я реализовал небольшой тестовый набор и два случая x ^{(5 ⁄ ₁₂)}, соответствующих приведенному выше.

#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 );
}

Выход:

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

Я подозреваю, что точность более точного 5/12 ограничена операцией rsqrt.

Другой ответ, потому что он сильно отличается от моего предыдущего ответа, и это очень быстро. Относительная погрешность 3e-8. Хотите большей точности? Добавьте еще пару чебычевских терминов. Лучше всего сохранить нечетный порядок, поскольку это приводит к небольшому разрыву между 2 ^ n-эпсилон и 2 ^ n + эпсилон.

#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);
}

Приложение: что здесь происходит?
В каждом запросе ниже объясняется, как работает приведенный выше код.

Обзор
Приведенный выше код определяет две функции: double pow512norm (double x) и double pow512 (double x). Последний является точкой входа в люкс; это функция, которую пользовательский код должен вызывать для вычисления x ^ (5/12). Функция pow512norm(x) использует полиномы Чебышева для аппроксимации x ^ (5/12), но только для x в диапазоне [1,2]. (Используйте pow512norm(x) для значений x вне этого диапазона, и результат будет мусором.)

Функция pow512(x) разбивает входящий x на пару (double s, int n), так что x = s * 2^n и такие, что 1≤s ‹2. Дальнейшее разделение n на (int q, unsigned int r) таким образом, что n = 12*q + r и r меньше 12, позволяет мне разделить проблему поиска x ^ (5/12) на части:

1. x^(5/12)=(s^(5/12))*((2^n)^(5/12)) через (u v) ^ a = (u ^ a) (v ^ a) для положительных u, v и вещественных a.
2. s^(5/12) рассчитывается через pow512norm(s).
3. (2^n)^(5/12)=(2^(12*q+r))^(5/12) путем подстановки.
4. 2^(12*q+r)=(2^(12*q))*(2^r) через u^(a+b)=(u^a)*(u^b) для положительного u, действительного a, b.
5. (2^(12*q+r))^(5/12)=(2^(5*q))*((2^r)^(5/12)) через еще несколько манипуляций.
6. (2^r)^(5/12) рассчитывается по таблице поиска pow2_512.
7. Вычислите pow512norm(s)*pow2_512[qr.rem], и мы почти у цели. Здесь qr.rem - значение r, вычисленное на шаге 3 выше. Все, что нужно, - это умножить это на 2^(5*q), чтобы получить желаемый результат.
8. Именно это и делает функция математической библиотеки ldexp.

Аппроксимация функций
Цель состоит в том, чтобы найти легко вычислимое приближение f (x) = x ^ (5/12), которое «достаточно хорошо» для рассматриваемой проблемы. Наше приближение должно быть в некотором смысле близким к f (x). Риторический вопрос: что означает «близко к»? Две конкурирующие интерпретации сводят к минимуму среднеквадратичную ошибку по сравнению с минимизацией максимальной абсолютной ошибки.

Я буду использовать аналогию с фондовым рынком, чтобы описать разницу между ними. Предположим, вы хотите накопить на свою возможную пенсию. Если вам за двадцать, лучше всего инвестировать в акции или фонды фондового рынка. Это потому, что в течение достаточно длительного периода времени фондовый рынок в среднем превосходит любую другую инвестиционную схему. Однако все мы видели случаи, когда вкладывать деньги в акции - это очень плохо. Если вам за пятьдесят или за шестьдесят (или за сорок, если вы хотите выйти на пенсию молодым), вам нужно инвестировать немного более консервативно. Эти спады могут серьезно отразиться на вашем пенсионном портфеле.

Вернемся к приближению функции: как потребитель некоторого приближения, вы обычно беспокоитесь о наихудшей ошибке, а не о производительности «в среднем». Используйте некоторое приближение, построенное для получения наилучшей производительности «в среднем» (например, методом наименьших квадратов), а закон Мерфи требует, чтобы ваша программа проводила много времени, используя приближение именно там, где производительность намного хуже, чем в среднем. Вам нужно минимаксное приближение, то есть то, что минимизирует максимальную абсолютную ошибку в некоторой области. Хорошая математическая библиотека будет использовать минимаксный подход, а не метод наименьших квадратов, потому что это позволяет авторам математической библиотеки обеспечить некоторую гарантированную производительность своей библиотеки.

Математические библиотеки обычно используют многочлен или рациональный многочлен для аппроксимации некоторой функции f (x) в некоторой области a≤x≤b. Предположим, что функция f (x) является аналитической в ​​этой области, и вы хотите аппроксимировать функцию некоторым многочленом p (x) степени N.Для данной степени N существует некоторый магический уникальный многочлен p (x) такой, что p ( x) -f (x) имеет N + 2 экстремума над [a, b] и такие, что абсолютные значения этих N + 2 экстремумов равны друг другу. Нахождение этого волшебного многочлена p (x) - это святой Грааль аппроксиматоров функций.

Я не нашел для тебя этого святого Грааля. Вместо этого я использовал приближение Чебышева. Многочлены Чебышева первого рода - это ортогональный (но не ортонормированный) набор многочленов с некоторыми очень хорошими особенностями, когда дело доходит до приближения функций. Приближение Чебышева часто очень близко к магическому многочлену p (x). (Фактически, алгоритм обмена Ремеза, который действительно находит этот полином святого Грааля, обычно начинается с приближения Чебышева.)

pow512norm (x)
Эта функция использует аппроксимацию Чебышева, чтобы найти некоторый многочлен p * (x), который аппроксимирует x ^ (5/12). Здесь я использую p * (x), чтобы отличить это приближение Чебышева от магического полинома p (x), описанного выше. Чебышёвское приближение p * (x) найти несложно; находка p (x) - медведь. Приближение Чебышева p * (x) - это sum_i Cn [i] * Tn (i, x), где Cn [i] - это коэффициенты Чебышева, а Tn (i, x) - полиномы Чебышева, вычисленные в x.

Я использовал Wolfram alpha, чтобы найти для меня коэффициенты Чебышева Cn. Например, вычисляется Cn[1]. В первом поле после поля ввода указан желаемый ответ, в данном случае 0,166658. Это не так много цифр, как хотелось бы. Нажмите на «больше цифр» и вуаля, вы получите намного больше цифр. Wolfram alpha бесплатна; существует ограничение на количество вычислений, которые он будет выполнять. Он достигает этого предела на условиях более высокого порядка. (Если вы купите программу mathematica или у вас есть доступ к ней, вы сможете вычислить эти старшие коэффициенты с высокой степенью точности.)

Многочлены Чебышева Tn (x) вычисляются в массиве Tn. Помимо предоставления чего-то очень близкого к магическому многочлену p (x), еще одна причина для использования приближения Чебышева состоит в том, что значения этих многочленов Чебышева легко вычисляются: начните с Tn[0]=1 и Tn[1]=x, а затем итеративно вычислите Tn[i]=2*x*Tn[i-1] - Tn[i-2]. (Я использовал в качестве индексной переменной «ii», а не «i» в своем коде. Я никогда не использую «i» в качестве имени переменной. Сколько слов в английском языке содержат букву «i»? последовательные 'i'?)

pow512 (x)
pow512 - функция, которую должен вызывать код пользователя. Я уже описал основы этой функции выше. Еще несколько деталей: функция математической библиотеки frexp(x) возвращает значение s и показатель iexp для входных данных x. (Незначительная проблема: я хочу s между 1 и 2 для использования с pow512norm, но frexp возвращает значение от 0,5 до 1.) Функция математической библиотеки div возвращает частное и остаток для целочисленного деления за один раз. Наконец, я использую математическую библиотечную функцию ldexp, чтобы соединить три части и сформировать окончательный ответ.

Ян Стивенсон написал этот код, который, по его словам, превосходит pow(). Он описывает идею следующим образом:

Pow в основном реализован с использованием журнала: pow(a,b)=x(logx(a)*b). поэтому нам нужен быстрый журнал и быстрая экспонента - не имеет значения, какой x, поэтому мы используем 2. Хитрость в том, что число с плавающей запятой уже находится в формате журнала:

a=M*2E

Взятие бревна с обеих сторон дает:

log2(a)=log2(M)+E

или проще:

log2(a)~=E

Другими словами, если мы возьмем представление числа с плавающей запятой и извлечем экспоненту, мы получим что-то, что является хорошей отправной точкой в ​​качестве его журнала. Оказывается, когда мы делаем это, массируя битовые шаблоны, мантисса дает хорошее приближение к ошибке, и это работает довольно хорошо.

Этого должно быть достаточно для простых расчетов освещения, но если вам нужно что-то получше, вы можете извлечь мантиссу и использовать ее для расчета коэффициента квадратичной коррекции, который довольно точен.

Во-первых, использование поплавков в настоящее время не принесет многого на большинстве машин. На самом деле, удвоение может быть быстрее. Ваша сила 1.0 / 2.4 равна 5/12 или 1/3 * (1 + 1/4). Несмотря на то, что это вызывает cbrt (один раз) и sqrt (дважды!), Он все равно в два раза быстрее, чем при использовании pow (). (Оптимизация: -O3, компилятор: 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));
}

Возможно, это не ответ на ваш вопрос.

2.4f и 1/2.4f вызывают у меня большие подозрения, потому что это как раз те возможности, которые используются для преобразования между sRGB и линейным цветовым пространством RGB. Так что, возможно, вы действительно пытаетесь оптимизировать именно это. Я не знаю, поэтому это может не ответить на ваш вопрос.

В этом случае попробуйте использовать справочную таблицу. Что-то вроде:

__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, ...

Если вы используете 16-битные данные, измените их соответствующим образом. Я бы в любом случае сделал таблицу 16-битной, чтобы вы могли при необходимости сглаживать результат при работе с 8-битными данными. Это, очевидно, не будет работать очень хорошо, если ваши данные с плавающей запятой для начала - но на самом деле нет смысла хранить данные sRGB с плавающей запятой, поэтому вы можете сначала преобразовать в 16-бит / 8-бит. а затем сделайте преобразование из линейного в sRGB.

(Причина, по которой sRGB не имеет смысла как числа с плавающей запятой, заключается в том, что HDR должен быть линейным, а sRGB удобен только для хранения на диске или отображения на экране, но не удобен для манипуляций.)

Биномиальный ряд учитывает постоянную экспоненту, но вы сможете использовать ее, только если сможете нормализовать весь ваш ввод до диапазона [1,2). (Обратите внимание, что он вычисляет (1 + x) ^ a). Вам нужно будет провести некоторый анализ, чтобы решить, сколько терминов вам нужно для желаемой точности.

Я отвечу на вопрос, который вы действительно хотели задать: как выполнить быстрое преобразование sRGB ‹-> в линейную RGB-подсветку. Чтобы сделать это точно и эффективно, мы можем использовать полиномиальные приближения. Следующие полиномиальные приближения были сгенерированы с помощью sollya и имеют относительную погрешность в худшем случае 0,0144%.

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);
}

И вход sollya, используемый для генерации полиномов:

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,4 вы можете либо создать таблицу поиска для всех значений 2,4 и lirp, либо, возможно, функцию более высокого порядка для заполнения значений in-betweem, если таблица недостаточно точна (в основном огромная таблица журнала).

Или, значение в квадрате * значение до 2/5, которое может взять начальное квадратное значение из первой половины функции, а затем 5-й корень из него. Для 5-го корня вы можете использовать Ньютон или какой-нибудь другой быстрый аппроксиматор, хотя, честно говоря, как только вы дойдете до этой точки, вам, вероятно, будет лучше просто выполнить функции exp и log с соответствующими сокращенными функциями ряда самостоятельно.

Следующее - идея, которую вы можете использовать с любым из быстрых методов расчета. Поможет ли это ускорить работу, зависит от того, как поступают ваши данные. Вы можете использовать тот факт, что, если вы знаете x и pow(x, n), вы можете использовать скорость изменения мощности для вычисления разумного приближения pow(x + delta, n) для малых delta с помощью одного умножения и сложения (более или менее). Если последовательные значения, которые вы вводите в свои функции мощности, достаточно близки друг к другу, это приведет к амортизации полной стоимости точных вычислений по нескольким вызовам функций. Обратите внимание, что вам не нужно дополнительное вычисление мощности, чтобы получить производную. Вы можете расширить это, чтобы использовать вторую производную, чтобы вы могли использовать квадратичную, что увеличило бы delta, которое вы могли бы использовать, и при этом получило бы ту же точность.

Таким образом, традиционно powf(x, p) = x^p решается путем переписывания x на x=2^(log2(x)), что делает powf(x,p) = 2^(p*log2(x)), что превращает проблему в два приближения exp2() и log2(). Это дает преимущество работы с большими мощностями p, однако недостатком является то, что это не оптимальное решение для постоянной мощности p и превышения указанной входной границы 0 ≤ x ≤ 1.

Когда степень p > 1, ответом является тривиальный минимаксный полином над границей 0 ≤ x ≤ 1, что имеет место для p = 12/5 = 2.4, как показано ниже:

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);
}

Однако, когда p < 1 минимаксное приближение за пределами 0 ≤ x ≤ 1 не сходится надлежащим образом к желаемой точности. Один из вариантов [не совсем] - это переписать задачу y=x^p=x^(p+m)/x^m, где m=1,2,3 - положительное целое число, делая новое приближение мощности p > 1, но это вводит деление, которое по своей сути является более медленным.

Однако есть еще один вариант, который состоит в том, чтобы разложить вход x на его показатель с плавающей запятой и форму мантиссы:

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)

Минимаксное приближение mx^(5/12) по 1 ≤ mx < 2 теперь сходится намного быстрее, чем раньше, без деления, но требует 12-точечного LUT для 2^(k/12). Код ниже:

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);
}