Войти
ФлеймФорумПрограммирование

Буст pow2

Страницы: 1 2 3 Следующая »
#0
6:06, 5 мар. 2018

Хочется написать аппроксимацию pow2 (в стандарте C известной как exp2f).
Быструю. На x86 CPU.

Быстрая - означает быстрее деления.

Желательно - tableless, branchless; чтобы хорошо SIMD'ифицировалась.

На поведение в левых областях - можно забить (считать |x|<50, например).

Первое, что приходит в голову:

inline float pow2(float x)
{
    return bits2float(float2bits(x+639.0f)<<9);
}
Но точность не радует - это тупо lerp между 2^n.
Но даже при этом - относительная погрешность меньше 7%.

Мысль по улучшению выглядит примерно так:

inline float pow2(float x)
{
    uint32_t bits=float2bits(x+383.0f)<<8;
    float y=bits2float(bits);
    float z=bits2float(~bits);
    float w=-0.1177399f*(12.4932976f+y*z);
    return y*w;
}
Но всё-таки многовато операций.

Код:
http://rextester.com/ISG99233

Т. к. на rextester скорость дико прыгает - просьба позапускать у себя, кому не жалко.

Идеи по альтернативным версиям - приветствуются.


#1
8:40, 5 мар. 2018

По здравом размышлении наверно так:

inline float pow2(float x)
{
    uint32_t bits=float2bits(x+383.0f)<<8;
    float y=bits2float(bits);
    float z=bits2float(~bits);
    return (-0.1177399f*y)*(12.4932976f+y*z);
}
чтобы цепочка зависимых умножений была длины 2, а не 3.
И на FMA хорошо ложится.

#2
12:24, 5 мар. 2018
float pow2(float x)
{
    uint32_t bits=float2bits(x+383.0f)<<8;
    float y=bits2float(bits);
    float z=bits2float(~bits);
    return (-0.1177399f*y)*(12.4932976f+y*z);
}
Блин.
pow2(float):
        addss   .LC0(%rip), %xmm0
        movd    %xmm0, %eax
        sall    $8, %eax
        movl    %eax, -4(%rsp)
        notl    %eax
        movss   -4(%rsp), %xmm1
        movl    %eax, -4(%rsp)
        movss   -4(%rsp), %xmm0
        mulss   %xmm1, %xmm0
        mulss   .LC2(%rip), %xmm1
        addss   .LC1(%rip), %xmm0
        mulss   %xmm1, %xmm0
        ret
Я всего-то хочу:
pow2_asm(float):
        addss pow2_asm(float)::c0(%rip),%xmm0
        pslld $8,%xmm0
        pcmpeqb %xmm1,%xmm1
        pxor %xmm0,%xmm1
        mulps %xmm0,%xmm1
        mulps pow2_asm(float)::c2(%rip),%xmm0
        addps pow2_asm(float)::c1(%rip),%xmm1
        mulps %xmm1,%xmm0

Код:
http://rextester.com/YLW88140

Просьба мерить - всё-ещё в силе.

#3
12:32, 5 мар. 2018
uint32_t bits;
float y;
float z;

#define pow2(x, res) \
{ \
    bits=float2bits(x+383.0f)<<8; \
    y=bits2float(bits); \
    z=bits2float(~bits); \
    res = (-0.1177399f*y)*(12.4932976f+y*z); \
}
#4
13:13, 5 мар. 2018

lookid
Гм, а это к чему? Думаешь, что компилятор не заинлайнит?

#5
13:47, 5 мар. 2018

умножил К на 10 чтоб меньше результаты плясали
i7-7700HQ, gcc version 7.2.0
gcc pow2.cpp -O3

Function      |         Tacts/call|   Max. rel. error
--------------+-------------------+------------------
movss         |                3.0|                 -
divss         |                3.3|                 -
std::exp2     |               86.2|                 0
pow2_try1     |                3.1|            0.0615
pow2_try2     |                4.1|            0.0054
pow2_try3     |                4.3|            0.0054
pow2_asm      |                4.3|            0.0054
gcc pow2.cpp -O3 -march=native
Function      |         Tacts/call|   Max. rel. error
--------------+-------------------+------------------
movss         |                3.0|                 -
divss         |                3.2|                 -
std::exp2     |               86.4|                 0
pow2_try1     |                3.3|            0.0615
pow2_try2     |                3.6|            0.0054
pow2_try3     |                3.7|            0.0054
pow2_asm      |                4.2|            0.0054

#6
13:56, 5 мар. 2018

kipar
Спасибо.

В 3-тактное деление поверить сложно.
Возможно -O3 включает ffast-math и автовекторизацию.

Можешь запостить дизасм?

#7
14:07, 5 мар. 2018

FordPerfect
с -O2 результаты те же самые.
дизасм эээ видимо вот так:

_Z9the_divssf:
.LFB5338:
  .seh_endprologue
  vmovss  .LC6(%rip), %xmm1
  vdivss  %xmm0, %xmm1, %xmm0
  ret
#8
14:12, 5 мар. 2018

Я в Хулионе для тумана целую экспоненту применял.

delta := 256 - (256 - (delta and $7F)) shr (delta shr 7);

#9
14:13, 5 мар. 2018

Хотя, похоже на правду.
Переписал divss() на асме, скорость не поменялась.
Так что видимо не автовекторизация.

Тогда занятно.

Что это? Очень быстрое деление?
Или их несколько параллельно выполняется? Я думал, что в ядре одна делилка (а складывалок/умножалок/модулей-FMA - больше).
Или это на rdtsc такой множитель?

#10
14:18, 5 мар. 2018

kipar
> с -O2 результаты те же самые.
Ясно. Подтверждает гипотезу, что не в автовекторизации дело.
Дизасм интересен именно цикла, а не функции, т. к. предположительно она заинлайнится и вызова именно этого того кода всё-равно не случилось.

1 frag / 2 deaths
Диапазон (и тип) delta?

#11
14:23, 5 мар. 2018

FordPerfect
целое положительное
только это exp(-x)

#12
14:28, 5 мар. 2018

FordPerfect

+ Показать
#13
14:34, 5 мар. 2018

kipar
Спасибо. Да, скалярные:

.L7:
  vdivss  %xmm0, %xmm4, %xmm2
  subl  $1, %eax
  vaddss  %xmm3, %xmm0, %xmm0
  vaddss  %xmm2, %xmm1, %xmm1
  jne  .L7
#14
14:48, 5 мар. 2018

1 frag / 2 deaths

delta := 256 - LUT[delta and $7F] shr (delta shr 7);
по идее точнее и в ту же скорость.
Страницы: 1 2 3 Следующая »
ФлеймФорумПрограммирование

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