В статье ««Магическая константа» 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 = X^P

Показатель степени 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 для того, чтобы на последнем шаге объединить результаты X^P = X^c×X^dp.

Целая степень X^c вычисляется в простом цикле for (int n = 1; n <= c; n++), и обсуждать здесь особенно нечего. Всё «волшебство» тех самых «магических чисел» начинается именно при вычислении дробной части степени X^dp.

И здесь мы сразу выполняем всё то колдунство, описанное в статье, которую я взял за основу для текущей. Алгоритм поиска первого приближения взят именно оттуда, а константу MN = 0x3f7a3bea, базирующуюся на σ = 0,0450465 я заменил на MN = 0x3f7a3с59, основанную на σ = 0,0450333. Где-то на просторах Интернета я прочёл, что так будет лучше, так что оставлю этот вопрос народу для весёлых математических холиваров.

Итак, первое приближение для X^dp выглядит следующим образом:

    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 имеем значение первого приближения X^dp.

Следующим этапом необходимо выполнить уточнение результата, для чего я использовал метод Ньютона, аналогично авторам прочих аналогичных алгоритмов.

Для этого пришлось вывести итерационную формулу для нашего случая.

Итак, исходную функцию Y = X^P перепишем в виде Y^1/P – X = 0. Кстати, уже отсюда следует, что случай с P = 0 следует обрабатывать до достижения этого участка кода. Таким образом, f(Y) = Y^1/P – X, и f’(Y) = (1/P)×Y^(1/P-1).

Итерационная формула рассчитывается как Y_n+1 = Y_n – f(Y_n)/f’(Y_n).

Подставив наши значения для 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 = X^P требуется вычислить g(Y) = ty^(1/dp).

Если вычислить q = 1/dp, то g(Y) = ty^q, и алгоритм переходит в рекурсию.

И на данном этапе, следует вспомнить цель всего мероприятия – быстрое вычисление приближённого значения степени числа, а это значит, что вместо выполнения длительных малоэффективных рекурсий достаточно остановиться на вычислении одного приблизительного значения ty^q аналогично тому, как проводилось вычисление X^P.

Разложим q на целую c и дробную inv часть. Целую степень ty^c вычислим, как и ранее, в цикле, а первое приближение дробной части степени ty^inv снова пользуясь «магическим числом»:

    y.iVal = (int)(MN + inv * (y.iVal - MN));
    inv = y.fVal;

В inv имеем первое приближение значения дробную часть степени 1/dp, на чём я и остановился, предположив, что дальнейшие уточнения сильно замедлят алгоритм, при этом большой выгоды не дадут. Проверку этого предположения я оставил на будущие тесты готового алгоритма, которые следуют в статье далее.

Остаётся только собрать воедино значение ty^(1/dp), далее вычислить по формуле (1) значение X^dp, и дале искомое X^P.

Однако, реализуя этот алгоритм возникла сложность. Вышеупомянутые циклы вычисления целых частей степеней X^P и ty^q действительны только для P > 0, так как используют действие умножение. Для вычисления тех же неизвестных при P < 0 требует либо перемножения обратных величин, либо деления. Во-первых, это усложняет алгоритм, разделяя его уже на начальном этапе на два отдельных участка. Во-вторых, использование деления в цикле серьёзно замедлит выполнение алгоритма для случаев с P < 0.

Я принял решение, искать X^c и X^dp с 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;
    }

Лучшее враг хорошего?

Ещё на этапе написания основной функции у меня появилась идея выделить участок вычисления значения ty^inv в отдельную функцию, подобную функции быстрого вычисления обратного квадратного корня. Никаких проверок исходных данных, «плохие» результаты для расчёта значений за границей области допустимых значений и области применения, зато гораздо выше скорость вычисления. По окончании написания основной функции эта идея была реализована. Таким образом, к этапу тестирования вышли две функции – основная 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 значений на каждый порядок интервала, я получил результаты вычисления функций Y_p = FastPower(Х, 0.75) и Y = Math.Pow(Х, 0.75), рассчитал величины относительной ошибки вычисления 100*|Y_p – 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, и перед началом вычисления X^dp проверку 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.

 Из полученных данных видно, что стандартная функция 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, хотя последний и обладает большей стабильностью скорости выполнения при любых значениях своих параметров.

В любом случае, решение об использования разработанных функций остаётся на усмотрение разработчиков, которые сами оценивают соответствие этих функций конкретной задаче или требованиям проекта.

Ссылки

1. «Магическая константа» 0x5f3759df

2. Алгоритмы быстрого возведения в степень