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

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

Страницы: 16 7 8 9 10 11
#150
(Правка: 1:36) 1:34, 27 янв 2022

Подобрал коэффициенты с помощью https://www.sollya.org/ :

float log2_v10(float x)
{
    float y=-127.0f;
    uint32_t bits,bias=0x00400000u;
    memcpy(&bits,&x,sizeof(float));
    if((bits>>23)-1u>253u)
    {
        if(!(x>0.0f)) y=(x==0.0f?-INFINITY:nanf(""));
        else if(x==+INFINITY) y=+INFINITY;
        else {x*=4294967296.0f;memcpy(&bits,&x,sizeof(float));y-=32.0f;}
    }
    bits+=bias;
    y+=float(int32_t(bits>>23));
    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.8853900432586669921875f
           +x2*(0.961808860301971435546875f
           +x2*(0.57604229450225830078125f
           +x2*(0.44434440135955810546875f))));
}

https://rextester.com/AVFHUR5660

  Linf: abs=1.09e-07, rel=2.44e-07, 3.397 ulp
  RMS : abs= 2.1e-08, rel=4.71e-08, 0.580 ulp
#151
(Правка: 7:46) 7:37, 27 янв 2022

FordPerfect, напомните, пожалуйста, значение RMS rel = 4.71e-08 получено для какого диапазона x?
Я для диапазона 0.5<=x<2 (по всем числам диапазона, кроме x=1) получил данную характеристику, равную 4.52359719E-8. Оптимизация проводилась для инструкций FMA, аппроксимация строится на промежутке 0.75<=x<1.5. Коэффициенты следующие:

    2.88539004f,                        // b0
    0.961808383f,                       // b1
    0.576109707f,                       // b2
    0.442826569f,                       // b3
#152
8:05, 27 янв 2022

Aenigma
> значение RMS rel = 4.71e-08 получено для какого диапазона x?
[0.5; 2.0).
Это без FMA. С FMA получается 4.52e-08 (точнее 4.52388242e-08).

#153
(Правка: 8:37) 8:22, 27 янв 2022

Сравнил графики теоретической относительной погрешности для этих двух наборов коэффициентов: они немного разные, но практически идентичные по качеству. Главные максимумы находятся в тех же местах и имеют такую же величину. RMS тоже почти совпала. Кстати, я для повышения точности RMS измерял, используя команду fyl2x:

double dlb(float x, float y)
  {
  _asm{
  fld1
  fld x
  fyl2x
  fld y
  fsub st(0),st(1)
  fdivrp st(1),st(0)
  fmul st(0),st(0)
  }}
#154
3:56, 28 янв 2022

Подставил коэффициенты из #151 в свою замерялку, получил:
Без FMA

  Linf: abs=1.08e-07, rel=2.46e-07, 3.451 ulp
  RMS : abs= 2.1e-08, rel=4.71e-08, 0.580 ulp

С FMA

  Linf: abs=1.03e-07, rel= 2.1e-07, 2.991 ulp
  RMS : abs=1.94e-08, rel=4.52e-08, 0.557 ulp

Да, результаты очень близкие, хотя коэффициенты разные:

       #151      bits                          #150      bits 
2.88539004  [4038AA3B]   2.8853900432586669921875   [4038AA3B]
0.961808383 [3F763913]   0.961808860301971435546875 [3F76391B]
0.576109707 [3F137BED]   0.57604229450225830078125  [3F137782]
0.442826569 [3EE2BA2A]   0.44434440135955810546875  [3EE3811C]
#155
(Правка: 2 фев 2022, 7:35) 18:55, 1 фев 2022

Допилил код, реализующий дробно-рациональную аппроксимацию. По точности немного уступает многочлену (этого не избежать), но по производительности оба варианта практически идентичны. Возможно, на каких-то процессорах будет преимущество по скорости, поэтому привожу код.
Среднеквадратическая относительная погрешность на 0.5<=x<2 составляет 4.97420936e-08.

