math: Refactor how to use libm-test-ulps

The current approach tracks math maximum supported errors by explicitly
setting them per function and architecture. On newer implementations or
new compiler versions, the file is updated with newer values if it
shows higher results. The idea is to track the maximum known error, to
update the manual with the obtained values.

The constant libm-test-ulps shows little value, where it is usually a
mechanical change done by the maintainer, for past releases it is
usually ignored whether the ulp change resulted from a compiler
regression, and the math tests already have a maximum ulp error that
triggers a regression.

It was shown by a recent update after the new acosf [1] implementation
that is correctly rounded, where the libm-test-ulps was indeed from a
compiler issue.

This patch removes all arch-specific libm-test-ulps, adds system generic
libm-test-ulps where applicable, and changes its semantics. The generic
files now track specific implementation constraints, like if it is
expected to be correctly rounded, or if the system-specific has
different error expectations.

Now multiple libm-test-ulps can be defined, and system-specific
overrides generic implementation.  This is for the case where
arch-specific implementation might show worse precision than generic
implementation, for instance, the cbrtf on i686.

Regressions are only reported if the implementation shows larger errors
than 9 ulps (13 for IBM long double) unless it is overridden by
libm-test-ulps and the maximum error is not printed at the end of tests.
The regen-ulps rule is also removed since it does not make sense to
update the libm-test-ulps automatically.

The manual error table is also removed, Paul Zimmermann and others have
been tracking libm precision with a more comprehensive analysis for some
releases; so link to his work instead.

[1] https://sourceware.org/git/?p=glibc.git;a=commit;h=9cc9f8e11e8fb8f54f1e84d9f024917634a78201

4 files changed, 104 insertions, 135 deletions

diff --git a/math/Makefile b/math/Makefile
index 7aabf35..b1d7c76 100644
--- a/math/Makefile
+++ b/math/Makefile
@@ -621,11 +621,13 @@ libm-vec-test-wrappers = $(addsuffix -wrappers, $(libm-vec-tests))
 test-extras += $(libm-vec-test-wrappers)
 extra-test-objs += $(addsuffix .o, $(libm-vec-test-wrappers))
 
-ulps-file = $(firstword $(wildcard $(sysdirs:%=%/libm-test-ulps)))
+ulps-file = $(wildcard $(sysdirs:%=%/libm-test-ulps))
+
+space := $(subst ,, )
 
 $(objpfx)libm-test-ulps.h: $(ulps-file) gen-libm-test.py
 	$(make-target-directory)
-	$(PYTHON) gen-libm-test.py -u $< -H $@
+	$(PYTHON) gen-libm-test.py -u $(subst $(space),:,$(ulps-file)) -H $@
 
 libm-test-funcs-auto = \
   acos \
@@ -1301,29 +1303,6 @@ endef
 object-suffixes-left := $(test-types)
 include $(o-iterator)
 
-# Run the math programs to automatically generate ULPs files.
-.PHONY: regen-ulps
-
-run-regen-ulps = $(test-wrapper-env) \
-	    $(run-program-env) \
-	    $($*-ENV) $(rtld-prefix) $(objpfx)$${run}
-
-regen-ulps: $(addprefix $(objpfx),$(libm-tests))
-	rm -f $(objpfx)ULPs; rm -f $(objpfx)NewUlps; \
-	cp $(ulps-file) $(objpfx)libm-test-ulps; \
-	for run in $(libm-tests); do \
-	  echo "Regenerating ULPs for $${run}"; \
-	  $(run-regen-ulps) -u -o $(objpfx); \
-	  cat $(objpfx)ULPs >> $(objpfx)libm-test-ulps; \
-	  rm $(objpfx)ULPs; \
-	done; \
-	$(PYTHON) gen-libm-test.py -n $(objpfx)NewUlps \
-				   -u $(objpfx)libm-test-ulps; \
-	echo "Automatic regeneration of ULPs complete."; \
-	echo "Difference between the current baseline and the new baseline is:";\
-	diff -urN $(ulps-file) $(objpfx)NewUlps; \
-	echo "Copy $(objpfx)NewUlps to $(ulps-file) (relative to source)."
-
 # The generated sysd-rules file defines rules like this for sources
 # coming from sysdeps/ directories.  These rules find the generic sources.
 define o-iterator-doit
