Por que o GCC não otimiza a * a * a * a * a * a para (a * a * a) * (a * a * a)?

Estou fazendo alguma otimização numérica em um aplicativo científico. Uma coisa que notei é que o GCC otimiza a chamada pow(a,2), compilando-a a*a, mas a chamada pow(a,6)não é otimizada e na verdade chama a função de biblioteca pow, o que diminui bastante o desempenho. (Por outro lado, o compilador Intel C ++ , executável icc, eliminará a chamada da biblioteca pow(a,6).)

O que eu estou curioso é que quando eu substituído pow(a,6)com a*a*a*a*a*autilizando GCC 4.5.1 e opções " -O3 -lm -funroll-loops -msse4", ele usa 5 mulsdinstruções:

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

enquanto se eu escrever (a*a*a)*(a*a*a), ele produzirá

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

que reduz o número de instruções de multiplicação para 3. icctem um comportamento semelhante.

Por que os compiladores não reconhecem esse truque de otimização?

Porque a matemática de ponto flutuante não é associativa . A maneira como você agrupa os operandos na multiplicação de ponto flutuante afeta a precisão numérica da resposta.

Como resultado, a maioria dos compiladores é muito conservadora ao reordenar os cálculos de ponto flutuante, a menos que eles possam ter certeza de que a resposta permanecerá a mesma ou que você diga a eles que não se importa com a precisão numérica. Por exemplo: a -fassociative-mathopção do gcc, que permite ao gcc reassociar operações de ponto flutuante, ou mesmo a -ffast-mathopção que permite compensações ainda mais agressivas de precisão em relação à velocidade.

Lambdageek indica corretamente que, como a associatividade não é válida para números de ponto flutuante, a "otimização" dea*a*a*a*a*apara(a*a*a)*(a*a*a)pode alterar o valor. É por isso que não é permitido pelo C99 (a menos que seja especificamente permitido pelo usuário, via sinalizador do compilador ou pragma). Geralmente, a suposição é de que o programador escreveu o que fez por uma razão, e o compilador deve respeitar isso. Se você quiser(a*a*a)*(a*a*a), escreva isso.

Isso pode ser uma dor de escrever; por que o compilador não pode fazer [o que você considera] a coisa certa quando você usa pow(a,6)? Porque seria a coisa errada a fazer. Em uma plataforma com uma boa biblioteca de matemática, pow(a,6)é significativamente mais preciso que um a*a*a*a*a*aou outro (a*a*a)*(a*a*a). Apenas para fornecer alguns dados, fiz um pequeno experimento no meu Mac Pro, medindo o pior erro na avaliação de um ^ 6 para todos os números flutuantes de precisão única 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 em powvez de uma árvore de multiplicação reduz o erro vinculado por um fator de 4 . Os compiladores não devem (e geralmente não fazem) "otimizações" que aumentam o erro, a menos que sejam licenciadas pelo usuário (por exemplo, via -ffast-math).

Observe que o GCC fornece __builtin_powi(x,n)uma alternativa para pow( ), o que deve gerar uma árvore de multiplicação em linha. Use isso se desejar trocar a precisão pelo desempenho, mas não deseja ativar a matemática rápida.

Um outro caso semelhante: a maioria dos compiladores não vai optimizar a + b + c + da (a + b) + (c + d)(esta é uma optimização desde a segunda expressão pode ser melhor pipeline) e avaliá-la como determinado (isto é, quanto (((a + b) + c) + d)). Isso também ocorre por causa dos casos de canto:

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

Isso gera 1.000000e-05 0.000000e+00

O Fortran (projetado para computação científica) possui um operador de energia embutido e, até onde eu sei, os compiladores do Fortran geralmente otimizam o aumento para potências inteiras de maneira semelhante à que você descreve. Infelizmente, o C / C ++ não possui um operador de energia, apenas a função de biblioteca pow(). Isso não impede que os compiladores inteligentes tratem powespecialmente e o computem de maneira mais rápida em casos especiais, mas parece que eles fazem isso com menos frequência ...

Alguns anos atrás, eu estava tentando torná-lo mais conveniente para calcular potências inteiras da maneira ideal, e criei o seguinte. É C ++, não C, e ainda depende do compilador ser um pouco inteligente sobre como otimizar / incorporar coisas. De qualquer forma, espero que você ache útil na prática:

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

_{Esclarecimento para os curiosos: isso não encontra a maneira ideal de calcular potências, mas como encontrar a solução ideal é um problema completo de NP e vale a pena fazer para potências pequenas de qualquer maneira (ao contrário de usar pow), não há razão para se preocupar com o detalhe.}

Então apenas use-o como power<6>(a).

