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

Question

Я занимаюсь численной оптимизацией научного приложения.Одна вещь, которую я заметил, заключается в том, что 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 ведет себя аналогично.

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

Answer 1

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

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

Answer 2

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

Answer 3

Другой аналогичный случай: большинство компиляторов не оптимизируют 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

Answer 4

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

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

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

EDIT:

Вот что я получаю от моего компилятора:

Для 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

Answer 5

GCC действительно оптимизирует 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  ; temp = x
imul eax, edi  ; temp = x * temp
imul eax, edi  ; temp = x * temp
imul eax, eax  ; temp = temp * temp

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

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

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

Answer 6

Поскольку 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 - но все равно вычисляет число в интервале ошибока * а * а * а * а * а.Это нормально, так как вы уже знаете, что работаете с интервалами, а не с идеальными числами.

Answer 7

Ни один из авторов еще не упомянул о сокращении выражений с плавающей запятой (стандарт 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

Answer 8

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

Answer 9

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

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

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

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

Answer 10

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

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

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) это может на самом деле снизить точность, потому что «a» было некоторым специальным значением, которое позволяет умножение без ошибок (степень 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". Это было бы единственной разумной вещью, которую должен делать поставщик компилятора.