The following test program computes cosh of 1, 1/2, 1/4, ..., 1/2^-128.
#include <stdio.h>
#include <math.h>
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-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:
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).
The following test program computes cosh of 1, 1/2, 1/4, ..., 1/2^-128.
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:
but a little bit later, the outputs unexpectedly become wrong:
It looks as if the early exit code path from this special case in
ld128/e_coshl.cis absent-mindedly returning the wrong variable:But
wis not 1: it's 1 + expm1(input), i.e. about exp(input), i.e. about 1+input (for small inputs).