ldexp: Work around OpenBSD/mips64 bug

[Bruno Haible]

On OpenBSD 7.3/mips64 I see a test failure:

FAIL: test-ldexp
================

../../gltests/test-ldexp.h:114: assertion 'y == x' failed
Abort trap (core dumped)
FAIL test-ldexp (exit status: 134)

A simplified test case is:
  ldexp (1.9269695883136991774e-308, 0) returns 7.32274e-245

This patch provides a workaround.


2023-08-19  Bruno Haible  <bruno@clisp.org>

        ldexp: Work around OpenBSD/mips64 bug.
        * lib/math.in.h (ldexp): New declaration.
        * lib/ldexp.c: New file, based on lib/ldexpl.c.
        * lib/ldexpl.c: Moved the implementation to lib/ldexp.c. Just include
        it.
        * m4/math_h.m4 (gl_MATH_H_REQUIRE_DEFAULTS): Initialize GNULIB_LDEXP.
        (gl_MATH_H_DEFAULTS): Initialize REPLACE_LDEXP.
        * m4/ldexp.m4 (gl_FUNC_LDEXP): Require gl_MATH_H_DEFAULTS and
        gl_FUNC_ISNAND. Invoke gl_FUNC_LDEXP_WORKS. Set REPLACE_LDEXP. Consider
        it when setting LDEXP_LIBM.
        (gl_FUNC_LDEXP_WORKS): New macro.
        * modules/math (Makefile.am): Substitute GNULIB_LDEXP, REPLACE_LDEXP.
        * modules/ldexp (Files): Add lib/ldexp.c.
        (Depends-on): Add math, isnand.
        (configure.ac): Set GL_COND_OBJ_LDEXP. Invoke gl_MATH_MODULE_INDICATOR.
        (Makefile.am): Conditionally compile ldexp.c.
        * modules/ldexpl (Files): Add lib/ldexp.c.
        * doc/posix-functions/ldexp.texi: Mention the OpenBSD bug.

diff --git a/doc/posix-functions/ldexp.texi b/doc/posix-functions/ldexp.texi
index 4698bbad6f..4bc5c0ef9f 100644
--- a/doc/posix-functions/ldexp.texi
+++ b/doc/posix-functions/ldexp.texi
@@ -8,6 +8,10 @@
 
 Portability problems fixed by Gnulib:
 @itemize
+@item
+This function produces wrong results on some platforms:
+@c ldexp(1.9269695883136991774e-308,0) returns 7.32274e-245.
+OpenBSD 7.3/mips64.
 @end itemize
 
 Portability problems not fixed by Gnulib:
