while utiliza una raiz que programa potencia numero elevar dev cómo cuadrada con calcule c++ c math optimization exponent

c++ - utiliza - Optimizaciones para pow() con const exponente no entero?



programa que calcule la potencia de un numero en c++ (10)

Tengo puntos calientes en mi código donde estoy haciendo pow() ocupando alrededor del 10-20% de mi tiempo de ejecución.

Mi entrada a pow(x,y) es muy específica, así que me pregunto si hay una manera de lanzar dos aproximaciones pow() (una para cada exponente) con un rendimiento superior:

  • Tengo dos exponentes constantes: 2.4 y 1 / 2.4.
  • Cuando el exponente es 2.4, x estará en el rango (0.090473935, 1.0).
  • Cuando el exponente es 1 / 2.4, x estará en el rango (0.0031308, 1.0).
  • Estoy usando vectores float SSE / AVX. Si se pueden aprovechar las especificaciones de la plataforma, ¡adelante!

Una tasa de error máxima de alrededor del 0,01% es ideal, aunque también me interesan los algoritmos de precisión completa (para float ).

Ya estoy usando una approximation pow() rápida, pero no tiene en cuenta estas limitaciones. ¿Es posible hacerlo mejor?

En la veta de pirateo IEEE 754, aquí hay otra solución que es más rápida y menos "mágica". Logra un margen de error de .08% en aproximadamente una docena de ciclos de reloj (para el caso de p = 2.4, en una CPU Intel Merom).

Los números de punto flotante se inventaron originalmente como una aproximación a los logaritmos, por lo que puede usar el valor entero como una aproximación de log2 . Esto es algo realizable de forma portátil aplicando la instrucción de conversión de un entero a un valor de punto flotante, para obtener otro valor de coma flotante.

Para completar el cálculo de pow , puede multiplicar por un factor constante y convertir el logaritmo en la instrucción de conversión a entero. En SSE, las instrucciones relevantes son cvtdq2ps y cvtps2dq .

Aunque no es tan simple. El campo de exponente en IEEE 754 está firmado, con un valor de sesgo de 127 que representa un exponente de cero. Este sesgo debe eliminarse antes de multiplicar el logaritmo y volverse a agregar antes de exponenciarse. Además, el ajuste de sesgo por sustracción no funcionará en cero. Afortunadamente, ambos ajustes se pueden lograr multiplicando por un factor constante de antemano.

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 ) es el factor constante. Esta función es bastante especializada: no funcionará con exponentes fraccionarios pequeños, porque el factor constante crece exponencialmente con el inverso del exponente y se desbordará. No funcionará con exponentes negativos. Los exponentes grandes conducen a un error alto, porque los bits de mantisa se mezclan con los bits de exponente por la multiplicación.

Pero, son solo 4 instrucciones rápidas de largo. Pre-multiplicar, convertir de "entero" (a logaritmo), multiplicar por potencia, convertir a "entero" (a partir del logaritmo). Las conversiones son muy rápidas en esta implementación de SSE. También podemos exprimir un coeficiente constante adicional en la primera multiplicación.

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

Algunos ensayos con exponente = 2.4 muestran que esto se sobreestima consistentemente en aproximadamente 5%. (Siempre se garantiza que la rutina sobreestimará). Simplemente podría multiplicar por 0.95, pero algunas instrucciones más nos darán unos 4 dígitos decimales de precisión, lo cual debería ser suficiente para los gráficos.

La clave es hacer coincidir la sobreestimación con una subestimación, y tomar el promedio.

  • Calcule x ^ 0.8: cuatro instrucciones, error ~ + 3%.
  • Calcule x ^ -0.4: una rsqrtps . (Esto es bastante preciso, pero sacrifica la capacidad de trabajar con cero).
  • Calcule x ^ 0.4: uno mulps .
  • Calcule x ^ -0.2: una rsqrtps .
  • Calcule x ^ 2: un mulps .
  • Calcule x ^ 3: un mulps .
  • x ^ 2.4 = x ^ 2 * x ^ 0.4: un mulps . Esta es la sobreestimacion
  • x ^ 2.4 = x ^ 3 * x ^ -0.4 * x ^ -0.2: dos mulps . Esta es la subestimacion
  • Promedio de lo anterior: un addps , un mulps .

Contador de instrucciones: catorce, incluidas dos conversiones con latencia = 5 y dos estimaciones recíprocas de raíz cuadrada con rendimiento = 4.

Para tomar correctamente el promedio, queremos ponderar las estimaciones por los errores esperados. La subestimación aumenta el error a una potencia de 0.6 frente a 0.4, por lo que esperamos que sea 1.5x errónea. La ponderación no agrega ninguna instrucción; se puede hacer en el pre-factor. Llamando al coeficiente a: a ^ 0.5 = 1.5 a ^ -0.75, y a = 1.38316186.

