Быстрый точный sin/cos из float32. (6 стр) / Общее / Форум / Программирование игр / GameDev.ru

На всякий случай, сохраню здесь:

Код на Maxima, считающий оптимальный в L²-смысле полином, приближающий синус, с коэффициентами, прибитыми гвоздями к представимым во float32 числам (в примере - на \(\left[ -\pi ;\; +\pi \right]\), коэффициент при x гарантированно 1):

+ Показать

− Скрыть

fpprec:256;
fpprintprec:17;
ratprint:false;

l:-%pi$
r:+%pi$

log2floor(x):=floor(log(abs(x))/log(2));
snap2float(x):=floor(x*2^(-log2floor(x)+23)+1/2)/2^(-log2floor(x)+23);

poly(x,n):=x+apply("+",makelist(a[i]*x^(2*i+1),i,1,n/2))$
pdst(f,n):=integrate((f(x)-poly(x,n))^2,x,l,r)$
coef(f,b,n):=flatten(solve(flatten([makelist(diff(pdst(f,n),a[i])=0,i,length(b)+1,n/2),makelist(a[i]=b[i],i,1,length(b))])))$
appr(f,b,n):=apply(ev,append([poly(x,n)],coef(f,b,n)))$
prec(f,b,n):=apply(ev,append([pdst(f,n)],coef(f,b,n)))$

snapappr(f,b,n):=if(length(b)+1<n/2) then snapappr(f,append(b,[snap2float(ev(a[length(b)+1],coef(f,b,n)))]),n) else appr(f,b,n);

p:snapappr(sin,[],9);
bfloat(p);

Код на Maxima, выделяющий часть (двоичных) знаков числа \(\pi\), в виде вещественного числа:

+ Показать

Точное (десятичное) значение float32-числа по его битам (исключая денормализованные и т. п.):

bits2decimal(n):=bfloat((1+mod(n,2^23)/2^23)*2^(mod(floor(n/2^23),2^8)-127));

В качестве онлайн-Maxima использую http://maxima.cesga.es/ .

FordPerfect
Постоялец

#76
18:28, 1 сен 2015

Ну и быстрые версии, для тех, кому плевать на точность:

http://rextester.com/RIP91276

Bishop
Постоялец

#77
20:08, 1 сен 2015

