Как на компьютере вычисляются корни и др мат. функции?

[photon]

Тогда доработанный вариант можно записать так:

Используется синтаксис C++

template <size_t N>
inline float newton_rsqrt(float const d0, float const d1)
{
    float a = d0;
    a *= d1;
    a *= d1;
    a += 1.5f;
    d1 *= a;
    return newton_rsqrt<N - 1>(d0, d1);
}

template <>
inline float newton_rsqrt<0>(float const, float const d1)
{
    return d1;
}


template <size_t N>
inline float rsqrt(float const d0)
{
    union
    {
        float f;
        int32_t i;
    } d1;

    d1.f = d0;
    d1.i = 0x5F3759DF - (d1.i >> 1);
    d0 *= -0.5f;
    return newton_rsqrt<N>(d0, d1.f);
}

template <size_t N>
inline float fast_sqrt(float const x)
{
    return x * rsqrt<N>(x);
}

inline float rsqrt(float const x)
{
    return rsqrt<2>(x);
}

inline float fast_sqrt(float const x)
{
    return x * rsqrt(x);
}

[=SSN=]

photon в сообщении #1481099 писал(а):

Тогда доработанный вариант можно записать так:

А зачем это надо? На современных процессорах sqrt вычисляется аппаратно и без всяких Ньютонов.
По скорости этот метод тоже проигрывает аппаратному sqrt.
Если только для обучения студентов..

[photon]

=SSN= в сообщении #1481104 писал(а):

По скорости этот метод тоже проигрывает аппаратному sqrt.

Вас не затруднит привести результаты замеров скорости?

-- Fri Aug 28, 2020 12:42:48 --

И даже если вы правы насчет скорости, тема посвящена алгоритмам вычисления. Не думаю, что ответ "используйте sqrt()" - в такой ситуации хорош.

Еще быстрее вычисления можно провести, если мы заранее что-то знаем об аргументе. Например, если нас интересуют корни только из чисел в диапазоне $[0, 1]$ - вполне возможны ситуации, когда данные как-то отнормированы. В таком случае, хороший результат дают интерполяционные методы - табулируем какое-то разумное число точек (например $4096$) в нужном диапазоне, и вычисляем результат просто как линейную интерполяцию по двум соседям.

[Dmitriy40]

Влезу с комментом про скорость. Цифры приведу для Haswell, мне так удобнее, а в более новых ничего кардинально вроде бы не изменилось.
Сначала для скаляров. Команда SQRTSS имеет latency 11 тактов. Каждая итерация Ньютона в коде выше требует трёх умножений и сложения, причём все операции линейно зависимы и параллельно выполниться не могут. Сложение учитывать не будем (вдруг оно оптимизируется в команды FMA с той же latency в 5 тактов, но уже суммарно со сложением), но три умножения суммарно дают latency уже 15 тактов на итерацию. Всё, дальше можно не считать, даже одна итерация Ньютона уже медленнее аппаратной команды. Если всё же прикинуть, то 2 итерации плюс финальное умножение плюс быстрый инверсный корень дают примерно 38 тактов latency. Вместо 11.
Теперь векторно и для AVX. Команда VSQRTPS для ymm регистра имеет latency в 19 тактов. Команды VMULPS имеют latency те же 5 тактов, значит примерно 38 тактов так и останутся. Замедление уменьшилось с трёх с лишним до двух раз.
Double точность рассматривать не буду так как двух итераций Ньютона недостаточно для получения соответствующей точности. А достаточные 4 итерации займут около 68 тактов против менее 30 тактов команды VSQRTPD.

Итого, во всех вариантах аппаратная команда вдвое-втрое быстрее.

При этом польза от алгоритмического варианта будет если отсутствует аппаратная команда, или неподдерживаемый формат чисел, или ради обучения. Т.е. далеко не нулевая.

[=SSN=]

Dmitriy40 в сообщении #1481118 писал(а):

Итого, во всех вариантах аппаратная команда вдвое-втрое быстрее.

Да, именно это и имелось ввиду! Спасибо за подробный ответ!

Ну, и для полноты картины, сколько тактов выполняется команда RSQRTSS ?

[Dmitriy40]

=SSN= в сообщении #1481133 писал(а):

Ну, и для полноты картины, сколько тактов выполняется команда RSQRTSS ?

Latency 5 тактов для xmm (и скалярная и векторная) и 7 для ymm (векторная). Сдвигом и вычитанием выгоднее, уложится в два такта, а точность вполне сравнима, 0.04% команда и 0.4% алгоритм.