diff --git a/math/README.libm-test b/math/README.libm-test
index f32c629..84b2979 100644
--- a/math/README.libm-test
+++ b/math/README.libm-test
@@ -30,8 +30,8 @@ libm.  These files, along with generated files named
 cases for each math function that are specific to the target platform
 but still independent of the real floating type.  The results of the
 processing are "libm-test-<func>.c" and a file "libm-test-ulps.h" with
-platform specific deltas by which the actual math function results may
-deviate from the expected results and still be considered correct.
+specific math results that can be either generic for the floating
+type or platform specific.
 
 The test drivers "test-double-<func>.c", "test-float-<func>.c", and
 "test-ldouble-<func>.c", generated by the Makefile, test the normal
@@ -42,13 +42,11 @@ of macros just before including the generic "libm-test.c" file. Each
 driver is compiled into a single executable test program with the
 corresponding name.
 
-As mentioned above, the "gen-libm-test.py" script looks for a file
-named "libm-test-ulps" in the platform specific sysdep directory (or
-its fpu or nofpu subdirectory) and for each variant (real floating
-type and rounding mode) of every tested function reads from it the
-maximum difference expressed as Units of Least Precision (ULP) the
-actual result of the function may deviate from the expected result
-before it's considered incorrect.
+The math tests do not report up to 9 Units of Least Precision (ULP)
+(13 for IBM long double format) difference between the obtained
+result and the expected one as a regression.  The "gen-libm-test.py"
+script looks for files named "libm-test-ulps" in the sysdep directories
+to generate the "libm-test-ulps.h" file.
 
 The "auto-libm-test-out-<func>" files contain sets of test cases to
 exercise, the conditions under which to exercise each, and the
@@ -57,29 +55,39 @@ expected results.  The files are generated by the
 the comments in gen-auto-libm-tests.c for details about the content
 and format of the -in and -out files.
 
-How can I generate "libm-test-ulps"?
+How can I use "libm-test-ulps"?
 ====================================
 
