4 ms·
For this specific case, use Python's expm1, see https://docs.python.org/3/library/math.html#math.expm1 https://docs.python.org/3/library/math.html#math.expm1 an
by eesmith 3y ago
For this specific case, use Python's expm1, see https://docs.python.org/3/library/math.html#math.expm1 https://docs.python.org/3/library/math.html#math.expm1 and history of expm1 at https://en.wikipedia.org/wiki/Exponential_function#expm1 https://en.wikipedia.org/wiki/Exponential_function#expm1
def f(x):
return math.expm1(x)/x
The expm1(x) means "exp(x) minus 1".
>>> f(1e-15)
1.0000000000000007
The technique described gives
1.0000000000000004
which is 1 step smaller than the value computed via expm1():
>>> math.nextafter(f(1e-15), 0)
1.0000000000000004
Both are within 1 ulp of the more precise value from WolframAlpha:
1.0000000000000005000000000000001666666666666667083333333333333416...
- Sharlin 3y agoTo be precise, it’s C’s `expm1`. Like other things in math.h, it’s just adopted as is to most popular languages, even including PHP and JS.
- eesmith 3y agoTo be precise, it's most likely Kahan's `expm1`, since Wikipedia attributes it to him. It adds that expm1 was first implemented on an HP calculator in the 1970s. It became part of BSD 4.3 due to Kahan and one of his students. The documentation at https://archive.org/details/unix-programmers-reference-manual-4.3-berkeley-software-distribution/page/n161/mode/2up?q=kahan+expm1 https://archive.org/details/unix-programmers-reference-manua... mentions "expm1" was available in HP BASIC, and mentions the Apple C compiler at the time used "exp1", which you can verify from the Apple Numerics Manual Second Edition 1988 at https://archive.org/details/apple-numerics-manual-second-edition-1988/page/n85/mode/2up?q=exp1+compiler https://archive.org/details/apple-numerics-manual-second-edi... That makes it BASIC's expm1 before it was C's. ;) While yes, Python's expm1 uses the platform libm's expm1 through C, Python's math module will avoid using the vendor libm for cases where the vendor implementation is known to have numerical issues. For one such example, on macOS the following: #include <math.h> #include <stdio.h> int main(void) { double x = -1.0; x = nextafter(x, 0.0); // -0.9999999999999999 printf("lgamma(%.19f) -> %.17f\n", x, lgamma(x)); return 0; } reports lgamma(-0.9999999999999998890) -> 36.25168935623486988 while Python, using its own lgamma() implementation: import math x = math.nextafter(-1.0, 0.0) print(f"lgamma({x}) -> {math.lgamma(x)}") reports: lgamma(-0.9999999999999999) -> 36.7368005696771 showing that Python's lgamma is not a pass-through to C's. FWIW, WolframAlpha says the value is 36.73680056... so Python's is more accurate. EDIT: Python's math.expm1 was added in version 3.2 back in February 20th, 2011, and did not depend on C's expm1. Instead, it used Kahan's method directly: /* Mathematically, expm1(x) = exp(x) - 1. The expm1 function is designed to avoid the significant loss of precision that arises from direct evaluation of the expression exp(x) - 1, for x near 0. */ double _Py_expm1(double x) { /* For abs(x) >= log(2), it's safe to evaluate exp(x) - 1 directly; this also works fine for infinities and nans. For smaller x, we can use a method due to Kahan that achieves close to full accuracy. */ if (fabs(x) < 0.7) { double u; u = exp(x); if (u == 1.0) return x; else return (u - 1.0) * x / log(u); } else return exp(x) - 1.0; }
- zokier 3y agoexpm1 is specified in IEEE 754-2008 (didn't bother to check old revisions)
- jjgreen 3y agoOfficially arrived in C99 if my memory serves ...