JuliaMath / JuliaMath/openlibm
ld128 coshl mishandles tiny inputs
- Dominant language
- C
- Stars
- 632
- Forks
- 163
- Avg merge
- 4m
- Merged PRs (30d)
- 1
Description
The following test program computes cosh of 1, 1/2, 1/4, ..., 1/2^-128.
```c
#include
#include
int main(int argc, char **argv) {
long double d = 1.0;
for (int i = 0; i < 128; i++) {
printf("coshl(%La) = %La\n", d, coshl(d));
d /= 2;
}
}
```
For small ε, cosh(ε) ≈ 1+ε^2/2. So if we run on AArch64, which has 128-bit long doubles with a 112-bit mantissa, we expect that by the time the input is about 2^-56, the output will be indistinguishable from 1. And indeed this looks sensible to start off with:
```
coshl(0x1p-54) = 0x1.0000000000000000000000000008p+0
coshl(0x1p-55) = 0x1.0000000000000000000000000002p+0
coshl(0x1p-56) = 0x1p+0
coshl(0x1p-57) = 0x1p+0
```
but a little bit later, the outputs unexpectedly become wrong:
```coshl(0x1p-70) = 0x1p+0
coshl(0x1p-71) = 0x1p+0
coshl(0x1p-72) = 0x1.000000000000000001p+0
coshl(0x1p-73) = 0x1.0000000000000000008p+0
coshl(0x1p-74) = 0x1.0000000000000000004p+0
coshl(0x1p-75) = 0x1.0000000000000000002p+0
coshl(0x1p-76) = 0x1.0000000000000000001p+0
```
It looks as if the early exit code path from this special case in `ld128/e_coshl.c` is absent-mindedly returning the wrong variable:
```c
if (ex < 0x3fb80000) /* |x| < 2^-116 */
return w; /* cosh(tiny) = 1 */
```
But `w` is not 1: it's 1 + expm1(input), i.e. about exp(input), i.e. about 1+input (for small inputs).
Contributor guide
No contributing guide indexed for this repository
Assessment
This issue has not been assessed yet.