_declspec(naked) float _vectorcall ln1(float x)
{
static const float ct[8] =          // Таблица констант
  {
    1.0f,                           // a0 = 1
    0.693147182f,                   // b0 = ln 2
    0.343713552f,                   // a2
    0.606786430f,                   // b2
    1.23005533f,                    // a1
    1.19918299f,                    // b1
    0.00755311921f,                 // a3
    0.0821869969f                   // b3
  };
  _asm
  {
    vmovd eax,xmm0                  // Прочитать в eax двоичное представление числа x
    mov edx,-126                    // В edx - смещение для коррекции порядка нормализованных float
    cmp eax,0x7F800000              // Сравнить число x со значением Inf
    jge ln_exit                     // Если x=+Inf или +NaN, вернуть аргумент x без изменений
    rol eax,9                       // Циклически сдвинуть число так, чтобы порядок оказался в al
    movzx ecx,al                    // Получить в ecx порядок числа со смещением +127
    jecxz ln_denorm                 // Перейти, если x=0 или x денормализованное (т.е. близко к 0)
      ln_calc:                      // Точка входа для вычисления логарифма после нормализации числа 
    cmp eax,0x55555600              // Сравнить мантиссу с 4/3
    setb al                         // al=1, если мантисса меньше 4/3; иначе al=0
    sub cl,al                       // ecx=k+126, где k - добавка к двоичному логарифму (целая часть)
    or al,126                       // Изменить порядок x так, чтобы привести к диапазону 2/3<=x<4/3
    add ecx,edx                     // ecx=k
    ror eax,9                       // Получить в eax двоичное представление числа и проверить его знак
    jc ln_nan                       // Перейти, если x<0 (аргумент некорректный); продолжить для x>0
    mov edx,offset ct               // В edx - адрес таблицы констант
    vmovd xmm0,eax                  // Перенести из eax в xmm0 приведённое значение x
    vmovups xmm1,[edx]              // xmm1 = b2 : a2 | b0 : a0
    vsubss xmm0,xmm0,xmm1           // xmm0 = x-1 = t
    vshufps xmm0,xmm0,xmm0,0        // xmm0 = t : t | t : t
    vcvtsi2ss xmm4,xmm4,ecx         // xmm4 = k
    vmulps xmm2,xmm0,xmm0           // xmm2 = t^2 : t^2 | t^2 : t^2
    vfmadd231ps xmm1,xmm0,[edx+16]  // xmm1 = b3*t+b2 : a3*t+a2 | b1*t+b0 : a1*t+a0
    vmovhlps xmm3,xmm3,xmm1         // xmm3 = b3*t+b2 : a3*t+a2
    vfmadd231ps xmm1,xmm2,xmm3      // xmm1 = b3*t^3+b2*t^2+b1*t+b0 : a3*t^3+a2*t^2+a1*t+a0 = Q : P
    vmovshdup xmm2,xmm1             // xmm2 = Q
    vdivss xmm1,xmm1,xmm2           // xmm1 = P/Q
    vfmadd213ss xmm0,xmm1,xmm4      // xmm1 = P/Q*t+k - готовый двоичный логарифм
    vmulss xmm0,xmm0,[edx+4]        // Перевести двоичный логарифм в натуральный
    ret                             // Вернуться
      ln_cont:                      // Продолжаем в штатном режиме, когда денормализованное x>0
    lea edx,[edx+2*ecx-31]          // Обновить смещение порядка для денормализованного числа
    neg ecx                         // 32-ecx - это необходимое количество сдвигов мантиссы влево
    shl eax,cl                      // Сформировать мантиссу нормализованного числа
    jmp ln_calc                     // Число нормализовано. Перейти для вычисления логарифма
      ln_denorm:                    // Точка входа для обработки денормализованного x, в т.ч. x=0
    btr eax,8                       // Прочитать в cf и сбросить знаковый бит числа
    bsr ecx,eax                     // Найти номер старшего единичного бита; установить zf, если x=0
    ja ln_cont                      // Перейти, если денормализованное x>0
    jz ln_inf                       // Если x=0, перейти для формирования результата -Inf
      ln_nan:                       // Точка входа для формирования результата NaN, если x<0 (или -NaN)
    sbb eax,eax                     // eax=-1 (значение NaN) в случае x<0
      ln_inf:                       // Точка входа для формирования результата -Inf в случае x=0
    or eax,0xFF800000               // Сформировать в eax значение -Inf, если x=0; или NaN, если x<0
    vmovd xmm0,eax                  // Перенести результат из eax в xmm0
      ln_exit:                      // В xmm0 готов результат для нештатного значения x
    ret                             // Вернуться
  }
}
Страницы: 16 7 8 9 10 11
ПрограммированиеФорумОбщее