Ускоряем pow21.10.2021 05:47

В этой статье я хочу поделиться несколькими нестандартными алгоритмами для быстрого возведения числа в степень, а также продемонстрировать их реализацию и сравнить их быстродействие в 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# и 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, т.е. находил значение:

\frac{1}{\sqrt{x}}Разработчики для этого написали такую функцию

float FastInvSqrt(float x) {
    float xhalf = 0.5f * x;
    int i = *(int*)&x;
    i = 0x5f3759df - (i >> 1);
    x = *(float*)&i;
    x = x*(1.5f-(xhalf*x*x));
    return x;
}

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

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

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

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

Увеличение скорости: в среднем в ~7.5 раз, преимущество сохраняется до возведения чисел в степень 134217728 в C++/C# и 4096 в Java.
Погрешность: нет
Ограничения: степень должна быть целым числом не меньше 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
Шаги цикла:

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

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

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

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

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

Увеличение скорости: в ~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 и о бинарном методе, мне захотелось объединить их для создания функции, которая могла бы быстро возводить число в любую степень. Для этого я придумал следующую формулу:

el = \left | \left \lceil e \right \rceil \right |\\ x^e = (x^{\frac{e}{el}})^{el}

Мы разбиваем степень на две части: 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;
}

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

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

el = \left \lfloor e \right \rfloor\\ x^e = x^{el}*x^{e - el}

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

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

Увеличение скорости: в ~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(finish - start).count();

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

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

image-loader.svgimage-loader.svgimage-loader.svg

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

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

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

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

image-loader.svgimage-loader.svg

Вывод

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

image-loader.svg

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

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

© Habrahabr.ru