El error final es aproximadamente .015%, o 2 órdenes de magnitud mejor que el resultado inicial de fastpow . El tiempo de ejecución es de alrededor de una docena de ciclos para un ciclo ocupado con variables de origen y destino volatile ... aunque se superpone a las iteraciones, el uso en el mundo real también verá el paralelismo a nivel de instrucción. Considerando SIMD, ¡eso es un rendimiento de un resultado escalar por 3 ciclos!

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

Bueno ... lo siento, no pude publicar esto antes. Y extenderlo a x ^ 1 / 2.4 queda como ejercicio; v).

Actualizar con estadísticas

Implementé un pequeño arnés de prueba y dos x ( 5/12 ) casos correspondientes a los anteriores.

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

Salida:

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

Sospecho que la precisión del 5/12 más preciso está siendo limitado por la operación rsqrt .

En primer lugar, usar flotadores no va a comprar mucho en la mayoría de las máquinas hoy en día. De hecho, los dobles pueden ser más rápidos. Su potencia, 1.0 / 2.4, es 5/12 o 1/3 * (1 + 1/4). Aunque esto llama a cbrt (una vez) y sqrt (¡dos veces!), Sigue siendo dos veces más rápido que usar pow (). (Optimización: -O3, compilador: 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)); }

Esto podría no responder su pregunta.

El 2.4f y el 1/2.4f me hacen muy sospechoso, porque esos son exactamente los poderes utilizados para convertir entre sRGB y un espacio de color RGB lineal. Entonces, en realidad podrías estar tratando de optimizar eso, específicamente. No lo sé, y es por eso que esto podría no responder su pregunta.

Si este es el caso, intente usar una tabla de búsqueda. Algo como:

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

Si está utilizando datos de 16 bits, cámbielos según corresponda. Haría la tabla de 16 bits de todos modos para que pueda oscilar el resultado si es necesario cuando se trabaja con datos de 8 bits. This obviously won''t work very well if your data is floating point to begin with -- but it doesn''t really make sense to store sRGB data in floating point, so you might as well convert to 16-bit / 8-bit first and then do the conversion from linear to sRGB.

(The reason sRGB doesn''t make sense as floating point is that HDR should be linear, and sRGB is only convenient for storing on disk or displaying on screen, but not convenient for manipulation.)

Ian Stephenson escribió este código que, según él, supera a pow() . Él describe la idea de la siguiente manera:

Pow básicamente se implementa utilizando log''s: pow(a,b)=x(logx(a)*b) . entonces necesitamos un registro rápido y un exponente rápido, no importa qué es x, entonces usamos 2. El truco es que un número de coma flotante ya está en un formato de estilo de registro:

a=M*2E

Tomar el registro de ambos lados da:

log2(a)=log2(M)+E

o más simplemente:

log2(a)~=E

En otras palabras, si tomamos la representación en coma flotante de un número y extraemos el exponente, tenemos algo que es un buen punto de partida para su registro. Resulta que cuando hacemos esto masajeando los patrones de bits, el Mantissa termina dando una buena aproximación al error, y funciona bastante bien.

Esto debería ser lo suficientemente bueno para cálculos simples de iluminación, pero si necesita algo mejor, puede extraer el Mantissa y usarlo para calcular un factor de corrección cuadrático que es bastante preciso.

Otra respuesta porque esto es muy diferente de mi respuesta anterior, y esto es increíblemente rápido. El error relativo es 3e-8. ¿Quieres más precisión? Agregue un par de términos más de Chebychev. Lo mejor es mantener el orden impar ya que esto hace que haya una pequeña discontinuidad entre 2 ^ n-epsilon y 2 ^ n + epsilon.

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

Adición: ¿Qué está pasando aquí?
Por solicitud, a continuación se explica cómo funciona el código anterior.

Visión de conjunto
El código anterior define dos funciones, double pow512norm (double x) y double pow512 (double x) . Este último es el punto de entrada al conjunto; esta es la función que el código de usuario debe llamar para calcular x ^ (5/12). La función pow512norm(x) usa polinomios de Chebyshev para aproximar x ^ (5/12), pero solo para x en el rango [1,2]. (Use pow512norm(x) para valores de x fuera de ese rango y el resultado será basura).

