Войти
ПрограммированиеФорумОбщее

Трюки с float: быстрое вычисление логарифмов. (комментарии) (5 стр)

Страницы: 14 5 6 710 Следующая »
#60
8:25, 10 янв. 2022

Заметил, что в <math.h> есть expm1, но нету cosm1 (даже в C23). В некоторых других библиотеках - есть.

Aenigma
> К сожалению, на VS компилятор генерирует не лучший код, поэтому не могу адекватно замерить производительность и сравнить со своим кодом.
По идее, код из https://gamedev.ru/code/forum/?id=132465&page=3&m=5503341#m35 не особо сложно перевести на GCC ассемблер, и сравнить хоть бы и прямо на godbolt (он умеет запускать код).


#61
8:29, 10 янв. 2022

Тест для #56:
https://rextester.com/MWKR48397

#62
(Правка: 10:51) 10:48, 10 янв. 2022

Авось не сломал:
https://godbolt.org/z/zEE9hxTza

Это перевод #35 на GCC асм.

#63
(Правка: 12:51) 12:39, 10 янв. 2022

FordPerfect, спасибо. У вас среднеквадратическая относительная погрешность на 1% ниже, чем у меня. Это вполне хорошее соответствие, учитывая, что мы по разным критериям определяли константы (я - на основе коэффициентов ряда Чебышёва). А быстродействие от запуска к запуску показывает большой разброс, вряд ли можно доверять полученным цифрам. Кстати, если у меня развернуть цикл и оптимизировать передачу параметра/результата (в примере это делается через стек), можно ещё несколько процентов быстродействия выгадать. В целом, можно констатировать, что мы очень близки к оптимуму.

P.S. Почему-то в моём коде есть команда vzeroupper, а в вашем нет. По идее, она и у меня не нужна, т.к. я работаю только с xmm (не трогаю старшие половины ymm).

#64
19:16, 10 янв. 2022

Что-то я дизасм показывал под 32 бита, а запускал под 64.
Вот оба под x64:
https://godbolt.org/z/443eYdoE9
Заодно от использования FPU стека избавился.

Aenigma
> А быстродействие от запуска к запуску показывает большой разброс, вряд ли можно доверять полученным цифрам.
Да, на онлайн-компиляторах это типично так.

> P.S. Почему-то в моём коде есть команда vzeroupper, а в вашем нет. По идее, она и у меня не нужна, т.к. я работаю только с xmm (не трогаю старшие половины ymm).
Подозреваю, GCC толком не знает, что внутри ассемблерной вставки, и перестраховался.

Я, кстати, статью писал, про ассемблерные вставки в GCC, может пригодится: https://gamedev.ru/code/articles/gcc_inline_asm .

#65
19:34, 10 янв. 2022

Там, кстати, баг был. xmm3 используется в коде но не указан в clobber-list. Исправил:
https://godbolt.org/z/KovK5sne9
Кстати, только GCC понимает "Yz", а clang - нет. Так что для портирования с этим надо бы будет что-то сделать.

#66
20:02, 10 янв. 2022

> понимает "Yz"
Переделал на "x":
https://godbolt.org/z/x5a4o97Gf
Вроде без проблем. Теперь clang тоже понимает.

#67
(Правка: 11:18) 6:43, 11 янв. 2022

FordPerfect, спасибо за информацию, может пригодится.

Для полноты картины выкладываю ещё несколько вариантов кода для логарифма с двойной точностью. Хотя это и не совсем по теме, но, думаю, представляет интерес. Производительность в зависимости от процессора может до 2-х раз превышать библиотечную функцию. По точности незначительно уступают ей.

1. Натуральный логарифм с двойной точностью на AVX-командах с использованием аппроксимации по методу Чебышёва. Имеет максимальное быстродействие из предложенных вариантов, но точность чуть пониже, чем по методу Паде-Чебышёва.

