m4/fmal.m4

   1 # fmal.m4 serial 7
   2 dnl Copyright (C) 2011-2020 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_FMAL],
   8 [
   9   AC_REQUIRE([gl_MATH_H_DEFAULTS])
  10   AC_REQUIRE([gl_LONG_DOUBLE_VS_DOUBLE])
  11
  12   dnl Persuade glibc <math.h> to declare fmal().
  13   AC_REQUIRE([gl_USE_SYSTEM_EXTENSIONS])
  14
  15   dnl Determine FMAL_LIBM.
  16   gl_MATHFUNC([fmal], [long double], [(long double, long double, long double)],
  17     [extern
  18      #ifdef __cplusplus
  19      "C"
  20      #endif
  21      long double fmal (long double, long double, long double);
  22     ])
  23   if test $gl_cv_func_fmal_no_libm = yes \
  24      || test $gl_cv_func_fmal_in_libm = yes; then
  25     dnl Also check whether it's declared.
  26     dnl IRIX 6.5 has fmal() in libm but doesn't declare it in <math.h>,
  27     dnl and the function is buggy.
  28     AC_CHECK_DECL([fmal], , [REPLACE_FMAL=1], [[#include <math.h>]])
  29     if test $REPLACE_FMAL = 0; then
  30       gl_FUNC_FMAL_WORKS
  31       case "$gl_cv_func_fmal_works" in
  32         *no) REPLACE_FMAL=1 ;;
  33       esac
  34     fi
  35   else
  36     HAVE_FMAL=0
  37   fi
  38   if test $HAVE_FMAL = 0 || test $REPLACE_FMAL = 1; then
  39     dnl Find libraries needed to link lib/fmal.c.
  40     if test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 1; then
  41       AC_REQUIRE([gl_FUNC_FMA])
  42       FMAL_LIBM="$FMA_LIBM"
  43     else
  44       AC_REQUIRE([gl_FUNC_FREXPL])
  45       AC_REQUIRE([gl_FUNC_LDEXPL])
  46       AC_REQUIRE([gl_FUNC_FEGETROUND])
  47       FMAL_LIBM=
  48       dnl Append $FREXPL_LIBM to FMAL_LIBM, avoiding gratuitous duplicates.
  49       case " $FMAL_LIBM " in
  50         *" $FREXPL_LIBM "*) ;;
  51         *) FMAL_LIBM="$FMAL_LIBM $FREXPL_LIBM" ;;
  52       esac
  53       dnl Append $LDEXPL_LIBM to FMAL_LIBM, avoiding gratuitous duplicates.
  54       case " $FMAL_LIBM " in
  55         *" $LDEXPL_LIBM "*) ;;
  56         *) FMAL_LIBM="$FMAL_LIBM $LDEXPL_LIBM" ;;
  57       esac
  58       dnl Append $FEGETROUND_LIBM to FMAL_LIBM, avoiding gratuitous duplicates.
  59       case " $FMAL_LIBM " in
  60         *" $FEGETROUND_LIBM "*) ;;
  61         *) FMAL_LIBM="$FMAL_LIBM $FEGETROUND_LIBM" ;;
  62       esac
  63     fi
  64   fi
  65   AC_SUBST([FMAL_LIBM])
  66 ])
  67
  68 dnl Test whether fmal() has any of the 15 known bugs of glibc 2.11.3 on x86_64.
  69 AC_DEFUN([gl_FUNC_FMAL_WORKS],
  70 [
  71   AC_REQUIRE([AC_PROG_CC])
  72   AC_REQUIRE([AC_CANONICAL_HOST]) dnl for cross-compiles
  73   AC_REQUIRE([gl_FUNC_LDEXPL])
  74   save_LIBS="$LIBS"
  75   LIBS="$LIBS $FMAL_LIBM $LDEXPL_LIBM"
  76   AC_CACHE_CHECK([whether fmal works], [gl_cv_func_fmal_works],
  77     [
  78       AC_RUN_IFELSE(
  79         [AC_LANG_SOURCE([[
  80 #include <float.h>
  81 #include <math.h>
  82 /* Override the values of <float.h>, like done in float.in.h.  */
  83 #if defined __i386__ && (defined __BEOS__ || defined __OpenBSD__)
  84 # undef LDBL_MANT_DIG
  85 # define LDBL_MANT_DIG   64
  86 # undef LDBL_MIN_EXP
  87 # define LDBL_MIN_EXP    (-16381)
  88 # undef LDBL_MAX_EXP
  89 # define LDBL_MAX_EXP    16384
  90 #endif
  91 #if defined __i386__ && (defined __FreeBSD__ || defined __DragonFly__)
  92 # undef LDBL_MANT_DIG
  93 # define LDBL_MANT_DIG   64
  94 # undef LDBL_MIN_EXP
  95 # define LDBL_MIN_EXP    (-16381)
  96 # undef LDBL_MAX_EXP
  97 # define LDBL_MAX_EXP    16384
  98 #endif
  99 #if (defined _ARCH_PPC || defined _POWER) && defined _AIX && (LDBL_MANT_DIG == 106) && defined __GNUC__
 100 # undef LDBL_MIN_EXP
 101 # define LDBL_MIN_EXP DBL_MIN_EXP
 102 #endif
 103 #if defined __sgi && (LDBL_MANT_DIG >= 106)
 104 # undef LDBL_MANT_DIG
 105 # define LDBL_MANT_DIG 106
 106 # if defined __GNUC__
 107 #  undef LDBL_MIN_EXP
 108 #  define LDBL_MIN_EXP DBL_MIN_EXP
 109 # endif
 110 #endif
 111 long double p0 = 0.0L;
 112 int main()
 113 {
 114   int failed_tests = 0;
 115   /* This test fails on glibc 2.11 powerpc.  */
 116   {
 117     volatile long double x = 1.5L; /* 3 * 2^-1 */
 118     volatile long double y = x;
 119     volatile long double z = ldexpl (1.0L, LDBL_MANT_DIG + 1); /* 2^65 */
 120     /* x * y + z with infinite precision: 2^65 + 9 * 2^-2.
 121        Lies between (2^63 + 0) * 2^2 and (2^63 + 1) * 2^2
 122        and is closer to (2^63 + 1) * 2^2, therefore the rounding
 123        must round up and produce (2^63 + 1) * 2^2.  */
 124     volatile long double expected = z + 4.0L;
 125     volatile long double result = fmal (x, y, z);
 126     if (result != expected)
 127       failed_tests |= 1;
 128   }
 129   /* This test fails on glibc 2.11 powerpc.  */
 130   {
 131     volatile long double x = 1.25L; /* 2^0 + 2^-2 */
 132     volatile long double y = - x;
 133     volatile long double z = ldexpl (1.0L, LDBL_MANT_DIG + 1); /* 2^65 */
 134     /* x * y + z with infinite precision: 2^65 - 2^0 - 2^-1 - 2^-4.
 135        Lies between (2^64 - 1) * 2^1 and 2^64 * 2^1
 136        and is closer to (2^64 - 1) * 2^1, therefore the rounding
 137        must round down and produce (2^64 - 1) * 2^1.  */
 138     volatile long double expected = (ldexpl (1.0L, LDBL_MANT_DIG) - 1.0L) * 2.0L;
 139     volatile long double result = fmal (x, y, z);
 140     if (result != expected)
 141       failed_tests |= 1;
 142   }
 143   /* This test fails on glibc 2.11 x86,x86_64,powerpc, glibc 2.7 hppa,sparc,
 144      OSF/1 5.1, mingw.  */
 145   {
 146     volatile long double x = 1.0L + ldexpl (1.0L, 1 - LDBL_MANT_DIG); /* 2^0 + 2^-63 */
 147     volatile long double y = x;
 148     volatile long double z = 4.0L; /* 2^2 */
 149     /* x * y + z with infinite precision: 2^2 + 2^0 + 2^-62 + 2^-126.
 150        Lies between (2^63 + 2^61) * 2^-61 and (2^63 + 2^61 + 1) * 2^-61
 151        and is closer to (2^63 + 2^61 + 1) * 2^-61, therefore the rounding
 152        must round up and produce (2^63 + 2^61 + 1) * 2^-61.  */
 153     volatile long double expected = 4.0L + 1.0L + ldexpl (1.0L, 3 - LDBL_MANT_DIG);
 154     volatile long double result = fmal (x, y, z);
 155     if (result != expected)
 156       failed_tests |= 2;
 157   }
 158   /* This test fails on glibc 2.11 x86,x86_64,powerpc glibc 2.7 hppa,sparc,
 159      OSF/1 5.1, mingw.  */
 160   {
 161     volatile long double x = 1.0L + ldexpl (1.0L, 1 - LDBL_MANT_DIG); /* 2^0 + 2^-63 */
 162     volatile long double y = - x;
 163     volatile long double z = 8.0L; /* 2^3 */
 164     /* x * y + z with infinite precision: 2^2 + 2^1 + 2^0 - 2^-62 - 2^-126.
 165        Lies between (2^63 + 2^62 + 2^61 - 1) * 2^-61 and
 166        (2^63 + 2^62 + 2^61) * 2^-61 and is closer to
 167        (2^63 + 2^62 + 2^61 - 1) * 2^-61, therefore the rounding
 168        must round down and produce (2^63 + 2^62 + 2^61 - 1) * 2^-61.  */
 169     volatile long double expected = 7.0L - ldexpl (1.0L, 3 - LDBL_MANT_DIG);
 170     volatile long double result = fmal (x, y, z);
 171     if (result != expected)
 172       failed_tests |= 2;
 173   }
 174   /* This test fails on glibc 2.11 powerpc.  */
 175   {
 176     volatile long double x = 1.25L; /* 2^0 + 2^-2 */
 177     volatile long double y = - 0.75L; /* - 2^0 + 2^-2 */
 178     volatile long double z = ldexpl (1.0L, LDBL_MANT_DIG); /* 2^64 */
 179     /* x * y + z with infinite precision: 2^64 - 2^0 + 2^-4.
 180        Lies between (2^64 - 2^0) and 2^64 and is closer to (2^64 - 2^0),
 181        therefore the rounding must round down and produce (2^64 - 2^0).  */
 182     volatile long double expected = ldexpl (1.0L, LDBL_MANT_DIG) - 1.0L;
 183     volatile long double result = fmal (x, y, z);
 184     if (result != expected)
 185       failed_tests |= 1;
 186   }
 187   if ((LDBL_MANT_DIG % 2) == 1)
 188     {
 189       /* These tests fail on glibc 2.7 hppa,sparc, OSF/1 5.1.  */
 190       {
 191         volatile long double x = 1.0L + ldexpl (1.0L, - (LDBL_MANT_DIG + 1) / 2); /* 2^0 + 2^-27 */
 192         volatile long double y = 1.0L - ldexpl (1.0L, - (LDBL_MANT_DIG + 1) / 2); /* 2^0 - 2^-27 */
 193         volatile long double z = - ldexpl (1.0L, LDBL_MIN_EXP - LDBL_MANT_DIG); /* - 2^-1074 */
 194         /* x * y + z with infinite precision: 2^0 - 2^-54 - 2^-1074.
 195            Lies between (2^53 - 1) * 2^-53 and 2^53 * 2^-53 and is closer to
 196            (2^53 - 1) * 2^-53, therefore the rounding must round down and
 197            produce (2^53 - 1) * 2^-53.  */
 198         volatile long double expected = 1.0L - ldexpl (1.0L, - LDBL_MANT_DIG);
 199         volatile long double result = fmal (x, y, z);
 200         if (result != expected)
 201           failed_tests |= 4;
 202       }
 203       {
 204         volatile long double x = 1.0L + ldexpl (1.0L, - (LDBL_MANT_DIG + 1) / 2); /* 2^0 + 2^-57 */
 205         volatile long double y = x;
 206         volatile long double z = ldexpl (1.0L, - LDBL_MANT_DIG); /* 2^-113 */
 207         /* x * y + z with infinite precision: 2^0 + 2^-56 + 2^-113 + 2^-114.
 208            Lies between (2^112 + 2^56) * 2^-112 and (2^112 + 2^56 + 1) * 2^-112
 209            and is closer to (2^112 + 2^56 + 1) * 2^-112, therefore the rounding
 210            must round up and produce (2^112 + 2^56 + 1) * 2^-112.  */
 211         volatile long double expected =
 212           1.0L + ldexpl (1.0L, - (LDBL_MANT_DIG - 1) / 2) + ldexpl (1.0L, 1 - LDBL_MANT_DIG);
 213         volatile long double result = fmal (x, y, z);
 214         if (result != expected)
 215           failed_tests |= 4;
 216       }
 217     }
 218   else
 219     {
 220       /* These tests fail on glibc 2.11 x86,x86_64,powerpc, mingw.  */
 221       {
 222         volatile long double x = 1.0L + ldexpl (1.0L, - LDBL_MANT_DIG / 2); /* 2^0 + 2^-32 */
 223         volatile long double y = x;
 224         volatile long double z = ldexpl (1.0L, LDBL_MIN_EXP - LDBL_MANT_DIG); /* 2^-16445 */
 225         /* x * y + z with infinite precision: 2^0 + 2^-31 + 2^-64 + 2^-16445.
 226            Lies between (2^63 + 2^32 + 0) * 2^-63 and (2^63 + 2^32 + 1) * 2^-63
 227            and is closer to (2^63 + 2^32 + 1) * 2^-63, therefore the rounding
 228            must round up and produce (2^63 + 2^32 + 1) * 2^-63.  */
 229         volatile long double expected =
 230           1.0L + ldexpl (1.0L, 1 - LDBL_MANT_DIG / 2) + ldexpl (1.0L, 1 - LDBL_MANT_DIG);
 231         volatile long double result = fmal (x, y, z);
 232         if (result != expected)
 233           failed_tests |= 8;
 234       }
 235       {
 236         volatile long double x = 1.0L + ldexpl (1.0L, - LDBL_MANT_DIG / 2); /* 2^0 + 2^-32 */
 237         volatile long double y = x;
 238         volatile long double z = ldexpl (1.0L, - LDBL_MANT_DIG); /* 2^-64 */
 239         /* x * y + z with infinite precision: 2^0 + 2^-31 + 2^-63.
 240            Rounding must return this value unchanged.  */
 241         volatile long double expected = 1.0L + ldexpl (1.0L, 1 - LDBL_MANT_DIG / 2) + ldexpl (1.0L, 1 - LDBL_MANT_DIG);
 242         volatile long double result = fmal (x, y, z);
 243         if (result != expected)
 244           failed_tests |= 8;
 245       }
 246       {
 247         volatile long double x = 1.0L + ldexpl (1.0L, - LDBL_MANT_DIG / 2); /* 2^0 + 2^-32 */
 248         volatile long double y = x;
 249         volatile long double z = ldexpl (1.0L, 1 - LDBL_MANT_DIG); /* 2^-63 */
 250         /* x * y + z with infinite precision: 2^0 + 2^-31 + 2^-63 + 2^-64.
 251            Lies between (2^63 + 2^32 + 1) * 2^-63 and (2^63 + 2^32 + 2) * 2^-63
 252            and is at the same distance from each.  According to the round-to-even
 253            rule, the rounding must round up and produce (2^63 + 2^32 + 2) * 2^-63.  */
 254         volatile long double expected = 1.0L + ldexpl (1.0L, -31) + ldexpl (1.0L, -62);
 255         volatile long double result = fmal (x, y, z);
 256         if (result != expected)
 257           failed_tests |= 8;
 258       }
 259       {
 260         volatile long double x = 1.0L + ldexpl (1.0L, - LDBL_MANT_DIG / 2); /* 2^0 + 2^-32 */
 261         volatile long double y = x;
 262         volatile long double z = ldexpl (1.0L, LDBL_MANT_DIG / 2 + 1); /* 2^33 */
 263         /* x * y + z with infinite precision: 2^33 + 2^0 + 2^-31 + 2^-64.
 264            Lies between (2^63 + 2^30) * 2^-30 and (2^63 + 2^30 + 1) * 2^-30
 265            and is closer to (2^63 + 2^30 + 1) * 2^-30, therefore the rounding
 266            must round up and produce (2^63 + 2^30 + 1) * 2^-30.  */
 267         volatile long double expected = z + 1.0L + ldexp (1.0L, 2 - LDBL_MANT_DIG / 2);
 268         volatile long double result = fmal (x, y, z);
 269         if (result != expected)
 270           failed_tests |= 8;
 271       }
 272       {
 273         volatile long double x = 1.0L + ldexpl (1.0L, - LDBL_MANT_DIG / 2); /* 2^0 + 2^-32 */
 274         volatile long double y = x;
 275         volatile long double z = - ldexpl (1.0, 1 - LDBL_MANT_DIG); /* - 2^-63 */
 276         /* x * y + z with infinite precision: 2^0 + 2^-31 - 2^-64.
 277            Lies between (2^63 + 2^32 - 1) * 2^-63 and (2^63 + 2^32) * 2^-63
 278            and is at the same distance from each.  According to the round-to-even
 279            rule, the rounding must round up and produce (2^63 + 2^32) * 2^-63.  */
 280         volatile long double expected = 1.0L + ldexpl (1.0L, 1 - LDBL_MANT_DIG / 2);
 281         volatile long double result = fmal (x, y, z);
 282         if (result != expected)
 283           failed_tests |= 8;
 284       }
 285       {
 286         volatile long double x = 1.0L + ldexpl (1.0L, - LDBL_MANT_DIG / 2); /* 2^0 + 2^-32 */
 287         volatile long double y = x;
 288         volatile long double z = - 1.0L; /* - 2^0 */
 289         /* x * y + z with infinite precision: 2^-31 + 2^-64.
 290            Rounding must return this value unchanged.  */
 291         volatile long double expected = ldexpl (1.0L, 1 - LDBL_MANT_DIG / 2) + ldexpl (1.0L, - LDBL_MANT_DIG);
 292         volatile long double result = fmal (x, y, z);
 293         if (result != expected)
 294           failed_tests |= 8;
 295       }
 296       {
 297         volatile long double x = 1.0L + ldexpl (1.0L, - LDBL_MANT_DIG / 2); /* 2^0 + 2^-32 */
 298         volatile long double y = - x;
 299         volatile long double z = 2.0L; /* 2^1 */
 300         /* x * y + z with infinite precision: 2^0 - 2^31 - 2^-64.
 301            Rounding must return this value unchanged.  */
 302         volatile long double expected = 1.0L - ldexpl (1.0L, 1 - LDBL_MANT_DIG / 2) - ldexpl (1.0L, - LDBL_MANT_DIG);
 303         volatile long double result = fmal (x, y, z);
 304         if (result != expected)
 305           failed_tests |= 8;
 306       }
 307       {
 308         volatile long double x = 1.0L + ldexpl (1.0L, - LDBL_MANT_DIG / 2); /* 2^0 + 2^-32 */
 309         volatile long double y = - x;
 310         volatile long double z = ldexpl (1.0L, LDBL_MANT_DIG / 2 + 2); /* 2^34 */
 311         /* x * y + z with infinite precision: 2^34 - (2^0 + 2^-31 + 2^-64).
 312            Lies between (2^64 - 2^30 - 1) * 2^-30 and (2^64 - 2^30) * 2^-30
 313            and is closer to (2^64 - 2^30 - 1) * 2^-30, therefore the rounding
 314            must round down and produce (2^64 - 2^30 - 1) * 2^-30.  */
 315         volatile long double expected = z - 1.0L - ldexpl (1.0L, 2 - LDBL_MANT_DIG / 2);
 316         volatile long double result = fmal (x, y, z);
 317         if (result != expected)
 318           failed_tests |= 8;
 319       }
 320       {
 321         volatile long double x = 1.0L + ldexpl (1.0L, - LDBL_MANT_DIG / 2 - 1); /* 2^0 + 2^-33 */
 322         volatile long double y = 1.0L - ldexpl (1.0L, - LDBL_MANT_DIG / 2 - 1); /* 2^0 - 2^-33 */
 323         volatile long double z = - ldexpl (1.0L, - LDBL_MANT_DIG - 1); /* 2^-65 */
 324         /* x * y + z with infinite precision: 2^0 - 2^-65 - 2^-66.
 325            Lies between (2^64 - 1) * 2^-64 and 2^64 * 2^-64 and is closer to
 326            (2^64 - 1) * 2^-64, therefore the rounding must round down and
 327            produce (2^64 - 1) * 2^-64.  */
 328         volatile long double expected = 1.0L - ldexpl (1.0L, - LDBL_MANT_DIG);
 329         volatile long double result = fmal (x, y, z);
 330         if (result != expected)
 331           failed_tests |= 8;
 332       }
 333       {
 334         volatile long double x = 1.0L + ldexpl (1.0L, - LDBL_MANT_DIG / 2 - 1); /* 2^0 + 2^-33 */
 335         volatile long double y = 1.0L - ldexpl (1.0L, - LDBL_MANT_DIG / 2 - 1); /* 2^0 - 2^-33 */
 336         volatile long double z = - 1.0L; /* 2^0 */
 337         /* x * y + z with infinite precision: - 2^-66.
 338            Rounding must return this value unchanged.  */
 339         volatile long double expected = - ldexpl (1.0L, - LDBL_MANT_DIG - 2);
 340         volatile long double result = fmal (x, y, z);
 341         if (result != expected)
 342           failed_tests |= 8;
 343       }
 344     }
 345   /* This test fails on glibc 2.11 x86,x86_64,powerpc, glibc 2.7 hppa,sparc,
 346      FreeBSD 6.4 x86, mingw.  */
 347   {
 348     long double minus_inf = -1.0L / p0;
 349     volatile long double x = ldexpl (1.0L, LDBL_MAX_EXP - 1);
 350     volatile long double y = ldexpl (1.0L, LDBL_MAX_EXP - 1);
 351     volatile long double z = minus_inf;
 352     volatile long double result = fmal (x, y, z);
 353     if (!(result == minus_inf))
 354       failed_tests |= 16;
 355   }
 356   /* This test fails on glibc 2.11 x86,x86_64,powerpc glibc 2.7 hppa,sparc,
 357      Mac OS X 10.5, FreeBSD 6.4 x86, OSF/1 5.1, mingw.  */
 358   {
 359     volatile long double x = ldexpl (1.0L, LDBL_MAX_EXP - 1);
 360     volatile long double y = 2.0L;
 361     volatile long double z =
 362       - ldexpl (ldexpl (1.0L, LDBL_MAX_EXP - 1) - ldexpl (1.0L, LDBL_MAX_EXP - LDBL_MANT_DIG - 1), 1);
 363     volatile long double expected = ldexpl (1.0L, LDBL_MAX_EXP - LDBL_MANT_DIG);
 364     volatile long double result = fmal (x, y, z);
 365     if (result != expected)
 366       failed_tests |= 32;
 367   }
 368   return failed_tests;
 369 }]])],
 370         [gl_cv_func_fmal_works=yes],
 371         [gl_cv_func_fmal_works=no],
 372         [dnl Guess yes on native Windows with MSVC.
 373          dnl Otherwise guess no, even on glibc systems.
 374          gl_cv_func_fmal_works="$gl_cross_guess_normal"
 375          case "$host_os" in
 376            mingw*)
 377              AC_EGREP_CPP([Known], [
 378 #ifdef _MSC_VER
 379  Known
 380 #endif
 381                ], [gl_cv_func_fmal_works="guessing yes"])
 382              ;;
 383          esac
 384         ])
 385     ])
 386   LIBS="$save_LIBS"
 387 ])
 388
 389 # Prerequisites of lib/fmal.c.
 390 AC_DEFUN([gl_PREREQ_FMAL], [:])