-To automatically generate a new "libm-test-ulps" run "make regen-ulps".
-This generates the file "math/NewUlps" in the build directory.  The file
-contains the sorted results of all the tests.  You can use the "NewUlps"
-file as the machine's updated "libm-test-ulps" file.  Copy "NewUlps" to
-"libm-test-ulps" in the appropriate machine sysdep directory.  Verify
-the changes, post your patch, and check it in after review.
-
-To manually generate a new "libm-test-ulps" file, first remove "ULPs"
-file in the current directory, then you can execute for example:
-    ./testrun.sh math/test-double -u --ignore-max-ulp=yes
-This generates a file "ULPs" with all double ULPs in it, ignoring any
-previously calculated ULPs, and running with the newly built dynamic
-loader and math library (assumes you didn't install your build).  Now
-generate the ULPs for all other formats, the tests will be appending the
-data to the "ULPs" file.  As final step run "gen-libm-test.py" with the
-file as input and ask to generate a pretty printed output in the file
-"NewUlps":
-  gen-libm-test.py -u ULPs -n NewUlps
-Copy "NewUlps" to "libm-test-ulps" in the appropriate machine sysdep
-directory.
+A "libm-test-ulps" is required only to test for extra constraints in
+the math tests.  The file contains lines for maximal errors of single
+functions, like:
+
+Function "yn":
+float: 2
+double: 6
+
+It means that if the "yn" shows error larger than 2 ULP for float
+or 6 ULP for double, the related test for "symbol" will fail.  It can
+be useful to check for correctly rounded implementation, where the
+expected ULP is 0.
+
+The function is tested with default FE_TONEAREST rounding mode.  To
+check with a different one, the function definition name should be
+prepended with an underline plus the rounding mode 'downward' (FE_DOWNWARD),
+'towardzero' (FE_TOWARDZERO), or 'upward' (FE_UPWARD).  For instance,
+
+Function "yn_downward":
+float: 3
+double: 7
+
+It means that 'yn' will be checked with FE_DOWNWARD rounding mode
+and any error larger than 3 ULPs for float or 7 ULPs for double will be
+reported as a regression.
+
+The keywords are float, double, ldouble, and float128.
+
+Also, multiple "libm-test-ulps" can be added, "gen-libm-test.py" will
+merge the input in only one table.
 
 Note that the test drivers have an option "-u" to output an unsorted
 list of all epsilons that the functions have.  The output can be read
@@ -87,16 +95,6 @@ in directly but it's better to pretty print it first.
 "gen-libm-test.py" has an option to generate a pretty-printed and
 sorted new ULPs file from the output of the test drivers.
 
-Contents of libm-test-ulps
-==========================
-
-Since libm-test-ulps can be generated automatically, just a few notes.
-The file contains lines for maximal errors of single functions, like:
-
-Function "yn":
-double: 6
-
-The keywords are float, double, and ldouble.
 
 Adding tests to libm-test-<func>.inc
 ====================================
diff --git a/math/gen-libm-test.py b/math/gen-libm-test.py
index 6c48372..a76fe87 100755
--- a/math/gen-libm-test.py
+++ b/math/gen-libm-test.py
@@ -121,19 +121,20 @@ class Ulps(object):
         """Initialize an Ulps object."""
         # normal[function][float_type] is the ulps value, and likewise
         # for real and imag.
-        self.normal = defaultdict(lambda: defaultdict(lambda: 0))
-        self.real = defaultdict(lambda: defaultdict(lambda: 0))
-        self.imag = defaultdict(lambda: defaultdict(lambda: 0))
+        self.normal = defaultdict(lambda: defaultdict(lambda: -1))
+        self.real = defaultdict(lambda: defaultdict(lambda: -1))
+        self.imag = defaultdict(lambda: defaultdict(lambda: -1))
         # List of ulps kinds, in the order in which they appear in
         # sorted ulps files.
         self.ulps_kinds = (('Real part of ', self.real),
                            ('Imaginary part of ', self.imag),
                            ('', self.normal))
+        self.ulps_file = []
         self
 
     def read(self, ulps_file):
         """Read ulps from a file into an Ulps object."""
-        self.ulps_file = ulps_file
+        self.ulps_file.append(ulps_file)
         with open(ulps_file, 'r') as f:
             ulps_dict = None
             ulps_fn = None
@@ -165,10 +166,7 @@ class Ulps(object):
                     if line_first not in ALL_FLOATS:
                         raise ValueError('bad ulps line: %s' % line)
                     ulps_val = int(line_second)
-                    if ulps_val > 0:
-                        ulps_dict[ulps_fn][line_first] = max(
-                            ulps_dict[ulps_fn][line_first],
-                            ulps_val)
+                    ulps_dict[ulps_fn][line_first] = ulps_val
 
     def all_functions(self):
         """Return the set of functions with ulps and whether they are
@@ -181,27 +179,6 @@ class Ulps(object):
                 complex[f] = True if k_prefix else False
         return funcs, complex
 
-    def write(self, ulps_file):
-        """Write ulps back out as a sorted ulps file."""
-        # Output is sorted first by function name, then by (real,
-        # imag, normal), then by float type.
-        out_data = {}
-        for order, (prefix, d) in enumerate(self.ulps_kinds):
-            for fn in d.keys():
-                fn_data = ['%s: %d' % (f, d[fn][f])
-                           for f in sorted(d[fn].keys())]
-                fn_text = 'Function: %s"%s":\n%s' % (prefix, fn,
-                                                     '\n'.join(fn_data))
-                out_data[(fn, order)] = fn_text
-        out_list = [out_data[fn_order] for fn_order in sorted(out_data.keys())]
-        out_text = ('# Begin of automatic generation\n\n'
-                    '# Maximal error of functions:\n'
-                    '%s\n\n'
-                    '# end of automatic generation\n'
-                    % '\n\n'.join(out_list))
-        with open(ulps_file, 'w') as f:
-            f.write(out_text)
-
     @staticmethod
     def ulps_table(name, ulps_dict):
         """Return text of a C table of ulps."""
@@ -227,7 +204,7 @@ class Ulps(object):
                          '  const char *name;\n'
                          '  FLOAT max_ulp[%d];\n'
                          '};'
-                         % (self.ulps_file, len(ALL_FLOATS)))
+                         % (', '.join(self.ulps_file), len(ALL_FLOATS)))
         macro_list = []
         for i, f in enumerate(ALL_FLOATS):
             if f.startswith('i'):
@@ -251,18 +228,6 @@ class Ulps(object):
             f.write(header_text)
 
 
-def read_all_ulps(srcdir):
-    """Read all platforms' libm-test-ulps files."""
-    all_ulps = {}
-    for dirpath, dirnames, filenames in os.walk(srcdir):
-        if 'libm-test-ulps' in filenames:
-            with open(os.path.join(dirpath, 'libm-test-ulps-name')) as f:
-                name = f.read().rstrip()
-            all_ulps[name] = Ulps()
-            all_ulps[name].read(os.path.join(dirpath, 'libm-test-ulps'))
-    return all_ulps
-
-
 def read_auto_tests(test_file):
     """Read tests from auto-libm-test-out-<function> (possibly None)."""
     auto_tests = defaultdict(lambda: defaultdict(dict))
@@ -653,12 +618,8 @@ def main():
                         help='input file with automatically generated tests')
     parser.add_argument('-c', dest='inc_input', metavar='FILE',
                         help='input file .inc file with tests')
