ld128 coshl mishandles tiny inputs

[statham-arm]

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

[ararslan]

Naive question: could this branch be amended to return one instead of w?

Also, I checked the OpenBSD source (where openlibm's 128-bit long double definition for coshl came from) and the code hasn't changed since the version imported here, so it has the same issue. When the issue is fixed here it could be worthwhile for someone to contribute the fix back upstream as well.