Here's what I don't get. Why isn't Math.expm1 just implemented in JavaScript?
Math.expm1 = (x) => Math.pow(Math.E, (x)) - 1;
With this monkey patched version, Math.expm1(-0) now returns 0. I've been writing JavaScript for a long time and I didn't even know Math.expm1 was a thing. Was it really necessary to implement this inside V8?
expm1 avoids catastrophic cancellation (due to finite precision floating point), and so is more accurate for x close to zero than naive(x) = exp(x) - 1, e.g. expm1(1e-100) != 0, but naive(1e-100) = 0.
Even good old C includes several functions that seem "redundant" at first glance [0] but are necessary for edge cases. I suppose this is similar to mechanics like fused multiply-add [1].
The reason expm1 exists is because that naive approach is numerically poor when x is small (then e^x is almost 1, and subtracting two almost equal values throws away significant digits and gives you an inaccurate result).
I don't know what the rules for JavaScript semantics are - this version will break if someone redefines `Math.pow`. Is that ok or does that break JavaScript semantics?
That's the wrong answer and your algorithm is bad. The point of expm1 is to be accurate for small arguments using the Taylor series or Chebeyshev polynomial representation.
I know what a Padé approximant is but I am not sure how it is relevant. expm1 is available as part of the math library in JavaScript, just as it is also available for C. You do not have to implement it yourself.
The exact implementation depends on your math library, but in glibc it appears to reduce the argument modulo log 2, and then uses an approximation found in part by the Remez algorithm for the exponential function on the range [0, 0.5 log2]. My assumption here is that the Remez algorithm will often be preferred over Padé, but I am not an expert on these things. The Padé approximant gives an approximation “near” a point but the Remez algorithm finds an approximation over an entire interval. In either case, the Remez algorithm is only used to construct one part of the function, and there is also some effort to propagate error from the argument reduction. The resulting algorithm claims <1 ULP error.
For matrix exponentiation you would need a library.