source/projects/vb.mi.rppls_tilde/vcvrack/include/simd/sse_mathfun.h
inline __m128 sse_mathfun_sin_ps(__m128 x) { // any x
	__m128 xmm1, xmm2 = _mm_setzero_ps(), xmm3, sign_bit, y;

	__m128i emm0, emm2;
	sign_bit = x;
	/* take the absolute value */
	const __m128 inv_sign_mask = _mm_castsi128_ps(_mm_set1_epi32(~0x80000000));
	x = _mm_and_ps(x, inv_sign_mask);
	/* extract the sign bit (upper one) */
	const __m128 sign_mask = _mm_castsi128_ps(_mm_set1_epi32(0x80000000));
	sign_bit = _mm_and_ps(sign_bit, sign_mask);

	/* scale by 4/Pi */
	const __m128 cephes_FOPI = _mm_set_ps1(1.27323954473516); // 4 / M_PI
	y = _mm_mul_ps(x, cephes_FOPI);

	/* store the integer part of y in mm0 */
	emm2 = _mm_cvttps_epi32(y);
	/* j=(j+1) & (~1) (see the cephes sources) */
	emm2 = _mm_add_epi32(emm2, _mm_set1_epi32(1));
	emm2 = _mm_and_si128(emm2, _mm_set1_epi32(~1));
	y = _mm_cvtepi32_ps(emm2);

	/* get the swap sign flag */
	emm0 = _mm_and_si128(emm2, _mm_set1_epi32(4));
	emm0 = _mm_slli_epi32(emm0, 29);
	/* get the polynom selection mask
	   there is one polynom for 0 <= x <= Pi/4
	   and another one for Pi/4<x<=Pi/2

	   Both branches will be computed.
	*/
	emm2 = _mm_and_si128(emm2, _mm_set1_epi32(2));
	emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());

	__m128 swap_sign_bit = _mm_castsi128_ps(emm0);
	__m128 poly_mask = _mm_castsi128_ps(emm2);
	sign_bit = _mm_xor_ps(sign_bit, swap_sign_bit);

	/* The magic pass: "Extended precision modular arithmetic"
	   x = ((x - y * DP1) - y * DP2) - y * DP3; */
	xmm1 = _mm_set_ps1(-0.78515625);
	xmm2 = _mm_set_ps1(-2.4187564849853515625e-4);
	xmm3 = _mm_set_ps1(-3.77489497744594108e-8);
	xmm1 = _mm_mul_ps(y, xmm1);
	xmm2 = _mm_mul_ps(y, xmm2);
	xmm3 = _mm_mul_ps(y, xmm3);
	x = _mm_add_ps(x, xmm1);
	x = _mm_add_ps(x, xmm2);
	x = _mm_add_ps(x, xmm3);

	/* Evaluate the first polynom  (0 <= x <= Pi/4) */
	y = _mm_set_ps1(2.443315711809948E-005);
	__m128 z = _mm_mul_ps(x, x);

	y = _mm_mul_ps(y, z);
	y = _mm_add_ps(y, _mm_set_ps1(-1.388731625493765E-003));
	y = _mm_mul_ps(y, z);
	y = _mm_add_ps(y, _mm_set_ps1(4.166664568298827E-002));
	y = _mm_mul_ps(y, z);
	y = _mm_mul_ps(y, z);
	__m128 tmp = _mm_mul_ps(z, _mm_set_ps1(0.5));
	y = _mm_sub_ps(y, tmp);
	y = _mm_add_ps(y, _mm_set_ps1(1.0));

	/* Evaluate the second polynom  (Pi/4 <= x <= 0) */

	__m128 y2 = _mm_set_ps1(-1.9515295891E-4);
	y2 = _mm_mul_ps(y2, z);
	y2 = _mm_add_ps(y2, _mm_set_ps1(8.3321608736E-3));
	y2 = _mm_mul_ps(y2, z);
	y2 = _mm_add_ps(y2, _mm_set_ps1(-1.6666654611E-1));
	y2 = _mm_mul_ps(y2, z);
	y2 = _mm_mul_ps(y2, x);
	y2 = _mm_add_ps(y2, x);

	/* select the correct result from the two polynoms */
	xmm3 = poly_mask;
	y2 = _mm_and_ps(xmm3, y2); //, xmm3);
	y = _mm_andnot_ps(xmm3, y);
	y = _mm_add_ps(y, y2);
	/* update the sign */
	y = _mm_xor_ps(y, sign_bit);
	return y;
}
/* evaluation of 4 sines at onces, using only SSE1+MMX intrinsics so
   it runs also on old athlons XPs and the pentium III of your grand
   mother.

   The code is the exact rewriting of the cephes sinf function.
   Precision is excellent as long as x < 8192 (I did not bother to
   take into account the special handling they have for greater values
   -- it does not return garbage for arguments over 8192, though, but
   the extra precision is missing).

   Note that it is such that sinf((float)M_PI) = 8.74e-8, which is the
   surprising but correct result.

   Performance is also surprisingly good, 1.33 times faster than the
   macos vsinf SSE2 function, and 1.5 times faster than the
   __vrs4_sinf of amd's ACML (which is only available in 64 bits). Not
   too bad for an SSE1 function (with no special tuning) !
   However the latter libraries probably have a much better handling of NaN,
   Inf, denormalized and other special arguments..

   On my core 1 duo, the execution of this function takes approximately 95 cycles.

   From what I have observed on the experiments with Intel AMath lib, switching to an
   SSE2 version would improve the perf by only 10%.

   Since it is based on SSE intrinsics, it has to be compiled at -O2 to
   deliver full speed.
*/
evaluation of 4 sines at onces, using only SSE1+MMX intrinsics so
it runs also on old athlons XPs and the pentium III of your grand
mother.

