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

Aenigma
> Я не совсем понял, по какой формуле вы аппроксимируете. У вас вижу почему-то 6 констант. У меня 8 констант, но одна — свободная, т.е. фактически 7.

    float z=x-1.0f;
    float P=1.44269504f+z*(1.41014671f+z*0.247352005f);
    float Q=x+z*(0.4774388f+z*(0.576839184f+z*0.045954678f));
    return y+z*(P/Q);

Я аппроксимирую log2(x)/(x-1) как P(x-1)/Q(x-1). Т. к. x=z+1, то по сути Q=1.0+z*(1.4774388+z*(0.576839184+z*0.045954678)). Через x чуть точнее.
И z*(P/Q) - чуть точнее, чем (z*P)/Q.

> И, кстати, точность у меня получилось лишь незначительно хуже (точные цифры пока не готов привести), чем через P((x-1)/(x+1)). У вас на сколько хуже точность?
P((x-1)/(x+1)) ( https://gamedev.ru/code/forum/?id=132465&page=6&m=5505593#m80 ):

  Linf: abs=1.11e-07, rel=2.46e-07, 3.479 ulp
  RMS : abs= 2.1e-08, rel=4.72e-08, 0.581 ulp

P(x)/Q(x) ( https://gamedev.ru/code/forum/?id=132465&page=9&m=5508853#m133 ):

  Linf: abs=1.29e-07, rel=3.03e-07, 4.232 ulp
  RMS : abs=2.24e-08, rel=5.25e-08, 0.644 ulp

Это без FMA.
С FMA (тупо включить -mfma в настройках, GCC по умолчанию автоматом преобразует a*b+c во fma(a,b,c)):
P((x-1)/(x+1)):

  Linf: abs=1.06e-07, rel=2.13e-07, 3.001 ulp
  RMS : abs=1.95e-08, rel=4.53e-08, 0.558 ulp

P(x)/Q(x):

  Linf: abs=1.14e-07, rel=2.28e-07, 3.396 ulp
  RMS : abs=2.05e-08, rel=   5e-08, 0.613 ulp

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

У меня среднеквадратическая относительная погрешность на 0.5<=x<2 составляет 4.59E-8 для аппроксимации многочленом и 4.96E-8 - для дробно-рациональной аппроксимации (степени P и Q равны 3). Различие 8%. У вас различие 11 %. Побольше, чем у меня. Либо у меня аппроксимация многочленом хуже вашей (кстати, в какой-то мере это может быть, потому что я не занимался тонкой подстройкой коэффициентов), либо у вас хуже аппроксимация дробно-рациональной функцией.

Aenigma
> Вы не смотрели теоретический график погрешности вашей аппроксимации?
Вот (формула из #130; в #133 числитель и знаменатель поделены на 0.693147180):
https://rextester.com/SPBCH54796
(нажать F8 чтобы увидеть график)

Aenigma
> Как я понял, у вас числитель аппроксимируется многочленом второй степени, знаменатель - третьей
Т. к. я домножаю P на z, то в каком-то смысле числитель - тоже 3-й.

Я имел в виду степень без умножения на (x-1), я тоже частное потом домножаю на (x-1). Похоже небольшая систематическая составляющая ошибки у вас имеется. Относительная погрешность почти доходит до 4E-8, это многовато. У меня относительная погрешность не превышает примерно 2.5E-9.

Удалось переплюнуть:
__m512 _mm512_log2_ps (__m512 a)
__m256 _mm256_log2_ps (__m256 a)
?

dds
Что-то они как-то не зажгли:
https://godbolt.org/z/WjhP6jvnY
Я криво мерял?
Да, там только скалярно, но даже если умножить на 16 - всё-равно медленнее.

> Я криво мерял?
Вижу в коде

call      *__svml_log2f16@GOTPCREL(%rip)                #52.7

а не

vlog2ps

Насколько я понимаю, она есть на Knights Corner/Xeon Phi. Но я даже не особо понимаю, как заставить компилятор её выдавать.

FordPerfect
Как то так:
https://godbolt.org/z/3Pv9jcref

dds
Забавно. Всё-равно call      *__svml_log2f16, но теперь в 100 раз быстрее.
0.7 ns на значение, т. е. примерно 11 ns на функцию.

Да, #144 выглядит более похожим на правду.
Всё-равно не особо понятно, почему #142 медленнее в 100 раз, а не в 16.

Если забить на отрицательные, денормалы, inf/nan - моё тоже можно легко SIMD'ифицировать. Если не забивать - сложнее, но тоже делается.
По идее если сделать особые случаи через y

float log2_v9(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.88539003e+00f
           +x2*(9.61808794e-01f
           +x2*(5.76110615e-01f
           +x2*(4.42826415e-01f))));
}

легче SIMD'ифицировать.

А, тест в #142 всё-равно чушь

    int N=100000;
    volatile float src[1024]={0.0f};
    volatile float dst[1024]={0.0f};
    double t=clock();
    for(int i=0;i<N;++i)
        dst[i&1023]=f(src[i&1023]);
    t=clock()-t;

т. к. там замеряется скорость log2(0.0).

Вот замер скорости для #146:
https://rextester.com/GIIOLT84405