To 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-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-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
#include
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)