Вопрос: Почему GCC не оптимизирует a * a * a * a * a * a to (a * a * a) * (a * a * a)?


Я делаю некоторую численную оптимизацию в научном приложении. Одна вещь, которую я заметил, это то, что GCC оптимизирует вызов pow(a,2)путем его компиляции в a*a, но вызов pow(a,6)не оптимизирован и фактически вызовет библиотечную функцию pow, что значительно замедляет производительность. (В противоположность, Компилятор Intel C ++ , исполняемый файл icc, исключит вызов библиотеки для pow(a,6).)

Мне любопытно, что когда я заменяю pow(a,6)с a*a*a*a*a*aиспользуя GCC 4.5.1 и опции " -O3 -lm -funroll-loops -msse4», он использует 5 mulsdинструкции:

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

а если я напишу (a*a*a)*(a*a*a), он будет производить

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

что уменьшает количество команд умножения до 3. iccимеет подобное поведение.

Почему компиляторы не признают этот трюк оптимизации?


1958


источник


Ответы:


Потому как Плавающая математическая точка не ассоциативна , Способ группировки операндов при умножении с плавающей запятой влияет на числовую точность ответа.

В результате большинство компиляторов очень консервативны в отношении переупорядочения вычислений с плавающей запятой, если они не могут быть уверены, что ответ останется неизменным или если вы не скажете им, что вам не нужна численная точность. Например: -fassociative-mathвариант gcc, который позволяет gcc перезаписывать операции с плавающей запятой или даже -ffast-mathвариант, который позволяет еще более агрессивные компромиссы точности против скорости.


2555



Lambdageek правильно указывает, что, поскольку ассоциативность не выполняется для чисел с плавающей запятой, «оптимизация» a*a*a*a*a*aв (a*a*a)*(a*a*a)может изменить значение. Вот почему он запрещен C99 (если это специально не разрешено пользователем, через флаг компилятора или прагма). Как правило, предполагается, что программист написал то, что сделал по какой-то причине, и компилятор должен это уважать. Если ты хочешь (a*a*a)*(a*a*a), напишите это.