La función pow512(x) divide la x entrante en un par (double s, int n) tal que x = s * 2^n y tal que 1≤ s <2. Una partición adicional de n en (int q, unsigned int r) tal que n = 12*q + r y r es menor que 12 me permite dividir el problema de encontrar x ^ (5/12) en partes:

  1. x^(5/12)=(s^(5/12))*((2^n)^(5/12)) vía (u v) ^ a = (u ^ a) (v ^ a) para positivo u, vy real a.
  2. s^(5/12) se calcula a través de pow512norm(s) .
  3. (2^n)^(5/12)=(2^(12*q+r))^(5/12) través de la sustitución.
  4. 2^(12*q+r)=(2^(12*q))*(2^r) través de u^(a+b)=(u^a)*(u^b) para positivo u, real a, b.
  5. (2^(12*q+r))^(5/12)=(2^(5*q))*((2^r)^(5/12)) mediante algunas manipulaciones adicionales.
  6. (2^r)^(5/12) se calcula mediante la tabla de búsqueda pow2_512 .
  7. Calcula pow512norm(s)*pow2_512[qr.rem] y ya casi llegamos. Aquí qr.rem es el valor de r calculado en el paso 3 anterior. Todo lo que se necesita es multiplicar esto por 2^(5*q) para obtener el resultado deseado.
  8. Eso es exactamente lo que hace la función de biblioteca matemática ldexp .

Aproximación de la función
El objetivo aquí es llegar a una aproximación fácilmente computable de f (x) = x ^ (5/12) que sea "lo suficientemente buena" para el problema en cuestión. Nuestra aproximación debería estar cerca de f (x) en algún sentido. Pregunta retórica: ¿Qué significa "cerca de"? Dos interpretaciones competitivas minimizan el error cuadrático medio frente a la minimización del error absoluto máximo.

Utilizaré una analogía del mercado de valores para describir la diferencia entre estos. Supongamos que quiere ahorrar para su eventual jubilación. Si tiene veintitantos años, lo mejor que puede hacer es invertir en acciones o en fondos del mercado de acciones. Esto se debe a que, durante un lapso de tiempo suficientemente largo, el mercado de valores en general supera a cualquier otro esquema de inversión. Sin embargo, todos hemos visto momentos en los que invertir dinero en acciones es algo muy malo. Si tiene entre 50 y 60 años (o 40 si quiere jubilarse joven), necesita invertir un poco más conservadoramente. Esos descensos pueden afectar su cartera de jubilación.

Volver a la aproximación de funciones: como consumidor de cierta aproximación, generalmente está preocupado por el peor de los casos, más que por el rendimiento "en promedio". Utilice una aproximación construida para obtener el mejor rendimiento "en promedio" (por ejemplo, mínimos cuadrados) y la ley de Murphy establece que su programa empleará mucho tiempo utilizando la aproximación exactamente donde el rendimiento es mucho peor que el promedio. Lo que quiere es una aproximación minimax, algo que minimice el error absoluto máximo sobre algún dominio. Una buena biblioteca de matemáticas tendrá un enfoque minimax en lugar de un enfoque de mínimos cuadrados, ya que esto permite a los autores de la biblioteca matemática dar un rendimiento garantizado de su biblioteca.

Las bibliotecas matemáticas generalmente usan un polinomio o un polinomio racional para aproximar alguna función f (x) sobre algún dominio a≤x≤b. Suponga que la función f (x) es analítica sobre este dominio y desea aproximar la función por un polinomio p (x) de grado N. Para un grado dado N existe un polinomio mágico único (p) x tal que p ( x) -f (x) tiene N + 2 extremos sobre [a, b] y tal que los valores absolutos de estos extremos N + 2 son todos iguales entre sí. Encontrar este polinomio mágico p (x) es el santo grial de los aproximadores de funciones.

No encontré ese Santo Grial para ti. En su lugar, utilicé una aproximación de Chebyshev. Los polinomios de Chebyshev del primer tipo son un conjunto de polinomios ortogonales (pero no ortonormales) con algunas características muy agradables cuando se trata de la aproximación de funciones. La aproximación de Chebyshev a menudo está muy cerca de ese polinomio mágico p (x). (De hecho, el algoritmo de intercambio de Remez que sí encuentra ese polinomio de santo grial generalmente comienza con una aproximación de Chebyshev).

pow512norm (x)
Esta función utiliza la aproximación de Chebyshev para encontrar un polinomio p * (x) que se aproxima a x ^ (5/12). Aquí estoy usando p * (x) para distinguir esta aproximación de Chebyshev del polinomio mágico p (x) descrito anteriormente. La aproximación de Chebyshev p * (x) es fácil de encontrar; encontrar p (x) es un oso. La aproximación de Chebyshev p * (x) es sum_i Cn [i] * Tn (i, x), donde Cn [i] son ​​los coeficientes de Chebyshev y Tn (i, x) son los polinomios de Chebyshev evaluados en x.

