В этой статье я хочу поделиться несколькими нестандартными алгоритмами для быстрого возведения числа в степень, а также продемонстрировать их реализацию и сравнить их быстродействие в C++, C# и Java.
Сравнить точность алгоритмов можно прямо сейчас на этой странице.
В конце будет краткая памятка по тому, где и когда лучше применять какой из методов. При правильном выборе можно добиться увеличения скорости вычислений в 5 раз при погрешности ~1%, а иногда и вовсе без неё.
Названия алгоритмам я придумал сам (за исключением бинарной степени), так как нигде не нашёл официальных версий, но вы можете называть их иначе.
Для расчета прироста скорости и погрешности будем сравнивать эти методы со стандартными функциями pow, Math.Pow и Math.pow в C++, C# и Java соответственно. О том, как производилось сравнение, будет сказано в частях “Сравнение производительности” и “Сравнение точности”.
Погрешность: <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# и Java
// 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 было подобрано вручную для достижения наибольшей точности вычислений.
Погрешность: нет, но стоит отметить, что операция умножения не ассоциативна для чисел с плавающей точкой, т.е. 1.21 * 1.21 не то же самое, что 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# и Java
// 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, … 2n (в коде этому соответствует b *= b), каждый раз сдвигая биты e вправо (e >>= 1) пока оно не равно 0 и тогда, когда последний бит e не равен нулю ((e & 1) != 0), домножать результат v на полученное b.
Пример: возвести 2 в степень 5.
v = 1, e = 5 = 1012, b = 2
Шаги цикла:
Погрешность: ~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, (long long)el);
if(e < 0.0) {
return 1.0 / result;
}
return result;
}
Реализация в C# и Java
// 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, (long)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 - бинарную степень.Таким образом, объединяя этих два узкоспециализированных метода, мы получили универсальный метод. Но эту идею можно реализовать по-другому.
Погрешность: ~0.7%
Реализация в C++:
double FastPowerFractional(double b, double e) {
if(b == 1.0 || e == 0.0) {
return 1.0;
}
double absExp = fabs(e);
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# и Java
// C#
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;
}
// 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;
}
По сути, любое число состоит из суммы двух частей: целой и дробной. Целую можно использовать для возведения основания в степень при помощи бинарного возведения, а дробную - при помощи “старой” аппроксимации.
В результате получаем следующую формулу:
Она, в отличии от формулы “делящего” метода, никак не искажает дробную часть. Это позволяет добиться намного большей точности.
Погрешность: <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;
}
Реализация в C# и Java
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#. Возможно, если использовать его в “делящей–” и “дробной быстрых степенях" вместо “старой” аппроксимации, можно достигнуть лучшей точности ценой немного меньшей скорости.
Псевдокод сравнения производительности в 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 страницу, на которой можно построить графики реальных и полученных одним из методов значений чисел, лежащих в заданном промежутке и возведенных в указанную степень:
Ниже на этой же странице можно самим провести сравнение точности каждого метода для всех степеней, лежащих в заданном промежутке:
Код всего проекта доступен на гитхабе. Там вы найдете и реализации всех алгоритмов в трех языках, и программы для сравнения точности для каждого языка, и код сайта.
Спасибо за внимание. Надеюсь эта статья поможет сделать ваш код немного быстрее.
habr.com
Сравнить точность алгоритмов можно прямо сейчас на этой странице.
В конце будет краткая памятка по тому, где и когда лучше применять какой из методов. При правильном выборе можно добиться увеличения скорости вычислений в 5 раз при погрешности ~1%, а иногда и вовсе без неё.
Содержание
- Алгоритмы (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# и Java
// 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, однако при сравнении со стандартными функциями погрешности, как уже сказано ранее, не возникает.
Ограничения: степень должна быть целым числом не меньше 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# и Java
// 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, … 2n (в коде этому соответствует b *= b), каждый раз сдвигая биты e вправо (e >>= 1) пока оно не равно 0 и тогда, когда последний бит e не равен нулю ((e & 1) != 0), домножать результат v на полученное b.
Пример: возвести 2 в степень 5.
v = 1, e = 5 = 1012, b = 2
Шаги цикла:
- e = 1012 - последний 1 → v *= b → v = 2
b *= b → b = 4
e >>= 1 → e = 102 = 2 - e = 102 - последний 0 → пропускаем
b *= b → b = 16
e >>= 1 → e = 1 - e = 12 - последний 0 → v *= b → v = 32
...
e = 0 → выход из цикла
Алгоритм: "Делящая быстрая степень"
Увеличение скорости: в ~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, (long long)el);
if(e < 0.0) {
return 1.0 / result;
}
return result;
}
Реализация в C# и Java
// 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, (long)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);
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# и Java
// C#
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;
}
// 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;
}
Реализация в C# и Java
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 страницу, на которой можно построить графики реальных и полученных одним из методов значений чисел, лежащих в заданном промежутке и возведенных в указанную степень:
Ниже на этой же странице можно самим провести сравнение точности каждого метода для всех степеней, лежащих в заданном промежутке:
Вывод
Каждый из перечисленных мною методов дает различную точность и скорость. Поэтому, прежде чем использовать стандартную функцию возведения в степень, стоит проанализировать какие у вас будут входные данные, и насколько вам важна точность. Выбор правильного метода может существенно ускорить работу вашей программы. Чтобы сделать этот процесс проще, я подготовил вот такую шпаргалку:
Код всего проекта доступен на гитхабе. Там вы найдете и реализации всех алгоритмов в трех языках, и программы для сравнения точности для каждого языка, и код сайта.
Спасибо за внимание. Надеюсь эта статья поможет сделать ваш код немного быстрее.
Ускоряем pow
В этой статье я хочу поделиться несколькими нестандартными алгоритмами для быстрого возведения числа в степень, а также продемонстрировать их реализацию и сравнить их быстродействие в C++, C# и Java....
habr.com