Почему 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 имеет схожее поведение.

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

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

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

Lambdageek правильно указывает, что, поскольку ассоциативность не выполняется для чисел с плавающей запятой, «оптимизация» a*a*a*a*a*a to (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 пользователем (например, через -ffast-math ).

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

Другой подобный случай: большинство компиляторов не будут оптимизировать a + b + c + d (a + b) + (c + d) (это оптимизация, так как второе выражение может быть конвейеризовано лучше) и оценивать его как заданное (т.е. as (((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 ++, к сожалению, нет оператора мощности, только библиотечная функция pow() . Это не мешает умным компиляторам обрабатывать pow специально и вычислять его быстрее для особых случаев, но, похоже, они делают это реже …

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

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

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

Затем просто используйте его как power<6>(a) .

Это позволяет легко набирать полномочия (не нужно указывать 6 a с помощью 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 

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

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

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

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 -xc - 

Есть много флагов 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 

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

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

Состояние по умолчанию для 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

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

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

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

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

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

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

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

  • Каков наиболее эффективный способ для плавающего и двойного сравнения?
  • Преобразовать десятичное число в двойное?
  • проблема с плавающей запятой в R?
  • извлечение мантиссы и экспонента из двойного в c #
  • Сколько двойных чисел существует между 0.0 и 1.0?
  • Детали реализации оборудования с плавающей запятой
  • Ассемблер x86: сравнение с плавающей запятой
  • Сколько значимых цифр имеет float и double в java?
  • Точное хранение больших целых чисел
  • Float и двойной тип данных в Java
  • Каково первое целое число, которое плавающий IEEE 754 не может точно представлять?
  • Interesting Posts

    Невозможно прикрепить файл * .mdf как базу данных

    Как заставить Notepad ++ связывать тип файла с языком?

    Почему этот пример FINDSTR с несколькими буквальными поисковыми строками не находит совпадения?

    Как я могу настроить фокус (и отображать клавиатуру) на моем EditText программно

    Есть ли способ взломать пароль администратора Windows в Linux (используя файл SAM) БЕЗ его сброса?

    Есть ли способ изменить значение android: windowSoftInputMode из classа java?

    Есть ли ll команда symlink в Windows 7

    Угловая ng-repeat добавляет бутстрап-строку каждые 3 или 4 столбца

    Создать случайное число в LESS CSS?

    Как запустить отдельный проект без отладки в Visual Studio?

    Зачем мне сначала устанавливать значения Autologon в реестре до того, как он работает, и могу ли я исправить это?

    Как установить java jdk 7 на Snow Leopard

    Изменение целевых настроек процессора в Visual Studio 2010 Express

    java.security.InvalidAlgorithmParameterException: параметр trustAnchors должен быть не пустым в Linux, или почему пустая

    Легкий способ загрузки установочного компакт-диска WinXP с USB-накопителя USB?

    Давайте будем гением компьютера.