Isso facilita a digitação de poderes (não é necessário especificar 6 as com parênteses) e permite que você tenha esse tipo de otimização sem -ffast-mathter algo dependente da precisão, como a soma compensada (um exemplo em que a ordem das operações é essencial) .

Você provavelmente também pode esquecer que este é C ++ e apenas usá-lo no programa C (se compilar com um compilador C ++).

Espero que isso possa ser útil.

EDITAR:

Isto é o que recebo do meu 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 optimizar a*a*a*a*a*aa (a*a*a)*(a*a*a)quando um é um número inteiro. Eu tentei com este comando:

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

Existem muitas bandeiras do gcc, mas nada sofisticado. Eles significam: Leia de stdin; use nível de otimização de O2; lista de linguagem assembly de saída em vez de um binário; a listagem deve usar a sintaxe da linguagem de montagem Intel; a entrada está no idioma C (geralmente o idioma é deduzido da extensão do arquivo de entrada, mas não há extensão de arquivo ao ler no stdin); e escreva para stdout.

Aqui está a parte importante da saída. Eu o anotei com alguns comentários indicando o que está acontecendo na linguagem assembly:

; 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

Estou usando o sistema GCC no Linux Mint 16 Petra, um derivado do Ubuntu. Aqui está a versão do gcc:

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

Como outros pôsteres observaram, essa opção não é possível no ponto flutuante, porque a aritmética do ponto flutuante não é associativa.

Como um número de ponto flutuante de 32 bits - como 1.024 - não é 1.024. Em um computador, 1.024 é um intervalo: de (1.024-e) a (1.024 + e), em que "e" representa um erro. Algumas pessoas não conseguem perceber isso e também acreditam que * em * a significa multiplicação de números de precisão arbitrária sem que haja erros associados a esses números. A razão pela qual algumas pessoas não percebem isso talvez seja o cálculo matemático que exercitavam nas escolas de ensino fundamental: trabalhando apenas com números ideais sem erros, e acreditando que não há problema em simplesmente ignorar "e" durante a multiplicação. Eles não veem o "e" implícito em "float a = 1.2", "a * a * a" e códigos C similares.

Caso a maioria dos programadores reconheça (e consiga executar) a idéia de que a expressão C a * a * a * a * a * a não esteja realmente trabalhando com números ideais, o compilador GCC estará livre para otimizar "a * a * a * a * a * a "digamos" t = (a * a); t * t * t "que requer um número menor de multiplicações. Infelizmente, porém, o compilador GCC não sabe se o programador que está escrevendo o código pensa que "a" é um número com ou sem erro. E assim o GCC fará apenas a aparência do código-fonte - porque é isso que o GCC vê com seu "olho nu".

... depois que você souber que tipo de programador você é, poderá usar a opção "-ffast-math" para informar ao GCC que "Ei, GCC, eu sei o que estou fazendo!". Isso permitirá que o GCC converta um * a * a * a * a * a em um pedaço de texto diferente - parece diferente de um * a * a * a * a * a - mas ainda calcula um número dentro do intervalo de erro de a * a * a * a * a * a. Tudo bem, já que você já sabe que está trabalhando com intervalos, e não com números ideais.

Nenhum pôster mencionou ainda a contração de expressões flutuantes (norma ISO C, 6.5p8 e 7.12.2). Se o FP_CONTRACTpragma estiver definido como ON, o compilador poderá considerar uma expressão comoa*a*a*a*a*a como uma única operação, como se avaliada exatamente com um único arredondamento. Por exemplo, um compilador pode substituí-lo por uma função de energia interna que é mais rápida e precisa. Isso é particularmente interessante, pois o comportamento é parcialmente controlado pelo programador diretamente no código-fonte, enquanto as opções do compilador fornecidas pelo usuário final às vezes podem ser usadas incorretamente.

O estado padrão do FP_CONTRACTpragma é definido pela implementação, para que um compilador possa fazer essas otimizações por padrão. Portanto, o código portátil que precisa seguir estritamente as regras da IEEE 754 deve explicitamente configurá-lo para OFF.

Se um compilador não suporta esse pragma, ele deve ser conservador, evitando qualquer otimização, caso o desenvolvedor tenha escolhido configurá-lo OFF.

O GCC não suporta esse pragma, mas com as opções padrão, ele assume que é ON; portanto, para destinos com uma FMA de hardware, se alguém quiser impedir a transformação a*b+cem fma (a, b, c), precisará fornecer uma opção como -ffp-contract=off(definir explicitamente o pragma como OFF) ou -std=c99(dizer ao GCC para se adequar a alguns Versão padrão C, aqui C99, siga o parágrafo acima). No passado, a última opção não estava impedindo a transformação, o que significa que o GCC não estava em conformidade com este ponto: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845

