m4/fma.m4

   1 # fma.m4 serial 2
   2 dnl Copyright (C) 2011-2014 Free Software Foundation, Inc.
   3 dnl This file is free software; the Free Software Foundation
   4 dnl gives unlimited permission to copy and/or distribute it,
   5 dnl with or without modifications, as long as this notice is preserved.
   6
   7 AC_DEFUN([gl_FUNC_FMA],
   8 [
   9   AC_REQUIRE([gl_MATH_H_DEFAULTS])
  10
  11   dnl Determine FMA_LIBM.
  12   gl_MATHFUNC([fma], [double], [(double, double, double)],
  13     [extern
  14      #ifdef __cplusplus
  15      "C"
  16      #endif
  17      double fma (double, double, double);
  18     ])
  19   if test $gl_cv_func_fma_no_libm = yes \
  20      || test $gl_cv_func_fma_in_libm = yes; then
  21     dnl Also check whether it's declared.
  22     dnl IRIX 6.5 has fma() in libm but doesn't declare it in <math.h>,
  23     dnl and the function is buggy.
  24     AC_CHECK_DECL([fma], , [REPLACE_FMA=1], [[#include <math.h>]])
  25     if test $REPLACE_FMA = 0; then
  26       gl_FUNC_FMA_WORKS
  27       case "$gl_cv_func_fma_works" in
  28         *no) REPLACE_FMA=1 ;;
  29       esac
  30     fi
  31   else
  32     HAVE_FMA=0
  33   fi
  34   if test $HAVE_FMA = 0 || test $REPLACE_FMA = 1; then
  35     dnl Find libraries needed to link lib/fmal.c.
  36     AC_REQUIRE([gl_FUNC_FREXP])
  37     AC_REQUIRE([gl_FUNC_LDEXP])
  38     AC_REQUIRE([gl_FUNC_FEGETROUND])
  39     FMA_LIBM=
  40     dnl Append $FREXP_LIBM to FMA_LIBM, avoiding gratuitous duplicates.
  41     case " $FMA_LIBM " in
  42       *" $FREXP_LIBM "*) ;;
  43       *) FMA_LIBM="$FMA_LIBM $FREXP_LIBM" ;;
  44     esac
  45     dnl Append $LDEXP_LIBM to FMA_LIBM, avoiding gratuitous duplicates.
  46     case " $FMA_LIBM " in
  47       *" $LDEXP_LIBM "*) ;;
  48       *) FMA_LIBM="$FMA_LIBM $LDEXP_LIBM" ;;
  49     esac
  50     dnl Append $FEGETROUND_LIBM to FMA_LIBM, avoiding gratuitous duplicates.
  51     case " $FMA_LIBM " in
  52       *" $FEGETROUND_LIBM "*) ;;
  53       *) FMA_LIBM="$FMA_LIBM $FEGETROUND_LIBM" ;;
  54     esac
  55   fi
  56   AC_SUBST([FMA_LIBM])
  57 ])
  58
  59 dnl Test whether fma() has any of the 7 known bugs of glibc 2.11.3 on x86_64.
  60 AC_DEFUN([gl_FUNC_FMA_WORKS],
  61 [
  62   AC_REQUIRE([AC_PROG_CC])
  63   AC_REQUIRE([AC_CANONICAL_HOST]) dnl for cross-compiles
  64   AC_REQUIRE([gl_FUNC_LDEXP])
  65   save_LIBS="$LIBS"
  66   LIBS="$LIBS $FMA_LIBM $LDEXP_LIBM"
  67   AC_CACHE_CHECK([whether fma works], [gl_cv_func_fma_works],
  68     [
  69       AC_RUN_IFELSE(
  70         [AC_LANG_SOURCE([[
  71 #include <float.h>
  72 #include <math.h>
  73 double p0 = 0.0;
  74 int main()
  75 {
  76   int failed_tests = 0;
  77   /* These tests fail with glibc 2.11.3 on x86_64.  */
  78   {
  79     volatile double x = 1.5; /* 3 * 2^-1 */
  80     volatile double y = x;
  81     volatile double z = ldexp (1.0, DBL_MANT_DIG + 1); /* 2^54 */
  82     /* x * y + z with infinite precision: 2^54 + 9 * 2^-2.
  83        Lies between (2^52 + 0) * 2^2 and (2^52 + 1) * 2^2
  84        and is closer to (2^52 + 1) * 2^2, therefore the rounding
  85        must round up and produce (2^52 + 1) * 2^2.  */
  86     volatile double expected = z + 4.0;
  87     volatile double result = fma (x, y, z);
  88     if (result != expected)
  89       failed_tests |= 1;
  90   }
  91   {
  92     volatile double x = 1.25; /* 2^0 + 2^-2 */
  93     volatile double y = - x;
  94     volatile double z = ldexp (1.0, DBL_MANT_DIG + 1); /* 2^54 */
  95     /* x * y + z with infinite precision: 2^54 - 2^0 - 2^-1 - 2^-4.
  96        Lies between (2^53 - 1) * 2^1 and 2^53 * 2^1
  97        and is closer to (2^53 - 1) * 2^1, therefore the rounding
  98        must round down and produce (2^53 - 1) * 2^1.  */
  99     volatile double expected = (ldexp (1.0, DBL_MANT_DIG) - 1.0) * 2.0;
 100     volatile double result = fma (x, y, z);
 101     if (result != expected)
 102       failed_tests |= 2;
 103   }
 104   {
 105     volatile double x = 1.0 + ldexp (1.0, 1 - DBL_MANT_DIG); /* 2^0 + 2^-52 */
 106     volatile double y = x;
 107     volatile double z = 4.0; /* 2^2 */
 108     /* x * y + z with infinite precision: 2^2 + 2^0 + 2^-51 + 2^-104.
 109        Lies between (2^52 + 2^50) * 2^-50 and (2^52 + 2^50 + 1) * 2^-50
 110        and is closer to (2^52 + 2^50 + 1) * 2^-50, therefore the rounding
 111        must round up and produce (2^52 + 2^50 + 1) * 2^-50.  */
 112     volatile double expected = 4.0 + 1.0 + ldexp (1.0, 3 - DBL_MANT_DIG);
 113     volatile double result = fma (x, y, z);
 114     if (result != expected)
 115       failed_tests |= 4;
 116   }
 117   {
 118     volatile double x = 1.0 + ldexp (1.0, 1 - DBL_MANT_DIG); /* 2^0 + 2^-52 */
 119     volatile double y = - x;
 120     volatile double z = 8.0; /* 2^3 */
 121     /* x * y + z with infinite precision: 2^2 + 2^1 + 2^0 - 2^-51 - 2^-104.
 122        Lies between (2^52 + 2^51 + 2^50 - 1) * 2^-50 and
 123        (2^52 + 2^51 + 2^50) * 2^-50 and is closer to
 124        (2^52 + 2^51 + 2^50 - 1) * 2^-50, therefore the rounding
 125        must round down and produce (2^52 + 2^51 + 2^50 - 1) * 2^-50.  */
 126     volatile double expected = 7.0 - ldexp (1.0, 3 - DBL_MANT_DIG);
 127     volatile double result = fma (x, y, z);
 128     if (result != expected)
 129       failed_tests |= 8;
 130   }
 131   {
 132     volatile double x = 1.25; /* 2^0 + 2^-2 */
 133     volatile double y = - 0.75; /* - 2^0 + 2^-2 */
 134     volatile double z = ldexp (1.0, DBL_MANT_DIG); /* 2^53 */
 135     /* x * y + z with infinite precision: 2^53 - 2^0 + 2^-4.
 136        Lies between (2^53 - 2^0) and 2^53 and is closer to (2^53 - 2^0),
 137        therefore the rounding must round down and produce (2^53 - 2^0).  */
 138     volatile double expected = ldexp (1.0, DBL_MANT_DIG) - 1.0;
 139     volatile double result = fma (x, y, z);
 140     if (result != expected)
 141       failed_tests |= 16;
 142   }
 143   if ((DBL_MANT_DIG % 2) == 1)
 144     {
 145       volatile double x = 1.0 + ldexp (1.0, - (DBL_MANT_DIG + 1) / 2); /* 2^0 + 2^-27 */
 146       volatile double y = 1.0 - ldexp (1.0, - (DBL_MANT_DIG + 1) / 2); /* 2^0 - 2^-27 */
 147       volatile double z = - ldexp (1.0, DBL_MIN_EXP - DBL_MANT_DIG); /* - 2^-1074 */
 148       /* x * y + z with infinite precision: 2^0 - 2^-54 - 2^-1074.
 149          Lies between (2^53 - 1) * 2^-53 and 2^53 * 2^-53 and is closer to
 150          (2^53 - 1) * 2^-53, therefore the rounding must round down and
 151          produce (2^53 - 1) * 2^-53.  */
 152       volatile double expected = 1.0 - ldexp (1.0, - DBL_MANT_DIG);
 153       volatile double result = fma (x, y, z);
 154       if (result != expected)
 155         failed_tests |= 32;
 156     }
 157   {
 158     double minus_inf = -1.0 / p0;
 159     volatile double x = ldexp (1.0, DBL_MAX_EXP - 1);
 160     volatile double y = ldexp (1.0, DBL_MAX_EXP - 1);
 161     volatile double z = minus_inf;
 162     volatile double result = fma (x, y, z);
 163     if (!(result == minus_inf))
 164       failed_tests |= 64;
 165   }
 166   return failed_tests;
 167 }]])],
 168         [gl_cv_func_fma_works=yes],
 169         [gl_cv_func_fma_works=no],
 170         [dnl Guess no, even on glibc systems.
 171          gl_cv_func_fma_works="guessing no"
 172         ])
 173     ])
 174   LIBS="$save_LIBS"
 175 ])
 176
 177 # Prerequisites of lib/fma.c.
 178 AC_DEFUN([gl_PREREQ_FMA], [:])