compile c++ gcc 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.

// Collection of sin/cos routines for float32.

// 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>

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

// Intrinsic version (scalar).

// |x|<65536:       rel. error < 1.9*10^-7 [ 2 ulp ]
// |x|<262144:      abs. error < 4*10^-6
// Garbage for |x|>8388608.
inline float sin_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);

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

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

// Intrinsic version (parallel).

// |x|<65536:       rel. error < 1.9*10^-7 [ 2 ulp ]
// |x|<262144:      abs. error < 4*10^-6
// Garbage for |x|>8388608.
inline __m128 sin_f32_sse2_x4_intrin(__m128 x)
{
	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);

__m128 y=_mm_mul_ps(c0,x);

// Round y to closest integer.
	y=_mm_add_ps(y,c1);
	__m128 sign=y;
	y=_mm_sub_ps(y,c1);

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

x=_mm_sub_ps(x,_mm_mul_ps(c2,y));
	x=_mm_sub_ps(x,_mm_mul_ps(c3,y));
	x=_mm_sub_ps(x,_mm_mul_ps(c4,y));
	x=_mm_sub_ps(x,_mm_mul_ps(c5,y));
	
    __m128 x2=_mm_mul_ps(x,x);

__m128 ret=_mm_mul_ps(c6,x2);
	ret=_mm_mul_ps(_mm_add_ps(c7,ret),x2);
	ret=_mm_mul_ps(_mm_add_ps(c8,ret),x2);
    x=_mm_xor_ps(x,sign);
	ret=_mm_mul_ps(_mm_add_ps(c9,ret),x2);
	ret=_mm_add_ps(x,_mm_mul_ps(ret,x));
    return ret;
}

// |x|<32768:       rel. error < 1.9*10^-7 [ 2 ulp ]
// |x|<131072:      abs. error < 1*10^-6
// Garbage for |x|>8388608.
inline __m128 cos_f32_sse2_x4_intrin(__m128 x)
{
	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);

__m128 y=_mm_mul_ps(c0,x);

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

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

x=_mm_sub_ps(x,_mm_mul_ps(c2,y));
	x=_mm_sub_ps(x,_mm_mul_ps(c3,y));
	x=_mm_sub_ps(x,_mm_mul_ps(c4,y));
	x=_mm_sub_ps(x,_mm_mul_ps(c5,y));
	
    __m128 x2=_mm_mul_ps(x,x);

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

// Fast imprecise scalar intrinsic version.

// |x|<50:          abs. error < 3*10^-6
// |x|<1000:        abs. error < 6*10^-5
// |x|<16000:       abs. error < 0.001
// |x|<1000000:     abs. error < 0.06
// Garbage for |x|>8388608.
inline float sin_f32_sse2_intrin_fast(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.14159265358979f);
	static const __m128 c7=_mm_set_ss(-1.8561107572168111e-4f);
	static const __m128 c8=_mm_set_ss(+8.3145750686526298e-3f);
	static const __m128 c9=_mm_set_ss(-1.6665840148925781e-1f);

__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));
	
    __m128 x2=_mm_mul_ss(x_sse,x_sse);

__m128 ret_sse=_mm_mul_ss(c7,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|<50:          abs. error < 3*10^-6
// |x|<1000:        abs. error < 6*10^-5
// |x|<8000:        abs. error < 0.001
// |x|<1000000:     abs. error < 0.06
// Garbage for |x|>8388608.    
inline float cos_f32_sse2_intrin_fast(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.14159265358979f);
	static const __m128 c7=_mm_set_ss(-1.8561107572168111e-4f);
	static const __m128 c8=_mm_set_ss(+8.3145750686526298e-3f);
	static const __m128 c9=_mm_set_ss(-1.6665840148925781e-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));
	
    __m128 x2=_mm_mul_ss(x_sse,x_sse);

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

// Fast imprecise parallel intrinsic version.

