Почему 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 ++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-ffast-math

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

$ 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

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

Потому что 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.

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