Трюки с float: быстрое вычисление логарифмов. (комментарии) (3 стр) / Общее / Форум / Программирование игр / GameDev.ru

Kloun

все-таки интересно. ктонибудь на see логарифм нарисует?

Для одинарной точности использую следующую функцию натурального логарифма собственной разработки на основе аппроксимации по методу Чебышёва. Она даёт 7-8 верных значащих цифр, т.е. почти полную точность в рамках формата float. По производительности превосходит известные мне аналоги (включая библиотечную функцию logf), за исключением совсем грубых приближений. Для ускорения работы функции в ней не реализована поддержка денормализованных значений x; такие значения считаются нулём. Код предназначен для 32-разрядных проектов на Visual Studio, но может быть адаптирован и под другие среды, где разрешены ассемблерные вставки. Процессор должен поддерживать расширения набора инструкций AVX+FMA. В противном случае используемые SIMD-команды следует заменить SSE-аналогами.

#define salc _emit 0xD6                 // Определяем команду salc её машинным кодом
_declspec(naked) float _vectorcall ln1(float x)
  {
  static const float ct[6] =            // Таблица констант
    {
    1.0f,                               // 1
    0.576395561f,                       // 5^5*b2
    0.961802196f,                       // 5^3*b1
    2.88539007f,                        // 5*b0
    0.439077714f,                       // 5^7*b3
    0.693147181f                        // ln(2)
    };
  _asm
    {
    vmovd eax, xmm0                     // Прочитать в eax двоичное представление числа x
    cmp eax, 0x7F800000                 // Сравнить число x со значением Inf
    jge ln_end                          // Если x=+Inf или +NaN, вернуть аргумент x без изменений
    ror eax, 23                         // Циклически сдвинуть число так, чтобы порядок оказался в al
    movzx ecx, al                       // Получить в ecx порядок числа со смещением +127
    jecxz ln_break                      // Перейти, если x=0 или x денормализованное (т.е. близко к 0)
    salc                                // al=-1, если установлен старший бит явной части мантиссы; иначе al=0
    sub cl, al                          // ecx=k+127, где k - добавка к двоичному логарифму (целая часть)
    add al, 127                         // Изменить порядок x так, чтобы привести его к диапазону 0.75<=x<1.5
    sub ecx, 127                        // ecx=k
    ror eax, 9                          // Получить в eax двоичное представление числа и проверить его знак
    vcvtsi2ss xmm3, xmm3, ecx           // xmm3=k, это добавка к двоичному логарифму (целая часть)
    sbb ecx, ecx                        // ecx=-1 (значение NaN), если x<0; или ecx=0, если x>0
    jnc ln_cont                         // Перейти, если x>0, т.е. аргумент корректный. Ошибка, если x<=0
      ln_break :                        // Точка входа для случая x=0 или денормализованного x
    xchg eax, ecx                       // Далее результат будет формироваться в eax
    or eax, 0xFF800000                  // Сформировать в eax значение -Inf, если x=0 или x денормализованное
      ln_cont :                         // Точка входа для обработки случая x>0
    vmovd xmm0, eax                     // Перенести из eax в xmm0 результат или приведённое значение x
    jnz ln_end                          // Перейти в конец функции, если ошибка, и продолжить для x>0
    mov eax, offset ct                  // В eax адрес таблицы констант
    vaddss xmm1, xmm0, [eax]            // xmm1=x+1
    vsubss xmm0, xmm0, [eax]            // xmm0=x-1
    vdivss xmm0, xmm0, xmm1             // xmm0=(x-1)/(x+1)=t/5
    vmovss xmm1, [eax+16]               // Инициализировать сумму старшим коэффициентом: xmm1=5^7*b3
    vmulss xmm2, xmm0, xmm0             // xmm2=t^2/25 - подготовить аргумент для вычисления многочлена
      ln_loop :                         // Цикл вычисления многочлена (ecx пробегает от 1 до 3)
    inc ecx                             // Увеличить счётчик цикла на 1. Флаг pf установится только при ecx=3
    vfmadd213ss xmm1, xmm2, [eax+ecx*4] // Выполнить очередной шаг вычислений по схеме Горнера
    jnp ln_loop                         // Учесть все 4 коэффициента. В xmm1 формируется 5*f(t)
    vfmadd213ss xmm0, xmm1, xmm3        // xmm0=t*f(t)+k - получили готовый двоичный логарифм
    vmulss xmm0, xmm0, [eax+20]         // Перевести двоичный логарифм в натуральный
      ln_end :                          // Результат в xmm0 готов
    ret                                 // Вернуться
    }
  }

