Почему GCC не оптимизирует a*a*a*a*a*a до (a*a*a)*(a*a*a)?

Я делаю некоторую числовую оптимизацию для научного приложения. Одна вещь, которую я заметил, это то, что GCC оптимизирует вызов pow(a,2) составив его в a*a, но звонок pow(a,6) не оптимизирован и будет вызывать библиотечную функцию pow, что сильно тормозит производительность. (Напротив, Intel C++ Compiler, исполняемый файл 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 имеет похожее поведение.

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

12 ответов

Решение

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

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

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( ), который должен генерировать встроенное дерево умножения. Используйте это, если вы хотите поменять точность на производительность, но не хотите включать быструю математику.

Еще один похожий случай: большинство компиляторов не оптимизируют 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

Fortran (разработанный для научных вычислений) имеет встроенный оператор питания, и, насколько я знаю, компиляторы Fortran обычно оптимизируют повышение до целочисленных степеней аналогично тому, что вы описываете. C/C++, к сожалению, не имеет оператора power, только библиотечная функция 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 a с паренами), и позволяет вам иметь такую ​​оптимизацию без -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

GCC действительно оптимизирует a * a * a * a * a * a до (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, но ничего особенного. Они имеют в виду: читать со стандартного ввода; использовать уровень оптимизации O2; выводит список ассемблера вместо двоичного файла; в листинге должен использоваться синтаксис языка ассемблера Intel; ввод осуществляется на языке C (обычно язык определяется по расширению входного файла, но при чтении из stdin расширение файла отсутствует); и напиши в стандартный вывод.

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

    ; 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

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

Потому что 32-разрядное число с плавающей запятой, такое как 1.024, не равно 1.024. В компьютере 1,024 - это интервал: от (1,024-е) до (1,024+ е), где "е" обозначает ошибку. Некоторые люди не понимают этого и также считают, что * в * a означает умножение чисел произвольной точности без каких-либо ошибок, связанных с этими числами. Причиной, по которой некоторые люди не понимают этого, возможно, являются математические вычисления, которые они выполняли в начальных школах: работа только с идеальными числами без ошибок и уверенность в том, что можно просто игнорировать "е" при выполнении умножения. Они не видят "e", подразумеваемое в "float a=1.2", "a*a*a" и подобных кодах C.

Если большинство программистов признают (и смогут выполнять) идею о том, что выражение C a * a * a * a * a * a на самом деле не работает с идеальными числами, компилятор GCC тогда БЕСПЛАТНО оптимизирует "a * a * a * a * a * a "в скажем"t=(a*a); t*t*t", что требует меньшего числа умножений. Но, к сожалению, компилятор GCC не знает, думает ли программист, пишущий код, что "a" - это число с ошибкой или без нее. И поэтому GCC будет делать только то, на что похож исходный код - потому что это то, что GCC видит "невооруженным глазом".

... как только вы узнаете, какой вы программист, вы можете использовать переключатель "-ffast-math", чтобы сообщить GCC: "Эй, GCC, я знаю, что я делаю!" Это позволит GCC преобразовать a * a * a * a * a * a в другой фрагмент текста - он выглядит иначе, чем a * a * a * a * a * a - но все равно вычисляет число в интервале ошибок а * а * а * а * а * а. Это нормально, так как вы уже знаете, что работаете с интервалами, а не с идеальными числами.

Никто из авторов еще не упомянул о сокращении выражений с плавающей запятой (стандарт ISO C, 6.5p8 и 7.12.2). Если FP_CONTRACT Прагма установлена ​​в ONкомпилятору разрешено рассматривать выражение, такое как a*a*a*a*a*a как одна операция, как если бы она была оценена точно с одним округлением. Например, компилятор может заменить его внутренней функцией power, которая является более быстрой и точной. Это особенно интересно, так как поведение частично контролируется программистом непосредственно в исходном коде, в то время как параметры компилятора, предоставляемые конечным пользователем, могут иногда использоваться неправильно.

Состояние по умолчанию FP_CONTRACT Прагма определяется реализацией, так что компилятору разрешено выполнять такую ​​оптимизацию по умолчанию. Таким образом, переносимый код, который должен строго следовать правилам IEEE 754, должен явно установить его в OFF,

Если компилятор не поддерживает эту прагму, он должен быть консервативным, избегая любой такой оптимизации, в случае, если разработчик решил установить ее на OFF,

GCC не поддерживает эту прагму, но с параметрами по умолчанию она предполагает, что ON; таким образом, для целей с аппаратным FMA, если кто-то хочет предотвратить преобразование a*b+c Для FMA (A, B, C), необходимо предоставить такую ​​опцию, как -ffp-contract=off (явно установить прагму OFF) или же -std=c99 (чтобы сообщить GCC о соответствии некоторой стандартной версии C, здесь C99, таким образом следуйте приведенному выше параграфу). В прошлом последний вариант не препятствовал преобразованию, а это означает, что GCC не соответствовал этому пункту: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845

Библиотечные функции, такие как "pow", обычно тщательно создаются для получения минимально возможной ошибки (в общем случае). Обычно это достигается аппроксимацией функций сплайнами (согласно комментарию Паскаля, наиболее распространенной реализацией, по-видимому, является использование алгоритма Ремеза)

Принципиально следующая операция:

pow(x,y);

имеет присущую ошибку приблизительно такой же величины, как ошибка при любом одиночном умножении или делении.

Пока следующая операция:

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

имеет собственную ошибку, которая более чем в 5 раз превышает ошибку одиночного умножения или деления (потому что вы комбинируете 5 умножений).

Компилятор должен быть очень внимателен к той оптимизации, которую он выполняет:

  1. если оптимизировать pow(a,6) в a*a*a*a*a*a это может улучшить производительность, но радикально снизить точность для чисел с плавающей запятой.
  2. если оптимизировать a*a*a*a*a*a в pow(a,6) это может на самом деле снизить точность, потому что "а" было некоторым специальным значением, которое позволяет умножение без ошибок (степень 2 или небольшое целое число)
  3. если оптимизировать pow(a,6) в (a*a*a)*(a*a*a) или же (a*a)*(a*a)*(a*a) все еще может быть потеря точности по сравнению с pow функция.

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

Единственное, что имеет смысл (личное мнение и, очевидно, выбор в GCC без какой-либо конкретной оптимизации или флага компилятора) для оптимизации, - это заменить "pow(a,2)" на "a*a". Это было бы единственной разумной вещью, которую должен сделать поставщик компилятора.

Как указал Lambdageek, умножение чисел с плавающей запятой не ассоциативно, и вы можете получить меньшую точность, но также, когда вы получите лучшую точность, вы можете поспорить с оптимизацией, потому что вам нужно детерминированное приложение. Например, в клиент-сервер симуляции игры, где каждый клиент должен симулировать тот же мир, который вы хотите, чтобы вычисления с плавающей запятой были детерминированными.

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

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

Нет ничего плохого в предоставлении подсказок компилятору, как вы сделали здесь. Это нормальная и ожидаемая часть процесса микрооптимизации - перестановка операторов и выражений, чтобы увидеть, какие различия они принесут.

Хотя компилятор может быть оправдан при рассмотрении двух выражений для получения противоречивых результатов (без надлежащих переключателей), вам не нужно ограничиваться этим ограничением. Разница будет невероятно мала - настолько, что, если разница для вас важна, вам не следует использовать стандартную арифметику с плавающей запятой.

На этот вопрос уже есть несколько хороших ответов, но для полноты картины я хотел бы отметить, что применимым разделом стандарта C является 5.1.2.2.3/15 (который совпадает с разделом 1.9/9 в Стандарт C++11). В этом разделе говорится, что операторы могут быть перегруппированы, только если они действительно ассоциативны или коммутативны.

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

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

становится

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

с -O -funsafe-math-optimizations, Это переупорядочение нарушает IEEE-754, поэтому требует флаг.

Целые числа со знаком, как отметил Питер Кордес в комментарии, могут сделать эту оптимизацию без -funsafe-math-optimizations так как он выполняется именно тогда, когда переполнения нет, а при переполнении вы получаете неопределенное поведение. Итак, вы получаете

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

просто -O, Для целых чисел без знака это даже проще, поскольку они работают с модами степеней 2 и поэтому могут свободно переупорядочиваться даже при переполнении.

Другие вопросы по тегам