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

[donhatch]

Description

expm1rel(x) prematurely returns Infinity for x in the interval [709.782712893384, 716.3568913878178].
It should return a finite number instead.

Related Issues

Original feature request: #115

Questions

No.

Demo

No response

Reproduction

The following javascript program demonstrates, comparing against an overflow-resistant alternative implementation. The magic input values were found by binary search.

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


// Overflow resistant expression for expm1rel(x):
//   expm1rel(x) = (exp(x)-1) / x
//               = (exp(x/2)-1) * (exp(x/2)+1) / x
//               = expm1(x/2) * (expm1(x/2) + 2) / x
//               = expm1(x/2) / x * (expm1(x/2) + 2)
const overflow_resistant_expm1rel = x => {
  if (x === 0) return 1;
  const temp = expm1(x/2);
  if (temp === Infinity) return Infinity;
  return temp / x * (temp + 2);
};
    
const data = {
  "Last value for which expm1rel correctly returns a finite number": 709.7827128933839,
  "First value for which expm1rel prematurely overflows":            709.782712893384,
  "A nice round number for which expm1rel prematurely overflows":    710,
  "Last value for which expm1rel prematurely overflows":             716.3568913878178,
  "First value for which expm1rel correctly returns Infinity":       716.3568913878179,
};
  
for (const key in data) {
  const x = data[key];
  console.log("  "+key+": "+x);
  console.log("                         expm1rel("+x+") = "+expm1rel(x));
  console.log("      overflow_resistant_expm1rel("+x+") = "+overflow_resistant_expm1rel(x));
  console.log();
}

Expected Results

expm1rel should return the same values as overflow_resistant_expm1rel,
up to floating-point roundoff error.

Actual Results

Output when run using nodejs:

  Last value for which expm1rel correctly returns a finite number: 709.7827128933839
                         expm1rel(709.7827128933839) = 2.532737276079728e+305
      overflow_resistant_expm1rel(709.7827128933839) = 2.532737276079728e+305

  First value for which expm1rel prematurely overflows: 709.782712893384
                         expm1rel(709.782712893384) = Infinity
      overflow_resistant_expm1rel(709.782712893384) = 2.5327372760800158e+305

  A nice round number for which expm1rel prematurely overflows: 710
                         expm1rel(710) = Infinity
      overflow_resistant_expm1rel(710) = 3.1464715016362133e+305

  Last value for which expm1rel prematurely overflows: 716.3568913878178
                         expm1rel(716.3568913878178) = Infinity
      overflow_resistant_expm1rel(716.3568913878178) = 1.7976931348621814e+308

  First value for which expm1rel correctly returns Infinity: 716.3568913878179
                         expm1rel(716.3568913878179) = Infinity
      overflow_resistant_expm1rel(716.3568913878179) = Infinity

Version

pnpm list reports: @stdlib/math-base-special-expm1rel@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.

[stdlib-bot]

👋 Hi there! 👋

And thank you for opening your first issue! We will get back to you shortly. 🏃 💨

If you have any further questions while waiting for a response, please join our Zulip community to chat with project maintainers and other community members.

[kgryte]

@donhatch Thank you for reporting this issue and thank you for suggesting the alternative formulation. I just pushed 0a659c4 which leverages that insight, while still preserving 1:1 results relative to expm1(x)/x for medium-sized x. The patch extends the range of computable values based on your suggestion above.

The patch is now available both in the main repo and the standalone repository (ref: https://github.com/stdlib-js/math-base-special-expm1rel/blob/4477702f9201e83666d913937b06a3667dc430f3/lib/main.js).

I am not sure our next standalone package release timeline. @Planeshifter May have more insight here.