_declspec(naked) float _vectorcall ln1(float x)
  {
  static const float ct[6] =            // Таблица констант
    {
    1.0f,                               // 1
    0.576395561f,                       // 5^5*b2
    0.961802196f,                       // 5^3*b1
    2.88539007f,                        // 5*b0
    0.439077714f,                       // 5^7*b3
    0.693147181f                        // ln(2)
    };
  _asm
    {
    movd eax, xmm0                      // Прочитать в eax двоичное представление числа x
    cmp eax, 0x7F800000                 // Сравнить число x со значением Inf
    jge ln_end                          // Если x=+Inf или +NaN, вернуть аргумент x без изменений
    ror eax, 23                         // Циклически сдвинуть число так, чтобы порядок оказался в al
    movzx ecx, al                       // Получить в ecx порядок числа со смещением +127
    jecxz ln_break                      // Перейти, если x=0 или x денормализованное (т.е. близко к 0)
    salc                                // al=-1, если установлен старший бит явной части мантиссы; иначе al=0
    sub cl, al                          // ecx=k+127, где k - добавка к двоичному логарифму
    add al, 127                         // Изменить порядок x так, чтобы привести его к диапазону 0.75<=x<1.5
    sub ecx, 127                        // ecx=k
    ror eax, 9                          // Получить в eax двоичное представление числа и проверить его знак
    cvtsi2ss xmm3, ecx                  // xmm3=k, это добавка к двоичному логарифму (целая часть)
    sbb ecx, ecx                        // ecx=-1 (значение NaN), если x<0; или ecx=0, если x>0
    jnc ln_cont                         // Перейти, если x>0, т.е. аргумент корректный. Ошибка, если x<=0
      ln_break :                        // Точка входа для случая x=0 или денормализованного x
    xchg eax, ecx                       // Далее результат будет формироваться в eax
    or eax, 0xFF800000                  // Сформировать в eax значение -Inf, если x=0 или x денормализованное
      ln_cont :                         // Точка входа для обработки случая x>0
    movd xmm0, eax                      // Перенести из eax в xmm0 результат или приведённое значение x
    jnz ln_end                          // Перейти в конец функции, если ошибка, и продолжить для x>0
    mov eax, offset ct                  // В eax адрес таблицы констант
    movaps xmm1, xmm0                   // xmm1=x (приведённое значение)
    subss xmm0, [eax]                   // xmm0=x-1
    addss xmm1, [eax]                   // xmm1=x+1
    divss xmm0, xmm1                    // xmm0=(x-1)/(x+1)=t/5
    movss xmm1, [eax+16]                // Инициализировать сумму старшим коэффициентом: xmm1=5^7*b3
    movaps xmm2, xmm0                   // xmm2=t/5
    mulss xmm2, xmm2                    // xmm2=t^2/25 - подготовить аргумент для вычисления многочлена
      ln_loop :                         // Цикл вычисления многочлена по схеме Горнера (ecx пробегает от 1 до 3)
    inc ecx                             // Увеличить счётчик цикла на 1. Флаг pf установится только при ecx=3
    mulss xmm1, xmm2                    // Умножить текущую сумму на t^2/25
    addss xmm1, [eax+ecx*4]             // Прибавить текущий коэффициент
    jnp ln_loop                         // Учесть все 4 коэффициента. В xmm1 формируется 5*f(t)
    mulss xmm0, xmm1                    // xmm0=t*f(t) - дробная часть двоичного логарифма
    addss xmm0, xmm3                    // xmm0=t*f(t)+k - готовый двоичный логарифм
    mulss xmm0, [eax+20]                // Перевести двоичный логарифм в натуральный
      ln_end :                          // Результат в xmm0 готов
    ret                                 // Вернуться
    }
  }

