mitsuba-af602c6fd98a/src/libcore/ssemath.cpp
/* natural logarithm computed for 4 simultaneous float
   return NaN for x <= 0
*/
natural logarithm computed for 4 simultaneous float
return NaN for x <= 0
natural logarithm computed for 4 simultaneous float return NaN for x < = 0
__m128 log_ps(__m128 x) {
    typedef __m128 v4sf;
    typedef __m128i v4si;
    v4si emm0;
    v4sf one = constants::ps_1.ps;
    v4sf invalid_mask = _mm_cmple_ps(x, _mm_setzero_ps());
    x = _mm_max_ps(x, constants::min_norm_pos.ps);  
    emm0 = _mm_srli_epi32(_mm_castps_si128(x), 23);
    x = _mm_and_ps(x, constants::inv_mant_mask.ps);
    x = _mm_or_ps(x,  constants::ps_0p5.ps);
    emm0 = _mm_sub_epi32(emm0, constants::pi32_0x7f.pi);
    v4sf e = _mm_cvtepi32_ps(emm0);
    e = _mm_add_ps(e, one);
    v4sf mask = _mm_cmplt_ps(x, constants::cephes_SQRTHF.ps);
    v4sf tmp = _mm_and_ps(x, mask);
    x = _mm_sub_ps(x, one);
    e = _mm_sub_ps(e, _mm_and_ps(one, mask));
    x = _mm_add_ps(x, tmp);
    v4sf z = _mm_mul_ps(x,x);
    v4sf y = constants::cephes_log_p0.ps;
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, constants::cephes_log_p1.ps);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, constants::cephes_log_p2.ps);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, constants::cephes_log_p3.ps);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, constants::cephes_log_p4.ps);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, constants::cephes_log_p5.ps);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, constants::cephes_log_p6.ps);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, constants::cephes_log_p7.ps);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, constants::cephes_log_p8.ps);
    y = _mm_mul_ps(y, x);
    y = _mm_mul_ps(y, z);
    tmp = _mm_mul_ps(e, constants::cephes_log_q1.ps);
    y = _mm_add_ps(y, tmp);
    tmp = _mm_mul_ps(z, constants::ps_0p5.ps);
    y = _mm_sub_ps(y, tmp);
    tmp = _mm_mul_ps(e, constants::cephes_log_q2.ps);
    x = _mm_add_ps(x, y);
    x = _mm_add_ps(x, tmp);
    x = _mm_or_ps(x, invalid_mask); 
    return x;
}
__m128 log_ps ( __m128 x ) { typedef __m128 v4sf ; typedef __m128i v4si ; v4si emm0 ; v4sf one = constants :: ps_1 . ps ; v4sf invalid_mask = _mm_cmple_ps ( x , _mm_setzero_ps ( ) ) ; x = _mm_max_ps ( x , constants :: min_norm_pos . ps ) ; emm0 = _mm_srli_epi32 ( _mm_castps_si128 ( x ) , 23 ) ; x = _mm_and_ps ( x , constants :: inv_mant_mask . ps ) ; x = _mm_or_ps ( x , constants :: ps_0p5 . ps ) ; emm0 = _mm_sub_epi32 ( emm0 , constants :: pi32_0x7f . pi ) ; v4sf e = _mm_cvtepi32_ps ( emm0 ) ; e = _mm_add_ps ( e , one ) ; v4sf mask = _mm_cmplt_ps ( x , constants :: cephes_SQRTHF . ps ) ; v4sf tmp = _mm_and_ps ( x , mask ) ; x = _mm_sub_ps ( x , one ) ; e = _mm_sub_ps ( e , _mm_and_ps ( one , mask ) ) ; x = _mm_add_ps ( x , tmp ) ; v4sf z = _mm_mul_ps ( x , x ) ; v4sf y = constants :: cephes_log_p0 . ps ; y = _mm_mul_ps ( y , x ) ; y = _mm_add_ps ( y , constants :: cephes_log_p1 . ps ) ; y = _mm_mul_ps ( y , x ) ; y = _mm_add_ps ( y , constants :: cephes_log_p2 . ps ) ; y = _mm_mul_ps ( y , x ) ; y = _mm_add_ps ( y , constants :: cephes_log_p3 . ps ) ; y = _mm_mul_ps ( y , x ) ; y = _mm_add_ps ( y , constants :: cephes_log_p4 . ps ) ; y = _mm_mul_ps ( y , x ) ; y = _mm_add_ps ( y , constants :: cephes_log_p5 . ps ) ; y = _mm_mul_ps ( y , x ) ; y = _mm_add_ps ( y , constants :: cephes_log_p6 . ps ) ; y = _mm_mul_ps ( y , x ) ; y = _mm_add_ps ( y , constants :: cephes_log_p7 . ps ) ; y = _mm_mul_ps ( y , x ) ; y = _mm_add_ps ( y , constants :: cephes_log_p8 . ps ) ; y = _mm_mul_ps ( y , x ) ; y = _mm_mul_ps ( y , z ) ; tmp = _mm_mul_ps ( e , constants :: cephes_log_q1 . ps ) ; y = _mm_add_ps ( y , tmp ) ; tmp = _mm_mul_ps ( z , constants :: ps_0p5 . ps ) ; y = _mm_sub_ps ( y , tmp ) ; tmp = _mm_mul_ps ( e , constants :: cephes_log_q2 . ps ) ; x = _mm_add_ps ( x , y ) ; x = _mm_add_ps ( x , tmp ) ; x = _mm_or_ps ( x , invalid_mask ) ; return x ; }
natural logarithm computed for 4 simultaneous float
return NaN for x <= 0
natural logarithm computed for 4 simultaneous float return NaN for x < = 0
// cut off denormalized stuff // keep only the fractional part /* part2:
       if( x < SQRTHF ) {
         e -= 1;
         x = x + x - 1.0;
       } else { x = x - 1.0; }
    */ // negative arg will be NAN