> On OpenBSD 5.1/SPARC64, the sqrtl() function returns values that have a
> relative error > 50000000 ulp. Obviously unusable.
Likewise for the hypotl() function. Witness:
#include <math.h>
#include <stdio.h>
#include <float.h>
volatile long double x;
volatile long double y;
volatile long double z;
int main ()
{
x = 2.5541394760659556563446062497337725156L;
y = 7.7893454113437840832487794525518765265L;
z = hypotl (x, y);
long double err = (z * z - (x * x + y * y));
int i;
for (i = 1; i <= LDBL_MANT_DIG; i++)
err = err * 2.0L;
printf ("%.36Lf\n%.36Lf\n%.36Lf\n%Lg\n", x, y, z, err);
return 0;
}
prints
2.554139476065955656344606249733772516
7.789345411343784083248779452551876527
8.197409981233154068081014187198178182
-1.66122e+09
So this is an error of ca. 200 million ulps! Here's the workaround:
2012-03-13 Bruno Haible <br...@clisp.org>
hypotl: Bypass broken implementation in OpenBSD 5.1/SPARC.
* m4/hypotl.m4 (gl_FUNC_HYPOTL_WORKS): New macro.
(gl_FUNC_HYPOTL): Invoke it. If the function does not work, set
REPLACE_HYPOTL to 1.
* doc/posix-functions/hypotl.texi: Mention the OpenBSD 5.1/SPARC bug.
--- doc/posix-functions/hypotl.texi.orig Wed Mar 14 03:23:54 2012
+++ doc/posix-functions/hypotl.texi Wed Mar 14 02:58:39 2012
@@ -11,6 +11,9 @@
@item
This function is missing on some platforms:
FreeBSD 6.0, NetBSD 5.0, OpenBSD 3.8, Minix 3.1.8, AIX 5.1, HP-UX 11, IRIX
6.5, Solaris 9, Cygwin, MSVC 9, Interix 3.5, BeOS.
+@item
+This function produces very imprecise results on some platforms:
+OpenBSD 5.1/SPARC.
@end itemize
Portability problems fixed by Gnulib module @code{hypotl-ieee}:
--- m4/hypotl.m4.orig Wed Mar 14 03:23:54 2012
+++ m4/hypotl.m4 Wed Mar 14 03:16:44 2012
@@ -1,4 +1,4 @@
-# hypotl.m4 serial 3
+# hypotl.m4 serial 4
dnl Copyright (C) 2012 Free Software Foundation, Inc.
dnl This file is free software; the Free Software Foundation
dnl gives unlimited permission to copy and/or distribute it,
@@ -21,6 +21,16 @@
LIBS="$save_LIBS"
if test $ac_cv_func_hypotl = yes; then
HYPOTL_LIBM="$HYPOT_LIBM"
+
+ save_LIBS="$LIBS"
+ LIBS="$LIBS $HYPOTL_LIBM"
+ gl_FUNC_HYPOTL_WORKS
+ LIBS="$save_LIBS"
+ case "$gl_cv_func_hypotl_works" in
+ *yes) ;;
+ *) REPLACE_HYPOTL=1 ;;
+ esac
+
m4_ifdef([gl_FUNC_HYPOTL_IEEE], [
if test $gl_hypotl_required = ieee && test $REPLACE_HYPOTL = 0; then
AC_CACHE_CHECK([whether hypotl works according to ISO C 99 with IEC
60559],
@@ -105,3 +115,55 @@
fi
AC_SUBST([HYPOTL_LIBM])
])
+
+dnl Test whether hypotl() works.
+dnl On OpenBSD 5.1/SPARC,
+dnl hypotl (2.5541394760659556563446062497337725156L,
7.7893454113437840832487794525518765265L)
+dnl has rounding errors that eat up the last 8 to 9 decimal digits.
+AC_DEFUN([gl_FUNC_HYPOTL_WORKS],
+[
+ AC_REQUIRE([AC_PROG_CC])
+ AC_REQUIRE([AC_CANONICAL_HOST]) dnl for cross-compiles
+ AC_CACHE_CHECK([whether hypotl works], [gl_cv_func_hypotl_works],
+ [
+ AC_RUN_IFELSE(
+ [AC_LANG_SOURCE([[
+#include <float.h>
+#include <math.h>
+static long double
+my_ldexpl (long double x, int d)
+{
+ for (; d > 0; d--)
+ x *= 2.0L;
+ for (; d < 0; d++)
+ x *= 0.5L;
+ return x;
+}
+volatile long double x;
+volatile long double y;
+volatile long double z;
+int main ()
+{
+ long double err;
+
+ x = 2.5541394760659556563446062497337725156L;
+ y = 7.7893454113437840832487794525518765265L;
+ z = hypotl (x, y);
+ err = z * z - (x * x + y * y);
+ err = my_ldexpl (err, LDBL_MANT_DIG);
+ if (err < 0)
+ err = - err;
+ if (err > 1000.0L)
+ return 1;
+ return 0;
+}
+]])],
+ [gl_cv_func_hypotl_works=yes],
+ [gl_cv_func_hypotl_works=no],
+ [case "$host_os" in
+ openbsd*) gl_cv_func_hypotl_works="guessing no";;
+ *) gl_cv_func_hypotl_works="guessing yes";;
+ esac
+ ])
+ ])
+])