Aenigma
Можешь сравнить с моим?
https://gamedev.ru/code/forum/?id=233033&page=8&m=4682379#m114
По точности и скорости.
Асм получается такой:
https://godbolt.org/z/Ee91EE6aj

FordPerfect, у вас побольше констант, и в 32-битной версии ресурсы тратятся впустую на пересылку данных через машинный стек и в стек сопроцессора. 64-битный код должен быть быстрее. Ну и ручная оптимизация на ассемблере иногда позволяет ускорить код (я, кстати, сознательно оставил цикл, потому что он мне нравится). Так что приведённые результаты тестирования не окончательные. В общем, собрал ваш код под 32-бита на Visual Studio 2013 с максимальной оптимизацией по скорости. У себя отключил умножение на ln(2), чтобы был двоичный логарифм, как у вас. Тестировал на AMD FX-8350. Получил следующие результаты.

Среднее время работы
Ваш код: 1513
Мой AVX-код: 952
Мой SSE-код: 1310

Среднеквадратическая относительная погрешность на промежутке [0.5; 2)
Ваш код: 0.001461
Мой AVX-код: 5.048E-8
Мой SSE-код: 5.253E-8

У вас, похоже, как-то не так константы аппроксимации определяются. Каким методом пользовались? Ещё вопрос: что делаете, если на вход поступает денормализованное число? Самому разбираться сейчас нет времени. Я в этом случае считаю такое число нулём.
И ещё один момент: должно быть lb(0)=-Inf, а не NaN.

Aenigma
Спасибо за тест.

> впустую на пересылку данных через машинный стек и в стек сопроцессора
Вроде это ABI такое, и чинится -msseregparm.

> Ваш код: 0.001461
Странно, у меня получается заметно другое (~1e-7):
https://rextester.com/MRRNP67539 (нажать F8 чтобы показало график)
https://rextester.com/VIUG71532

> Каким методом пользовались?
Код старый, это ещё вспомнить. Но судя по графику на https://rextester.com/MRRNP67539 - там что-то совсем тупое, типа Тейлора.

> Ещё вопрос: что делаете, если на вход поступает денормализованное число?
Делаю вид, что оно 1.M...M*2^-127, а не 0.M...M*2^-126.

> И ещё один момент: должно быть lb(0)=-Inf, а не NaN.
Точно.

> я, кстати, сознательно оставил цикл, потому что он мне нравится
Гм, и по производительности не бьёт?

Ну и MSVC ассемблер с FMA выглядит чуть странно, с учётом того, что FMA типично доступно только на 64-битных машинах (хотя там - и из 32-битного кода тоже доступно), а под x64 MSVC ассемблер вроде не умеет.

Aenigma
> Среднеквадратическая относительная погрешность на промежутке
А, понял. Я абсолютную замерял. Да, относительная у меня отстой.

Тейлора вообще не надо использовать в компьютерной арифметике, от него мало толку. Нужно использовать Чебышёва или Паде-Чебышёва (смотреть по ситуации, что лучше). В идеале - равномерное приближение, но его коэффициенты сложно определяются, а по точности почти то же, что и для Чебышёва.

Да, цикл бьёт по производительности, но не критично, а на новом Ryzen 5600X - вообще почти незаметно. Он или круто предсказывает переходы, или умеет разворачивать небольшие циклы (что в сущности одно и то же).

Самое главное, вам нужно выбрать хорошую аппроксимацию. Судя по всему, у вас аппроксимация строится на промежутке [1; 2), а это очень плохо для точности. Вот, представьте, что x=0.999, логарифм которого примерно равен -0.001443 (до 4-х значащих цифр). После изменения порядка для приведения в диапазон [1; 2) число превратится в 1.998, а его логарифм равен 0.9986 (тоже до 4-х значащих цифр). После вычитания 1 получаем -0.0014. Две значащие цифры потеряли. Чувствуете, в чём источник проблем? Промежуток аппроксимации должен выбираться вокруг 1.

