Start generating libm tests automatically with MPFR.

6 files changed, 2110 insertions, 19 deletions

diff --git a/math/Makefile b/math/Makefile
index fcccab2051..d178789fb5 100644
--- a/math/Makefile
+++ b/math/Makefile
@@ -112,7 +112,8 @@ ulps-file = $(firstword $(wildcard $(sysdirs:%=%/libm-test-ulps)))
 
 $(addprefix $(objpfx), $(libm-tests-generated)): $(objpfx)libm-test.stmp
 
-$(objpfx)libm-test.stmp: $(ulps-file) libm-test.inc gen-libm-test.pl
+$(objpfx)libm-test.stmp: $(ulps-file) libm-test.inc gen-libm-test.pl \
+			 auto-libm-test-out
 	$(make-target-directory)
 	$(PERL) gen-libm-test.pl -u $< -o "$(objpfx)"
 	@echo > $@
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
new file mode 100644
index 0000000000..5dc790c692
--- /dev/null
+++ b/math/auto-libm-test-in
@@ -0,0 +1,27 @@
+# libm test inputs for gen-auto-libm-tests.c.
+# Copyright (C) 1997-2013 Free Software Foundation, Inc.
+# This file is part of the GNU C Library.
+#
+# The GNU C Library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# The GNU C Library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with the GNU C Library; if not, see
+# <http://www.gnu.org/licenses/>.  */
+
+sqrt 0
+sqrt -0
+sqrt 2209
+sqrt 4
+sqrt 2
+sqrt 0.25
+sqrt 6642.25
+sqrt 15190.5625
+sqrt 0.75
diff --git a/math/auto-libm-test-out b/math/auto-libm-test-out
new file mode 100644
index 0000000000..bdb6bdc8b8
--- /dev/null
+++ b/math/auto-libm-test-out
@@ -0,0 +1,225 @@
+sqrt 0
+= sqrt downward flt-32 0x0p+0f : 0x0p+0f :
+= sqrt tonearest flt-32 0x0p+0f : 0x0p+0f :
+= sqrt towardzero flt-32 0x0p+0f : 0x0p+0f :
+= sqrt upward flt-32 0x0p+0f : 0x0p+0f :
+= sqrt downward dbl-64 0x0p+0 : 0x0p+0 :
+= sqrt tonearest dbl-64 0x0p+0 : 0x0p+0 :
+= sqrt towardzero dbl-64 0x0p+0 : 0x0p+0 :
+= sqrt upward dbl-64 0x0p+0 : 0x0p+0 :
+= sqrt downward ldbl-96-intel 0x0p+0L : 0x0p+0L :
+= sqrt tonearest ldbl-96-intel 0x0p+0L : 0x0p+0L :
+= sqrt towardzero ldbl-96-intel 0x0p+0L : 0x0p+0L :
+= sqrt upward ldbl-96-intel 0x0p+0L : 0x0p+0L :
+= sqrt downward ldbl-96-m68k 0x0p+0L : 0x0p+0L :
+= sqrt tonearest ldbl-96-m68k 0x0p+0L : 0x0p+0L :
+= sqrt towardzero ldbl-96-m68k 0x0p+0L : 0x0p+0L :
+= sqrt upward ldbl-96-m68k 0x0p+0L : 0x0p+0L :
+= sqrt downward ldbl-128 0x0p+0L : 0x0p+0L :
+= sqrt tonearest ldbl-128 0x0p+0L : 0x0p+0L :
+= sqrt towardzero ldbl-128 0x0p+0L : 0x0p+0L :
+= sqrt upward ldbl-128 0x0p+0L : 0x0p+0L :
+= sqrt downward ldbl-128ibm 0x0p+0L : 0x0p+0L :
+= sqrt tonearest ldbl-128ibm 0x0p+0L : 0x0p+0L :
+= sqrt towardzero ldbl-128ibm 0x0p+0L : 0x0p+0L :
+= sqrt upward ldbl-128ibm 0x0p+0L : 0x0p+0L :
+sqrt -0
+= sqrt downward flt-32 -0x0p+0f : -0x0p+0f :
+= sqrt tonearest flt-32 -0x0p+0f : -0x0p+0f :
+= sqrt towardzero flt-32 -0x0p+0f : -0x0p+0f :
+= sqrt upward flt-32 -0x0p+0f : -0x0p+0f :
+= sqrt downward dbl-64 -0x0p+0 : -0x0p+0 :
+= sqrt tonearest dbl-64 -0x0p+0 : -0x0p+0 :
+= sqrt towardzero dbl-64 -0x0p+0 : -0x0p+0 :
+= sqrt upward dbl-64 -0x0p+0 : -0x0p+0 :
+= sqrt downward ldbl-96-intel -0x0p+0L : -0x0p+0L :
+= sqrt tonearest ldbl-96-intel -0x0p+0L : -0x0p+0L :
+= sqrt towardzero ldbl-96-intel -0x0p+0L : -0x0p+0L :
+= sqrt upward ldbl-96-intel -0x0p+0L : -0x0p+0L :
+= sqrt downward ldbl-96-m68k -0x0p+0L : -0x0p+0L :
+= sqrt tonearest ldbl-96-m68k -0x0p+0L : -0x0p+0L :
+= sqrt towardzero ldbl-96-m68k -0x0p+0L : -0x0p+0L :
+= sqrt upward ldbl-96-m68k -0x0p+0L : -0x0p+0L :
+= sqrt downward ldbl-128 -0x0p+0L : -0x0p+0L :
+= sqrt tonearest ldbl-128 -0x0p+0L : -0x0p+0L :
+= sqrt towardzero ldbl-128 -0x0p+0L : -0x0p+0L :
+= sqrt upward ldbl-128 -0x0p+0L : -0x0p+0L :
+= sqrt downward ldbl-128ibm -0x0p+0L : -0x0p+0L :
+= sqrt tonearest ldbl-128ibm -0x0p+0L : -0x0p+0L :
+= sqrt towardzero ldbl-128ibm -0x0p+0L : -0x0p+0L :
+= sqrt upward ldbl-128ibm -0x0p+0L : -0x0p+0L :
+sqrt 2209
+= sqrt downward flt-32 0x8.a1p+8f : 0x2.fp+4f :
+= sqrt tonearest flt-32 0x8.a1p+8f : 0x2.fp+4f :
+= sqrt towardzero flt-32 0x8.a1p+8f : 0x2.fp+4f :
+= sqrt upward flt-32 0x8.a1p+8f : 0x2.fp+4f :
+= sqrt downward dbl-64 0x8.a1p+8 : 0x2.fp+4 :
+= sqrt tonearest dbl-64 0x8.a1p+8 : 0x2.fp+4 :
+= sqrt towardzero dbl-64 0x8.a1p+8 : 0x2.fp+4 :
+= sqrt upward dbl-64 0x8.a1p+8 : 0x2.fp+4 :
+= sqrt downward ldbl-96-intel 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt tonearest ldbl-96-intel 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt towardzero ldbl-96-intel 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt upward ldbl-96-intel 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt downward ldbl-96-m68k 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt tonearest ldbl-96-m68k 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt towardzero ldbl-96-m68k 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt upward ldbl-96-m68k 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt downward ldbl-128 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt tonearest ldbl-128 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt towardzero ldbl-128 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt upward ldbl-128 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt downward ldbl-128ibm 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt tonearest ldbl-128ibm 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt towardzero ldbl-128ibm 0x8.a1p+8L : 0x2.fp+4L :
+= sqrt upward ldbl-128ibm 0x8.a1p+8L : 0x2.fp+4L :
+sqrt 4
+= sqrt downward flt-32 0x4p+0f : 0x2p+0f :
+= sqrt tonearest flt-32 0x4p+0f : 0x2p+0f :
+= sqrt towardzero flt-32 0x4p+0f : 0x2p+0f :
+= sqrt upward flt-32 0x4p+0f : 0x2p+0f :
+= sqrt downward dbl-64 0x4p+0 : 0x2p+0 :
+= sqrt tonearest dbl-64 0x4p+0 : 0x2p+0 :
+= sqrt towardzero dbl-64 0x4p+0 : 0x2p+0 :
+= sqrt upward dbl-64 0x4p+0 : 0x2p+0 :
+= sqrt downward ldbl-96-intel 0x4p+0L : 0x2p+0L :
+= sqrt tonearest ldbl-96-intel 0x4p+0L : 0x2p+0L :
+= sqrt towardzero ldbl-96-intel 0x4p+0L : 0x2p+0L :
+= sqrt upward ldbl-96-intel 0x4p+0L : 0x2p+0L :
+= sqrt downward ldbl-96-m68k 0x4p+0L : 0x2p+0L :
+= sqrt tonearest ldbl-96-m68k 0x4p+0L : 0x2p+0L :
+= sqrt towardzero ldbl-96-m68k 0x4p+0L : 0x2p+0L :
+= sqrt upward ldbl-96-m68k 0x4p+0L : 0x2p+0L :
+= sqrt downward ldbl-128 0x4p+0L : 0x2p+0L :
+= sqrt tonearest ldbl-128 0x4p+0L : 0x2p+0L :
+= sqrt towardzero ldbl-128 0x4p+0L : 0x2p+0L :
+= sqrt upward ldbl-128 0x4p+0L : 0x2p+0L :
+= sqrt downward ldbl-128ibm 0x4p+0L : 0x2p+0L :
+= sqrt tonearest ldbl-128ibm 0x4p+0L : 0x2p+0L :
+= sqrt towardzero ldbl-128ibm 0x4p+0L : 0x2p+0L :
+= sqrt upward ldbl-128ibm 0x4p+0L : 0x2p+0L :
+sqrt 2
+= sqrt downward flt-32 0x2p+0f : 0x1.6a09e6p+0f : inexact
+= sqrt tonearest flt-32 0x2p+0f : 0x1.6a09e6p+0f : inexact
+= sqrt towardzero flt-32 0x2p+0f : 0x1.6a09e6p+0f : inexact
+= sqrt upward flt-32 0x2p+0f : 0x1.6a09e8p+0f : inexact
+= sqrt downward dbl-64 0x2p+0 : 0x1.6a09e667f3bccp+0 : inexact
+= sqrt tonearest dbl-64 0x2p+0 : 0x1.6a09e667f3bcdp+0 : inexact
+= sqrt towardzero dbl-64 0x2p+0 : 0x1.6a09e667f3bccp+0 : inexact
+= sqrt upward dbl-64 0x2p+0 : 0x1.6a09e667f3bcdp+0 : inexact
+= sqrt downward ldbl-96-intel 0x2p+0L : 0x1.6a09e667f3bcc908p+0L : inexact
+= sqrt tonearest ldbl-96-intel 0x2p+0L : 0x1.6a09e667f3bcc908p+0L : inexact
+= sqrt towardzero ldbl-96-intel 0x2p+0L : 0x1.6a09e667f3bcc908p+0L : inexact
+= sqrt upward ldbl-96-intel 0x2p+0L : 0x1.6a09e667f3bcc90ap+0L : inexact
+= sqrt downward ldbl-96-m68k 0x2p+0L : 0x1.6a09e667f3bcc908p+0L : inexact
+= sqrt tonearest ldbl-96-m68k 0x2p+0L : 0x1.6a09e667f3bcc908p+0L : inexact
+= sqrt towardzero ldbl-96-m68k 0x2p+0L : 0x1.6a09e667f3bcc908p+0L : inexact
+= sqrt upward ldbl-96-m68k 0x2p+0L : 0x1.6a09e667f3bcc90ap+0L : inexact
+= sqrt downward ldbl-128 0x2p+0L : 0x1.6a09e667f3bcc908b2fb1366ea95p+0L : inexact
+= sqrt tonearest ldbl-128 0x2p+0L : 0x1.6a09e667f3bcc908b2fb1366ea95p+0L : inexact
+= sqrt towardzero ldbl-128 0x2p+0L : 0x1.6a09e667f3bcc908b2fb1366ea95p+0L : inexact
+= sqrt upward ldbl-128 0x2p+0L : 0x1.6a09e667f3bcc908b2fb1366ea96p+0L : inexact
+= sqrt downward ldbl-128ibm 0x2p+0L : 0x1.6a09e667f3bcc908b2fb1366ea8p+0L : inexact
+= sqrt tonearest ldbl-128ibm 0x2p+0L : 0x1.6a09e667f3bcc908b2fb1366ea8p+0L : inexact
+= sqrt towardzero ldbl-128ibm 0x2p+0L : 0x1.6a09e667f3bcc908b2fb1366ea8p+0L : inexact
+= sqrt upward ldbl-128ibm 0x2p+0L : 0x1.6a09e667f3bcc908b2fb1366ebp+0L : inexact
+sqrt 0.25
+= sqrt downward flt-32 0x4p-4f : 0x8p-4f :
+= sqrt tonearest flt-32 0x4p-4f : 0x8p-4f :
+= sqrt towardzero flt-32 0x4p-4f : 0x8p-4f :
+= sqrt upward flt-32 0x4p-4f : 0x8p-4f :
+= sqrt downward dbl-64 0x4p-4 : 0x8p-4 :
+= sqrt tonearest dbl-64 0x4p-4 : 0x8p-4 :
+= sqrt towardzero dbl-64 0x4p-4 : 0x8p-4 :
+= sqrt upward dbl-64 0x4p-4 : 0x8p-4 :
+= sqrt downward ldbl-96-intel 0x4p-4L : 0x8p-4L :
+= sqrt tonearest ldbl-96-intel 0x4p-4L : 0x8p-4L :
+= sqrt towardzero ldbl-96-intel 0x4p-4L : 0x8p-4L :
+= sqrt upward ldbl-96-intel 0x4p-4L : 0x8p-4L :
+= sqrt downward ldbl-96-m68k 0x4p-4L : 0x8p-4L :
+= sqrt tonearest ldbl-96-m68k 0x4p-4L : 0x8p-4L :
+= sqrt towardzero ldbl-96-m68k 0x4p-4L : 0x8p-4L :
+= sqrt upward ldbl-96-m68k 0x4p-4L : 0x8p-4L :
+= sqrt downward ldbl-128 0x4p-4L : 0x8p-4L :
+= sqrt tonearest ldbl-128 0x4p-4L : 0x8p-4L :
+= sqrt towardzero ldbl-128 0x4p-4L : 0x8p-4L :
+= sqrt upward ldbl-128 0x4p-4L : 0x8p-4L :
+= sqrt downward ldbl-128ibm 0x4p-4L : 0x8p-4L :
+= sqrt tonearest ldbl-128ibm 0x4p-4L : 0x8p-4L :
+= sqrt towardzero ldbl-128ibm 0x4p-4L : 0x8p-4L :
+= sqrt upward ldbl-128ibm 0x4p-4L : 0x8p-4L :
+sqrt 6642.25
+= sqrt downward flt-32 0x1.9f24p+12f : 0x5.18p+4f :
+= sqrt tonearest flt-32 0x1.9f24p+12f : 0x5.18p+4f :
+= sqrt towardzero flt-32 0x1.9f24p+12f : 0x5.18p+4f :
+= sqrt upward flt-32 0x1.9f24p+12f : 0x5.18p+4f :
+= sqrt downward dbl-64 0x1.9f24p+12 : 0x5.18p+4 :
+= sqrt tonearest dbl-64 0x1.9f24p+12 : 0x5.18p+4 :
+= sqrt towardzero dbl-64 0x1.9f24p+12 : 0x5.18p+4 :
+= sqrt upward dbl-64 0x1.9f24p+12 : 0x5.18p+4 :
+= sqrt downward ldbl-96-intel 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt tonearest ldbl-96-intel 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt towardzero ldbl-96-intel 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt upward ldbl-96-intel 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt downward ldbl-96-m68k 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt tonearest ldbl-96-m68k 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt towardzero ldbl-96-m68k 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt upward ldbl-96-m68k 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt downward ldbl-128 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt tonearest ldbl-128 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt towardzero ldbl-128 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt upward ldbl-128 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt downward ldbl-128ibm 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt tonearest ldbl-128ibm 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt towardzero ldbl-128ibm 0x1.9f24p+12L : 0x5.18p+4L :
+= sqrt upward ldbl-128ibm 0x1.9f24p+12L : 0x5.18p+4L :
+sqrt 15190.5625
+= sqrt downward flt-32 0x3.b569p+12f : 0x7.b4p+4f :
+= sqrt tonearest flt-32 0x3.b569p+12f : 0x7.b4p+4f :
+= sqrt towardzero flt-32 0x3.b569p+12f : 0x7.b4p+4f :
+= sqrt upward flt-32 0x3.b569p+12f : 0x7.b4p+4f :
+= sqrt downward dbl-64 0x3.b569p+12 : 0x7.b4p+4 :
+= sqrt tonearest dbl-64 0x3.b569p+12 : 0x7.b4p+4 :
+= sqrt towardzero dbl-64 0x3.b569p+12 : 0x7.b4p+4 :
+= sqrt upward dbl-64 0x3.b569p+12 : 0x7.b4p+4 :
+= sqrt downward ldbl-96-intel 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt tonearest ldbl-96-intel 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt towardzero ldbl-96-intel 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt upward ldbl-96-intel 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt downward ldbl-96-m68k 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt tonearest ldbl-96-m68k 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt towardzero ldbl-96-m68k 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt upward ldbl-96-m68k 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt downward ldbl-128 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt tonearest ldbl-128 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt towardzero ldbl-128 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt upward ldbl-128 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt downward ldbl-128ibm 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt tonearest ldbl-128ibm 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt towardzero ldbl-128ibm 0x3.b569p+12L : 0x7.b4p+4L :
+= sqrt upward ldbl-128ibm 0x3.b569p+12L : 0x7.b4p+4L :
+sqrt 0.75
+= sqrt downward flt-32 0xcp-4f : 0xd.db3d7p-4f : inexact
+= sqrt tonearest flt-32 0xcp-4f : 0xd.db3d7p-4f : inexact
+= sqrt towardzero flt-32 0xcp-4f : 0xd.db3d7p-4f : inexact
+= sqrt upward flt-32 0xcp-4f : 0xd.db3d8p-4f : inexact
+= sqrt downward dbl-64 0xcp-4 : 0xd.db3d742c2655p-4 : inexact
+= sqrt tonearest dbl-64 0xcp-4 : 0xd.db3d742c2655p-4 : inexact
+= sqrt towardzero dbl-64 0xcp-4 : 0xd.db3d742c2655p-4 : inexact
+= sqrt upward dbl-64 0xcp-4 : 0xd.db3d742c26558p-4 : inexact
+= sqrt downward ldbl-96-intel 0xcp-4L : 0xd.db3d742c265539dp-4L : inexact
+= sqrt tonearest ldbl-96-intel 0xcp-4L : 0xd.db3d742c265539ep-4L : inexact
+= sqrt towardzero ldbl-96-intel 0xcp-4L : 0xd.db3d742c265539dp-4L : inexact
+= sqrt upward ldbl-96-intel 0xcp-4L : 0xd.db3d742c265539ep-4L : inexact
+= sqrt downward ldbl-96-m68k 0xcp-4L : 0xd.db3d742c265539dp-4L : inexact
+= sqrt tonearest ldbl-96-m68k 0xcp-4L : 0xd.db3d742c265539ep-4L : inexact
+= sqrt towardzero ldbl-96-m68k 0xcp-4L : 0xd.db3d742c265539dp-4L : inexact
+= sqrt upward ldbl-96-m68k 0xcp-4L : 0xd.db3d742c265539ep-4L : inexact
+= sqrt downward ldbl-128 0xcp-4L : 0xd.db3d742c265539d92ba16b83c5cp-4L : inexact
+= sqrt tonearest ldbl-128 0xcp-4L : 0xd.db3d742c265539d92ba16b83c5cp-4L : inexact
+= sqrt towardzero ldbl-128 0xcp-4L : 0xd.db3d742c265539d92ba16b83c5cp-4L : inexact
+= sqrt upward ldbl-128 0xcp-4L : 0xd.db3d742c265539d92ba16b83c5c8p-4L : inexact
+= sqrt downward ldbl-128ibm 0xcp-4L : 0xd.db3d742c265539d92ba16b83c4p-4L : inexact
+= sqrt tonearest ldbl-128ibm 0xcp-4L : 0xd.db3d742c265539d92ba16b83c4p-4L : inexact
+= sqrt towardzero ldbl-128ibm 0xcp-4L : 0xd.db3d742c265539d92ba16b83c4p-4L : inexact
+= sqrt upward ldbl-128ibm 0xcp-4L : 0xd.db3d742c265539d92ba16b83c8p-4L : inexact
diff --git a/math/gen-auto-libm-tests.c b/math/gen-auto-libm-tests.c
new file mode 100644
index 0000000000..57053ad420
--- /dev/null
+++ b/math/gen-auto-libm-tests.c
@@ -0,0 +1,1626 @@
+/* Generate expected output for libm tests with MPFR and MPC.
+   Copyright (C) 2013 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+/* Compile this program as:
+
+   gcc -std=gnu99 -O2 -Wall -Wextra gen-auto-libm-tests.c -lmpc -lmpfr -lgmp \
+     -o gen-auto-libm-tests
+
+   (use of current MPC and MPFR versions recommended) and run it as:
+
+   gen-auto-libm-tests auto-libm-test-in auto-libm-test-out
+
+   The input file auto-libm-test-in contains three kinds of lines:
+
+   Lines beginning with "#" are comments, and are ignored, as are
+   empty lines.
+
+   Other lines are test lines, of the form "function input1 input2
+   ... [flag1 flag2 ...]".  Inputs are either finite real numbers or
+   integers, depending on the function under test.  Real numbers may
+   be in any form acceptable to mpfr_strtofr (base 0); integers in any
+   form acceptable to mpz_set_str (base 0).  In addition, real numbers
+   may be certain special strings such as "pi", as listed in the
+   special_real_inputs array.
+
+   Each flag is a flag name possibly followed by a series of
+   ":condition".  Conditions may be any of the names of floating-point
+   formats in the floating_point_formats array, "long32" and "long64"
+   to indicate the number of bits in the "long" type, or other strings
+   for which libm-test.inc defines a TEST_COND_<condition> macro (with
+   "-"- changed to "_" in the condition name) evaluating to nonzero
+   when the condition is true and zero when the condition is false.
+   The meaning is that the flag applies to the test if all the listed
+   conditions are true.  "flag:cond1:cond2 flag:cond3:cond4" means the
+   flag applies if ((cond1 && cond2) || (cond3 && cond4)).
+
+   A real number specified as an input is considered to represent the
+   set of real numbers arising from rounding the given number in any
+   direction for any supported floating-point format; any roundings
+   that give infinity are ignored.  Each input on a test line has all
+   the possible roundings considered independently.  Each resulting
+   choice of the tuple of inputs to the function is ignored if the
+   mathematical result of the function involves a NaN or an exact
+   infinity, and is otherwise considered for each floating-point
+   format for which all those inputs are exactly representable.  Thus
+   tests may result in "overflow", "underflow" and "inexact"
+   exceptions; "invalid" may arise only when the final result type is
+   an integer type and it is the conversion of a mathematically
+   defined finite result to integer type that results in that
+   exception.
+
+   By default, it is assumed that "overflow" and "underflow"
+   exceptions should be correct, but that "inexact" exceptions should
+   only be correct for functions listed as exactly determined.  For
+   such functions, "underflow" exceptions should respect whether the
+   machine has before-rounding or after-rounding tininess detection.
+   For other functions, it is considered that if the exact result is
+   somewhere between the greatest magnitude subnormal of a given sign
+   (exclusive) and the least magnitude normal of that sign
+   (inclusive), underflow exceptions are permitted but optional on all
+   machines, and they are also permitted but optional for smaller
+   subnormal exact results for functions that are not exactly
+   determined.  errno setting is expected for overflow to infinity and
+   underflow to zero (for real functions), and for out-of-range
+   conversion of a finite result to integer type, and is considered
+   permitted but optional for all other cases where overflow
+   exceptions occur, and where underflow exceptions occur or are
+   permitted.  In other cases (where no overflow or underflow is
+   permitted), errno is expected to be left unchanged.
+
+   The flag "no-test-inline" indicates a test is disabled for inline
+   function testing; "xfail" indicates the test is disabled as
+   expected to produce incorrect results.  Otherwise, test flags are
+   of the form "spurious-<exception>" and "missing-<exception>", for
+   any exception ("overflow", "underflow", "inexact", "invalid",
+   "divbyzero"), "spurious-errno" and "missing-errno", to indicate
+   when tests are expected to deviate from the exception and errno
+   settings corresponding to the mathematical results.  "xfail",
+   "spurious-" and "missing-" flags should be accompanied by a comment
+   referring to an open bug in glibc Bugzilla.
+
+   The output file auto-libm-test-out contains the test lines from
+   auto-libm-test-in, and, after the line for a given test, some
+   number of output test lines.  An output test line is of the form "=
+   function rounding-mode format input1 input2 ... : output1 output2
+   ... : flags".  rounding-mode is "tonearest", "towardzero", "upward"
+   or "downward".  format is a name from the floating_point_formats
+   array, possibly followed by a sequence of ":flag" for flags from
+   "long32", "long64", "before-rounding" and "after-rounding" (the
+   last two indicating tests where expectations for underflow
+   exceptions depend on how the architecture detects tininess).
+   Inputs and outputs are specified as hex floats with the required
+   suffix for the floating-point type, or plus_infty or minus_infty
+   for infinite expected results, or as integer constant expressions
+   (not necessarily with the right type) or IGNORE for integer inputs
+   and outputs.  Flags are "no-test-inline", "xfail", "<exception>",
+   "<exception>-ok", "errno-<value>", "errno-<value>-ok", where
+   "<exception>" and "errno-<value>" are unconditional, indicating
+   that a correct result means the given exception should be raised or
+   errno should be set to the given value, and other settings may be
+   conditional or unconditional; "-ok" means not to test for the given
+   exception or errno value (whether because it was marked as possibly
+   missing or spurious, or because the calculation of correct results
+   indicated it was optional).  */
+
+#define _GNU_SOURCE
+
+#include <assert.h>
+#include <ctype.h>
+#include <errno.h>
+#include <error.h>
+#include <stdbool.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <gmp.h>
+#include <mpfr.h>
+#include <mpc.h>
+
+#define ARRAY_SIZE(A) (sizeof (A) / sizeof ((A)[0]))
+
+/* The supported floating-point formats.  */
+typedef enum
+  {
+    fp_flt_32,
+    fp_dbl_64,
+    fp_ldbl_96_intel,
+    fp_ldbl_96_m68k,
+    fp_ldbl_128,
+    fp_ldbl_128ibm,
+    fp_num_formats,
+    fp_first_format = 0
+  } fp_format;
+
+/* Structure describing a single floating-point format.  */
+typedef struct
+{
+  /* The name of the format.  */
+  const char *name;
+  /* The suffix to use on floating-point constants with this
+     format.  */
+  const char *suffix;
+  /* A string for the largest normal value, or NULL for IEEE formats
+     where this can be determined automatically.  */
+  const char *max_string;
+  /* The number of mantissa bits.  */
+  int mant_dig;
+  /* The least N such that 2^N overflows.  */
+  int max_exp;
+  /* One more than the least N such that 2^N is normal.  */
+  int min_exp;
+  /* The largest normal value.  */
+  mpfr_t max;
+  /* The least positive normal value, 2^(MIN_EXP-1).  */
+  mpfr_t min;
+  /* The greatest positive subnormal value.  */
+  mpfr_t subnorm_max;
+  /* The least positive subnormal value, 2^(MIN_EXP-MANT_DIG).  */
+  mpfr_t subnorm_min;
+} fp_format_desc;
+
+/* List of floating-point formats, in the same order as the fp_format
+   enumeration.  */
+static fp_format_desc fp_formats[fp_num_formats] =
+  {
+    { "flt-32", "f", NULL, 24, 128, -125, {}, {}, {}, {} },
+    { "dbl-64", "", NULL, 53, 1024, -1021, {}, {}, {}, {} },
+    { "ldbl-96-intel", "L", NULL, 64, 16384, -16381, {}, {}, {}, {} },
+    { "ldbl-96-m68k", "L", NULL, 64, 16384, -16382, {}, {}, {}, {} },
+    { "ldbl-128", "L", NULL, 113, 16384, -16381, {}, {}, {}, {} },
+    { "ldbl-128ibm", "L", "0x1.fffffffffffff7ffffffffffff8p+1023",
+      106, 1024, -968, {}, {}, {}, {} },
+  };
+
+/* The supported rounding modes.  */
+typedef enum
+  {
+    rm_downward,
+    rm_tonearest,
+    rm_towardzero,
+    rm_upward,
+    rm_num_modes,
+    rm_first_mode = 0
+  } rounding_mode;
+
+/* Structure describing a single rounding mode.  */
+typedef struct
+{
+  /* The name of the rounding mode.  */
+  const char *name;
+  /* The MPFR rounding mode.  */
+  mpfr_rnd_t mpfr_mode;
+} rounding_mode_desc;
+
+/* List of rounding modes, in the same order as the rounding_mode
+   enumeration.  */
+static const rounding_mode_desc rounding_modes[rm_num_modes] =
+  {
+    { "downward", MPFR_RNDD },
+    { "tonearest", MPFR_RNDN },
+    { "towardzero", MPFR_RNDZ },
+    { "upward", MPFR_RNDU },
+  };
+
+/* The supported exceptions.  */
+typedef enum
+  {
+    exc_divbyzero,
+    exc_inexact,
+    exc_invalid,
+    exc_overflow,
+    exc_underflow,
+    exc_num_exceptions,
+    exc_first_exception = 0
+  } fp_exception;
+
+/* List of exceptions, in the same order as the fp_exception
+   enumeration.  */
+static const char *const exceptions[exc_num_exceptions] =
+  {
+    "divbyzero",
+    "inexact",
+    "invalid",
+    "overflow",
+    "underflow",
+  };
+
+/* The internal precision to use for most MPFR calculations, which
+   must be at least 2 more than the greatest precision of any
+   supported floating-point format.  */
+static int internal_precision;
+
+/* A value that overflows all supported floating-point formats.  */
+static mpfr_t global_max;
+
+/* A value that is at most half the least subnormal in any
+   floating-point format and so is rounded the same way as all
+   sufficiently small positive values.  */
+static mpfr_t global_min;
+
+/* The maximum number of (real or integer) arguments to a function
+   handled by this program (complex arguments count as two real
+   arguments).  */
+#define MAX_NARGS 4
+
+/* The maximum number of (real or integer) return values from a
+   function handled by this program.  */
+#define MAX_NRET 2
+
+/* A type of a function argument or return value.  */
+typedef enum
+  {
+    /* No type (not a valid argument or return value).  */
+    type_none,
+    /* A floating-point value with the type corresponding to that of
+       the function.  */
+    type_fp,
+    /* An integer value of type int.  */
+    type_int,
+    /* An integer value of type long.  */
+    type_long,
+    /* An integer value of type long long.  */
+    type_long_long,
+  } arg_ret_type;
+
+/* A type of a generic real or integer value.  */
+typedef enum
+  {
+    /* No type.  */
+    gtype_none,
+    /* Floating-point (represented with MPFR).  */
+    gtype_fp,
+    /* Integer (represented with GMP).  */
+    gtype_int,
+  } generic_value_type;
+
+/* A generic value (argument or result).  */
+typedef struct
+{
+  /* The type of this value.  */
+  generic_value_type type;
+  /* Its value.  */
+  union
+  {
+    mpfr_t f;
+    mpz_t i;
+  } value;
+} generic_value;
+
+/* A type of input flag.  */
+typedef enum
+  {
+    flag_no_test_inline,
+    flag_xfail,
+    /* The "spurious" and "missing" flags must be in the same order as
+       the fp_exception enumeration.  */
+    flag_spurious_divbyzero,
+    flag_spurious_inexact,
+    flag_spurious_invalid,
+    flag_spurious_overflow,
+    flag_spurious_underflow,
+    flag_spurious_errno,
+    flag_missing_divbyzero,
+    flag_missing_inexact,
+    flag_missing_invalid,
+    flag_missing_overflow,
+    flag_missing_underflow,
+    flag_missing_errno,
+    num_input_flag_types,
+    flag_first_flag = 0,
+    flag_spurious_first = flag_spurious_divbyzero,
+    flag_missing_first = flag_missing_divbyzero
+  } input_flag_type;
+
+/* List of flags, in the same order as the input_flag_type
+   enumeration.  */
+static const char *const input_flags[num_input_flag_types] =
+  {
+    "no-test-inline",
+    "xfail",
+    "spurious-divbyzero",
+    "spurious-inexact",
+    "spurious-invalid",
+    "spurious-overflow",
+    "spurious-underflow",
+    "spurious-errno",
+    "missing-divbyzero",
+    "missing-inexact",
+    "missing-invalid",
+    "missing-overflow",
+    "missing-underflow",
+    "missing-errno",
+  };
+
+/* An input flag, possibly conditional.  */
+typedef struct
+{
+  /* The type of this flag.  */
+  input_flag_type type;
+  /* The conditions on this flag, as a string ":cond1:cond2..." or
+     NULL.  */
+  const char *cond;
+} input_flag;
+
+/* Structure describing a single test from the input file (which may
+   expand into many tests in the output).  The choice of function,
+   which implies the numbers and types of arguments and results, is
+   implicit rather than stored in this structure (except as part of
+   the source line).  */
+typedef struct
+{
+  /* The text of the input line describing the test, including the
+     trailing newline.  */
+  const char *line;
+  /* The number of combinations of interpretations of input values for
+     different floating-point formats and rounding modes.  */
+  size_t num_input_cases;
+  /* The corresponding lists of inputs.  */
+  generic_value **inputs;
+  /* The number of flags for this test.  */
+  size_t num_flags;
+  /* The corresponding list of flags.  */
+  input_flag *flags;
+  /* The old output for this test.  */
+  const char *old_output;
+} input_test;
+
+/* Ways to calculate a function.  */
+typedef enum
+  {
+    /* MPFR function with a single argument and result.  */
+    mpfr_f_f,
+  } func_calc_method;
+
+/* Description of how to calculate a function.  */
+typedef struct
+{
+  /* Which method is used to calculate the function.  */
+  func_calc_method method;
+  /* The specific function called.  */
+  union
+  {
+    int (*mpfr_f_f) (mpfr_t, const mpfr_t, mpfr_rnd_t);
+  } func;
+} func_calc_desc;
+
+/* Structure describing a function handled by this program.  */
+typedef struct
+{
+  /* The name of the function.  */
+  const char *name;
+  /* The number of arguments.  */
+  size_t num_args;
+  /* The types of the arguments.  */
+  arg_ret_type arg_types[MAX_NARGS];
+  /* The number of return values.  */
+  size_t num_ret;
+  /* The types of the return values.  */
+  arg_ret_type ret_types[MAX_NRET];
+  /* Whether the function has exactly determined results and
+     exceptions.  */
+  bool exact;
+  /* Whether the function is a complex function, so errno setting is
+     optional.  */
+  bool complex_fn;
+  /* How to calculate this function.  */
+  func_calc_desc calc;
+  /* The number of tests allocated for this function.  */
+  size_t num_tests_alloc;
+  /* The number of tests for this function.  */
+  size_t num_tests;
+  /* The tests themselves.  */
+  input_test *tests;
+} test_function;
+
+#define FUNC_mpfr_f_f(NAME, MPFR_FUNC, EXACT)		\
+  {							\
+    NAME, 1, { type_fp }, 1, { type_fp }, EXACT, false,	\
+    { mpfr_f_f, { .mpfr_f_f = MPFR_FUNC } }, 0, 0, NULL \
+  }
+
+/* List of functions handled by this program.  */
+static test_function test_functions[] =
+  {
+    FUNC_mpfr_f_f ("sqrt", mpfr_sqrt, true),
+  };
+
+/* Allocate memory, with error checking.  */
+
+static void *
+xmalloc (size_t n)
+{
+  void *p = malloc (n);
+  if (p == NULL)
+    error (EXIT_FAILURE, errno, "xmalloc failed");
+  return p;
+}
+
+static void *
+xrealloc (void *p, size_t n)
+{
+  p = realloc (p, n);
+  if (p == NULL)
+    error (EXIT_FAILURE, errno, "xrealloc failed");
+  return p;
+}
+
+static char *
+xstrdup (const char *s)
+{
+  char *p = strdup (s);
+  if (p == NULL)
+    error (EXIT_FAILURE, errno, "xstrdup failed");
+  return p;
+}
+
+/* Assert that the result of an MPFR operation was exact; that is,
+   that the returned ternary value was 0.  */
+
+static void
+assert_exact (int i)
+{
+  assert (i == 0);
+}
+
+/* Return the generic type of an argument or return value type T.  */
+
+static generic_value_type
+generic_arg_ret_type (arg_ret_type t)
+{
+  switch (t)
+    {
+    case type_fp:
+      return gtype_fp;
+
+    case type_int:
+    case type_long:
+    case type_long_long:
+      return gtype_int;
+
+    default:
+      abort ();
+    }
+}
+
+/* Free a generic_value *V.  */
+
+static void
+generic_value_free (generic_value *v)
+{
+  switch (v->type)
+    {
+    case gtype_fp:
+      mpfr_clear (v->value.f);
+      break;
+
+    case gtype_int:
+      mpz_clear (v->value.i);
+      break;
+
+    default:
+      abort ();
+    }
+}
+
+/* Copy a generic_value *SRC to *DEST.  */
+
+static void
+generic_value_copy (generic_value *dest, const generic_value *src)
+{
+  dest->type = src->type;
+  switch (src->type)
+    {
+    case gtype_fp:
+      mpfr_init (dest->value.f);
+      assert_exact (mpfr_set (dest->value.f, src->value.f, MPFR_RNDN));
+      break;
+
+    case gtype_int:
+      mpz_init (dest->value.i);
+      mpz_set (dest->value.i, src->value.i);
+      break;
+
+    default:
+      abort ();
+    }
+}
+
+/* Initialize data for floating-point formats.  */
+
+static void
+init_fp_formats ()
+{
+  int global_max_exp = 0, global_min_subnorm_exp = 0;
+  for (fp_format f = fp_first_format; f < fp_num_formats; f++)
+    {
+      if (fp_formats[f].mant_dig + 2 > internal_precision)
+	internal_precision = fp_formats[f].mant_dig + 2;
+      if (fp_formats[f].max_exp > global_max_exp)
+	global_max_exp = fp_formats[f].max_exp;
+      int min_subnorm_exp = fp_formats[f].min_exp - fp_formats[f].mant_dig;
+      if (min_subnorm_exp < global_min_subnorm_exp)
+	global_min_subnorm_exp = min_subnorm_exp;
+      mpfr_init2 (fp_formats[f].max, fp_formats[f].mant_dig);
+      if (fp_formats[f].max_string != NULL)
+	{
+	  char *ep = NULL;
+	  assert_exact (mpfr_strtofr (fp_formats[f].max,
+				      fp_formats[f].max_string,
+				      &ep, 0, MPFR_RNDN));
+	  assert (*ep == 0);
+	}
+      else
+	{
+	  assert_exact (mpfr_set_ui_2exp (fp_formats[f].max, 1,
+					  fp_formats[f].max_exp,
+					  MPFR_RNDN));
+	  mpfr_nextbelow (fp_formats[f].max);
+	}
+      mpfr_init2 (fp_formats[f].min, fp_formats[f].mant_dig);
+      assert_exact (mpfr_set_ui_2exp (fp_formats[f].min, 1,
+				      fp_formats[f].min_exp - 1,
+				      MPFR_RNDN));
+      mpfr_init2 (fp_formats[f].subnorm_max, fp_formats[f].mant_dig);
+      assert_exact (mpfr_set (fp_formats[f].subnorm_max, fp_formats[f].min,
+			      MPFR_RNDN));
+      mpfr_nextbelow (fp_formats[f].subnorm_max);
+      mpfr_nextbelow (fp_formats[f].subnorm_max);
+      mpfr_init2 (fp_formats[f].subnorm_min, fp_formats[f].mant_dig);
+      assert_exact (mpfr_set_ui_2exp (fp_formats[f].subnorm_min, 1,
+				      min_subnorm_exp, MPFR_RNDN));
+    }
+  mpfr_set_default_prec (internal_precision);
+  mpfr_init (global_max);
+  assert_exact (mpfr_set_ui_2exp (global_max, 1, global_max_exp, MPFR_RNDN));
+  mpfr_init (global_min);
+  assert_exact (mpfr_set_ui_2exp (global_min, 1, global_min_subnorm_exp - 1,
+				  MPFR_RNDN));
+}
+
+/* Fill in mpfr_t values for special strings in input arguments.  */
+
+static size_t
+special_fill_max (mpfr_t res0, mpfr_t res1 __attribute__ ((unused)),
+		  fp_format format)
+{
+  mpfr_init2 (res0, fp_formats[format].mant_dig);
+  assert_exact (mpfr_set (res0, fp_formats[format].max, MPFR_RNDN));
+  return 1;
+}
+
+static size_t
+special_fill_minus_max (mpfr_t res0, mpfr_t res1 __attribute__ ((unused)),
+			fp_format format)
+{
+  mpfr_init2 (res0, fp_formats[format].mant_dig);
+  assert_exact (mpfr_neg (res0, fp_formats[format].max, MPFR_RNDN));
+  return 1;
+}
+
+static size_t
+special_fill_pi (mpfr_t res0, mpfr_t res1, fp_format format)
+{
+  mpfr_init2 (res0, fp_formats[format].mant_dig);
+  mpfr_const_pi (res0, MPFR_RNDU);
+  mpfr_init2 (res1, fp_formats[format].mant_dig);
+  mpfr_const_pi (res1, MPFR_RNDD);
+  return 2;
+}
+
+static size_t
+special_fill_minus_pi (mpfr_t res0, mpfr_t res1, fp_format format)
+{
+  mpfr_init2 (res0, fp_formats[format].mant_dig);
+  mpfr_const_pi (res0, MPFR_RNDU);
+  assert_exact (mpfr_neg (res0, res0, MPFR_RNDN));
+  mpfr_init2 (res1, fp_formats[format].mant_dig);
+  mpfr_const_pi (res1, MPFR_RNDD);
+  assert_exact (mpfr_neg (res1, res1, MPFR_RNDN));
+  return 2;
+}
+
+/* A special string accepted in input arguments.  */
+typedef struct
+{
+  /* The string.  */
+  const char *str;
+  /* The function that interprets it for a given floating-point
+     format, filling in up to two mpfr_t values and returning the
+     number of values filled.  */
+  size_t (*func) (mpfr_t, mpfr_t, fp_format);
+} special_real_input;
+
+/* List of special strings accepted in input arguments.  */
+
+static const special_real_input special_real_inputs[] =
+  {
+    { "max", special_fill_max },
+    { "-max", special_fill_minus_max },
+    { "pi", special_fill_pi },
+    { "-pi", special_fill_minus_pi },
+  };
+
+/* Given a real number R computed in round-to-zero mode, set the
+   lowest bit as a sticky bit if INEXACT, and saturate the exponent
+   range for very large or small values.  */
+
+static void
+adjust_real (mpfr_t r, bool inexact)
+{
+  if (!inexact)
+    return;
+  /* NaNs are exact, as are infinities in round-to-zero mode.  */
+  assert (mpfr_regular_p (r));
+  if (mpfr_cmpabs (r, global_min) < 0)
+    assert_exact (mpfr_copysign (r, global_min, r, MPFR_RNDN));
+  else if (mpfr_cmpabs (r, global_max) > 0)
+    assert_exact (mpfr_copysign (r, global_max, r, MPFR_RNDN));
+  else
+    {
+      mpz_t tmp;
+      mpz_init (tmp);
+      mpfr_exp_t e = mpfr_get_z_2exp (tmp, r);
+      mpz_setbit (tmp, 0);
+      assert_exact (mpfr_set_z_2exp (r, tmp, e, MPFR_RNDN));
+      mpz_clear (tmp);
+    }
+}
+
+/* Given a finite real number R with sticky bit, compute the roundings
+   to FORMAT in each rounding mode, storing the results in RES, the
+   before-rounding exceptions in EXC_BEFORE and the after-rounding
+   exceptions in EXC_AFTER.  */
+
+static void
+round_real (mpfr_t res[rm_num_modes],
+	    unsigned int exc_before[rm_num_modes],
+	    unsigned int exc_after[rm_num_modes],
+	    mpfr_t r, fp_format format)
+{
+  assert (mpfr_number_p (r));
+  for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
+    {
+      mpfr_init2 (res[m], fp_formats[format].mant_dig);
+      exc_before[m] = exc_after[m] = 0;
+      bool inexact = mpfr_set (res[m], r, rounding_modes[m].mpfr_mode);
+      if (mpfr_cmpabs (res[m], fp_formats[format].max) > 0)
+	{
+	  inexact = true;
+	  exc_before[m] |= 1U << exc_overflow;
+	  exc_after[m] |= 1U << exc_overflow;
+	  bool overflow_inf;
+	  switch (m)
+	    {
+	    case rm_tonearest:
+	      overflow_inf = true;
+	      break;
+	    case rm_towardzero:
+	      overflow_inf = false;
+	      break;
+	    case rm_downward:
+	      overflow_inf = mpfr_signbit (res[m]);
+	      break;
+	    case rm_upward:
+	      overflow_inf = !mpfr_signbit (res[m]);
+	      break;
+	    default:
+	      abort ();
+	    }
+	  if (overflow_inf)
+	    mpfr_set_inf (res[m], mpfr_signbit (res[m]) ? -1 : 1);
+	  else
+	    assert_exact (mpfr_copysign (res[m], fp_formats[format].max,
+					 res[m], MPFR_RNDN));
+	}
+      if (mpfr_cmpabs (r, fp_formats[format].min) < 0)
+	{
+	  /* Tiny before rounding; may or may not be tiny after
+	     rounding, and underflow applies only if also inexact
+	     around rounding to a possibly subnormal value.  */
+	  bool tiny_after_rounding
+	    = mpfr_cmpabs (res[m], fp_formats[format].min) < 0;
+	  /* To round to a possibly subnormal value, and determine
+	     inexactness as a subnormal in the process, scale up and
+	     round to integer, then scale back down.  */
+	  mpfr_t tmp;
+	  mpfr_init (tmp);
+	  assert_exact (mpfr_mul_2si (tmp, r, (fp_formats[format].mant_dig
+					       - fp_formats[format].min_exp),
+				      MPFR_RNDN));
+	  int rint_res = mpfr_rint (tmp, tmp, rounding_modes[m].mpfr_mode);
+	  /* The integer must be representable.  */
+	  assert (rint_res == 0 || rint_res == 2 || rint_res == -2);
+	  /* If rounding to full precision was inexact, so must
+	     rounding to subnormal precision be inexact.  */
+	  if (inexact)
+	    assert (rint_res != 0);
+	  else
+	    inexact = rint_res != 0;
+	  assert_exact (mpfr_mul_2si (res[m], tmp,
+				      (fp_formats[format].min_exp
+				       - fp_formats[format].mant_dig),
+				      MPFR_RNDN));
+	  mpfr_clear (tmp);
+	  if (inexact)
+	    {
+	      exc_before[m] |= 1U << exc_underflow;
+	      if (tiny_after_rounding)
+		exc_after[m] |= 1U << exc_underflow;
+	    }
+	}
+      if (inexact)
+	{
+	  exc_before[m] |= 1U << exc_inexact;
+	  exc_after[m] |= 1U << exc_inexact;
+	}
+    }
+}
+
+/* Handle the input argument at ARG (NUL-terminated), updating the
+   lists of test inputs in IT accordingly.  NUM_PREV_ARGS arguments
+   are already in those lists.  The argument, of type GTYPE, comes
+   from file FILENAME, line LINENO.  */
+
+static void
+handle_input_arg (const char *arg, input_test *it, size_t num_prev_args,
+		  generic_value_type gtype,
+		  const char *filename, unsigned int lineno)
+{
+  size_t num_values = 0;
+  generic_value values[2 * fp_num_formats];
+  switch (gtype)
+    {
+    case gtype_fp:
+      for (fp_format f = fp_first_format; f < fp_num_formats; f++)
+	{
+	  mpfr_t extra_values[2];
+	  size_t num_extra_values = 0;
+	  for (size_t i = 0; i < ARRAY_SIZE (special_real_inputs); i++)
+	    {
+	      if (strcmp (arg, special_real_inputs[i].str) == 0)
+		{
+		  num_extra_values
+		    = special_real_inputs[i].func (extra_values[0],
+						   extra_values[1], f);
+		  assert (num_extra_values > 0
+			  && num_extra_values <= ARRAY_SIZE (extra_values));
+		  break;
+		}
+	    }
+	  if (num_extra_values == 0)
+	    {
+	      mpfr_t tmp;
+	      char *ep;
+	      mpfr_init (tmp);
+	      bool inexact = mpfr_strtofr (tmp, arg, &ep, 0, MPFR_RNDZ);
+	      if (*ep != 0 || !mpfr_number_p (tmp))
+		error_at_line (EXIT_FAILURE, 0, filename, lineno,
+			       "bad floating-point argument: '%s'", arg);
+	      adjust_real (tmp, inexact);
+	      mpfr_t rounded[rm_num_modes];
+	      unsigned int exc_before[rm_num_modes];
+	      unsigned int exc_after[rm_num_modes];
+	      round_real (rounded, exc_before, exc_after, tmp, f);
+	      mpfr_clear (tmp);
+	      if (mpfr_number_p (rounded[rm_upward]))
+		{
+		  mpfr_init2 (extra_values[num_extra_values],
+			      fp_formats[f].mant_dig);
+		  assert_exact (mpfr_set (extra_values[num_extra_values],
+					  rounded[rm_upward], MPFR_RNDN));
+		  num_extra_values++;
+		}
+	      if (mpfr_number_p (rounded[rm_downward]))
+		{
+		  mpfr_init2 (extra_values[num_extra_values],
+			      fp_formats[f].mant_dig);
+		  assert_exact (mpfr_set (extra_values[num_extra_values],
+					  rounded[rm_downward], MPFR_RNDN));
+		  num_extra_values++;
+		}
+	      for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
+		mpfr_clear (rounded[m]);
+	    }
+	  for (size_t i = 0; i < num_extra_values; i++)
+	    {
+	      bool found = false;
+	      for (size_t j = 0; j < num_values; j++)
+		{
+		  if (mpfr_equal_p (values[j].value.f, extra_values[i])
+		      && ((mpfr_signbit (values[j].value.f) != 0)
+			  == (mpfr_signbit (extra_values[i]) != 0)))
+		    {
+		      found = true;
+		      break;
+		    }
+		}
+	      if (!found)
+		{
+		  assert (num_values < ARRAY_SIZE (values));
+		  values[num_values].type = gtype_fp;
+		  mpfr_init2 (values[num_values].value.f,
+			      fp_formats[f].mant_dig);
+		  assert_exact (mpfr_set (values[num_values].value.f,
+					  extra_values[i], MPFR_RNDN));
+		  num_values++;
+		}
+	      mpfr_clear (extra_values[i]);
+	    }
+	}
+      break;
+
+    case gtype_int:
+      num_values = 1;
+      values[0].type = gtype_int;
+      int ret = mpz_init_set_str (values[0].value.i, arg, 0);
+      if (ret != 0)
+	error_at_line (EXIT_FAILURE, 0, filename, lineno,
+		       "bad integer argument: '%s'", arg);
+      break;
+
+    default:
+      abort ();
+    }
+  assert (num_values > 0 && num_values <= ARRAY_SIZE (values));
+  if (it->num_input_cases >= SIZE_MAX / num_values)
+    error_at_line (EXIT_FAILURE, 0, filename, lineno, "too many input cases");
+  generic_value **old_inputs = it->inputs;
+  size_t new_num_input_cases = it->num_input_cases * num_values;
+  generic_value **new_inputs = xmalloc (new_num_input_cases
+					* sizeof (new_inputs[0]));
+  for (size_t i = 0; i < it->num_input_cases; i++)
+    {
+      for (size_t j = 0; j < num_values; j++)
+	{
+	  size_t idx = i * num_values + j;
+	  new_inputs[idx] = xmalloc ((num_prev_args + 1)
+				     * sizeof (new_inputs[idx][0]));
+	  for (size_t k = 0; k < num_prev_args; k++)
+	    generic_value_copy (&new_inputs[idx][k], &old_inputs[i][k]);
+	  generic_value_copy (&new_inputs[idx][num_prev_args], &values[j]);
+	}
+      for (size_t j = 0; j < num_prev_args; j++)
+	generic_value_free (&old_inputs[i][j]);
+      free (old_inputs[i]);
+    }
+  free (old_inputs);
+  for (size_t i = 0; i < num_values; i++)
+    generic_value_free (&values[i]);
+  it->inputs = new_inputs;
+  it->num_input_cases = new_num_input_cases;
+}
+
+/* Handle the input flag ARG (NUL-terminated), storing it in *FLAG.
+   The flag comes from file FILENAME, line LINENO.  */
+
+static void
+handle_input_flag (char *arg, input_flag *flag,
+		   const char *filename, unsigned int lineno)
+{
+  char *ep = strchr (arg, ':');
+  if (ep == NULL)
+    {
+      ep = strchr (arg, 0);
+      assert (ep != NULL);
+    }
+  char c = *ep;
+  *ep = 0;
+  bool found = false;
+  for (input_flag_type i = flag_first_flag; i <= num_input_flag_types; i++)
+    {
+      if (strcmp (arg, input_flags[i]) == 0)
+	{
+	  found = true;
+	  flag->type = i;
+	  break;
+	}
+    }
+  if (!found)
+    error_at_line (EXIT_FAILURE, 0, filename, lineno, "unknown flag: '%s'",
+		   arg);
+  *ep = c;
+  if (c == 0)
+    flag->cond = NULL;
+  else
+    flag->cond = xstrdup (ep);
+}
+
+/* Add the test LINE (file FILENAME, line LINENO) to the test
+   data.  */
+
+static void
+add_test (char *line, const char *filename, unsigned int lineno)
+{
+  size_t num_tokens = 1;
+  char *p = line;
+  while ((p = strchr (p, ' ')) != NULL)
+    {
+      num_tokens++;
+      p++;
+    }
+  if (num_tokens < 2)
+    error_at_line (EXIT_FAILURE, 0, filename, lineno,
+		   "line too short: '%s'", line);
+  p = strchr (line, ' ');
+  size_t func_name_len = p - line;
+  for (size_t i = 0; i < ARRAY_SIZE (test_functions); i++)
+    {
+      if (func_name_len == strlen (test_functions[i].name)
+	  && strncmp (line, test_functions[i].name, func_name_len) == 0)
+	{
+	  test_function *tf = &test_functions[i];
+	  if (num_tokens < 1 + tf->num_args)
+	    error_at_line (EXIT_FAILURE, 0, filename, lineno,
+			   "line too short: '%s'", line);
+	  if (tf->num_tests == tf->num_tests_alloc)
+	    {
+	      tf->num_tests_alloc = 2 * tf->num_tests_alloc + 16;
+	      tf->tests
+		= xrealloc (tf->tests,
+			    tf->num_tests_alloc * sizeof (tf->tests[0]));
+	    }
+	  input_test *it = &tf->tests[tf->num_tests];
+	  it->line = line;
+	  it->num_input_cases = 1;
+	  it->inputs = xmalloc (sizeof (it->inputs[0]));
+	  it->inputs[0] = NULL;
+	  it->old_output = NULL;
+	  p++;
+	  for (size_t j = 0; j < tf->num_args; j++)
+	    {
+	      char *ep = strchr (p, ' ');
+	      if (ep == NULL)
+		{
+		  ep = strchr (p, '\n');
+		  assert (ep != NULL);
+		}
+	      if (ep == p)
+		error_at_line (EXIT_FAILURE, 0, filename, lineno,
+			       "empty token in line: '%s'", line);
+	      for (char *t = p; t < ep; t++)
+		if (isspace ((unsigned char) *t))
+		  error_at_line (EXIT_FAILURE, 0, filename, lineno,
+				 "whitespace in token in line: '%s'", line);
+	      char c = *ep;
+	      *ep = 0;
+	      handle_input_arg (p, it, j,
+				generic_arg_ret_type (tf->arg_types[j]),
+				filename, lineno);
+	      *ep = c;
+	      p = ep + 1;
+	    }
+	  it->num_flags = num_tokens - 1 - tf->num_args;
+	  it->flags = xmalloc (it->num_flags * sizeof (it->flags[0]));
+	  for (size_t j = 0; j < it->num_flags; j++)
+	    {
+	      char *ep = strchr (p, ' ');
+	      if (ep == NULL)
+		{
+		  ep = strchr (p, '\n');
+		  assert (ep != NULL);
+		}
+	      if (ep == p)
+		error_at_line (EXIT_FAILURE, 0, filename, lineno,
+			       "empty token in line: '%s'", line);
+	      for (char *t = p; t < ep; t++)
+		if (isspace ((unsigned char) *t))
+		  error_at_line (EXIT_FAILURE, 0, filename, lineno,
+				 "whitespace in token in line: '%s'", line);
+	      char c = *ep;
+	      *ep = 0;
+	      handle_input_flag (p, &it->flags[j], filename, lineno);
+	      *ep = c;
+	      p = ep + 1;
+	    }
+	  assert (*p == 0);
+	  tf->num_tests++;
+	  return;
+	}
+    }
+  error_at_line (EXIT_FAILURE, 0, filename, lineno,
+		 "unknown function in line: '%s'", line);
+}
+
+/* Read in the test input data from FILENAME.  */
+
+static void
+read_input (const char *filename)
+{
+  FILE *fp = fopen (filename, "r");
+  if (fp == NULL)
+    error (EXIT_FAILURE, errno, "open '%s'", filename);
+  unsigned int lineno = 0;
+  for (;;)
+    {
+      size_t size = 0;
+      char *line = NULL;
+      ssize_t ret = getline (&line, &size, fp);
+      if (ret == -1)
+	break;
+      lineno++;
+      if (line[0] == '#' || line[0] == '\n')
+	continue;
+      add_test (line, filename, lineno);
+    }
+  if (ferror (fp))
+    error (EXIT_FAILURE, errno, "read from '%s'", filename);
+  if (fclose (fp) != 0)
+    error (EXIT_FAILURE, errno, "close '%s'", filename);
+}
+
+/* Calculate the generic results (round-to-zero with sticky bit) for
+   the function described by CALC, with inputs INPUTS.  */
+
+static void
+calc_generic_results (generic_value *outputs, generic_value *inputs,
+		      const func_calc_desc *calc)
+{
+  bool inexact;
+  switch (calc->method)
+    {
+    case mpfr_f_f:
+      assert (inputs[0].type == gtype_fp);
+      outputs[0].type = gtype_fp;
+      mpfr_init (outputs[0].value.f);
+      inexact = calc->func.mpfr_f_f (outputs[0].value.f, inputs[0].value.f,
+				     MPFR_RNDZ);
+      adjust_real (outputs[0].value.f, inexact);
+      break;
+
+    default:
+      abort ();
+    }
+}
+
+/* Return the number of bits for integer type TYPE, where "long" has
+   LONG_BITS bits (32 or 64).  */
+
+static int
+int_type_bits (arg_ret_type type, int long_bits)
+{
+  assert (long_bits == 32 || long_bits == 64);
+  switch (type)
+    {
+    case type_int:
+      return 32;
+      break;
+
+    case type_long:
+      return long_bits;
+      break;
+
+    case type_long_long:
+      return 64;
+      break;
+
+    default:
+      abort ();
+    }
+}
+
+/* Check whether an integer Z fits a given type TYPE, where "long" has
+   LONG_BITS bits (32 or 64).  */
+
+static bool
+int_fits_type (mpz_t z, arg_ret_type type, int long_bits)
+{
+  int bits = int_type_bits (type, long_bits);
+  bool ret = true;
+  mpz_t t;
+  mpz_init (t);
+  mpz_ui_pow_ui (t, 2, bits - 1);
+  if (mpz_cmp (z, t) >= 0)
+    ret = false;
+  mpz_neg (t, t);
+  if (mpz_cmp (z, t) < 0)
+    ret = false;
+  mpz_clear (t);
+  return ret;
+}
+
+/* Print a generic value V to FP (name FILENAME), preceded by a space,
+   for type TYPE, floating-point format FORMAT, LONG_BITS bits per
+   long, printing " IGNORE" instead if IGNORE.  */
+
+static void
+output_generic_value (FILE *fp, const char *filename, const generic_value *v,
+		      bool ignore, arg_ret_type type, fp_format format,
+		      int long_bits)
+{
+  if (ignore)
+    {
+      if (fputs (" IGNORE", fp) < 0)
+	error (EXIT_FAILURE, errno, "write to '%s'", filename);
+      return;
+    }
+  assert (v->type == generic_arg_ret_type (type));
+  const char *suffix;
+  switch (type)
+    {
+    case type_fp:
+      suffix = fp_formats[format].suffix;
+      break;
+
+    case type_int:
+      suffix = "";
+      break;
+
+    case type_long:
+      suffix = "L";
+      break;
+
+    case type_long_long:
+      suffix = "LL";
+      break;
+
+    default:
+      abort ();
+    }
+  switch (v->type)
+    {
+    case gtype_fp:
+      if (mpfr_inf_p (v->value.f))
+	{
+	  if (fputs ((mpfr_signbit (v->value.f)
+		      ? " minus_infty" : " plus_infty"), fp) < 0)
+	    error (EXIT_FAILURE, errno, "write to '%s'", filename);
+	}
+      else
+	{
+	  assert (mpfr_number_p (v->value.f));
+	  if (mpfr_fprintf (fp, " %Ra%s", v->value.f, suffix) < 0)
+	    error (EXIT_FAILURE, errno, "mpfr_fprintf to '%s'", filename);
+	}
+      break;
+
+    case gtype_int: ;
+      int bits = int_type_bits (type, long_bits);
+      mpz_t tmp;
+      mpz_init (tmp);
+      mpz_ui_pow_ui (tmp, 2, bits - 1);
+      mpz_neg (tmp, tmp);
+      if (mpz_cmp (v->value.i, tmp) == 0)
+	{
+	  mpz_add_ui (tmp, tmp, 1);
+	  if (mpfr_fprintf (fp, " (%Zd%s-1)", tmp, suffix) < 0)
+	    error (EXIT_FAILURE, errno, "mpfr_fprintf to '%s'", filename);
+	}
+      else
+	{
+	  if (mpfr_fprintf (fp, " %Zd%s", v->value.i, suffix) < 0)
+	    error (EXIT_FAILURE, errno, "mpfr_fprintf to '%s'", filename);
+	}
+      mpz_clear (tmp);
+      break;
+
+    default:
+      abort ();
+    }
+}
+
+/* Generate test output to FP (name FILENAME) for test function TF,
+   input test IT, choice of input values INPUTS.  */
+
+static void
+output_for_one_input_case (FILE *fp, const char *filename, test_function *tf,
+			   input_test *it, generic_value *inputs)
+{
+  bool long_bits_matters = false;
+  bool fits_long32 = true;
+  for (size_t i = 0; i < tf->num_args; i++)
+    {
+      generic_value_type gtype = generic_arg_ret_type (tf->arg_types[i]);
+      assert (inputs[i].type == gtype);
+      if (gtype == gtype_int)
+	{
+	  bool fits_64 = int_fits_type (inputs[i].value.i, tf->arg_types[i],
+					64);
+	  if (!fits_64)
+	    return;
+	  if (tf->arg_types[i] == type_long
+	      && !int_fits_type (inputs[i].value.i, tf->arg_types[i], 32))
+	    {
+	      long_bits_matters = true;
+	      fits_long32 = false;
+	    }
+	}
+    }
+  generic_value generic_outputs[MAX_NRET];
+  calc_generic_results (generic_outputs, inputs, &tf->calc);
+  bool ignore_output_long32[MAX_NRET] = { false };
+  bool ignore_output_long64[MAX_NRET] = { false };
+  for (size_t i = 0; i < tf->num_ret; i++)
+    {
+      assert (generic_outputs[i].type
+	      == generic_arg_ret_type (tf->ret_types[i]));
+      switch (generic_outputs[i].type)
+	{
+	case gtype_fp:
+	  if (!mpfr_number_p (generic_outputs[i].value.f))
+	    goto out; /* Result is NaN or exact infinity.  */
+	  break;
+
+	case gtype_int:
+	  ignore_output_long32[i] = !int_fits_type (generic_outputs[i].value.i,
+						    tf->ret_types[i], 32);
+	  ignore_output_long64[i] = !int_fits_type (generic_outputs[i].value.i,
+						    tf->ret_types[i], 64);
+	  if (ignore_output_long32[i] != ignore_output_long64[i])
+	    long_bits_matters = true;
+	  break;
+
+	default:
+	  abort ();
+	}
+    }
+  /* Iterate over relevant sizes of long and floating-point formats.  */
+  for (int long_bits = 32; long_bits <= 64; long_bits += 32)
+    {
+      if (long_bits == 32 && !fits_long32)
+	continue;
+      if (long_bits == 64 && !long_bits_matters)
+	continue;
+      const char *long_cond;
+      if (long_bits_matters)
+	long_cond = (long_bits == 32 ? ":long32" : ":long64");
+      else
+	long_cond = "";
+      bool *ignore_output = (long_bits == 32
+			     ? ignore_output_long32
+			     : ignore_output_long64);
+      for (fp_format f = fp_first_format; f < fp_num_formats; f++)
+	{
+	  bool fits = true;
+	  mpfr_t res[rm_num_modes];
+	  unsigned int exc_before[rm_num_modes];
+	  unsigned int exc_after[rm_num_modes];
+	  for (size_t i = 0; i < tf->num_args; i++)
+	    {
+	      if (inputs[i].type == gtype_fp)
+		round_real (res, exc_before, exc_after, inputs[i].value.f, f);
+	      if (!mpfr_equal_p (res[rm_tonearest], inputs[i].value.f))
+		fits = false;
+	      for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
+		mpfr_clear (res[m]);
+	      if (!fits)
+		break;
+	    }
+	  if (!fits)
+	    continue;
+	  /* The inputs fit this type, so compute the ideal outputs
+	     and exceptions.  */
+	  mpfr_t all_res[MAX_NRET][rm_num_modes];
+	  unsigned int all_exc_before[MAX_NRET][rm_num_modes];
+	  unsigned int all_exc_after[MAX_NRET][rm_num_modes];
+	  unsigned int merged_exc_before[rm_num_modes] = { 0 };
+	  unsigned int merged_exc_after[rm_num_modes] = { 0 };
+	  /* For functions not exactly determined, track whether
+	     underflow is required (some result is inexact, and
+	     magnitude does not exceed the greatest magnitude
+	     subnormal), and permitted (not an exact zero, and
+	     magnitude does not exceed the least magnitude
+	     normal).  */
+	  bool must_underflow = false;
+	  bool may_underflow = false;
+	  for (size_t i = 0; i < tf->num_ret; i++)
+	    {
+	      switch (generic_outputs[i].type)
+		{
+		case gtype_fp:
+		  round_real (all_res[i], all_exc_before[i], all_exc_after[i],
+			      generic_outputs[i].value.f, f);
+		  for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
+		    {
+		      merged_exc_before[m] |= all_exc_before[i][m];
+		      merged_exc_after[m] |= all_exc_after[i][m];
+		      if (!tf->exact)
+			{
+			  must_underflow
+			    |= ((all_exc_before[i][m]
+				 & (1U << exc_inexact)) != 0
+				&& (mpfr_cmpabs (generic_outputs[i].value.f,
+						fp_formats[f].subnorm_max)
+				    <= 0));
+			  may_underflow
+			    |= (!mpfr_zero_p (generic_outputs[i].value.f)
+				&& mpfr_cmpabs (generic_outputs[i].value.f,
+						fp_formats[f].min) <= 0);
+			}
+		    }
+		  break;
+
+		case gtype_int:
+		  if (ignore_output[i])
+		    for (rounding_mode m = rm_first_mode;
+			 m < rm_num_modes;
+			 m++)
+		      {
+			merged_exc_before[m] |= 1U << exc_invalid;
+			merged_exc_after[m] |= 1U << exc_invalid;
+		      }
+		  break;
+
+		default:
+		  abort ();
+		}
+	    }
+	  assert (may_underflow || !must_underflow);
+	  for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
+	    {
+	      bool before_after_matters
+		= tf->exact && merged_exc_before[m] != merged_exc_after[m];
+	      for (int after = 0; after <= 1; after++)
+		{
+		  if (after == 1 && !before_after_matters)
+		    continue;
+		  const char *after_cond;
+		  if (before_after_matters)
+		    after_cond = (after
+				  ? ":after-rounding"
+				  : ":before-rounding");
+		  else
+		    after_cond = "";
+		  unsigned int merged_exc = (after
+					     ? merged_exc_after[m]
+					     : merged_exc_before[m]);
+		  if (fprintf (fp, "= %s %s %s%s%s", tf->name,
+			       rounding_modes[m].name, fp_formats[f].name,
+			       long_cond, after_cond) < 0)
+		    error (EXIT_FAILURE, errno, "write to '%s'", filename);
+		  /* Print inputs.  */
+		  for (size_t i = 0; i < tf->num_args; i++)
+		    output_generic_value (fp, filename, &inputs[i], false,
+					  tf->arg_types[i], f, long_bits);
+		  if (fputs (" :", fp) < 0)
+		    error (EXIT_FAILURE, errno, "write to '%s'", filename);
+		  /* Print outputs.  */
+		  bool must_erange = false;
+		  for (size_t i = 0; i < tf->num_ret; i++)
+		    {
+		      generic_value g;
+		      g.type = generic_outputs[i].type;
+		      switch (g.type)
+			{
+			case gtype_fp:
+			  if (mpfr_inf_p (all_res[i][m])
+			      && (all_exc_before[i][m]
+				  & (1U << exc_overflow)) != 0)
+			    must_erange = true;
+			  if (mpfr_zero_p (all_res[i][m])
+			      && (all_exc_before[i][m]
+				  & (1U << exc_underflow)) != 0)
+			    must_erange = true;
+			  mpfr_init2 (g.value.f, fp_formats[f].mant_dig);
+			  assert_exact (mpfr_set (g.value.f, all_res[i][m],
+						  MPFR_RNDN));
+			  break;
+
+			case gtype_int:
+			  mpz_init (g.value.i);
+			  mpz_set (g.value.i, generic_outputs[i].value.i);
+			  break;
+
+			default:
+			  abort ();
+			}
+		      output_generic_value (fp, filename, &g, ignore_output[i],
+					    tf->ret_types[i], f, long_bits);
+		      generic_value_free (&g);
+		    }
+		  if (fputs (" :", fp) < 0)
+		    error (EXIT_FAILURE, errno, "write to '%s'", filename);
+		  /* Print miscellaneous flags (passed through from
+		     input).  */
+		  for (size_t i = 0; i < it->num_flags; i++)
+		    switch (it->flags[i].type)
+		      {
+		      case flag_no_test_inline:
+		      case flag_xfail:
+			if (fprintf (fp, " %s%s",
+				     input_flags[it->flags[i].type],
+				     (it->flags[i].cond
+				      ? it->flags[i].cond
+				      : "")) < 0)
+			  error (EXIT_FAILURE, errno, "write to '%s'",
+				 filename);
+			break;
+		      default:
+			break;
+		      }
+		  /* Print exception flags and compute errno
+		     expectations where not already computed.  */
+		  bool may_edom = false;
+		  bool must_edom = false;
+		  bool may_erange = must_erange || may_underflow;
+		  for (fp_exception e = exc_first_exception;
+		       e < exc_num_exceptions;
+		       e++)
+		    {
+		      bool expect_e = (merged_exc & (1U << e)) != 0;
+		      bool e_optional = false;
+		      switch (e)
+			{
+			case exc_divbyzero:
+			  if (expect_e)
+			    may_erange = must_erange = true;
+			  break;
+
+			case exc_inexact:
+			  if (!tf->exact)
+			    e_optional = true;
+			  break;
+
+			case exc_invalid:
+			  if (expect_e)
+			    may_edom = must_edom = true;
+			  break;
+
+			case exc_overflow:
+			  if (expect_e)
+			    may_erange = true;
+			  break;
+
+			case exc_underflow:
+			  if (expect_e)
+			    may_erange = true;
+			  if (must_underflow)
+			    assert (expect_e);
+			  if (may_underflow && !must_underflow)
+			    e_optional = true;
+			  break;
+
+			default:
+			  abort ();
+			}
+		      if (e_optional)
+			{
+			  if (fprintf (fp, " %s-ok", exceptions[e]) < 0)
+			    error (EXIT_FAILURE, errno, "write to '%s'",
+				   filename);
+			}
+		      else
+			{
+			  if (expect_e)
+			    if (fprintf (fp, " %s", exceptions[e]) < 0)
+			      error (EXIT_FAILURE, errno, "write to '%s'",
+				     filename);
+			  input_flag_type okflag;
+			  okflag = (expect_e
+				    ? flag_missing_first
+				    : flag_spurious_first) + e;
+			  for (size_t i = 0; i < it->num_flags; i++)
+			    if (it->flags[i].type == okflag)
+			      if (fprintf (fp, " %s-ok%s",
+					   exceptions[e],
+					   (it->flags[i].cond
+					    ? it->flags[i].cond
+					    : "")) < 0)
+				error (EXIT_FAILURE, errno, "write to '%s'",
+				       filename);
+			}
+		    }
+		  /* Print errno expectations.  */
+		  if (tf->complex_fn)
+		    {
+		      must_edom = false;
+		      must_erange = false;
+		    }
+		  if (may_edom && !must_edom)
+		    {
+		      if (fputs (" errno-edom-ok", fp) < 0)
+			error (EXIT_FAILURE, errno, "write to '%s'",
+			       filename);
+		    }
+		  else
+		    {
+		      if (must_edom)
+			if (fputs (" errno-edom", fp) < 0)
+			  error (EXIT_FAILURE, errno, "write to '%s'",
+				 filename);
+		      input_flag_type okflag = (must_edom
+						? flag_missing_errno
+						: flag_spurious_errno);
+		      for (size_t i = 0; i < it->num_flags; i++)
+			if (it->flags[i].type == okflag)
+			  if (fprintf (fp, " errno-edom-ok%s",
+				       (it->flags[i].cond
+					? it->flags[i].cond
+					: "")) < 0)
+			    error (EXIT_FAILURE, errno, "write to '%s'",
+				   filename);
+		    }
+		  if (may_erange && !must_erange)
+		    {
+		      if (fputs (" errno-erange-ok", fp) < 0)
+			error (EXIT_FAILURE, errno, "write to '%s'",
+			       filename);
+		    }
+		  else
+		    {
+		      if (must_erange)
+			if (fputs (" errno-erange", fp) < 0)
+			  error (EXIT_FAILURE, errno, "write to '%s'",
+				 filename);
+		      input_flag_type okflag = (must_erange
+						? flag_missing_errno
+						: flag_spurious_errno);
+		      for (size_t i = 0; i < it->num_flags; i++)
+			if (it->flags[i].type == okflag)
+			  if (fprintf (fp, " errno-erange-ok%s",
+				       (it->flags[i].cond
+					? it->flags[i].cond
+					: "")) < 0)
+			    error (EXIT_FAILURE, errno, "write to '%s'",
+				   filename);
+		    }
+		  if (putc ('\n', fp) < 0)
+		    error (EXIT_FAILURE, errno, "write to '%s'", filename);
+		}
+	    }
+	  for (size_t i = 0; i < tf->num_ret; i++)
+	    {
+	      if (generic_outputs[i].type == gtype_fp)
+		for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
+		  mpfr_clear (all_res[i][m]);
+	    }
+	}
+    }
+ out:
+  for (size_t i = 0; i < tf->num_ret; i++)
+    generic_value_free (&generic_outputs[i]);
+}
+
+/* Generate test output data to FILENAME.  */
+
+static void
+generate_output (const char *filename)
+{
+  FILE *fp = fopen (filename, "w");
+  if (fp == NULL)
+    error (EXIT_FAILURE, errno, "open '%s'", filename);
+  for (size_t i = 0; i < ARRAY_SIZE (test_functions); i++)
+    {
+      test_function *tf = &test_functions[i];
+      for (size_t j = 0; j < tf->num_tests; j++)
+	{
+	  input_test *it = &tf->tests[j];
+	  if (fputs (it->line, fp) < 0)
+	    error (EXIT_FAILURE, errno, "write to '%s'", filename);
+	  for (size_t k = 0; k < it->num_input_cases; k++)
+	    output_for_one_input_case (fp, filename, tf, it, it->inputs[k]);
+	}
+    }
+  if (fclose (fp) != 0)
+    error (EXIT_FAILURE, errno, "close '%s'", filename);
+}
+
+int
+main (int argc, char **argv)
+{
+  if (argc != 3)
+    error (EXIT_FAILURE, 0, "usage: gen-auto-libm-tests <input> <output>");
+  const char *input_filename = argv[1];
+  const char *output_filename = argv[2];
+  init_fp_formats ();
+  read_input (input_filename);
+  generate_output (output_filename);
+  exit (EXIT_SUCCESS);
+}
diff --git a/math/gen-libm-test.pl b/math/gen-libm-test.pl
index bd701d7d84..81f5fabdd7 100755
--- a/math/gen-libm-test.pl
+++ b/math/gen-libm-test.pl
@@ -39,10 +39,11 @@ use Getopt::Std;
 
 use strict;
 
