¿Por qué GCC no optimiza a * a * a * a * a * a to (a * a * a) * (a * a * a)?

Estoy haciendo una optimización numérica en una aplicación científica. Una cosa que noté es que GCC optimizará la llamada pow(a,2)compilándola a*a, pero la llamada pow(a,6)no está optimizada y realmente llamará a la función de biblioteca pow, lo que ralentiza enormemente el rendimiento. (En contraste, el compilador Intel C ++ , ejecutable icc, eliminará la llamada a la biblioteca pow(a,6)).

Lo que me interesa es que cuando lo reemplacé pow(a,6)con a*a*a*a*a*aGCC 4.5.1 y las opciones " -O3 -lm -funroll-loops -msse4", utiliza 5 mulsdinstrucciones:

movapd  %xmm14, %xmm13
mulsd   %xmm14, %xmm13
mulsd   %xmm14, %xmm13
mulsd   %xmm14, %xmm13
mulsd   %xmm14, %xmm13
mulsd   %xmm14, %xmm13

mientras que si escribo (a*a*a)*(a*a*a), producirá

movapd  %xmm14, %xmm13
mulsd   %xmm14, %xmm13
mulsd   %xmm14, %xmm13
mulsd   %xmm13, %xmm13

lo que reduce el número de instrucciones de multiplicación a 3. icctiene un comportamiento similar.

¿Por qué los compiladores no reconocen este truco de optimización?

Porque Floating Point Math no es asociativo . La forma en que agrupa los operandos en la multiplicación de coma flotante tiene un efecto en la precisión numérica de la respuesta.

Como resultado, la mayoría de los compiladores son muy conservadores al reordenar los cálculos de coma flotante a menos que puedan estar seguros de que la respuesta seguirá siendo la misma, o a menos que les diga que no le importa la precisión numérica. Por ejemplo: la -fassociative-mathopción de gcc que le permite a gcc reasociar operaciones de coma flotante, o incluso la -ffast-mathopción que permite compensaciones aún más agresivas de precisión contra velocidad.

Lambdageek señala correctamente que debido a que la asociatividad no se cumple para los números de coma flotante, la "optimización" dea*a*a*a*a*ato(a*a*a)*(a*a*a)puede cambiar el valor. Es por eso que C99 no lo permite (a menos que el usuario lo permita específicamente, a través del indicador del compilador o pragma). En general, se supone que el programador escribió lo que hizo por una razón, y el compilador debería respetar eso. Si quieres(a*a*a)*(a*a*a), escribe eso.

Sin embargo, puede ser un dolor de escribir; ¿Por qué el compilador no puede hacer [lo que consideras que es] lo correcto cuando lo usas pow(a,6)? Porque sería lo incorrecto hacer. En una plataforma con una buena biblioteca de matemáticas, pow(a,6)es significativamente más preciso que sea a*a*a*a*a*ao (a*a*a)*(a*a*a). Solo para proporcionar algunos datos, realicé un pequeño experimento en mi Mac Pro, midiendo el peor error al evaluar un ^ 6 para todos los números flotantes de precisión simple entre [1,2):

worst relative error using    powf(a, 6.f): 5.96e-08
worst relative error using (a*a*a)*(a*a*a): 2.94e-07
worst relative error using     a*a*a*a*a*a: 2.58e-07

Usar en powlugar de un árbol de multiplicación reduce el error limitado por un factor de 4 . Los compiladores no deben (y generalmente no lo hacen) "optimizaciones" que aumenten el error a menos que el usuario lo autorice (por ejemplo, a través de -ffast-math).

Tenga en cuenta que GCC proporciona __builtin_powi(x,n)como alternativa a pow( ), que debería generar un árbol de multiplicación en línea. Úselo si desea cambiar la precisión por el rendimiento, pero no desea habilitar las matemáticas rápidas.

Otro caso similar: la mayoría de los compiladores no optimiza a + b + c + da (a + b) + (c + d)(esto es una optimización desde la segunda expresión se puede pipeline mejor) y evaluarlo como dado (es decir, como (((a + b) + c) + d)). Esto también se debe a casos de esquina:

float a = 1e35, b = 1e-5, c = -1e35, d = 1e-5;
printf("%e %e\n", a + b + c + d, (a + b) + (c + d));