// |x|<50:          abs. error < 3*10^-6
// |x|<1000:        abs. error < 6*10^-5
// |x|<16000:       abs. error < 0.001
// |x|<1000000:     abs. error < 0.06
// Garbage for |x|>8388608.
inline __m128 sin_f32_sse2_x4_intrin_fast(__m128 x)
{
    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.14159265358979f);
	static const __m128 c7=_mm_set_ps1(-1.8561107572168111e-4f);
	static const __m128 c8=_mm_set_ps1(+8.3145750686526298e-3f);
	static const __m128 c9=_mm_set_ps1(-1.6665840148925781e-1f);

__m128 y=_mm_mul_ps(c0,x);

// Round y to closest integer.
	y=_mm_add_ps(y,c1);
	__m128 sign=y;
	y=_mm_sub_ps(y,c1);

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

x=_mm_sub_ps(x,_mm_mul_ps(c2,y));
	
    __m128 x2=_mm_mul_ps(x,x);

__m128 ret=_mm_mul_ps(c7,x2);
	ret=_mm_mul_ps(_mm_add_ps(c8,ret),x2);
    x=_mm_xor_ps(x,sign);
	ret=_mm_mul_ps(_mm_add_ps(c9,ret),x2);
	ret=_mm_add_ps(x,_mm_mul_ps(ret,x));
    return ret;
}

// |x|<50:          abs. error < 3*10^-6
// |x|<1000:        abs. error < 6*10^-5
// |x|<8000:        abs. error < 0.001
// |x|<1000000:     abs. error < 0.06
// Garbage for |x|>8388608.    
inline __m128 cos_f32_sse2_x4_intrin_fast(__m128 x)
{
    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.14159265358979f);
	static const __m128 c7=_mm_set_ps1(-1.8561107572168111e-4f);
	static const __m128 c8=_mm_set_ps1(+8.3145750686526298e-3f);
	static const __m128 c9=_mm_set_ps1(-1.6665840148925781e-1f);
    static const __m128 c10=_mm_set_ps1(0.5f);

__m128 y=_mm_mul_ps(c0,x);

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

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

x=_mm_sub_ps(x,_mm_mul_ps(c2,y));
	
    __m128 x2=_mm_mul_ps(x,x);

__m128 ret=_mm_mul_ps(c7,x2);
	ret=_mm_mul_ps(_mm_add_ps(c8,ret),x2);
    x=_mm_xor_ps(x,sign);
	ret=_mm_mul_ps(_mm_add_ps(c9,ret),x2);
	ret=_mm_add_ps(x,_mm_mul_ps(ret,x));
    return ret;
}

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

// Tests.

void test_range(
    float (*f)(float),
    double (*f_ref)(double),
    float a,float b,
    float step // 0 -> exhaustive.
)
{
    float x;
    float arg_abs=a;
    float err_abs=0.0;
    float arg_rel=a;
    float err_rel=0.0;
    uint32_t cnt=0;
    x=a;
    while(x<b)
    {
        double val=double(f(x));
        double ref=f_ref(double(x));
        double d=std::abs(val-ref);
        if((cnt&0xFFFFF)==0)
            printf("\rx=%.7g. %d values tested...                  ",x,cnt);
        ++cnt;
        if(x==a||d>err_abs)
        {
            arg_abs=x;
            err_abs=d;
        }
        if(d!=0.0) d/=std::abs(ref);
        if(x==a||d>err_rel)
        {
            arg_rel=x;
            err_rel=d;
        }
        if(step==0.0f) x=std::nextafter(x,b);
        else x+=step;
    }
    printf("\r[%.7g ; %.7g]:                 \n",a,b);
    printf("Max absolute error: %.7g (at %.7g).\n",err_abs,arg_abs);
    printf("Max relative error: %.7g (at %.7g).\n",err_rel,arg_rel);
}

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

void measure_speed(float (*f)(float))
{
    static const int M=1<<10;
    static const int N=1<<25;
    static const float D=1000.0f;
    static_assert(M>0&&(M&(M-1))==0,"M is not a power of 2.");
    static float x[M];
    volatile float sum=0.0;
    for(int i=0;i<M;++i)
        x[i]=D*float(sin(double(i)));
    uint64_t t=clock();
    for(int i=0;i<N;++i)
        sum+=f(x[i&(M-1)]);
    t=clock()-t;
    printf("~%.1f ns per call.\n",1.0e+9*(double(t)/double(N))/double(CLOCKS_PER_SEC));
}

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

