expm1l: Work around a NetBSD 10.0/i386 bug

expm1l: Work around a NetBSD 10.0/i386 bug

[Bruno Haible]

On NetBSD 10.0/i386 I see this test failure:

FAIL: test-expm1l
=================

../../gltests/test-expm1.h:41: assertion 'y >= x' failed
FAIL test-expm1l (exit status: 134)

The cause is that for arguments such as x = 2^-1075, expm1l(x) returns 0,
which it never should. (exp(x) - 1 is always ≥ x.)

This patch provides a workaround.


2024-04-06  Bruno Haible  <bruno@clisp.org>

	expm1l: Work around a NetBSD 10.0/i386 bug.
	* m4/expm1l.m4 (gl_FUNC_EXPM1L): Test the value of
	expm1l(2^LDBL_MIN_EXP).
	* doc/posix-functions/expm1l.texi: Mention the NetBSD 10.0/i386 bug.
	* tests/test-expm1.h (test_function): Allow more deviation on
	NetBSD/i386.

diff --git a/doc/posix-functions/expm1l.texi b/doc/posix-functions/expm1l.texi
index f8c68c9d45..f8a015b3da 100644
--- a/doc/posix-functions/expm1l.texi
+++ b/doc/posix-functions/expm1l.texi
@@ -15,9 +15,12 @@
 This function is not declared on some platforms:
 IRIX 6.5.
 @item
+This function returns 0 instead of small positive values on some platforms:
+NetBSD 10.0/i386.
+@item
 This function produces results which are accurate to only 16 digits on some
 platforms:
-musl libc 1.2.2/arm64, musl libc 1.2.2/s390x, Mac OS X 10.5, NetBSD 10.0.
+musl libc 1.2.2/arm64, musl libc 1.2.2/s390x, Mac OS X 10.5, NetBSD 10.0/x86_64.
 @end itemize
 
 Portability problems not fixed by Gnulib:
diff --git a/m4/expm1l.m4 b/m4/expm1l.m4
index cdd51ae201..3e75891ac3 100644
--- a/m4/expm1l.m4
+++ b/m4/expm1l.m4
@@ -1,5 +1,5 @@
 # expm1l.m4
-# serial 13
+# serial 14
 dnl Copyright (C) 2010-2024 Free Software Foundation, Inc.
 dnl This file is free software; the Free Software Foundation
 dnl gives unlimited permission to copy and/or distribute it,
@@ -125,7 +125,7 @@ AC_DEFUN([gl_FUNC_EXPM1L]
   long double (* volatile my_expm1l) (long double) = argc ? expm1l : dummy;
   int result = 0;
   /* This test fails on musl 1.2.2/arm64, musl 1.2.2/s390x, Mac OS X 10.5,
-     NetBSD 10.0.  */
+     NetBSD 10.0/x86_64.  */
   {
     const long double TWO_LDBL_MANT_DIG = /* 2^LDBL_MANT_DIG */
       (long double) (1U << ((LDBL_MANT_DIG - 1) / 5))
@@ -141,6 +141,18 @@ AC_DEFUN([gl_FUNC_EXPM1L]
     if (!(err >= -100.0L && err <= 100.0L))
       result |= 1;
   }
+  /* This test fails on NetBSD 10.0/i386.  */
+  {
+    int i;
+    long double x;
+    volatile long double y;
+    for (i = -1, x = 0.5L; i > LDBL_MIN_EXP; i--, x *= 0.5L)
+      ;
+    /* Here i = LDBL_MIN_EXP, x = 2^LDBL_MIN_EXP.  */
+    y = my_expm1l (x);
+    if (!(y >= x))
+      result |= 2;
+  }
   return result;
 }
             ]])],
diff --git a/tests/test-expm1.h b/tests/test-expm1.h
index 01128c491d..fc92617a05 100644
--- a/tests/test-expm1.h
+++ b/tests/test-expm1.h
@@ -52,7 +52,9 @@ test_function (void)
     /* Error bound, in ulps.  */
     const DOUBLE err_bound =
       (sizeof (DOUBLE) > sizeof (double) ?
-#if defined __i386__ && defined __FreeBSD__
+#if defined __i386__ && defined __NetBSD__
+       L_(1025.0)
+#elif defined __i386__ && defined __FreeBSD__
        /* On FreeBSD/x86 6.4, the 'long double' type really has only 53 bits of
           precision in the compiler but 64 bits of precision at runtime.  See
           <https://lists.gnu.org/r/bug-gnulib/2008-07/msg00063.html>.
@@ -83,6 +85,8 @@ test_function (void)
       (sizeof (DOUBLE) > sizeof (double) ?
 #if defined __i386__ && defined __FreeBSD__
        L_(1536.0)
+#elif defined __i386__ && defined __NetBSD__
+       L_(1026.0)
 #else
        L_(11.0)
 #endif