m4/fmal.m4

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