m4/fma.m4

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