double _declspec(naked) _vectorcall ln21(double x)
  {
  static const double ct[11] = // Таблица констант
    {                            
    -1.0,                      // -1
    1.41421356237309505,       // 2^0.5
    0.262334391706879229,      // b5
    0.320598534127666242,      // b4
    0.412198585848994984,      // b3
    0.577078016345480558,      // b2
    0.961796693925989802,      // b1
    2.88539008177792681,       // b0
    0.213668367342346557,      // b7
    0.220912115142181838,      // b6
    0.693147180559945309       // ln 2
    };
  _asm
    {
    mov cx,0x3FF0              // Выставить в ch страший байт порядка и подготовить маску в cl
    vpextrw eax,xmm0,3         // Прочитать в ax старшее слово числа, содержащее знаковый бит и порядок
    mov edx,offset ct          // В edx - адрес таблицы констант
    cmp ax,0x7FF0              // Проверить условие, когда на вход подано +Inf или +NaN
    jge ln_end                 // В этом случае функция просто возвращает исходное значение
    or cl,al                   // Выставить младший байт порядка x
    sar ax,4                   // Получить в eax порядок числа плюс смещение 1023
    jg ln_cont                 // Перейти, если x>0 и число не является денормализованным
    vxorpd xmm0,xmm0,xmm0      // xmm0 = 0
    test ax,0x07FF             // zf=1, если x=0 или денормализованное
    mov cx,0xFFF0              // Подготовить в cx старшее слово значения -Inf
    jz ln_cont                 // Перейти, если x=0 или x денормализованное
    or cl,ch                   // Теперь cx=FFFFh - старшее слово значения NaN; sf=1
      ln_cont:                 // Точка входа для формирования приведённого значения x или результата
    vpinsrw xmm0,xmm0,ecx,3    // Заменить порядок числа. В штатной ветке 1<=x<2
    jle ln_end                 // Перейти в конец функции, если результат -Inf или NaN
    sub eax,1023               // Вычислить в eax правильный порядок, это будет целая часть логарифма
    vcomisd xmm0,[edx+8]       // Сравнить x и корень из 2
    jb ln_x_done               // Перейти, если 1 <= x < 2^0.5, т.е. не требуется коррекция порядка
    and ecx,-17                // Вычесть из порядка 1, чтобы стало 2^0.5/2 <= x < 1
    inc eax                    // Целая часть логарифма соотвтетственно увеличивается на 1
    vpinsrw xmm0,xmm0,ecx,3    // Обновляем значение xmm0
      ln_x_done:               // В xmm0 - приведённое значение x, такое что 2^0.5/2 <= x < 2^0.5
    vaddsd xmm4,xmm0,[edx]     // xmm4 = x-1
    vsubsd xmm2,xmm0,[edx]     // xmm2 = x+1
    vcvtsi2sd xmm3,xmm3,eax    // В xmm3 - целая часть двоичного логарифма
    vdivsd xmm0,xmm4,xmm2      // xmm0 = (x-1)/(x+1) = u
    vmovupd xmm1,[edx+64]      // Инициализировать суммы старшими коэффициентами: xmm1 = b6 : b7
    vmulsd xmm2,xmm0,xmm0      // xmm2 = u^2
    xor eax,eax                // eax = 0 - инициализируем счётчик цикла
    vmulsd xmm4,xmm2,xmm2      // xmm4 = u^4
    vunpcklpd xmm4,xmm4,xmm4   // xmm4 = u^4 : u^4
      ln_loop:                 // Цикл вычисления многочлена по схемам Эстрина и Горнера
    add al,16                  // Получить в eax смещение очередной пары коэффициентов относительно ct
    vfmadd213pd xmm1,xmm4,[edx+eax] // Умножить текущие суммы на u^4 и прибавить коэффициенты
    jnp ln_loop                // Повторить цикл 3 раза (учесть все 4 пары коэффициентов)
    vunpckhpd xmm4,xmm1,xmm1   // В xmm1 - сумма по нечётным членам, в xmm4 - по чётным
    vfmadd213sd xmm1,xmm2,xmm4 // В xmm1 - готовая частичная сумма ряда Чебышёва
    vfmadd213sd xmm0,xmm1,xmm3 // В xmm0 - готовый двоичный логарифм
    vmulsd xmm0,xmm0,[edx+80]  // Перевести двоичный логарифм в натуральный
      ln_end:                  // Результат в xmm0 готов
    ret                        // Вернуться
    }
  }
#68
(Правка: 11:50) 6:43, 11 янв. 2022

2. Натуральный логарифм с двойной точностью на AVX-командах с использованием аппроксимации по методу Паде-Чебышёва. Имеет максимальную точность из предложенных вариантов, но быстродействие чуть пониже, чем по методу Чебышёва.

