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

Suslik
Не знаю, я для искусства считаю.

foxes

Ой.

Спасибо.

FordPerfect
Попробуй заменить memcpy на обычное присваивание. и 2ULL*m на (m<<1)
Я сделал стало по шустрей но не намного.
Во первых преобразование между целой и дробной частью идет через память это много жрет.
Вот такое вот выражение можно сделать короче на один int32
n=(n&0x000FFFFFFFFFFFFF);
n=(n|0x3FF0000000000000);
Если выполнять эту операцию только над старшей частью так как младшая остается неизменной.

Алгоритм от intel работает за 40-50 тактов c полным соответствием FPU. только там таблица используется на 8 килобайт как ее считать я не знаю, а вот таблицу выложить могу.
intelLog2Table

Код из #12 - фигня.

Вот более внятная версия:

double f_exp(double x)
{
    double y=((((((
        2.021099337273173939e-014)*x+
        7.761021455128987523e-012)*x+
        2.483526865641275904e-009)*x+
        6.357828776041666314e-007)*x+
        0.0001220703125)*x+
        0.015625)*x+
        1.0;
    y*=y;y*=y;y*=y;
    y*=y;y*=y;y*=y;
    return y;
}

double f_log(double x)
{
    double y,z;
    y=(((((((
        +0.09008290449235264)*x
        -0.8232520933634452)*x
        +3.268347597049261)*x
        -7.395442213971012)*x
        +10.54199156531427)*x
        -9.979689072844229)*x
        +6.882403741291307)*x
        -2.584443170188027;
    z=f_exp(-y);
    y-=1.0-x*z;
    z=f_exp(-y);
    y-=1.0-x*z;
    return y;
}

Тоже не подарок.

Точность:
~1e-16 - exp
~5e-17 - log.

Время:
80 тактов exp,
370 тактов log.

Но это на моей машине.

Спасибо.

Итак в чем собственно суть алгоритма Intel, для числа 1.35355339059327376220042218105242 результат 0.43675179543982401
Код для того чтобы посчитать числа очень близкие к 1.35355339059327376220042218105242 и их аналоги по мантиссе:

FordPerfect
> Код из #12
Не знаю в чем фишка твоего алгоритма но на значение x=1.5678 выдает 0.44967336277248604, а должно быть 0.64874153049172767 при этом считает за 120 тактов.

foxes

>Не знаю в чем фишка твоего алгоритма но на значение x=1.5678 выдает 0.44967336277248604, а должно быть >0.64874153049172767 при этом считает за 120 тактов.

Там ln(), не log2().

>при этом считает за 120 тактов

а это радует.

FordPerfect
> а это радует.
ага а если поставить  x=14.5678 то за 9000 тактов. а так у тебя ж в диапазоне... точность кстати высокая даже сравнимо с FPU полностью 1-2 бита разницы.

foxes
Ага, и результат inf.
Там сказано, что функция заточена на значения x∈[e^(-1/2),e^(+1/2)] ([-1/2;+1/2] для exp).
Экспоненту надо откусывать отдельно, из битового представления числа, например.

Да я понял уже

foxes
>Итак в чем собственно суть алгоритма Intel

Так, я торможу, или table lookup не параллелится (ну не считая AVX2 gathered loads)?
Т. е. _mm_log2_pd внутре скалярная?
(Или она AVX2-only?)

>Во первых преобразование между целой и дробной частью идет через память это много жрет.
Переделал эту функцию в целых, всё-равно хреново:

inline double f(double x)
{
    double ret;
    uint64_t n,m,m0;
    uint64_t n0,n1,k;
    memcpy(&n,&x,8);
    m0=((n&0x7FF0000000000000)>>11ULL);
    n=(n&0x000FFFFFFFFFFFFF);
    n<<=10ULL;
    m=0ULL;
    for(int i=0;i<53;++i)
    {
        k=(((n>>62ULL)|(n>>63ULL))&1ULL);
        m=m+m+k;
        n=(n>>k)-(k<<61ULL);
        n1=(n>>32ULL);
        n0=(n&0xFFFFFFFF);
        n=n+n+(((n1*n1)<<2)+((n1*n0)>>29)+((n0*n0)>>62));
    }
    m=(m&0x000FFFFFFFFFFFFF);
    m=(m|0x3FF0000000000000);
    memcpy(&x,&m,8);
    ret=x;
    m0=(m0|0x40A0000000000000);
    memcpy(&x,&m0,8);
    ret=ret+x-3072.0;
    return ret;
}

Замеришь, ради интереса?

+((n0*n0)>>62) можно выкинуть без особой потери для точности.

FordPerfect
> Т. е. _mm_log2_pd внутре скалярная?
> (Или она AVX2-only?)
Да, реализацию которую я разбирал на AVX, но я сделал для себя на SSE тоже, хотя точно утверждать что то сложно поскольку прежде чем она вызовется происходит определение типа процессора. Там заточка под вычисление сразу двух логарифмов, а таблицы не параллельно считываются, поскольку там из всей этой таблицы всего две пары значений считываются. Разобраться бы как они таблицу вычислили, загружать такую в код жирно будет.

FordPerfect
> Переделал эту функцию в целых, всё-равно хреново:
Ну не ты тут только хуже сделал целочисленное умножение долго, я имел ввиду через юнионы копировать - если для x64 или SSE2 такой алгоритм будет без использования памяти работать.
2000 тактов.

foxes
>2000 тактов.

Спасибо.

Здесь точно всё правильно?
Зачем 2 поля в таблице, если оба с ответом тупо суммируются? Можно было прямо в таблице просуммировать.

FordPerfect
Тут очень много нюансов связанных с точностью
Допустим строка:

  double s6=((((((-0.026810595806536702)*s3+0.046373904173250588)*s3 
                        -0.083554326761742681)*s3+0.16058097174555816)*s3 
                        -0.34719362445817492)*s3+0.00089412050833233693)*s3;

вычисляется правильно только на FPU, благодаря большой разрядности, если те же вычисления (в той же последовательности) реализовывать на SSE double 64, то результат разъедется. Если как следует приглядеться тут много парадоксов, которые я постарался переложить в обычный код, но без них не считается правильно.

А можешь сравнить вот этот код

double f_sqrt(double x)
{
    double y;
    float z;
    uint32_t zn;
    z=x;
    memcpy(&zn,&z,4);
    zn=0x5F3504F3-(zn>>1);
    memcpy(&z,&zn,4);
    y=0.5*z*(3.0-z*z*x);
    y=0.5*y*(3.0-y*y*x);
    return 0.5*x*y*(3.0-y*y*x);
}

со своим (из #8 который)?
Точность, скорость.
На своей машине (на моей не очень показательно).

memcpy можешь поменять на union.