math: sync with GNU libc

The format of the ULPS files have changed, non-glibc architecture
files needs to be updated later.
Add all math tests from latest GNU libc and allow to compile
and run them on uClibc-ng and GNU libc systems.

1 files changed, 11329 insertions, 3624 deletions

diff --git a/test/math/libm-test.inc b/test/math/libm-test.inc
index f50b48b..2fc97ee 100644
--- a/test/math/libm-test.inc
+++ b/test/math/libm-test.inc
@@ -1,4 +1,4 @@
-/* Copyright (C) 1997-2002, 2003, 2004, 2005, 2006 Free Software Foundation, Inc.
+/* Copyright (C) 1997-2016 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Andreas Jaeger <aj@suse.de>, 1997.
 
@@ -24,40 +24,46 @@
    Macros:
    FUNC(function): converts general function name (like cos) to
    name with correct suffix (e.g. cosl or cosf)
-   MATHCONST(x):   like FUNC but for constants (e.g convert 0.0 to 0.0L)
-   FLOAT:          floating point type to test
-   - TEST_MSG:     informal message to be displayed
-   CHOOSE(Clongdouble,Cdouble,Cfloat,Cinlinelongdouble,Cinlinedouble,Cinlinefloat):
+   FLOAT:	   floating point type to test
+   - TEST_MSG:	   informal message to be displayed
    chooses one of the parameters as delta for testing
    equality
-   PRINTF_EXPR     Floating point conversion specification to print a variable
+   PRINTF_EXPR	   Floating point conversion specification to print a variable
    of type FLOAT with printf.  PRINTF_EXPR just contains
    the specifier, not the percent and width arguments,
    e.g. "f".
    PRINTF_XEXPR	   Like PRINTF_EXPR, but print in hexadecimal format.
-   PRINTF_NEXPR    Like PRINTF_EXPR, but print nice.  */
+   PRINTF_NEXPR Like PRINTF_EXPR, but print nice.
+   PREFIX A macro which defines the prefix for common macros for the
+   type (i.e LDBL, DBL, or FLT).
+   LIT A function which appends the correct suffix to a literal.
+   TYPE_STR A macro which defines a stringitized name of the type.
+   FTOSTR This macro defines a function similar in type to snprintf
+   which converts a FLOAT to a string.  */
 
 /* This testsuite has currently tests for:
    acos, acosh, asin, asinh, atan, atan2, atanh,
-   cbrt, ceil, copysign, cos, cosh, erf, erfc, exp, exp10, exp2, expm1,
-   fabs, fdim, floor, fma, fmax, fmin, fmod, fpclassify,
+   cbrt, ceil, copysign, cos, cosh, drem, erf, erfc, exp, exp10, exp2, expm1,
+   fabs, fdim, finite, floor, fma, fmax, fmin, fmod, fpclassify,
    frexp, gamma, hypot,
-   ilogb, isfinite, isinf, isnan, isnormal,
+   ilogb, isfinite, isinf, isnan, isnormal, issignaling,
    isless, islessequal, isgreater, isgreaterequal, islessgreater, isunordered,
    j0, j1, jn,
    ldexp, lgamma, log, log10, log1p, log2, logb,
-   modf, nearbyint, nextafter,
-   pow, remainder, remquo, rint, lrint, llrint,
+   modf, nearbyint, nextafter, nexttoward,
+   pow, pow10, remainder, remquo, rint, lrint, llrint,
    round, lround, llround,
    scalb, scalbn, scalbln, signbit, sin, sincos, sinh, sqrt, tan, tanh, tgamma, trunc,
    y0, y1, yn, significand
 
    and for the following complex math functions:
    cabs, cacos, cacosh, carg, casin, casinh, catan, catanh,
-   ccos, ccosh, cexp, clog, cpow, cproj, csin, csinh, csqrt, ctan, ctanh.
+   ccos, ccosh, cexp, cimag, clog, clog10, conj, cpow, cproj, creal,
+   csin, csinh, csqrt, ctan, ctanh.
 
-   At the moment the following functions aren't tested:
-   drem, nan
+   At the moment the following functions and macros aren't tested:
+   lgamma_r,
+   nan.
 
    Parameter handling is primitive in the moment:
    --verbose=[0..3] for different levels of output:
@@ -78,13 +84,7 @@
    against.  These implemented tests should check all cases that are
    specified in ISO C99.
 
-   Exception testing: At the moment only divide-by-zero and invalid
-   exceptions are tested.  Overflow/underflow and inexact exceptions
-   aren't checked at the moment.
-
-   NaN values: There exist signalling and quiet NaNs.  This implementation
-   only uses quiet NaN as parameter but does not differenciate
-   between the two kinds of NaNs as result.
+   NaN values: The payload of NaNs is not examined.
 
    Inline functions: Inlining functions should give an improvement in
    speed - but not in precission.  The inlined functions return
@@ -115,14 +115,14 @@
 # define _GNU_SOURCE
 #endif
 
-#undef __CHK_COMPLEX_STUFF
-#define __CHK_COMPLEX_STUFF 0
-
-#include "libm-test-ulps.h"
 #include <complex.h>
 #include <math.h>
 #include <float.h>
+#if defined(__GLIBC__) || defined(__UCLIBC__) && defined(__UCLIBC_HAS_FENV__)
+#include <fenv.h>
+#else
 #include "fenv.h"
+#endif
 #include <limits.h>
 
 #include <errno.h>
@@ -131,93 +131,311 @@
 #include <string.h>
 #include <getopt.h>
 
+#define TININESS_AFTER_ROUNDING 0
+
+#include "math-tests.h"
+#include "math-tests-arch.h"
+
+/* This header defines func_ulps, func_real_ulps and func_imag_ulps
+   arrays.  */
+#include "libm-test-ulps.h"
+
+/* Allow platforms without all rounding modes to test properly,
+   assuming they provide an __FE_UNDEFINED in <bits/fenv.h> which
+   causes fesetround() to return failure.  */
+#ifndef FE_TONEAREST
+# define FE_TONEAREST	__FE_UNDEFINED
+#endif
+#ifndef FE_TOWARDZERO
+# define FE_TOWARDZERO	__FE_UNDEFINED
+#endif
+#ifndef FE_UPWARD
+# define FE_UPWARD	__FE_UNDEFINED
+#endif
+#ifndef FE_DOWNWARD
+# define FE_DOWNWARD	__FE_UNDEFINED
+#endif
+
 /* Possible exceptions */
 #define NO_EXCEPTION			0x0
 #define INVALID_EXCEPTION		0x1
 #define DIVIDE_BY_ZERO_EXCEPTION	0x2