diff --git a/lib/ldexp.c b/lib/ldexp.c
new file mode 100644
index 0000000000..dd09ca3814
--- /dev/null
+++ b/lib/ldexp.c
@@ -0,0 +1,98 @@
+/* Multiply a 'float' by a power of 2.
+   Copyright 2002-2003, 2007-2023 Free Software Foundation, Inc.
+
+   This file is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as
+   published by the Free Software Foundation; either version 2.1 of the
+   License, or (at your option) any later version.
+
+   This file is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with this program.  If not, see <https://www.gnu.org/licenses/>.  */
+
+/* Written by Paolo Bonzini and Bruno Haible.  */
+
+#include <config.h>
+
+/* Specification.  */
+#include <math.h>
+
+# include <float.h>
+#ifdef USE_LONG_DOUBLE
+# include "fpucw.h"
+#endif
+
+/* Avoid some warnings from "gcc -Wshadow".
+   This file doesn't use the exp() function.  */
+#undef exp
+#define exp exponent
+
+#ifdef USE_LONG_DOUBLE
+# define FUNC ldexpl
+# define DOUBLE long double
+# define ISNAN isnanl
+# define DECL_ROUNDING DECL_LONG_DOUBLE_ROUNDING
+# define BEGIN_ROUNDING() BEGIN_LONG_DOUBLE_ROUNDING ()
+# define END_ROUNDING() END_LONG_DOUBLE_ROUNDING ()
+# define L_(literal) literal##L
+#else
+# define FUNC ldexp
+# define DOUBLE double
+# define ISNAN isnand
+# define DECL_ROUNDING
+# define BEGIN_ROUNDING()
+# define END_ROUNDING()
+# define L_(literal) literal
+#endif
+
+DOUBLE
+FUNC (DOUBLE x, int exp)
+{
+  DECL_ROUNDING
+
+  BEGIN_ROUNDING ();
+
+  /* Check for zero, nan and infinity. */
+  if (!(ISNAN (x) || x + x == x))
+    {
+      unsigned int uexp;
+      DOUBLE factor;
+
+      if (exp < 0)
+        {
+          /* Avoid signed integer overflow when exp == INT_MIN.  */
+          uexp = (unsigned int) (-1 - exp) + 1;
+          factor = L_(0.5);
+        }
+      else
+        {
+          uexp = exp;
+          factor = L_(2.0);
+        }
+
+      if (uexp > 0)
+        {
+          unsigned int bit;
+
+          for (bit = 1;;)
+            {
+              /* Invariant: Here bit = 2^i, factor = 2^-2^i or = 2^2^i,
+                 and bit <= uexp.  */
+              if (uexp & bit)
+                x *= factor;
+              bit <<= 1;
+              if (bit > uexp)
+                break;
+              factor = factor * factor;
+            }
+        }
+    }
+
+  END_ROUNDING ();
+
+  return x;
+}
diff --git a/lib/ldexpl.c b/lib/ldexpl.c
index 89917e7a10..186e963173 100644
--- a/lib/ldexpl.c
+++ b/lib/ldexpl.c
@@ -1,10 +1,6 @@
-/* Emulation for ldexpl.
-   Contributed by Paolo Bonzini
-
+/* Multiply a 'float' by a power of 2.
    Copyright 2002-2003, 2007-2023 Free Software Foundation, Inc.
 
-   This file is part of gnulib.
-
    This file is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as
    published by the Free Software Foundation; either version 2.1 of the
@@ -18,6 +14,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with this program.  If not, see <https://www.gnu.org/licenses/>.  */
 
+/* Written by Paolo Bonzini and Bruno Haible.  */
+
 #include <config.h>
 
 /* Specification.  */
@@ -38,52 +36,8 @@ ldexpl (long double x, int exp)
 
 #else
 
-# include <float.h>
-# include "fpucw.h"
-
-long double
-ldexpl (long double x, int exp)
-{
-  unsigned int uexp;
-  long double factor;
-  unsigned int bit;
-  DECL_LONG_DOUBLE_ROUNDING
-
-  BEGIN_LONG_DOUBLE_ROUNDING ();
-
-  /* Check for zero, nan and infinity. */
-  if (!(isnanl (x) || x + x == x))
-    {
-      if (exp < 0)
-        {
-          /* Avoid signed integer overflow when exp == INT_MIN.  */
-          uexp = (unsigned int) (-1 - exp) + 1;
-          factor = 0.5L;
-        }
-      else
-        {
-          uexp = exp;
-          factor = 2.0L;
-        }
-
-      if (uexp > 0)
-        for (bit = 1;;)
-          {
-            /* Invariant: Here bit = 2^i, factor = 2^-2^i or = 2^2^i,
-               and bit <= uexp.  */
-            if (uexp & bit)
-              x *= factor;
-            bit <<= 1;
-            if (bit > uexp)
-              break;
-            factor = factor * factor;
-          }
-    }
-
-  END_LONG_DOUBLE_ROUNDING ();
-
-  return x;
-}
+# define USE_LONG_DOUBLE
+# include "ldexp.c"
 
 #endif
 
