ld128 coshl mishandles tiny inputs · Issue #285 · JuliaMath/openlibm · GitHub
Skip to content

ld128 coshl mishandles tiny inputs #285

Open

Description

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

Activity

  1. ararslan commented on Jan 17, 2024

    @ararslan
    Member

    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.

  2. added
    bugIncorrect or unexpected results, crashes, etc.
    on Jan 19, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugIncorrect or unexpected results, crashes, etc.

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions