¿Optimizaciones para pow () con un exponente const no entero?

Otra respuesta porque esto es muy diferente de mi respuesta anterior, y esto es muy rápido. El error relativo es 3e-8. Quieres más precisión? Añade un par de términos más de Chebychev. Es mejor mantener el orden impar ya que esto hace que para 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);
}

Anexo: ¿Qué está pasando aquí?
Por solicitud, lo siguiente explica cómo funciona el código anterior.

Sinopsis
El código anterior define dos funciones, double pow512norm (double x) y double pow512 (double x). Este último es el punto de entrada a la suite; esta es la función que el código de usuario debe llamar para calcular x^(5/12). La función pow512norm(x) utiliza 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 el entrante x en un par (double s, int n) tal que x = s * 2^n y tal que 1≤sn 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)) via (u v)^a=(u^a) (v^a) para u positivo,v y a real.
2. s^(5/12) se calcula a través de pow512norm(s).
3. (2^n)^(5/12)=(2^(12*q+r))^(5/12) mediante sustitución.
4. 2^(12*q+r)=(2^(12*q))*(2^r) vía u^(a+b)=(u^a)*(u^b) para u positivo, a real, b.
5. (2^(12*q+r))^(5/12)=(2^(5*q))*((2^r)^(5/12)) a través de algunas manipulaciones más.
6. (2^r)^(5/12) se calcula mediante la tabla de búsqueda pow2_512.
7. Calcular pow512norm(s)*pow2_512[qr.rem] y estamos casi allí. Aqui qr.rem es el valor r calculado en el paso 3 anterior. Todo lo que se necesita es multiplicar esto por 2^(5*q) para producir el resultado deseado.
8. Eso es exactamente lo que hace la función de biblioteca matemática ldexp.

Aproximación De Funciones
El objetivo aquí es llegar a una aproximación fácilmente computable de f(x) = x^(5/12) que es 'lo suficientemente bueno' para el problema en cuestión. Nuestra aproximación debe estar cerca de f (x) en algún sentido. Pregunta retórica: ¿Qué hace 'cerca de' ¿mala? Dos interpretaciones competidoras son minimizar el error cuadrado medio versus minimizar el error absoluto máximo.

Voy a utilizar una analogía del mercado de valores para describir la diferencia entre estos. Supongamos que desea ahorrar para su eventual jubilación. Si usted está en sus veinte años, lo mejor que puede hacer es invertir en acciones o fondos del mercado de valores. Esto se debe a que durante un período de tiempo lo suficientemente largo, el mercado de valores en promedio supera a cualquier otro esquema de inversión. Sin embargo, todos hemos visto veces cuando poner dinero en acciones es una cosa muy mala que hacer. Si usted está en sus cincuenta o sesenta (o cuarenta años si desea retirarse joven) que necesita para invertir un poco más conservador. Esas bajadas pueden tener en tu cartera de jubilación.

Volver a la aproximación de funciones: Como consumidor de alguna aproximación, normalmente le preocupa el error en el peor de los casos en lugar del rendimiento "en promedio". Utilice alguna aproximación construida para dar el mejor rendimiento " en promedio " (por ejemplo, mínimos cuadrados) y la ley de Murphy dicta que su programa pasará mucho tiempo usando la aproximación exactamente donde el rendimiento es mucho peor que el promedio. Lo que quieres es una aproximación minimax, algo que minimice el error absoluto máximo sobre algún dominio. Una buena biblioteca de matemáticas tomará un enfoque minimax en lugar de un enfoque de mínimos cuadrados porque esto permite a los autores de la biblioteca de matemáticas dar un rendimiento garantizado de su biblioteca.

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

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 ortogonal (pero no ortonormal) de polinomios 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 Remez que encuentra ese santo grial el polinomio comienza típicamente con una aproximación de Chebyshev.)

Pow512norm(x)
Esta función utiliza la aproximación de Chebyshev para encontrar algún 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 los Cn [i] son los Chebyshev los coeficientes y Tn (i, x) son los polinomios de Chebyshev evaluados en x.

Usé 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 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 computación que hará. Llega a ese límite en términos de orden más altos. (Si usted compra o tiene acceso a mathematica, podrá calcular esos coeficientes de alto orden con un alto grado de precisión.)

Los polinomios de Chebyshev Tn(x) se calculan en el array 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 calcule iterativamente Tn[i]=2*x*Tn[i-1] - Tn[i-2]. (Yo usé 'ii' como variable de índice en lugar de' i ' en mi código. Nunca uso 'i' como nombre de variable. ¿Cuántas palabras en el idioma inglés tienen una " i " en la palabra? ¿Cuántos tienen dos i consecutivas?)

Pow512 (x)
pow512 es la función que el código de usuario debe llamar. Ya he descrito los conceptos básicos de esta función anteriormente. Algunos detalles más: La función math library frexp(x) devuelve el significante s y el exponente iexp para la entrada x. (Tema menor: quiero s entre 1 y 2 para usar con pow512norm pero frexp devuelve un valor entre 0.5 y 1.) La función math library div devuelve el cociente y el resto para la división de enteros en un pie de página hinchado. Finalmente, utilizo la función math library ldexp para juntar las tres partes para formar la respuesta final.

Ian Stephenson escribió este código que él afirma supera pow(). Él describe la idea de la siguiente manera:

Pow se implementa básicamente usando log's: pow(a,b)=x(logx(a)*b). así que necesita un registro rápido y exponente rápido-it no importa lo que x es por lo que usamos 2. El truco es que un punto flotante el número ya está en un estilo de registro formato:

a=M*2E

Tomando el registro de ambos lados da:

log2(a)=log2(M)+E

O más simplemente:

log2(a)~=E

En otros palabras si tomamos la flotación representación puntual de un número, y extraer el Exponente que tenemos algo que es un buen punto de partida como su tronco. Resulta que cuando haga esto masajeando los patrones de bits, la Mantisa termina dando un buen aproximación al error, y funciona bastante bien.

Esto debería ser lo suficientemente bueno para simple cálculos de iluminación, pero si necesita algo mejor, a continuación, se puede extraer la Mantisa, y usa eso para calcular un factor de corrección cuadrática lo cual es bastante exacto.

En primer lugar, el uso de flotadores no va a comprar mucho en la mayoría de las máquinas de hoy en día. De hecho, los dobles pueden ser más rápidos. Su poder, 1.0 / 2.4, es 5/12 o 1/3*(1 + 1/4). A pesar de que esto está llamando cbrt (una vez) y sqrt (dos veces!) sigue siendo el doble de rápido que usando 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 a su pregunta.

El 2.4f y 1/2.4f me hacen sospechar mucho, porque esos son exactamente los poderes utilizados para convertir entre sRGB y un espacio de color RGB lineal. Por lo que en realidad podría estar tratando de optimizar que, específicamente. No lo sé, por lo que esto podría no responder a tu 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, cambie según corresponda. Yo haría la tabla de 16 bits de todos modos para que pueda tramar el resultado si es necesario cuando se trabaja con datos de 8 bits. Esto obviamente no funcionará muy bien si sus datos son de punto flotante para empezar -- pero realmente no tiene sentido almacenar datos sRGB en punto flotante, por lo que también podría convertir a 16 bits / 8 bits primero y luego hacer la conversión de lineal a sRGB.

(La razón por la que sRGB no tiene sentido como punto flotante es que HDR debe ser lineal, y sRGB solo es conveniente para almacenar en disco o visualización en pantalla, pero no conveniente para la manipulación.)

La serie binomial da cuenta de un exponente constante, pero solo podrá usarlo si puede normalizar toda su entrada al rango [1,2). (Tenga en cuenta que calcula (1 + x)^a). Tendrá que hacer un análisis para decidir cuántos términos necesita para su precisión deseada.

Para exponentes de 2.4, puede hacer una tabla de búsqueda para todos sus valores de 2.4 y lirp o tal vez una función de orden superior para completar los valores intermedios si la tabla no era lo suficientemente precisa (básicamente una tabla de registro enorme.)

O, valor cuadrado * valor a los 2/5s que podría tomar el valor cuadrado inicial de la primera mitad de la función y luego la 5a raíz. Para la 5a raíz, podrías Newton o hacer algún otro aproximador rápido, aunque honestamente una vez que llegues a este punto, probablemente sea mejor que solo haga las funciones exp y log con las funciones de serie abreviada adecuadas usted mismo.

La siguiente es una idea que puede usar con cualquiera de los métodos de cálculo rápidos. Si ayuda a que las cosas vayan más rápido depende de cómo lleguen sus datos. Puedes usar el hecho de que si sabes x y pow(x, n), puedes usar la tasa de cambio de la potencia para calcular una aproximación razonable de pow(x + delta, n) para delta pequeños, con una sola multiplicación y suma (más o menos). Si los valores sucesivos que alimenta sus funciones de potencia están lo suficientemente cerca entre sí, esto amortizaría el costo total de la precisión cálculo sobre múltiples llamadas a funciones. Tenga en cuenta que no necesita un cálculo pow extra para obtener la derivada. Usted podría extender esto para usar la segunda derivada para que pueda usar una cuadrática, lo que aumentaría el delta que podría usar y aún así obtener la misma precisión.

Voy a responder a la pregunta que realmente quería preguntar, que es cómo hacer una rápida conversión RGB lineal sRGB. Para hacer esto con precisión y eficiencia podemos usar aproximaciones polinómicas. Las siguientes aproximaciones polinómicas se han generado con sollya, y tienen un error relativo en el peor de los casos de 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);
}

Y la entrada sollya utilizada para generar los polinomios:

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

Así que tradicionalmente el powf(x, p) = x^p se resuelve reescribiendo x como x=2^(log2(x)) haciendo powf(x,p) = 2^(p*log2(x)), lo que transforma el problema en dos aproximaciones exp2() & log2(). Esto tiene la ventaja de trabajar con potencias mayores p, sin embargo, la desventaja es que esta no es la solución óptima para una potencia constante p y sobre una entrada determinada 0 ≤ x ≤ 1.

Cuando la potencia p > 1, la respuesta es un polinomio minimax trivial sobre el límite 0 ≤ x ≤ 1, que es el caso de p = 12/5 = 2.4 como se puede ver abajo:

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

Sin embargo, cuando p < 1 la aproximación de minimax sobre el límite 0 ≤ x ≤ 1 no converge apropiadamente a la precisión deseada. Una opción [no realmente] es reescribir el problema y=x^p=x^(p+m)/x^m donde m=1,2,3 es un entero positivo, haciendo la nueva aproximación de potencia p > 1 pero esto introduce la división que es inherentemente más lenta.

Sin embargo, hay otra opción que es descomponer la entrada x como su exponente de punto flotante y forma mantisa:

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)

La aproximación minimax de mx^(5/12) sobre 1 ≤ mx < 2 ahora converge mucho más rápido que antes, sin división, pero requiere LUT de 12 puntos para el 2^(k/12). El código es el siguiente:

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