+#define OVERFLOW_EXCEPTION		0x4
+#define UNDERFLOW_EXCEPTION		0x8
+#define INEXACT_EXCEPTION		0x10
 /* The next flags signals that those exceptions are allowed but not required.   */
-#define INVALID_EXCEPTION_OK		0x4
-#define DIVIDE_BY_ZERO_EXCEPTION_OK	0x8
+#define INVALID_EXCEPTION_OK		0x20
+#define DIVIDE_BY_ZERO_EXCEPTION_OK	0x40
+#define OVERFLOW_EXCEPTION_OK		0x80
+#define UNDERFLOW_EXCEPTION_OK		0x100
+/* For "inexact" exceptions, the default is allowed but not required
+   unless INEXACT_EXCEPTION or NO_INEXACT_EXCEPTION is specified.  */
+#define NO_INEXACT_EXCEPTION		0x200
 #define EXCEPTIONS_OK INVALID_EXCEPTION_OK+DIVIDE_BY_ZERO_EXCEPTION_OK
-/* Some special test flags, passed togther with exceptions.  */
-#define IGNORE_ZERO_INF_SIGN		0x10
-
-/* Various constants (we must supply them precalculated for accuracy).  */
-#define M_PI_6l			.52359877559829887307710723054658383L
-#define M_E2l			7.389056098930650227230427460575008L
-#define M_E3l			20.085536923187667740928529654581719L
-#define M_2_SQRT_PIl		3.5449077018110320545963349666822903L	/* 2 sqrt (M_PIl)  */
-#define M_SQRT_PIl		1.7724538509055160272981674833411451L	/* sqrt (M_PIl)  */
-#define M_LOG_SQRT_PIl		0.57236494292470008707171367567652933L	/* log(sqrt(M_PIl))  */
-#define M_LOG_2_SQRT_PIl	1.265512123484645396488945797134706L	/* log(2*sqrt(M_PIl))  */
-#define M_PI_34l		(M_PIl - M_PI_4l)		/* 3*pi/4 */
-#define M_PI_34_LOG10El		(M_PIl - M_PI_4l) * M_LOG10El
-#define M_PI2_LOG10El		M_PI_2l * M_LOG10El
-#define M_PI4_LOG10El		M_PI_4l * M_LOG10El
-#define M_PI_LOG10El		M_PIl * M_LOG10El
-#define M_SQRT_2_2		0.70710678118654752440084436210484903L /* sqrt (2) / 2 */
-
-static FILE *ulps_file;	/* File to document difference.  */
-static int output_ulps;	/* Should ulps printed?  */
+/* Some special test flags, passed together with exceptions.  */
+#define IGNORE_ZERO_INF_SIGN		0x400
+#define TEST_NAN_SIGN			0x800
+#define NO_TEST_INLINE			0x1000
+#define XFAIL_TEST			0x2000
+/* Indicate errno settings required or disallowed.  */
+#define ERRNO_UNCHANGED			0x4000
+#define ERRNO_EDOM			0x8000
+#define ERRNO_ERANGE			0x10000
+/* Flags generated by gen-libm-test.pl, not entered here manually.  */
+#define IGNORE_RESULT			0x20000
+#define NON_FINITE			0x40000
+#define TEST_SNAN			0x80000
+
+#define __CONCATX(a,b) __CONCAT(a,b)
+
+#define TYPE_MIN __CONCATX (PREFIX, _MIN)
+#define TYPE_TRUE_MIN __CONCATX (PREFIX, _TRUE_MIN)
+#define TYPE_MAX __CONCATX (PREFIX, _MAX)
+#define MIN_EXP __CONCATX (PREFIX, _MIN_EXP)
+#define MAX_EXP __CONCATX (PREFIX, _MAX_EXP)
+#define MANT_DIG __CONCATX (PREFIX, _MANT_DIG)
+
+/* Maximum character buffer to store a stringitized FLOAT value.  */
+#define FSTR_MAX (128)
+
+#if TEST_INLINE
+# define ULP_IDX __CONCATX (ULP_I_, PREFIX)
+# define QTYPE_STR "i" TYPE_STR
+#else
+# define ULP_IDX __CONCATX (ULP_, PREFIX)
+# define QTYPE_STR TYPE_STR
+#endif
+
+/* Format specific test macros.  */
+#define TEST_COND_binary32 (MANT_DIG == 24	\
+			    && MIN_EXP == -125	\
+			    && MAX_EXP == 128)
+
+#define TEST_COND_binary64 (MANT_DIG == 53	\
+			    && MIN_EXP == -1021	\
+			    && MAX_EXP == 1024)
+
+#define TEST_COND_binary128 (MANT_DIG == 113		\
+			     && MIN_EXP == -16381	\
+			     && MAX_EXP == 16384)
+
+#define TEST_COND_ibm128 (MANT_DIG == 106)
+
+#define TEST_COND_intel96 (MANT_DIG == 64	\
+			   && MIN_EXP == -16381	\
+			   && MAX_EXP == 16384)
+
+#define TEST_COND_m68k96 (MANT_DIG == 64	\
+			  && MIN_EXP == -16382	\
+			  && MAX_EXP == 16384)
+
+/* Values underflowing only for float.  */
+#if TEST_COND_binary32
+# define UNDERFLOW_EXCEPTION_FLOAT	UNDERFLOW_EXCEPTION
+# define UNDERFLOW_EXCEPTION_OK_FLOAT	UNDERFLOW_EXCEPTION_OK
+#else
+# define UNDERFLOW_EXCEPTION_FLOAT	0
+# define UNDERFLOW_EXCEPTION_OK_FLOAT	0
+#endif
+
+/* Values underflowing only for double or types with a larger least
+   positive normal value.  */
+#if TEST_COND_binary32 || TEST_COND_binary64 || TEST_COND_ibm128
+# define UNDERFLOW_EXCEPTION_DOUBLE	UNDERFLOW_EXCEPTION
+# define UNDERFLOW_EXCEPTION_OK_DOUBLE	UNDERFLOW_EXCEPTION_OK
+#else
+# define UNDERFLOW_EXCEPTION_DOUBLE	0
+# define UNDERFLOW_EXCEPTION_OK_DOUBLE	0
+#endif
+
+/* Values underflowing only for IBM long double or types with a larger least
+   positive normal value.  */
+#if TEST_COND_binary32 || TEST_COND_ibm128
+# define UNDERFLOW_EXCEPTION_LDOUBLE_IBM	UNDERFLOW_EXCEPTION
+#else
+# define UNDERFLOW_EXCEPTION_LDOUBLE_IBM	0
+#endif
+
+/* Values underflowing on architectures detecting tininess before
+   rounding, but not on those detecting tininess after rounding.  */
+#define UNDERFLOW_EXCEPTION_BEFORE_ROUNDING	(TININESS_AFTER_ROUNDING \
+						 ? 0			\
+						 : UNDERFLOW_EXCEPTION)
+
+#if LONG_MAX == 0x7fffffff
+# define TEST_COND_long32	1
+# define TEST_COND_long64	0
+#else
+# define TEST_COND_long32	0
+# define TEST_COND_long64	1
+#endif
+#define TEST_COND_before_rounding	(!TININESS_AFTER_ROUNDING)
+#define TEST_COND_after_rounding	TININESS_AFTER_ROUNDING
+
+/* Various constants derived from pi.  We must supply them precalculated for
+   accuracy.  They are written as a series of postfix operations to keep
+   them concise yet somewhat readable.  */
+
+/* (pi * 3) / 4 */
+#define lit_pi_3_m_4_d		LIT (2.356194490192344928846982537459627163)
+/* pi * 3 / (4 * ln(10)) */
+#define lit_pi_3_m_4_ln10_m_d	LIT (1.023282265381381010614337719073516828)
+/* pi / (2 * ln(10)) */
+#define lit_pi_2_ln10_m_d	LIT (0.682188176920920673742891812715677885)
+/* pi / (4 * ln(10)) */
+#define lit_pi_4_ln10_m_d	LIT (0.341094088460460336871445906357838943)
+/* pi / ln(10) */
+#define lit_pi_ln10_d		LIT (1.364376353841841347485783625431355770)
+/* pi / 2 */
+#define lit_pi_2_d		LITM (M_PI_2)
+/* pi / 4 */
+#define lit_pi_4_d		LITM (M_PI_4)
+/* pi */
+#define lit_pi			LITM (M_PI)
+
+/* Other useful constants.  */
+
+/* e */
+#define lit_e			LITM (M_E)
+
+#define ulps_file_name "ULPs"	/* Name of the ULPs file.  */
+static FILE *ulps_file;		/* File to document difference.  */
+static int output_ulps;		/* Should ulps printed?  */
+static char *output_dir;	/* Directory where generated files will be written.  */
 
 static int noErrors;	/* number of errors */
 static int noTests;	/* number of tests (without testing exceptions) */
 static int noExcTests;	/* number of tests for exception flags */