static void test_sse2_x4(__m128 (*f)(__m128 x),double (*f_ref)(double))
{
    static const uint32_t N=1<<25;
    __m128 x=_mm_set_ps(0.0f,1.0f,500.0f,17000.0f);
    __m128 y;
    uint64_t t=clock();
    __m128 s;
    for(uint32_t i=0;i<N;++i)
    {
        y=f(x);
		s=_mm_add_ps(s,y);
    }
	float tmp0;
	_mm_store_ss(&tmp0,s);
	volatile float tmp1=tmp0;
    t=clock()-t;
    printf("~%.1f ns per call.\n",1.0e+9*(double(t)/double(N))/double(CLOCKS_PER_SEC));

printf("[ Abs. error: ");
    static float fx[4];
	static float fy[4];
    float e=0.0f;
    for(unsigned int i=0;i<1000;++i)
    {
        for(unsigned int j=0;j<4;++j)
            fx[j]=1.0f*sin(float(4*i+j));
		x=_mm_loadu_ps(fx);
        y=f(x);
		_mm_storeu_ps(fy,y);
        for(unsigned int j=0;j<4;++j)
            if(fabs(fy[j]-float(f_ref(fx[j])))>e)
                e=fabs(fy[j]-float(f_ref(fx[j])));

}
    printf("~%.7g ]\n",e);
}

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

int main()
{
    static const float pi=3.14159265358979323846f;

printf("========          Speed:          ========\n\n");
    printf("sinf: ");
    measure_speed(sinf);
	printf("cosf: ");
    measure_speed(cosf);
	printf("sin_f32_sse2_intrin: ");
    measure_speed(sin_f32_sse2_intrin);
	printf("cos_f32_sse2_intrin: ");
    measure_speed(cos_f32_sse2_intrin);
    printf("sin_f32_sse2_intrin_fast: ");
    measure_speed(sin_f32_sse2_intrin_fast);
    printf("cos_f32_sse2_intrin_fast: ");
    measure_speed(cos_f32_sse2_intrin_fast);
	printf("\nsin_f32_sse2_x4_intrin: \n");
	test_sse2_x4(sin_f32_sse2_x4_intrin,sin);
	printf("\ncos_f32_sse2_x4_intrin: \n");
	test_sse2_x4(cos_f32_sse2_x4_intrin,cos);
    printf("\nsin_f32_sse2_x4_intrin_fast: \n");
	test_sse2_x4(sin_f32_sse2_x4_intrin_fast,sin);
	printf("\ncos_f32_sse2_x4_intrin_fast: \n");
	test_sse2_x4(cos_f32_sse2_x4_intrin_fast,cos);
    printf("\n========        Precision:        ========\n");
	float (* const s)(float)=sin_f32_sse2_intrin_fast;
	double (* const s_ref)(double)=sin;
	float (* const c)(float)=cos_f32_sse2_intrin_fast;
	double (* const c_ref)(double)=cos;
	const float step=1.0f/1000.0f; // Relative step, 0.0f for exhaustive.
    printf("\nsin_f32_sse2_intrin_fast: \n");
    for(float d=pi/8.0f;d<1.0e+7f;d*=2.0f)
        test_range(s,s_ref,+d,+2.0f*d,d*step);
    printf("\ncos_f32_sse_intrin: \n");
    for(float d=pi/8.0f;d<1.0e+7f;d*=2.0f)
        test_range(c,c_ref,+d,+2.0f*d,d*step);
    return 0;
}

Show compiler warnings [ + ] Compiler args [ + ] Show input
Compilation time: 0.42 sec, absolute running time: 3.94 sec, cpu time: 3.8 sec, memory peak: 4 Mb, absolute service time: 4.37 sec