2 dnl Copyright (C) 2011-2012 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.
7 AC_DEFUN([gl_FUNC_FMA],
9 AC_REQUIRE([gl_MATH_H_DEFAULTS])
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
16 case "$gl_cv_func_fma_works" in
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])
28 dnl Append $FREXP_LIBM to FMA_LIBM, avoiding gratuitous duplicates.
31 *) FMA_LIBM="$FMA_LIBM $FREXP_LIBM" ;;
33 dnl Append $LDEXP_LIBM to FMA_LIBM, avoiding gratuitous duplicates.
36 *) FMA_LIBM="$FMA_LIBM $LDEXP_LIBM" ;;
38 dnl Append $FEGETROUND_LIBM to FMA_LIBM, avoiding gratuitous duplicates.
40 *" $FEGETROUND_LIBM "*) ;;
41 *) FMA_LIBM="$FMA_LIBM $FEGETROUND_LIBM" ;;
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],
50 AC_REQUIRE([AC_PROG_CC])
51 AC_REQUIRE([AC_CANONICAL_HOST]) dnl for cross-compiles
52 AC_REQUIRE([gl_FUNC_LDEXP])
54 LIBS="$LIBS $FMA_LIBM $LDEXP_LIBM"
55 AC_CACHE_CHECK([whether fma works], [gl_cv_func_fma_works],
65 /* These tests fail with glibc 2.11.3 on x86_64. */
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)
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)
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)
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)
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)
131 if ((DBL_MANT_DIG % 2) == 1)
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)
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))
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"
165 # Prerequisites of lib/fma.c.
166 AC_DEFUN([gl_PREREQ_FMA], [:])