-static int noXFails;	/* number of expected failures.  */
-static int noXPasses;	/* number of unexpected passes.  */
+static int noErrnoTests;/* number of tests for errno values */
 
 static int verbose;
 static int output_max_error;	/* Should the maximal errors printed?  */
 static int output_points;	/* Should the single function results printed?  */
 static int ignore_max_ulp;	/* Should we ignore max_ulp?  */
 
-static FLOAT minus_zero, plus_zero;
-static FLOAT plus_infty, minus_infty, nan_value, max_value, min_value;
+#define plus_zero	LIT (0.0)
+#define minus_zero	LIT (-0.0)
+#define plus_infty	FUNC (__builtin_inf) ()
+#define minus_infty	-(FUNC (__builtin_inf) ())
+#define qnan_value	FUNC (__builtin_nan) ("")
+#define snan_value	FUNC (__builtin_nans) ("")
+#define max_value	TYPE_MAX
+#define min_value	TYPE_MIN
+#define min_subnorm_value TYPE_TRUE_MIN
+
+/* For nexttoward tests.  */
+#define snan_value_ld	__builtin_nansl ("")
 
 static FLOAT max_error, real_max_error, imag_max_error;
 
+static FLOAT prev_max_error, prev_real_max_error, prev_imag_max_error;
 
-#define BUILD_COMPLEX(real, imag) \
-  ({ __complex__ FLOAT __retval;					      \
-     __real__ __retval = (real);					      \
-     __imag__ __retval = (imag);					      \
-     __retval; })
+static FLOAT max_valid_error;
 
