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

Log2 (_mm_log2_pd intel - ламеры :)) (3 стр)

Страницы: 1 2 3
#30
0:00, 26 ноя 2014

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.

#31
1:34, 26 ноя 2014

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);
}

У меня выиграша не даёт (меньше умножений, больше загрузок констант), но вдруг?

#32
1:41, 26 ноя 2014

также.

#33
1:59, 26 ноя 2014

И так магическая строка:

s3=s3*s1-1.44140625+(s2-s3)*s1;

которая может быт упрощена до:

s3=s2*s1-1.44140625;

но смысл в том что в S3 хранится старшая часть мантиссы, а в (S2-S3) младшая, S2 содержит мантиссу полностью.
При таком вычислении результат получается совсем с другой точностью.

s3=s3*s1+(s2-s3)*s1;
#34
2:08, 26 ноя 2014

foxes
Это понятно.
http://en.wikipedia.org/wiki/Kahan_summation_algorithm

Правка:

http://ru.wikipedia.org/wiki/Алгоритм_Кэхэна

#35
2:15, 26 ноя 2014

Ну вот с парными значениями в таблице такая же песня, там они разного порядка (-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

+ Показать
Страницы: 1 2 3
ПрограммированиеФорумОбщее

Тема в архиве.

Тема закрыта.