New module 'fma'.

* lib/math.in.h (fma): New declaration.
* lib/fma.c: New file.
* m4/fma.m4: New file.
* m4/fegetround.m4: New file.
* m4/math_h.m4 (gl_MATH_H): Test whethern fma is declared.
(gl_MATH_H_DEFAULTS): Initialize GNULIB_FMA, HAVE_FMA, REPLACE_FMA.
* modules/math (Makefile.am): Substitute GNULIB_FMA, HAVE_FMA,
REPLACE_FMA.
* modules/fma: New file.
* doc/posix-functions/fma.texi: Mention the new module and the various
bugs.

3 files changed, 191 insertions, 2 deletions

diff --git a/m4/fegetround.m4 b/m4/fegetround.m4
new file mode 100644
index 0000000000..37db5c8858
--- /dev/null
+++ b/m4/fegetround.m4
@@ -0,0 +1,20 @@
+# fegetround.m4 serial 1
+dnl Copyright (C) 2011 Free Software Foundation, Inc.
+dnl This file is free software; the Free Software Foundation
+dnl gives unlimited permission to copy and/or distribute it,
+dnl with or without modifications, as long as this notice is preserved.
+
+AC_DEFUN([gl_FUNC_FEGETROUND],
+[
+  dnl Determine FEGETROUND_LIBM.
+  gl_MATHFUNC([fegetround], [int], [(void)], [#include <fenv.h>])
+  if test $gl_cv_func_fegetround_no_libm = no \
+     && test $gl_cv_func_fegetround_in_libm = no; then
+    HAVE_FEGETROUND=0
+  else
+    HAVE_FEGETROUND=1
+    AC_DEFINE([HAVE_FEGETROUND], [1],
+      [Define to 1 if you have the 'fegetround' function.])
+  fi
+  AC_SUBST([FEGETROUND_LIBM])
+])
diff --git a/m4/fma.m4 b/m4/fma.m4
new file mode 100644
index 0000000000..3ce3f9a05c
--- /dev/null
+++ b/m4/fma.m4
@@ -0,0 +1,166 @@
+# fma.m4 serial 1
+dnl Copyright (C) 2011 Free Software Foundation, Inc.
+dnl This file is free software; the Free Software Foundation
+dnl gives unlimited permission to copy and/or distribute it,
+dnl with or without modifications, as long as this notice is preserved.
+
+AC_DEFUN([gl_FUNC_FMA],
+[
+  AC_REQUIRE([gl_MATH_H_DEFAULTS])
+
+  dnl Determine FMA_LIBM.
+  gl_MATHFUNC([fma], [double], [(double, double, double)])
+  if test $gl_cv_func_fma_no_libm = yes \
+     || test $gl_cv_func_fma_in_libm = yes; then
+    gl_FUNC_FMA_WORKS
+    case "$gl_cv_func_fma_works" in
+      *no) REPLACE_FMA=1 ;;
+    esac
+  else
+    HAVE_FMA=0
+  fi
+  if test $HAVE_FMA = 0 || test $REPLACE_FMA = 1; then
+    dnl Find libraries needed to link lib/fmal.c.
+    AC_REQUIRE([gl_FUNC_FREXP])
+    AC_REQUIRE([gl_FUNC_LDEXP])
+    AC_REQUIRE([gl_FUNC_FEGETROUND])
+    FMA_LIBM=
+    dnl Append $FREXP_LIBM to FMA_LIBM, avoiding gratuitous duplicates.
+    case " $FMA_LIBM " in
+      *" $FREXP_LIBM "*) ;;
+      *) FMA_LIBM="$FMA_LIBM $FREXP_LIBM" ;;
+    esac
+    dnl Append $LDEXP_LIBM to FMA_LIBM, avoiding gratuitous duplicates.
+    case " $FMA_LIBM " in
+      *" $LDEXP_LIBM "*) ;;
+      *) FMA_LIBM="$FMA_LIBM $LDEXP_LIBM" ;;
+    esac
+    dnl Append $FEGETROUND_LIBM to FMA_LIBM, avoiding gratuitous duplicates.
+    case " $FMA_LIBM " in
+      *" $FEGETROUND_LIBM "*) ;;
+      *) FMA_LIBM="$FMA_LIBM $FEGETROUND_LIBM" ;;
+    esac
+  fi
+  AC_SUBST([FMA_LIBM])
+])
+
+dnl Test whether fma() has any of the 7 known bugs of glibc 2.11.3 on x86_64.
+AC_DEFUN([gl_FUNC_FMA_WORKS],
+[
+  AC_REQUIRE([AC_PROG_CC])
+  AC_REQUIRE([AC_CANONICAL_HOST]) dnl for cross-compiles
+  AC_REQUIRE([gl_FUNC_LDEXP])
+  save_LIBS="$LIBS"
+  LIBS="$LIBS $FMA_LIBM $LDEXP_LIBM"
+  AC_CACHE_CHECK([whether fma works], [gl_cv_func_fma_works],
+    [
+      AC_RUN_IFELSE(
+        [AC_LANG_SOURCE([[
+#include <float.h>
+#include <math.h>
+double p0 = 0.0;
+int main()
+{
+  int failed_tests = 0;
+  /* These tests fail with glibc 2.11.3 on x86_64.  */
+  {
+    volatile double x = 1.5; /* 3 * 2^-1 */
+    volatile double y = x;
+    volatile double z = ldexp (1.0, DBL_MANT_DIG + 1); /* 2^54 */
+    /* x * y + z with infinite precision: 2^54 + 9 * 2^-2.
+       Lies between (2^52 + 0) * 2^2 and (2^52 + 1) * 2^2
+       and is closer to (2^52 + 1) * 2^2, therefore the rounding
+       must round up and produce (2^52 + 1) * 2^2.  */
+    volatile double expected = z + 4.0;
+    volatile double result = fma (x, y, z);
+    if (result != expected)
+      failed_tests |= 1;
+  }
+  {
+    volatile double x = 1.25; /* 2^0 + 2^-2 */
+    volatile double y = - x;
+    volatile double z = ldexp (1.0, DBL_MANT_DIG + 1); /* 2^54 */
+    /* x * y + z with infinite precision: 2^54 - 2^0 - 2^-1 - 2^-4.
+       Lies between (2^53 - 1) * 2^1 and 2^53 * 2^1
+       and is closer to (2^53 - 1) * 2^1, therefore the rounding
+       must round down and produce (2^53 - 1) * 2^1.  */
+    volatile double expected = (ldexp (1.0, DBL_MANT_DIG) - 1.0) * 2.0;
+    volatile double result = fma (x, y, z);
+    if (result != expected)
+      failed_tests |= 2;
+  }
+  {
+    volatile double x = 1.0 + ldexp (1.0, 1 - DBL_MANT_DIG); /* 2^0 + 2^-52 */
+    volatile double y = x;
+    volatile double z = 4.0; /* 2^2 */
+    /* x * y + z with infinite precision: 2^2 + 2^0 + 2^-51 + 2^-104.
+       Lies between (2^52 + 2^50) * 2^-50 and (2^52 + 2^50 + 1) * 2^-50
+       and is closer to (2^52 + 2^50 + 1) * 2^-50, therefore the rounding
+       must round up and produce (2^52 + 2^50 + 1) * 2^-50.  */
+    volatile double expected = 4.0 + 1.0 + ldexp (1.0, 3 - DBL_MANT_DIG);
+    volatile double result = fma (x, y, z);
+    if (result != expected)
+      failed_tests |= 4;
+  }
+  {
+    volatile double x = 1.0 + ldexp (1.0, 1 - DBL_MANT_DIG); /* 2^0 + 2^-52 */
+    volatile double y = - x;
+    volatile double z = 8.0; /* 2^3 */
+    /* x * y + z with infinite precision: 2^2 + 2^1 + 2^0 - 2^-51 - 2^-104.
+       Lies between (2^52 + 2^51 + 2^50 - 1) * 2^-50 and
+       (2^52 + 2^51 + 2^50) * 2^-50 and is closer to
+       (2^52 + 2^51 + 2^50 - 1) * 2^-50, therefore the rounding
+       must round down and produce (2^52 + 2^51 + 2^50 - 1) * 2^-50.  */
+    volatile double expected = 7.0 - ldexp (1.0, 3 - DBL_MANT_DIG);
+    volatile double result = fma (x, y, z);
+    if (result != expected)
+      failed_tests |= 8;
+  }
+  {
+    volatile double x = 1.25; /* 2^0 + 2^-2 */
+    volatile double y = - 0.75; /* - 2^0 + 2^-2 */
+    volatile double z = ldexp (1.0, DBL_MANT_DIG); /* 2^53 */
+    /* x * y + z with infinite precision: 2^53 - 2^0 + 2^-4.
+       Lies between (2^53 - 2^0) and 2^53 and is closer to (2^53 - 2^0),
+       therefore the rounding must round down and produce (2^53 - 2^0).  */
+    volatile double expected = ldexp (1.0, DBL_MANT_DIG) - 1.0;
+    volatile double result = fma (x, y, z);
+    if (result != expected)
+      failed_tests |= 16;
+  }
+  if ((DBL_MANT_DIG % 2) == 1)
+    {
+      volatile double x = 1.0 + ldexp (1.0, - (DBL_MANT_DIG + 1) / 2); /* 2^0 + 2^-27 */
+      volatile double y = 1.0 - ldexp (1.0, - (DBL_MANT_DIG + 1) / 2); /* 2^0 - 2^-27 */
+      volatile double z = - ldexp (1.0, DBL_MIN_EXP - DBL_MANT_DIG); /* - 2^-1074 */
+      /* x * y + z with infinite precision: 2^0 - 2^-54 - 2^-1074.
+         Lies between (2^53 - 1) * 2^-53 and 2^53 * 2^-53 and is closer to
+         (2^53 - 1) * 2^-53, therefore the rounding must round down and
+         produce (2^53 - 1) * 2^-53.  */
+      volatile double expected = 1.0 - ldexp (1.0, - DBL_MANT_DIG);
+      volatile double result = fma (x, y, z);
+      if (result != expected)
+        failed_tests |= 32;
+    }
+  {
+    double minus_inf = -1.0 / p0;
+    volatile double x = ldexp (1.0, DBL_MAX_EXP - 1);
+    volatile double y = ldexp (1.0, DBL_MAX_EXP - 1);
+    volatile double z = minus_inf;
+    volatile double result = fma (x, y, z);
+    if (!(result == minus_inf))
+      failed_tests |= 64;
+  }
+  return failed_tests;
+}]])],
+        [gl_cv_func_fma_works=yes],
+        [gl_cv_func_fma_works=no],
+        [dnl Guess no, even on glibc systems.
+         gl_cv_func_fma_works="guessing no"
+        ])
+    ])
+  LIBS="$save_LIBS"
+])
+
+# Prerequisites of lib/fma.c.
+AC_DEFUN([gl_PREREQ_FMA], [:])
diff --git a/m4/math_h.m4 b/m4/math_h.m4
index e5a2892481..878aaff3a1 100644
--- a/m4/math_h.m4
+++ b/m4/math_h.m4
@@ -1,4 +1,4 @@
-# math_h.m4 serial 53
+# math_h.m4 serial 54
 dnl Copyright (C) 2007-2011 Free Software Foundation, Inc.
 dnl This file is free software; the Free Software Foundation
 dnl gives unlimited permission to copy and/or distribute it,