diff --git a/lib/math.in.h b/lib/math.in.h
index f841a1356e..96abf10b81 100644
--- a/lib/math.in.h
+++ b/lib/math.in.h
@@ -1434,6 +1434,29 @@ _GL_WARN_ON_USE (ldexpf, "ldexpf is unportable - "
 # endif
 #endif
 
+/* Return x * 2^exp.  */
+#if @GNULIB_LDEXP@
+# if @REPLACE_LDEXP@
+#  if !(defined __cplusplus && defined GNULIB_NAMESPACE)
+#   undef ldexp
+#   define ldexp rpl_ldexp
+#  endif
+_GL_FUNCDECL_RPL (ldexp, double, (double x, int exp));
+_GL_CXXALIAS_RPL (ldexp, double, (double x, int exp));
+# else
+/* Assume ldexp is always declared.  */
+_GL_CXXALIAS_SYS (ldexp, double, (double x, int exp));
+# endif
+# if __GLIBC__ >= 2
+_GL_CXXALIASWARN (ldexp);
+# endif
+#elif defined GNULIB_POSIXCHECK
+# undef ldexp
+/* Assume ldexp is always declared.  */
+_GL_WARN_ON_USE (ldexp, "ldexp is unportable - "
+                 "use gnulib module ldexp for portability");
+#endif
+
 /* Return x * 2^exp.  */
 #if @GNULIB_LDEXPL@ && @REPLACE_LDEXPL@
 # if !(defined __cplusplus && defined GNULIB_NAMESPACE)
diff --git a/m4/ldexp.m4 b/m4/ldexp.m4
index ef2c4c2825..513d3807c1 100644
--- a/m4/ldexp.m4
+++ b/m4/ldexp.m4
@@ -1,4 +1,4 @@
-# ldexp.m4 serial 1
+# ldexp.m4 serial 2
 dnl Copyright (C) 2010-2023 Free Software Foundation, Inc.
 dnl This file is free software; the Free Software Foundation
 dnl gives unlimited permission to copy and/or distribute it,
@@ -6,6 +6,9 @@
 
 AC_DEFUN([gl_FUNC_LDEXP],
 [
+  AC_REQUIRE([gl_MATH_H_DEFAULTS])
+  AC_REQUIRE([gl_FUNC_ISNAND]) dnl for ISNAND_LIBM
+
   AC_REQUIRE([gl_CHECK_LDEXP_NO_LIBM])
   LDEXP_LIBM=
   if test $gl_cv_func_ldexp_no_libm = no; then
@@ -30,6 +33,20 @@ AC_DEFUN([gl_FUNC_LDEXP]
       LDEXP_LIBM=-lm
     fi
   fi
+
+  save_LIBS="$LIBS"
+  LIBS="$LIBS $LDEXP_LIBM"
+  gl_FUNC_LDEXP_WORKS
+  LIBS="$save_LIBS"
+  case "$gl_cv_func_ldexp_works" in
+    *yes) ;;
+    *) REPLACE_LDEXP=1 ;;
+  esac
+
+  if test $REPLACE_LDEXP = 1; then
+    dnl Find libraries needed to link lib/ldexp.c.
+    LDEXP_LIBM="$ISNAND_LIBM"
+  fi
   AC_SUBST([LDEXP_LIBM])
 ])
 
@@ -52,3 +69,36 @@ AC_DEFUN([gl_CHECK_LDEXP_NO_LIBM]
         [gl_cv_func_ldexp_no_libm=no])
     ])
 ])
