Function sseExp2

src/OpenColorIO/SSE.h:235–307

exp2 function in SSE version 2 The function exp2() is evaluated by performing argument reduction and then using Chebyshev polynomials to evaluate the function over a restricted range.

// reduction and then using Chebyshev polynomials to evaluate the function 
// over a restricted range.
inline __m128 sseExp2(__m128 x)
{
    // y = exp2( x ) = exp2(integer + fraction) 
    //               = exp2(integer) * exp2(fraction) 
    //               = zf * mexp

    // Compute the largest integer not greater than x, i.e., floor(x)
    // Note: cvttps_epi32 simply cast the float value to int. That means cvttps_epi32(-2.7) = -2
    // rather than -3, hence for negative numbers we need to add -1. This ensures that "fraction"
    // is always in the range [0, 1).  Note that _mm_castps_si128(0xFFFFFFFF) is -1.
    // If x is outside the INT_MIN to INT_MAX range, _mm_cvttps_epi32 will return 0x80000000
    // (i.e. INT_MIN, just the sign bit set), which Intel calls the "integer indefinite" value. 
    // When 1 is subtracted from INT_MIN, it gives INT_MAX.  So floor_x is wrong for values
    // outside [INT_MIN, INT_MAX] but it's ignored thanks to the checks at the bottom.
    // It's also wrong for x=NaN, but again it's ok since the polynomial returns NaN and
    // hence the output is NaN, regardless of floor_x.
    __m128i floor_x
        = _mm_add_epi32(                                // add a pair of integer arguments
            _mm_cvttps_epi32(x),                        // convert float to int via truncation
            _mm_castps_si128(                           // reinterpret cast float to int
                _mm_cmpnle_ps(EZERO, x)));              // NOT( EZERO <= x ) ? 0xFFFFFFFF : 0

    // Compute exp2(floor_x) by moving floor_x to the exponent bits of the floating-point number.
    __m128 zf
        = _mm_castsi128_ps(                             // reinterpret cast int to float
            _mm_slli_epi32(                             // left shift by EXP_SHIFT
                _mm_add_epi32(floor_x, EBIAS),          // add a pair of integer arguments
                EXP_SHIFT));

    __m128 iexp = _mm_cvtepi32_ps(floor_x);             // convert floor_x to float
    __m128 fraction = _mm_sub_ps(x, iexp);              // x - iexp

    // Compute exp2(fraction) using a polynomial approximation.
    __m128 mexp
        = _mm_add_ps(
            _mm_mul_ps(
                _mm_add_ps(
                    _mm_mul_ps(
                        _mm_add_ps(
                            _mm_mul_ps(
                                _mm_add_ps(
                                    _mm_mul_ps(PNEXP4, fraction),
                                    PNEXP3),
                                fraction),
                            PNEXP2),
                        fraction),
                    PNEXP1),
                fraction),
            PNEXP0);

    __m128 exp2 = _mm_mul_ps(zf, mexp);                 // zf * mexp

    // Handle underflow:
    // If the (unbiased) exponent of zf is less than -126, the result is smaller than
    // the smallest representable floating-point number and an underflow computation is
    // potentially happening. When this happens, force the result to zero.
    // Note that as described above, floor_x is inaccurate, so the test here uses x.
    exp2 = _mm_andnot_ps(                               // NOT(...) AND exp2