double _declspec(naked) _vectorcall ln22(double x)
  {
  static const double ct[11] = // Таблица констант
    {                           
    -1.0,                      // -1
    1.41421356237309505,       // 2^0.5
    1.00565869936870511,       // a2
    0.644659714443160503,      // b2
    -3.23215475233190253,      // a1
    -1.40993519139999135,      // b1
    2.50817479464989201,       // a0
    0.869267143631546224,      // b0
    -0.0438956451157205899,    // a3
    -0.0722936249724093899,    // b3
    0.693147180559945309       // ln 2
    };
  _asm
    {
    mov cx,0x3FF0              // Выставить в ch страший байт порядка и подготовить маску в cl
    vpextrw eax,xmm0,3         // Прочитать в ax старшее слово числа, содержащее знаковый бит и порядок
    mov edx,offset ct          // В edx - адрес таблицы констант
    cmp ax,0x7FF0              // Проверить условие, когда на вход подано +Inf или +NaN
    jge ln_end                 // В этом случае функция просто возвращает исходное значение
    or cl,al                   // Выставить младший байт порядка x
    sar ax,4                   // Получить в eax порядок числа плюс смещение 1023
    jg ln_cont                 // Перейти, если x>0 и число не является денормализованным
    vxorpd xmm0,xmm0,xmm0      // xmm0 = 0
    test ax,0x07FF             // zf=1, если x=0 или денормализованное
    mov cx,0xFFF0              // Подготовить в cx старшее слово значения -Inf
    jz ln_cont                 // Перейти, если x=0 или x денормализованное
    or cl,ch                   // Теперь cx=FFFFh - старшее слово значения NaN; sf=1
      ln_cont:                 // Точка входа для формирования приведённого значения x или результата
    vpinsrw xmm0,xmm0,ecx,3    // Заменить порядок числа. В штатной ветке 1<=x<2
    jle ln_end                 // Перейти в конец функции, если результат -Inf или NaN
    sub eax,1023               // Вычислить в eax правильный порядок, это будет целая часть логарифма
    vcomisd xmm0,[edx+8]       // Сравнить x и корень из 2
    jb ln_x_done               // Перейти, если 1 <= x < 2^0.5, т.е. не требуется коррекция порядка
    and ecx,-17                // Вычесть из порядка 1, чтобы стало 2^0.5/2 <= x < 1
    inc eax                    // Целая часть логарифма соотвтетственно увеличивается на 1
    vpinsrw xmm0,xmm0,ecx,3    // Обновляем значение xmm0
      ln_x_done:               // В xmm0 - приведённое значение x, такое что 2^0.5/2 <= x < 2^0.5
    vsubsd xmm2,xmm0,[edx]     // xmm2 = x+1
    vaddsd xmm4,xmm0,[edx]     // xmm4 = x-1
    vdivsd xmm2,xmm4,xmm2      // xmm2 = (x-1)/(x+1) = u
    vcvtsi2sd xmm3,xmm3,eax    // В xmm3 - целая часть двоичного логарифма
    xor eax,eax                // eax = 0 - инициализируем счётчик цикла
    vmulsd xmm2,xmm2,xmm2      // xmm2 = u^2
    vmovupd xmm1,[edx+64]      // Инициализировать суммы старшими коэффициентами: xmm1 = b3 : a3
    vunpcklpd xmm2,xmm2,xmm2   // xmm2 = u^2 : u^2
      ln_loop:                 // Цикл вычисления многочленов по схеме Горнера
    add al,16                  // Получить в eax смещение очередной пары коэффициентов относительно ct
    vfmadd213pd xmm1,xmm2,[edx+eax] // Умножить текущие суммы на u^2 и прибавить коэффициенты
    jnp ln_loop                // Повторить цикл 3 раза (учесть все 4 пары коэффициентов)
    vunpckhpd xmm2,xmm1,xmm1   // xmm1 = P; xmm2 = Q
    vfmadd213sd xmm0,xmm2,xmm2 // xmm0 = Q*(x+1)
    vdivsd xmm0,xmm1,xmm0      // xmm0 = P/Q/(x+1)
    vfmadd213sd xmm0,xmm4,xmm3 // xmm0 = u*P/Q+k - готовый двоичный логарифм
    vmulsd xmm0,xmm0,[edx+80]  // Перевести двоичный логарифм в натуральный
      ln_end:                  // Результат в xmm0 готов
    ret                        // Вернуться
    }
  }

3. Натуральный логарифм с двойной точностью на SSE-командах с использованием аппроксимации по методу Чебышёва.

