# Faster computation of $ke^{-(x - h)^2}$

The question is quite simple; almost every computer language today provides the $$\exp(x)$$ function in their standard library to compute expressions like $$ke^{-(x - h)^2}.$$ However, I would like to know whether this function is the fastest way to compute the above expression. In other words, is there some way to compute $$ke^{-(x - h)^2}$$ faster than $$\exp(x)$$ in standard libraries while keeping the result very accurate?

I would like to specify that Taylor series will not work for my application, nor will any other polynomial approximations.

In the end, I wish to compute a series of these expressions, $$\sum^n_{i = 0} k_ie^{-(x - h_i)^2}$$ for some integral $$n \geq 0.$$

• "nor will any other polynomial approximations" Eh what?
– orlp
Jul 25 '20 at 20:45
• The way $\exp$ is often computed is doing Pade' approximations. These are rational approximations.
– plop
Jul 26 '20 at 3:23
• @Venkataram Sivaram Does the hardware platform you work on support fused multiply-add (FMA) operation, and is this functionality exposed in the programming language you program in? Jul 27 '20 at 4:40
• @Venkataram Sivaram Is the computation performed in single precision or double precision? Jul 27 '20 at 7:03
• Please don't use "Edit:". Instead, revise your question to what it should have been from the start, so it reads well for someone who encounters it for the first time. See cs.meta.stackexchange.com/q/657/755.
– D.W.
Jul 28 '20 at 1:19

There are several reasons why your standard library implementation of exp is slower than you might hope:

