Which algorithm is used to compute exponential functions the GNU C++ Standard Library?

Question

Please consider std::exp defined in the header cmath in C++ numerics library. Now, please consider an implementation of the C++ Standard Library, say libstdc++.

Considering there are various algorithms to compute the elementary functions, such as arithmetic-geometric mean iteration algorithm to compute the exponential function and three others shown here;

Could you please name the particular algorithm being used to compute the exponential function in libstdc++, if possible?

PS: I could not pinpoint either the correct tarballs containing the std::exp implementation or comprehend the relevant file contents, I'm afraid.

Hoad

libstdc++ just forwards to the compiler intrinsic __builtin_exp, which will have a varying implementation by platform and compiler. – Hippocrates

__builtin_exp is just a fancy name for ::exp, the C function declared in math.h. Gcc does not contain any implementation of exp, you want to look for it in your C library (e.g. glibc). – Enviable

dug up the sources for you. See the updated answer – Mcnalley

Answer 1

It doesn't use any intricate algorithm at all. Note that std::exp is only defined for a very limited number of types: float, double and long double + any Integral type that is castable to double. That makes it not necessary to implement complicated maths.

Currently, it uses the builtin __builtin_expf as can be verified from the source code. This compiles to a call to expf on my machine which is a call into libm coming from glibc. Let's see what we find in their source code. When we search for expf we find that this internally calls __ieee754_expf which is a system-dependant implementation. Both i686 and x86_64 just include a glibc/sysdeps/ieee754/flt-32/e_expf.c which finally gives us an implementation (reduced for brevity, the look into the sources

It is basically a order 3 polynomial approximation for floats:

static inline uint32_t
top12 (float x)
{
  return asuint (x) >> 20;
}

float
__expf (float x)
{
  uint64_t ki, t;
  /* double_t for better performance on targets with FLT_EVAL_METHOD==2.  */
  double_t kd, xd, z, r, r2, y, s;

  xd = (double_t) x;
  // [...] skipping fast under/overflow handling

  /* x*N/Ln2 = k + r with r in [-1/2, 1/2] and int k.  */
  z = InvLn2N * xd;

  /* Round and convert z to int, the result is in [-150*N, 128*N] and
     ideally ties-to-even rule is used, otherwise the magnitude of r
     can be bigger which gives larger approximation error.  */
  kd = roundtoint (z);
  ki = converttoint (z);
  r = z - kd;

  /* exp(x) = 2^(k/N) * 2^(r/N) ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */
  t = T[ki % N];
  t += ki << (52 - EXP2F_TABLE_BITS);
  s = asdouble (t);
  z = C[0] * r + C[1];
  r2 = r * r;
  y = C[2] * r + 1;
  y = z * r2 + y;
  y = y * s;
  return (float) y;
}

Similarly, for 128-bit long double, it's an order 7 approximation and for double they use more complicated algorithm that I can't make sense of right now.

Mcnalley

You are looking at x87 instructions, which are only remotely relevant for long double nowadays. If you look at asm, __builtin_expf becomes a call to expf (provided by the C library). – Enviable

@MarcGlisse I've tried getting gcc to inline that call but didn't succeed. – Mcnalley

Of course it didn't inline it. How often have you seen functions from external libraries inlined in your code? – Enviable

I found this answer really helpful! Oddly enough, it looks like this is not the actual algorithm which is used in current releases of glibc. The source for version 2.26 (the most recent release, from Feb 1 2018) is different - see sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/flt-32/… . The source was recently modified (sourceware.org/git/…) but this change seems to have not made it through to release – Armful

I got my timelines wrong in the previous comment. Glibc version 2.26 is from Aug 2017 and has the old algorithm (different from that posted in this answer) but version 2.27 (from Feb 1 2018) uses the algorithm presented above. Still, if your glibc isn't super up-to-date, then the algorithm your machine uses is probably different. – Armful

the mathematical formula on which this is based is the taylor series of exp(x) – Docile