Почему pow (int, int) так медленно?

Question

Я работаю над несколькими упражнениями Эйлера для улучшения моих знаний о C++.

Я написал следующую функцию:

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) 

в

if (c == sqrt((a*a)+(b*b)) && b < c) 

вычисление происходит в 2 миллисекунды. Есть ли какая-то очевидная деталь реализации pow(int, int), которую я пропускаю, что делает первое выражение вычисляемым так медленнее?

Answer 1

pow() работает с числами реальной плавающей точкой и использует под капотом формулу

pow(x,y) = e^(y log(x)) 

для расчета x^y. int преобразуются в double перед вызовом pow. (log это натуральный логарифм, электронной основе) Поэтому

x^2 использованием pow() происходит медленнее, чем x*x.

Редактировать на основе соответствующих комментариев

• Использования pow даже с целыми показателями могут давать неправильные результаты (PaulMcKenzie)
• В дополнении к использованию математической функции с двойной типа, pow является функцией звонок (в то время как x*x нет) (jtbandes)
• Многие современные компиляторы фактически оптимизируют силу wi го постоянных целых аргументов, но на это не следует полагаться.

Answer 2

Вы выбрали один из самых медленных возможных способов проверить

c*c == a*a + b*b // assuming c is non-negative 

Это компилирует три целочисленных умножений (один из которых может быть поднят из петли). Даже без pow() вы все еще конвертируете в double и получаете квадратный корень, что ужасно для пропускной способности. (А также латентность, но предсказание ветвей + спекулятивное выполнение на современных CPU означает, что латентность здесь не является фактором).

инструкция SQRTSD Intel Haswell имеет пропускную способность одного за 8-14 циклов (source: Agner Fog's instruction tables), так что даже если ваша версия 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), чтобы кодировать с этим внутренним контуром. См. Полный код на the Godbolt compiler explorer.

# a*a is hoisted out of the loop. It's 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's 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, which can issue at one iteration per 2.5 cycles.

Haswell запускает imul r32, r32 с пропускной способностью одной итерации за такт, поэтому два умножения внутри внутреннего контура не насыщают порт 1 при двух умножениях на 2,5 c. Это оставляет место для того, чтобы впитать неизбежные конфликты ресурсов с ADD и SUB, крадующим порт 1.

Мы даже не близки к каким-либо другим узким местам исполнения, поэтому является узким местом переднего плана, и это должен работать на одной итерации за 2,5 цикла на Intel Haswell и позже.

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

---

Вы могли бы также векторизации это с SSE или AVX, проверка 4 или 8 последовательных b значения параллельно. Поскольку хиты действительно редки, вы просто проверяете, был ли у кого-то из 8 хит, а затем разобраться, какой из них был в редком случае, когда был матч.

См. Также ссылку на wiki для x86.