Hyperbolic functions: move the libm values 3 floats outward, or bound them without the libm?

[Jordan08]

The tests added to this fork (tests/elementary.cpp, commit 0e2504f) found bounds of GAOL's hyperbolic functions that do not enclose the exact values. On Linux x86_64 with glibc 2.31, it happens for asinh(-0x1.ee84df02a8766p-4) and acosh(0x1.01fd62fff333fp+0). Two fixes are ready, one on master and one on the branch hyperbolic-rigorous. This issue is to choose between them.

Why the bounds are wrong

mathlib, which bounds GAOL's other elementary functions, has no hyperbolic function; crlibm has sinh and cosh only. GAOL therefore takes sinh, cosh, tanh, asinh, acosh and atanh from the libm of the system, rounded to nearest, and moves the value one float outward (gaol_double_op_apmathlib.h, gaol_double_op_crlibm.h).

That gives a bound only when the libm returns one of the two doubles around the exact value, and libms are not always that accurate. On 2000 random arguments of each function (exact values computed by mpmath with 2000 bits), the libm returned a double beyond those two this many times:

libm	sinh	cosh	tanh	asinh	acosh	atanh
glibc 2.31, x86_64	2	0	86	0	9	2
musl, x86_64 and aarch64	2	3	80	0	9	2
MinGW-w64, x64 (under Wine)	2	0	86	0	9	2
MinGW-w64, x86 (under Wine)	0	0	0	0	0	0
glibc 2.44, aarch64 and armv7	0	0	0	0	0	0

The libm was never more than one float beyond, but each of these times GAOL's bound on that side did not enclose the exact value.

Codac calls gaol::sinh(), gaol::cosh(), etc. for its intervals, and builds GAOL from the official repository, so it gets the same bounds.

Solution 1: move the libm values 3 floats outward (master, commit fdabc12)

• Change: two helpers, called instead of previous_float() and next_float() in the functions taken from the libm (gaol_double_op_apmathlib.h, and tanh, acosh, asinh and atanh in gaol_double_op_crlibm.h).
• Guarantee: the bounds enclose the exact values as long as the libm is within two floats of them. That is still an assumption about the libm. The tests check it on each platform of the continuous integration, but it is not proven.
• Tightness: bounds 3 floats from the tightest ones, sometimes 2 or 4.
• Speed: 30 to 60 ns more per call, for four more nextafter().

Solution 2: bound them without the libm (branch hyperbolic-rigorous, commit ead154d)

• Change: about 300 lines in gaol_interval.cpp, computed in interval arithmetic with exp_dn(), exp_up(), log_dn() and log_up():
  • sinh, cosh: Taylor series below 1 (the terms after the tenth are bounded), (e^x ∓ e^-x)/2 up to 20, then e^x/2 and one float beyond;
  • tanh: sinh/cosh below 1, 1 - 2/(e^2x + 1) above;
  • asinh, acosh, atanh: the libm value is only a starting point. Each bound is the nearest double that the enclosures of sinh, of cosh - 1 = 2 sinh(x/2)^2 and of tanh prove to be a bound, searched double by double. A formula with log() is used when the libm is more than 16 doubles away.
  • Near 0 and beyond 2^30, the first terms of the series and of the asymptotic expansions give the bounds.
• Guarantee: it rests only on mathlib (or crlibm) bounding exp and log, as GAOL's other functions do.
• Tightness: bounds 0 to 3 floats from the tightest ones (4 once for sinh), 0.3 to 1.9 on average.
• Speed: 6 to 45 times slower.

Both pass all the tests of tests/. Neither changes the variant that uses only the libm (gaol_double_op_m.h, without mathlib or crlibm).

Measurements

2000 random arguments of each function; Linux x86_64, glibc 2.31, GCC 9, -O2.

Bounds that do not enclose the exact value:

	sinh	cosh	tanh	asinh	acosh	atanh
GAOL (libm, 1 float)	2	0	86	0	9	2
Solution 1 (libm, 3 floats)	0	0	0	0	0	0
Solution 2 (without the libm)	0	0	0	0	0	0

Mean number of floats between the bounds and the tightest ones (largest in parentheses):

	sinh	cosh	tanh	asinh	acosh	atanh