-use vars qw ($input $output);
+use vars qw ($input $output $auto_input);
 use vars qw (%results);
 use vars qw (%beautify @all_floats);
 use vars qw ($output_dir $ulps_file);
+use vars qw (%auto_tests);
 
 # all_floats is sorted and contains all recognised float types
 @all_floats = ('double', 'float', 'idouble',
@@ -51,6 +52,12 @@ use vars qw ($output_dir $ulps_file);
 %beautify =
   ( "minus_zero" => "-0",
     "plus_zero" => "+0",
+    "-0x0p+0f" => "-0",
+    "-0x0p+0" => "-0",
+    "-0x0p+0L" => "-0",
+    "0x0p+0f" => "+0",
+    "0x0p+0" => "+0",
+    "0x0p+0L" => "+0",
     "minus_infty" => "-inf",
     "plus_infty" => "inf",
     "qnan_value" => "qNaN",
@@ -99,9 +106,11 @@ $ulps_file = $opt_u if ($opt_u);
 $output_dir = $opt_o if ($opt_o);
 
 $input = "libm-test.inc";
+$auto_input = "auto-libm-test-out";
 $output = "${output_dir}libm-test.c";
 
 &parse_ulps ($ulps_file);
+&parse_auto_input ($auto_input);
 &generate_testfile ($input, $output) unless ($opt_n);
 &output_ulps ("${output_dir}libm-test-ulps.h", $ulps_file) unless ($opt_n);
 &print_ulps_file ("${output_dir}NewUlps") if ($opt_n);
@@ -121,6 +130,9 @@ sub beautify {
       return '-' . $beautify{$tmp};
     }
   }
+  if ($arg =~ /^-?0x[0-9a-f.]*p[-+][0-9]+f$/) {
+    $arg =~ s/f$//;
+  }
   if ($arg =~ /[0-9]L$/) {
     $arg =~ s/L$//;
   }
@@ -218,7 +230,7 @@ sub parse_args {
   # consistency check
   if ($current_arg == $#args) {
     die ("wrong number of arguments")
-      unless ($args[$current_arg] =~ /EXCEPTION|ERRNO|IGNORE_ZERO_INF_SIGN|TEST_NAN_SIGN|NO_TEST_INLINE/);
+      unless ($args[$current_arg] =~ /EXCEPTION|ERRNO|IGNORE_ZERO_INF_SIGN|TEST_NAN_SIGN|NO_TEST_INLINE|XFAIL_TEST/);
   } elsif ($current_arg < $#args) {
     die ("wrong number of arguments");
   } elsif ($current_arg > ($#args+1)) {
@@ -309,17 +321,165 @@ sub parse_args {
   print $file "    $cline },\n";
 }
 
+# Convert a condition from auto-libm-test-out to C form.
+sub convert_condition {
+  my ($cond) = @_;
+  my (@conds, $ret);
+  @conds = split /:/, $cond;
+  foreach (@conds) {
+    s/-/_/g;
+    s/^/TEST_COND_/;
+  }
+  $ret = join " && ", @conds;
+  return "($ret)";
+}
+
+# Return text to OR a value into an accumulated flags string.
+sub or_value {
+  my ($cond) = @_;
+  if ($cond eq "0") {
+    return "";
+  } else {
+    return " | $cond";
+  }
+}
+
+# Return text to OR a conditional expression between two values into
+# an accumulated flags string.
+sub or_cond_value {
+  my ($cond, $if, $else) = @_;
+  if ($cond eq "1") {
+    return or_value ($if);
+  } elsif ($cond eq "0") {
+    return or_value ($else);
+  } else {
+    return or_value ("($cond ? $if : $else)");
+  }
+}
+
 # Generate libm-test.c
 sub generate_testfile {
   my ($input, $output) = @_;
-  my ($lasttext);
-  my (@args, $i);
 
   open INPUT, $input or die ("Can't open $input: $!");
   open OUTPUT, ">$output" or die ("Can't open $output: $!");
 
   # Replace the special macros
   while (<INPUT>) {
+    # AUTO_TESTS (function, mode),
+    if (/^\s*AUTO_TESTS_/) {
+      my ($descr, $func, $mode, $auto_test, $num_auto_tests);
+      ($descr, $func, $mode) = ($_ =~ /AUTO_TESTS_(\w+)\s*\((\w+),\s*(\w+)\)/);
+      $num_auto_tests = 0;
+      foreach $auto_test (sort keys %{$auto_tests{$func}{$mode}}) {
+	my ($finputs, $format, $inputs, $outputs, $flags);
+	my ($format_conv, $flags_conv, @flags, %flag_cond);
+	$num_auto_tests++;
+	($finputs, $outputs, $flags) = split / : */, $auto_test;
+	($format, $inputs) = split / /, $finputs, 2;
+	$inputs =~ s/ /, /g;
+	$outputs =~ s/ /, /g;
+	$format_conv = convert_condition ($format);
+	print OUTPUT "#if $format_conv\n";
+	@flags = split / /, $flags;
+	foreach (@flags) {
+	  if (/^([^:]*):(.*)$/) {
+	    my ($flag, $cond);
+	    $flag = $1;
+	    $cond = convert_condition ($2);
+	    if (defined ($flag_cond{$flag})) {
+	      if ($flag_cond{$flag} ne "1") {
+		$flag_cond{$flag} .= " || $cond";
+	      }
+	    } else {
+	      $flag_cond{$flag} = $cond;
+	    }
+	  } else {
+	    $flag_cond{$_} = "1";
+	  }
+	}
+	$flags_conv = "";
+	if (defined ($flag_cond{"no-test-inline"})) {
+	  $flags_conv .= or_cond_value ($flag_cond{"no-test-inline"},
+					"NO_TEST_INLINE", "0");
+	}
+	if (defined ($flag_cond{"xfail"})) {
+	  $flags_conv .= or_cond_value ($flag_cond{"xfail"},
+					"XFAIL_TEST", "0");
+	}
+	my (@exc_list) = qw(divbyzero inexact invalid overflow underflow);
+	my ($exc);
+	foreach $exc (@exc_list) {
+	  my ($exc_expected, $exc_ok, $no_exc);
+	  $exc_expected = "\U$exc\E_EXCEPTION";
+	  $exc_ok = "\U$exc\E_EXCEPTION_OK";
+	  $no_exc = "0";
+	  if ($exc eq "inexact") {
+	    $exc_ok = "0";
+	    $no_exc = "NO_INEXACT_EXCEPTION";
+	  }
+	  if (defined ($flag_cond{$exc})) {
+	    if ($flag_cond{$exc} ne "1") {
+	      die ("unexpected condition for $exc\n");
+	    }
+	    if (defined ($flag_cond{"$exc-ok"})) {
+	      $flags_conv .= or_cond_value ($flag_cond{"$exc-ok"},
+					    $exc_ok, $exc_expected);
+	    } else {
+	      $flags_conv .= or_value ($exc_expected);
+	    }
+	  } else {
+	    if (defined ($flag_cond{"$exc-ok"})) {
+	      $flags_conv .= or_cond_value ($flag_cond{"$exc-ok"},
+					    $exc_ok, $no_exc);
+	    } else {
+	      $flags_conv .= or_value ($no_exc);
+	    }
+	  }
+	}
+	my ($errno_expected, $errno_unknown_cond);
+	if (defined ($flag_cond{"errno-edom"})) {
+	  if ($flag_cond{"errno-edom"} ne "1") {
+	    die ("unexpected condition for errno-edom");
+	  }
+	  if (defined ($flag_cond{"errno-erange"})) {
+	    die ("multiple errno values expected");
+	  }
+	  $errno_expected = "ERRNO_EDOM";
+	} elsif (defined ($flag_cond{"errno-erange"})) {
+	  if ($flag_cond{"errno-erange"} ne "1") {
+	    die ("unexpected condition for errno-erange");
+	  }
+	  $errno_expected = "ERRNO_ERANGE";
+	} else {
+	  $errno_expected = "ERRNO_UNCHANGED";
+	}
+	if (defined ($flag_cond{"errno-edom-ok"})) {
+	  if (defined ($flag_cond{"errno-erange-ok"})
+	      && $flag_cond{"errno-erange-ok"} ne $flag_cond{"errno-edom-ok"}) {
+	    $errno_unknown_cond = "($flag_cond{\"errno-edom-ok\"} || $flag_cond{\"errno-erange-ok\"})";
+	  } else {
+	    $errno_unknown_cond = $flag_cond{"errno-edom-ok"};
+	  }
+	} elsif (defined ($flag_cond{"errno-erange-ok"})) {
+	  $errno_unknown_cond = $flag_cond{"errno-erange-ok"};
+	} else {
+	  $errno_unknown_cond = "0";
+	}
+	$flags_conv .= or_cond_value ($errno_unknown_cond,
+				      "0", $errno_expected);
+	if ($flags_conv ne "") {
+	  $flags_conv =~ s/^ \|/,/;
+	}
+	&parse_args (\*OUTPUT, $descr,
+		     "$func, $inputs, $outputs$flags_conv");
+	print OUTPUT "#endif\n";
+      }
+      if ($num_auto_tests == 0) {
+	die ("no automatic tests for $func, $mode\n");
+      }
+      next;
+    }
 
     # TEST_...
     if (/^\s*TEST_/) {
@@ -573,3 +733,20 @@ sub output_ulps {
   print ULP "  };\n";
   close ULP;
 }
+
+# Parse auto-libm-test-out.
+sub parse_auto_input {
+  my ($file) = @_;
+  open AUTO, $file or die ("Can't open $file: $!");
+  while (<AUTO>) {
+    chop;
+    next if !/^= /;
+    s/^= //;
+    if (/^(\S+) (\S+) (.*)$/) {
+      $auto_tests{$1}{$2}{$3} = 1;
+    } else {
+      die ("bad automatic test line: $_\n");
+    }
+  }
+  close AUTO;
+}
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 0c88abc3db..017fa63122 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -175,12 +175,13 @@ struct ulp_data
 #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			0x2000
-#define ERRNO_EDOM			0x4000
-#define ERRNO_ERANGE			0x8000
+#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			0x10000
+#define IGNORE_RESULT			0x20000
 
 /* Values underflowing only for float.  */
 #ifdef TEST_FLOAT
@@ -225,6 +226,47 @@ struct ulp_data
 # define NO_TEST_INLINE_DOUBLE	0
 #endif
 
+/* Conditions used by tests generated by gen-auto-libm-tests.c.  */
+#ifdef TEST_FLOAT
+# define TEST_COND_flt_32	1
+#else
+# define TEST_COND_flt_32	0
+#endif
+#if defined TEST_DOUBLE || (defined TEST_LDOUBLE && LDBL_MANT_DIG == 53)
+# define TEST_COND_dbl_64	1
+#else
+# define TEST_COND_dbl_64	0
+#endif
+#if defined TEST_LDOUBLE && LDBL_MANT_DIG == 64 && LDBL_MIN_EXP == -16381
+# define TEST_COND_ldbl_96_intel	1
+#else
+# define TEST_COND_ldbl_96_intel	0
+#endif
+#if defined TEST_LDOUBLE && LDBL_MANT_DIG == 64 && LDBL_MIN_EXP == -16382
+# define TEST_COND_ldbl_96_m68k	1
+#else
+# define TEST_COND_ldbl_96_m68k	0
+#endif
+#if defined TEST_LDOUBLE && LDBL_MANT_DIG == 113
+# define TEST_COND_ldbl_128	1
+#else
+# define TEST_COND_ldbl_128	0
+#endif
+#if defined TEST_LDOUBLE && LDBL_MANT_DIG == 106
+# define TEST_COND_ldbl_128ibm	1
+#else
+# define TEST_COND_ldbl_128ibm	0
+#endif
+#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 (we must supply them precalculated for accuracy).  */
 #define M_PI_6l			.52359877559829887307710723054658383L
 #define M_E2l			7.389056098930650227230427460575008L
@@ -997,6 +1039,8 @@ check_longlong (const char *test_name, long long int computed,
 static int
 enable_test (int exceptions)
 {
+  if (exceptions & XFAIL_TEST)
+    return 0;
 #ifdef TEST_INLINE
   if (exceptions & NO_TEST_INLINE)
     return 0;
@@ -13707,24 +13751,15 @@ sinh_test_upward (void)
 
 static const struct test_f_f_data sqrt_test_data[] =
   {
-    TEST_f_f (sqrt, 0, 0, NO_INEXACT_EXCEPTION),
     TEST_f_f (sqrt, qnan_value, qnan_value, NO_INEXACT_EXCEPTION),
     TEST_f_f (sqrt, plus_infty, plus_infty, NO_INEXACT_EXCEPTION),
 
-    TEST_f_f (sqrt, minus_zero, minus_zero, NO_INEXACT_EXCEPTION),
-
     /* sqrt (x) == qNaN plus invalid exception for x < 0.  */
     TEST_f_f (sqrt, -1, qnan_value, NO_INEXACT_EXCEPTION|INVALID_EXCEPTION|ERRNO_EDOM),
     TEST_f_f (sqrt, -max_value, qnan_value, NO_INEXACT_EXCEPTION|INVALID_EXCEPTION|ERRNO_EDOM),
     TEST_f_f (sqrt, minus_infty, qnan_value, NO_INEXACT_EXCEPTION|INVALID_EXCEPTION|ERRNO_EDOM),
 
-    TEST_f_f (sqrt, 2209, 47, NO_INEXACT_EXCEPTION),
-    TEST_f_f (sqrt, 4, 2, NO_INEXACT_EXCEPTION),
-    TEST_f_f (sqrt, 2, M_SQRT2l, INEXACT_EXCEPTION),
-    TEST_f_f (sqrt, 0.25, 0.5, NO_INEXACT_EXCEPTION),
-    TEST_f_f (sqrt, 6642.25, 81.5, NO_INEXACT_EXCEPTION),
-    TEST_f_f (sqrt, 15190.5625L, 123.25L, NO_INEXACT_EXCEPTION),
-    TEST_f_f (sqrt, 0.75L, 0.866025403784438646763723170752936183L, INEXACT_EXCEPTION),
+    AUTO_TESTS_f_f (sqrt, tonearest),
   };
 
 static void