В статье ««Магическая константа» 0x5f3759df» проводится подробный анализ алгоритма быстрого вычисления обратного квадратного корня, изучение вопроса появления «магической константы», а также её модификации в целях дальнейшего использования в алгоритме вычисления чисел в любой степени. Была получена новая «магическая константа» и обозначено направление для разработки алгоритма с её использованием.
В этой статье я покажу алгоритм функции вычисления любой степени положительного числа, использующий новую «магическую константу». Кроме того, я приведу результаты её сравнения с исходной функцией вычисления обратного квадратного корня, а также со стандартной функцией вычисления степени Math.Pow.
На затравку
Алгоритм вычисления степени можно разделить на пять основных частей:
— предварительные проверки начальных условий;
— вычисление целой части значения степени заданного аргумента;
— вычисление первого приближения дробной части значения степени;
— уточнение результата дробной части значения степени;
— перемножение результатов вычисления целой и дробной частей степени для получения окончательного результата.
Рассмотрим исходную функцию вычисления обратного квадратного корня.
float FastInvSqrt(float x) { float xhalf = 0.5f * x; int i = *(int*)&x; // представим биты float в виде целого числа i = 0x5f3759df - (i >> 1); // какого черта здесь происходит ? x = *(float*)&i; x = x*(1.5f-(xhalf*x*x)); return x; }
Вероятно, для ускорения алгоритма, предварительных проверок значения аргумента в функции не предусмотрено. Хотя, строго говоря, область допустимых значений функции включает строго положительные числа, а математическая основа дальнейшего использования «магического числа» 0x5f3759df требует дополнительного ограничения для аргумента не превышать значение +1.
Вычисления целой части степени в алгоритме не предусмотрено, за неимением таковой для P = -1/2.
Первое приближение вычисляется с использованием вышеупомянутого «магического числа».
И наконец, уточнение результата проводится с использованием метода Ньютона. Генерация результата здесь происходит на этой стадии, за неимением целой части значения степени P, как говорилось выше.
Собираем «свой Лунапарк»
Приступим к построению функции вычисления степени. Для этого определимся с требованием к исходным данным. Искомая математическая функция имеет вид
Y = XP
Показатель степени P может принимать любое значение из множества действительных чисел R. Однако, область допустимых значений аргумента Х зависит от P. При P, принадлежащему множеству целых чисел Z, Х Є (-∞; +∞). Но для прочих P из множества R, область допустимых значений X сужается до (0; +∞). Результатом для всех сочетаний аргументов и показателей степени, не попадающих в область допустимых значений, будет float.NaN – программный аналог математической неопределённости.
Кроме того, при P = 0 для всех действительных Х результатом будет 1, и дальнейшие расчёты не имеют смысла. Аналогично, при P = 1 результатом будет X.
Отдельно следует отметить значение аргумента X = 0. В этом случае, для любых P < 0 значение функции является вопросом далеко не тривиальным. В дебрях матана и функана легко потерять исходный посыл своих действий, по этой причине примем результатом возведения 0 в отрицательную степень значение float.PositiveInfinity, то есть +∞.
При X = 0 для P > 0 дела несколько проще, и здесь результатом примем 0.
Учитывая вышеизложенное, я поставил следующие условия в начале (и одно в теле) функции:
— результат 1 при P = 0;
— результат X при P = 1;
— результат float.PositiveInfinity при X = 0 и P < 0;
— результат 0 при X = 0 и P > 0;
— результат float.NaN при X < 0 для всех P, не принадлежащих множеству целых чисел Z.
В ходе реализации алгоритма расчёта значения функции, проверку последних двух условий я реализовал отдельно для P > 0 и P < 0. Об этом я расскажу далее.
Следующим этапом, разделим степень P на целую часть c и дробную часть dp для того, чтобы на последнем шаге объединить результаты XP = Xc×Xdp.
Целая степень Xc вычисляется в простом цикле for (int n = 1; n <= c; n++), и обсуждать здесь особенно нечего. Всё «волшебство» тех самых «магических чисел» начинается именно при вычислении дробной части степени Xdp.
И здесь мы сразу выполняем всё то колдунство, описанное в статье, которую я взял за основу для текущей. Алгоритм поиска первого приближения взят именно оттуда, а константу MN = 0x3f7a3bea, базирующуюся на σ = 0,0450465 я заменил на MN = 0x3f7a3с59, основанную на σ = 0,0450333. Где-то на просторах Интернета я прочёл, что так будет лучше, так что оставлю этот вопрос народу для весёлых математических холиваров.
Итак, первое приближение для Xdp выглядит следующим образом:
FIBUnion y = new FIBUnion(); y.fVal = x; y.iVal = (int)(MN + dp * (y.iVal - MN)); float ty = y.fVal;
где структура FIBUnion задана следующим образом:
public struct FIBUnion { [FieldOffset(0)] public float fVal; [FieldOffset(0)] public int iVal; [FieldOffset(0)] public bool bVal; }
Она позволяет обращаться к одной и той же ячейке памяти размером 32 бита как к разным типам данных – float, int или bool.
В итоге, в ty имеем значение первого приближения Xdp.
Следующим этапом необходимо выполнить уточнение результата, для чего я использовал метод Ньютона, аналогично авторам прочих аналогичных алгоритмов.
Для этого пришлось вывести итерационную формулу для нашего случая.
Итак, исходную функцию Y = XP перепишем в виде Y1/P – X = 0. Кстати, уже отсюда следует, что случай с P = 0 следует обрабатывать до достижения этого участка кода. Таким образом, f(Y) = Y1/P – X, и f’(Y) = (1/P)×Y(1/P-1).
Итерационная формула рассчитывается как Yn+1 = Yn – f(Yn)/f’(Yn).
Подставив наши значения для f(Y) и f’(Y), и упростив выражение, получим рабочую версию формулы для уточнения результата:
Res = Y×(1+P×(X/Y(1/P)-1)).
Исходными данными для этой формулы будут Y = ty, и P = dp, и окончательное выражение примет вид:
Res = ty×(1+dp×(X/ty(1/dp)-1)) (1)
И здесь можно заметить забавный факт, что для поиска значения Y = XP требуется вычислить g(Y) = ty(1/dp).
Если вычислить q = 1/dp, то g(Y) = tyq, и алгоритм переходит в рекурсию.
И на данном этапе, следует вспомнить цель всего мероприятия – быстрое вычисление приближённого значения степени числа, а это значит, что вместо выполнения длительных малоэффективных рекурсий достаточно остановиться на вычислении одного приблизительного значения tyq аналогично тому, как проводилось вычисление XP.
Разложим q на целую c и дробную inv часть. Целую степень tyc вычислим, как и ранее, в цикле, а первое приближение дробной части степени tyinv снова пользуясь «магическим числом»:
y.iVal = (int)(MN + inv * (y.iVal - MN)); inv = y.fVal;
В inv имеем первое приближение значения дробную часть степени 1/dp, на чём я и остановился, предположив, что дальнейшие уточнения сильно замедлят алгоритм, при этом большой выгоды не дадут. Проверку этого предположения я оставил на будущие тесты готового алгоритма, которые следуют в статье далее.
Остаётся только собрать воедино значение ty(1/dp), далее вычислить по формуле (1) значение Xdp, и дале искомое XP.
Однако, реализуя этот алгоритм возникла сложность. Вышеупомянутые циклы вычисления целых частей степеней XP и tyq действительны только для P > 0, так как используют действие умножение. Для вычисления тех же неизвестных при P < 0 требует либо перемножения обратных величин, либо деления. Во-первых, это усложняет алгоритм, разделяя его уже на начальном этапе на два отдельных участка. Во-вторых, использование деления в цикле серьёзно замедлит выполнение алгоритма для случаев с P < 0.
Я принял решение, искать Xc и Xdp с c и dp , взятыми по модулю, а затем, при P < 0, взять обратное значение полученного результата.
Это, не отменило необходимость раздельного вычисления по двум ветвям для P < 0 и P > 0, но передвинуло его ближе к окончанию алгоритма. Эти две ветви имеют протяжённые повторяющиеся участки кода, однако, выделять их в отдельную функцию, или обкладывать условными переходами я не стал, исходя из принципа «скорость > размер».
В ходе оптимизации кода и поиска возможностей для ускорения вычислений я наткнулся на статью Алгоритмы быстрого возведения в степень. Из приведённых вариантов я выбрал для применения бинарный алгоритм вычисления целой степени, заменив использованные изначально циклы умножения. Так как операции умножения чисел с плавающей запятой являются весьма затратными по процессорному времени, я скорректировал приведённый в статье бинарный алгоритм для дополнительного уменьшения операций умножения. Для степеней 0 – 5 разница скорости вычисления несущественна, но с ростом степени преимущество бинарного алгоритма резко возрастает. Полученный алгоритм показал более чем двукратное увеличение скорости вычисления степени 15 произвольных действительных чисел в сравнении с методом цикла умножений.
После проведения оптимизации кода в целях увеличения скорости его выполнения свет увидела функция Power(float x, float p):
// быстрое вычисление степени p аргумента х c проверками аргументов public static float Power(float x, float p) { const float eLim = 0.01f; // наименьшая степень с допустимой ошибкой const int MN = 0x3f7a3c59; // проверки начальных условий if (p == 0) return 1; // уточнение для p = 0 if (p == 1) return x; // уточнение для p = 1 //if (x < 0) return float.NaN; // для целых p условие не верно, см. далее if (x == 0) // корректировка для x = 0 { if (p < 0) return float.PositiveInfinity; else return 0; } float dp; // c + dp = |p|, где c - целое, а dp - дробный остаток int c = (int)p; if (p < 0) { dp = c - p; // dp >= 0 c = ~(--c); // |c| >= 0 } else dp = p - c; // dp, c >= 0 //if (dp != 0 && x < 0) return float.NaN; // для ускорения будет внутри if, см. далее // yn^|p| = yn^dp*yn^c // вычислим fRes = yn^c // через циклы умножений /*float fRes = 1; for (int n = 1; n <= c; n++) { fRes *= x; }*/ // ускоренное вычисление float fRes; float inv = x; if ((c & 1) == 0) // проверка чётности { fRes = 1; } else { fRes = inv; } c = c >> 1; while (c != 0) { inv *= inv; if ((c & 1) == 1) { fRes *= inv; } c = c >> 1; } // окончание вычисления целой части степени if (p < 0) { //if (dp == 0) return 1 / fRes; // дальше вычислять нечего if (dp < eLim) return 1 / fRes; // для снижения отклонения if (x < 0) return float.NaN; // для нецелых p при x < 0 // вычислим yn^dp, алгоритм см. FastPower float inv = 1 / dp; c = (int)inv; inv = inv - c; FIBUnion y = new FIBUnion(); y.fVal = x; y.iVal = (int)(MN + dp * (y.iVal - MN)); float ty = y.fVal; y.iVal = (int)(MN + inv * (y.iVal - MN)); inv = y.fVal; /*for (int n = 1; n <= c; n++) { inv *= ty; }*/ float v = ty; if ((c & 1) == 1) { inv *= v; } c = c >> 1; while (c != 0) { v *= v; if ((c & 1) == 1) { inv *= v; } c = c >> 1; } inv = x / inv; return 1 / (fRes * ty * (1 + (dp * (--inv)))); } else { // аналогично для p>0 //if (dp == 0) return fRes; if (dp < eLim) return fRes; if (x < 0) return float.NaN; float inv = 1 / dp; c = (int)inv; inv = inv - c; FIBUnion y = new FIBUnion(); y.fVal = x; y.iVal = (int)(MN + dp * (y.iVal - MN)); float ty = y.fVal; y.iVal = (int)(MN + inv * (y.iVal - MN)); inv = y.fVal; /*for (int n = 1; n <= c; n++) { inv *= ty; }*/ float v = ty; if ((c & 1) == 1) { inv *= v; } c = c >> 1; while (c != 0) { v *= v; if ((c & 1) == 1) { inv *= v; } c = c >> 1; } // Формула метода Ньютона немного отличается для p>0 inv = x / inv; return fRes * ty * (1 + (dp * (--inv))); } }
И не забудьте структуру FIBUnion:
public struct FIBUnion { [FieldOffset(0)] public float fVal; [FieldOffset(0)] public int iVal; [FieldOffset(0)] public bool bVal; }
Лучшее враг хорошего?
Ещё на этапе написания основной функции у меня появилась идея выделить участок вычисления значения tyinv в отдельную функцию, подобную функции быстрого вычисления обратного квадратного корня. Никаких проверок исходных данных, «плохие» результаты для расчёта значений за границей области допустимых значений и области применения, зато гораздо выше скорость вычисления. По окончании написания основной функции эта идея была реализована. Таким образом, к этапу тестирования вышли две функции – основная Power(float x, float p) и её упрощённый и быстрый вариант FastPower(float x, float p):
// быстрое вычисление степени a Э (0, 1) аргумента х > 0 без проверок public static float FastPower(float x, float a) { /* "магическое" число MN = L(B-s) для типа "float" х = [(-1)^sign]*[2^(E-B)]*M, где int Е - 8 бит экспоненты, B = 127, int M - 23 бита мантиссы, L = 2^23 если представить как х = (1+m)2^e, где e = E-B, а m = M/L а его целочисленная интерпретация i = M+LE настоящее магическое число s = 0.0450465 заменено более удачным s = 0.0450333 */ const int MN = 0x3f7a3c59; // до замены s было MN = 0x3f7a3bea /* метод Ньютона. yn+1 = yn - f(yn)/f'(yn) y = x^a; f(y) = y^(1/a)-x; f'(y) = (1/a)*1/(y^(1/a-1)) yn+1 = yn*(1+a*(x/yn^(1/a)-1)) */ // найдём yn^(1/a) float inv = 1 / a; int c = (int)inv; inv = inv - c; // теперь c + inv = 1/a, где c - целое, а inv - дробный остаток // первое приближение в ty FIBUnion y = new FIBUnion(); y.fVal = x; y.iVal = (int)(MN + a * (y.iVal - MN)); float ty = y.fVal; // yn^(1/a) = yn^inv*yn^c // вычисляем приблизительно yn^inv y.iVal = (int)(MN + inv * (y.iVal - MN)); inv = y.fVal; // вычисляем yn^c и соберём yn^(|1/a|) => inv // через циклы умножений /*float fRes = 1; for (int n = 1; n <= c; n++) { fRes *= x; }*/ // ускоренное вычисление float fRes; float inv = x; if ((c & 1) == 0) // проверка чётности { fRes = 1; } else { fRes = inv; } c = c >> 1; while (c != 0) { inv *= inv; if ((c & 1) == 1) { fRes *= inv; } c = c >> 1; } // окончание вычисления целой части степени //yn+1 = yn*(1+a*(x/yn^(1/a)-1)) inv = x / inv; return ty * (1 + (a * (--inv))); }
Это кто тут такой лохматенький?
Начнём выполнять проверки, исходя из принципа программирования «рабочий код -> быстрый код -> красивый код». Следовательно, первым этапом следует оценить результаты выполнения разработанных функций. Естественно, корректность результатов я проверял ещё на этапе разработки, однако, так как функции вычисляют приближённое значение, следует исследовать вопрос получаемых ими отклонений.
Вычисление целой части степени производится с точностью, определённой для операции умножения чисел типа float. Основной вклад в возникающее в процессе вычисления отклонение обусловлен выбранным способом вычисления дробной части степени. Эта часть алгоритма идентична для обеих разработанных функций. По этой причине оценку отклонений вычислений я проводил только для функции FastPower.
Для тестирования я генерировал значения аргумента из интервала (0; +1) в массив для дальнейшего вычисления функций, значения которых сохранялись в другой массив:
float[] fX = new float[iIt]; float[] fY = new float[iIt]; Random Rnd = new Random(iRnd); // инициализация массива аргументов for (int i = 0; i < iIt; i++) { fX[i] = ((float)Rnd.NextDouble())*flEnd+flBeg; }
Контрольные значения, относительно которых вычислялись отклонения вычислений я получал стандартной функцией Math.Pow.
Для 4 000 случайных значений Х в интервале (+0,0001; +1), по 1 000 значений на каждый порядок интервала, я получил результаты вычисления функций Yp = FastPower(Х, 0.75) и Y = Math.Pow(Х, 0.75), рассчитал величины относительной ошибки вычисления 100*|Yp – Y|/Y, в процентах, и построил диаграмму зависимости относительной ошибки вычисления функции FastPower от значения аргумента (рис. 1):
Из диаграммы очевидна периодичность относительной ошибки вычисления, по этой причине, я разбил диаграмму на три участка по одному порядку в интервале (+0,0001; +1) и вычислил смещение, требуемое для совпадения этих участков. Результат представлен на диаграмме (рис. 2):
На диаграмме легко видеть, что зависимость относительной ошибки функции FastPower от X при P = 0,75 логопериодическая, с логарифмическим периодом 16.
Аналогичным образом, я построил диаграмму зависимости относительной ошибки функции FastPower от X при P = 0,55 (рис. 3):
Очевидно, что для степени 0,55 периодичности относительной ошибки не наблюдается. Характер диаграммы сильно отличается от степени 0,75.
В дальнейшем, собрав данные для степеней от 0,05 до 0,95 с шагом 0,1, я построил общую диаграмму, с целью попытаться определить динамику изменения зависимости относительной ошибки от значения степени P (рис. 4):
Диаграммы непрерывны, хотя и содержат в себе множество особых точек, в которых происходят изломы. Однако, даже для «соседних» взятых значений P диаграммы резко отличаются друг от друга. На этом этапе появилось предположение, что зависимость относительной ошибки от значения степени носит характер, схожий с её зависимостью от аргумента.
С целью проверки этого предположения, для 4 000 случайных значений P в интервале (0; +1), я получил результаты вычисления функций FastPower, при Х = 0,05, Х = 0,30, Х = 0,55 и Х = 0,95, и построил диаграмму зависимости относительной ошибки вычисления функции FastPower от значения аргумента (рис. 5):
Полученные диаграммы подтвердили и сделанное предположение о характере зависимости относительной ошибки вычисления функции FastPower от значения степени, а также обозначили «слабое место» в алгоритме разработанной функции. Из диаграмм рис. 1 – 5 видно, что в большинстве случаев относительная ошибка вычисления не превышает 3%. Однако, в области малых значений Р, она резко возрастает, и, очевидно, стремится к +∞ (рис. 5). Сложившуюся ситуацию можно оставить без внимания для быстрой небезопасной функции FastPower, но в функции Power необходимо предусмотреть защитный механизм. Для этого я объявил константу eLim = 0.01f, и перед началом вычисления Xdp проверку dp = 0 заменил на dp < eLim. Таким образом, в области малых dp программа проводит вычисления как для dp = 0, и не допускает роста ошибки до бесконечности.
Весёлые старты
Настало время узнать, ради чего проделана вся работа. Для измерения времени выполнения функций я использовал Stopwatch из пространства имен System.Diagnostics C#. Для сравнения я использовал функцию Math.Pow. Каждый тест состоял из циклов, в каждом из которых генерировался массив из аргументов для вычислений, случайно выбранных из заданного диапазона. Из всех циклов выбирался цикл с наименьшим значением времени выполнения (самый быстрый). Для тестов с разными функциями на одном и том же диапазоне, в целях создания максимально близких условий, генерация случайных чисел проводилась с одним и тем же начальным значением. Измерение времени выполнения функции, на примере Power организовано в виде следующего алгоритма:
Stopwatch sw = new Stopwatch(); Random Rnd = new Random(iRnd); for (int cc = 0; cc<iCyc; cc++) { sw.Reset(); // инициализация массива аргументов for (int i = 0; i < iIt; i++) { fX[i] = ((float)Rnd.NextDouble())*flEnd+flBeg; } sw.Start(); for (int i = 0; i < iIt; i++) { fY[i] = Func.Power(fX[i], flP); } sw.Stop(); }
Для тестов с аргументом в диапазоне Х > 0 я также использовал функции FastPower, и FastInvSqr для быстрого вычисления обратного квадратного корня. Все результаты получены со Stopwatch в режиме использования счетчика производительности высокого разрешения (IsHighResolution == True), и выражаются в тактах таймера (ElapsedTicks). Получены следующие результаты измерения (таб. 1):
Таблица 1. Количество циклов 300, количество вычислений в цикле 2000.
№ пп |
Диапазон Х |
P |
Math.Pow |
Power |
FastPower |
FastInvSqr |
1 |
(-2; +2) |
-2,75 |
363 |
295 |
|
|
2 |
(-10; +10) |
-2,75 |
356 |
296 |
|
|
3 |
(-2; +2) |
2,75 |
363 |
273 |
|
|
4 |
(-2; +2) |
-15,13 |
358 |
493 |
|
|
5 |
(0; +2) |
-15,13 |
359 |
613 |
|
|
6 |
(0; +2) |
15,13 |
359 |
583 |
|
|
7 |
(-5; -1) |
-15,13 |
362 |
150 |
|
|
8 |
(0; +2) |
-0,11 |
472 |
482 |
|
|
9 |
(0; +2) |
-0,83 |
391 |
361 |
|
|
10 |
(0; +2) |
0,11 |
473 |
453 |
326 |
|
11 |
(0; +2) |
0,83 |
392 |
327 |
204 |
|
12 |
(0; +2) |
0,01 |
480 |
1839 |
1713 |
|
13 |
(0; +2) |
0,009 |
480 |
137 |
1885 |
|
14 |
(0; +2) |
-0,5 |
422 |
372 |
188 |
91 |
Из полученных данных видно, что стандартная функция Math.Pow отличается довольно высокой стабильностью по скорости исполнения, слабо зависящей от величины и знака как аргументов, так и степеней. Заметное замедление наблюдается в области вычисления малых степеней, близких к нулю. Изучение скорости выполнения функции Math.Pow не является предметом данной статьи, и этот вопрос рассматривается исключительно в рамках сравнения со скоростями исполнения других функций.
На основании пар пп. 2, 3; 5, 6; 8, 10 и 9, 11 можно видеть, что функция Power выполняется на (5 ÷ 10)% быстрее для положительных степеней, что, вероятно, объясняется наличием действий деления в алгоритме обработки отрицательных степеней.
Результаты пп. 4, 5 и 7 показывают, что быстрая обработка случаев с X < 0 при P не принадлежащем множеству целых чисел, в результате которой функции выдают результат float.NaN ускоряет выполнение функции Power на 50% и более. Функция Math.Pow же, очевидно, обрабатывает такие случаи с той же скоростью, что и прочие.
Пары результатов пп. 8,9 и 10, 11, а также п. 12 показывают замедление алгоритма функций Power и FastPower вследствие выполнения циклов умножения при расчёте дробной части степени. Замедление достигает 5,6 раза и ограничено только условием dp < eLim для функции Power. Время же выполнения функции FastPower будет расти и далее, при приближении степени к 0 как видно из п. 13. В случаях, когда целая часть степени больше 1, выполняются циклы умножения в начальной части алгоритма, что заметно при сравнении троек результатов пп. 6, 10, 11 и 5, 8, 9. Таким образом, наличие модуля целой части степени больше 1, и дробной части – меньше 0,5(0)1 замедляет выполнение алгоритма функций Power и FastPower минимум на (27 ÷ 30)% каждое.
И вполне предсказуемо, функция FastInvSqr показывает лучший результат, что показано в п. 14. Скорость её выполнения выше даже случая, когда в функции Power не выполняется почти никаких математических операций (п. 7). Интересно отметить, что скорость выполнения функции Math.Pow в случае п. 14 имеет необычную величину, в сравнении с прочими степенями в интервале (0; +1). Однако, интереснейший вопрос исследования скорости выполнения функции Math.Pow выходит за рамки данной статьи.
В заключение
Реализован алгоритм вычисления степени на основе использования «магического числа» 0x3f7a3c59 в виде функции Power с областью допустимых значений обоих параметров (-∞; +∞), и более быстрой функции FastPower, с областью допустимых значений параметра X (0; +∞), и параметра P (0; +1).
Погрешность вычислений алгоритма, как правило, не превышают 3%.
В области небольших по абсолютной величине, и не слишком близко лежащих к 0 значениях аргумента и показателя степени, функции Power и FastPower показывают хорошую скорость вычисления относительно Math.Pow, хотя последний и обладает большей стабильностью скорости выполнения при любых значениях своих параметров.
В любом случае, решение об использования разработанных функций остаётся на усмотрение разработчиков, которые сами оценивают соответствие этих функций конкретной задаче или требованиям проекта.
Ссылки
ссылка на оригинал статьи https://habr.com/ru/articles/820713/
Добавить комментарий