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
Тема в архиве.