Если экспонента 5, то число должно быть от 1*2^(5-15) до 1.99(9)*2^(5-15), а в случае 22 это не так.
Опытным путём удалось обнаружить, что при value == 22.0 :
exp (without bias) == 19
exp (with bias) == 5
mantissa (normalized) == 0.6875
Осталось понять как мантисса такой получилась.
поправка: mantissa == 0.6875
Если exp = 5, а мантисса 0.6875, то как раз будет
0.6875*(2^5) = 22
kipar
Да не, всё так, формула подходит как к флоат так и халф, но с одним маленьким ньюансом, хер поймёшь когда 0 когда 1 перед мантисой,
если взять флоат например то там мантиса равна 1.375, например 1.375*(2^4) = 22
можно найти добиться
Отредактируйте пожалуйста.Предпоследний абзац.
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 // Вернуться } }
Далее привожу аналогичную функцию на SSE-командах. По производительности уступает предыдущему варианту, но всё равно быстрее библиотечной logf.
_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 соответственно.
Тема в архиве.