mitsuba-af602c6fd98a/src/libcore/ssemath.cpp
/* evaluation of 4 sines at onces, using only SSE2.

   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 SSE2.
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 SSE2 . 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 .
__m128 sin_ps(__m128 x) { 
    typedef __m128 v4sf;
    typedef __m128i v4si;
    v4sf xmm1, xmm2 = _mm_setzero_ps(), xmm3, sign_bit, y;
    v4si emm0, emm2;
    sign_bit = x;
    x = _mm_and_ps(x, constants::inv_mant_mask.ps);
    sign_bit = _mm_and_ps(sign_bit, constants::sign_mask.ps);
    y = _mm_mul_ps(x, constants::cephes_FOPI.ps);
    emm2 = _mm_cvttps_epi32(y);
    emm2 = _mm_add_epi32(emm2, constants::pi32_1.pi);
    emm2 = _mm_and_si128(emm2, constants::pi32_inv1.pi);
    y = _mm_cvtepi32_ps(emm2);
    emm0 = _mm_and_si128(emm2, constants::pi32_4.pi);
    emm0 = _mm_slli_epi32(emm0, 29);
    emm2 = _mm_and_si128(emm2, constants::pi32_2.pi);
    emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
    v4sf swap_sign_bit = _mm_castsi128_ps(emm0);
    v4sf poly_mask = _mm_castsi128_ps(emm2);
    sign_bit = _mm_xor_ps(sign_bit, swap_sign_bit);
    xmm1 = constants::minus_cephes_DP1.ps;
    xmm2 = constants::minus_cephes_DP2.ps;
    xmm3 = constants::minus_cephes_DP3.ps;
    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 = constants::coscof_p0.ps;
    v4sf z = _mm_mul_ps(x,x);
    y = _mm_mul_ps(y, z);
    y = _mm_add_ps(y, constants::coscof_p1.ps);
    y = _mm_mul_ps(y, z);
    y = _mm_add_ps(y, constants::coscof_p2.ps);
    y = _mm_mul_ps(y, z);
    y = _mm_mul_ps(y, z);
    v4sf tmp = _mm_mul_ps(z, constants::ps_0p5.ps);
    y = _mm_sub_ps(y, tmp);
    y = _mm_add_ps(y, constants::ps_1.ps);
    v4sf y2 = constants::sincof_p0.ps;
    y2 = _mm_mul_ps(y2, z);
    y2 = _mm_add_ps(y2, constants::sincof_p1.ps);
    y2 = _mm_mul_ps(y2, z);
    y2 = _mm_add_ps(y2, constants::sincof_p2.ps);
    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;
}
__m128 sin_ps ( __m128 x ) { typedef __m128 v4sf ; typedef __m128i v4si ; v4sf xmm1 , xmm2 = _mm_setzero_ps ( ) , xmm3 , sign_bit , y ; v4si emm0 , emm2 ; sign_bit = x ; x = _mm_and_ps ( x , constants :: inv_mant_mask . ps ) ; sign_bit = _mm_and_ps ( sign_bit , constants :: sign_mask . ps ) ; y = _mm_mul_ps ( x , constants :: cephes_FOPI . ps ) ; emm2 = _mm_cvttps_epi32 ( y ) ; emm2 = _mm_add_epi32 ( emm2 , constants :: pi32_1 . pi ) ; emm2 = _mm_and_si128 ( emm2 , constants :: pi32_inv1 . pi ) ; y = _mm_cvtepi32_ps ( emm2 ) ; emm0 = _mm_and_si128 ( emm2 , constants :: pi32_4 . pi ) ; emm0 = _mm_slli_epi32 ( emm0 , 29 ) ; emm2 = _mm_and_si128 ( emm2 , constants :: pi32_2 . pi ) ; emm2 = _mm_cmpeq_epi32 ( emm2 , _mm_setzero_si128 ( ) ) ; v4sf swap_sign_bit = _mm_castsi128_ps ( emm0 ) ; v4sf poly_mask = _mm_castsi128_ps ( emm2 ) ; sign_bit = _mm_xor_ps ( sign_bit , swap_sign_bit ) ; xmm1 = constants :: minus_cephes_DP1 . ps ; xmm2 = constants :: minus_cephes_DP2 . ps ; xmm3 = constants :: minus_cephes_DP3 . ps ; 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 = constants :: coscof_p0 . ps ; v4sf z = _mm_mul_ps ( x , x ) ; y = _mm_mul_ps ( y , z ) ; y = _mm_add_ps ( y , constants :: coscof_p1 . ps ) ; y = _mm_mul_ps ( y , z ) ; y = _mm_add_ps ( y , constants :: coscof_p2 . ps ) ; y = _mm_mul_ps ( y , z ) ; y = _mm_mul_ps ( y , z ) ; v4sf tmp = _mm_mul_ps ( z , constants :: ps_0p5 . ps ) ; y = _mm_sub_ps ( y , tmp ) ; y = _mm_add_ps ( y , constants :: ps_1 . ps ) ; v4sf y2 = constants :: sincof_p0 . ps ; y2 = _mm_mul_ps ( y2 , z ) ; y2 = _mm_add_ps ( y2 , constants :: sincof_p1 . ps ) ; y2 = _mm_mul_ps ( y2 , z ) ; y2 = _mm_add_ps ( y2 , constants :: sincof_p2 . ps ) ; 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 SSE2.
evaluation of 4 sines at onces using only SSE2 .
// any x /* take the absolute value */ /* extract the sign bit (upper one) */ /* scale by 4/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 */