# A clean way to compute exp(-1/x^2) near x=0 in C?

I am looking for a clean way to compute preciseley, when $x$ is very close to zero: $$\exp(-1/x^2)$$ using C. What is the best way (speed, precision, etc.) ?

• You probably already know this, but this function is not analytic at x=0 and its Taylor expansion there is just identically f(x)=0. Therefore the common trick of using a low order polynomial approximation is not helpful. – k20 Dec 12 '13 at 18:34
• @k20 Indeed.... – mellow Dec 12 '13 at 18:52
• @mellow do you expect a problem by just implementing it straightforward in C using IEEE arithmetic? – GertVdE Dec 12 '13 at 20:23
• @mellow what are you doing with the result? Taken on its own, you're unlikely to end up with useful precision 'near 0', but the scale may be able to be cancelled in part of a larger expression. – Max Hutchinson Dec 12 '13 at 21:09
• It's not a sensible thing to calculate using floating point arithmetic, since if x<.03 (assuming 11 bit exponent) you'll underflow: en.wikipedia.org/wiki/Arithmetic_underflow – p.s. Dec 14 '13 at 9:42

Use $\exp\left(-\tfrac1{x^2}\right) =\frac1{\exp\left(\tfrac1{x^2}\right)}$ and expand in the denominator, for example as

$\exp\left(-\tfrac1{x^2}\right) \approx \frac1{1+1/x^2+1/(2x^4)}=\frac{4x^4}{(1+2x^2)^2+1}$

which is correct for $1/x^6$ below machine precision, i.e. for $x>1000$ and has the qualitatively correct behavior for smaller x.

You can also use the halving and squaring technique to get the cross-over from numerically to qualitatively exact closer to zero,

$\exp\left(-\tfrac1{x^2}\right)=\exp\left(-\tfrac1{nx^2}\right)^n \approx \left[\frac1{1+1/(nx^2)+1/(2n^2x^4)}\right]^{n} = \left[\frac{4n^2x^4}{1+(1+2nx^2)^2}\right]^{n}$

which is numerically exact for $1/(nx)^6$ below machine precision, $nx^2>10^5$, for example with $n=10$ for $x>100$, and qualitatively exact below that.

For the numerator in the last example to give a value above machine precision, $400x^4>10^{-2}$ or $x>0.02$ is necessary. Since this approximation is an upper bound for the true value, this is as close as it gets, i.e., one could alternatively return $0$ for $x<0.02$ and $\exp(-1/x^2)$ for $x\ge 0.02$.

• My question was about x close to zero... – mellow Dec 14 '13 at 19:22
• See my last sentence. For $|x|<1/7$, $1+\exp(-\tfrac1{x^2})$ will be numerically indistinguishable from $1$. For $|x|<1/28\approx 0.036$, evaluation of $\exp(-\tfrac1{x^2})$ will produce underflow for double floating point numbers. So using double f(double x) { return abs(x)<0.02?0:exp(-1/(x*x));} will be quite safe. – Lutz Lehmann Dec 14 '13 at 19:37
• My problem is precisely to distinguish $1+\exp(-1/x^2)$ from $1$ numerically. I want to fight with the underflow problem. – mellow Dec 15 '13 at 11:47
• Depending on the intended application, you may achieve that by keeping the exponents separate with a data structure for sums of terms $c_k\cdot \exp(a_k)$ instead of instantly carrying out the operations on doubles. This leads to (semi-)symbolic computations. One application of such a concept that I know of is Malajovic' extension of the Dandelin-Graeffe iteration for polynomial roots. – Lutz Lehmann Dec 15 '13 at 11:55

As has been pointed out in the comments on the question, the function is not analytic at $x=0$ and a Taylor expansion around this point is not possible. However, the function is analytic everywhere else. If you have an idea of the size of your (small) values for $y=x^2$, then you can of course still do one of the following:

• A Taylor expansion around a representative value for $y$.
• A low-degree polynomial approximation for the range of values $y$ you care about.
• A low-degree rational approximation.

Because the function is very smooth, you will likely good a fairly good approximation using any of these techniques.

$$y(x) = \exp(-1/x^2)$$ in C is exp(-1/(x*x)). Since most floating point in C is based on powers of 2, simple keep y in two parts: an integer power of 2 and a significand. Reconstitute y as needed with ldexp().

Calculation of z (see below) may lose up to 2 bits. y_expo2 will not lose anymore. y_signifcand may lose 1 bit. Reconstituting y should yield no more than a 4 bit loss of precision. But each x is halved, when less than 1e-3, another 2 bits is lost.

The helps maintain usable y_expo2 and y_significand down to |x| of about 1e-7. But smaller than that, y_expo2 has all the precision and y_signifcand little or none.

double f(double x) {
return exp(-1 / (x * x));
}

void fmellow(double x, double *y_expo2, double *y_signifcand) {
static const double ln2 = 0.69314718055994530941723212145818; // log(2.0);
double z = -1.0 / (x * x) / ln2;
// y = exp(z), but this rapidly rounds to 0.0
// Instead, break the floating-point number into a normalized fraction
//   and an integral power of 2.
double fraction;
fraction = modf(z, y_expo2);
*y_signifcand = pow(2.0, fraction);
}

