Comments 40
double[] powtable = {1,2,4,8,16,32, ..., 8.98846567431157954e307};
double pow2(double x)
{
bool sign = x < 0.0;
x = sign ? -x : x;
int xf = (int)x;
x -= xf;
double yt = powtable[xf];
double y = x * (0.6931471805599465398 + x * (0.2402265069590422222 +
x * (0.05550410866590987068 + x * (0.00961812909717555934 +
x * (0.0013333558738165095559 + x * (0.00015403509189194102748 +
x * (0.000015253232908458899497 + x * (1.3207676270599404858e-6 +
x * (1.0258347084283025531e-7 + (6.537941907237267033e-9 +
6.3026908837748924689e-10 * x) * x)))))))));
return sign ? 1 / (y * yt + yt) : y * yt + yt;
}UPD: Если ограничиться точностью single, то будет в 6 раз быстрее.
double y = 2.5868889777e-9 + x * (0.693146928693029174 +
x * (0.240230502044993504 +x * (0.0554804263257944027 +
x * (0.00968458045229538239 + (0.0012387821479302309 +
0.000218775047688226791 * x) * x))));int xf = (int)x;не помешает добавить строчку
if (xf >= 1024) return sign ? 0.0 : double.PositiveInfinity;
Пишите статью "Уточняем pow" :-)
double) просто сдвиг показателя степени на 52 бита влево? Для этого не нужно FPU.Цитируемый код Q_rsqrt неполон без авторских комментариев, сохраняющих актуальность и при возведении в произвольную степень - "evil floating point bit level hacking" и "what the fuck?"
Странно заниматься ускорением функции pow без разбора exp и log.
У бинарного возведения в степень есть погрешность, потому что умножение не ассоциативно в поле чисел с плавающей точкой. x * x * x * x это не то же самое, что y = x * x, y * y.
Кстати, есть более эффективные небинарные возведения, например x ^ 6 — это y = x * x * x, y * y. В бинарном воведении переменных будет на 1 больше: y = x * x, z = y * y, z * y.
Число переменных может влиять на эффективность. Они могут не влезть во все регистры, особенно актуально для больших степеней.
Более того, число умножений тоже может быть меньше, например для x^27:
Тернарное умножение: 6 умножений, 2 переменных:
y = x * x * x // 3
z = y * y * y // 9
result = z * z * z // 27Бинарное умножение: 7 умножений, 6 переменных:
x1 = x * x // 2
x2 = x1 * x1 // 4
x3 = x2 * x2 // 8
x4 = x3 * x3 // 16
x5 = x4 * x3 // 24
x6 = x5 * x1 // 26
result = x6 * x // 27
Интересный вопрос. Наверняка есть, подозреваю, что можно использовать динамическое программирование. Но такое лучше спрашивать у математиков, алгоритмистов. Попробую подумать на досуге.
Формулы нет, оптимальная схема возведения в степень находится только поиском. Во всяком случае так было у Кнута, и я не слышал новостей на эту тему.
Спасибо, что написали про ассоциативность, добавлю пометку.
Однако на практике погрешности не возникает. Видимо, встроенные методы тоже подвержены этому эффекту (на скриншотах ниже погрешность отображается в столбце Sum difference).
Тесты в C++
Тесты в C#
Тесты в Java
Про способ с возведением через кубы я тоже думал. Было бы интересно посмотреть на программную реализацию троичного возведения в степень и сравнить его с другими методами.
Однако на практике погрешности не возникает. Видимо, встроенные методы тоже подвержены этому
А если сравнить с тупым перемножением? Его все равно не помешало бы добавить в таблицу для сравнения.
Про способ с возведением через кубы я тоже думал. Было бы интересно посмотреть на программную реализацию троичного возведения в степень и сравнить его с другими методами.
Там не совсем кубы, а в теории вообще обобщенное количество, которое на больших степенях может давать лучший результат. Но вроде такое разложение NP полное — его имеет смысл использовать разве что в генерации кода, но не в исполнении.
А если сравнить с тупым перемножением? Его все равно не помешало бы добавить в таблицу для сравнения.По идее, у «тупого перемножения» больше погрешность, чем у метода деления степени пополам, потому что намного больше операций (умножений), каждая из которых выполняется с погрешностью.
Погрешность у них одинаковая, как это ни странно. Возведение в квадрат погрешность удваивает. Умножение на точно известное число добавляет к погрешности маленькую константу. Можно показать что погрешность почти полностью определяется показателем степени, а порядок умножения влияет мало.
<ошибся веткой>
Алгоритм: Бинарное возведение в степень ...
Погрешность: нет
Я бы не стал так категорично...
$ echo '
#include <stdio.h>
#include <assert.h>
#include <math.h>
double BinaryPower(double b, unsigned long long e) {
double v = 1.0;
while (e) {
if (e & 1) v *= b;
b *= b;
e >>= 1;
}
return v;
}
int main()
{
unsigned e = 134217728;
double barr[] = {1.000001, 1.0000001, 1.00000001, 0}, *b = barr;
while (*b) {
printf("calc %.10g ** %u:\n fast pow ==> %.16f\n native ==> %.16f\n",
*b, e, BinaryPower(*b, e), pow(*b,e));
b++;
}
return 0;
}
// ' > test.c; gcc -O2 -Wall -Wextra test.c -o test; ./test
calc 1.000001 ** 134217728:
- fast pow ==> 19497974326417947731259943008814618865863786077838628618240.0000000000000000
+ native ==> 19497974384856317387011039256973190612600604089375799640064.0000000000000000
calc 1.0000001 ** 134217728:
- fast pow ==> 674530.4760286132805049
+ native ==> 674530.4755217875353992
calc 1.00000001 ** 134217728:
- fast pow ==> 3.8273676342423024
+ native ==> 3.8273676116718129Как видим ошибка может быть довольно значительной, я уж не говорю про -Ofast (без -fno-fast-math).
Поправка: пример компилился в -m32toolchain... для -m64 погрешности действительно нет (по крайней мере на этих числах).
1.94979745417325751802166478351188792876322264823289192517468 × 10^58
674530.47074108455938268917802974681284444414341034203174237732783
3.8273676654778624379533125217738601044713095101002538124954655334
Т.е. разница между двумя реализациями
pow не столь значительна, как разница между ними и точным результатом. В последнем случае «быстрая» реализация ещё и ближе к точному результату, чем «родная».Ну у Wolfram precision и accuracy настраиваемые (если не ошибаюсь оно умеет длинную FP-арифметику, exact quantities и всё такое)...
Пример алгоритма был для double (aka 64 bits IEEE normalized double-precision floating-point number) и C/C++...
Я это о чем, собственно - что-то математически правильное (и доказуемо верное) в конкретной реализации на конкретном языке для какой-либо платформы может вылиться в неслабую такую погрешность (из-за переполнений, недостаточной precision промежуточных результатов и т.п.)...
Я про другое: что «родная»
pow даёт столь же неточный результат.Несомненно, что бинарное возведение в степень на практике даёт неточный результат.
Ну дак а я о чем, "Погрешность: нет" - как бы не совсем верно, если про этот алгоритм...
Я про другое: что «родная»
powдаёт столь же неточный результат.
В рамках double - конечно.
Я вам больше скажу, например тот же bigfloat.pow (даже с precision 400) также выдает "неточный" результат (и он менее точен чем нативный pow):
bigfloat.pow, precision 400 ...
>>> pow(1.000001, 134217728, precision(400))
BigFloat.exact('19497974326443151952384768395025354130242292944980784774315.431771682120918607912604365932634709150271941538088343395908095', precision=400)
>>> pow(1.0000001, 134217728, precision(400))
BigFloat.exact('674530.47602706407027793197206421845263783142259940138702038470170721140355151884268337259494331459943481248761388507182459', precision=400)
>>> pow(1.00000001, 134217728, precision(400))
BigFloat.exact('3.8273676342578585041456476752437466436493872281452027009938444315573508186002376274873069800877205918834944719451685524823', precision=400)Для сравнения возьмем корни от обоих (этот алгоритм у bigfloat насколько знаю точнее чем pow):
bigfloat.root в сравнении с Wolfram ...
>>> setcontext(precision(220)); p = getcontext().precision;
>>> root(pow(1.000001, 134217728), 134217728)
BigFloat.exact('1.0000009999999999177333620536956004798412322998046875000000000000000', precision=220)
>>> root(BigFloat.exact('1.94979745417325751802166478351188792876322264823289192517468e58', precision=p), 134217728)
BigFloat.exact('1.0000010000000000000000000000000000000000000000000000000000000000002', precision=220)
>>>
>>> root(pow(1.0000001, 134217728), 134217728)
BigFloat.exact('1.0000001000000000583867176828789524734020233154296875000000000000000', precision=220)
>>> root(BigFloat.exact('674530.47074108455938268917802974681284444414341034203174237732783', precision=p), 134217728)
BigFloat.exact('1.0000001000000000000000000000000000000000000000000000000000000000003', precision=220)
>>>
>>> root(pow(1.00000001, 134217728), 134217728)
BigFloat.exact('1.0000000099999999392252902907785028219223022460937500000000000000000', precision=220)
>>> root(BigFloat.exact('3.8273676654778624379533125217738601044713095101002538124954655334', precision=p), 134217728)
BigFloat.exact('1.0000000100000000000000000000000000000000000000000000000000000000001', precision=220)Но если округлить результаты до оригинальной accuracy - то все результаты верны.
Что не отменяет факт, что результат pow может сильно отличатся и погрешность есть, хоть и небольшая - для bigfloat.pow (и BinaryPower) около 11.04e-7%, а для pow - 08.04e-7%, что очень неплохо для степени 134217728.
│ algorithm │ error % │
├──────────────┼───────────────┤
│ BinaryPower | 0.0000011043% │
│ bigfloat.pow │ 0.0000011042% │
│ native pow │ 0.0000008045% │
-m32 означает что все промежуточные результаты приводятся к float. Не мудрено получить погрешность.
Нет, не приводятся - double он и в Африке double, и что там, что тут - 64 бит (размеры floating-point data types не зависит от разрядности и одинаковы и для x86 и для x64). При этом даже для x86, внутренние FPU регистры - 80 бит. Различие может быть лишь в оптимизации, да в используемых инструкциях (mulsd vs fmul).
А можно ли улучшить "грубый" результат каким-либо алгоритмом уточнения корня (как это делается для быстрого обратного квадратного корня)?
Боюсь, что уточнение результата так же, как и в алгоритме быстрого обратного корня, для произвольной степени будет не оптимально. В быстром обратном корне используется метод Ньютона, который использует производную от изначальной функции. А у нас функция вида 
, так что её производная будет 
, что также содержит в себе степень. Получается, что нам придется помимо того, что возводить
в степень
, также возводить
в степень
, для чего снова придется использовать один из быстрых алгоритмов, который будет снова использовать себя, и так далее, пока не дойдем до какой-то степени, которую можно посчитать за один такт.
По поводу других методов уточнения результата ничего сказать не могу, не изучал эту тему.
Ускоряем pow