• It should to be accurate to 0.5ulp (or 1ulp at worst), no matter what number you give it. This means that special care may need to be taken, say, in the subnormal range.
• It needs to work for all floating point numbers, and return a reasonable answer for all of them. Assuming that you're on an IEEE-754 platform, this includes NaN, Inf, etc. The standard explains what the "correct" answer is for all of those cases, and the library function needs to get it right.
• It may not be specialised for your specific instruction set (i.e. your CPU may have some advanced floating-point instructions that libc doesn't use by default).

This is what it means to be general-purpose.

Chances are that your standard library has had a lot of effort put into it, and you cannot beat it and still be general-purpose. To get more performance, you need to relax some restrictions. So we're going to need some details. How much accuracy is enough? Can you assume anything about the range of inputs?

I've gone through a typical implementation elsewhere (see this previous answer), and that should help you get an idea of how these things can be implemented.

• I would like to keep the same amount of precision that you mention, but suppose that the range of $h$ and $k$ are both $[0, 1000].$ Would this allow computation to be made slightly faster? Also, the application I wish to use this for involves computing the sum of several terms of the specified form. Is there a way to simply this sum? Jul 27 '20 at 1:28
• Probably not on the first point, unless your standard library only supplies a double-precision exp and you are working in single-precision. On the second point, one option is to use vector instructions to compute multiple values at the same time. Although anything table-driven is going to have problems there because of the dependent loads. Jul 27 '20 at 2:57
• @VenkataramSivaram, please edit the question to state your requirements -- don't leave them in comments here. Thank you!
– D.W.
Jul 27 '20 at 3:00
• Now that I think of it, you can probably compute exp(-x^2) more efficiently than using the built-in exp. Let me think about it. This might be a good question for Computational Science. Jul 27 '20 at 4:15
• Thank you, I have posted this question in the Computational Science stack (scicomp.stackexchange.com/questions/35634/…). Perhaps my wording in this question will make my intentions a bit clearer... Jul 27 '20 at 16:52

TL;DR: Using the exp function of your library is likely the fastest way to compute $$ke^{-(x-h)^{2}}$$, but it's usually worthwhile to do some limited experiments to double-check that assumption.

Generally speaking, when programming for a reasonably mature platform, standard math libraries tend to be highly optimized by both domain experts and microarchitecture experts. While not all functions may be optimized equally well, the most commonly used math functions (exp, log, sin, cos, atan{2}) tend to be the most heavily optimized.

I am assuming you have already profiled your code to establish that the multiple calls to exp are a bottleneck in your code, and that you have double-checked your algorithm(s) to minimize the calls to this function. I further assume that you have already established that you cannot perform the computation at lower precision (say, float instead of double), which results in a significant performance increase on most platforms.

Are you using the latest compiler and libraries available for your platform? Performance improvements are incorporated all the time, so recent tool chains with their associated libraries tend to offer the highest performance. Are you targeting the compiler's code generation to the architecture that most closely reflects your processor's architecture? Newer processors tend to add performance enhancing hardware, such as fused multiply-add (FMA) units and wider SIMD operations and compiler often need to be instructed to use them via compiler flags, e.g. -march=core-avx2.

Also, make sure you are maxing out compiler optimizations. Some advanced optimizations may require adding compiler switches by hand as they are not subsumed under -O3. Examples could be auto-vectorization, whole-program optimization (by use of an optimizing linker), or profile-guided optimizations. Your math library may offer multiple levels of performance / accuracy trade-offs. For example, Intel's MKL provides three modes: high accuracy (maximum error < 1 ulp), lower accuracy (maximum eror < 4 ulp), enhanced performance. The lower the accuracy requirement, the higher the performance.

Note that the overall numerical error in the evaluation of the expression will very likely be dominated by the error in the exp argument magnified through exponentiation. Depending on the magnitude of the argument, a 1 ulp error in the input can well turn into a 1000 ulp error in the output. In light of that, the exp function itself does not have to be extremely accurate.

Standard math library functions need to follow the relevant language specification exactly, which includes overhead for the handling of special cases and detection of errors. Standards may also mandate certain accuracy requirements. If your use case allows the elimination of special case handling and a reduction in accuracy, you could try to roll your own function, like the exemplary C implementation below, which requires hardware support for FMA. It is usually a good idea to use the tool chain specific attributes to force inlining of any custom function to eliminate function call overhead and improve instruction scheduling flexibility.

#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>

double uint64_as_double (uint64_t a)
{
double r;
memcpy (&r, &a, sizeof r);
return r;
}

uint64_t double_as_uint64 (double a)
{
uint64_t r;
memcpy (&r, &a, sizeof r);
return r;
}

/* Compute exponential function e**x. Maximum error found in testing: < 0.9 ulp */
double my_exp (double a)
{
const double ln2_hi = 6.9314718055829871e-01;
const double ln2_lo = 1.6465949582897082e-12;
const double l2e = 1.4426950408889634; // log2(e)
const double cvt = 6755399441055744.0; // 3 * 2**51
double f, j, p, r;
uint64_t i;

// exp(a) = exp2(i) * exp(f); i = rint (a / log(2))
j = fma (l2e, a, cvt);
i = double_as_uint64 (j);
j = j - cvt;
f = fma (j, -ln2_hi, a);
f = fma (j, -ln2_lo, f);
// approximate p = exp(f) on interval [-log(2)/2, +log(2)/2]
p = fma (p, f, 2.7630903491116071e-7); // 0x1.28af3fcaa8f70p-22
p = fma (p, f, 2.7557514543681978e-6); // 0x1.71dee62382584p-19
p = fma (p, f, 2.4801491039342422e-5); // 0x1.a01997c8b03e6p-16
p = fma (p, f, 1.9841269589067952e-4); // 0x1.a01a01475dae0p-13
p = fma (p, f, 1.3888888945916467e-3); // 0x1.6c16c1852b7d7p-10
p = fma (p, f, 8.3333333334557717e-3); // 0x1.11111111224c6p-7
p = fma (p, f, 4.1666666666519782e-2); // 0x1.55555555502a5p-5
p = fma (p, f, 1.6666666666666477e-1); // 0x1.5555555555511p-3
p = fma (p, f, 5.0000000000000122e-1); // 0x1.000000000000bp-1
p = fma (p, f, 1.0000000000000000e+0); // 0x1.0000000000000p+0
p = fma (p, f, 1.0000000000000000e+0); // 0x1.0000000000000p+0

// exp(a) = 2**i * exp(f);
uint64_t ri = (double_as_uint64 (p) + (i << 52));
r = uint64_as_double (ri);

// handle special cases
double fa = fabs (a);
if (! (fa < 708.0)) { // |a| >= 708 requires double scaling
i = (a > 0.0) ? 0ULL : 0x8030000000000000ULL;
r = uint64_as_double (0x7fe0000000000000ULL + i);
r = r * uint64_as_double (ri - i - 0x3ff0000000000000ULL);
if (! (fa < 746.0)) { // |a| >= 746 severe overflow / underflow
r = (a > 0.0) ? INFINITY : 0.0;
if (isnan (a)) {
r = a + a;
}
}
}
return r;
}