ld128 coshl mishandles tiny inputs

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).

JuliaMath / openlibm

ld128 coshl mishandles tiny inputs #285