double _declspec(naked) _vectorcall ln23(double x)
  {
  _declspec(align(16)) static const double ct[11] = // Таблица констант
    {                            
    -1.0,                    // -1
    1.41421356237309505,     // 2^0.5
    0.262334391706879229,    // b5
    0.320598534127666242,    // b4
    0.412198585848994984,    // b3
    0.577078016345480558,    // b2
    0.961796693925989802,    // b1
    2.88539008177792681,     // b0
    0.213668367342346557,    // b7
    0.220912115142181838,    // b6
    0.693147180559945309     // ln 2
    };
  _asm
    {
    mov cx,0x3FF0            // Выставить в ch страший байт порядка, и подготовить маску в cl
    pextrw eax,xmm0,3        // Прочитать в ax старшее слово числа, содержащее знаковый бит и порядок
    cmp ax,0x7FF0            // Проверить условие, когда на вход подано +Inf или +NaN
    jge ln_end               // В этом случае функция просто возвращает исходное значение
    or cl,al                 // Выставить младший байт порядка x
    sar ax,4                 // Получить в eax порядок числа плюс смещение 1023
    jg ln_cont               // Перейти, если x>0 и число не является денормализованным
    xorpd xmm0,xmm0          // xmm0 = 0
    test ax,0x07FF           // zf=1, если x=0 или денормализованное
    mov cx,0xFFF0            // Подготовить в cx старшее слово значения -Inf
    jz ln_cont               // Перейти, если x=0 или x денормализованное
    or cl,ch                 // Теперь cx=FFFFh - старшее слово значения NaN; sf=1
      ln_cont:               // Точка входа для формирования приведённого значения x или результата
    pinsrw xmm0,ecx,3        // Заменить порядок числа. В штатной ветке 1<=x<2
    jle ln_end               // Перейти в конец функции, если результат -Inf или NaN
    mov edx,offset ct        // В edx - адрес таблицы констант
    sub eax,1023             // Вычислить в eax правильный порядок, это будет целая часть логарифма
    comisd xmm0,[edx+8]      // Сравнить x и корень из 2
    jb ln_x_done             // Перейти, если 1 <= x < 2^0.5, т.е. не требуется коррекция порядка
    and ecx,-17              // Вычесть из порядка 1, чтобы стало 2^0.5/2 <= x < 1
    inc eax                  // Целая часть логарифма соотвтетственно увеличивается на 1
    pinsrw xmm0,ecx,3        // Обновляем значение xmm0
      ln_x_done:             // В xmm0 - приведённое значение x, такое что 2^0.5/2 <= x < 2^0.5
    movapd xmm1,xmm0         // xmm1 = xmm0 = x (приведённое значение)
    subsd xmm0,[edx]         // xmm0 = x+1
    addsd xmm1,[edx]         // xmm1 = x-1
    cvtsi2sd xmm3,eax        // В xmm3 - целая часть двоичного логарифма
    xor eax,eax              // eax = 0 - инициализируем счётчик цикла
    divsd xmm1,xmm0          // xmm1 = (x-1)/(x+1) = u
    movapd xmm0,[edx+64]     // Инициализировать суммы старшими коэффициентами: xmm0 = b6 : b7
    movapd xmm2,xmm1         // xmm2 = u
    mulsd xmm2,xmm2          // xmm2 = u^2
    movddup xmm4,xmm2        // xmm4 = u^2 : u^2
    mulpd xmm4,xmm4          // xmm4 = u^4 : u^4
      ln_loop:               // Цикл вычисления многочлена по схемам Эстрина и Горнера
    add al,16                // Получить в eax смещение очередной пары коэффициентов относительно ct
    mulpd xmm0,xmm4          // Умножить текущие суммы на u^4
    addpd xmm0,[edx+eax]     // Прибавить текущие коэффициенты
    jnp ln_loop              // Повторить цикл 3 раза (учесть все 4 пары коэффициентов)
    mulsd xmm2,xmm0          // Сумму по нечётным членам умножить на u^2 и поместить в xmm2
    unpckhpd xmm0,xmm0       // В xmm0 - сумма по чётным членам
    addsd xmm0,xmm2          // В xmm0 - готовая частичная сумма ряда Чебышёва
    mulsd xmm0,xmm1          // Получить в xmm0 дробную часть двоичного логарифма
    addsd xmm0,xmm3          // Получить в xmm0 готовый двоичный логарифм
    mulsd xmm0,[edx+80]      // Перевести двоичный логарифм в натуральный
      ln_end:                // Результат в xmm0 готов
    ret                      // Вернуться
    }
  }