-    parser.add_argument('-u', dest='ulps_file', metavar='FILE',
-                        help='input file with ulps')
-    parser.add_argument('-s', dest='srcdir', metavar='DIR',
-                        help='input source directory with all ulps')
-    parser.add_argument('-n', dest='ulps_output', metavar='FILE',
-                        help='generate sorted ulps file FILE')
+    parser.add_argument('-u', dest='ulps_file', metavar='list',
+                        help='input files with ulps (multiple input separated by colon')
     parser.add_argument('-C', dest='c_output', metavar='FILE',
                         help='generate output C file FILE from .inc file')
     parser.add_argument('-H', dest='ulps_header', metavar='FILE',
@@ -668,12 +629,11 @@ def main():
     args = parser.parse_args()
     ulps = Ulps()
     if args.ulps_file is not None:
-        ulps.read(args.ulps_file)
+        # Iterate in reverse order so arch specific definitions can override
+        # the generic ones.
+        for ulp_file in reversed(args.ulps_file.split(':')):
+            ulps.read(ulp_file)
     auto_tests = read_auto_tests(args.auto_input)
-    if args.srcdir is not None:
-        all_ulps = read_all_ulps(args.srcdir)
-    if args.ulps_output is not None:
-        ulps.write(args.ulps_output)
     if args.ulps_header is not None:
         ulps.write_header(args.ulps_header)
     if args.c_output is not None:
diff --git a/math/libm-test-support.c b/math/libm-test-support.c
index 3a7b096..31ae990 100644
--- a/math/libm-test-support.c
+++ b/math/libm-test-support.c
@@ -109,10 +109,12 @@
 
 #include "libm-test-support.h"
 
+#include <array_length.h>
 #include <argp.h>
 #include <errno.h>
 #include <string.h>
 #include <assert.h>
+#include <stdbool.h>
 
 /* This header defines func_ulps, func_real_ulps and func_imag_ulps
    arrays.  */
@@ -146,6 +148,8 @@ static FLOAT prev_max_error, prev_real_max_error, prev_imag_max_error;
 
 static FLOAT max_valid_error;
 
+static bool is_complex = false;
+
 /* 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
@@ -186,6 +190,20 @@ fmt_ftostr (char *dest, size_t size, int precision, const char *conversion,
   FTOSTR (dest, size, format, value);
 }
 
+
+static FLOAT
+default_max_valid_error (int exact, int testing_ibm128)
+{
+  if (testing_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.  */
+    return max_valid_error = exact ? 3 : 16;
+  else
+    return max_valid_error = exact ? 0 : 9;
+}
+
 /* 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.  */
