[BUG]: `math/base/special/expm1` overflows prematurely

[donhatch]

Description

Stdlib's expm1(709.782712893384) incorrectly returns Infinity,
whereas five other similar expressions all correctly return 1.7976931348622732e+308
(which is the correctly-rounded-to-nearest-representable answer for all 6 expressions,
as I'll show below, and is safely hundreds of ulps away from overflowing):

  Math.exp(709.782712893384)   = 1.7976931348622732e+308
       exp(709.782712893384)   = 1.7976931348622732e+308
  Math.exp(709.782712893384)-1 = 1.7976931348622732e+308
       exp(709.782712893384)-1 = 1.7976931348622732e+308
Math.expm1(709.782712893384)   = 1.7976931348622732e+308
     expm1(709.782712893384)   = Infinity

This is especially surprising and noticeable since it causes expm1's overflow threshold to be smaller than that of the mathematically-always-larger function exp.

Why 1.7976931348622732e+308 is the correctly-rounded-to-nearest answer

First, let's find the exact mathematical value of the input javascript/C/C++/python literal 709.782712893384.
In python3:

  >>> '%.10000g' % 709.782712893384
  '709.7827128933839730962063185870647430419921875'

Calculating exp (or expm1, the distinction is negligible) of that in high precision in worlframalpha:

  exp(709.7827128933839730962063185870647430419921875)

yields:

  1.797693134862273217839649630900041649872882... × 10^308

Using that as a literal in python3 should produce the correctly-rounded answer:

  >>> 1.797693134862273217839649630900041649872882e+308
  1.7976931348622732e+308

And javascript agrees:

  > 1.797693134862273217839649630900041649872882e+308
  1.7976931348622732e+308

Conclusion: stdlib's exp(), exp()-1, Math.exp(), Math.exp()-1, Math.expm1() all give the correctly-rounded-to-nearest answer 1.7976931348622732e+308,
and stdlib's expm1() does not.

Is Infinity a correct answer?

Despite 1.7976931348622732e+308 being the correctly-rounded-to-nearest answer as shown above, it is still conceivable that Infinity could be a correct answer.
For that to be the case, it would have to be that the mathematically-exact answer is within 1.5 (or so... I think) ulps of overflowing.
But that is not the case: Number.MAX_VALUE is 1.7976931348623157e+308, which is many (213 or so) ulps larger, as the reproducer program below shows.
So, Infinity is not a correct answer.

Other considerations

I noticed that the expectation that expm1's first overflow is (incorrectly) at exactly literal 709.782712893384 (which is the same number as literal 7.09782712893383973096e+02) is currently hard-coded here in the source code of the expm1rel function at https://github.com/stdlib-js/math-base-special-expm1rel/blob/4477702f9201e83666d913937b06a3667dc430f3/lib/main.js :

var HUGE_VALUE = 7.09782712893383973096e+02; // 0x40862E42 0xFEFA39EF
...
if ( x < HUGE_VALUE ) {
    return expm1( x ) / x;
}
... alternative overflow-resistant implementation ...

So if the implementation of expm1 gets fixed so that it no longer overflows prematurely for this value, then presumably this definition of HUGE_VALUE then needs to be updated to the next higher floating-point number (at which expm1 correctly overflows).

In practice, failing to update HUGE_VALUE wouldn't be a disaster; it would just mean expm1rel(709.782712893384) would still use the slightly-less-accurate-in-general overflow-resistant alternative formula, which happens to return exactly the same answer in this case. But it seems bad form to depend on this accidental fact.

It seems to me that the strategy of depending on this magic HUGE_VALUE constant is perhaps not the best strategy to begin with.
A possibly-more-robust strategy would be to get rid of this constant, and instead make the code simply evaluate expm1(x) and test whether it returns Infinity, instead.
That would free expm1rel's implementation from such tight dependence on accidental properties of the expm1 function and implementation.

Related Issues

Related issue #12436 fixed premature overflow of expm1rel, whose implementation is still perhaps overly-entangled with the overflow threshold of expm1, as noted above.

Questions

No.

Demo

No response

Reproduction

For bare minimal reproduction, one can simply print out the values of the 6 expressions shown at the beginning of this issue report.

The following script provides that and more context, confirming that other nearby values are behaving correctly, and confirming other claims made above.

import exp from '@stdlib/math-base-special-exp';  // pnpm install @stdlib/math-base-special-exp
import expm1 from '@stdlib/math-base-special-expm1';  // pnpm install @stdlib/math-base-special-expm1
import nextafter from 'nextafter';  // pnpm install nextafter

const pred = x => nextafter(x, -Infinity);  // predecessor
const succ = x => nextafter(x, Infinity);   // successor

const data = {
  "-Infinity": -Infinity,
  "zero": 0,
  "One": 1,
  "Last value for which both Math.expm1 and stdlib's expm1 correctly returns a finite number": 709.7827128933839,
  "The one value for which Math.exmp1 returns a finite answer but stdlib's expm1 overflows":   709.782712893384,
  "HUGE_VALUE 7.09782712893383973096e+02 in expm1rel source code (same as previous)":          7.09782712893383973096e+02,
  "First value for which both Math.expm1 and stdlib's expm1 return Infinity":                  709.7827128933841,
  "Infinity": Infinity,
  "NaN": NaN,
};

for (const key in data) {
  const x = data[key];
  console.log("  "+key+": "+x);
  if (x !== -Infinity) console.log("      pred("+x+") = "+pred(x));
  if (x !== Infinity) console.log("      succ("+x+") = "+succ(x));

  console.log("        Math.exp("+x+")   = "+Math.exp(x));
  console.log("             exp("+x+")   = "+exp(x));
  console.log("        Math.exp("+x+")-1 = "+(Math.exp(x)-1));
  console.log("             exp("+x+")-1 = "+(exp(x)-1));
  console.log("      Math.expm1("+x+")   = "+Math.expm1(x));
  console.log("           expm1("+x+")   = "+expm1(x));
  if (!Object.is(expm1(x), Math.expm1(x))) {  // use Object.is instead of === so NaN will be considered equal to NaN
    console.warn("HEY! Math.expm1 and stdlib's expm1 gave different answers!");
  }
  console.log();
}

console.log("  HUGE_VAL = 7.09782712893383973096e+02 = "+7.09782712893383973096e+02);  // confirm that it's 709.782712893384 as claimed

const number_of_ulps_to_infinity = x => {
  let answer = 0;
  while (x < Infinity) {
    answer++;
    x = succ(x);
  }
  return answer;
};
console.log("  number_of_ulps_to_infinity(1.7976931348622732e+308) = "+number_of_ulps_to_infinity(1.7976931348622732e+308));  // 214

Expected Results

...
  Math.exp(709.782712893384)   = 1.7976931348622732e+308
       exp(709.782712893384)   = 1.7976931348622732e+308
  Math.exp(709.782712893384)-1 = 1.7976931348622732e+308
       exp(709.782712893384)-1 = 1.7976931348622732e+308
Math.expm1(709.782712893384)   = 1.7976931348622732e+308
     expm1(709.782712893384)   = 1.7976931348622732e+308
...

Actual Results

...
  Math.exp(709.782712893384)   = 1.7976931348622732e+308
       exp(709.782712893384)   = 1.7976931348622732e+308
  Math.exp(709.782712893384)-1 = 1.7976931348622732e+308
       exp(709.782712893384)-1 = 1.7976931348622732e+308
Math.expm1(709.782712893384)   = 1.7976931348622732e+308
     expm1(709.782712893384)   = Infinity
HEY! Math.expm1 and stdlib's expm1 gave different answers!
...

Version

pnpm list says: @stdlib/math-base-special-expm1@0.2.4

Environments

N/A

Browser Version

No response

Node.js / npm Version

No response

Platform

No response

Checklist

• Read and understood the Code of Conduct.
• Searched for existing issues and pull requests.

[donhatch]

changing title following lead of #12436 (comment)

[donhatch]

Possibly-relevant observation...

I see in https://github.com/stdlib-js/stdlib/blob/develop/lib/node_modules/%40stdlib/math/base/special/expm1/lib/main.js the following:

var OVERFLOW_THRESHOLD = 7.09782712893383973096e+02; // 0x40862E42 0xFEFA39EF
...
/**
...
* ## Notes
...
* -   For IEEE 754 double, if \\(x > 7.09782712893383973096\mbox{e+}02\\), then \\(\operatorname{expm1}(x)\\) will overflow.
...
*/
...
                if ( y >= OVERFLOW_THRESHOLD ) { // if |x| >= 709.78...
                        return PINF;
                }

(recall that literal 7.09782712893383973096e+02 is the same value as literal 709.782712893384, i.e. it's the one input value for which expm1 is misbehaving).

My reading of this is that this note in the comment is correct, but that the code incorrectly uses it.
This could be fixed by simply changing the '>=' to '>' in the code.
I tried this, in a scratch local copy of stdlib-js, and it appears to fix the problem.

It also looks like the same comment and incorrect code occurs in the C implementation in https://github.com/stdlib-js/stdlib/blob/develop/lib/node_modules/%40stdlib/math/base/special/expm1/src/main.c (although I have never tested that, as far as I know).