#69
11:15, 11 янв. 2022

Попробовал взять коэффициенты из https://sleef.org/ :
https://github.com/shibatch/sleef/blob/master/src/libm/sleefsimdsp.c#L2746
И они несколько получше моих.

#70
(Правка: 12:14) 12:12, 11 янв. 2022

FordPerfect
Ну, понятие "лучше" имеет изрядную долю неопределённости. Оно зависит от выбранной метрики. Можно целенаправленно оптимизировать среднеквадратическую относительную погрешность и получать идеальные коэффициенты с этой точки зрения.

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

#71
13:56, 11 янв. 2022

Aenigma
> Во-первых, почти любая математическая операция над наилучшим равномерным приближением сразу же лишает его этого статуса.
Важен не статус, а возможность делать анализ погрешностей: зная максимальную погрешность до операции можно оценить максимальную погрешность после.

> Во-вторых, для конечного пользователя вполне разумным является, например, требование не минимаксной ошибки, а минимальной средней или среднеквадратической ошибки.
А вот в этом случае ничего оценить нельзя: расчет может залипнуть в особой точке, в которой погрешность в 100 раз больше средней и вся точность улетит в трубу. В общем, все эти средние погрешности применимы в достаточно редких случаях (ибо понятие "среднего" зависит от выбора распределения, в отличие от "максимума").

#72
15:10, 11 янв. 2022

}:+()___ [Smile]

Важен не статус, а возможность делать анализ погрешностей

Когда строится универсальная аппроксимация функции, которая будет применяться в самых разных областях знаний, то важна как раз уверенность, что аппроксимация хорошая, т.е. обладает некоторым свойством оптимальности, можно сказать, "статусом".

расчет может залипнуть в особой точке, в которой погрешность в 100 раз больше средней и вся точность улетит в трубу

Люди, которые занимаются аппроксимациями, обладают достаточным опытом, чтобы понимать, как должна вести себя аппроксимируемая функция и как лучше выбрать промежуток аппроксимации. Именно поэтому, например, вместо логарифма строится аппроксимация функции f(t)=lb((a+t)/(a-t))/t и т.д.

#73
(Правка: 19:55) 19:55, 11 янв. 2022

Aenigma
По тому параметру, который меня больше всего интересует (worst-case погрешность в ULP) #35, #56 и https://github.com/shibatch/sleef/blob/master/src/libm/sleefsimdsp.c#L2746 одинаковые (3 ULP).
Погрешность в ULP для log2 в каком-то смысле абсолютная, т. к. на [1;2) ULP константный.

Но, походу, я криво мерял.
Перемерил. Тупо подставив константы из sleef в свой код. Для [0.5; 2.0):

#56:
  Linf: abs=1.19e-07, rel=2.47e-07, 3.000 ulp
  RMS : abs=2.28e-08, rel=5.21e-08, 0.632 ulp
sleef:
  Linf: abs=1.19e-07, rel=2.47e-07, 3.000 ulp
  RMS : abs=2.29e-08, rel=5.22e-08, 0.634 ulp

https://rextester.com/IGOS68851

#74
5:33, 12 янв. 2022

Чуть причесал код. Получилось так:
https://rextester.com/XDLJ1747

float log2_v4(float x)
{
    uint32_t bits,bias=0x00400000u;
    memcpy(&bits,&x,sizeof(float));
    if((bits>>23)-1u>253u)
    {
        if(!(x>0.0f)) return (x==0.0f?-INFINITY:nanf(""));
        if(x==+INFINITY) return +INFINITY;
        return log2_v4(4294967296.0f*x)-32.0f;
    }
    bits+=bias;
    float y=float(int32_t(bits>>23))-127.0f;
    bits=((bits&0x007FFFFFu)|0x3F800000u)-bias;
    memcpy(&x,&bits,sizeof(float));
    x=(x-1.0f)/(x+1.0f);
    float x2=x*x;
    return y+x*(2.88539003e+00f
           +x2*(9.61808794e-01f
           +x2*(5.76110615e-01f
           +x2*(4.42826415e-01f))));
}

Так норм?

Страницы: 14 5 6 710 Следующая »
ПрограммированиеФорумОбщее