Utilicé Wolfram alpha para encontrar los coeficientes de Chebyshev Cn para mí. Por ejemplo, esto calcula Cn[1] . El primer cuadro después del cuadro de entrada tiene la respuesta deseada, 0.166658 en este caso. No son tantos los dígitos como me gustaría. Haga clic en "más dígitos" y listo, obtendrá muchos más dígitos. Wolfram alpha es gratis; hay un límite en la cantidad de cálculo que hará. Llega a ese límite en términos de orden superior. (Si compra o tiene acceso a mathematica, podrá calcular esos coeficientes de orden alto con un alto grado de precisión).

Los polinomios de Chebyshev Tn (x) se calculan en la matriz Tn . Más allá de dar algo muy cercano al polinomio mágico p (x), otra razón para usar la aproximación de Chebyshev es que los valores de esos polinomios de Chebyshev se calculan fácilmente: Comience con Tn[0]=1 y Tn[1]=x , y luego iterativamente calcule Tn[i]=2*x*Tn[i-1] - Tn[i-2] . (Utilicé ''ii'' como la variable de índice en lugar de ''i'' en mi código. Nunca uso ''i'' como nombre de variable. ¿Cuántas palabras en inglés tienen una ''i'' en la palabra? ¿Cuántas tienen dos? ''i'' consecutivos?)

pow512 (x)
pow512 es la función que debería estar llamando el código de usuario. Ya describí los conceptos básicos de esta función más arriba. Algunos detalles más: la función de biblioteca matemática frexp(x) devuelve el significado y el exponente iexp para la entrada x . ( pow512norm menor: quiero s entre 1 y 2 para usar con pow512norm pero frexp devuelve un valor entre 0.5 y 1.) La función de la biblioteca matemática div devuelve el cociente y el resto para la división de enteros en un solo swop. Finalmente, utilizo la función de biblioteca matemática ldexp para juntar las tres partes para formar la respuesta final.

Binomial series does account for a constant exponent, but you will be able to use it only if you can normalize all your input to the range [1,2). (Note that it computes (1+x)^a). You''ll have to do some analysis to decide how many terms you need for your desired accuracy.

For exponents of 2.4, you could either make a lookup table for all your 2.4 values and lirp or perhaps higher-order function to fill in the in-betweem values if the table wasn''t accurate enough (basically a huge log table.)

Or, value squared * value to the 2/5s which could take the initial square value from the first half of the function and then 5th root it. For the 5th root, you could Newton it or do some other fast approximator, though honestly once you get to this point, your probably better off just doing the exp and log functions with the appropriate abbreviated series functions yourself.

I shall answer the question you really wanted to ask, which is how to do fast sRGB <-> linear RGB conversion. To do this precisely and efficiently we can use polynomial approximations. The following polynomial approximations have been generated with sollya, and have a worst case relative error of 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); }

And the sollya input used to generate the polynomials:

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

So traditionally the powf(x, p) = x^p is solved by rewriting x as x=2^(log2(x)) making powf(x,p) = 2^(p*log2(x)) , which transforms the problem into two approximations exp2() & log2() . This has the advantage of working with larger powers p , however the downside is that this is not the optimal solution for a constant power p and over a specified input bound 0 ≤ x ≤ 1 .

When the power p > 1 , the answer is a trivial minimax polynomial over the bound 0 ≤ x ≤ 1 , which is the case for p = 12/5 = 2.4 as can be seen below:

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

However when p < 1 the minimax approximation over the bound 0 ≤ x ≤ 1 does not appropriately converge to the desired accuracy. One option [not really] is to rewrite the problem y=x^p=x^(p+m)/x^m where m=1,2,3 is a positive integer, making the new power approximation p > 1 but this introduces division which is inherently slower.

There''s however another option which is to decompose the input x as its floating point exponent and mantissa form:

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)

The minimax approximation of mx^(5/12) over 1 ≤ mx < 2 now converges much faster than before, without division, but requires 12 point LUT for the 2^(k/12) . The code is below:

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

The following is an idea you can use with any of the fast calculation methods. Whether it helps things go faster depends on how your data arrives. You can use the fact that if you know x and pow(x, n) , you can use the rate of change of the power to compute a reasonable approximation of pow(x + delta, n) for small delta , with a single multiply and add (more or less). If successive values you feed your power functions are close enough together, this would amortize the full cost of the accurate calculation over multiple function calls. Note that you don''t need an extra pow calculation to get the derivative. You could extend this to use the second derivative so you can use a quadratic, which would increase the delta you could use and still get the same accuracy.