FordPerfect
На инстриктах читаемость кода упала в разы :(

FordPerfect
Постоялец

#78
20:32, 1 сен 2015

Bishop
По сравнению с C++ версией или с asm?

Ну а вообще "любители MSVC должны страдать".
Портабельность, всё такое.

PaulSh
Постоялец

#79
20:34, 1 сен 2015

FordPerfect
> Просьба потестить.

Win7 x64, i7-4710HQ
VS2015, Release x64, макс оптимизация (под спойлером),

sin_f32_sse2_intrin на 20% медленнее предыдущего С++ аналога (sin_f32)
И это я ещё вынес объявления статиков из тела функции, а то скорость была 2 раза хуже sin_f32.

Дисассемблер и оптимизации:

+ Показать

Я так понимаю, VS забивает на твои интринзики, переписывает по-своему и, в итоге, не может хорошо оптимизировать. Может я что-то напутал в параметрах проекта?
Код тот же, что и раньше:

+ Показать

− Скрыть

#include <stdio.h>
#include <iostream>
#include <chrono>
#include <conio.h>
#include <string>
#include <sstream>
#include <cmath>
#include <vector>

#include "Eigen/Dense"
#include "test.h"

#include <Windows.h>

#include "QuickTrig.h"




//==============================================================================

// C++ versions.

// |x|<65536:       rel. error < 1.9*10^-7 [ 2 ulp ]
// |x|<262144:      abs. error < 4*10^-6
// Garbage for |x|>8388608.
// Values are for SSE case.
// FPU has both better range and accuracy.
inline float sin_f32(float x)
{
    float y = 0.318309886f*x; // 1/(pi)

                              // Round y to closest integer.
    y = y + 12582912.0f; // 1.5*2^23

    uint32_t bits;
    memcpy(&bits, &y, 4);
    float sign;
    bits = (bits << 31) | 0x3F800000;
    memcpy(&sign, &bits, 4);
#if 0 // Not necessary if floating-point math is guaranteed  to be handled by SSE ( -mfpmath=sse ).
    // Fixup for FPU carrying calculation
    // at 80 bits.
    asm("fstp %[bits]\n\t"
        "fld %[bits]\n\t"
        :"+t"(y)
        : [bits]"m"(bits)
        : "cc");
#endif

    y = y - 12582912.0f;

    // Reduced argument.
    // pi is split into 4 parts, Kahan-summation-style.
    // Cf. http://www.netlib.org/cephes/

    x = (((x
        - y*3.140625f)                       // 12 bits, 9 non-zeros
        - y*9.670257568359375e-4f)           // 10 bits, 9 non-zeros
        - y*6.2771141529083251953125e-7f)    //  9 bits
        - y*1.2154201013012385202950e-10f;   // 23 bits

    float x2 = x*x;

    return sign*(
        x
        + x*x2*(-1.6666658222675323e-1f    // -11184805*2^-26
            + x2*(+8.3330469205975532e-3f      // +8947541*2^-30
                + x2*(-1.9808833894785493e-4f      // -13612527*2^-36
                    + x2*(2.6050554015241175e-6f)))));
}

// |x|<32768:       rel. error < 1.9*10^-7 [ 2 ulp ]
// |x|<131072:      abs. error < 1*10^-6
// Garbage for |x|>8388608.
// Values are for SSE case.
// FPU has both better range and accuracy.
inline float cos_f32(float x)
{
    float y = 0.318309886f*x + 0.5f; // 1/(pi)

                                     // Round y to closest integer.
    y = y + 12582912.0f; // 1.5*2^23

    uint32_t bits;
    memcpy(&bits, &y, 4);
    float sign;
    bits = (bits << 31) | 0x3F800000;
    memcpy(&sign, &bits, 4);
#if 0 // Not necessary if floating-point math is guaranteed  to be handled by SSE ( -mfpmath=sse ).
    // Fixup for FPU carrying calculation
    // at 80 bits.
    asm("fstp %[bits]\n\t"
        "fld %[bits]\n\t"
        :"+t"(y)
        : [bits]"m"(bits)
        : "cc");
#endif

    y = y - 12582912.0f;
    y -= 0.5f;

    x = (((x
        - y*3.140625f)                       // 12 bits, 9 non-zeros
        - y*9.670257568359375e-4f)           // 10 bits, 9 non-zeros
        - y*6.2771141529083251953125e-7f)    //  9 bits
        - y*1.2154201013012385202950e-10f;   // 23 bits

    float x2 = x*x;

    return sign*(
        x
        + x*x2*(-1.6666658222675323e-1f    // -11184805*2^-26
            + x2*(+8.3330469205975532e-3f      // +8947541*2^-30
                + x2*(-1.9808833894785493e-4f      // -13612527*2^-36
                    + x2*(2.6050554015241175e-6f)))));
}/**/



//==============================================================================

// Intrinsic version (scalar).

// |x|<65536:       rel. error < 1.9*10^-7 [ 2 ulp ]
// |x|<262144:      abs. error < 4*10^-6
// Garbage for |x|>8388608.

static const __m128 c0 = _mm_set_ss(0.318309886f);
static const __m128 c1 = _mm_set_ss(12582912.0f);
static const __m128 c2 = _mm_set_ss(3.140625f);
static const __m128 c3 = _mm_set_ss(9.670257568359375e-4f);
static const __m128 c4 = _mm_set_ss(6.2771141529083251953125e-7f);
static const __m128 c5 = _mm_set_ss(1.2154201013012385202950e-10f);
static const __m128 c6 = _mm_set_ss(+2.6050554015241175e-6f);
static const __m128 c7 = _mm_set_ss(-1.9808833894785493e-4f);
static const __m128 c8 = _mm_set_ss(+8.3330469205975532e-3f);
static const __m128 c9 = _mm_set_ss(-1.6666658222675323e-1f);

inline float sin_f32_sse2_intrin(float x)
{
    __m128 x_sse = _mm_load_ss(&x);
    __m128 y = _mm_mul_ss(c0, x_sse);

    // Round y to closest integer.
    y = _mm_add_ss(y, c1);
    __m128 sign = y;
    y = _mm_sub_ss(y, c1);

    sign = _mm_castsi128_ps(_mm_slli_epi32(_mm_castps_si128(sign), 31));

    x_sse = _mm_sub_ss(x_sse, _mm_mul_ss(c2, y));
    x_sse = _mm_sub_ss(x_sse, _mm_mul_ss(c3, y));
    x_sse = _mm_sub_ss(x_sse, _mm_mul_ss(c4, y));
    x_sse = _mm_sub_ss(x_sse, _mm_mul_ss(c5, y));

    __m128 x2 = _mm_mul_ss(x_sse, x_sse);

    __m128 ret_sse = _mm_mul_ss(c6, x2);
    ret_sse = _mm_mul_ss(_mm_add_ss(c7, ret_sse), x2);
    ret_sse = _mm_mul_ss(_mm_add_ss(c8, ret_sse), x2);
    x_sse = _mm_xor_ps(x_sse, sign);
    ret_sse = _mm_mul_ss(_mm_add_ss(c9, ret_sse), x2);
    ret_sse = _mm_add_ss(x_sse, _mm_mul_ss(ret_sse, x_sse));
    float ret;

    _mm_store_ss(&ret, ret_sse);
    return ret;
}

/*
// |x|<32768:       rel. error < 1.9*10^-7 [ 2 ulp ]
// |x|<131072:      abs. error < 1*10^-6
// Garbage for |x|>8388608.
inline float cos_f32_sse2_intrin(float x)
{
    static const __m128 c0 = _mm_set_ss(0.318309886f);
    static const __m128 c1 = _mm_set_ss(12582912.0f);
    static const __m128 c2 = _mm_set_ss(3.140625f);
    static const __m128 c3 = _mm_set_ss(9.670257568359375e-4f);
    static const __m128 c4 = _mm_set_ss(6.2771141529083251953125e-7f);
    static const __m128 c5 = _mm_set_ss(1.2154201013012385202950e-10f);
    static const __m128 c6 = _mm_set_ss(+2.6050554015241175e-6f);
    static const __m128 c7 = _mm_set_ss(-1.9808833894785493e-4f);
    static const __m128 c8 = _mm_set_ss(+8.3330469205975532e-3f);
    static const __m128 c9 = _mm_set_ss(-1.6666658222675323e-1f);
    static const __m128 c10 = _mm_set_ss(0.5f);

    __m128 x_sse = _mm_load_ss(&x);
    __m128 y = _mm_mul_ss(c0, x_sse);

    // Round y to closest integer.
    y = _mm_add_ss(y, c10);
    y = _mm_add_ss(y, c1);
    __m128 sign = y;
    y = _mm_sub_ss(y, c1);
    y = _mm_sub_ss(y, c10);

    sign = _mm_castsi128_ps(_mm_slli_epi32(_mm_castps_si128(sign), 31));

    x_sse = _mm_sub_ss(x_sse, _mm_mul_ss(c2, y));
    x_sse = _mm_sub_ss(x_sse, _mm_mul_ss(c3, y));
    x_sse = _mm_sub_ss(x_sse, _mm_mul_ss(c4, y));
    x_sse = _mm_sub_ss(x_sse, _mm_mul_ss(c5, y));

    __m128 x2 = _mm_mul_ss(x_sse, x_sse);

    __m128 ret_sse = _mm_mul_ss(c6, x2);
    ret_sse = _mm_mul_ss(_mm_add_ss(c7, ret_sse), x2);
    ret_sse = _mm_mul_ss(_mm_add_ss(c8, ret_sse), x2);
    x_sse = _mm_xor_ps(x_sse, sign);
    ret_sse = _mm_mul_ss(_mm_add_ss(c9, ret_sse), x2);
    ret_sse = _mm_add_ss(x_sse, _mm_mul_ss(ret_sse, x_sse));
    float ret;

    _mm_store_ss(&ret, ret_sse);
    return ret;
}/**/



int main()
{
    CQuickTrig trig;

    volatile float aval = 0.0;
    volatile float vals[40000000];
    float delta = 0.0f, pdelta = 0.0f;


    for (int i = 0; i < 40000000; i++)
    {
        vals[i] = -1000.f + 1.0f * i / 20000;
    }

    
    LARGE_INTEGER Freq, T1, T2;
    double TicksToSecKoef = 0.0;
    

    QueryPerformanceFrequency(&Freq);
  TicksToSecKoef = 1.0 / Freq.QuadPart;
    
    QueryPerformanceCounter(&T1);
    
    
    for (int j = 0; j < 20; j++)
    for (int i = 0; i < 40000000; i++)
    {
        //aval = sin_f32_sse2_intrin(vals[i]);
        aval = sin_f32(vals[i]);        
        //aval = trig.QuickSin(vals[i]);     
        //aval = cos(vals[i]);

        /*float val = sin(vals[i]);
        float lde = fabs(val - aval);

        if (lde > delta)
        {
            delta = lde;
            pdelta = fabs(delta * 100.0f / val);
        }/**/
    }
    
    QueryPerformanceCounter(&T2);
    

    
    std::stringstream ss;
    ss << "t = " << (double)(T2.QuadPart - T1.QuadPart) * TicksToSecKoef << ", d = " << delta << ", " << pdelta << "%";
    std::string str = ss.str();
    std::cout << str;

    _getch();

    return 0;
}

Bishop
Постоялец

#80
21:18, 1 сен 2015

FordPerfect
> По сравнению с C++ версией или с asm?
С asm конечно.

FordPerfect
Постоялец

#81
21:43, 1 сен 2015

PaulSh
Как интере~есно...
Спасибо, кстати.

>И это я ещё вынес объявления статиков из тела функции, а то скорость была 2 раза хуже sin_f32.
Ох ты блин. И в gcc тоже. Т. е. оно не умеет вынести _mm_set в compile-time и честно пишет вот такой ужас:
http://www.everfall.com/paste/id.php?g7yldd8go81p

То, что компиляторы интринсики уровня <emmintrin.h> могут реализовать через AVX - для меня новость (это как бы намекает на уровень знания мной интринсиков). Проверил - gcc, icc и clang делают то же самое.

И, похоже все компиляторы решают константы сначала загрузить в регистр, а уже потом делать арифметику. В смысле, в случае интринсик. И для данного конкретного кода это, похоже, получается тормознее.

Вот, ради интереса, что пишут gcc и icc, если вынести константы отдельно (код инициализации констант - там же):
gcc: http://www.everfall.com/paste/id.php?qlhorvlyhs8k
icc: http://www.everfall.com/paste/id.php?p2gi754bx6h8

Для C++ случая, без интринсик (но с SSE2-арифметикой) - они вполне себе складывают/множат прямо из памяти:
http://www.everfall.com/paste/id.php?6sp5nwfmo0am (gcc)

Тут впрочем, часть операций вообще идёт не через SSE.
Но для 4-ствольной версии такой вариает не идёт же?

Bishop
Дело вкуса, наверно.
Мне asm читать, пожалуй, тоже симпатичнее.
Механически переносить C++ код на интринсики похоже, что проще, чем на asm.

FordPerfect
Постоялец

#82
23:20, 1 сен 2015

А если у тебя много КЦ синусов, то можно написать вот так:

+ Показать

− Скрыть

inline void sin_f32_sse2_intrin_stream(const float *xs,float *ys,int n)
{
    if(n<=0) return;
    int n4=n/4;
    int nd=n-4*n4;
    for(int i=0;i<n4;++i)
    {
        __m128 x=_mm_loadu_ps(xs+4*i);
        __m128 y=_mm_mul_ps(c0,x);
        y=_mm_add_ps(y,c1);
        __m128 sign=y;
        y=_mm_sub_ps(y,c1);
        x=_mm_sub_ps(x,_mm_mul_ps(c2,y));
        x=_mm_sub_ps(x,_mm_mul_ps(c3,y));
        sign=_mm_castsi128_ps(_mm_slli_epi32(_mm_castps_si128(sign),31));
        x=_mm_sub_ps(x,_mm_mul_ps(c4,y));
        x=_mm_sub_ps(x,_mm_mul_ps(c5,y));
        x=_mm_xor_ps(x,sign);
        _mm_storeu_ps(ys+4*i,x);
    }
    for(int i=0;i<n4;++i)
    {
        __m128 x=_mm_loadu_ps(ys+4*i);
        __m128 x2=_mm_mul_ps(x,x);
        __m128 y=_mm_mul_ps(c6,x2);
        y=_mm_mul_ps(_mm_add_ps(c7,y),x2);
        y=_mm_mul_ps(_mm_add_ps(c8,y),x2);
      y=_mm_mul_ps(_mm_add_ps(c9,y),x2);
        y=_mm_add_ps(x,_mm_mul_ps(y,x));
        _mm_storeu_ps(ys+4*i,y);
    }
    for(int i=0;i<nd;++i)
        ys[4*n4+i]=sin_f32_sse2_intrin(xs[4*n4+i]);
}

И тогда константы в циклах замечательно лягут в регистры (у gcc, по крайней мере), и будет шустрее (по моим замерам).

http://rextester.com/SGW57300 - GCC.
http://rextester.com/WOJGS76842 - MSVC.

Полторы наносекунды на синус, хм.

FordPerfect
Постоялец

#83
0:08, 2 сен 2015

PaulSh
Ну... я целюсь в SSE2.
А так да: __m256-версия - логичный следующий шаг.

FordPerfect
Постоялец

#84
0:35, 2 сен 2015

PaulSh

Написал такое:

+ Показать

− Скрыть

// This code is in the public domain. Where that dedication is not
// recognized, you are granted a perpetual, irrevocable license to copy
// and modify it as you see fit.

// WARNING 1: the code relies on rounding mode being round-to-nearest!
// WARNING 2: the error bounds below are not rigorously tested.

#include <cstdio>
#include <cmath>
#include <cstdint>
#include <cstring>
#include <ctime>
#include <immintrin.h>

static const __m256 c0=_mm256_set1_ps(0.318309886f);
static const __m256 c1=_mm256_set1_ps(12582912.0f);
static const __m256 c2=_mm256_set1_ps(3.140625f);
static const __m256 c3=_mm256_set1_ps(9.670257568359375e-4f);
static const __m256 c4=_mm256_set1_ps(6.2771141529083251953125e-7f);
static const __m256 c5=_mm256_set1_ps(1.2154201013012385202950e-10f);
static const __m256 c6=_mm256_set1_ps(+2.6050554015241175e-6f);
static const __m256 c7=_mm256_set1_ps(-1.9808833894785493e-4f);
static const __m256 c8=_mm256_set1_ps(+8.3330469205975532e-3f);
static const __m256 c9=_mm256_set1_ps(-1.6666658222675323e-1f);
static const __m256 c10=_mm256_set1_ps(0.5f);

inline float sin_f32_avx_intrin(float x)
{
  __m256 x_avx=_mm256_set1_ps(x);
  __m256 y=_mm256_mul_ps(c0,x_avx);

    // Round y to closest integer.
  y=_mm256_add_ps(y,c1);
  __m256 sign=y;
  y=_mm256_sub_ps(y,c1);

  sign=_mm256_castsi256_ps(_mm256_slli_epi32(_mm256_castps_si256(sign),31));

  x_avx=_mm256_sub_ps(x_avx,_mm256_mul_ps(c2,y));
  x_avx=_mm256_sub_ps(x_avx,_mm256_mul_ps(c3,y));
  x_avx=_mm256_sub_ps(x_avx,_mm256_mul_ps(c4,y));
  x_avx=_mm256_sub_ps(x_avx,_mm256_mul_ps(c5,y));
  
    __m256 x2=_mm256_mul_ps(x_avx,x_avx);

  __m256 ret_avx=_mm256_mul_ps(c6,x2);
  ret_avx=_mm256_mul_ps(_mm256_add_ps(c7,ret_avx),x2);
  ret_avx=_mm256_mul_ps(_mm256_add_ps(c8,ret_avx),x2);
    x_avx=_mm256_xor_ps(x_avx,sign);
  ret_avx=_mm256_mul_ps(_mm256_add_ps(c9,ret_avx),x2);
  ret_avx=_mm256_add_ps(x_avx,_mm256_mul_ps(ret_avx,x_avx));
  float ret[8];
  
  _mm256_storeu_ps(ret,ret_avx);
    return ret[0];
}

inline void sin_f32_avx_intrin_stream(const float *xs,float *ys,int n)
{
    if(n<=0) return;
    int n8=n/8;
    int nd=n-8*n8;
    for(int i=0;i<n8;++i)
    {
        __m256 x=_mm256_loadu_ps(xs+8*i);
        __m256 y=_mm256_mul_ps(c0,x);
        y=_mm256_add_ps(y,c1);
        __m256 sign=y;
        y=_mm256_sub_ps(y,c1);
        x=_mm256_sub_ps(x,_mm256_mul_ps(c2,y));
        x=_mm256_sub_ps(x,_mm256_mul_ps(c3,y));
        sign=_mm256_castsi256_ps(_mm256_slli_epi32(_mm256_castps_si256(sign),31));
        x=_mm256_sub_ps(x,_mm256_mul_ps(c4,y));
        x=_mm256_sub_ps(x,_mm256_mul_ps(c5,y));
        x=_mm256_xor_ps(x,sign);
        _mm256_storeu_ps(ys+8*i,x);
    }
    for(int i=0;i<n8;++i)
    {
        __m256 x=_mm256_loadu_ps(ys+8*i);
        __m256 x2=_mm256_mul_ps(x,x);
        __m256 y=_mm256_mul_ps(c6,x2);
        y=_mm256_mul_ps(_mm256_add_ps(c7,y),x2);
        y=_mm256_mul_ps(_mm256_add_ps(c8,y),x2);
      y=_mm256_mul_ps(_mm256_add_ps(c9,y),x2);
        y=_mm256_add_ps(x,_mm256_mul_ps(y,x));
        _mm256_storeu_ps(ys+8*i,y);
    }
    for(int i=0;i<nd;++i)
        ys[8*n8+i]=sin_f32_avx_intrin(xs[8*n8+i]);
}

int main()
{
    static const int M=1024;
    static const int N=1000*1024;
    static float x[M];
    static float y[M];
    for(int i=0;i<M;++i)
        x[i]=float(i);
    double t=double(clock());
    for(int i=0;i<N;++i)
    {
        sin_f32_avx_intrin_stream(x,y,M);
    }
    t=double(clock())-t;
    float b=0.0f;
    for(int i=0;i<M;++i)
    {
        if(fabs(sinf(x[i])-y[i])>b) b=fabs(sinf(x[i])-y[i]);
    }
    printf("Max. abs. error: %.7g\n",b);
    printf("%.1f ns/value\n",1.0e+9*t/(double(N)*double(M)*double(CLOCKS_PER_SEC)));
    return 0;
}

Но и меня его и проверить-то негде.

Оно вообще работает?

FordPerfect
Постоялец

#85
0:40, 2 сен 2015

Да, цивилизованные люди, видимо, использовали бы FMA, но речь не об этом, пока.

FordPerfect
Постоялец

#86
19:24, 5 сен 2015

FordPerfect
>Та как-то не зажгла.
>Вопрос: так и должно быть, или у меня руки кривые?
Понял, я хреново мерил.
Если убрать

x[1]=s;

в test_sse_x4, то будет примерно та же скорость, что и у одноствольной.

И да, есть желающие потестить код из #88? Я не знаю даже, запускается ли он.