@@ -41,7 +41,7 @@ AC_DEFUN([gl_MATH_H],
   gl_WARN_ON_USE_PREPARE([[#include <math.h>]],
     [acosf acosl asinf asinl atanf atanl
      ceilf ceill copysign copysignf copysignl cosf cosl coshf
-     expf expl fabsf floorf floorl fmodf frexpf frexpl
+     expf expl fabsf floorf floorl fma fmodf frexpf frexpl
      ldexpf ldexpl logb logf logl log10f modff powf
      rint rintf rintl round roundf roundl sinf sinl sinhf sqrtf sqrtl
      tanf tanl tanhf trunc truncf truncl])
@@ -80,6 +80,7 @@ AC_DEFUN([gl_MATH_H_DEFAULTS],
   GNULIB_FLOOR=0;     AC_SUBST([GNULIB_FLOOR])
   GNULIB_FLOORF=0;    AC_SUBST([GNULIB_FLOORF])
   GNULIB_FLOORL=0;    AC_SUBST([GNULIB_FLOORL])
+  GNULIB_FMA=0;       AC_SUBST([GNULIB_FMA])
   GNULIB_FMODF=0;     AC_SUBST([GNULIB_FMODF])
   GNULIB_FREXPF=0;    AC_SUBST([GNULIB_FREXPF])
   GNULIB_FREXP=0;     AC_SUBST([GNULIB_FREXP])
@@ -133,6 +134,7 @@ AC_DEFUN([gl_MATH_H_DEFAULTS],
   HAVE_EXPF=1;                 AC_SUBST([HAVE_EXPF])
   HAVE_EXPL=1;                 AC_SUBST([HAVE_EXPL])
   HAVE_FABSF=1;                AC_SUBST([HAVE_FABSF])
+  HAVE_FMA=1;                  AC_SUBST([HAVE_FMA])
   HAVE_FMODF=1;                AC_SUBST([HAVE_FMODF])
   HAVE_FREXPF=1;               AC_SUBST([HAVE_FREXPF])
   HAVE_ISNANF=1;               AC_SUBST([HAVE_ISNANF])
@@ -183,6 +185,7 @@ AC_DEFUN([gl_MATH_H_DEFAULTS],
   REPLACE_FLOOR=0;             AC_SUBST([REPLACE_FLOOR])
   REPLACE_FLOORF=0;            AC_SUBST([REPLACE_FLOORF])
   REPLACE_FLOORL=0;            AC_SUBST([REPLACE_FLOORL])
+  REPLACE_FMA=0;               AC_SUBST([REPLACE_FMA])
   REPLACE_FREXPF=0;            AC_SUBST([REPLACE_FREXPF])
   REPLACE_FREXP=0;             AC_SUBST([REPLACE_FREXP])
   REPLACE_FREXPL=0;            AC_SUBST([REPLACE_FREXPL])