Тёмная тема
- foxes
- Постоялец
FordPerfect
> Зачем 2 поля в таблице, если оба с ответом тупо суммируются?
Вот например:
s2-=s3; s2=s2*s1; s3=s3*s1-1.44140625; s3=s3+s2;
Сразу понятно что s1 можно за скобку вынести, а вычитание вообще убрать - а вот и нельзя, результат будет другим.
FordPerfect
> со своим (из #8 который)?
Скорость 190-65 тактов (memcpy) 26-31 такт (union), точность -9 бит, у меня 40-50 тактов, точность -4 бита (количество бит с конца от 52 бит мантиссы не совпадающие с результатом). Если сделать еще одну итерацию методом Ньютона, то будет точности -0 бит и скорость 40-50 тактов.
double _STDCALL _DsqrtDEF(double _a)
{
double y;
conv4b z;
unsigned _int32 zn;
z.valfloat=_a;
zn=0x5F3504F3-(z.valuint>>1);
z.valuint=zn;
y=0.5*z.valfloat*(3.0-z.valfloat*z.valfloat*_a);
y=0.5*y*(3.0-y*y*_a);
y=0.5*y*(3.0-y*y*_a);
return 0.5*_a*y*(3.0-y*y*_a);
}в принципе я также делал но решил заменить 3 итерации Нъютона на Герона.
double _STDCALL _DsqrtDEF(double _a)
{
if (_a==0) return 0;
conv4b c;
conv4b al;
c.valfloat=_a;
al.valuint=0x5f3504f3-(c.valuint>>1);
float g=al.valfloat;
float mr=g*(3.0f-c.valfloat*g*g);
double r=2.0/mr+_a*mr*0.5;
return r*0.25+_a/r;
}тут за счет меньшего количества операций при использовании в коде получается быстрее работает, а когда много умножений и процессор без конвейера получается медленнее.
Чтобы посчитать количество неверных бит нужно взять у чисел с плавающей точкой их бинарный вид и вычесть их друг из друга, двоичный логарифм от абсолютной разницы покажет на сколько бит они отличаются.
biterror=ilog2(*(unsigned __int64*)&val1-*(unsigned __int64*)&val2); //для double biterror=ilog2(*(unsigned __int32*)&val1-*(unsigned __int32*)&val2); //для float
не велика конечно потеря, но если делать 3 итерации Нъютона для абсолютной точности в x64 то скорость 60-70 тактов становиться, а для 2-х 50-60 (вместо 26-31), хотя это можно компенсировать если первые итерации делать на float там все равно только для первой части мантиссы биты вычисляются. У меня по прежнему 40-50.
foxes
Спасибо.
А можешь тогда ещё это сравнить:
double f_sqrt(double x) { double y,w; float z; uint32_t zn; z=x; memcpy( &zn,&z,4); zn=0x5F3504F3-( zn>>1); memcpy( &z,&zn,4); y=z*( 3.0-z*z*x); y=y*( 12.0-y*y*x); return 0.0001220703125*x*y*( 768.0-y*y*x); //y=y*(768.0-y*y*x); //return 9.094947017729282379e-013*x*y*(201326592.0-y*y*x); }
У меня выиграша не даёт (меньше умножений, больше загрузок констант), но вдруг?
также.
И так магическая строка:
s3=s3*s1-1.44140625+(s2-s3)*s1;
которая может быт упрощена до:
s3=s2*s1-1.44140625;
но смысл в том что в S3 хранится старшая часть мантиссы, а в (S2-S3) младшая, S2 содержит мантиссу полностью.
При таком вычислении результат получается совсем с другой точностью.
s3=s3*s1+(s2-s3)*s1;
foxes
Это понятно.
http://en.wikipedia.org/wiki/Kahan_summation_algorithm
Правка:
Ну вот с парными значениями в таблице такая же песня, там они разного порядка (-013 и +003 степени) как раз разница 52 бита.
Вот как раз то что не дает реализовать точный алгоритм вычисления:
> //Алгебраически, c всегда должна равняться нулю. Берегитесь слишком оптимизирующих компиляторов!
И теперь чтобы посчитать эту таблицу все что нужно это правильный алгоритм вычисления логарифма, где:
tablehi=TrueLog2-myLog2; tablelow=(TrueLog2-myLog2)-tablehi;
Замкнутый круг. Получается что алгоритма генерации такой таблицы нет, а лежит она у них также прямо в исходном коде.
Другими словами, по логике Intel, можно свести вариативность значений разницы между реальными значениями Log2(x) и некой F(x) к минимуму, записав эту разницу в таблицу. Шаг мантиссы X для этой таблицы составляет 0.001953125 в диапазоне от 1 до 2 всего 512 значений пар, общий размер 8208 байт.
В результате имеем:
Log2(x)=F(x)+Log2Table[IdTable(x)];
так можно любую функцию выразить с абсолютной точностью без F(x).
intelLog2Table - usable
алгоритм _mm_log2_pd
Код Intel AVX _mm_log2_pd
Свой вариант для одного значения X на SSE
Тема в архиве.
Тема закрыта.