Почему pow (int, int) так медленно?
Я работаю над несколькими упражнениями Эйлера для улучшения моих знаний о С++.
Я написал следующую функцию:
int a = 0,b = 0,c = 0;
for (a = 1; a <= SUMTOTAL; a++)
{
for (b = a+1; b <= SUMTOTAL-a; b++)
{
c = SUMTOTAL-(a+b);
if (c == sqrt(pow(a,2)+pow(b,2)) && b < c)
{
std::cout << "a: " << a << " b: " << b << " c: "<< c << std::endl;
std::cout << a * b * c << std::endl;
}
}
}
Это вычисляется за 17 миллисекунд.
Однако, если я изменяю строку
if (c == sqrt(pow(a,2)+pow(b,2)) && b < c)
to
if (c == sqrt((a*a)+(b*b)) && b < c)
вычисление происходит через 2 миллисекунды. Есть ли некоторые очевидные детали реализации pow(int, int)
, которые я пропускаю, что делает первое выражение вычисляемым так медленнее?
Ответы
Ответ 1
pow()
работает с реальными числами с плавающей запятой и использует под капотом формулу
pow(x,y) = e^(y log(x))
для вычисления x^y
. int
преобразуются в double
перед вызовом pow
. (log
- натуральный логарифм, основанный на e)
x^2
, используя pow()
, поэтому медленнее, чем x*x
.
Изменить на основе соответствующих комментариев
- Использование
pow
даже с целыми показателями может привести к неправильным результатам (PaulMcKenzie)
- В дополнение к использованию математической функции с двойным типом,
pow
является вызовом функции (в то время как x*x
не является) (jtbandes)
- Многие современные компиляторы фактически оптимизируют pow с постоянными целыми аргументами, но на это не следует полагаться.
Ответ 2
Вы выбрали один из самых медленных способов проверки
c*c == a*a + b*b // assuming c is non-negative
Это компилируется с тремя целыми умножениями (один из которых может быть выведен из цикла). Даже без pow()
вы все еще конвертируете в double
и получаете квадратный корень, что ужасно для пропускной способности. (А также латентность, но предсказание ветвей + спекулятивное выполнение на современных CPU означает, что латентность здесь не является фактором).
Инструкция Intel Haswell SQRTSD имеет пропускную способность одного на 8-14 циклов (источник: таблицы инструкций Agner Fog), поэтому даже если ваш sqrt()
поддерживает насыщенность исполняемого блока FP sqrt, он все еще примерно в 4 раза медленнее, чем то, что я получил gcc для испускания (ниже).
Вы также можете оптимизировать условие цикла для выхода из цикла, когда часть условия b < c
становится ложной, поэтому компилятор должен выполнить только одну версию этой проверки.
void foo_optimized()
{
for (int a = 1; a <= SUMTOTAL; a++) {
for (int b = a+1; b < SUMTOTAL-a-b; b++) {
// int c = SUMTOTAL-(a+b); // gcc won't always transform signed-integer math, so this prevents hoisting (SUMTOTAL-a) :(
int c = (SUMTOTAL-a) - b;
// if (b >= c) break; // just changed the loop condition instead
// the compiler can hoist a*a out of the loop for us
if (/* b < c && */ c*c == a*a + b*b) {
// Just print a newline. std::endl also flushes, which bloats the asm
std::cout << "a: " << a << " b: " << b << " c: "<< c << '\n';
std::cout << a * b * c << '\n';
}
}
}
}
Скомпилирует (с gcc6.2 -O3 -mtune=haswell
) код с этим внутренним циклом. Полный код проводник компилятора Godbolt.
# a*a is hoisted out of the loop. It in r15d
.L6:
add ebp, 1 # b++
sub ebx, 1 # c--
add r12d, r14d # ivtmp.36, ivtmp.43 # not sure what this is or why it in the loop, would have to look again at the asm outside
cmp ebp, ebx # b, _39
jg .L13 ## This is the loop-exit branch, not-taken until the end
## .L13 is the rest of the outer loop.
## It sets up for the next entry to this inner loop.
.L8:
mov eax, ebp # multiply a copy of the counters
mov edx, ebx
imul eax, ebp # b*b
imul edx, ebx # c*c
add eax, r15d # a*a + b*b
cmp edx, eax # tmp137, tmp139
jne .L6
## Fall-through into the cout print code when we find a match
## extremely rare, so should predict near-perfectly
В Intel Haswell все эти инструкции по 1 мкп каждый. (И макро-предохранитель cmp/jcc попадает в сопоставления и ветвления.) Таким образом, 10 fused-domain uops, которые могут выдаваться на одной итерации за 2,5 цикла.
Haswell запускает imul r32, r32
с пропускной способностью одной итерации за такт, поэтому два умножения внутри внутреннего цикла не насыщают порт 1 с двумя умножениями на 2.5c. Это оставляет место, чтобы впитать неизбежные конфликты ресурсов из ADD и SUB, крадущего порт 1.
Мы даже не близки к каким-либо другим узким местам исполнения, поэтому узкое место переднего плана является единственной проблемой, и это должно выполняться на одной итерации за 2,5 цикла на Intel Haswell, а затем.
Loop-unrolling может помочь здесь уменьшить количество попыток на проверку. например используйте lea ecx, [rbx+1]
для вычисления b + 1 для следующей итерации, поэтому мы можем imul ebx, ebx
не использовать MOV, чтобы сделать его неразрушающим.
Возможно также снижение прочности. Учитывая b*b
, мы могли бы попытаться вычислить (b-1) * (b-1)
без IMUL. (b-1) * (b-1) = b*b - 2*b + 1
, поэтому, возможно, мы можем сделать lea ecx, [rbx*2 - 1]
, а затем вычесть из b*b
. (Нет меток адресации, которые вычитают вместо добавления. Хм, возможно, мы могли бы сохранить -b
в регистре и подсчитать до нуля, поэтому мы могли бы использовать lea ecx, [rcx + rbx*2 - 1]
для обновления b*b
в ECX, учитывая -b
в EBX).
Если вы на самом деле не столкнулись с пропускной способностью IMUL, это может привести к большему количеству ошибок, а не победе. Было бы забавно видеть, насколько хорошо компилятор будет делать это с уменьшением силы в источнике С++.
Возможно, вы могли бы также векторизовать это с помощью SSE или AVX, одновременно проверяя 4 или 8 последовательных значений b
. Поскольку хиты действительно редки, вы просто проверяете, был ли у кого-то из 8 ударов, а затем разобрал, какой из них был в редком случае, когда был матч.
См. также x86 tag wiki для большей информации по оптимизации.