The code is the exact rewriting of the cephes sinf function.
Precision is excellent as long as x < 8192 (I did not bother to
take into account the special handling they have for greater values
it does not return garbage for arguments over 8192, though, but
the extra precision is missing).

Note that it is such that sinf((float)M_PI) = 8.74e-8, which is the
surprising but correct result.

Performance is also surprisingly good, 1.33 times faster than the
macos vsinf SSE2 function, and 1.5 times faster than the
vrs4_sinf of amd's ACML (which is only available in 64 bits). Not
too bad for an SSE1 function (with no special tuning) .
However the latter libraries probably have a much better handling of NaN,
Inf, denormalized and other special arguments

On my core 1 duo, the execution of this function takes approximately 95 cycles.

From what I have observed on the experiments with Intel AMath lib, switching to an
SSE2 version would improve the perf by only 10%.

Since it is based on SSE intrinsics, it has to be compiled at -O2 to
deliver full speed.
evaluation of 4 sines at onces using only SSE1 + MMX intrinsics so it runs also on old athlons XPs and the pentium III of your grand mother . The code is the exact rewriting of the cephes sinf function . Precision is excellent as long as x < 8192 ( I did not bother to take into account the special handling they have for greater values it does not return garbage for arguments over 8192 though but the extra precision is missing ) . Note that it is such that sinf (( float ) M_PI ) = 8 . 74e - 8 which is the surprising but correct result . Performance is also surprisingly good 1 . 33 times faster than the macos vsinf SSE2 function and 1 . 5 times faster than the vrs4_sinf of amd ' s ACML ( which is only available in 64 bits ) . Not too bad for an SSE1 function ( with no special tuning ) . However the latter libraries probably have a much better handling of NaN Inf denormalized and other special arguments On my core 1 duo the execution of this function takes approximately 95 cycles . From what I have observed on the experiments with Intel AMath lib switching to an SSE2 version would improve the perf by only 10% . Since it is based on SSE intrinsics it has to be compiled at - O2 to deliver full speed .
inline __m128 sse_mathfun_sin_ps(__m128 x) { 
	__m128 xmm1, xmm2 = _mm_setzero_ps(), xmm3, sign_bit, y;
	__m128i emm0, emm2;
	sign_bit = x;
	const __m128 inv_sign_mask = _mm_castsi128_ps(_mm_set1_epi32(~0x80000000));
	x = _mm_and_ps(x, inv_sign_mask);
	const __m128 sign_mask = _mm_castsi128_ps(_mm_set1_epi32(0x80000000));
	sign_bit = _mm_and_ps(sign_bit, sign_mask);
	const __m128 cephes_FOPI = _mm_set_ps1(1.27323954473516); 
	y = _mm_mul_ps(x, cephes_FOPI);
	emm2 = _mm_cvttps_epi32(y);
	emm2 = _mm_add_epi32(emm2, _mm_set1_epi32(1));
	emm2 = _mm_and_si128(emm2, _mm_set1_epi32(~1));
	y = _mm_cvtepi32_ps(emm2);
	emm0 = _mm_and_si128(emm2, _mm_set1_epi32(4));
	emm0 = _mm_slli_epi32(emm0, 29);
	emm2 = _mm_and_si128(emm2, _mm_set1_epi32(2));
	emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
	__m128 swap_sign_bit = _mm_castsi128_ps(emm0);
	__m128 poly_mask = _mm_castsi128_ps(emm2);
	sign_bit = _mm_xor_ps(sign_bit, swap_sign_bit);
	xmm1 = _mm_set_ps1(-0.78515625);
	xmm2 = _mm_set_ps1(-2.4187564849853515625e-4);
	xmm3 = _mm_set_ps1(-3.77489497744594108e-8);
	xmm1 = _mm_mul_ps(y, xmm1);
	xmm2 = _mm_mul_ps(y, xmm2);
	xmm3 = _mm_mul_ps(y, xmm3);
	x = _mm_add_ps(x, xmm1);
	x = _mm_add_ps(x, xmm2);
	x = _mm_add_ps(x, xmm3);
	y = _mm_set_ps1(2.443315711809948E-005);
	__m128 z = _mm_mul_ps(x, x);
	y = _mm_mul_ps(y, z);
	y = _mm_add_ps(y, _mm_set_ps1(-1.388731625493765E-003));
	y = _mm_mul_ps(y, z);
	y = _mm_add_ps(y, _mm_set_ps1(4.166664568298827E-002));
	y = _mm_mul_ps(y, z);
	y = _mm_mul_ps(y, z);
	__m128 tmp = _mm_mul_ps(z, _mm_set_ps1(0.5));
	y = _mm_sub_ps(y, tmp);
	y = _mm_add_ps(y, _mm_set_ps1(1.0));
	__m128 y2 = _mm_set_ps1(-1.9515295891E-4);
	y2 = _mm_mul_ps(y2, z);
	y2 = _mm_add_ps(y2, _mm_set_ps1(8.3321608736E-3));
	y2 = _mm_mul_ps(y2, z);
	y2 = _mm_add_ps(y2, _mm_set_ps1(-1.6666654611E-1));
	y2 = _mm_mul_ps(y2, z);
	y2 = _mm_mul_ps(y2, x);
	y2 = _mm_add_ps(y2, x);
	xmm3 = poly_mask;
	y2 = _mm_and_ps(xmm3, y2); 
	y = _mm_andnot_ps(xmm3, y);
	y = _mm_add_ps(y, y2);
	y = _mm_xor_ps(y, sign_bit);
	return y;
}
inline __m128 sse_mathfun_sin_ps ( __m128 x ) { __m128 xmm1 , xmm2 = _mm_setzero_ps ( ) , xmm3 , sign_bit , y ; __m128i emm0 , emm2 ; sign_bit = x ; const __m128 inv_sign_mask = _mm_castsi128_ps ( _mm_set1_epi32 ( ~ 0x80000000 ) ) ; x = _mm_and_ps ( x , inv_sign_mask ) ; const __m128 sign_mask = _mm_castsi128_ps ( _mm_set1_epi32 ( 0x80000000 ) ) ; sign_bit = _mm_and_ps ( sign_bit , sign_mask ) ; const __m128 cephes_FOPI = _mm_set_ps1 ( 1.27323954473516 ) ; y = _mm_mul_ps ( x , cephes_FOPI ) ; emm2 = _mm_cvttps_epi32 ( y ) ; emm2 = _mm_add_epi32 ( emm2 , _mm_set1_epi32 ( 1 ) ) ; emm2 = _mm_and_si128 ( emm2 , _mm_set1_epi32 ( ~ 1 ) ) ; y = _mm_cvtepi32_ps ( emm2 ) ; emm0 = _mm_and_si128 ( emm2 , _mm_set1_epi32 ( 4 ) ) ; emm0 = _mm_slli_epi32 ( emm0 , 29 ) ; emm2 = _mm_and_si128 ( emm2 , _mm_set1_epi32 ( 2 ) ) ; emm2 = _mm_cmpeq_epi32 ( emm2 , _mm_setzero_si128 ( ) ) ; __m128 swap_sign_bit = _mm_castsi128_ps ( emm0 ) ; __m128 poly_mask = _mm_castsi128_ps ( emm2 ) ; sign_bit = _mm_xor_ps ( sign_bit , swap_sign_bit ) ; xmm1 = _mm_set_ps1 ( -0.78515625 ) ; xmm2 = _mm_set_ps1 ( -2.4187564849853515625e-4 ) ; xmm3 = _mm_set_ps1 ( -3.77489497744594108e-8 ) ; xmm1 = _mm_mul_ps ( y , xmm1 ) ; xmm2 = _mm_mul_ps ( y , xmm2 ) ; xmm3 = _mm_mul_ps ( y , xmm3 ) ; x = _mm_add_ps ( x , xmm1 ) ; x = _mm_add_ps ( x , xmm2 ) ; x = _mm_add_ps ( x , xmm3 ) ; y = _mm_set_ps1 ( 2.443315711809948E-005 ) ; __m128 z = _mm_mul_ps ( x , x ) ; y = _mm_mul_ps ( y , z ) ; y = _mm_add_ps ( y , _mm_set_ps1 ( -1.388731625493765E-003 ) ) ; y = _mm_mul_ps ( y , z ) ; y = _mm_add_ps ( y , _mm_set_ps1 ( 4.166664568298827E-002 ) ) ; y = _mm_mul_ps ( y , z ) ; y = _mm_mul_ps ( y , z ) ; __m128 tmp = _mm_mul_ps ( z , _mm_set_ps1 ( 0.5 ) ) ; y = _mm_sub_ps ( y , tmp ) ; y = _mm_add_ps ( y , _mm_set_ps1 ( 1.0 ) ) ; __m128 y2 = _mm_set_ps1 ( -1.9515295891E-4 ) ; y2 = _mm_mul_ps ( y2 , z ) ; y2 = _mm_add_ps ( y2 , _mm_set_ps1 ( 8.3321608736E-3 ) ) ; y2 = _mm_mul_ps ( y2 , z ) ; y2 = _mm_add_ps ( y2 , _mm_set_ps1 ( -1.6666654611E-1 ) ) ; y2 = _mm_mul_ps ( y2 , z ) ; y2 = _mm_mul_ps ( y2 , x ) ; y2 = _mm_add_ps ( y2 , x ) ; xmm3 = poly_mask ; y2 = _mm_and_ps ( xmm3 , y2 ) ; y = _mm_andnot_ps ( xmm3 , y ) ; y = _mm_add_ps ( y , y2 ) ; y = _mm_xor_ps ( y , sign_bit ) ; return y ; }
evaluation of 4 sines at onces, using only SSE1+MMX intrinsics so
it runs also on old athlons XPs and the pentium III of your grand
mother.
evaluation of 4 sines at onces using only SSE1 + MMX intrinsics so it runs also on old athlons XPs and the pentium III of your grand mother .
// any x /* take the absolute value */ /* extract the sign bit (upper one) */ /* scale by 4/Pi */ // 4 / M_PI /* store the integer part of y in mm0 */ /* j=(j+1) & (~1) (see the cephes sources) */ /* get the swap sign flag */ /* get the polynom selection mask
	   there is one polynom for 0 <= x <= Pi/4
	   and another one for Pi/4<x<=Pi/2

	   Both branches will be computed.
	*/ /* The magic pass: "Extended precision modular arithmetic"
	   x = ((x - y * DP1) - y * DP2) - y * DP3; */ /* Evaluate the first polynom  (0 <= x <= Pi/4) */ /* Evaluate the second polynom  (Pi/4 <= x <= 0) */ /* select the correct result from the two polynoms */ //, xmm3); /* update the sign */