Esto salidas 1.000000e-05 0.000000e+00

Fortran (diseñado para computación científica) tiene un operador de potencia incorporado, y que yo sepa, los compiladores de Fortran comúnmente optimizarán la elevación a potencias enteras de manera similar a lo que usted describe. Desafortunadamente, C / C ++ no tiene un operador de energía, solo la función de biblioteca pow(). Esto no evita que los compiladores inteligentes lo traten powespecialmente y lo computen de una manera más rápida para casos especiales, pero parece que lo hacen con menos frecuencia ...

Hace algunos años, estaba tratando de hacer más conveniente el cálculo de potencias enteras de manera óptima, y ​​se me ocurrió lo siguiente. Sin embargo, es C ++, no C, y aún depende de que el compilador sea algo inteligente sobre cómo optimizar / alinear las cosas. De todos modos, espero que les resulte útil en la práctica:

template<unsigned N> struct power_impl;

template<unsigned N> struct power_impl {
    template<typename T>
    static T calc(const T &x) {
        if (N%2 == 0)
            return power_impl<N/2>::calc(x*x);
        else if (N%3 == 0)
            return power_impl<N/3>::calc(x*x*x);
        return power_impl<N-1>::calc(x)*x;
    }
};

template<> struct power_impl<0> {
    template<typename T>
    static T calc(const T &) { return 1; }
};

template<unsigned N, typename T>
inline T power(const T &x) {
    return power_impl<N>::calc(x);
}

_{Aclaración para los curiosos: esto no encuentra la forma óptima de calcular las potencias, pero dado que encontrar la solución óptima es un problema NP-completo y esto solo vale la pena hacerlo para las potencias pequeñas de todos modos (en lugar de usar pow), no hay razón para preocuparse Con el detalle.}

Entonces solo úsalo como power<6>(a).

Esto facilita la escritura de poderes (no es necesario deletrear 6 as con parens), y le permite tener este tipo de optimización sin -ffast-mathtener que depender de la precisión, como la suma compensada (un ejemplo donde el orden de las operaciones es esencial) .

Probablemente también pueda olvidar que se trata de C ++ y simplemente usarlo en el programa C (si se compila con un compilador de C ++).

Espero que esto pueda ser útil.

EDITAR:

Esto es lo que obtengo de mi compilador:

para a*a*a*a*a*a,

    movapd  %xmm1, %xmm0
    mulsd   %xmm1, %xmm0
    mulsd   %xmm1, %xmm0
    mulsd   %xmm1, %xmm0
    mulsd   %xmm1, %xmm0
    mulsd   %xmm1, %xmm0

para (a*a*a)*(a*a*a),

    movapd  %xmm1, %xmm0
    mulsd   %xmm1, %xmm0
    mulsd   %xmm1, %xmm0
    mulsd   %xmm0, %xmm0

para power<6>(a),

    mulsd   %xmm0, %xmm0
    movapd  %xmm0, %xmm1
    mulsd   %xmm0, %xmm1
    mulsd   %xmm0, %xmm1

GCC realmente se optimiza a*a*a*a*a*apara (a*a*a)*(a*a*a)cuando a es un número entero. Intenté con este comando:

$ echo 'int f(int x) { return x*x*x*x*x*x; }' | gcc -o - -O2 -S -masm=intel -x c -

Hay muchas banderas de gcc pero nada lujoso. Significan: leer de stdin; utilizar el nivel de optimización de O2; salida del listado del lenguaje ensamblador en lugar de un binario; la lista debe usar la sintaxis del lenguaje ensamblador Intel; la entrada está en lenguaje C (generalmente el idioma se infiere de la extensión del archivo de entrada, pero no hay extensión de archivo cuando se lee desde stdin); y escribe a stdout.

Aquí está la parte importante de la salida. Lo he anotado con algunos comentarios que indican lo que está sucediendo en el lenguaje ensamblador:

; x is in edi to begin with.  eax will be used as a temporary register.
mov  eax, edi  ; temp = x
imul eax, edi  ; temp = x * temp
imul eax, edi  ; temp = x * temp
imul eax, eax  ; temp = temp * temp

Estoy usando el sistema GCC en Linux Mint 16 Petra, un derivado de Ubuntu. Aquí está la versión de gcc:

$ gcc --version
gcc (Ubuntu/Linaro 4.8.1-10ubuntu9) 4.8.1

Como han señalado otros carteles, esta opción no es posible en coma flotante, porque la aritmética de coma flotante no es asociativa.

Porque un número de coma flotante de 32 bits, como 1.024, no es 1.024. En una computadora, 1.024 es un intervalo: de (1.024-e) a (1.024 + e), donde "e" representa un error. Algunas personas no se dan cuenta de esto y también creen que * en a * a representa la multiplicación de números de precisión arbitraria sin que haya ningún error asociado a esos números. La razón por la cual algunas personas no se dan cuenta de esto es quizás los cálculos matemáticos que ejercieron en las escuelas primarias: trabajar solo con números ideales sin errores adjuntos y creer que está bien simplemente ignorar "e" mientras se realiza la multiplicación. No ven la "e" implícita en "float a = 1.2", "a * a * a" y códigos C similares.

Si la mayoría de los programadores reconocen (y pueden ejecutar) la idea de que la expresión C a * a * a * a * a * a en realidad no funciona con números ideales, el compilador GCC sería GRATIS para optimizar "a * a * a * a * a * a "en decir" t = (a * a); t * t * t "que requiere un número menor de multiplicaciones. Pero desafortunadamente, el compilador de GCC no sabe si el programador que escribe el código piensa que "a" es un número con o sin error. Y así, GCC solo hará lo que parece el código fuente, porque eso es lo que GCC ve a simple vista.

... una vez que sepa qué tipo de programador que se encuentre, puede utilizar el interruptor "-ffast-matemáticas" para decirle a gcc que "Hey, GCC, sé lo que estoy haciendo!". Esto permitirá que GCC convierta a * a * a * a * a * a en un texto diferente; se ve diferente de a * a * a * a * a * a, pero aún calcula un número dentro del intervalo de error de a * a * a * a * a * a. Esto está bien, ya que ya sabe que está trabajando con intervalos, no con números ideales.

Ningún póster ha mencionado todavía la contracción de las expresiones flotantes (estándar ISO C, 6.5p8 y 7.12.2). Si el FP_CONTRACTpragma se establece en ON, el compilador puede considerar una expresión a*a*a*a*a*acomo una sola operación, como si se evaluara exactamente con un solo redondeo. Por ejemplo, un compilador puede reemplazarlo por una función de potencia interna que es más rápida y más precisa. Esto es particularmente interesante ya que el comportamiento está parcialmente controlado por el programador directamente en el código fuente, mientras que las opciones del compilador proporcionadas por el usuario final a veces se pueden usar incorrectamente.

El estado predeterminado del FP_CONTRACTpragma está definido por la implementación, de modo que un compilador puede realizar tales optimizaciones de manera predeterminada. Por lo tanto, el código portátil que debe seguir estrictamente las reglas IEEE 754 debe establecerlo explícitamente OFF.

Si un compilador no admite este pragma, debe ser conservador evitando dicha optimización, en caso de que el desarrollador haya elegido configurarlo OFF.

GCC no admite este pragma, pero con las opciones predeterminadas, asume que es así ON; por lo tanto, para objetivos con un FMA de hardware, si se quiere evitar la transformación a*b+ca fma (a, b, c), se debe proporcionar una opción como -ffp-contract=off(establecer explícitamente el pragma en OFF) o -std=c99(para indicarle a GCC que se ajuste a algunos Versión estándar C, aquí C99, por lo tanto, siga el párrafo anterior). En el pasado, la última opción no impedía la transformación, lo que significa que GCC no se ajustaba a este punto: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845

Como Lambdageek señaló, la multiplicación por flotación no es asociativa y se puede obtener menos precisión, pero también cuando se obtiene una mayor precisión se puede argumentar en contra de la optimización, porque se desea una aplicación determinista. Por ejemplo, en la simulación de juego cliente / servidor, donde cada cliente tiene que simular el mismo mundo en el que desea que los cálculos de coma flotante sean deterministas.

Las funciones de la biblioteca como "pow" generalmente están cuidadosamente diseñadas para producir el mínimo error posible (en caso genérico). Esto generalmente se logra aproximando funciones con splines (según el comentario de Pascal, la implementación más común parece estar usando el algoritmo Remez )

fundamentalmente la siguiente operación:

pow(x,y);