void testf() {
double x = 1.0;
int i;
for (i = 0; i < 24; i++) {
double y = f(x);
double y_expo2;
double y_signifcand;
fmellow(x, &y_expo2, &y_signifcand);
double y_reconstituted = ldexp(y_signifcand, y_expo2);
double error = (y - y_reconstituted) / y_reconstituted;
printf("x=%8.2e  y=%9.2e  y_expo=%9.2e  y_signifcand=%9.2e  y_reconstituted=%9.2e  diff=%9.2e\n",
x, y, y_expo2, y_signifcand, y_reconstituted, error);
x *= 0.5;
}
}

int main() {
testf();
return 0;
}

x=1.00e+00  y= 3.68e-01  y_expo=-1.00e+00  y_signifcand= 7.36e-01  y_reconstituted= 3.68e-01  diff= 0.00e+00
x=5.00e-01  y= 1.83e-02  y_expo=-5.00e+00  y_signifcand= 5.86e-01  y_reconstituted= 1.83e-02  diff=-1.89e-16
x=2.50e-01  y= 1.13e-07  y_expo=-2.30e+01  y_signifcand= 9.44e-01  y_reconstituted= 1.13e-07  diff=-2.35e-16
x=1.25e-01  y= 1.60e-28  y_expo=-9.20e+01  y_signifcand= 7.94e-01  y_reconstituted= 1.60e-28  diff=-8.39e-16
x=6.25e-02  y=6.62e-112  y_expo=-3.69e+02  y_signifcand= 7.96e-01  y_reconstituted=6.62e-112  diff=-3.63e-15
x=3.12e-02  y= 0.00e+00  y_expo=-1.48e+03  y_signifcand= 8.01e-01  y_reconstituted= 0.00e+00  diff=      nan
x=1.56e-02  y= 0.00e+00  y_expo=-5.91e+03  y_signifcand= 8.24e-01  y_reconstituted= 0.00e+00  diff=      nan
x=7.81e-03  y= 0.00e+00  y_expo=-2.36e+04  y_signifcand= 9.23e-01  y_reconstituted= 0.00e+00  diff=      nan
x=3.91e-03  y= 0.00e+00  y_expo=-9.45e+04  y_signifcand= 7.26e-01  y_reconstituted= 0.00e+00  diff=      nan
x=1.95e-03  y= 0.00e+00  y_expo=-3.78e+05  y_signifcand= 5.55e-01  y_reconstituted= 0.00e+00  diff=      nan
x=9.77e-04  y= 0.00e+00  y_expo=-1.51e+06  y_signifcand= 7.60e-01  y_reconstituted= 0.00e+00  diff=      nan
x=4.88e-04  y= 0.00e+00  y_expo=-6.05e+06  y_signifcand= 6.69e-01  y_reconstituted= 0.00e+00  diff=      nan
x=2.44e-04  y= 0.00e+00  y_expo=-2.42e+07  y_signifcand= 7.99e-01  y_reconstituted= 0.00e+00  diff=      nan
x=1.22e-04  y= 0.00e+00  y_expo=-9.68e+07  y_signifcand= 8.16e-01  y_reconstituted= 0.00e+00  diff=      nan
x=6.10e-05  y= 0.00e+00  y_expo=-3.87e+08  y_signifcand= 8.89e-01  y_reconstituted= 0.00e+00  diff=      nan
x=3.05e-05  y= 0.00e+00  y_expo=-1.55e+09  y_signifcand= 6.24e-01  y_reconstituted= 0.00e+00  diff=      nan
x=1.53e-05  y= 0.00e+00  y_expo=-6.20e+09  y_signifcand= 6.07e-01  y_reconstituted= 0.00e+00  diff=      nan
x=7.63e-06  y= 0.00e+00  y_expo=-2.48e+10  y_signifcand= 5.44e-01  y_reconstituted= 0.00e+00  diff=      nan
x=3.81e-06  y= 0.00e+00  y_expo=-9.91e+10  y_signifcand= 7.01e-01  y_reconstituted= 0.00e+00  diff=      nan
x=1.91e-06  y= 0.00e+00  y_expo=-3.97e+11  y_signifcand= 9.68e-01  y_reconstituted= 0.00e+00  diff=      nan
x=9.54e-07  y= 0.00e+00  y_expo=-1.59e+12  y_signifcand= 8.78e-01  y_reconstituted= 0.00e+00  diff=      nan
x=4.77e-07  y= 0.00e+00  y_expo=-6.35e+12  y_signifcand= 5.95e-01  y_reconstituted= 0.00e+00  diff=      nan
x=2.38e-07  y= 0.00e+00  y_expo=-2.54e+13  y_signifcand= 5.03e-01  y_reconstituted= 0.00e+00  diff=      nan
x=1.19e-07  y= 0.00e+00  y_expo=-1.02e+14  y_signifcand= 5.11e-01  y_reconstituted= 0.00e+00  diff=      nan