GAOL (libm, 1 float)	0.99 (1)	0.99 (1)	0.98 (1)	1 (1)	1 (1)	1 (1)
Solution 1 (libm, 3 floats)	2.99 (4)	2.99 (3)	3.02 (4)	3 (3)	3.00 (4)	3.00 (4)
Solution 2 (without the libm)	0.29 (4)	0.29 (2)	1.02 (3)	1.88 (3)	1.63 (3)	1.08 (3)

Time per call on a point interval, in ns (best of 3 runs of 10^6 calls). Arguments are in [-5, 5] for sinh, cosh and tanh, [-100, 100] for asinh, [1, 100] for acosh, and [-0.99, 0.99] for atanh.

	sinh	cosh	tanh	asinh	acosh	atanh
GAOL (libm, 1 float)	87	59	90	77	63	103
Solution 1 (libm, 3 floats)	145	89	144	121	108	144
Solution 2 (without the libm)	519	513	728	1131	1491	4522

To decide

• Solution 1: keep master as it is
• Solution 2: replace the commit of solution 1 on master by the commit of hyperbolic-rigorous

[Jordan08]

The first run of the continuous integration (commit dde06ad) found libms much less accurate than the ones measured above.

Solution 1 does not give bounds with MinGW-w64 11, 12 and 13 (x86 and x64)

• acosh near 1. The exact value of acosh(0x1.0000000000001p+0) lies between 0x1.6a09e667f3bccp-26 and 0x1.6a09e667f3bcdp-26.
  • MinGW-w64 12.2 x86: the libm returns about 15,000 floats below it, and GAOL gives [0x1.6a09e667efffdp-26, 0x1.6a09e667f0003p-26], whose upper bound is below the exact value.
  • MinGW-w64 12.2 x64: the libm returns about 25 million floats above it, and GAOL gives [0x1.6a09e67fffffcp-26, 0x1.6a09e68000002p-26], whose lower bound is above the exact value.
  • The same happens for other arguments near 1, with 11.2, 12.2 and 13.2.
• asinh of large negative numbers. With MinGW-w64 11.2 the libm returns NaN, e.g. asinh(-0x1.fffffffffffffp+1023) = [nan, nan], and so does asinh_rel().

MinGW-w64 14.2 and 15.2 pass, as do Visual Studio 2022 and 2026 (x86, x64, arm64), MSYS2, macOS, Linux (glibc and musl) and the other architectures. No number of floats is enough for these libms.

Solution 2 on the same platforms

The branch hyperbolic-rigorous-ci is master with the commit of solution 2 in place of the commit of solution 1 (532c9a4 reverts fdabc12, 1b6e0f5 is ead154d). It runs the same workflows.

• Everything passes where master passes.
• asinh no longer fails with MinGW-w64 11 to 13.
• acosh near 1 still fails with MinGW-w64 11, 12 and 13, on x86 and x64:
  • Loose bounds. When the libm is more than 16 floats away, the bounds fall back on log(x + sqrt((x-1)(x+1))), which encloses the exact value but is loose near 1: acosh(0x1.0000000000001p+0) = [0x1.6a09e63ffffffp-26, 0x1.6a09e67ffffefp-26] with MinGW-w64 x64, some 2^26 floats wide, where the tests allow 8.
  • One wrong result. acosh(0x1.0000008f56f9ap+0) gives the inverted interval [0x1.0ee7e22f1025ep-12, 0x1.0ee7e22f1024fp-12], while the exact value is 0x1.0ee7e22f10569p-12. It happens with MinGW-w64 11.2, 12.2 and 13.2 x64 and 11.2 x86.
  • Not reproduced on Linux. Moving the value of the libm by as many floats there (−793, ±20, ±25 million) does not reproduce it: the bounds stay enclosures, looser when the fallback is taken. Something else in these MinGW-w64 versions is involved, and it is not understood yet.

Where this leaves the choice

As they stand, neither solution is right with MinGW-w64 11 to 13:

• solution 1, because their libms are millions of floats away;
• solution 2, because of its loose fallback near 1 and a wrong result not explained yet.

[Jordan08]

MinGW-w64 older than version 12 is now refused (commit 62d0107)

The branch diag-mingw-acosh traced the acosh() of solution 2 with MinGW-w64 13.2 (x64, Release and Debug), 11.2 (x86), and 15.2 (x64) for comparison.

What the libm returns. With MinGW-w64 13.2:

