Что быстрее, чем std:: pow?

Ответ 1

Похоже, что у Мартина Анкерла есть несколько статей по этому поводу, Оптимизированный Approximative pow() в C/С++ является одним, и он имеет две быстрых версии, одна из них выглядит следующим образом:

inline double fastPow(double a, double b) {
  union {
    double d;
    int x[2];
  } u = { a };
  u.x[1] = (int)(b * (u.x[1] - 1072632447) + 1072632447);
  u.x[0] = 0;
  return u.d;
}

который опирается на тип punning через объединение, которое является undefined поведением в С++, из проекта стандартного раздела 9.5 [class.union]:

В объединении не более чем один из нестатических членов данных может быть активным в любое время, то есть значение at большинство из нестатических элементов данных могут быть сохранены в союзе в любое время. [...]

но большинство компиляторов, включая gcc поддерживают это с четко определенным поведением:

Практика чтения из другого члена профсоюза, чем тот, который был недавно написан (называемый "пингом типа" ), является обычным явлением. Даже при использовании -fstrict-aliasing допускается использование пула типа, при условии, что доступ к памяти осуществляется через тип объединения

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

Он также ссылается на второй Оптимизированное приближение pow() для Java, C/С++ и С#.

Первая статья также ссылается на его микрообъекты здесь

Ответ 2

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

Ответ 3

Насколько велики ваши целые числа? Известны ли они во время компиляции? Гораздо лучше вычислить x^2 как x*x, а не pow(x,2). Примечание. Почти все приложения pow() к целочисленной мощности включают в себя повышение числа до второй или третьей мощности (или мультипликативного обратного в случае отрицательных показателей). Использование pow() является избыточным в таких случаях. Используйте шаблон для этих малых целых степеней или просто используйте x*x.

Если целые числа малы, но неизвестны во время компиляции, скажем между -12 и +12, умножение по-прежнему будет бить pow() и не потеряет точность. Вам не нужно одиннадцать умножений для вычисления x ^ 12. Четверо будут делать. Используйте тот факт, что x ^ (2n) = (x ^ n) ^ 2 и x ^ (2n + 1) = x * ((x ^ n) ^ 2). Например, x ^ 12 есть ((x * x * x) ^ 2) ^ 2. Два умножения для вычисления x ^ 3 (x * x * x), еще один для вычисления x ^ 6 и один последний для вычисления x ^ 12.

Ответ 4

ДА! Очень быстро, если вам нужно только 'y'/'n' как long/int, что позволяет вам избежать медленной функции FPU FSCALE. Это версия Agner Fog x86, оптимизированная для рук, если вам нужны результаты только с 'y'/'n' в качестве INT. Я обновил его до __fastcall/__declspec (голый) для скорости/размера, использовал ECX для передачи 'n' (числа с плавающей запятой всегда передаются в стек для 32-битных MSVC++), так что с моей стороны очень незначительные изменения, в основном это Агнер работа. Он был протестирован/отлажен/скомпилирован на MS Visual VC++ 2005 Express/Pro, поэтому в новых версиях все должно быть в порядке. Точность по отношению к универсальной функции CRT pow() очень хорошая.

extern double __fastcall fs_power(double x, long n);

// Raise 'x' to the power 'n' (INT-only) in ASM by the great Agner Fog!

__declspec(naked) double __fastcall fs_power(double x, long n) { __asm {
    MOV  EAX, ECX     ;// Move 'n' to eax
;// abs(n) is calculated by inverting all bits and adding 1 if n < 0:
    CDQ               ;// Get sign bit into all bits of edx
    XOR  EAX, EDX     ;// Invert bits if negative
    SUB  EAX, EDX     ;// Add 1 if negative. Now eax = abs(n)
    JZ   RETZERO      ;// End if n = 0
    FLD1              ;// ST(0) = 1.0 (FPU push1)
    FLD  QWORD PTR [ESP+4] ;// Load 'x' : ST(0) = 'x', ST(1) = 1.0 (FPU push2)
    JMP  L2           ;// Jump into loop
L1: ;// Top of loop
    FMUL ST(0), ST(0) ;// Square x
L2: ;// Loop entered here
    SHR  EAX, 1       ;// Get each bit of n into carry flag
    JNC  L1           ;// No carry. Skip multiplication, goto next
    FMUL ST(1), ST(0) ;// Multiply by x squared i times for bit # i
    JNZ  L1           ;// End of loop. Stop when nn = 0
    FSTP ST(0)        ;// Discard ST(0) (FPU Pop1)
    TEST EDX, EDX     ;// Test if 'n' was negative
    JNS  RETPOS       ;// Finish if 'n' was positive
    FLD1              ;// ST(0) = 1.0, ST(1) = x^abs(n)
    FDIVR             ;// Reciprocal
RETPOS:    ;// Finish, success!
    RET  4 ;//(FPU Pop2 occurs by compiler on assignment

RETZERO:
    FLDZ   ;// Ret 0.0, fail, if n was 0
    RET  4
}}