Это может быть болью писать; почему компилятор не может сделать [то, что вы считаете] правильным, когда вы используете pow(a,6)? Потому что это будет неправильно вещь которую нужно сделать. На платформе с хорошей математической библиотекой, pow(a,6)значительно более точен, чем a*a*a*a*a*aили (a*a*a)*(a*a*a), Чтобы предоставить некоторые данные, я провел небольшой эксперимент на своем Mac Pro, измеряя худшую ошибку при оценке ^ 6 для всех чисел с плавающей точкой с одной точностью между [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

С помощью powвместо дерева умножения уменьшает ошибку, связанную с коэффициент 4 , Компиляторы не должны (и, как правило, не делать), делать «оптимизации», которые увеличивают ошибку, если лицензия для этого не сделана пользователем (например, через -ffast-math).

Обратите внимание, что GCC предоставляет __builtin_powi(x,n)в качестве альтернативы pow( ), который должен генерировать встроенное дерево умножения. Используйте это, если вы хотите скомпрометировать точность для производительности, но не хотите включать ускоренную математику.


608



Другой подобный случай: большинство компиляторов не будут оптимизировать a + b + c + dв (a + b) + (c + d)(это оптимизация, поскольку второе выражение может быть конвейеризовано лучше) и оценивать его как заданное (т. е. как (((a + b) + c) + d)). Это тоже из-за угловых случаев:

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

Эти результаты 1.000000e-05 0.000000e+00


150



Fortran (предназначенный для научных вычислений) имеет встроенный оператор мощности, и насколько я знаю, компиляторы Fortran обычно оптимизируют повышение до целых полномочий аналогично тому, что вы описываете. К сожалению, у C / C ++ нет оператора мощности, только функция библиотеки pow(), Это не мешает умным компиляторам powспециально и вычисляя его быстрее для особых случаев, но кажется, что они делают это реже ...

Несколько лет назад я пытался сделать его более удобным для вычисления целочисленных мощностей оптимальным образом и придумал следующее. Это C ++, а не C, хотя и все еще зависит от того, как компилятор немного соображает, как оптимизировать / встроить вещи. Во всяком случае, надеюсь, что вы найдете это полезным на практике:

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

Уточнение для любопытных: это не находит оптимального способа вычисления полномочий, но поскольку нахождение оптимального решения является NP-полной задачей и в любом случае это стоит делать только для небольших держав (в отличие от использования pow), нет никаких оснований суетиться с деталями.

Тогда просто используйте его как power<6>(a),

Это позволяет легко набирать полномочия (не нужно указывать 6 as с parens), и позволяет вам проводить такую ​​оптимизацию без -ffast-mathв случае, если у вас есть что-то точно зависимое, например, компенсированное суммирование (пример, где порядок операций является существенным).

Возможно, вы также можете забыть, что это C ++ и просто использовать его в программе C (если он компилируется с помощью компилятора C ++).

Надеюсь, это может быть полезно.

РЕДАКТИРОВАТЬ:

Это то, что я получаю от своего компилятора:

Для a*a*a*a*a*a,

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

Для (a*a*a)*(a*a*a),

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

Для power<6>(a),

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

74



GCC фактически оптимизирует a * a * a * a * a * a to (a * a * a) * (a * a * a), когда a - целое число. Я пробовал эту команду:

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

Есть много флагов gcc, но ничего необычного. Они означают: Читайте от stdin; использовать уровень оптимизации O2; выводить список языков ассемблера вместо двоичного; в листинге должен использоваться синтаксис языка ассемблера Intel; вход на языке C (обычно язык выводится из расширения входного файла, но при чтении из stdin нет расширения файла); и напишите в stdout.

Вот важная часть вывода. Я комментировал это с некоторыми комментариями, указывающими, что происходит на ассемблере:

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

Я использую систему GCC на Linux Mint 16 Petra, производном Ubuntu. Вот версия gcc:

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

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


49



Because a 32-bit floating-point number - such as 1.024 - is not 1.024. In a computer, 1.024 is an interval: from (1.024-e) to (1.024+e), where "e" represents an error. Some people fail to realize this and also believe that * in a*a stands for multiplication of arbitrary-precision numbers without there being any errors attached to those numbers. The reason why some people fail to realize this is perhaps the math computations they exercised in elementary schools: working only with ideal numbers without errors attached, and believing that it is OK to simply ignore "e" while performing multiplication. They do not see the "e" implicit in "float a=1.2", "a*a*a" and similar C codes.

Should majority of programmers recognize (and be able to execute on) the idea that C expression a*a*a*a*a*a is not actually working with ideal numbers, the GCC compiler would then be FREE to optimize "a*a*a*a*a*a" into say "t=(a*a); t*t*t" which requires a smaller number of multiplications. But unfortunately, the GCC compiler does not know whether the programmer writing the code thinks that "a" is a number with or without an error. And so GCC will only do what the source code looks like - because that is what GCC sees with its "naked eye".

... once you know what kind of programmer you are, you can use the "-ffast-math" switch to tell GCC that "Hey, GCC, I know what I am doing!". This will allow GCC to convert a*a*a*a*a*a into a different piece of text - it looks different from a*a*a*a*a*a - but still computes a number within the error interval of a*a*a*a*a*a. This is OK, since you already know you are working with intervals, not ideal numbers.


48



I would not have expected this case to be optimized at all. It can't be very often where an expression contains subexpressions that can be regrouped to remove entire operations. I would expect compiler writers to invest their time in areas which would be more likely to result in noticeable improvements, rather than covering a rarely encountered edge case.

I was surprised to learn from the other answers that this expression could indeed be optimized with the proper compiler switches. Either the optimization is trivial, or it is an edge case of a much more common optimization, or the compiler writers were extremely thorough.

There's nothing wrong with providing hints to the compiler as you've done here. It's a normal and expected part of the micro-optimization process to rearrange statements and expressions to see what differences they will bring.

While the compiler may be justified in considering the two expressions to deliver inconsistent results (without the proper switches), there's no need for you to be bound by that restriction. The difference will be incredibly tiny - so much so that if the difference matters to you, you should not be using standard floating point arithmetic in the first place.


27



No posters have mentioned the contraction of floating expressions yet (ISO C standard, 6.5p8 and 7.12.2). If the FP_CONTRACT pragma is set to ON, the compiler is allowed to regard an expression such as a*a*a*a*a*a as a single operation, as if evaluated exactly with a single rounding. For instance, a compiler may replace it by an internal power function that is both faster and more accurate. This is particularly interesting as the behavior is partly controlled by the programmer directly in the source code, while compiler options provided by the end user may sometimes be used incorrectly.

The default state of the FP_CONTRACT pragma is implementation-defined, so that a compiler is allowed to do such optimizations by default. Thus portable code that needs to strictly follow the IEEE 754 rules should explicitly set it to OFF.

If a compiler doesn't support this pragma, it must be conservative by avoiding any such optimization, in case the developer has chosen to set it to OFF.

GCC doesn't support this pragma, but with the default options, it assumes it to be ON; thus for targets with a hardware FMA, if one wants to prevent the transformation a*b+c to fma(a,b,c), one needs to provide an option such as -ffp-contract=off (to explicitly set the pragma to OFF) or -std=c99 (to tell GCC to conform to some C standard version, here C99, thus follow the above paragraph). In the past, the latter option was not preventing the transformation, meaning that GCC was not conforming on this point: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845


27



As Lambdageek pointed out float multiplication is not associative and you can get less accuracy, but also when get better accuracy you can argue against optimisation, because you want a deterministic application. For example in game simulation client/server, where every client has to simulate the same world you want floating point calculations to be deterministic.


26



Library functions like "pow" are usually carefully crafted to yield the minimum possible error (in generic case). This is usually achieved approximating functions with splines (according to Pascal's comment the most common implementation seems to be using Remez algorithm)

fundamentally the following operation:

pow(x,y);

has a inherent error of approximately the same magnitude as the error in any single multiplication or division.

While the following operation:

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

has a inherent error that is greater more than 5 times the error of a single multiplication or division (because you are combining 5 multiplications).

The compiler should be really carefull to the kind of optimization it is doing:

  1. if optimizing pow(a,6) to a*a*a*a*a*a it may improve performance, but drastically reduce the accuracy for floating point numbers.
  2. if optimizing a*a*a*a*a*a to pow(a,6) it may actually reduce the accuracy because "a" was some special value that allows multiplication without error (a power of 2 or some small integer number)
  3. if optimizing pow(a,6) to (a*a*a)*(a*a*a) or (a*a)*(a*a)*(a*a) there still can be a loss of accuracy compared to pow function.

In general you know that for arbitrary floating point values "pow" has better accuracy than any function you could eventually write, but in some special cases multiple multiplications may have better accuracy and performance, it is up to the developer choosing what is more appropriate, eventually commenting the code so that noone else would "optimize" that code.

The only thing that make sense (personal opinion, and apparently a choice in GCC wichout any particular optimization or compiler flag) to optimize should be replacing "pow(a,2)" with "a*a". That would be the only sane thing a compiler vendor should do.


21