Como Lambdageek apontou, a multiplicação de flutuador não é associativa e você pode obter menos precisão, mas também quando obtém uma melhor precisão, pode argumentar contra a otimização, porque deseja uma aplicação determinística. Por exemplo, em cliente / servidor de simulação de jogo, em que todo cliente precisa simular o mesmo mundo em que você deseja que os cálculos de ponto flutuante sejam determinísticos.

Funções de biblioteca como "pow" são geralmente criadas com cuidado para gerar o mínimo erro possível (no caso genérico). Isso geralmente é alcançado aproximando funções com splines (de acordo com o comentário de Pascal, a implementação mais comum parece estar usando o algoritmo Remez )

fundamentalmente a seguinte operação:

pow(x,y);

tem um erro inerente de aproximadamente a mesma magnitude que o erro em qualquer multiplicação ou divisão única .

Enquanto a seguinte operação:

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

tem um erro inerente que é maior que 5 vezes o erro de uma única multiplicação ou divisão (porque você está combinando 5 multiplicações).

O compilador deve ter muito cuidado com o tipo de otimização que está fazendo:

1. se otimizar pow(a,6)a a*a*a*a*a*aque pode melhorar o desempenho, mas reduzir drasticamente a precisão de números de ponto flutuante.
2. se otimizar a*a*a*a*a*a para pow(a,6)ele pode realmente reduzir a precisão, pois "a" era algum valor especial que permite a multiplicação sem erro (uma potência de 2 ou algum número inteiro pequeno)
3. se otimizar pow(a,6)para (a*a*a)*(a*a*a)ou (a*a)*(a*a)*(a*a)ainda houver perda de precisão em comparação à powfunção.

Em geral, você sabe que para valores arbitrários de ponto flutuante "pow" tem melhor precisão do que qualquer função que você possa escrever, mas em alguns casos especiais várias multiplicações podem ter melhor precisão e desempenho, cabe ao desenvolvedor escolher o que é mais apropriado, eventualmente comentando o código para que ninguém mais "otimize" esse código.

A única coisa que faz sentido (opinião pessoal e, aparentemente, uma escolha no GCC que não inclua nenhum sinalizador de otimização ou compilador) a ser otimizado deve ser substituir "pow (a, 2)" por "a * a". Essa seria a única coisa sã que um fornecedor de compilador deve fazer.

Eu não esperava que esse caso fosse otimizado. Não pode ser muito frequente que uma expressão contenha subexpressões que possam ser reagrupadas para remover operações inteiras. Eu esperaria que os escritores de compiladores investissem seu tempo em áreas que provavelmente resultariam em melhorias visíveis, em vez de cobrir um caso extremo raramente encontrado.

Fiquei surpreso ao saber das outras respostas que essa expressão poderia realmente ser otimizada com as opções apropriadas do compilador. A otimização é trivial, ou é um caso extremo de uma otimização muito mais comum, ou os escritores do compilador foram extremamente minuciosos.

Não há nada errado em fornecer dicas para o compilador, como você fez aqui. É uma parte normal e esperada do processo de micro-otimização reorganizar declarações e expressões para ver quais diferenças elas trarão.

Embora o compilador possa ser justificado ao considerar as duas expressões para fornecer resultados inconsistentes (sem as opções apropriadas), não há necessidade de você ficar vinculado por essa restrição. A diferença será incrivelmente pequena - tanto que, se a diferença importa para você, você não deve usar aritmética padrão de ponto flutuante em primeiro lugar.

Já existem algumas boas respostas para essa pergunta, mas, para fins de completude, eu gostaria de salientar que a seção aplicável do padrão C é 5.1.2.2.3 / 15 (que é a mesma que a seção 1.9 / 9 no Padrão C ++ 11). Esta seção afirma que os operadores só podem ser reagrupados se forem realmente associativos ou comutativos.

Na verdade, o gcc pode fazer essa otimização, mesmo para números de ponto flutuante. Por exemplo,

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

torna-se

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

com -O -funsafe-math-optimizations. Essa reordenação viola a IEEE-754, no entanto, portanto, requer a bandeira.

Inteiros assinados, como Peter Cordes apontou em um comentário, podem fazer essa otimização sem, -funsafe-math-optimizationsuma vez que ela é mantida exatamente quando não há excesso e, se houver excesso, você obtém um comportamento indefinido. Então você recebe

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

com apenas -O. Para números inteiros não assinados, é ainda mais fácil, pois eles trabalham com potências mod de 2 e, portanto, podem ser reordenados livremente, mesmo diante do estouro.