-#define BUILD_COMPLEX_INT(real, imag) \
-  ({ __complex__ int __retval;						      \
-     __real__ __retval = (real);					      \
-     __imag__ __retval = (imag);					      \
-     __retval; })
+/* Sufficient numbers of digits to represent any floating-point value
+   unambiguously (for any choice of the number of bits in the first
+   hex digit, in the case of TYPE_HEX_DIG).  When used with printf
+   formats where the precision counts only digits after the point, 1
+   is subtracted from these values. */
+#define TYPE_DECIMAL_DIG __CONCATX (PREFIX, _DECIMAL_DIG)
+#define TYPE_HEX_DIG ((MANT_DIG + 6) / 4)
+
+/* Compare KEY (a string, with the name of a function) with ULP (a
+   pointer to a struct ulp_data structure), returning a value less
+   than, equal to or greater than zero for use in bsearch.  */
+
+static int
+compare_ulp_data (const void *key, const void *ulp)
+{
+  const char *keystr = key;
+  const struct ulp_data *ulpdat = ulp;
+  return strcmp (keystr, ulpdat->name);
+}
 
+/* Return the ulps for NAME in array DATA with NMEMB elements, or 0 if
+   no ulps listed.  */
 
-#define MANT_DIG CHOOSE ((LDBL_MANT_DIG-1), (DBL_MANT_DIG-1), (FLT_MANT_DIG-1),  \
-                         (LDBL_MANT_DIG-1), (DBL_MANT_DIG-1), (FLT_MANT_DIG-1))
+static FLOAT
+find_ulps (const char *name, const struct ulp_data *data, size_t nmemb)
+{
+  const struct ulp_data *entry = bsearch (name, data, nmemb, sizeof (*data),
+					  compare_ulp_data);
+  if (entry == NULL)
+    return 0;
+  else
+    return entry->max_ulp[ULP_IDX];
+}
 
 static void
-init_max_error (void)
+init_max_error (const char *name, int exact)
 {
   max_error = 0;
   real_max_error = 0;
   imag_max_error = 0;
+  prev_max_error = find_ulps (name, func_ulps,
+			      sizeof (func_ulps) / sizeof (func_ulps[0]));
+  prev_real_max_error = find_ulps (name, func_real_ulps,
+				   (sizeof (func_real_ulps)
+				    / sizeof (func_real_ulps[0])));
+  prev_imag_max_error = find_ulps (name, func_imag_ulps,
+				   (sizeof (func_imag_ulps)
+				    / sizeof (func_imag_ulps[0])));
+#if TEST_COND_ibm128
+  /* The documented accuracy of IBM long double division is 3ulp (see
+     libgcc/config/rs6000/ibm-ldouble-format), so do not require
+     better accuracy for libm functions that are exactly defined for
+     other formats.  */
+  max_valid_error = exact ? 3 : 14;
+#else
+  max_valid_error = exact ? 0 : 9;
+#endif
+  prev_max_error = (prev_max_error <= max_valid_error
+		    ? prev_max_error
+		    : max_valid_error);
+  prev_real_max_error = (prev_real_max_error <= max_valid_error
+			 ? prev_real_max_error
+			 : max_valid_error);
+  prev_imag_max_error = (prev_imag_max_error <= max_valid_error
+			 ? prev_imag_max_error
+			 : max_valid_error);
   feclearexcept (FE_ALL_EXCEPT);
+  errno = 0;
 }
 
 static void
 set_max_error (FLOAT current, FLOAT *curr_max_error)
 {
-  if (current > *curr_max_error)
+  if (current > *curr_max_error && current <= max_valid_error)
     *curr_max_error = current;
 }
 
 
+/* Print a FLOAT.  */
+static void
+print_float (FLOAT f)
+{
+  /* As printf doesn't differ between a sNaN and a qNaN, do this manually.  */
+  if (issignaling (f))
+    printf ("sNaN\n");
+  else if (isnan (f))
+    printf ("qNaN\n");
+  else
+    {
+      char fstrn[FSTR_MAX], fstrx[FSTR_MAX];
+      FTOSTR (fstrn, FSTR_MAX, "% .*" PRINTF_EXPR, TYPE_DECIMAL_DIG - 1, f);
+      FTOSTR (fstrx, FSTR_MAX, "% .*" PRINTF_XEXPR, TYPE_HEX_DIG - 1, f);
+      printf ("%s  %s\n", fstrn, fstrx);
+    }
+}
+
 /* Should the message print to screen?  This depends on the verbose flag,
    and the test status.  */
 static int
-print_screen (int ok, int xfail)
+print_screen (int ok)
 {
   if (output_points
       && (verbose > 1
-	  || (verbose == 1 && ok == xfail)))
+	  || (verbose == 1 && ok == 0)))
     return 1;
   return 0;
 }
@@ -226,51 +444,33 @@ print_screen (int ok, int xfail)
 /* Should the message print to screen?  This depends on the verbose flag,
    and the test status.  */
 static int
-print_screen_max_error (int ok, int xfail)
+print_screen_max_error (int ok)
 {
   if (output_max_error
       && (verbose > 1
-	  || ((verbose == 1) && (ok == xfail))))
+	  || ((verbose == 1) && (ok == 0))))
     return 1;
   return 0;
 }
 
 /* Update statistic counters.  */
 static void
-update_stats (int ok, int xfail)
+update_stats (int ok)
 {
   ++noTests;
-  if (ok && xfail)
-    ++noXPasses;
-  else if (!ok && xfail)
-    ++noXFails;
-  else if (!ok && !xfail)
+  if (!ok)
     ++noErrors;
 }
 
 static void
