В статье ««Магическая константа» 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 задана следующим образом:
[StructLayout(LayoutKind.Explicit)]
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:
[StructLayout(LayoutKind.Explicit)]
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, хотя последний и обладает большей стабильностью скорости выполнения при любых значениях своих параметров.
В любом случае, решение об использования разработанных функций остаётся на усмотрение разработчиков, которые сами оценивают соответствие этих функций конкретной задаче или требованиям проекта.