+
+dnl Test whether ldexp() works (this fails on OpenBSD 7.3/mips64).
+AC_DEFUN([gl_FUNC_LDEXP_WORKS],
+[
+  AC_REQUIRE([AC_PROG_CC])
+  AC_REQUIRE([AC_CANONICAL_HOST]) dnl for cross-compiles
+  AC_CACHE_CHECK([whether ldexp works], [gl_cv_func_ldexp_works],
+    [
+      AC_RUN_IFELSE(
+        [AC_LANG_SOURCE([[
+#include <math.h>
+int main()
+{
+  int result = 0;
+  {
+    volatile double x = 1.9269695883136991774e-308;
+    volatile double y = ldexp (x, 0);
+    if (y != x)
+      result |= 1;
+  }
+  return result;
+}]])],
+        [gl_cv_func_ldexp_works=yes],
+        [gl_cv_func_ldexp_works=no],
+        [case "$host_os" in
+           openbsd*)          gl_cv_func_ldexp_works="guessing no" ;;
+                              # Guess yes on native Windows.
+           mingw* | windows*) gl_cv_func_ldexp_works="guessing yes" ;;
+           *)                 gl_cv_func_ldexp_works="guessing yes" ;;
+         esac
+        ])
+    ])
+])
diff --git a/m4/math_h.m4 b/m4/math_h.m4
index d2e90ff1eb..c214f8efa8 100644
--- a/m4/math_h.m4
+++ b/m4/math_h.m4
@@ -1,4 +1,4 @@
-# math_h.m4 serial 125
+# math_h.m4 serial 126
 dnl Copyright (C) 2007-2023 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,6 +125,7 @@ AC_DEFUN([gl_MATH_H_REQUIRE_DEFAULTS]
     gl_MODULE_INDICATOR_INIT_VARIABLE([GNULIB_ISNANF])
     gl_MODULE_INDICATOR_INIT_VARIABLE([GNULIB_ISNAND])
     gl_MODULE_INDICATOR_INIT_VARIABLE([GNULIB_ISNANL])
+    gl_MODULE_INDICATOR_INIT_VARIABLE([GNULIB_LDEXP])
     gl_MODULE_INDICATOR_INIT_VARIABLE([GNULIB_LDEXPF])
     gl_MODULE_INDICATOR_INIT_VARIABLE([GNULIB_LDEXPL])
     gl_MODULE_INDICATOR_INIT_VARIABLE([GNULIB_LOG])
@@ -319,6 +320,7 @@ AC_DEFUN([gl_MATH_H_DEFAULTS]
   REPLACE_ISFINITE=0;               AC_SUBST([REPLACE_ISFINITE])
   REPLACE_ISINF=0;                  AC_SUBST([REPLACE_ISINF])
   REPLACE_ISNAN=0;                  AC_SUBST([REPLACE_ISNAN])
+  REPLACE_LDEXP=0;                  AC_SUBST([REPLACE_LDEXP])
   REPLACE_LDEXPL=0;                 AC_SUBST([REPLACE_LDEXPL])
   REPLACE_LOG=0;                    AC_SUBST([REPLACE_LOG])
   REPLACE_LOGF=0;                   AC_SUBST([REPLACE_LOGF])
diff --git a/modules/ldexp b/modules/ldexp
index 4e6ab1bb4d..2f55db6416 100644
--- a/modules/ldexp
+++ b/modules/ldexp
@@ -2,14 +2,22 @@ Description:
 ldexp() function: multiply a 'double' by a power of 2.
 
 Files:
+lib/ldexp.c
 m4/ldexp.m4
 
 Depends-on:
+math
+isnand          [test $REPLACE_LDEXP = 1]
 
 configure.ac:
 gl_FUNC_LDEXP
+gl_CONDITIONAL([GL_COND_OBJ_LDEXP], [test $REPLACE_LDEXP = 1])
+gl_MATH_MODULE_INDICATOR([ldexp])
 
 Makefile.am:
+if GL_COND_OBJ_LDEXP
+lib_SOURCES += ldexp.c
+endif
 
 Include:
 <math.h>
diff --git a/modules/ldexpl b/modules/ldexpl
index 5eae816b38..515ce713cd 100644
--- a/modules/ldexpl
+++ b/modules/ldexpl
@@ -3,6 +3,7 @@ ldexpl() function: multiply a 'long double' by a power of 2.
 
 Files:
 lib/ldexpl.c
+lib/ldexp.c
 m4/ldexpl.m4
 
 Depends-on:
diff --git a/modules/math b/modules/math
index f34fc8a978..f107141533 100644
--- a/modules/math
+++ b/modules/math
@@ -89,6 +89,7 @@ math.h: math.in.h $(top_builddir)/config.status $(CXXDEFS_H) 
$(ARG_NONNULL_H) $(
              -e 's/@''GNULIB_ISNANF''@/$(GNULIB_ISNANF)/g' \
              -e 's/@''GNULIB_ISNAND''@/$(GNULIB_ISNAND)/g' \
              -e 's/@''GNULIB_ISNANL''@/$(GNULIB_ISNANL)/g' \
+             -e 's/@''GNULIB_LDEXP''@/$(GNULIB_LDEXP)/g' \
              -e 's/@''GNULIB_LDEXPF''@/$(GNULIB_LDEXPF)/g' \
              -e 's/@''GNULIB_LDEXPL''@/$(GNULIB_LDEXPL)/g' \
              -e 's/@''GNULIB_LOG''@/$(GNULIB_LOG)/g' \
@@ -283,6 +284,7 @@ math.h: math.in.h $(top_builddir)/config.status 
$(CXXDEFS_H) $(ARG_NONNULL_H) $(
              -e 's|@''REPLACE_ITOLD''@|$(REPLACE_ITOLD)|g' \
              < $@-t4 > $@-t5
        $(AM_V_at)sed \
+             -e 's|@''REPLACE_LDEXP''@|$(REPLACE_LDEXP)|g' \
              -e 's|@''REPLACE_LDEXPL''@|$(REPLACE_LDEXPL)|g' \
              -e 's|@''REPLACE_LOG''@|$(REPLACE_LOG)|g' \
              -e 's|@''REPLACE_LOGF''@|$(REPLACE_LOGF)|g' \