Я после глубокого погружения в задачу пришёл к выводу, что для логарифма с одинарной точностью лучше всего использовать аппроксимацию многочленами Чебышёва на промежутке x из [2/3; 3/2]. Этот промежуток несколько избыточен, зато правая граница у него удобная. Фактически задействуется более узкий промежуток [3/4; 3/2] - его достаточно, т.к. при этом правая граница в 2 раза больше левой. Многочленами Чебышёва аппроксимируется функция

        1      5+t                 x-1
f(t) = --- lb ------ ,  где t = 5 ----- ;  -1<=t<=1.
        t      5-t                 x+1

При таком подходе для достижения одинарной точности достаточно взять 4 члена разложения в ряд Чебышёва.

Aenigma
> Он или круто предсказывает переходы, или умеет разворачивать небольшие циклы (что в сущности одно и то же).
По идее, т. к. длина цикла почти всегда одна и та же, относительно новый предсказатель (который помнит историю) обработает его с примерно 100% точностью.

О, да там условие цикла по parity flag.
Не знаю, насколько оно оправдано, но занятно!

> Вот, представьте, что x=0.999, логарифм которого примерно равен -0.001443 (до 4-х значащих цифр). После изменения порядка для приведения в диапазон [1; 2) число превратится в 1.998, а его логарифм равен 0.9986 (тоже до 4-х значащих цифр). После вычитания 1 получаем -0.0014.
Есть такое. Если тестировать на [1;2) - относительная погрешность кажется намного лучше, чем на самом деле.
Да, конечно, переделать стоит.

Код GCC для моей функции выглядит приятнее, чем MSVC:
https://godbolt.org/z/f6efqjY4Y

Насчёт предсказателя переходов познавательно, спасибо.

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

Нафига деление, что в одном варианте, что в другом? Не будет ли достаточно одного полинома?

Aenigma
> В идеале - равномерное приближение, но его коэффициенты сложно определяются, а по точности почти то же, что и для Чебышёва.
Лично у меня наоборот: алгоритм подбора коэффициентов по \(L_\infty\) давно написан, а Чебышева я никогда не заморачивался реализовывать.

> Самое главное, вам нужно выбрать хорошую аппроксимацию.
Скорее, нужно перестать пользоваться однобоким определением точности, когда рассматривается только погрешность по y, а погрешность по x почему-то считается равной нулю. Для практических целей это просто бесполезная потеря производительности на ровном месте. Есть задачи, где, действительно, нужна точность в районе единицы, но в них надо использовать не log, а log1p.

}:+()___ [Smile]
> Нафига деление, что в одном варианте, что в другом? Не будет ли достаточно одного полинома?
У меня вроде не получалось нормальную точность с полиномом.

> когда рассматривается только погрешность по y, а погрешность по x почему-то считается равной нулю
Кэхэн и компания как раз ратуют за разделение этих ошибок.

Кстати, в виде оффтопика, приведу кэхэновский пример о том, как внесение ошибок может (систематически!) повышать точность.

Итак, есть функция

double f(double x)
{
    if(x==0.0) return 1.0;
    else return (exp(x)-1.0)/x;
}

Ну и у неё всё численно плохо в районе 0.0.
А теперь, давайте внесём в неё ещё погрешностей:

double f(double x)
{
    double y=exp(x);
    double z=log(y);
    if(z==0.0) return 1.0;
    else return (y-1.0)/z;
}

И вот у неё точность порядка пары ULP примерно на всём диапазоне.
Пример, естественно, не про то что так надо делать (существует expm1). А про то, что гипотеза "внесение погрешности не повышает точность" - бывает ошибочна. Эта гипотеза может сбить с толку при дебаге, например.
Точность в данном примере закладывается на то, что exp и log реализованы достаточно точно возле 0.0 и 1.0 соответственно.