Ускоряем pow

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

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

В конце будет краткая памятка по тому, где и когда лучше применять какой из методов. При правильном выборе можно добиться увеличения скорости вычислений в 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 ^[2]. В ней подробно объясняется как работает этот код и как его можно улучшить для работы с любой степенью и double’ми вместо float’ов. В моих кодах также есть магическая константа 4606853616395542500L. Нашёл я её по следующей формуле (она описана в статье выше ^[2]):

//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#
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₂ - последний 0 → 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, (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, (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#
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;
}

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# ^[3]. Возможно, если использовать его в “делящей–” и “дробной быстрых степенях" вместо “старой” аппроксимации, можно достигнуть лучшей точности ценой немного меньшей скорости.

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

Сравнение производительности производилось следующим образом: генерируем 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++ ^[4], C# ^[5] и Java ^[6].

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

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

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

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

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

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

Вывод

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

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

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