-print_ulps (const char *test_name, FLOAT ulp)
-{
-  if (output_ulps)
-    {
-      fprintf (ulps_file, "Test \"%s\":\n", test_name);
-      fprintf (ulps_file, "%s: %.0" PRINTF_NEXPR "\n",
-	       CHOOSE("ldouble", "double", "float",
-		      "ildouble", "idouble", "ifloat"),
-	       FUNC(ceil) (ulp));
-    }
-}
-
-static void
 print_function_ulps (const char *function_name, FLOAT ulp)
 {
   if (output_ulps)
     {
+      char ustrn[FSTR_MAX];
+      FTOSTR (ustrn, FSTR_MAX, "%.0" PRINTF_NEXPR, FUNC (ceil) (ulp));
       fprintf (ulps_file, "Function: \"%s\":\n", function_name);
-      fprintf (ulps_file, "%s: %.0" PRINTF_NEXPR "\n",
-	       CHOOSE("ldouble", "double", "float",
-		      "ildouble", "idouble", "ifloat"),
-	       FUNC(ceil) (ulp));
+      fprintf (ulps_file, QTYPE_STR ": %s\n", ustrn);
     }
 }
 
@@ -281,21 +481,20 @@ print_complex_function_ulps (const char *function_name, FLOAT real_ulp,
 {
   if (output_ulps)
     {
+      char fstrn[FSTR_MAX];
       if (real_ulp != 0.0)
 	{
+	  FTOSTR (fstrn, FSTR_MAX, "%.0" PRINTF_NEXPR,
+	            FUNC (ceil) (real_ulp));
 	  fprintf (ulps_file, "Function: Real part of \"%s\":\n", function_name);
-	  fprintf (ulps_file, "%s: %.0" PRINTF_NEXPR "\n",
-		   CHOOSE("ldouble", "double", "float",
-			  "ildouble", "idouble", "ifloat"),
-		   FUNC(ceil) (real_ulp));
+	  fprintf (ulps_file, QTYPE_STR ": %s\n", fstrn);
 	}
       if (imag_ulp != 0.0)
 	{
+	  FTOSTR (fstrn, FSTR_MAX, "%.0" PRINTF_NEXPR,
+	            FUNC (ceil) (imag_ulp));
 	  fprintf (ulps_file, "Function: Imaginary part of \"%s\":\n", function_name);
-	  fprintf (ulps_file, "%s: %.0" PRINTF_NEXPR "\n",
-		   CHOOSE("ldouble", "double", "float",
-			  "ildouble", "idouble", "ifloat"),
-		   FUNC(ceil) (imag_ulp));
+	  fprintf (ulps_file, QTYPE_STR ": %s\n", fstrn);
 	}
 
 
@@ -308,11 +507,11 @@ print_complex_function_ulps (const char *function_name, FLOAT real_ulp,
 static void
 fpstack_test (const char *test_name)
 {
-#ifdef i386
+#if defined (__i386__) || defined (__x86_64__)
   static int old_stack;
   int sw;
 
-  __asm__ ("fnstsw" : "=a" (sw));
+  asm ("fnstsw" : "=a" (sw));
   sw >>= 11;
   sw &= 7;
 
@@ -328,11 +527,11 @@ fpstack_test (const char *test_name)
 
 
 static void
-print_max_error (const char *func_name, FLOAT allowed, int xfail)
+print_max_error (const char *func_name)
 {
   int ok = 0;
 
-  if (max_error == 0.0 || (max_error <= allowed && !ignore_max_ulp))
+  if (max_error == 0.0 || (max_error <= prev_max_error && !ignore_max_ulp))
     {
       ok = 1;
     }
@@ -341,53 +540,69 @@ print_max_error (const char *func_name, FLOAT allowed, int xfail)
     print_function_ulps (func_name, max_error);
 
 
-  if (print_screen_max_error (ok, xfail))
+  if (print_screen_max_error (ok))
     {
+      char mestr[FSTR_MAX], pmestr[FSTR_MAX];
+      FTOSTR (mestr, FSTR_MAX, "%.0" PRINTF_NEXPR, FUNC (ceil) (max_error));
+      FTOSTR (pmestr, FSTR_MAX, "%.0" PRINTF_NEXPR, FUNC (ceil) (prev_max_error));
       printf ("Maximal error of `%s'\n", func_name);
-      printf (" is      : %.0" PRINTF_NEXPR " ulp\n", FUNC(ceil) (max_error));
-      printf (" accepted: %.0" PRINTF_NEXPR " ulp\n", FUNC(ceil) (allowed));
+      printf (" is      : %s ulp\n", mestr);
+      printf (" accepted: %s ulp\n", pmestr);
     }
 
-  update_stats (ok, xfail);
+  update_stats (ok);
 }
 
 
 static void
-print_complex_max_error (const char *func_name, __complex__ FLOAT allowed,
-			 __complex__ int xfail)
+print_complex_max_error (const char *func_name)
 {
-  int ok = 0;
+  int real_ok = 0, imag_ok = 0, ok;
 
-  if ((real_max_error == 0 && imag_max_error == 0)
-      || (real_max_error <= __real__ allowed
-	  && imag_max_error <= __imag__ allowed
-	  && !ignore_max_ulp))
+  if (real_max_error == 0
+      || (real_max_error <= prev_real_max_error && !ignore_max_ulp))
     {
-      ok = 1;
+      real_ok = 1;
     }
 
-  if (!ok)
-    print_complex_function_ulps (func_name, real_max_error, imag_max_error);
+  if (imag_max_error == 0
+      || (imag_max_error <= prev_imag_max_error && !ignore_max_ulp))
+    {
+      imag_ok = 1;
+    }
 
+  ok = real_ok && imag_ok;
 
-  if (print_screen_max_error (ok, xfail))
+  if (!ok)
+    print_complex_function_ulps (func_name,
+				 real_ok ? 0 : real_max_error,
+				 imag_ok ? 0 : imag_max_error);
+
+  if (print_screen_max_error (ok))
     {
+      char rmestr[FSTR_MAX], prmestr[FSTR_MAX];
+      char imestr[FSTR_MAX], pimestr[FSTR_MAX];
+      FTOSTR (rmestr, FSTR_MAX, "%.0" PRINTF_NEXPR,
+		FUNC (ceil) (real_max_error));
+      FTOSTR (prmestr, FSTR_MAX, "%.0" PRINTF_NEXPR,
+		FUNC (ceil) (prev_real_max_error));
+      FTOSTR (imestr, FSTR_MAX, "%.0" PRINTF_NEXPR,
+		FUNC (ceil) (imag_max_error));
+      FTOSTR (pimestr, FSTR_MAX, "%.0" PRINTF_NEXPR,
+		FUNC (ceil) (prev_imag_max_error));
       printf ("Maximal error of real part of: %s\n", func_name);
-      printf (" is      : %.0" PRINTF_NEXPR " ulp\n",
-	      FUNC(ceil) (real_max_error));
-      printf (" accepted: %.0" PRINTF_NEXPR " ulp\n",
-	      FUNC(ceil) (__real__ allowed));
+      printf (" is      : %s ulp\n", rmestr);
+      printf (" accepted: %s ulp\n", prmestr);
       printf ("Maximal error of imaginary part of: %s\n", func_name);
-      printf (" is      : %.0" PRINTF_NEXPR " ulp\n",
-	      FUNC(ceil) (imag_max_error));
-      printf (" accepted: %.0" PRINTF_NEXPR " ulp\n",
-	      FUNC(ceil) (__imag__ allowed));
+      printf (" is      : %s ulp\n", imestr);
+      printf (" accepted: %s ulp\n", pimestr);
     }
 
-  update_stats (ok, xfail);
+  update_stats (ok);
 }
 
 
+#if FE_ALL_EXCEPT
 /* Test whether a given exception was raised.  */
 static void
 test_single_exception (const char *test_name,
@@ -396,19 +611,18 @@ test_single_exception (const char *test_name,
 		       int fe_flag,
 		       const char *flag_name)
 {
-#ifndef TEST_INLINE
   int ok = 1;
   if (exception & exc_flag)
     {
       if (fetestexcept (fe_flag))
 	{
-	  if (print_screen (1, 0))
+	  if (print_screen (1))
 	    printf ("Pass: %s: Exception \"%s\" set\n", test_name, flag_name);
 	}
       else
 	{
 	  ok = 0;
-	  if (print_screen (0, 0))
+	  if (print_screen (0))
 	    printf ("Failure: %s: Exception \"%s\" not set\n",
 		    test_name, flag_name);
 	}
@@ -418,23 +632,21 @@ test_single_exception (const char *test_name,
       if (fetestexcept (fe_flag))
 	{
 	  ok = 0;
-	  if (print_screen (0, 0))
+	  if (print_screen (0))
 	    printf ("Failure: %s: Exception \"%s\" set\n",
 		    test_name, flag_name);
 	}
       else
 	{
-	  if (print_screen (1, 0))
+	  if (print_screen (1))
 	    printf ("%s: Exception \"%s\" not set\n", test_name,
 		    flag_name);
 	}
     }
   if (!ok)
     ++noErrors;
-
-#endif
 }
-
+#endif
 
 /* Test whether exceptions given by EXCEPTION are raised.  Ignore thereby
    allowed but not required exceptions.
@@ -442,35 +654,162 @@ test_single_exception (const char *test_name,
 static void
 test_exceptions (const char *test_name, int exception)
 {
-  ++noExcTests;
+  if (TEST_EXCEPTIONS && EXCEPTION_TESTS (FLOAT))
+    {
+      ++noExcTests;
 #ifdef FE_DIVBYZERO
-  if ((exception & DIVIDE_BY_ZERO_EXCEPTION_OK) == 0)
-    test_single_exception (test_name, exception,
-			   DIVIDE_BY_ZERO_EXCEPTION, FE_DIVBYZERO,
-			   "Divide by zero");
+      if ((exception & DIVIDE_BY_ZERO_EXCEPTION_OK) == 0)
+	test_single_exception (test_name, exception,
+			       DIVIDE_BY_ZERO_EXCEPTION, FE_DIVBYZERO,
+			       "Divide by zero");
 #endif
 #ifdef FE_INVALID
-  if ((exception & INVALID_EXCEPTION_OK) == 0)
-    test_single_exception (test_name, exception, INVALID_EXCEPTION, FE_INVALID,
-			 "Invalid operation");
+      if ((exception & INVALID_EXCEPTION_OK) == 0)
+	test_single_exception (test_name, exception,
+			       INVALID_EXCEPTION, FE_INVALID,
+			       "Invalid operation");
+#endif
+#ifdef FE_OVERFLOW
+      if ((exception & OVERFLOW_EXCEPTION_OK) == 0)
+	test_single_exception (test_name, exception, OVERFLOW_EXCEPTION,
+			       FE_OVERFLOW, "Overflow");
 #endif
+      /* Spurious "underflow" and "inexact" exceptions are always
+	 allowed for IBM long double, in line with the underlying
+	 arithmetic.  */
+#ifdef FE_UNDERFLOW
+      if ((exception & UNDERFLOW_EXCEPTION_OK) == 0
+	  && !(TEST_COND_ibm128
+	       && (exception & UNDERFLOW_EXCEPTION) == 0))
+	test_single_exception (test_name, exception, UNDERFLOW_EXCEPTION,
+			       FE_UNDERFLOW, "Underflow");
+#endif
+#ifdef FE_INEXACT
+      if ((exception & (INEXACT_EXCEPTION | NO_INEXACT_EXCEPTION)) != 0
+	  && !(TEST_COND_ibm128
+	       && (exception & NO_INEXACT_EXCEPTION) != 0))
+	test_single_exception (test_name, exception, INEXACT_EXCEPTION,
+			       FE_INEXACT, "Inexact");
+#endif
+    }
   feclearexcept (FE_ALL_EXCEPT);
 }
 
+/* Test whether errno for TEST_NAME, set to ERRNO_VALUE, has value
+   EXPECTED_VALUE (description EXPECTED_NAME).  */
+static void
+test_single_errno (const char *test_name, int errno_value,
+		   int expected_value, const char *expected_name)
+{
+  if (errno_value == expected_value)
+    {
+      if (print_screen (1))
+	printf ("Pass: %s: errno set to %d (%s)\n", test_name, errno_value,
+		expected_name);
+    }
+  else
+    {
+      ++noErrors;
+      if (print_screen (0))
+	printf ("Failure: %s: errno set to %d, expected %d (%s)\n",
+		test_name, errno_value, expected_value, expected_name);
+    }
+}
+
+/* Test whether errno (value ERRNO_VALUE) has been for TEST_NAME set
+   as required by EXCEPTIONS.  */
+static void
+test_errno (const char *test_name, int errno_value, int exceptions)
+{
+  if (TEST_ERRNO)
+    {
+      ++noErrnoTests;
+      if (exceptions & ERRNO_UNCHANGED)
+	test_single_errno (test_name, errno_value, 0, "unchanged");
+      if (exceptions & ERRNO_EDOM)
+	test_single_errno (test_name, errno_value, EDOM, "EDOM");
+      if (exceptions & ERRNO_ERANGE)
+	test_single_errno (test_name, errno_value, ERANGE, "ERANGE");
+    }
+}
+
+/* Returns the number of ulps that GIVEN is away from EXPECTED.  */
+#define ULPDIFF(given, expected) \
+	(FUNC(fabs) ((given) - (expected)) / ulp (expected))
+
+/* Returns the size of an ulp for VALUE.  */
+static FLOAT
+ulp (FLOAT value)
+{
+  FLOAT ulp;
+
+  switch (fpclassify (value))
+    {
+      case FP_ZERO:
+	/* We compute the distance to the next FP which is the same as the
+	   value of the smallest subnormal number. Previously we used
+	   2^-(MANT_DIG - 1) which is too large a value to be useful. Note that we
+	   can't use ilogb(0), since that isn't a valid thing to do. As a point
+	   of comparison Java's ulp returns the next normal value e.g.
+	   2^(1 - MAX_EXP) for ulp(0), but that is not what we want for
+	   glibc.  */
+	/* Fall through...  */
+      case FP_SUBNORMAL:
+        /* The next closest subnormal value is a constant distance away.  */
+	ulp = FUNC(ldexp) (1.0, MIN_EXP - MANT_DIG);
+	break;
+
+      case FP_NORMAL:
+	ulp = FUNC(ldexp) (1.0, FUNC(ilogb) (value) - MANT_DIG + 1);
+	break;
+
+      default:
+	/* It should never happen. */
+	abort ();
+	break;
+    }
+  return ulp;
+}
 
 static void
 check_float_internal (const char *test_name, FLOAT computed, FLOAT expected,
-		      FLOAT max_ulp, int xfail, int exceptions,
-		      FLOAT *curr_max_error)
+		      int exceptions,
+		      FLOAT *curr_max_error, FLOAT max_ulp)
 {
   int ok = 0;
   int print_diff = 0;
   FLOAT diff = 0;
-  FLOAT ulp = 0;
+  FLOAT ulps = 0;
+  int errno_value = errno;
 
   test_exceptions (test_name, exceptions);
-  if (isnan (computed) && isnan (expected))
-    ok = 1;
+  test_errno (test_name, errno_value, exceptions);
+  if (exceptions & IGNORE_RESULT)
+    goto out;
+  if (issignaling (computed) && issignaling (expected))
+    {
+      if ((exceptions & TEST_NAN_SIGN) != 0
+	  && signbit (computed) != signbit (expected))
+	{
+	  ok = 0;
+	  printf ("signaling NaN has wrong sign.\n");
+	}
+      else
+	ok = 1;
+    }
+  else if (issignaling (computed) || issignaling (expected))
+    ok = 0;
+  else if (isnan (computed) && isnan (expected))
+    {
+      if ((exceptions & TEST_NAN_SIGN) != 0
+	  && signbit (computed) != signbit (expected))
+	{
+	  ok = 0;
+	  printf ("quiet NaN has wrong sign.\n");
+	}
+      else
+	ok = 1;
+    }
   else if (isinf (computed) && isinf (expected))
     {
       /* Test for sign of infinities.  */
@@ -483,115 +822,119 @@ check_float_internal (const char *test_name, FLOAT computed, FLOAT expected,
       else
 	ok = 1;
     }
-  /* Don't calc ulp for NaNs or infinities.  */
-  else if (isinf (computed) || isnan (computed) || isinf (expected) || isnan (expected))
+  /* Don't calculate ULPs for infinities or any kind of NaNs.  */
+  else if (isinf (computed) || isnan (computed)
+	   || isinf (expected) || isnan (expected))
     ok = 0;
   else
     {
       diff = FUNC(fabs) (computed - expected);
-      /* ilogb (0) isn't allowed.  */
-      if (expected == 0.0)
-	ulp = diff / FUNC(ldexp) (1.0, - MANT_DIG);
-      else
-	ulp = diff / FUNC(ldexp) (1.0, FUNC(ilogb) (expected) - MANT_DIG);
-      set_max_error (ulp, curr_max_error);
+      ulps = ULPDIFF (computed, expected);
+      set_max_error (ulps, curr_max_error);
       print_diff = 1;
       if ((exceptions & IGNORE_ZERO_INF_SIGN) == 0
 	  && computed == 0.0 && expected == 0.0
 	  && signbit(computed) != signbit (expected))
 	ok = 0;
-      else if (ulp <= 0.5 || (ulp <= max_ulp && !ignore_max_ulp))
+      else if (ulps <= max_ulp && !ignore_max_ulp)
 	ok = 1;
       else
-	{
-	  ok = 0;
-	  print_ulps (test_name, ulp);
-	}
-
+	ok = 0;
     }
-  if (print_screen (ok, xfail))
+  if (print_screen (ok))
     {
       if (!ok)
 	printf ("Failure: ");
       printf ("Test: %s\n", test_name);
       printf ("Result:\n");
-      printf (" is:         % .20" PRINTF_EXPR "  % .20" PRINTF_XEXPR "\n",
-	      computed, computed);
-      printf (" should be:  % .20" PRINTF_EXPR "  % .20" PRINTF_XEXPR "\n",
-	      expected, expected);
+      printf (" is:         ");
+      print_float (computed);
+      printf (" should be:  ");
+      print_float (expected);
       if (print_diff)
 	{
-	  printf (" difference: % .20" PRINTF_EXPR "  % .20" PRINTF_XEXPR
-		  "\n", diff, diff);
-	  printf (" ulp       : % .4" PRINTF_NEXPR "\n", ulp);
-	  printf (" max.ulp   : % .4" PRINTF_NEXPR "\n", max_ulp);
+	  char dstrn[FSTR_MAX], dstrx[FSTR_MAX];
+	  char ustrn[FSTR_MAX], mustrn[FSTR_MAX];
+	  FTOSTR (dstrn, FSTR_MAX, "% .*" PRINTF_EXPR,
+		  TYPE_DECIMAL_DIG - 1, diff);
+	  FTOSTR (dstrx, FSTR_MAX, "% .*" PRINTF_XEXPR,
+		  TYPE_HEX_DIG - 1, diff);
+	  FTOSTR (ustrn, FSTR_MAX, "% .4" PRINTF_NEXPR, ulps);
+	  FTOSTR (mustrn, FSTR_MAX, "% .4" PRINTF_NEXPR, max_ulp);
+	  printf (" difference: %s  %s\n", dstrn, dstrx);
+	  printf (" ulp       : %s\n", ustrn);
+	  printf (" max.ulp   : %s\n", mustrn);
 	}
     }
-  update_stats (ok, xfail);
+  update_stats (ok);
 
+ out:
   fpstack_test (test_name);
+  errno = 0;
 }
 
 
 static void
 check_float (const char *test_name, FLOAT computed, FLOAT expected,
-	     FLOAT max_ulp, int xfail, int exceptions)
+	     int exceptions)
 {
-  check_float_internal (test_name, computed, expected, max_ulp, xfail,
-			exceptions, &max_error);
+  check_float_internal (test_name, computed, expected,
+			exceptions, &max_error, prev_max_error);
 }
 
 
 static void
 check_complex (const char *test_name, __complex__ FLOAT computed,
 	       __complex__ FLOAT expected,
-	       __complex__ FLOAT max_ulp, __complex__ int xfail,
 	       int exception)
 {
-  FLOAT part_comp, part_exp, part_max_ulp;
-  int part_xfail;
-  char str[200];
+  FLOAT part_comp, part_exp;
+  char *str;
+
+  if (asprintf (&str, "Real part of: %s", test_name) == -1)
+    abort ();
 
-  sprintf (str, "Real part of: %s", test_name);
   part_comp = __real__ computed;
   part_exp = __real__ expected;
-  part_max_ulp = __real__ max_ulp;
-  part_xfail = __real__ xfail;
 
-  check_float_internal (str, part_comp, part_exp, part_max_ulp, part_xfail,
-			exception, &real_max_error);
+  check_float_internal (str, part_comp, part_exp,
+			exception, &real_max_error, prev_real_max_error);
+  free (str);
+
+  if (asprintf (&str, "Imaginary part of: %s", test_name) == -1)
+    abort ();
 
-  sprintf (str, "Imaginary part of: %s", test_name);
   part_comp = __imag__ computed;
   part_exp = __imag__ expected;
-  part_max_ulp = __imag__ max_ulp;
-  part_xfail = __imag__ xfail;
 
-  /* Don't check again for exceptions, just pass through the
-     zero/inf sign test.  */
-  check_float_internal (str, part_comp, part_exp, part_max_ulp, part_xfail,
-			exception & IGNORE_ZERO_INF_SIGN,
-			&imag_max_error);
+  /* Don't check again for exceptions or errno, just pass through the
+     other relevant flags.  */
+  check_float_internal (str, part_comp, part_exp,
+			exception & (IGNORE_ZERO_INF_SIGN
+				     | TEST_NAN_SIGN
+				     | IGNORE_RESULT),
+			&imag_max_error, prev_imag_max_error);
+  free (str);
 }
 
 
 /* Check that computed and expected values are equal (int values).  */
 static void
-check_int (const char *test_name, int computed, int expected, int max_ulp,
-	   int xfail, int exceptions)
+check_int (const char *test_name, int computed, int expected,
+	   int exceptions)
 {
-  int diff = computed - expected;
   int ok = 0;
+  int errno_value = errno;
 
   test_exceptions (test_name, exceptions);
+  test_e