19

I'm looking for an approximation of the natural exponential function operating on SSE element. Namely - __m128 exp( __m128 x ).

I have an implementation which is quick but seems to be very low in accuracy:

static inline __m128 FastExpSse(__m128 x)
{
    __m128 a = _mm_set1_ps(12102203.2f); // (1 << 23) / ln(2)
    __m128i b = _mm_set1_epi32(127 * (1 << 23) - 486411);
    __m128  m87 = _mm_set1_ps(-87);
    // fast exponential function, x should be in [-87, 87]
    __m128 mask = _mm_cmpge_ps(x, m87);

    __m128i tmp = _mm_add_epi32(_mm_cvtps_epi32(_mm_mul_ps(a, x)), b);
    return _mm_and_ps(_mm_castsi128_ps(tmp), mask);
}

Could anybody have an implementation with better accuracy yet as fast (Or faster)?

I'd be happy if it is written in C Style.

Thank You.

7
  • 1
    What is the accuracy of your current implementation (maximum relative error)? And what accuracy are you targeting for an improved version? Commented Oct 30, 2017 at 22:55
  • It is really bad (I'd even suspect there is an error, if someone could spot it). So probably anything which makes sense will beat it. Something below 1% would be great! Commented Oct 30, 2017 at 22:56
  • See: sse_mathfun and this answer (relates to log2 but most of the suggestions also include an exp function). Commented Oct 31, 2017 at 7:41
  • @PaulR, How different it is from njuffa's solution? Commented Oct 31, 2017 at 9:16
  • @Royi: there are several different options in the linked answer, but in general they are probably somewhat more accurate and somewhat slower. So it depends on whereabouts on the accuracy-versus-performance curve you want to be. Commented Oct 31, 2017 at 9:23

6 Answers 6

35

The C code below is a translation into SSE intrinsics of an algorithm I used in a previous answer to a similar question.

The basic idea is to transform the computation of the standard exponential function into computation of a power of 2: expf (x) = exp2f (x / logf (2.0f)) = exp2f (x * 1.44269504). We split t = x * 1.44269504 into an integer i and a fraction f, such that t = i + f and 0 <= f <= 1. We can now compute 2f with a polynomial approximation, then scale the result by 2i by adding i to the exponent field of the single-precision floating-point result.

One problem that exists with an SSE implementation is that we want to compute i = floorf (t), but there is no fast way to compute the floor() function. However, we observe that for positive numbers, floor(x) == trunc(x), and that for negative numbers, floor(x) == trunc(x) - 1, except when x is a negative integer. However, since the core approximation can handle an f value of 1.0f, using the approximation for negative arguments is harmless. SSE provides an instruction to convert single-precision floating point operands to integers with truncation, so this solution is efficient.

Peter Cordes points out that SSE4.1 supports a fast floor function _mm_floor_ps(), so a variant using SSE4.1 is also shown below. Not all toolchains automatically predefine the macro __SSE4_1__ when SSE 4.1 code generation is enabled, but gcc does.

Compiler Explorer (Godbolt) shows that gcc 7.2 compiles the code below into sixteen instructions for plain SSE and twelve instructions for SSE 4.1.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <emmintrin.h>
#ifdef __SSE4_1__
#include <smmintrin.h>
#endif

/* max. rel. error = 1.72863156e-3 on [-87.33654, 88.72283] */
__m128 fast_exp_sse (__m128 x)
{
    __m128 t, f, e, p, r;
    __m128i i, j;
    __m128 l2e = _mm_set1_ps (1.442695041f);  /* log2(e) */
    __m128 c0  = _mm_set1_ps (0.3371894346f);
    __m128 c1  = _mm_set1_ps (0.657636276f);
    __m128 c2  = _mm_set1_ps (1.00172476f);

    /* exp(x) = 2^i * 2^f; i = floor (log2(e) * x), 0 <= f <= 1 */   
    t = _mm_mul_ps (x, l2e);             /* t = log2(e) * x */
#ifdef __SSE4_1__
    e = _mm_floor_ps (t);                /* floor(t) */
    i = _mm_cvtps_epi32 (e);             /* (int)floor(t) */
#else /* __SSE4_1__*/
    i = _mm_cvttps_epi32 (t);            /* i = (int)t */
    j = _mm_srli_epi32 (_mm_castps_si128 (x), 31); /* signbit(t) */
    i = _mm_sub_epi32 (i, j);            /* (int)t - signbit(t) */
    e = _mm_cvtepi32_ps (i);             /* floor(t) ~= (int)t - signbit(t) */
#endif /* __SSE4_1__*/
    f = _mm_sub_ps (t, e);               /* f = t - floor(t) */
    p = c0;                              /* c0 */
    p = _mm_mul_ps (p, f);               /* c0 * f */
    p = _mm_add_ps (p, c1);              /* c0 * f + c1 */
    p = _mm_mul_ps (p, f);               /* (c0 * f + c1) * f */
    p = _mm_add_ps (p, c2);              /* p = (c0 * f + c1) * f + c2 ~= 2^f */
    j = _mm_slli_epi32 (i, 23);          /* i << 23 */
    r = _mm_castsi128_ps (_mm_add_epi32 (j, _mm_castps_si128 (p))); /* r = p * 2^i*/
    return r;
}

int main (void)
{
    union {
        float f[4];
        unsigned int i[4];
    } arg, res;
    double relerr, maxrelerr = 0.0;
    int i, j;
    __m128 x, y;

    float start[2] = {-0.0f, 0.0f};
    float finish[2] = {-87.33654f, 88.72283f};

    for (i = 0; i < 2; i++) {

        arg.f[0] = start[i];
        arg.i[1] = arg.i[0] + 1;
        arg.i[2] = arg.i[0] + 2;
        arg.i[3] = arg.i[0] + 3;
        do {
            memcpy (&x, &arg, sizeof(x));
            y = fast_exp_sse (x);
            memcpy (&res, &y, sizeof(y));
            for (j = 0; j < 4; j++) {
                double ref = exp ((double)arg.f[j]);
                relerr = fabs ((res.f[j] - ref) / ref);
                if (relerr > maxrelerr) {
                    printf ("arg=% 15.8e  res=%15.8e  ref=%15.8e  err=%15.8e\n", 
                            arg.f[j], res.f[j], ref, relerr);
                    maxrelerr = relerr;
                }
            }   
            arg.i[0] += 4;
            arg.i[1] += 4;
            arg.i[2] += 4;
            arg.i[3] += 4;
        } while (fabsf (arg.f[3]) < fabsf (finish[i]));
    }
    printf ("maximum relative errror = %15.8e\n", maxrelerr);
    return EXIT_SUCCESS;
}

An alternative design for fast_sse_exp() extracts the integer portion of the adjusted argument x / log(2) in round-to-nearest mode, using the well-known technique of adding the "magic" conversion constant 1.5 * 223 to force rounding in the correct bit position, then subtracting out the same number again. This requires that the SSE rounding mode in effect during the addition is "round to nearest or even", which is the default. wim pointed out in comments that some compilers may optimize out the addition and subtraction of the conversion constant cvt as redundant when aggressive optimization is used, interfering with the functionality of this code sequence, so it is recommended to inspect the machine code generated. The approximation interval for computation of 2f is now centered around zero, since -0.5 <= f <= 0.5, requiring a different core approximation.

/* max. rel. error <= 1.72860465e-3 on [-87.33654, 88.72283] */
__m128 fast_exp_sse (__m128 x)
{
    __m128 t, f, p, r;
    __m128i i, j;

    const __m128 l2e = _mm_set1_ps (1.442695041f); /* log2(e) */
    const __m128 cvt = _mm_set1_ps (12582912.0f);  /* 1.5 * (1 << 23) */
    const __m128 c0 =  _mm_set1_ps (0.238428936f);
    const __m128 c1 =  _mm_set1_ps (0.703448006f);
    const __m128 c2 =  _mm_set1_ps (1.000443142f);

    /* exp(x) = 2^i * 2^f; i = rint (log2(e) * x), -0.5 <= f <= 0.5 */
    t = _mm_mul_ps (x, l2e);             /* t = log2(e) * x */
    r = _mm_sub_ps (_mm_add_ps (t, cvt), cvt); /* r = rint (t) */
    f = _mm_sub_ps (t, r);               /* f = t - rint (t) */
    i = _mm_cvtps_epi32 (t);             /* i = (int)t */
    p = c0;                              /* c0 */
    p = _mm_mul_ps (p, f);               /* c0 * f */
    p = _mm_add_ps (p, c1);              /* c0 * f + c1 */
    p = _mm_mul_ps (p, f);               /* (c0 * f + c1) * f */
    p = _mm_add_ps (p, c2);              /* p = (c0 * f + c1) * f + c2 ~= exp2(f) */
    j = _mm_slli_epi32 (i, 23);          /* i << 23 */
    r = _mm_castsi128_ps (_mm_add_epi32 (j, _mm_castps_si128 (p))); /* r = p * 2^i*/
    return r;
}

The algorithm for the code in the question appears to be taken from the work of Nicol N. Schraudolph, which cleverly exploits the semi-logarithmic nature of IEEE-754 binary floating-point formats:

N. N. Schraudolph. "A fast, compact approximation of the exponential function." Neural Computation, 11(4), May 1999, pp.853-862.

After removal of the argument clamping code, it reduces to just three SSE instructions. The "magical" correction constant 486411 is not optimal for minimizing maximum relative error over the entire input domain. Based on simple binary search, the value 298765 seems to be superior, reducing maximum relative error for FastExpSse() to 3.56e-2 vs. maximum relative error of 1.73e-3 for fast_exp_sse().

/* max. rel. error = 3.55959567e-2 on [-87.33654, 88.72283] */
__m128 FastExpSse (__m128 x)
{
    __m128 a = _mm_set1_ps (12102203.0f); /* (1 << 23) / log(2) */
    __m128i b = _mm_set1_epi32 (127 * (1 << 23) - 298765);
    __m128i t = _mm_add_epi32 (_mm_cvtps_epi32 (_mm_mul_ps (a, x)), b);
    return _mm_castsi128_ps (t);
}

Schraudolph's algorithm basically uses the linear approximation 2f ~= 1.0 + f for f in [0,1], and its accuracy could be improved by adding a quadratic term. The clever part of Schraudolph's approach is computing 2i * 2f without explicitly separating the integer portion i = floor(x * 1.44269504) from the fraction. I see no way to extend that trick to a quadratic approximation, but one can certainly combine the floor() computation from Schraudolph with the quadratic approximation used above:

/* max. rel. error <= 1.72886892e-3 on [-87.33654, 88.72283] */
__m128 fast_exp_sse (__m128 x)
{
    __m128 f, p, r;
    __m128i t, j;
    const __m128 a = _mm_set1_ps (12102203.0f); /* (1 << 23) / log(2) */
    const __m128i m = _mm_set1_epi32 (0xff800000); /* mask for integer bits */
    const __m128 ttm23 = _mm_set1_ps (1.1920929e-7f); /* exp2(-23) */
    const __m128 c0 = _mm_set1_ps (0.3371894346f);
    const __m128 c1 = _mm_set1_ps (0.657636276f);
    const __m128 c2 = _mm_set1_ps (1.00172476f);

    t = _mm_cvtps_epi32 (_mm_mul_ps (a, x));
    j = _mm_and_si128 (t, m);            /* j = (int)(floor (x/log(2))) << 23 */
    t = _mm_sub_epi32 (t, j);
    f = _mm_mul_ps (ttm23, _mm_cvtepi32_ps (t)); /* f = (x/log(2)) - floor (x/log(2)) */
    p = c0;                              /* c0 */
    p = _mm_mul_ps (p, f);               /* c0 * f */
    p = _mm_add_ps (p, c1);              /* c0 * f + c1 */
    p = _mm_mul_ps (p, f);               /* (c0 * f + c1) * f */
    p = _mm_add_ps (p, c2);              /* p = (c0 * f + c1) * f + c2 ~= 2^f */
    r = _mm_castsi128_ps (_mm_add_epi32 (j, _mm_castps_si128 (p))); /* r = p * 2^i*/
    return r;
}
Sign up to request clarification or add additional context in comments.

22 Comments

@Royi Note that I posted my test framework along with fast_exp_sse() iself, so you should be able to verify my accuracy claim and test your existing function as well, for a side-by-side comparison.
For code-style, I'd highly recommend _mm_set1_ps() for the constant vectors in the first code block as well. Initializing a __m128 with a braced initializer is not even guaranteed to be portable, but I think it does work on gcc/clang and I assume MSVC. But repeating each constant 4 times in the source is not nice.
SSE4.1 roundps provides a fast _mm_floor_ps (), so you could #ifdef __SSE4_1__ that part to use a more efficient floor when it's enable at compile time. (MSVC doesn't define that macro, so you need other checks there...)
@PeterCordes, could you post your code which is optimized for SSE4? Thank You.
@PeterCordes Note that my code requires floor(t) both as an integer i (to be added to the exponent field later) and a floating-point number e (for the computation of the "fraction" f). So that would give us: e = _mm_floor_ps (t); i = _mm_cvtps_epi32 (e);.
|
17

A good increase in accuracy in my algorithm (implementation FastExpSse in the answer above) can be obtained at the cost of an integer subtraction and floating-point division by using FastExpSse(x/2)/FastExpSse(-x/2) instead of FastExpSse(x). The trick here is to set the shift parameter (298765 above) to zero so that the piecewise linear approximations in the numerator and denominator line up to give you substantial error cancellation. Roll it into a single function:

__m128 BetterFastExpSse (__m128 x)
{
  const __m128 a = _mm_set1_ps ((1 << 22) / float(M_LN2));  // to get exp(x/2)
  const __m128i b = _mm_set1_epi32 (127 * (1 << 23));       // NB: zero shift!
  __m128i r = _mm_cvtps_epi32 (_mm_mul_ps (a, x));
  __m128i s = _mm_add_epi32 (b, r);
  __m128i t = _mm_sub_epi32 (b, r);
  return _mm_div_ps (_mm_castsi128_ps (s), _mm_castsi128_ps (t));
}

(I'm not a hardware guy - how bad a performance killer is the division here?)

If you need exp(x) just to get y = tanh(x) (e.g. for neural networks), use FastExpSse with zero shift as follows:

a = FastExpSse(x);
b = FastExpSse(-x);
y = (a - b)/(a + b);

to get the same type of error cancellation benefit. The logistic function works similarly, using FastExpSse(x/2)/(FastExpSse(x/2) + FastExpSse(-x/2)) with zero shift. (This is just to show the principle - you obviously don't want to evaluate FastExpSse multiple times here, but roll it into a single function along the lines of BetterFastExpSse above.)

I did develop a series of higher-order approximations from this, ever more accurate but also slower. Unpublished but happy to collaborate if anyone wants to give them a spin.

And finally, for some fun: use in reverse gear to get FastLogSse. Chaining that with FastExpSse gives you both operator and error cancellation, and out pops a blazingly fast power function...

8 Comments

One division mixed in with a lot of multiplies is the best case. divps has worse latency than a multiply, but starting a division doesn't block the multiply / FMA units, so if you're doing multiple exp() calls on independent data, it's going to be about as cheap as another multiply. (And thus much better than using on large polynomial or double-precision or w/e.) The main cost here is from evaluating the same polynomial twice, not the division at the end, on modern x86.
@technosaurus I've removed auto. This is generic code, I wasn't meaning to imply any particular language.
@Peter no need to evaluate a polynomial twice - since the argument is merely negated, it's just an extra integer subtraction. I've expanded the code snippet to clarify. Good news about the division!
Oh, I hadn't realized how simple the rest of the code was. This probably will bottleneck on divps throughput on most CPUs, especially a 256-bit AVX version, if you're doing just this over an array. But to feed another calculation it might still be fine. @njuffa's 2nd version fast_exp_ssewith relative error 1.72886892e-3 might be faster, especially on Haswell (where FMA is available for very fast polylnomials, but divps isn't as fast as Skylake).
Nice. For completeness: After changing the last line of BetterFastExpSse() to return _mm_div_ps (_mm_castsi128_ps (s), _mm_castsi128_ps (t)); I observed a maximum relative error of 1.049e-2 for arguments in [-87.33654, 88.72283].
|
8

Going back through my notes from way back then, I did explore ways to improve the accuracy without using division. I used the same reinterpret-as-float trick but applied a polynomial correction to the mantissa which was essentially calculated in 16-bit fixed-point arithmetic (the only way to do it fast back then).

The cubic resp. quartic versions give you 4 resp. 5 significant digits of accuracy. There was no point increasing the order beyond that, as the noise of the low-precision arithmetic then starts to drown out the error of the polynomial approximation. Here are the plain C versions:

#include <stdint.h>

float fastExp3(register float x)  // cubic spline approximation
{
    union { float f; int32_t i; } reinterpreter;

    reinterpreter.i = (int32_t)(12102203.0f*x) + 127*(1 << 23);
    int32_t m = (reinterpreter.i >> 7) & 0xFFFF;  // copy mantissa
    // empirical values for small maximum relative error (8.34e-5):
    reinterpreter.i +=
         ((((((((1277*m) >> 14) + 14825)*m) >> 14) - 79749)*m) >> 11) - 626;
    return reinterpreter.f;
}

float fastExp4(register float x)  // quartic spline approximation
{
    union { float f; int32_t i; } reinterpreter;

    reinterpreter.i = (int32_t)(12102203.0f*x) + 127*(1 << 23);
    int32_t m = (reinterpreter.i >> 7) & 0xFFFF;  // copy mantissa
    // empirical values for small maximum relative error (1.21e-5):
    reinterpreter.i += (((((((((((3537*m) >> 16)
        + 13668)*m) >> 18) + 15817)*m) >> 14) - 80470)*m) >> 11);
    return reinterpreter.f;
}

The quartic one obeys (fastExp4(0f) == 1f), which can be important for fixed-point iteration algorithms.

How efficient are these integer multiply-shift-add sequences in SSE? On architectures where float arithmetic is just as fast, one could use that instead, reducing the arithmetic noise. This would essentially yield cubic and quartic extensions of @njuffa's answer above.

9 Comments

M_LN2 is a double, so the resulting asm has to convert float to double. With (float)((1 << 23)/M_LN2), gcc7.3 auto-vectorizes these pretty reasonably in a loop (godbolt.org/g/vEWofG), with SSE4.1 or AVX2 (vcvttps2dq packed FP->int conversion, vpmulld packed 32-bit multiply / vpsrad packed arithmetic right shift), but clang doesn't autovectorize at all; it still converts to scalar. So you'd need to manually vectorize for this to be good with clang, I guess.
I haven't looked at the speed vs. precision tradeoff or this vs. others, but vpmulld has one per 2 clock throughput on Haswell (1/4 the throughput of FP mul or FMA). It's still 2 uops on Skylake, but SKL can run it on p0 / p1 for a throughput of one per clock. But it takes up more execution resources. This might be good on Nehalem or maybe Sandybridge where 128-bit pmulld is 1 uop, and FP division is more expensive. Static analysis of these with IACA (What is IACA and how do I use it?) would be interesting and not difficult
@Royi sorry I don't know enough SSE/AVX for that.
After replacing (1 << 23)/M_LN2 with 12102203.0f, an exhaustive test of all IEEE-754 binary32 floating-point numbers in [-87.33654, 88.72283] shows: max. rel. error fastExp3() = 8.34e-5; max. rel. error fastExp4() = 1.21e-5
@njuffa I suspect the higher errors you find are due to the inaccurate constant - can you try the exact value 12102203.161561f? I've edited my answer to use that, in order to avoid the problems with (1 << 23)/M_LN2 that you and @Peter describe. Cheers
|
2

There is a paper about creating fast versions of these equations (tanh, cosh, artanh, sinh, etc):

http://ijeais.org/wp-content/uploads/2018/07/IJAER180702.pdf "Creating a Compiler Optimized Inlineable Implementation of Intel Svml Simd Intrinsics"

their tanh equation 6, on page 9 is very similar to @NicSchraudolph answer

1 Comment

Normally link-only answers aren't welcome on SO. This arguably belongs as a comment, but I'm not going to flag it because it's potentially useful enough.
2

I have developed for my purposes the following function that calculates quickly and accurately the natural exponent with single precision. The function works in the entire range of float values. The code is written under Visual Studio (x86). AVX is used instead of SSE, but that shouldn't be a problem. The accuracy of this function is almost the same as standard expf function, but significantly faster. Used approximation is based on the Chebyshev series expansion of the function f(t)=t/(2^(t/2)-1)+t/2 for t from the [-1; 1]. I thank Peter Cordes for his good advice.

_declspec(naked) float _vectorcall fexp(float x)
{
  static const float ct[7] =       // Constants table
  {
    1.44269502f,                   // lb(e)
    1.92596299E-8f,                // Correction to the value lb(e)
    -9.21120925E-4f,               // 16*b2
    0.115524396f,                  // 4*b1
    2.88539004f,                   // b0
    2.0f,                          // 2
    4.65661287E-10f                // 2^-31
  };
  _asm
  {
    mov ecx,offset ct              // ecx contains the address of constants tables
    vmulss xmm1,xmm0,[ecx]         // xmm1 = x*lb(e)
    vcvtss2si eax,xmm1             // eax = round(x*lb(e)) = k
    cdq                            // edx=-1, if x<0 or overflow, otherwise edx=0
    vmovss xmm3,[ecx+8]            // Initialize the sum with highest coefficient 16*b2
    and edx,4                      // edx=4, if x<0 or overflow, otherwise edx=0
    vcvtsi2ss xmm1,xmm1,eax        // xmm1 = k
    lea eax,[eax+8*edx]            // Add 32 to exponent, if x<0
    vfmsub231ss xmm1,xmm0,[ecx]    // xmm1 = x*lb(e)-k = t/2 in the range from -0,5 to 0,5
    add eax,126                    // The exponent of 2^(k-1) or 2^(k+31) with bias 127
    jle exp_low                    // Jump if x<<0 or overflow (|x| too large or x=NaN)
    vfmadd132ss xmm0,xmm1,[ecx+4]  // xmm0 = t/2 (corrected value)
    cmp eax,254                    // Check that the exponent is not too large
    jg exp_inf                     // Jump to set Inf if overflow
    vmulss xmm2,xmm0,xmm0          // xmm2 = t^2/4 - the argument of the polynomial
    shl eax,23                     // The bits of the float value 2^(k-1) or 2^(k+31)
    vfmadd213ss xmm3,xmm2,[ecx+12] // xmm3 = 4*b1+4*b2*t^2
    vmovd xmm1,eax                 // xmm1 = 2^(k-1) или 2^(k+31)
    vfmsub213ss xmm3,xmm2,xmm0     // xmm3 = -t/2+b1*t^2+b2*t^4
    vaddss xmm0,xmm0,xmm0          // xmm0 = t
    vaddss xmm3,xmm3,[ecx+16]      // xmm3 = b0-t/2+b1*t^2+b2*t^4 = f(t)-t/2
    vdivss xmm0,xmm0,xmm3          // xmm0 = t/(f(t)-t/2)
    vfmadd213ss xmm0,xmm1,xmm1     // xmm0 = e^x with shifted exponent of -1 or 31
    vmulss xmm0,xmm0,[ecx+edx+20]  // xmm0 = e^x
    ret                            // Return
      exp_low:                     // Handling the case of x<<0 or overflow
    vucomiss xmm0,[ecx]            // Check the sign of x and a condition x=NaN
    jp exp_end                     // Complete with NaN result, if x=NaN
      exp_inf:                     // Entry point for processing large x
    vxorps xmm0,xmm0,xmm0          // xmm0 = 0
    jc exp_end                     // Ready, if x<<0
    vrcpss xmm0,xmm0,xmm0          // xmm0 = Inf in case x>>0
      exp_end:                     // The result at xmm0 is ready
    ret                            // Return
  }
}

Below I post a simplified algorithm. Support for denormalized numbers in the result is removed here.

_declspec(naked) float _vectorcall fexp(float x)
{
  static const float ct[5] =       // Constants table
  {
    1.44269502f,                   // lb(e)
    1.92596299E-8f,                // Correction to the value lb(e)
    -9.21120925E-4f,               // 16*b2
    0.115524396f,                  // 4*b1
    2.88539004f                    // b0
  };
  _asm
  {
    mov edx,offset ct              // edx contains the address of constants tables
    vmulss xmm1,xmm0,[edx]         // xmm1 = x*lb(e)
    vcvtss2si eax,xmm1             // eax = round(x*lb(e)) = k
    vmovss xmm3,[edx+8]            // Initialize the sum with highest coefficient 16*b2
    vcvtsi2ss xmm1,xmm1,eax        // xmm1 = k
    cmp eax,127                    // Check that the exponent is not too large
    jg exp_break                   // Jump to set Inf if overflow
    vfmsub231ss xmm1,xmm0,[edx]    // xmm1 = x*lb(e)-k = t/2 in the range from -0,5 to 0,5
    add eax,127                    // Receive the exponent of 2^k with the bias 127
    jle exp_break                  // The result is 0, if x<<0
    vfmadd132ss xmm0,xmm1,[edx+4]  // xmm0 = t/2 (corrected value)
    vmulss xmm2,xmm0,xmm0          // xmm2 = t^2/4 - the argument of polynomial
    shl eax,23                     // eax contains the bits of 2^k
    vfmadd213ss xmm3,xmm2,[edx+12] // xmm3 = 4*b1+4*b2*t^2
    vmovd xmm1,eax                 // xmm1 = 2^k
    vfmsub213ss xmm3,xmm2,xmm0     // xmm3 = -t/2+b1*t^2+b2*t^4
    vaddss xmm0,xmm0,xmm0          // xmm0 = t
    vaddss xmm3,xmm3,[edx+16]      // xmm3 = b0-t/2+b1*t^2+b2*t^4 = f(t)-t/2
    vdivss xmm0,xmm0,xmm3          // xmm0 = t/(f(t)-t/2)
    vfmadd213ss xmm0,xmm1,xmm1     // xmm0 = 2^k*(t/(f(t)-t/2)+1) = e^x
    ret                            // Return
      exp_break:                   // Get 0 for x<0 or Inf for x>>0
    vucomiss xmm0,[edx]            // Check the sign of x and a condition x=NaN
    jp exp_end                     // Complete with NaN result, if x=NaN
    vxorps xmm0,xmm0,xmm0          // xmm0 = 0
    jc exp_end                     // Ready, if x<<0
    vrcpss xmm0,xmm0,xmm0          // xmm0 = Inf, if x>>0
      exp_end:                     // The result at xmm0 is ready
    ret                            // Return
  }
}

5 Comments

Some interesting ideas here. Looks like some room left for some optimization: Convert to/from integer is usually not faster than vroundss to implement nearbyint. Especially on AMD CPUs where it's 1 uop, 3 cycle latency. You aren't using the integer result for anything except a large-magnitude check, which you could do on the FP bit-pattern you're later getting with vmovd.
e.g. add eax,eax to shift out the sign bit, then compare the rest of the bit pattern for being above the bit-pattern of the largest magnitude input you handle. cmp eax, 0x4f000000 << 1 / ja abs_gt_2_31. or 0x4effffff for 2^31-1. Note that NaN and Inf have bit-patterns where this jumps, potentially needing only one (hopefully not-taken) branch on the fast path to catch all the weird cases and get them sorted out. If the normal path through your code produces +Inf or zero correctly for big magnitude positive or negative exponents up to some range, you could use that as the limit.
h-schmidt.net/FloatConverter/IEEE754.html is handy for getting the hex bit-patterns since you can't use std::bit_cast<> as an inline asm constant (I don't think), and MSVC inline assembly probably doesn't have NASM __float32__(1<<31). The bit-pattern for +128 is 0x43000000 (so that's just above the log2 of FLT_MAX, 2^127 x 1.999..., and similarly log2(FLT_MIN) is just about -128, although subnormal floats go down to 2^-149. Flushing would-be subnormal results to zero is a valid option, though.)
If you have code that sometimes produces exp(n) = Infinity, for n that's a bit too large, probably best to let that case continue into the normal path of execution, instead of causing a branch miss though, for values where that code still works, instead of going through the special case return 0 or Inf path. That looks like what you're doing, using float->int conversion overflow as the cutoff. Oh, but it looks like you later cmp eax,254/ja after unbiasing, so maybe we can delay some checking until then. (If we assume the special cases are rare, it's fine to not leave the main path ASAP)
Added a simplified and optimized version into response.
2

For softmax use, I'm envisioning the flow as:

auto a = _mm_mul_ps(x, _mm_set1_ps(12102203.2f));
auto b = _mm_castsi128_ps(_mm_cvtps_epi32(a)); // so far as in other variants

// copy 9 MSB from 0x3f800000 over 'b' so that 1 <= c < 2
//  - also  1 <= poly_eval(...) < 2
auto c = replace_exponent(b, _mm_set1_ps(1.0f));
auto d = poly_eval(c, kA, kB, kC);  // 2nd degree polynomial
auto e = replace_exponent(d, b);    // restore exponent : 2^i * 2^f

The exponent copying can be done as bitwise select using a proper mask (AVX-512 has vpternlogd, and I'm using actually Arm Neon vbsl).

All the input values x must be negative and clamped between -17-f(N) <= x <= -f(N), so that when scaled by (1<<23)/log(2), the maximum sum of the N resulting floating point values do not reach infinity and that the reciprocal does not become denormal. For N=3, f(N) = 4. Larger f(N) will trade off input precision.

The polyeval coefficients are generated for example by polyfit([1 1.5 2],[1 sqrt(2) 2]), with kA=0.343146, kB=-0.029437, kC=0.68292, producing strictly values smaller than 2 and preventing discontinuities. The maximum average error can be diminished by evaluating the polynomial at x=[1+max_err 1.5-eps 2], y=[1 2^(.5-eps) 2-max_err].

For strictly SSE/AVX, exponent replacement for 1.0f can be done by (x & 0x007fffff) | 0x3f800000). A two instruction sequence for the latter exponent replacement can be found by ensuring that poly_eval(x) evaluates to a range, which can be directly ored with b & 0xff800000.


Revisiting the problem, I found a superior method for 2nd degree polynomial compensation: the equation for 2 ^ x, 0 <= x <= 1 can be approximated with (1+x)(1+x)/3 + 2/3 with the largest relative error of 3.4130e-03.

The non-zero initial estimate of the exponential is already in the format of e' = 2^k(1+x). By simply two bit operations (assuming that we were computing negative values -- a bit more if we are computing generic exponentiations), we can

  1. remove the fractional bits: i = e' & 0x3f800000
  2. convert to format (1+x): o = e' | 0x3f800000

Then we blend: 3.0f * e'' = e' * o + i * 2.0f

  • if we are going to use the values for softmax, it doesn't matter that we have an extra scaling by 3.0f. If it bothers, we just compute i * 0.3333333f, o * 0.6666666f instead.
  // initial estimation -- see bottom of page for alternative
  auto a = _mm_mul_ps(x, _mm_set1_ps(12102203.2f));
  auto b = _mm_castsi128_ps(_mm_cvtps_epi32(a));

  // 2nd degree improvement
  auto i = _mm_and_ps(b, _mm_set1_ps(1.0f));
  auto o = _mm_or_ps(b, _mm_set1_ps(1.0f));
  auto i2 = _mm_add_ps(i, i);
  b = _mm_mul_ps(b, o);
  b = _mm_add_ps(b, i2); // OR FMA with b += 2*i

Additionally I found a more suitable conversion for the first problem (of generating the initial approximation b) in ARM64 architecture, which often does not instruction level parallelise the int<->float conversions:

Here one does

   auto a = _mm_add_ps(x, _mm_set1_ps(639.0f));
   a = _mm_castsi128_ps(_mm_shli_epi32(_mm_castps_si128(a), 9));

The maximum relative error in this method is ~6%, however just like Schraudolph's method the argument can be shifted to slightly better range to have mean relative error close to 0. The magic constant would be closer to 638.95f then.

In Intel with 3 different formats for MLA, this form of `638.95 + x * scale_log_2_log_e` would be a single instruction; that followed by integer shift by 9 would then make the exp approximation computable in 2 instructions.

enter image description here

(Here's the relative error approx(2^x) / 2^x for -3 < x < 3)

See also https://stackoverflow.com/a/78509807/1716339

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.