m4/fma.m4

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