A Double Precision Puzzle with the Gaussian

Some library computes $$e^{-\frac{x^2}{2}}$$ the following way:

xsq = fint(x * 1.6) / 1.6;
del = (x - xsq) * (x + xsq);
result = exp(-xsq * xsq * 0.5) * exp(-del * 0.5);

where fint(z) computes the floor of z.

Basically, x*x is rewritten as xsq*xsq+del. I have seen that trick once before, but I just can’t figure out where and why (except that it is probably related to high accuracy issues).

The answer is in the next post.

comments powered by Disqus
Tweet Submit to reddit
© 2006-16 Fabien Creative Commons License This work is licensed under a Creative Commons Attribution 4.0 International License.