Ускоряем pow

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

Сравнить точность алгоритмов можно прямо сейчас на этой странице.

В конце будет краткая памятка по тому, где и когда лучше применять какой из методов. При правильном выборе можно добиться увеличения скорости вычислений в 5 раз при погрешности ~1%, а иногда и вовсе без неё.

Содержание

1. Алгоритмы (5 штук)

2. Сравнение производительности

3. Сравнение точности

4. Вывод

На повестке дня у нас есть 5 алгоритмов: «Старая аппроксимация», «Бинарная степень», «Делящая быстрая степень», «Дробная” быстрая степень» и «Другая” аппроксимация».
Названия алгоритмам я придумал сам (за исключением бинарной степени), так как нигде не нашёл официальных версий, но вы можете называть их иначе.

Для расчета прироста скорости и погрешности будем сравнивать эти методы со стандартными функциями pow, Math.Pow и Math.pow в C++, C# и Java соответственно. О том, как производилось сравнение, будет сказано в частях “Сравнение производительности” и “Сравнение точности”.

Алгоритм: «Старая аппроксимация»

Увеличение скорости: в ~11 раз
Погрешность: <2%
Ограничения: приемлемая точность только для степеней от -1 до 1

Реализация в C++:

double OldApproximatePower(double b, double e) {     union {         double d;         long long i;     } u = { b };     u.i = (long long)(4606853616395542500L + e * (u.i - 4606853616395542500L));     return u.d; }

// C# double OldApproximatePower(double b, double e) {     long i = BitConverter.DoubleToInt64Bits(b);     i = (long)(4606853616395542500L + e * (i - 4606853616395542500L));     return BitConverter.Int64BitsToDouble(i); }

// Java double OldApproximatePower(double b, double e) {     long i = Double.doubleToLongBits(b);     i = (long)(4606853616395542500L + e * (i - 4606853616395542500L));     return Double.longBitsToDouble(i); }

Этот метод основан на алгоритме, использованном в игре Quake III Arena 2005 года. Он возводил число x в степень -0.5, т.е. находил значение:

float FastInvSqrt(float x) {     float xhalf = 0.5f * x;     int i = *(int*)&x;          // evil floating point bit level hacking     i = 0x5f3759df - (i >> 1);  // what the fuck?     x = *(float*)&i;      x = x*(1.5f-(xhalf*x*x));     return x; }

Узнал я об этом методе из статьи «Магическая константа» 0x5f3759df. В ней подробно объясняется как работает этот код и как его можно улучшить для работы с любой степенью и double’ми вместо float’ов. В моих кодах также есть магическая константа 4606853616395542500L. Нашёл я её по следующей формуле (она описана в статье выше):

//C# or Java long doubleApproximator = (long)((1L << 52) * ((1L << 10) - 1.0730088));

Число 1.0730088 было подобрано вручную для достижения наибольшей точности вычислений.

Алгоритм: Бинарное возведение в степень

Увеличение скорости: в среднем в ~7.5 раз, преимущество сохраняется до возведения чисел в степень 134217728 в C++/C# и 4096 в Java.
Погрешность: нет, но стоит отметить, что операция умножения не ассоциативна для чисел с плавающей точкой, т.е. 1.21 * 1.21 не то же самое, что 1.1 * 1.1 * 1.1 * 1.1, однако при сравнении со стандартными функциями погрешности, как уже сказано ранее, не возникает.
Ограничения: степень должна быть целым числом не меньше 0

Реализация в C++:

double BinaryPower(double b, unsigned long long e) {        double v = 1.0;        while(e != 0) {               if((e & 1) != 0) {                      v *= b;               }               b *= b;               e >>= 1;        }        return v; }

// C# double BinaryPower(double b, UInt64 e) {     double v = 1d;     while(e != 0) {         if((e & 1) != 0) {             v *= b;         }         b *= b;         e >>= 1;     }     return v; }

// Java double BinaryPower(double b, long e) {        double v = 1d;        while(e > 0) {               if((e & 1) != 0) {                      v *= b;               }               b *= b;               e >>= 1;        }        return v; }

Широко известный алгоритм для возведения любого числа в целую степень с абсолютной точностью. Принцип действия прост: есть целая степень e, чтобы получить число b в этой степени нужно возвести это число во все степени 1, 2, 4, … 2ⁿ (в коде этому соответствует b *= b), каждый раз сдвигая биты e вправо (e >>= 1) пока оно не равно 0 и тогда, когда последний бит e не равен нулю ((e & 1) != 0), домножать результат v на полученное b.
Пример: возвести 2 в степень 5.
v = 1, e = 5 = 101₂, b = 2
Шаги цикла:

1. e = 101₂ — последний 1 → v *= b → v = 2
   b *= b → b = 4
   e >>= 1 → e = 10₂ = 2

2. e = 10₂ — последний 0 → пропускаем
   b *= b → b = 16
   e >>= 1 → e = 1

3. e = 1₂ — последний 1 → v *= b → v = 32
   …
   e = 0 → выход из цикла

Результат: v = 32, что и есть 2⁵.

Алгоритм: «Делящая быстрая степень»

Увеличение скорости: в ~3.5 раз
Погрешность: ~13%
Примечание: в коде ниже присутствуют проверки для особых входных данных. Без них код работает всего на 10% быстрее, но погрешность возрастает в десятки раз (особенно при использовании отрицательных степеней).

Реализация в C++:

double FastPowerDividing(double b, double e) {        if(b == 1.0 || e == 0.0) {               return 1.0;        }          double eAbs = fabs(e);        double el = ceil(eAbs);        double basePart = OldApproximatePower(b, eAbs / el);        double result = BinaryPower(basePart, (unsigned long long)el);          if(e < 0.0) {               return 1.0 / result;        }        return result; }

// C# double FastPowerDividing(double b, double e) {        if(b == 1d || e == 0d) {               return 1d;        }        var eAbs = Math.Abs(e);        var el = Math.Ceiling(eAbs);        var basePart = OldApproximatePower(b, eAbs / el);        var result = BinaryPower(basePart, (UInt64)el);          if(e < 0d) {               return 1d / result;        }        return result; }

// Java double FastPowerDividing(double b, double e) {        if(b == 1d || e == 0d) {               return 1d;        }        var eAbs = Math.abs(e);        var el = Math.ceil(eAbs);        var basePart = OldApproximatePower(b, eAbs / el);        var result = BinaryPower(basePart, (long)el);          if(e < 0d) {               return 1d / result;        }        return result; }

Узнав о методе аппроксимации чисел в степенях от -1 до 1 и о бинарном методе, мне захотелось объединить их для создания функции, которая могла бы быстро возводить число в любую степень. Для этого я придумал следующую формулу:

Мы разбиваем степень на две части: e / el, которая всегда меньше или равна 1, и el, которая является целым числом. Теперь для расчета x^(e / el) мы можем использовать “старую” аппроксимацию, а для x^el — бинарную степень.Таким образом, объединяя этих два узкоспециализированных метода, мы получили универсальный метод. Но эту идею можно реализовать по-другому.

Алгоритм: «Дробная быстрая степень»

Увеличение скорости: в ~4.4 раза
Погрешность: ~0.7%

Реализация в C++:

double FastPowerFractional(double b, double e) {        if(b == 1.0 || e == 0.0) {               return 1.0;        }          double absExp = fabs(e);        unsigned long long eIntPart = (long long)absExp;        double eFractPart = absExp - eIntPart;        double result = OldApproximatePower(b, eFractPart) * BinaryPower(b, eIntPart);        if(e < 0.0) {               return 1.0 / result;        }        return result; }

// C# double FastPowerFractional(double b, double e) {        if(b == 1d || e == 0d) {               return 1d;        }          double absExp = Math.Abs(e);        UInt64 eIntPart = (UInt64)absExp;        double eFractPart = absExp - eIntPart;        double result = OldApproximatePower(b, eFractPart) * BinaryPower(b, eIntPart);        if(e < 0d) {               return 1d / result;        }        return result; }

// Java double FastPowerFractional(double b, double e) {        if(b == 1d || e == 0d) {               return 1d;        }          double absExp = Math.abs(e);        long eIntPart = (long)absExp;        double eFractPart = absExp - eIntPart;        double result = OldApproximatePower(b, eFractPart) * BinaryPower(b, eIntPart);        if(e < 0d) {               return 1d / result;        }        return result; }

По сути, любое число состоит из суммы двух частей: целой и дробной. Целую можно использовать для возведения основания в степень при помощи бинарного возведения, а дробную — при помощи “старой” аппроксимации.

В результате получаем следующую формулу:

Она, в отличии от формулы “делящего” метода, никак не искажает дробную часть. Это позволяет добиться намного большей точности.

Алгоритм: «Другая аппроксимация»

Увеличение скорости: в ~9 раз
Погрешность: <1.5%
Ограничения: точность стремительно падает при повышении абсолютного значения степени и остается приемлемой в промежутке [-10, 10]

Реализация в C++:

double AnotherApproximatePower(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; }

double AnotherApproxPower(double a, double b) {        int tmp = (int)(BitConverter.DoubleToInt64Bits(a) >> 32);        int tmp2 = (int)(b * (tmp - 1072632447) + 1072632447);        return BitConverter.Int64BitsToDouble(((long)tmp2) << 32); }

double AnotherApproxPower(double a, double b) {        int tmp = (int)(Double.doubleToLongBits(a) >> 32);        int tmp2 = (int)(b * (tmp - 1072632447) + 1072632447);        return Double.longBitsToDouble(((long)tmp2) << 32); }

Про историю этого алгоритма я ничего не знаю, я просто нашёл его тут: Optimized pow() approximation for Java, C / C++, and C#. Возможно, если использовать его в “делящей–” и “дробной быстрых степенях» вместо “старой” аппроксимации, можно достигнуть лучшей точности ценой немного меньшей скорости.

Сравнение производительности

Сравнение производительности производилось следующим образом: генерируем 500000 чисел-оснований в промежутке от 0.0 до 99999.0 и 500000 чисел-степеней в промежутке от A до B. Запоминаем текущее время, запускаем цикл на 500000 итераций, вычисляем значение основания в степени через функцию f и результат суммируем в calculationResult. По окончанию цикла снова замеряем время, разница во времени и есть время выполнения. Данная процедура повторяется 20 раз, конечный результат — усредненный за все 20 тестов.

Псевдокод сравнения производительности в C++:

(long long iterationsCount = 500000, double* bases, double* exps)   double calculationResult = 0.0;   double* base = bases;   double* exp = exps;   double* baseEnd = base + iterationsCount;   auto start = std::chrono::high_resolution_clock::now();   while(base < baseEnd) {          calculationResult += f(*base++, *exp++);   }   auto finish = std::chrono::high_resolution_clock::now();   auto time = std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count();

Аналогично измерялась скорость в C# и Java. В репозитории проекта можно посмотреть на реальный код для сравнения производительности в C++, C# и Java.

Тесты производились в каждом языке для степеней в промежутках [-10.5, 0], [0, 2], [0, 10.5], [0, 25.75], [0, 55.5]. Прирост скорости каждого метода по сравнению со стандартным в каждом языке для каждого сета степеней изображен на графиках ниже:

Рассмотреть подробнее результаты тестов можно посмотреть в этой таблице.

Тесты проводились на i5-10300H, 19.8 DDR4 GB usable RAM, 64-битная платформа.
C++: MSVC + /O2 + /Oi + /Ot
C#: optimize code

Сравнение точности

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

Ниже на этой же странице можно самим провести сравнение точности каждого метода для всех степеней, лежащих в заданном промежутке:

Вывод

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

Код всего проекта доступен на гитхабе. Там вы найдете и реализации всех алгоритмов в трех языках, и программы для сравнения точности для каждого языка, и код сайта.

Спасибо за внимание. Надеюсь эта статья поможет сделать ваш код немного быстрее.