[=SSN=]

Dmitriy40 в сообщении #1481137 писал(а):

Latency 5 тактов для xmm (и скалярная и векторная) и 7 для ymm (векторная). Сдвигом и вычитанием выгоднее, уложится в два такта, а точность вполне сравнима, 0.04% команда и 0.4% алгоритм.

Я про то, что SQRT() можно также вычислить через RSQRTSS и MULPS. Получится 10 тактов. Т.е., еще один такт можно выиграть.

[GAA]

Корень нельзя вычислить через RSQRTSS и MULPS т.к. RSQRTSS находит приближённое значение, хотя и с высокой относительной точностью: [относительная погрешность] не больше $1.5 \cdot 2^{-12}$ (см. Intel® 64 and IA-32 Architectures Software Developer’s Manual Volume 2).

[Dmitriy40]

=SSN= в сообщении #1481140 писал(а):

Я про то, что SQRT() можно также вычислить через RSQRTSS и MULPS. Получится 10 тактов. Т.е., еще один такт можно выиграть.

Нельзя, команда RSQRTxx чуть ли не единственная имеющая просто ну никакую точность (относительная ошибка не превышает $1.5 \cdot 2^{-12} \approx 0.0366\%$). Если такая точность устраивает (как в компьютерной графике), то ОК, но для обычных расчётов это всего 3 верных цифры, это простите ни о чём.
Другое дело что кажется её точности хватает для лишь одной итерации Ньютона после для получения нормальной single точности (точно анализировать лень) (возможно это и было целью увеличения точности команды по отношению к сдвигу и вычитанию) и тогда команда становится выгоднее сдвига с вычитанием и лишней итерацией Ньютона.

(Чисто техническое замечание)

Вместо MULPS в пару к RSQRTSS надо использовать MULSS: либо обе векторные (PS), либо обе скалярные (SS), смешивать не стоит.

[=SSN=]

OK. Спасибо!

[TOTAL]

=SSN= в сообщении #1481104 писал(а):

А зачем это надо? На современных процессорах sqrt вычисляется аппаратно и без всяких Ньютонов.

Аппаратное вычисление - это не волшебная палочка. Это тоже реализация метода, пусть и не Ньютона.

[=SSN=]

TOTAL в сообщении #1481156 писал(а):

Аппаратное вычисление - это не волшебная палочка. Это тоже реализация метода, пусть и не Ньютона.

Тогда следует повторить исходный вопрос:

IvanX в сообщении #1480468 писал(а):

Как на компьютере вычисляются корни, показательная функция, логарифмы, тригонометрические функции?

И если аппаратное вычисление sqrt не Ньютоном, то тогда чем?

PS. Табличные методы не предлагать, ибо для даблов таблицы будут слишком толстые.. ;)

PPS. А вот вариант с двумя таблицами, похоже, сработает.. :)

[worm2]

=SSN= в сообщении #1481159 писал(а):

PS. Табличные методы не предлагать, ибо для даблов таблицы будут слишком толстые.. ;)

Кстати, для точности в 12 битов таблица — в самый раз, должна в L1 кэш влезать.

Цитата:

PPS. А вот вариант с двумя таблицами, похоже, сработает.. :)

Две таблицы — это как?

[Dmitriy40]

По моему точность в 12 битов результата достигается уже при таблице из 8-ми (рассматриваем только степени двух) точек с линейной интерполяцией между точками (а это всего одно умножение на малобитное число). А при табличке из 512 точек с линейной интерполяцией достигается полная single точность в 24 бита.

[=SSN=]

worm2 в сообщении #1481224 писал(а):

Две таблицы — это как?

Поясню.

Для начала, нужно воздать должное студентам и таки вывести формулу Ньютона для вычисления обратного корня. :)

Допустим, мы хотим найти значение "x" при котором обращается в нуль функция f(x) заданная выражением:

$f(x)=a-\frac{1}{x^2}=0$

Очевидно, что эта функция равна нулю как раз при $x=\frac{1}{\sqrt{a}}$

Чтобы воспользоваться формулой Ньютона, вычислим производную функции f(x):

$f'(x)=\frac{2}{x^3}$

Откуда находим формулу для итераций:

$x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}=x_n-\frac{{x_n}^3}{2}(a-\frac{1}{{x_n}^2})=\frac{x_n}{2}\cdot(3-a\cdot{x_n}^2)$

Пусть вещественное положительное число "a" представлено в нормализованном формате с плавающей запятой:

$a=a_0,a_1a_2\dots a_{12}a_{13}\dots a_{24}\cdot2^N$, причем, нам известно, что $a_0=1$

Представим число "a" в виде:

$a=(1,a_1a_2\dots a_{12}+0,a_{13}\dots a_{24}\cdot2^{-12})\cdot2^N$

или:

$a=1,a_1a_2\dots a_{12}\cdot2^N+0,a_{13}\dots a_{24}\cdot 2^{-12}\cdot2^N=b+\Delta$

Причем:

$\frac{\Delta}{a}\leqslant\frac{\Delta}{b}< 2^{-12}$

Выберем начальное значение $x_0$ для итерационной формулы Ньютона в виде:

$x_0=\frac{1}{\sqrt{b}}=\frac{1}{\sqrt{a-\Delta}}$

Тогда для первой итерации находим значение $x_1$ в виде:

$x_1=\frac{x_0}{2} (3-a\cdot{x_0}^2)=\frac{1}{2\sqrt{a-\Delta}} (3-\frac{a}{a-\Delta})$

Или:

$x_1\approx\frac{1}{2\sqrt{a}}\cdot(1+\frac{1}{2}\cdot\frac{\Delta}{a}+\frac{3}{8}\cdot(\frac{\Delta}{a})^2)\cdot(2-\frac{\Delta}{a}-(\frac{\Delta}{a})^2)$

Или:

$x_1\approx\frac{1}{\sqrt{a}}\cdot(1-\frac{3}{8}(\frac{\Delta}{a})^2)$

Тогда для относительной ошибки $\delta x$ первой итерации находим оценку:

$\frac{x_1}{x}=\frac{x-\delta x}{x}=1-\frac{3}{8}(\frac{\Delta}{a})^2$

Или:

$\frac{\delta x}{x}=\frac{3}{8}(\frac{\Delta}{a})^2 < 2^{-25}$

То есть, для выбранной величины $\Delta$ одной итерации Ньютона достаточно для вычисления корня с одинарной точностью float..

...

Теперь снова вернемся к итерационной формуле Ньютона и воспользуемся двумя таблицами для вычисления начальных значений величин $\frac{3}{2}x_0$ и $\frac{1}{2}{x_0}^3$

Рассмотрим сначала случай четного "N".

Для вычисления $\frac{3}{2}x_0$ нам нужна таблица значений величин:

float f1[4096];

for (int n = 0; n < 4096; n++) f1 = 1.5f / sqrt(1.0f + (float)n/4096);

Для индексирования элементов в этой таблице мы будем использовать двоичные числа:

$a_{1 \dots {12}}=a_1a_2\dots a_{12}$

Само значение $\frac{3}{2}x_0$ при этом вычисляется по формуле:

$g_1=\frac{3}{2}x_0=f_1\cdot2^{-\frac{N}{2}}$

Для вычисления $\frac{1}{2}{x_0}^3$ нам нужна таблица значений величин:

float f2[4096];

for (int n = 0; n < 4096; n++) f2 = 0.5f / sqrt((1.0f + (float)n/4096)*(1.0f + (float)n/4096)*(1.0f + (float)n/4096));

Для индексирования элементов в этой таблице мы также будем использовать двоичные числа:

$a_{1 \dots {12}}=a_1a_2\dots a_{12}$

Само значение $\frac{1}{2}{x_0}^3$ при этом вычисляется по формуле:

$g_2=\frac{1}{2}{x_0}^3=f_2\cdot2^{-\frac{3N}{2}}$

Теперь формула для вычисления $x_1$ сводится к одному умножению и одному сложению:

$x_1=g_1+a\cdot g_2$

Еще одно умножение потребуется чтобы из $x_1$ получить $\sqrt{a}$:

$\sqrt{a}=x_1\cdot a$

...

Для случая нечетных "N" вычисления делают аналогично, с той лишь разницей, что нужно использовать индексы:

$a_{2 \dots {13}}=a_2a_3\dots a_{13}$

и вместо массивов f1 и f2 использовать массивы:

float f3[4096];

for (int n = 0; n < 4096; n++) f3 = 1.5f / sqrt(0.5f + (float)n/8192);

float f4[4096];

for (int n = 0; n < 4096; n++) f4 = 0.5f / sqrt((0.5f + (float)n/8192)*(0.5f + (float)n/8192)*(0.5f + (float)n/8192));

Ну и массивы f1 и f3 можно, КМК, объединить в один.
И точно также массивы f2 и f4 можно, КМК, тоже объединить в один.