• acosh(0x1.0000008f56f9ap+0) is 0x1.0ee7e22f1024ep-12, 793 floats below the exact value 0x1.0ee7e22f10569p-12;
• acosh(0x1.0000000000001p+0) is 0x1.6a09e67ffffffp-26, 25 million floats above the exact value.

The rounding direction is right around the call: upward before and after it, MXCSR 0x5fa0.

What solution 2 gives.

• Traced, it gives enclosures, loose when the fallback on log() is taken: [0x1.0ee7e22f1025ep-12, 0x1.0ee7e22f1124ep-12].
• Not traced, the same code gave the inverted [0x1.0ee7e22f1025ep-12, 0x1.0ee7e22f1024fp-12]. Its search for the upper bound stopped one float above the value of the libm, as if the check had passed there, which the trace shows it does not.
• The difference comes from how these compilers build the code. It was not narrowed down further.

Only these versions. The problems of both solutions only happen with the mingw-w64 of these versions:

• MinGW-Builds with runtime rt_v10 (GCC 12.2) and rt_v11 (GCC 13.2);
• the WinLibs build of GCC 11.2.

MinGW-w64 14.2 (rt_v12) and 15.2 (rt_v13) pass.

The refusal. The CMake build of master now stops with a message when mingw-w64 is older than version 12 (__MINGW64_VERSION_MAJOR), naming the MinGW-w64 to use instead. The continuous integration checks that the MinGW-w64 11.2, 12.2 and 13.2 of Chocolatey are refused.

Apart from these versions, both solutions passed on all the platforms of the continuous integration, so the choice rests on the comparison above: guarantee, tightness and speed.

[Jordan08]

Work on solution 1 since then

Solution 1 is still the one on master, and the choice between the two solutions is still open. On top of it:

• Exact and tightest values at special points (7d99807).
  • sinh, tanh, asinh and atanh of 0 are 0, cosh(0) is 1 and acosh(1) is 0. The libm values moved three floats outward gave [-3·2^-1074, 3·2^-1074] there.
  • cosh and sinh are [MAX, +oo] from 711 on, and the lower bound of tanh is the double below 1 from 20 on.
  • The three assertions of GAOL's own check trigonometric that solution 1 broke pass again: cosh([MAX, +inf]), sinh([-inf, -MAX]) and tanh([MAX, +inf]). make check passes its 15 checks.
• Speed (3b63390). The rounding direction is now set once for both bounds rather than for each bound. On an Intel i7-1185G7 (GCC 9.4, glibc 2.31), per interval:

	before	after
sinh()	140 ns	128 ns
cosh()	81 ns	72 ns
tanh()	137 ns	126 ns
asinh()	122 ns	111 ns
acosh()	106 ns	96 ns
atanh()	138 ns	127 ns

• Documentation. doc/accuracy.md (0536f71) documents the tightness of these functions with solution 1: valid as long as the libm is within two floats of the exact value.

[Jordan08]

Where this stands (master, 17 September)

Solution 1 is still the one on master, and the choice is still open. Nothing changed today in the hyperbolic functions; what changed around them:

• mathlib is in the sources (3rd/mathlib, 49e0801): the builds no longer download it, and the 110 jobs of the continuous integration pass with solution 1, Windows (Visual C++, MinGW-w64 15.2, MSYS2), macOS, musl and the qemu architectures included.
• doc/continuous-integration.md now lists the configurations refused or left out with their reasons, MinGW-w64 11 to 13 for this issue among them.

A third way, to be looked at

Both solutions pay for the same lack: no correctly rounded hyperbolic function in mathlib nor in CRlibm (but sinh and cosh). CORE-MATH has correctly rounded sinh, cosh, tanh, asinh, acosh and atanh in double precision, under the MIT licence, one source file per function. Taken into 3rd/ as mathlib was:

• Guarantee: the one of solution 2, without an assumption on the libm, and without the libm at all: MinGW-w64 11 to 13 would no longer be refused for their acosh().
• Tightness: one double from the tightest bounds (the value moved one double outward), or the tightest ones if called rounding downward and upward, which CORE-MATH supports.
• Speed: to be measured; its functions are meant to be as fast as a libm.
• To check first: portability. Some of its functions use unsigned __int128 or fma(), which Visual C++ and the 32-bit targets (i386, armhf) do not have, or not fast: whether the six functions compile and stay correctly rounded on every platform of the continuous integration is the question.

Say if this should be tried on a branch, as solution 2 was.