tiene un error inherente de aproximadamente la misma magnitud que el error en cualquier multiplicación o división .

Mientras que la siguiente operación:

float a=someValue;
float b=a*a*a*a*a*a;

tiene un error inherente que es mayor que 5 veces el error de una sola multiplicación o división (porque está combinando 5 multiplicaciones).

El compilador debe tener mucho cuidado con el tipo de optimización que está haciendo:

1. Si la optimización pow(a,6)de a*a*a*a*a*aque puede mejorar el rendimiento, pero reducir drásticamente la exactitud de los números de punto flotante.
2. si la optimización a*a*a*a*a*a a pow(a,6)que en realidad puede reducir la precisión porque "a" era algún valor especial que permite la multiplicación sin error (una potencia de 2 o un número entero pequeño)
3. Si la optimización pow(a,6)de (a*a*a)*(a*a*a)o (a*a)*(a*a)*(a*a)que todavía puede ser una pérdida de precisión en comparación con powla función.

En general, usted sabe que para valores arbitrarios de coma flotante "pow" tiene mejor precisión que cualquier función que eventualmente podría escribir, pero en algunos casos especiales, las multiplicaciones múltiples pueden tener una mejor precisión y rendimiento, depende del desarrollador elegir qué es más apropiado, eventualmente comentando el código para que nadie más "optimice" ese código.

Lo único que tiene sentido (opinión personal, y aparentemente una elección en GCC sin ninguna optimización particular o indicador de compilación) para optimizar debería ser reemplazar "pow (a, 2)" con "a * a". Eso sería lo único sensato que un vendedor de compiladores debería hacer.

No hubiera esperado que este caso se optimizara en absoluto. No puede ser muy frecuente que una expresión contenga subexpresiones que pueden reagruparse para eliminar operaciones completas. Esperaría que los escritores de compiladores inviertan su tiempo en áreas que tendrían más probabilidades de generar mejoras notables, en lugar de cubrir un caso marginal que rara vez se encuentra.

Me sorprendió saber de las otras respuestas que esta expresión podría optimizarse con los modificadores de compilador adecuados. O la optimización es trivial, o es un caso extremo de una optimización mucho más común, o los escritores del compilador fueron extremadamente minuciosos.

No hay nada de malo en proporcionar pistas al compilador como lo ha hecho aquí. Es una parte normal y esperada del proceso de microoptimización reorganizar las declaraciones y expresiones para ver qué diferencias traerán.

Si bien el compilador puede estar justificado al considerar las dos expresiones para entregar resultados inconsistentes (sin los modificadores adecuados), no es necesario que esté sujeto a esa restricción. La diferencia será increíblemente pequeña, tanto que si la diferencia es importante para usted, no debería usar la aritmética de coma flotante estándar en primer lugar.

Ya hay algunas buenas respuestas a esta pregunta, pero en aras de la exhaustividad, quería señalar que la sección correspondiente del estándar C es 5.1.2.2.3 / 15 (que es lo mismo que la sección 1.9 / 9 en el C ++ 11 estándar). Esta sección establece que los operadores solo pueden reagruparse si son realmente asociativos o conmutativos.

En realidad, gcc puede hacer esta optimización, incluso para números de punto flotante. Por ejemplo,

double foo(double a) {
  return a*a*a*a*a*a;
}

se convierte

foo(double):
    mulsd   %xmm0, %xmm0
    movapd  %xmm0, %xmm1
    mulsd   %xmm0, %xmm1
    mulsd   %xmm1, %xmm0
    ret

con -O -funsafe-math-optimizations. Sin embargo, este reordenamiento viola IEEE-754, por lo que requiere la bandera.

Los enteros firmados, como señaló Peter Cordes en un comentario, pueden hacer esta optimización sin, -funsafe-math-optimizationsya que se mantiene exactamente cuando no hay desbordamiento y si hay desbordamiento, obtienes un comportamiento indefinido. Entonces obtienes

foo(long):
    movq    %rdi, %rax
    imulq   %rdi, %rax
    imulq   %rdi, %rax
    imulq   %rax, %rax
    ret

con sólo -O. Para enteros sin signo, es aún más fácil ya que funcionan con potencias mod de 2 y, por lo tanto, se pueden reordenar libremente incluso ante el desbordamiento.