compile visual studio c++ online

Language: Layout:

// 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 <emmintrin.h>

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

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

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

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

Show compiler warnings [ + ] Compiler args [ + ] Show input
Compilation time: 0.94 sec, absolute running time: 3 sec, absolute service time: 3.97 sec