@@ -204,14 +222,17 @@ static const int ulp_idx = ULP_IDX;
    no ulps listed.  */
 
 static FLOAT
-find_ulps (const char *name, const struct ulp_data *data, size_t nmemb)
+find_ulps (const char *name, const struct ulp_data *data, size_t nmemb,
+	   int exact, int testing_ibm128)
 {
   const struct ulp_data *entry = bsearch (name, data, nmemb, sizeof (*data),
 					  compare_ulp_data);
   if (entry == NULL)
-    return 0;
+    return default_max_valid_error (exact, testing_ibm128);
   else
-    return entry->max_ulp[ulp_idx];
+    return entry->max_ulp[ulp_idx] == -1
+	   ? default_max_valid_error (exact, testing_ibm128)
+	   : entry->max_ulp[ulp_idx];
 }
 
 void
@@ -221,22 +242,15 @@ init_max_error (const char *name, int exact, int testing_ibm128)
   real_max_error = 0;
   imag_max_error = 0;
   test_ibm128 = testing_ibm128;
-  prev_max_error = find_ulps (name, func_ulps,
-			      sizeof (func_ulps) / sizeof (func_ulps[0]));
+  prev_max_error = find_ulps (name, func_ulps, array_length (func_ulps),
+			      exact, testing_ibm128);
   prev_real_max_error = find_ulps (name, func_real_ulps,
-				   (sizeof (func_real_ulps)
-				    / sizeof (func_real_ulps[0])));
+				   array_length (func_real_ulps), exact,
+				   testing_ibm128);
   prev_imag_max_error = find_ulps (name, func_imag_ulps,
-				   (sizeof (func_imag_ulps)
-				    / sizeof (func_imag_ulps[0])));
-  if (testing_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 : 16;
-  else
-    max_valid_error = exact ? 0 : 9;
+				   array_length (func_imag_ulps),
+				   exact, testing_ibm128);
+  max_valid_error = default_max_valid_error (exact, testing_ibm128);
   prev_max_error = (prev_max_error <= max_valid_error
 		    ? prev_max_error
 		    : max_valid_error);
@@ -476,6 +490,8 @@ check_complex_max_error (const char *func_name)
 
   update_stats (ok, TEST_MAXERROR);
   print_test_end (thisTest, func_name, TEST_MAXERROR);
+
+  is_complex = true;
 }
 
 
@@ -1322,6 +1338,22 @@ libm_test_finish (void)
     fclose (ulps_file);
 
   printf ("\nTest suite completed:\n");
+  printf ("  Maximum error found of ");
+  if (is_complex)
+    {
+      char rmestr[FSTR_MAX];
+      char imestr[FSTR_MAX];
+      FTOSTR (rmestr, FSTR_MAX, "%.0f", FUNC (ceil) (real_max_error));
+      FTOSTR (imestr, FSTR_MAX, "%.0f", FUNC (ceil) (imag_max_error));
+      printf ("`%s' ulp for real part and `%s' ulp for imaginary part\n",
+	      rmestr, imestr);
+    }
+  else
+    {
+      char mestr[FSTR_MAX];
+      FTOSTR (mestr, FSTR_MAX, "%.0f", FUNC (ceil) (max_error));
+      printf ("`%s' ulp\n", mestr);
+    }
   printf ("  %d max error test cases,\n", noMaxErrorTests);
   printf ("  %d input tests,\n", noTests);
   printf ("  - with %d tests for exception flags,\n", noExcTests);