From 83975da8b13b867dc5f741b9506799cd0529ef1b Mon Sep 17 00:00:00 2001 From: Andreas Enge Date: Tue, 30 Jun 2020 16:08:05 +0200 Subject: mpcheck: Move the tool from tests/ to its own subdirectory inside tools/. * tests/mpcheck-*: Move to... * tools/mpcheck/mpcheck-*: ...this. * tests/Makefile.am: Drop all references to mpcheck. * tools/mpcheck/Makefile.am, tools/mpcheck/README: New files. * tools/Makefile.am: Reference the new subdirectory. * configure.ac: Reference tools/mpcheck/Makefile. * Makefile.am: Add a target "mpcheck". * NEWS: Add an item. * doc/version.texi: Spurious update. --- Makefile.am | 5 +- NEWS | 1 + configure.ac | 4 +- doc/version.texi | 8 +- tests/Makefile.am | 7 +- tests/mpcheck-common.c | 193 ---------------------------- tests/mpcheck-double.c | 241 ----------------------------------- tests/mpcheck-float.c | 246 ------------------------------------ tests/mpcheck-float128.c | 251 ------------------------------------ tests/mpcheck-longdouble.c | 234 ---------------------------------- tests/mpcheck-template1.c | 155 ----------------------- tests/mpcheck-template2.c | 130 ------------------- tests/mpcheck-template3.c | 172 ------------------------- tools/Makefile.am | 3 +- tools/mpcheck/Makefile.am | 31 +++++ tools/mpcheck/README | 3 + tools/mpcheck/mpcheck-common.c | 193 ++++++++++++++++++++++++++++ tools/mpcheck/mpcheck-double.c | 242 +++++++++++++++++++++++++++++++++++ tools/mpcheck/mpcheck-float.c | 247 ++++++++++++++++++++++++++++++++++++ tools/mpcheck/mpcheck-float128.c | 252 +++++++++++++++++++++++++++++++++++++ tools/mpcheck/mpcheck-longdouble.c | 235 ++++++++++++++++++++++++++++++++++ tools/mpcheck/mpcheck-template1.c | 155 +++++++++++++++++++++++ tools/mpcheck/mpcheck-template2.c | 130 +++++++++++++++++++ tools/mpcheck/mpcheck-template3.c | 172 +++++++++++++++++++++++++ 24 files changed, 1675 insertions(+), 1635 deletions(-) delete mode 100644 tests/mpcheck-common.c delete mode 100644 tests/mpcheck-double.c delete mode 100644 tests/mpcheck-float.c delete mode 100644 tests/mpcheck-float128.c delete mode 100644 tests/mpcheck-longdouble.c delete mode 100644 tests/mpcheck-template1.c delete mode 100644 tests/mpcheck-template2.c delete mode 100644 tests/mpcheck-template3.c create mode 100644 tools/mpcheck/Makefile.am create mode 100644 tools/mpcheck/README create mode 100644 tools/mpcheck/mpcheck-common.c create mode 100644 tools/mpcheck/mpcheck-double.c create mode 100644 tools/mpcheck/mpcheck-float.c create mode 100644 tools/mpcheck/mpcheck-float128.c create mode 100644 tools/mpcheck/mpcheck-longdouble.c create mode 100644 tools/mpcheck/mpcheck-template1.c create mode 100644 tools/mpcheck/mpcheck-template2.c create mode 100644 tools/mpcheck/mpcheck-template3.c diff --git a/Makefile.am b/Makefile.am index a286ce7..be2ade6 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,6 +1,6 @@ ## Makefile.am -- Process this file with automake to produce Makefile.in ## -## Copyright (C) 2008, 2010, 2011, 2012, 2013, 2014 INRIA +## Copyright (C) 2008, 2010, 2011, 2012, 2013, 2014, 2020 INRIA ## ## This file is part of GNU MPC. ## @@ -32,3 +32,6 @@ EXTRA_DIST = doc/fdl-1.3.texi src/mpc-log.h Makefile.vc bench : cd tools/bench && $(MAKE) $(AM_MAKEFLAGS) bench +mpcheck : + cd tools/mpcheck && $(MAKE) $(AM_MAKEFLAGS) mpcheck + diff --git a/NEWS b/NEWS index 528ad71..b0c2aa6 100644 --- a/NEWS +++ b/NEWS @@ -2,6 +2,7 @@ Changes in version xxx: - New functions: mpc_sum, mpc_dot - Several functions are more robust with a reduced exponent range (for example corresponding to IEEE 754 binary formats) + - New tool mpcheck. Changes in version 1.1.0: - Minimally required library versions: GMP 5.0.0 and MPFR 3.0.0 diff --git a/configure.ac b/configure.ac index a97485c..9f8069b 100644 --- a/configure.ac +++ b/configure.ac @@ -1,4 +1,4 @@ -# Copyright (C) 2008, 2009, 2010, 2011, 2012, 2014, 2016, 2017, 2018 INRIA +# Copyright (C) 2008, 2009, 2010, 2011, 2012, 2014, 2016, 2017, 2018, 2020 INRIA # # This file is part of GNU MPC. # @@ -253,5 +253,5 @@ AC_DEFINE_UNQUOTED([MPC_GCC_VERSION], ["$GCC_VERSION"], [Version of gcc]) # Looks for short git hash if the version string contains "dev" MPC_GITVERSION -AC_CONFIG_FILES([Makefile src/Makefile tests/Makefile doc/Makefile tools/Makefile tools/bench/Makefile]) +AC_CONFIG_FILES([Makefile src/Makefile tests/Makefile doc/Makefile tools/Makefile tools/bench/Makefile tools/mpcheck/Makefile]) AC_OUTPUT diff --git a/doc/version.texi b/doc/version.texi index 56a933c..2fdcf8b 100644 --- a/doc/version.texi +++ b/doc/version.texi @@ -1,4 +1,4 @@ -@set UPDATED 11 January 2018 -@set UPDATED-MONTH January 2018 -@set EDITION 1.1.0 -@set VERSION 1.1.0 +@set UPDATED 16 January 2019 +@set UPDATED-MONTH January 2019 +@set EDITION 1.1.1dev +@set VERSION 1.1.1dev diff --git a/tests/Makefile.am b/tests/Makefile.am index 6f6ec61..f211d5c 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -1,6 +1,6 @@ ## tests/Makefile.am -- Process this file with automake to produce Makefile.in ## -## Copyright (C) 2008, 2009, 2010, 2011, 2012, 2013, 2016 INRIA +## Copyright (C) 2008, 2009, 2010, 2011, 2012, 2013, 2016, 2020 INRIA ## ## This file is part of GNU MPC. ## @@ -65,10 +65,7 @@ DATA_SETS = abs.dat acos.dat acosh.dat add.dat add_fr.dat arg.dat \ pow_fr.dat pow_si.dat pow_ui.dat pow_z.dat proj.dat rootofunity.dat \ sin.dat sinh.dat \ sqr.dat sqrt.dat strtoc.dat sub.dat sub_fr.dat tan.dat tanh.dat -EXTRA_DIST = data_check.tpl tgeneric.tpl $(DATA_SETS) $(DESCRIPTIONS) \ - mpcheck-float.c mpcheck-double.c mpcheck-longdouble.c \ - mpcheck-float128.c mpcheck-common.c \ - mpcheck-template1.c mpcheck-template2.c mpcheck-template3.c +EXTRA_DIST = data_check.tpl tgeneric.tpl $(DATA_SETS) $(DESCRIPTIONS) LOG_COMPILER = $(VALGRIND) AM_LOG_FLAGS = $(VALGRIND_OPTS) diff --git a/tests/mpcheck-common.c b/tests/mpcheck-common.c deleted file mode 100644 index 0133f2a..0000000 --- a/tests/mpcheck-common.c +++ /dev/null @@ -1,193 +0,0 @@ -/* shared routines for mpcheck-xxx - -Copyright (C) 2020 INRIA - -This file is part of GNU MPC. - -GNU MPC 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 3 of the License, or (at your -option) any later version. - -GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . -*/ - -#include -#include - -unsigned long Max_err_re = 0, Max_err_im = 0; -unsigned long Max_err = 0; /* for real-valued functions like abs */ -int recheck = 0; - -/* y is the reference value computed by MPC */ -static unsigned long -ulp_error (mpfr_t x, mpfr_t y) -{ - mpfr_t z; - mpfr_prec_t p = mpfr_get_prec (y); - unsigned long n, special = 0; - mpfr_exp_t e; - - if (mpfr_cmp (x, y) == 0) - return 0; - - assert (!mpfr_nan_p (x)); - assert (!mpfr_nan_p (y)); - - if (mpfr_inf_p (x) || mpfr_inf_p (y)) - { - if (mpfr_number_p (x) && mpfr_inf_p (y)) - { - /* only MPFR returns Inf: we replace y by the maximal - floating-point number */ - if (mpfr_sgn (y) > 0) /* +Inf */ - mpfr_nextbelow (y); - else - mpfr_nextabove (y); - special = 1; - } - else - { - /* only glibc returns Inf */ - if (mpfr_sgn (x) > 0) /* +Inf */ - mpfr_nextbelow (x); - else - mpfr_nextabove (x); - special = 2; - } - } - - mpfr_init2 (z, p); - mpfr_sub (z, x, y, MPFR_RNDN); - mpfr_abs (z, z, MPFR_RNDN); - /* divide by ulp(y) = 2^(EXP(y) - p) */ - e = mpfr_get_exp (y) - p; - /* take into account subnormals: smallest number is 2^(emin-1) */ - if (e < emin - 1) - e = emin - 1; - mpfr_div_2si (z, z, e, MPFR_RNDN); - n = mpfr_get_ui (z, MPFR_RNDZ); - mpfr_clear (z); - if (special == 1) /* restore original y */ - { - if (mpfr_sgn (y) > 0) - mpfr_nextabove (y); - else - mpfr_nextbelow (y); - } - if (special == 2) /* restore original x */ - { - if (mpfr_sgn (x) > 0) - mpfr_nextabove (x); - else - mpfr_nextbelow (x); - } - if (special == 0) - return n; - else - return (n < ULONG_MAX) ? n + 1 : n; -} - -static void -mpcheck_random (mpc_t z) -{ - unsigned long erange; - - mpc_urandom (z, state); - - /* set sign and exponent of the real part */ - mpz_set_ui (expz, 2 * (emax - emin + 1)); - mpz_urandomm (expz, state, expz); - erange = mpz_get_ui (expz); - /* bit 0 of erange gives the sign */ - if (erange & 1) - mpfr_neg (mpc_realref (z), mpc_realref (z), MPFR_RNDN); - mpfr_set_exp (mpc_realref (z), emin + (erange >> 1)); - /* make sure the numbers are representable in IEEE format */ - mpfr_subnormalize (mpc_realref (z), 0, MPFR_RNDN); - - /* set sign and exponent of the imaginary part */ - mpz_set_ui (expz, 2 * (emax - emin + 1)); - mpz_urandomm (expz, state, expz); - erange = mpz_get_ui (expz); - /* bit 0 of erange gives the sign */ - if (erange & 1) - mpfr_neg (mpc_imagref (z), mpc_imagref (z), MPFR_RNDN); - mpfr_set_exp (mpc_imagref (z), emin + (erange >> 1)); - /* make sure the numbers are representable in IEEE format */ - mpfr_subnormalize (mpc_imagref (z), 0, MPFR_RNDN); -} - -/* y is the reference value computed with MPC: - * return 0 if x and y are equal (including both NaN or Inf) - * return 1 if x and y disagree (neither NaN nor Inf) - * return 2 if x is NaN but not y - * return 4 if x and y are infinities of different signs - */ -static int -mpcheck_mpfr_cmp (mpfr_t x, mpfr_t y) -{ - if (mpfr_number_p (x) && mpfr_number_p (y)) - return mpfr_cmp (x, y); - /* now either x or y is NaN or Inf */ - if (mpfr_nan_p (x) || mpfr_nan_p (y)) - { - if (mpfr_nan_p (x) && mpfr_nan_p (y)) - return 0; - if (mpfr_nan_p (x)) - return 2; - fprintf (stderr, "Error, NaN and not NaN\n"); - exit (1); - } - /* now either x or y is Inf */ - if (mpfr_inf_p (x) && mpfr_inf_p (y)) - { - if (mpfr_sgn (x) * mpfr_sgn (y) > 0) - return 0; - return 4; - } - return 1; /* Inf and non-Inf: differ, will be treated by ulp_error */ -} - -/* y is the reference value computed with MPC */ -static int -mpcheck_mpc_cmp (mpc_t x, mpc_t y) -{ - return mpcheck_mpfr_cmp (mpc_realref (x), mpc_realref (y)) | - mpcheck_mpfr_cmp (mpc_imagref (x), mpc_imagref (y)); -} - -static void -report_maximal_errors (void) -{ - printf ("Maximal errors: %lu (real), %lu (imag), %lu (abs/arg)\n", - Max_err_re, Max_err_im, Max_err); -} - -/* return non-zero if both a and b are NaN, or if both are not NaN - and are equal */ -static int -mpfr_agree (mpfr_t a, mpfr_t b, int inex) -{ - if (mpfr_nan_p (a) || mpfr_nan_p (b)) - return mpfr_nan_p (a) && mpfr_nan_p (b); - else - /* in the case of even rounding, we can't tell for sure */ - return mpfr_equal_p (a, b) || inex == 2 || inex == -2; -} - -static int -cputime (void) -{ - struct rusage rus; - - getrusage (0, &rus); - return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000; -} - diff --git a/tests/mpcheck-double.c b/tests/mpcheck-double.c deleted file mode 100644 index ba85952..0000000 --- a/tests/mpcheck-double.c +++ /dev/null @@ -1,241 +0,0 @@ -/* mpcheck-double -- compare mpc functions against "double complex" - from the GNU libc implementation - -Copyright (C) 2020 INRIA - -This file is part of GNU MPC. - -GNU MPC 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 3 of the License, or (at your -option) any later version. - -GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . -*/ - -/* the GNU libc provides the following functions (as of 2.31), - with 'f' suffix for the float/binary32 version, with no suffix - for the double/binary64 version, with 'l' suffix for the long double - version, and with 'f128' suffix for the __float128 version: - - cabs casinh cexp csinh - cacos catan clog csqrt - cacosh catanh clog10 ctan - carg ccos cpow ctanh - casin ccosh csin -*/ - -#define _GNU_SOURCE /* for clog10 */ -#include -#include -#include -#include -#include -#include -#ifdef __GNUC__ -#include -#endif -#include "mpc-tests.h" - -#define PRECISION 53 -#define EMAX 1024 -#define TYPE double -#define SUFFIX - -#define mpc_get_type mpc_get_dc -#define mpc_set_type mpc_set_dc -#define mpfr_set_type mpfr_set_d - -gmp_randstate_t state; -mpz_t expz; /* global variable used in mpcheck_random */ -unsigned long seed = 0; -int verbose = 0; -mpfr_exp_t emin, emax; - -#include "mpcheck-common.c" - -#define FOO add -#define CFOO(x,y) (x+y) -#include "mpcheck-template3.c" - -#define FOO sub -#define CFOO(x,y) (x-y) -#include "mpcheck-template3.c" - -#define FOO mul -#define CFOO(x,y) (x*y) -#include "mpcheck-template3.c" - -#define FOO div -#define CFOO(x,y) (x/y) -#include "mpcheck-template3.c" - -#define FOO pow -#include "mpcheck-template3.c" - -#define FOO abs -#include "mpcheck-template2.c" - -#define FOO arg -#include "mpcheck-template2.c" - -#define FOO sqrt -#include "mpcheck-template1.c" - -#define FOO acos -#include "mpcheck-template1.c" - -#define FOO acosh -#include "mpcheck-template1.c" - -#define FOO asin -#include "mpcheck-template1.c" - -#define FOO asinh -#include "mpcheck-template1.c" - -#define FOO atan -#include "mpcheck-template1.c" - -#define FOO atanh -#include "mpcheck-template1.c" - -#define FOO cos -#include "mpcheck-template1.c" - -#define FOO cosh -#include "mpcheck-template1.c" - -#define FOO exp -#include "mpcheck-template1.c" - -#define FOO log -#include "mpcheck-template1.c" - -#define FOO log10 -#include "mpcheck-template1.c" - -#define FOO sin -#include "mpcheck-template1.c" - -#define FOO sinh -#include "mpcheck-template1.c" - -/* use reduced exponent range for tan and tanh */ -#define FOO_EMIN -8 -#define FOO_EMAX 8 - -#define FOO tan -#include "mpcheck-template1.c" - -#define FOO tanh -#include "mpcheck-template1.c" - -#undef FOO_EMIN -#undef FOO_EMAX - -int -main (int argc, char *argv[]) -{ - mpfr_prec_t p = PRECISION; /* precision of 'double' */ - unsigned long n = 1000000; /* default number of random tests per function */ - - while (argc >= 2 && argv[1][0] == '-') - { - if (argc >= 3 && strcmp (argv[1], "-p") == 0) - { - p = atoi (argv[2]); - argc -= 2; - argv += 2; - } - else if (argc >= 3 && strcmp (argv[1], "-seed") == 0) - { - seed = atoi (argv[2]); - argc -= 2; - argv += 2; - } - else if (argc >= 3 && strcmp (argv[1], "-num") == 0) - { - n = atoi (argv[2]); - argc -= 2; - argv += 2; - } - else if (strcmp (argv[1], "-v") == 0) - { - verbose ++; - argc --; - argv ++; - } - else if (strcmp (argv[1], "-check") == 0) - { - recheck = 1; - argc --; - argv ++; - } - else - { - fprintf (stderr, "Unknown option %s\n", argv[1]); - exit (1); - } - } - - /* set exponent range */ - emin = -EMAX - PRECISION + 4; /* should be -1073 */ - emax = EMAX; - mpfr_set_emin (emin); - mpfr_set_emax (emax); - - gmp_randinit_default (state); - mpz_init (expz); - -#ifdef __GNUC__ - printf ("GNU libc version: %s\n", gnu_get_libc_version ()); - printf ("GNU libc release: %s\n", gnu_get_libc_release ()); -#endif - - if (seed == 0) - seed = getpid (); - printf ("Using random seed %lu\n", seed); - - /* (complex,complex) -> complex */ - test_add (p, n); - test_sub (p, n); - test_mul (p, n); - test_div (p, n); - test_pow (p, n); - - /* complex -> real */ - test_abs (p, n); - test_arg (p, n); - - /* complex -> complex */ - test_sqrt (p, n); - test_acos (p, n); - test_acosh (p, n); - test_asin (p, n); - test_asinh (p, n); - test_atan (p, n); - test_atanh (p, n); - test_cos (p, n); - test_cosh (p, n); - test_exp (p, n); - test_log (p, n); - test_log10 (p, n); - test_sin (p, n); - test_sinh (p, n); - test_tan (p, n); - test_tanh (p, n); - - gmp_randclear (state); - mpz_clear (expz); - - report_maximal_errors (); - - return 0; -} diff --git a/tests/mpcheck-float.c b/tests/mpcheck-float.c deleted file mode 100644 index fe3cf77..0000000 --- a/tests/mpcheck-float.c +++ /dev/null @@ -1,246 +0,0 @@ -/* mpcheck-float -- compare mpc functions against "float complex" - from the GNU libc implementation - -Copyright (C) 2020 INRIA - -This file is part of GNU MPC. - -GNU MPC 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 3 of the License, or (at your -option) any later version. - -GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . -*/ - -/* the GNU libc provides the following functions (as of 2.31), - with 'f' suffix for the float/binary32 version, with no suffix - for the double/binary64 version, with 'l' suffix for the long double - version, and with 'f128' suffix for the __float128 version: - - cabs casinh cexp csinh - cacos catan clog csqrt - cacosh catanh clog10 ctan - carg ccos cpow ctanh - casin ccosh csin -*/ - -#define _GNU_SOURCE /* for clog10 */ -#include -#include -#include -#include -#include -#include -#include "mpc-tests.h" -#ifdef __GNUC__ -#include -#endif - -#define PRECISION 24 -#define EMAX 128 -#define TYPE float -#define SUFFIX f - -#define mpfr_set_type mpfr_set_flt - -static TYPE complex -mpc_get_type (mpc_t x, mpc_rnd_t rnd) -{ - /* there is no mpc_get_fltc function */ - return (TYPE complex) mpc_get_dc (x, rnd); -} - -static int -mpc_set_type (mpc_t x, TYPE complex y, mpc_rnd_t rnd) -{ - /* there is no mpc_set_fltc function */ - return mpc_set_dc (x, (double complex) y, rnd); -} - -gmp_randstate_t state; -mpz_t expz; /* global variable used in mpcheck_random */ -unsigned long seed = 0; -int verbose = 0; -mpfr_exp_t emin, emax; - -#include "mpcheck-common.c" - -#define FOO add -#define CFOO(x,y) (x+y) -#include "mpcheck-template3.c" - -#define FOO sub -#define CFOO(x,y) (x-y) -#include "mpcheck-template3.c" - -#define FOO mul -#define CFOO(x,y) (x*y) -#include "mpcheck-template3.c" - -#define FOO div -#define CFOO(x,y) (x/y) -#include "mpcheck-template3.c" - -#define FOO pow -#include "mpcheck-template3.c" - -#define FOO abs -#include "mpcheck-template2.c" - -#define FOO arg -#include "mpcheck-template2.c" - -#define FOO sqrt -#include "mpcheck-template1.c" - -#define FOO acos -#include "mpcheck-template1.c" - -#define FOO acosh -#include "mpcheck-template1.c" - -#define FOO asin -#include "mpcheck-template1.c" - -#define FOO asinh -#include "mpcheck-template1.c" - -#define FOO atan -#include "mpcheck-template1.c" - -#define FOO atanh -#include "mpcheck-template1.c" - -#define FOO cos -#include "mpcheck-template1.c" - -#define FOO cosh -#include "mpcheck-template1.c" - -#define FOO exp -#include "mpcheck-template1.c" - -#define FOO log -#include "mpcheck-template1.c" - -#define FOO log10 -#include "mpcheck-template1.c" - -#define FOO sin -#include "mpcheck-template1.c" - -#define FOO sinh -#include "mpcheck-template1.c" - -#define FOO tan -#include "mpcheck-template1.c" - -#define FOO tanh -#include "mpcheck-template1.c" - -int -main (int argc, char *argv[]) -{ - mpfr_prec_t p = PRECISION; /* precision of 'float' */ - unsigned long n = 1000000; /* default number of random tests per function */ - - while (argc >= 2 && argv[1][0] == '-') - { - if (argc >= 3 && strcmp (argv[1], "-p") == 0) - { - p = atoi (argv[2]); - argc -= 2; - argv += 2; - } - else if (argc >= 3 && strcmp (argv[1], "-seed") == 0) - { - seed = atoi (argv[2]); - argc -= 2; - argv += 2; - } - else if (argc >= 3 && strcmp (argv[1], "-num") == 0) - { - n = atoi (argv[2]); - argc -= 2; - argv += 2; - } - else if (strcmp (argv[1], "-v") == 0) - { - verbose ++; - argc --; - argv ++; - } - else if (strcmp (argv[1], "-check") == 0) - { - recheck = 1; - argc --; - argv ++; - } - else - { - fprintf (stderr, "Unknown option %s\n", argv[1]); - exit (1); - } - } - - /* set exponent range */ - emin = -EMAX - PRECISION + 4; /* should be -148 */ - emax = EMAX; - mpfr_set_emin (emin); - mpfr_set_emax (emax); - - gmp_randinit_default (state); - mpz_init (expz); - -#ifdef __GNUC__ - printf ("GNU libc version: %s\n", gnu_get_libc_version ()); - printf ("GNU libc release: %s\n", gnu_get_libc_release ()); -#endif - - if (seed == 0) - seed = getpid (); - printf ("Using random seed %lu\n", seed); - - /* (complex,complex) -> complex */ - test_add (p, n); - test_sub (p, n); - test_mul (p, n); - test_div (p, n); - test_pow (p, n); - - /* complex -> real */ - test_abs (p, n); - test_arg (p, n); - - /* complex -> complex */ - test_sqrt (p, n); - test_acos (p, n); - test_acosh (p, n); - test_asin (p, n); - test_asinh (p, n); - test_atan (p, n); - test_atanh (p, n); - test_cos (p, n); - test_cosh (p, n); - test_exp (p, n); - test_log (p, n); - test_log10 (p, n); - test_sin (p, n); - test_sinh (p, n); - test_tan (p, n); - test_tanh (p, n); - - gmp_randclear (state); - mpz_clear (expz); - - report_maximal_errors (); - - return 0; -} diff --git a/tests/mpcheck-float128.c b/tests/mpcheck-float128.c deleted file mode 100644 index 7b8f495..0000000 --- a/tests/mpcheck-float128.c +++ /dev/null @@ -1,251 +0,0 @@ -/* mpcheck-float128 -- compare mpc functions against "__float128 complex" - from the GNU libc implementation - -Copyright (C) 2020 INRIA - -This file is part of GNU MPC. - -GNU MPC 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 3 of the License, or (at your -option) any later version. - -GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . -*/ - -/* the GNU libc provides the following functions (as of 2.31), - with 'f' suffix for the float/binary32 version, with no suffix - for the double/binary64 version, with 'l' suffix for the long double - version, and with 'f128' suffix for the __float128 version: - - cabs casinh cexp csinh - cacos catan clog csqrt - cacosh catanh clog10 ctan - carg ccos cpow ctanh - casin ccosh csin -*/ - -#define _GNU_SOURCE /* for clog10 */ -#include -#include -#include -#include -#include -#include -#define MPFR_WANT_FLOAT128 -#include "mpc-tests.h" -#ifdef __GNUC__ -#include -#endif - -#define PRECISION 113 -#define EMAX 16384 -#define TYPE _Float128 -#define SUFFIX f128 - -#define mpfr_set_type mpfr_set_float128 - -static TYPE complex -mpc_get_type (mpc_t z, mpc_rnd_t rnd) -{ - TYPE x, y; - /* there is no mpc_get_float128c function */ - x = mpfr_get_float128 (mpc_realref (z), MPC_RND_RE(rnd)); - y = mpfr_get_float128 (mpc_imagref (z), MPC_RND_IM(rnd)); - return x + I * y; -} - -static int -mpc_set_type (mpc_t x, TYPE complex y, mpc_rnd_t rnd) -{ - /* there is no mpc_set_float128c function */ - mpfr_set_float128 (mpc_realref (x), crealf128 (y), MPC_RND_RE(rnd)); - mpfr_set_float128 (mpc_imagref (x), cimagf128 (y), MPC_RND_IM(rnd)); -} - -gmp_randstate_t state; -mpz_t expz; /* global variable used in mpcheck_random */ -unsigned long seed = 0; -int verbose = 0; -mpfr_exp_t emin, emax; - -#include "mpcheck-common.c" - -#define FOO add -#define CFOO(x,y) (x+y) -#include "mpcheck-template3.c" - -#define FOO sub -#define CFOO(x,y) (x-y) -#include "mpcheck-template3.c" - -#define FOO mul -#define CFOO(x,y) (x*y) -#include "mpcheck-template3.c" - -#define FOO div -#define CFOO(x,y) (x/y) -#include "mpcheck-template3.c" - -#define FOO pow -#include "mpcheck-template3.c" - -#define FOO abs -#include "mpcheck-template2.c" - -#define FOO arg -#include "mpcheck-template2.c" - -#define FOO sqrt -#include "mpcheck-template1.c" - -#define FOO acos -#include "mpcheck-template1.c" - -#define FOO acosh -#include "mpcheck-template1.c" - -#define FOO asin -#include "mpcheck-template1.c" - -#define FOO asinh -#include "mpcheck-template1.c" - -#define FOO atan -#include "mpcheck-template1.c" - -#define FOO atanh -#include "mpcheck-template1.c" - -#define FOO cos -#include "mpcheck-template1.c" - -#define FOO cosh -#include "mpcheck-template1.c" - -#define FOO exp -#include "mpcheck-template1.c" - -#define FOO log -#include "mpcheck-template1.c" - -#define FOO log10 -#include "mpcheck-template1.c" - -#define FOO sin -#include "mpcheck-template1.c" - -#define FOO sinh -#include "mpcheck-template1.c" - -#define FOO tan -#include "mpcheck-template1.c" - -#define FOO tanh -#include "mpcheck-template1.c" - -int -main (int argc, char *argv[]) -{ - mpfr_prec_t p = PRECISION; /* precision of 'double' */ - unsigned long n = 1000000; /* default number of random tests per function */ - - while (argc >= 2 && argv[1][0] == '-') - { - if (argc >= 3 && strcmp (argv[1], "-p") == 0) - { - p = atoi (argv[2]); - argc -= 2; - argv += 2; - } - else if (argc >= 3 && strcmp (argv[1], "-seed") == 0) - { - seed = atoi (argv[2]); - argc -= 2; - argv += 2; - } - else if (argc >= 3 && strcmp (argv[1], "-num") == 0) - { - n = atoi (argv[2]); - argc -= 2; - argv += 2; - } - else if (strcmp (argv[1], "-v") == 0) - { - verbose ++; - argc --; - argv ++; - } - else if (strcmp (argv[1], "-check") == 0) - { - recheck = 1; - argc --; - argv ++; - } - else - { - fprintf (stderr, "Unknown option %s\n", argv[1]); - exit (1); - } - } - - /* set exponent range */ - emin = -EMAX - 64 + 4; /* should be -16444 like for long double */ - emax = EMAX; - mpfr_set_emin (emin); - mpfr_set_emax (emax); - - gmp_randinit_default (state); - mpz_init (expz); - -#ifdef __GNUC__ - printf ("GNU libc version: %s\n", gnu_get_libc_version ()); - printf ("GNU libc release: %s\n", gnu_get_libc_release ()); -#endif - - if (seed == 0) - seed = getpid (); - printf ("Using random seed %lu\n", seed); - - /* (complex,complex) -> complex */ - test_add (p, n); - test_sub (p, n); - test_mul (p, n); - test_div (p, n); - test_pow (p, n); - - /* complex -> real */ - test_abs (p, n); - test_arg (p, n); - - /* complex -> complex */ - test_sqrt (p, n); - test_acos (p, n); - test_acosh (p, n); - test_asin (p, n); - test_asinh (p, n); - test_atan (p, n); - test_atanh (p, n); - test_cos (p, n); - test_cosh (p, n); - test_exp (p, n); - test_log (p, n); - test_log10 (p, n); - test_sin (p, n); - test_sinh (p, n); - test_tan (p, n); - test_tanh (p, n); - - gmp_randclear (state); - mpz_clear (expz); - - report_maximal_errors (); - - return 0; -} diff --git a/tests/mpcheck-longdouble.c b/tests/mpcheck-longdouble.c deleted file mode 100644 index 3dc1688..0000000 --- a/tests/mpcheck-longdouble.c +++ /dev/null @@ -1,234 +0,0 @@ -/* mpcheck-longdouble -- compare mpc functions against "long double complex" - from the GNU libc implementation - -Copyright (C) 2020 INRIA - -This file is part of GNU MPC. - -GNU MPC 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 3 of the License, or (at your -option) any later version. - -GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . -*/ - -/* the GNU libc provides the following functions (as of 2.31), - with 'f' suffix for the float/binary32 version, with no suffix - for the double/binary64 version, with 'l' suffix for the long double - version, and with 'f128' suffix for the __float128 version: - - cabs casinh cexp csinh - cacos catan clog csqrt - cacosh catanh clog10 ctan - carg ccos cpow ctanh - casin ccosh csin -*/ - -#define _GNU_SOURCE /* for clog10 */ -#include -#include -#include -#include -#include -#include -#include "mpc-tests.h" -#ifdef __GNUC__ -#include -#endif - -#define PRECISION 64 -#define EMAX 16384 -#define TYPE long double -#define SUFFIX l - -#define mpc_get_type mpc_get_ldc -#define mpc_set_type mpc_set_ldc -#define mpfr_set_type mpfr_set_ld - -gmp_randstate_t state; -mpz_t expz; /* global variable used in mpcheck_random */ -unsigned long seed = 0; -int verbose = 0; -mpfr_exp_t emin, emax; - -#include "mpcheck-common.c" - -#define FOO add -#define CFOO(x,y) (x+y) -#include "mpcheck-template3.c" - -#define FOO sub -#define CFOO(x,y) (x-y) -#include "mpcheck-template3.c" - -#define FOO mul -#define CFOO(x,y) (x*y) -#include "mpcheck-template3.c" - -#define FOO div -#define CFOO(x,y) (x/y) -#include "mpcheck-template3.c" - -#define FOO pow -#include "mpcheck-template3.c" - -#define FOO abs -#include "mpcheck-template2.c" - -#define FOO arg -#include "mpcheck-template2.c" - -#define FOO sqrt -#include "mpcheck-template1.c" - -#define FOO acos -#include "mpcheck-template1.c" - -#define FOO acosh -#include "mpcheck-template1.c" - -#define FOO asin -#include "mpcheck-template1.c" - -#define FOO asinh -#include "mpcheck-template1.c" - -#define FOO atan -#include "mpcheck-template1.c" - -#define FOO atanh -#include "mpcheck-template1.c" - -#define FOO cos -#include "mpcheck-template1.c" - -#define FOO cosh -#include "mpcheck-template1.c" - -#define FOO exp -#include "mpcheck-template1.c" - -#define FOO log -#include "mpcheck-template1.c" - -#define FOO log10 -#include "mpcheck-template1.c" - -#define FOO sin -#include "mpcheck-template1.c" - -#define FOO sinh -#include "mpcheck-template1.c" - -#define FOO tan -#include "mpcheck-template1.c" - -#define FOO tanh -#include "mpcheck-template1.c" - -int -main (int argc, char *argv[]) -{ - mpfr_prec_t p = PRECISION; /* precision of 'double' */ - unsigned long n = 1000000; /* default number of random tests per function */ - - while (argc >= 2 && argv[1][0] == '-') - { - if (argc >= 3 && strcmp (argv[1], "-p") == 0) - { - p = atoi (argv[2]); - argc -= 2; - argv += 2; - } - else if (argc >= 3 && strcmp (argv[1], "-seed") == 0) - { - seed = atoi (argv[2]); - argc -= 2; - argv += 2; - } - else if (argc >= 3 && strcmp (argv[1], "-num") == 0) - { - n = atoi (argv[2]); - argc -= 2; - argv += 2; - } - else if (strcmp (argv[1], "-v") == 0) - { - verbose ++; - argc --; - argv ++; - } - else if (strcmp (argv[1], "-check") == 0) - { - recheck = 1; - argc --; - argv ++; - } - else - { - fprintf (stderr, "Unknown option %s\n", argv[1]); - exit (1); - } - } - - /* set exponent range */ - emin = -EMAX - PRECISION + 4; /* should be -16444 */ - emax = EMAX; - mpfr_set_emin (emin); - mpfr_set_emax (emax); - - gmp_randinit_default (state); - mpz_init (expz); - -#ifdef __GNUC__ - printf ("GNU libc version: %s\n", gnu_get_libc_version ()); - printf ("GNU libc release: %s\n", gnu_get_libc_release ()); -#endif - - if (seed == 0) - seed = getpid (); - printf ("Using random seed %lu\n", seed); - - /* (complex,complex) -> complex */ - test_add (p, n); - test_sub (p, n); - test_mul (p, n); - test_div (p, n); - test_pow (p, n); - - /* complex -> real */ - test_abs (p, n); - test_arg (p, n); - - /* complex -> complex */ - test_sqrt (p, n); - test_acos (p, n); - test_acosh (p, n); - test_asin (p, n); - test_asinh (p, n); - test_atan (p, n); - test_atanh (p, n); - test_cos (p, n); - test_cosh (p, n); - test_exp (p, n); - test_log (p, n); - test_log10 (p, n); - test_sin (p, n); - test_sinh (p, n); - test_tan (p, n); - test_tanh (p, n); - - gmp_randclear (state); - mpz_clear (expz); - - report_maximal_errors (); - - return 0; -} diff --git a/tests/mpcheck-template1.c b/tests/mpcheck-template1.c deleted file mode 100644 index ee56e01..0000000 --- a/tests/mpcheck-template1.c +++ /dev/null @@ -1,155 +0,0 @@ -/* template for complex -> complex function, for example acos - -Copyright (C) 2020 INRIA - -This file is part of GNU MPC. - -GNU MPC 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 3 of the License, or (at your -option) any later version. - -GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . -*/ - -#define CAT1(X,Y) X ## Y -#define CAT2(X,Y) CAT1(X,Y) -#define STRINGIFY(x) #x -#define TOSTRING(x) STRINGIFY(x) - -#define FUN CAT2(test_,FOO) -#define CHECK CAT2(check_,FOO) -#define MPC_FOO CAT2(mpc_,FOO) -#define CFOO CAT2(CAT2(c,FOO),SUFFIX) -#define BAR TOSTRING(FOO) - -#define EXTRA 20 - -/* routine to double-check result from MPC with more precision and without - reduced exponent range */ -static void -CHECK (mpc_t x, mpc_t z) -{ - mpfr_exp_t saved_emin = mpfr_get_emin (); - mpfr_exp_t saved_emax = mpfr_get_emax (); - mpfr_prec_t p = mpfr_get_prec (mpc_realref (x)), pp = p + EXTRA; - mpc_t t; - int ok, inex_re, inex_im; - - /* enlarge exponent range to the largest possible */ - mpfr_set_emin (mpfr_get_emin_min ()); - mpfr_set_emax (mpfr_get_emax_max ()); - - mpc_init2 (t, pp); - MPC_FOO (t, x, MPC_RNDNN); - inex_re = mpfr_prec_round (mpc_realref (t), p, MPFR_RNDN); - inex_im = mpfr_prec_round (mpc_imagref (t), p, MPFR_RNDN); - /* restore exponent range */ - mpfr_set_emin (saved_emin); - mpfr_set_emax (saved_emax); - inex_re = mpfr_check_range (mpc_realref (t), inex_re, MPFR_RNDN); - inex_re = mpfr_subnormalize (mpc_realref (t), inex_re, MPFR_RNDN); - inex_im = mpfr_check_range (mpc_imagref (t), inex_im, MPFR_RNDN); - inex_im = mpfr_subnormalize (mpc_imagref (t), inex_im, MPFR_RNDN); - - /* check real parts agree */ - ok = mpfr_agree (mpc_realref (z), mpc_realref (t), inex_re) - && mpfr_agree (mpc_imagref (z), mpc_imagref (t), inex_im); - - if (!ok) - { - mpfr_printf ("Potential bug in mpc_%s for x=(%Re,%Re)\n", BAR, - mpc_realref (x), mpc_imagref (x)); - mpfr_printf (" mpc_%s to precision %lu gives (%Re,%Re)\n", BAR, p, - mpc_realref (z), mpc_imagref (z)); - mpfr_printf (" mpc_%s to precision %lu gives (%Re,%Re)\n", BAR, pp, - mpc_realref (t), mpc_imagref (t)); - exit (1); - } - - mpc_clear (t); -} - -static void -FUN (mpfr_prec_t p, unsigned long n) -{ - unsigned long i = 0; - mpc_t x, z, t; - TYPE complex xx, zz; - int inex, cmp; - unsigned long errors = 0, max_err_re = 0, max_err_im = 0; - /* allow reduced exponent range */ -#if defined(FOO_EMIN) || defined(FOO_EMAX) - mpfr_exp_t saved_emin = emin, saved_emax = emax; - emin = FOO_EMIN; - emax = FOO_EMAX; -#endif - - if (verbose > 0) - printf ("Testing function %s\n", BAR); - - gmp_randseed_ui (state, seed); - - mpc_init2 (x, p); - mpc_init2 (z, p); - mpc_init2 (t, p); - for (i = 0; i < n; i++) - { - mpcheck_random (x); - inex = MPC_FOO (z, x, MPC_RNDNN); - mpfr_subnormalize (mpc_realref (z), MPC_INEX_RE(inex), MPFR_RNDN); - mpfr_subnormalize (mpc_imagref (z), MPC_INEX_IM(inex), MPFR_RNDN); - xx = mpc_get_type (x, MPC_RNDNN); - zz = CFOO (xx); - mpc_set_type (t, zz, MPFR_RNDN); - cmp = mpcheck_mpc_cmp (t, z); - if (cmp != 0) - { - unsigned long err_re = ulp_error (mpc_realref (t), mpc_realref (z)); - unsigned long err_im = ulp_error (mpc_imagref (t), mpc_imagref (z)); - if (verbose > 0 && (err_re > max_err_re || err_im > max_err_im)) - { - mpfr_printf (" mpc_%s and c%s differ by (%lu,%lu) ulp(s)\n for x=(%Re,%Re)\n", - BAR, BAR, err_re, err_im, - mpc_realref (x), mpc_imagref (x)); - mpfr_printf (" mpc_%s gives (%Re,%Re)\n", BAR, - mpc_realref (z), mpc_imagref (z)); - mpfr_printf (" c%s gives (%Re,%Re)\n", BAR, - mpc_realref (t), mpc_imagref (t)); - } - if (recheck) - CHECK (x, z); - errors ++; - if (err_re > max_err_re) - { - max_err_re = err_re; - if (max_err_re > Max_err_re) - Max_err_re = max_err_re; - } - if (err_im > max_err_im) - { - max_err_im = err_im; - if (max_err_im > Max_err_im) - Max_err_im = max_err_im; - } - } - } - mpc_clear (x); - mpc_clear (z); - mpc_clear (t); - printf ("Errors for %s: %lu/%lu (max %lu,%lu) [seed %lu]\n", BAR, - errors, n, max_err_re, max_err_im, seed); -#if defined(FOO_EMIN) || defined(FOO_EMAX) - emin = saved_emin; - emax = saved_emax; -#endif -} - -#undef FOO -#undef CFOO diff --git a/tests/mpcheck-template2.c b/tests/mpcheck-template2.c deleted file mode 100644 index 7a2a8c0..0000000 --- a/tests/mpcheck-template2.c +++ /dev/null @@ -1,130 +0,0 @@ -/* template for complex -> real function, for example abs - -Copyright (C) 2020 INRIA - -This file is part of GNU MPC. - -GNU MPC 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 3 of the License, or (at your -option) any later version. - -GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . -*/ - -#define CAT1(X,Y) X ## Y -#define CAT2(X,Y) CAT1(X,Y) -#define STRINGIFY(x) #x -#define TOSTRING(x) STRINGIFY(x) - -#define FUN CAT2(test_,FOO) -#define CHECK CAT2(check_,FOO) -#define MPC_FOO CAT2(mpc_,FOO) -#define CFOO CAT2(CAT2(c,FOO),SUFFIX) -#define BAR TOSTRING(FOO) - -#define EXTRA 20 - -/* routine to double-check result from MPC with more precision and without - reduced exponent range */ -static void -CHECK (mpc_t x, mpfr_t z) -{ - mpfr_exp_t saved_emin = mpfr_get_emin (); - mpfr_exp_t saved_emax = mpfr_get_emax (); - mpfr_prec_t p = mpfr_get_prec (mpc_realref (x)), pp = p + EXTRA; - mpfr_t t; - int inex, ok; - - /* enlarge exponent range to the largest possible */ - mpfr_set_emin (mpfr_get_emin_min ()); - mpfr_set_emax (mpfr_get_emax_max ()); - - mpfr_init2 (t, pp); - MPC_FOO (t, x, MPFR_RNDN); - inex = mpfr_prec_round (t, p, MPFR_RNDN); - /* restore exponent range */ - mpfr_set_emin (saved_emin); - mpfr_set_emax (saved_emax); - inex = mpfr_check_range (t, inex, MPFR_RNDN); - inex = mpfr_subnormalize (t, inex, MPFR_RNDN); - - /* check real parts agree */ - ok = mpfr_agree (z, t, inex); - - if (!ok) - { - mpfr_printf ("Potential bug in mpc_%s for x=(%Re,%Re)\n", BAR, - mpc_realref (x), mpc_imagref (x)); - mpfr_printf (" mpc_%s to precision %lu gives %Re\n", BAR, p, z); - mpfr_printf (" mpc_%s to precision %lu gives %Re\n", BAR, pp, t); - exit (1); - } - - mpfr_clear (t); -} - -static void -FUN (mpfr_prec_t p, unsigned long n) -{ - unsigned long i = 0; - mpc_t x; - mpfr_t z, t; - TYPE complex xx; - TYPE zz; - int inex; - unsigned long errors = 0, max_err = 0; - - if (verbose > 0) - printf ("Testing function %s\n", BAR); - - gmp_randseed_ui (state, seed); - - mpc_init2 (x, p); - mpfr_init2 (z, p); - mpfr_init2 (t, p); - for (i = 0; i < n; i++) - { - mpcheck_random (x); - inex = MPC_FOO (z, x, MPC_RNDNN); - mpfr_subnormalize (z, inex, MPFR_RNDN); - xx = mpc_get_type (x, MPC_RNDNN); - zz = CFOO (xx); - mpfr_set_type (t, zz, MPFR_RNDN); - if (mpcheck_mpfr_cmp (t, z) != 0) - { - unsigned long err = ulp_error (t, z); - if (verbose > 0 && err > max_err) - { - mpfr_printf (" mpc_%s and c%s differ by %lu ulp(s)\n for x=(%Re,%Re)\n", - BAR, BAR, err, - mpc_realref (x), mpc_imagref (x)); - mpfr_printf (" mpc_%s gives %Re\n", BAR, z); - mpfr_printf (" c%s gives %Re\n", BAR, t); - } - if (recheck) - CHECK (x, z); - errors++; - if (err > max_err) - { - max_err = err; - if (max_err > Max_err) - Max_err = max_err; - } - } - } - mpc_clear (x); - mpfr_clear (z); - mpfr_clear (t); - printf ("Errors for %s: %lu/%lu (max %lu) [seed %lu]\n", - BAR, errors, n, max_err, seed); -} - -#undef FOO -#undef CFOO diff --git a/tests/mpcheck-template3.c b/tests/mpcheck-template3.c deleted file mode 100644 index e6718a8..0000000 --- a/tests/mpcheck-template3.c +++ /dev/null @@ -1,172 +0,0 @@ -/* template for (complex,complex) -> complex function, for example add - -Copyright (C) 2020 INRIA - -This file is part of GNU MPC. - -GNU MPC 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 3 of the License, or (at your -option) any later version. - -GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . -*/ - -#define CAT1(X,Y) X ## Y -#define CAT2(X,Y) CAT1(X,Y) -#define STRINGIFY(x) #x -#define TOSTRING(x) STRINGIFY(x) - -#define FUN CAT2(test_,FOO) -#define CHECK CAT2(check_,FOO) -#define MPC_FOO CAT2(mpc_,FOO) -#ifndef CFOO -#define CFOO CAT2(CAT2(c,FOO),SUFFIX) -#endif -#define BAR TOSTRING(FOO) - -#define EXTRA 20 - -/* routine to double-check result from MPC with more precision and without - reduced exponent range */ -static void -CHECK (mpc_t x, mpc_t y, mpc_t z) -{ - mpfr_exp_t saved_emin = mpfr_get_emin (); - mpfr_exp_t saved_emax = mpfr_get_emax (); - mpfr_prec_t p = mpfr_get_prec (mpc_realref (x)), pp = p + EXTRA; - mpc_t t; - int ok, inex_re, inex_im; - - /* enlarge exponent range to the largest possible */ - mpfr_set_emin (mpfr_get_emin_min ()); - mpfr_set_emax (mpfr_get_emax_max ()); - - mpc_init2 (t, pp); - MPC_FOO (t, x, y, MPC_RNDNN); - inex_re = mpfr_prec_round (mpc_realref (t), p, MPFR_RNDN); - inex_im = mpfr_prec_round (mpc_imagref (t), p, MPFR_RNDN); - /* restore exponent range */ - mpfr_set_emin (saved_emin); - mpfr_set_emax (saved_emax); - inex_re = mpfr_check_range (mpc_realref (t), inex_re, MPFR_RNDN); - inex_re = mpfr_subnormalize (mpc_realref (t), inex_re, MPFR_RNDN); - inex_im = mpfr_check_range (mpc_imagref (t), inex_im, MPFR_RNDN); - inex_im = mpfr_subnormalize (mpc_imagref (t), inex_im, MPFR_RNDN); - - /* check real parts agree */ - ok = mpfr_agree (mpc_realref (z), mpc_realref (t), inex_re) - && mpfr_agree (mpc_imagref (z), mpc_imagref (t), inex_im); - - if (!ok) - { - mpfr_printf ("Potential bug in mpc_%s for\n", BAR); - mpfr_printf (" x=(%Re,%Re)\n", mpc_realref (x), mpc_imagref (x)); - mpfr_printf (" y=(%Re,%Re)\n", mpc_realref (y), mpc_imagref (y)); - mpfr_printf (" mpc_%s to precision %lu gives (%Re,%Re)\n", BAR, p, - mpc_realref (z), mpc_imagref (z)); - mpfr_printf (" mpc_%s to precision %lu gives (%Re,%Re)\n", BAR, pp, - mpc_realref (t), mpc_imagref (t)); - exit (1); - } - - mpc_clear (t); -} - -static void -FUN (mpfr_prec_t p, unsigned long n) -{ - unsigned long i = 0; - mpc_t x, y, z, t; - TYPE complex xx, yy, zz; - int inex; - unsigned long errors = 0, max_err_re = 0, max_err_im = 0; - unsigned long special_errors = 0; - int cmp; - - if (verbose > 0) - printf ("Testing function %s\n", BAR); - - gmp_randseed_ui (state, seed); - - mpc_init2 (x, p); - mpc_init2 (y, p); - mpc_init2 (z, p); - mpc_init2 (t, p); - for (i = 0; i < n; i++) - { - mpcheck_random (x); - mpcheck_random (y); - inex = MPC_FOO (z, x, y, MPC_RNDNN); - mpfr_subnormalize (mpc_realref (z), MPC_INEX_RE(inex), MPFR_RNDN); - mpfr_subnormalize (mpc_imagref (z), MPC_INEX_IM(inex), MPFR_RNDN); - xx = mpc_get_type (x, MPC_RNDNN); - yy = mpc_get_dc (y, MPC_RNDNN); - zz = CFOO(xx, yy); - mpc_set_type (t, zz, MPFR_RNDN); - cmp = mpcheck_mpc_cmp (t, z); - if (cmp > 1) /* NaN or Inf */ - { - if (verbose > 1 || (verbose == 1 && special_errors == 0)) - { - printf (" mpc_%s and c%s differ\n", BAR, BAR); - mpfr_printf (" for x=(%Re,%Re)\n y=(%Re,%Re)\n", - mpc_realref (x), mpc_imagref (x), - mpc_realref (y), mpc_imagref (y)); - mpfr_printf (" mpc_%s gives (%Re,%Re)\n", BAR, - mpc_realref (z), mpc_imagref (z)); - mpfr_printf (" c%s gives (%Re,%Re)\n", BAR, - mpc_realref (t), mpc_imagref (t)); - special_errors ++; - } - errors ++; - } - else if (cmp != 0) - { - unsigned long err_re = ulp_error (mpc_realref (t), mpc_realref (z)); - unsigned long err_im = ulp_error (mpc_imagref (t), mpc_imagref (z)); - if (verbose > 0 && (err_re > max_err_re || err_im > max_err_im)) - { - printf (" mpc_%s and c%s differ by (%lu,%lu) ulp(s)\n", - BAR, BAR, err_re, err_im); - mpfr_printf (" for x=(%Re,%Re)\n y=(%Re,%Re)\n", - mpc_realref (x), mpc_imagref (x), - mpc_realref (y), mpc_imagref (y)); - mpfr_printf (" mpc_%s gives (%Re,%Re)\n", BAR, - mpc_realref (z), mpc_imagref (z)); - mpfr_printf (" c%s gives (%Re,%Re)\n", BAR, - mpc_realref (t), mpc_imagref (t)); - } - errors ++; - if (err_re > max_err_re) - { - max_err_re = err_re; - if (max_err_re > Max_err_re) - Max_err_re = max_err_re; - } - if (err_im > max_err_im) - { - max_err_im = err_im; - if (max_err_im > Max_err_im) - Max_err_im = max_err_im; - } - } - if (recheck && cmp != 0) - CHECK (x, y, z); - } - mpc_clear (x); - mpc_clear (y); - mpc_clear (z); - mpc_clear (t); - printf ("Errors for %s: %lu/%lu (max %lu,%lu) [seed %lu]\n", BAR, errors, n, - max_err_re, max_err_im, seed); -} - -#undef FOO -#undef CFOO diff --git a/tools/Makefile.am b/tools/Makefile.am index 5730d54..86257e1 100644 --- a/tools/Makefile.am +++ b/tools/Makefile.am @@ -1,6 +1,7 @@ ## tools/Makefile.am -- Process this file with automake to produce Makefile.in ## ## Copyright (C) 2014 CNRS +## Copyright (C) 2020 INRIA ## ## This file is part of GNU MPC. ## @@ -17,5 +18,5 @@ ## You should have received a copy of the GNU Lesser General Public License ## along with this program. If not, see http://www.gnu.org/licenses/ . -SUBDIRS = bench +SUBDIRS = bench mpcheck diff --git a/tools/mpcheck/Makefile.am b/tools/mpcheck/Makefile.am new file mode 100644 index 0000000..77ed6b0 --- /dev/null +++ b/tools/mpcheck/Makefile.am @@ -0,0 +1,31 @@ +## tools/mpcheck/Makefile.am -- Process this file with automake to produce Makefile.in +## +## Copyright (C) 2014 CNRS +## Copyright (C) 2020 INRIA +## +## This file is part of GNU MPC. +## +## GNU MPC 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 3 of the License, or (at your +## option) any later version. +## +## GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . + +AM_CPPFLAGS = -I$(top_srcdir)/src +AM_DEFAULT_SOURCE_EXT = .c + +LDADD = $(top_builddir)/src/libmpc.la -lm + +BINARIES = mpcheck-float mpcheck-double mpcheck-longdouble mpcheck-float128 +EXTRA_PROGRAMS = $(BINARIES) +CLEANFILES = $(BINARIES) + +mpcheck : $(BINARIES) + diff --git a/tools/mpcheck/README b/tools/mpcheck/README new file mode 100644 index 0000000..f35ff80 --- /dev/null +++ b/tools/mpcheck/README @@ -0,0 +1,3 @@ +The mpcheck tool provides tests to check the MPC library against the +C library. It is meant as a tool for developers. + diff --git a/tools/mpcheck/mpcheck-common.c b/tools/mpcheck/mpcheck-common.c new file mode 100644 index 0000000..0133f2a --- /dev/null +++ b/tools/mpcheck/mpcheck-common.c @@ -0,0 +1,193 @@ +/* shared routines for mpcheck-xxx + +Copyright (C) 2020 INRIA + +This file is part of GNU MPC. + +GNU MPC 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 3 of the License, or (at your +option) any later version. + +GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . +*/ + +#include +#include + +unsigned long Max_err_re = 0, Max_err_im = 0; +unsigned long Max_err = 0; /* for real-valued functions like abs */ +int recheck = 0; + +/* y is the reference value computed by MPC */ +static unsigned long +ulp_error (mpfr_t x, mpfr_t y) +{ + mpfr_t z; + mpfr_prec_t p = mpfr_get_prec (y); + unsigned long n, special = 0; + mpfr_exp_t e; + + if (mpfr_cmp (x, y) == 0) + return 0; + + assert (!mpfr_nan_p (x)); + assert (!mpfr_nan_p (y)); + + if (mpfr_inf_p (x) || mpfr_inf_p (y)) + { + if (mpfr_number_p (x) && mpfr_inf_p (y)) + { + /* only MPFR returns Inf: we replace y by the maximal + floating-point number */ + if (mpfr_sgn (y) > 0) /* +Inf */ + mpfr_nextbelow (y); + else + mpfr_nextabove (y); + special = 1; + } + else + { + /* only glibc returns Inf */ + if (mpfr_sgn (x) > 0) /* +Inf */ + mpfr_nextbelow (x); + else + mpfr_nextabove (x); + special = 2; + } + } + + mpfr_init2 (z, p); + mpfr_sub (z, x, y, MPFR_RNDN); + mpfr_abs (z, z, MPFR_RNDN); + /* divide by ulp(y) = 2^(EXP(y) - p) */ + e = mpfr_get_exp (y) - p; + /* take into account subnormals: smallest number is 2^(emin-1) */ + if (e < emin - 1) + e = emin - 1; + mpfr_div_2si (z, z, e, MPFR_RNDN); + n = mpfr_get_ui (z, MPFR_RNDZ); + mpfr_clear (z); + if (special == 1) /* restore original y */ + { + if (mpfr_sgn (y) > 0) + mpfr_nextabove (y); + else + mpfr_nextbelow (y); + } + if (special == 2) /* restore original x */ + { + if (mpfr_sgn (x) > 0) + mpfr_nextabove (x); + else + mpfr_nextbelow (x); + } + if (special == 0) + return n; + else + return (n < ULONG_MAX) ? n + 1 : n; +} + +static void +mpcheck_random (mpc_t z) +{ + unsigned long erange; + + mpc_urandom (z, state); + + /* set sign and exponent of the real part */ + mpz_set_ui (expz, 2 * (emax - emin + 1)); + mpz_urandomm (expz, state, expz); + erange = mpz_get_ui (expz); + /* bit 0 of erange gives the sign */ + if (erange & 1) + mpfr_neg (mpc_realref (z), mpc_realref (z), MPFR_RNDN); + mpfr_set_exp (mpc_realref (z), emin + (erange >> 1)); + /* make sure the numbers are representable in IEEE format */ + mpfr_subnormalize (mpc_realref (z), 0, MPFR_RNDN); + + /* set sign and exponent of the imaginary part */ + mpz_set_ui (expz, 2 * (emax - emin + 1)); + mpz_urandomm (expz, state, expz); + erange = mpz_get_ui (expz); + /* bit 0 of erange gives the sign */ + if (erange & 1) + mpfr_neg (mpc_imagref (z), mpc_imagref (z), MPFR_RNDN); + mpfr_set_exp (mpc_imagref (z), emin + (erange >> 1)); + /* make sure the numbers are representable in IEEE format */ + mpfr_subnormalize (mpc_imagref (z), 0, MPFR_RNDN); +} + +/* y is the reference value computed with MPC: + * return 0 if x and y are equal (including both NaN or Inf) + * return 1 if x and y disagree (neither NaN nor Inf) + * return 2 if x is NaN but not y + * return 4 if x and y are infinities of different signs + */ +static int +mpcheck_mpfr_cmp (mpfr_t x, mpfr_t y) +{ + if (mpfr_number_p (x) && mpfr_number_p (y)) + return mpfr_cmp (x, y); + /* now either x or y is NaN or Inf */ + if (mpfr_nan_p (x) || mpfr_nan_p (y)) + { + if (mpfr_nan_p (x) && mpfr_nan_p (y)) + return 0; + if (mpfr_nan_p (x)) + return 2; + fprintf (stderr, "Error, NaN and not NaN\n"); + exit (1); + } + /* now either x or y is Inf */ + if (mpfr_inf_p (x) && mpfr_inf_p (y)) + { + if (mpfr_sgn (x) * mpfr_sgn (y) > 0) + return 0; + return 4; + } + return 1; /* Inf and non-Inf: differ, will be treated by ulp_error */ +} + +/* y is the reference value computed with MPC */ +static int +mpcheck_mpc_cmp (mpc_t x, mpc_t y) +{ + return mpcheck_mpfr_cmp (mpc_realref (x), mpc_realref (y)) | + mpcheck_mpfr_cmp (mpc_imagref (x), mpc_imagref (y)); +} + +static void +report_maximal_errors (void) +{ + printf ("Maximal errors: %lu (real), %lu (imag), %lu (abs/arg)\n", + Max_err_re, Max_err_im, Max_err); +} + +/* return non-zero if both a and b are NaN, or if both are not NaN + and are equal */ +static int +mpfr_agree (mpfr_t a, mpfr_t b, int inex) +{ + if (mpfr_nan_p (a) || mpfr_nan_p (b)) + return mpfr_nan_p (a) && mpfr_nan_p (b); + else + /* in the case of even rounding, we can't tell for sure */ + return mpfr_equal_p (a, b) || inex == 2 || inex == -2; +} + +static int +cputime (void) +{ + struct rusage rus; + + getrusage (0, &rus); + return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000; +} + diff --git a/tools/mpcheck/mpcheck-double.c b/tools/mpcheck/mpcheck-double.c new file mode 100644 index 0000000..df216aa --- /dev/null +++ b/tools/mpcheck/mpcheck-double.c @@ -0,0 +1,242 @@ +/* mpcheck-double -- compare mpc functions against "double complex" + from the GNU libc implementation + +Copyright (C) 2020 INRIA + +This file is part of GNU MPC. + +GNU MPC 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 3 of the License, or (at your +option) any later version. + +GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . +*/ + +/* the GNU libc provides the following functions (as of 2.31), + with 'f' suffix for the float/binary32 version, with no suffix + for the double/binary64 version, with 'l' suffix for the long double + version, and with 'f128' suffix for the __float128 version: + + cabs casinh cexp csinh + cacos catan clog csqrt + cacosh catanh clog10 ctan + carg ccos cpow ctanh + casin ccosh csin +*/ + +#define _GNU_SOURCE /* for clog10 */ +#include +#include +#include +#include +#include +#include +#include +#ifdef __GNUC__ +#include +#endif +#include "mpc.h" + +#define PRECISION 53 +#define EMAX 1024 +#define TYPE double +#define SUFFIX + +#define mpc_get_type mpc_get_dc +#define mpc_set_type mpc_set_dc +#define mpfr_set_type mpfr_set_d + +gmp_randstate_t state; +mpz_t expz; /* global variable used in mpcheck_random */ +unsigned long seed = 0; +int verbose = 0; +mpfr_exp_t emin, emax; + +#include "mpcheck-common.c" + +#define FOO add +#define CFOO(x,y) (x+y) +#include "mpcheck-template3.c" + +#define FOO sub +#define CFOO(x,y) (x-y) +#include "mpcheck-template3.c" + +#define FOO mul +#define CFOO(x,y) (x*y) +#include "mpcheck-template3.c" + +#define FOO div +#define CFOO(x,y) (x/y) +#include "mpcheck-template3.c" + +#define FOO pow +#include "mpcheck-template3.c" + +#define FOO abs +#include "mpcheck-template2.c" + +#define FOO arg +#include "mpcheck-template2.c" + +#define FOO sqrt +#include "mpcheck-template1.c" + +#define FOO acos +#include "mpcheck-template1.c" + +#define FOO acosh +#include "mpcheck-template1.c" + +#define FOO asin +#include "mpcheck-template1.c" + +#define FOO asinh +#include "mpcheck-template1.c" + +#define FOO atan +#include "mpcheck-template1.c" + +#define FOO atanh +#include "mpcheck-template1.c" + +#define FOO cos +#include "mpcheck-template1.c" + +#define FOO cosh +#include "mpcheck-template1.c" + +#define FOO exp +#include "mpcheck-template1.c" + +#define FOO log +#include "mpcheck-template1.c" + +#define FOO log10 +#include "mpcheck-template1.c" + +#define FOO sin +#include "mpcheck-template1.c" + +#define FOO sinh +#include "mpcheck-template1.c" + +/* use reduced exponent range for tan and tanh */ +#define FOO_EMIN -8 +#define FOO_EMAX 8 + +#define FOO tan +#include "mpcheck-template1.c" + +#define FOO tanh +#include "mpcheck-template1.c" + +#undef FOO_EMIN +#undef FOO_EMAX + +int +main (int argc, char *argv[]) +{ + mpfr_prec_t p = PRECISION; /* precision of 'double' */ + unsigned long n = 1000000; /* default number of random tests per function */ + + while (argc >= 2 && argv[1][0] == '-') + { + if (argc >= 3 && strcmp (argv[1], "-p") == 0) + { + p = atoi (argv[2]); + argc -= 2; + argv += 2; + } + else if (argc >= 3 && strcmp (argv[1], "-seed") == 0) + { + seed = atoi (argv[2]); + argc -= 2; + argv += 2; + } + else if (argc >= 3 && strcmp (argv[1], "-num") == 0) + { + n = atoi (argv[2]); + argc -= 2; + argv += 2; + } + else if (strcmp (argv[1], "-v") == 0) + { + verbose ++; + argc --; + argv ++; + } + else if (strcmp (argv[1], "-check") == 0) + { + recheck = 1; + argc --; + argv ++; + } + else + { + fprintf (stderr, "Unknown option %s\n", argv[1]); + exit (1); + } + } + + /* set exponent range */ + emin = -EMAX - PRECISION + 4; /* should be -1073 */ + emax = EMAX; + mpfr_set_emin (emin); + mpfr_set_emax (emax); + + gmp_randinit_default (state); + mpz_init (expz); + +#ifdef __GNUC__ + printf ("GNU libc version: %s\n", gnu_get_libc_version ()); + printf ("GNU libc release: %s\n", gnu_get_libc_release ()); +#endif + + if (seed == 0) + seed = getpid (); + printf ("Using random seed %lu\n", seed); + + /* (complex,complex) -> complex */ + test_add (p, n); + test_sub (p, n); + test_mul (p, n); + test_div (p, n); + test_pow (p, n); + + /* complex -> real */ + test_abs (p, n); + test_arg (p, n); + + /* complex -> complex */ + test_sqrt (p, n); + test_acos (p, n); + test_acosh (p, n); + test_asin (p, n); + test_asinh (p, n); + test_atan (p, n); + test_atanh (p, n); + test_cos (p, n); + test_cosh (p, n); + test_exp (p, n); + test_log (p, n); + test_log10 (p, n); + test_sin (p, n); + test_sinh (p, n); + test_tan (p, n); + test_tanh (p, n); + + gmp_randclear (state); + mpz_clear (expz); + + report_maximal_errors (); + + return 0; +} diff --git a/tools/mpcheck/mpcheck-float.c b/tools/mpcheck/mpcheck-float.c new file mode 100644 index 0000000..9e4bdea --- /dev/null +++ b/tools/mpcheck/mpcheck-float.c @@ -0,0 +1,247 @@ +/* mpcheck-float -- compare mpc functions against "float complex" + from the GNU libc implementation + +Copyright (C) 2020 INRIA + +This file is part of GNU MPC. + +GNU MPC 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 3 of the License, or (at your +option) any later version. + +GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . +*/ + +/* the GNU libc provides the following functions (as of 2.31), + with 'f' suffix for the float/binary32 version, with no suffix + for the double/binary64 version, with 'l' suffix for the long double + version, and with 'f128' suffix for the __float128 version: + + cabs casinh cexp csinh + cacos catan clog csqrt + cacosh catanh clog10 ctan + carg ccos cpow ctanh + casin ccosh csin +*/ + +#define _GNU_SOURCE /* for clog10 */ +#include +#include +#include +#include +#include +#include +#include +#include "mpc.h" +#ifdef __GNUC__ +#include +#endif + +#define PRECISION 24 +#define EMAX 128 +#define TYPE float +#define SUFFIX f + +#define mpfr_set_type mpfr_set_flt + +static TYPE complex +mpc_get_type (mpc_t x, mpc_rnd_t rnd) +{ + /* there is no mpc_get_fltc function */ + return (TYPE complex) mpc_get_dc (x, rnd); +} + +static int +mpc_set_type (mpc_t x, TYPE complex y, mpc_rnd_t rnd) +{ + /* there is no mpc_set_fltc function */ + return mpc_set_dc (x, (double complex) y, rnd); +} + +gmp_randstate_t state; +mpz_t expz; /* global variable used in mpcheck_random */ +unsigned long seed = 0; +int verbose = 0; +mpfr_exp_t emin, emax; + +#include "mpcheck-common.c" + +#define FOO add +#define CFOO(x,y) (x+y) +#include "mpcheck-template3.c" + +#define FOO sub +#define CFOO(x,y) (x-y) +#include "mpcheck-template3.c" + +#define FOO mul +#define CFOO(x,y) (x*y) +#include "mpcheck-template3.c" + +#define FOO div +#define CFOO(x,y) (x/y) +#include "mpcheck-template3.c" + +#define FOO pow +#include "mpcheck-template3.c" + +#define FOO abs +#include "mpcheck-template2.c" + +#define FOO arg +#include "mpcheck-template2.c" + +#define FOO sqrt +#include "mpcheck-template1.c" + +#define FOO acos +#include "mpcheck-template1.c" + +#define FOO acosh +#include "mpcheck-template1.c" + +#define FOO asin +#include "mpcheck-template1.c" + +#define FOO asinh +#include "mpcheck-template1.c" + +#define FOO atan +#include "mpcheck-template1.c" + +#define FOO atanh +#include "mpcheck-template1.c" + +#define FOO cos +#include "mpcheck-template1.c" + +#define FOO cosh +#include "mpcheck-template1.c" + +#define FOO exp +#include "mpcheck-template1.c" + +#define FOO log +#include "mpcheck-template1.c" + +#define FOO log10 +#include "mpcheck-template1.c" + +#define FOO sin +#include "mpcheck-template1.c" + +#define FOO sinh +#include "mpcheck-template1.c" + +#define FOO tan +#include "mpcheck-template1.c" + +#define FOO tanh +#include "mpcheck-template1.c" + +int +main (int argc, char *argv[]) +{ + mpfr_prec_t p = PRECISION; /* precision of 'float' */ + unsigned long n = 1000000; /* default number of random tests per function */ + + while (argc >= 2 && argv[1][0] == '-') + { + if (argc >= 3 && strcmp (argv[1], "-p") == 0) + { + p = atoi (argv[2]); + argc -= 2; + argv += 2; + } + else if (argc >= 3 && strcmp (argv[1], "-seed") == 0) + { + seed = atoi (argv[2]); + argc -= 2; + argv += 2; + } + else if (argc >= 3 && strcmp (argv[1], "-num") == 0) + { + n = atoi (argv[2]); + argc -= 2; + argv += 2; + } + else if (strcmp (argv[1], "-v") == 0) + { + verbose ++; + argc --; + argv ++; + } + else if (strcmp (argv[1], "-check") == 0) + { + recheck = 1; + argc --; + argv ++; + } + else + { + fprintf (stderr, "Unknown option %s\n", argv[1]); + exit (1); + } + } + + /* set exponent range */ + emin = -EMAX - PRECISION + 4; /* should be -148 */ + emax = EMAX; + mpfr_set_emin (emin); + mpfr_set_emax (emax); + + gmp_randinit_default (state); + mpz_init (expz); + +#ifdef __GNUC__ + printf ("GNU libc version: %s\n", gnu_get_libc_version ()); + printf ("GNU libc release: %s\n", gnu_get_libc_release ()); +#endif + + if (seed == 0) + seed = getpid (); + printf ("Using random seed %lu\n", seed); + + /* (complex,complex) -> complex */ + test_add (p, n); + test_sub (p, n); + test_mul (p, n); + test_div (p, n); + test_pow (p, n); + + /* complex -> real */ + test_abs (p, n); + test_arg (p, n); + + /* complex -> complex */ + test_sqrt (p, n); + test_acos (p, n); + test_acosh (p, n); + test_asin (p, n); + test_asinh (p, n); + test_atan (p, n); + test_atanh (p, n); + test_cos (p, n); + test_cosh (p, n); + test_exp (p, n); + test_log (p, n); + test_log10 (p, n); + test_sin (p, n); + test_sinh (p, n); + test_tan (p, n); + test_tanh (p, n); + + gmp_randclear (state); + mpz_clear (expz); + + report_maximal_errors (); + + return 0; +} diff --git a/tools/mpcheck/mpcheck-float128.c b/tools/mpcheck/mpcheck-float128.c new file mode 100644 index 0000000..ba08efd --- /dev/null +++ b/tools/mpcheck/mpcheck-float128.c @@ -0,0 +1,252 @@ +/* mpcheck-float128 -- compare mpc functions against "__float128 complex" + from the GNU libc implementation + +Copyright (C) 2020 INRIA + +This file is part of GNU MPC. + +GNU MPC 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 3 of the License, or (at your +option) any later version. + +GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . +*/ + +/* the GNU libc provides the following functions (as of 2.31), + with 'f' suffix for the float/binary32 version, with no suffix + for the double/binary64 version, with 'l' suffix for the long double + version, and with 'f128' suffix for the __float128 version: + + cabs casinh cexp csinh + cacos catan clog csqrt + cacosh catanh clog10 ctan + carg ccos cpow ctanh + casin ccosh csin +*/ + +#define _GNU_SOURCE /* for clog10 */ +#include +#include +#include +#include +#include +#include +#include +#define MPFR_WANT_FLOAT128 +#include "mpc.h" +#ifdef __GNUC__ +#include +#endif + +#define PRECISION 113 +#define EMAX 16384 +#define TYPE _Float128 +#define SUFFIX f128 + +#define mpfr_set_type mpfr_set_float128 + +static TYPE complex +mpc_get_type (mpc_t z, mpc_rnd_t rnd) +{ + TYPE x, y; + /* there is no mpc_get_float128c function */ + x = mpfr_get_float128 (mpc_realref (z), MPC_RND_RE(rnd)); + y = mpfr_get_float128 (mpc_imagref (z), MPC_RND_IM(rnd)); + return x + I * y; +} + +static int +mpc_set_type (mpc_t x, TYPE complex y, mpc_rnd_t rnd) +{ + /* there is no mpc_set_float128c function */ + mpfr_set_float128 (mpc_realref (x), crealf128 (y), MPC_RND_RE(rnd)); + mpfr_set_float128 (mpc_imagref (x), cimagf128 (y), MPC_RND_IM(rnd)); +} + +gmp_randstate_t state; +mpz_t expz; /* global variable used in mpcheck_random */ +unsigned long seed = 0; +int verbose = 0; +mpfr_exp_t emin, emax; + +#include "mpcheck-common.c" + +#define FOO add +#define CFOO(x,y) (x+y) +#include "mpcheck-template3.c" + +#define FOO sub +#define CFOO(x,y) (x-y) +#include "mpcheck-template3.c" + +#define FOO mul +#define CFOO(x,y) (x*y) +#include "mpcheck-template3.c" + +#define FOO div +#define CFOO(x,y) (x/y) +#include "mpcheck-template3.c" + +#define FOO pow +#include "mpcheck-template3.c" + +#define FOO abs +#include "mpcheck-template2.c" + +#define FOO arg +#include "mpcheck-template2.c" + +#define FOO sqrt +#include "mpcheck-template1.c" + +#define FOO acos +#include "mpcheck-template1.c" + +#define FOO acosh +#include "mpcheck-template1.c" + +#define FOO asin +#include "mpcheck-template1.c" + +#define FOO asinh +#include "mpcheck-template1.c" + +#define FOO atan +#include "mpcheck-template1.c" + +#define FOO atanh +#include "mpcheck-template1.c" + +#define FOO cos +#include "mpcheck-template1.c" + +#define FOO cosh +#include "mpcheck-template1.c" + +#define FOO exp +#include "mpcheck-template1.c" + +#define FOO log +#include "mpcheck-template1.c" + +#define FOO log10 +#include "mpcheck-template1.c" + +#define FOO sin +#include "mpcheck-template1.c" + +#define FOO sinh +#include "mpcheck-template1.c" + +#define FOO tan +#include "mpcheck-template1.c" + +#define FOO tanh +#include "mpcheck-template1.c" + +int +main (int argc, char *argv[]) +{ + mpfr_prec_t p = PRECISION; /* precision of 'double' */ + unsigned long n = 1000000; /* default number of random tests per function */ + + while (argc >= 2 && argv[1][0] == '-') + { + if (argc >= 3 && strcmp (argv[1], "-p") == 0) + { + p = atoi (argv[2]); + argc -= 2; + argv += 2; + } + else if (argc >= 3 && strcmp (argv[1], "-seed") == 0) + { + seed = atoi (argv[2]); + argc -= 2; + argv += 2; + } + else if (argc >= 3 && strcmp (argv[1], "-num") == 0) + { + n = atoi (argv[2]); + argc -= 2; + argv += 2; + } + else if (strcmp (argv[1], "-v") == 0) + { + verbose ++; + argc --; + argv ++; + } + else if (strcmp (argv[1], "-check") == 0) + { + recheck = 1; + argc --; + argv ++; + } + else + { + fprintf (stderr, "Unknown option %s\n", argv[1]); + exit (1); + } + } + + /* set exponent range */ + emin = -EMAX - 64 + 4; /* should be -16444 like for long double */ + emax = EMAX; + mpfr_set_emin (emin); + mpfr_set_emax (emax); + + gmp_randinit_default (state); + mpz_init (expz); + +#ifdef __GNUC__ + printf ("GNU libc version: %s\n", gnu_get_libc_version ()); + printf ("GNU libc release: %s\n", gnu_get_libc_release ()); +#endif + + if (seed == 0) + seed = getpid (); + printf ("Using random seed %lu\n", seed); + + /* (complex,complex) -> complex */ + test_add (p, n); + test_sub (p, n); + test_mul (p, n); + test_div (p, n); + test_pow (p, n); + + /* complex -> real */ + test_abs (p, n); + test_arg (p, n); + + /* complex -> complex */ + test_sqrt (p, n); + test_acos (p, n); + test_acosh (p, n); + test_asin (p, n); + test_asinh (p, n); + test_atan (p, n); + test_atanh (p, n); + test_cos (p, n); + test_cosh (p, n); + test_exp (p, n); + test_log (p, n); + test_log10 (p, n); + test_sin (p, n); + test_sinh (p, n); + test_tan (p, n); + test_tanh (p, n); + + gmp_randclear (state); + mpz_clear (expz); + + report_maximal_errors (); + + return 0; +} diff --git a/tools/mpcheck/mpcheck-longdouble.c b/tools/mpcheck/mpcheck-longdouble.c new file mode 100644 index 0000000..6f7aefc --- /dev/null +++ b/tools/mpcheck/mpcheck-longdouble.c @@ -0,0 +1,235 @@ +/* mpcheck-longdouble -- compare mpc functions against "long double complex" + from the GNU libc implementation + +Copyright (C) 2020 INRIA + +This file is part of GNU MPC. + +GNU MPC 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 3 of the License, or (at your +option) any later version. + +GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . +*/ + +/* the GNU libc provides the following functions (as of 2.31), + with 'f' suffix for the float/binary32 version, with no suffix + for the double/binary64 version, with 'l' suffix for the long double + version, and with 'f128' suffix for the __float128 version: + + cabs casinh cexp csinh + cacos catan clog csqrt + cacosh catanh clog10 ctan + carg ccos cpow ctanh + casin ccosh csin +*/ + +#define _GNU_SOURCE /* for clog10 */ +#include +#include +#include +#include +#include +#include +#include +#include "mpc.h" +#ifdef __GNUC__ +#include +#endif + +#define PRECISION 64 +#define EMAX 16384 +#define TYPE long double +#define SUFFIX l + +#define mpc_get_type mpc_get_ldc +#define mpc_set_type mpc_set_ldc +#define mpfr_set_type mpfr_set_ld + +gmp_randstate_t state; +mpz_t expz; /* global variable used in mpcheck_random */ +unsigned long seed = 0; +int verbose = 0; +mpfr_exp_t emin, emax; + +#include "mpcheck-common.c" + +#define FOO add +#define CFOO(x,y) (x+y) +#include "mpcheck-template3.c" + +#define FOO sub +#define CFOO(x,y) (x-y) +#include "mpcheck-template3.c" + +#define FOO mul +#define CFOO(x,y) (x*y) +#include "mpcheck-template3.c" + +#define FOO div +#define CFOO(x,y) (x/y) +#include "mpcheck-template3.c" + +#define FOO pow +#include "mpcheck-template3.c" + +#define FOO abs +#include "mpcheck-template2.c" + +#define FOO arg +#include "mpcheck-template2.c" + +#define FOO sqrt +#include "mpcheck-template1.c" + +#define FOO acos +#include "mpcheck-template1.c" + +#define FOO acosh +#include "mpcheck-template1.c" + +#define FOO asin +#include "mpcheck-template1.c" + +#define FOO asinh +#include "mpcheck-template1.c" + +#define FOO atan +#include "mpcheck-template1.c" + +#define FOO atanh +#include "mpcheck-template1.c" + +#define FOO cos +#include "mpcheck-template1.c" + +#define FOO cosh +#include "mpcheck-template1.c" + +#define FOO exp +#include "mpcheck-template1.c" + +#define FOO log +#include "mpcheck-template1.c" + +#define FOO log10 +#include "mpcheck-template1.c" + +#define FOO sin +#include "mpcheck-template1.c" + +#define FOO sinh +#include "mpcheck-template1.c" + +#define FOO tan +#include "mpcheck-template1.c" + +#define FOO tanh +#include "mpcheck-template1.c" + +int +main (int argc, char *argv[]) +{ + mpfr_prec_t p = PRECISION; /* precision of 'double' */ + unsigned long n = 1000000; /* default number of random tests per function */ + + while (argc >= 2 && argv[1][0] == '-') + { + if (argc >= 3 && strcmp (argv[1], "-p") == 0) + { + p = atoi (argv[2]); + argc -= 2; + argv += 2; + } + else if (argc >= 3 && strcmp (argv[1], "-seed") == 0) + { + seed = atoi (argv[2]); + argc -= 2; + argv += 2; + } + else if (argc >= 3 && strcmp (argv[1], "-num") == 0) + { + n = atoi (argv[2]); + argc -= 2; + argv += 2; + } + else if (strcmp (argv[1], "-v") == 0) + { + verbose ++; + argc --; + argv ++; + } + else if (strcmp (argv[1], "-check") == 0) + { + recheck = 1; + argc --; + argv ++; + } + else + { + fprintf (stderr, "Unknown option %s\n", argv[1]); + exit (1); + } + } + + /* set exponent range */ + emin = -EMAX - PRECISION + 4; /* should be -16444 */ + emax = EMAX; + mpfr_set_emin (emin); + mpfr_set_emax (emax); + + gmp_randinit_default (state); + mpz_init (expz); + +#ifdef __GNUC__ + printf ("GNU libc version: %s\n", gnu_get_libc_version ()); + printf ("GNU libc release: %s\n", gnu_get_libc_release ()); +#endif + + if (seed == 0) + seed = getpid (); + printf ("Using random seed %lu\n", seed); + + /* (complex,complex) -> complex */ + test_add (p, n); + test_sub (p, n); + test_mul (p, n); + test_div (p, n); + test_pow (p, n); + + /* complex -> real */ + test_abs (p, n); + test_arg (p, n); + + /* complex -> complex */ + test_sqrt (p, n); + test_acos (p, n); + test_acosh (p, n); + test_asin (p, n); + test_asinh (p, n); + test_atan (p, n); + test_atanh (p, n); + test_cos (p, n); + test_cosh (p, n); + test_exp (p, n); + test_log (p, n); + test_log10 (p, n); + test_sin (p, n); + test_sinh (p, n); + test_tan (p, n); + test_tanh (p, n); + + gmp_randclear (state); + mpz_clear (expz); + + report_maximal_errors (); + + return 0; +} diff --git a/tools/mpcheck/mpcheck-template1.c b/tools/mpcheck/mpcheck-template1.c new file mode 100644 index 0000000..ee56e01 --- /dev/null +++ b/tools/mpcheck/mpcheck-template1.c @@ -0,0 +1,155 @@ +/* template for complex -> complex function, for example acos + +Copyright (C) 2020 INRIA + +This file is part of GNU MPC. + +GNU MPC 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 3 of the License, or (at your +option) any later version. + +GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . +*/ + +#define CAT1(X,Y) X ## Y +#define CAT2(X,Y) CAT1(X,Y) +#define STRINGIFY(x) #x +#define TOSTRING(x) STRINGIFY(x) + +#define FUN CAT2(test_,FOO) +#define CHECK CAT2(check_,FOO) +#define MPC_FOO CAT2(mpc_,FOO) +#define CFOO CAT2(CAT2(c,FOO),SUFFIX) +#define BAR TOSTRING(FOO) + +#define EXTRA 20 + +/* routine to double-check result from MPC with more precision and without + reduced exponent range */ +static void +CHECK (mpc_t x, mpc_t z) +{ + mpfr_exp_t saved_emin = mpfr_get_emin (); + mpfr_exp_t saved_emax = mpfr_get_emax (); + mpfr_prec_t p = mpfr_get_prec (mpc_realref (x)), pp = p + EXTRA; + mpc_t t; + int ok, inex_re, inex_im; + + /* enlarge exponent range to the largest possible */ + mpfr_set_emin (mpfr_get_emin_min ()); + mpfr_set_emax (mpfr_get_emax_max ()); + + mpc_init2 (t, pp); + MPC_FOO (t, x, MPC_RNDNN); + inex_re = mpfr_prec_round (mpc_realref (t), p, MPFR_RNDN); + inex_im = mpfr_prec_round (mpc_imagref (t), p, MPFR_RNDN); + /* restore exponent range */ + mpfr_set_emin (saved_emin); + mpfr_set_emax (saved_emax); + inex_re = mpfr_check_range (mpc_realref (t), inex_re, MPFR_RNDN); + inex_re = mpfr_subnormalize (mpc_realref (t), inex_re, MPFR_RNDN); + inex_im = mpfr_check_range (mpc_imagref (t), inex_im, MPFR_RNDN); + inex_im = mpfr_subnormalize (mpc_imagref (t), inex_im, MPFR_RNDN); + + /* check real parts agree */ + ok = mpfr_agree (mpc_realref (z), mpc_realref (t), inex_re) + && mpfr_agree (mpc_imagref (z), mpc_imagref (t), inex_im); + + if (!ok) + { + mpfr_printf ("Potential bug in mpc_%s for x=(%Re,%Re)\n", BAR, + mpc_realref (x), mpc_imagref (x)); + mpfr_printf (" mpc_%s to precision %lu gives (%Re,%Re)\n", BAR, p, + mpc_realref (z), mpc_imagref (z)); + mpfr_printf (" mpc_%s to precision %lu gives (%Re,%Re)\n", BAR, pp, + mpc_realref (t), mpc_imagref (t)); + exit (1); + } + + mpc_clear (t); +} + +static void +FUN (mpfr_prec_t p, unsigned long n) +{ + unsigned long i = 0; + mpc_t x, z, t; + TYPE complex xx, zz; + int inex, cmp; + unsigned long errors = 0, max_err_re = 0, max_err_im = 0; + /* allow reduced exponent range */ +#if defined(FOO_EMIN) || defined(FOO_EMAX) + mpfr_exp_t saved_emin = emin, saved_emax = emax; + emin = FOO_EMIN; + emax = FOO_EMAX; +#endif + + if (verbose > 0) + printf ("Testing function %s\n", BAR); + + gmp_randseed_ui (state, seed); + + mpc_init2 (x, p); + mpc_init2 (z, p); + mpc_init2 (t, p); + for (i = 0; i < n; i++) + { + mpcheck_random (x); + inex = MPC_FOO (z, x, MPC_RNDNN); + mpfr_subnormalize (mpc_realref (z), MPC_INEX_RE(inex), MPFR_RNDN); + mpfr_subnormalize (mpc_imagref (z), MPC_INEX_IM(inex), MPFR_RNDN); + xx = mpc_get_type (x, MPC_RNDNN); + zz = CFOO (xx); + mpc_set_type (t, zz, MPFR_RNDN); + cmp = mpcheck_mpc_cmp (t, z); + if (cmp != 0) + { + unsigned long err_re = ulp_error (mpc_realref (t), mpc_realref (z)); + unsigned long err_im = ulp_error (mpc_imagref (t), mpc_imagref (z)); + if (verbose > 0 && (err_re > max_err_re || err_im > max_err_im)) + { + mpfr_printf (" mpc_%s and c%s differ by (%lu,%lu) ulp(s)\n for x=(%Re,%Re)\n", + BAR, BAR, err_re, err_im, + mpc_realref (x), mpc_imagref (x)); + mpfr_printf (" mpc_%s gives (%Re,%Re)\n", BAR, + mpc_realref (z), mpc_imagref (z)); + mpfr_printf (" c%s gives (%Re,%Re)\n", BAR, + mpc_realref (t), mpc_imagref (t)); + } + if (recheck) + CHECK (x, z); + errors ++; + if (err_re > max_err_re) + { + max_err_re = err_re; + if (max_err_re > Max_err_re) + Max_err_re = max_err_re; + } + if (err_im > max_err_im) + { + max_err_im = err_im; + if (max_err_im > Max_err_im) + Max_err_im = max_err_im; + } + } + } + mpc_clear (x); + mpc_clear (z); + mpc_clear (t); + printf ("Errors for %s: %lu/%lu (max %lu,%lu) [seed %lu]\n", BAR, + errors, n, max_err_re, max_err_im, seed); +#if defined(FOO_EMIN) || defined(FOO_EMAX) + emin = saved_emin; + emax = saved_emax; +#endif +} + +#undef FOO +#undef CFOO diff --git a/tools/mpcheck/mpcheck-template2.c b/tools/mpcheck/mpcheck-template2.c new file mode 100644 index 0000000..7a2a8c0 --- /dev/null +++ b/tools/mpcheck/mpcheck-template2.c @@ -0,0 +1,130 @@ +/* template for complex -> real function, for example abs + +Copyright (C) 2020 INRIA + +This file is part of GNU MPC. + +GNU MPC 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 3 of the License, or (at your +option) any later version. + +GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . +*/ + +#define CAT1(X,Y) X ## Y +#define CAT2(X,Y) CAT1(X,Y) +#define STRINGIFY(x) #x +#define TOSTRING(x) STRINGIFY(x) + +#define FUN CAT2(test_,FOO) +#define CHECK CAT2(check_,FOO) +#define MPC_FOO CAT2(mpc_,FOO) +#define CFOO CAT2(CAT2(c,FOO),SUFFIX) +#define BAR TOSTRING(FOO) + +#define EXTRA 20 + +/* routine to double-check result from MPC with more precision and without + reduced exponent range */ +static void +CHECK (mpc_t x, mpfr_t z) +{ + mpfr_exp_t saved_emin = mpfr_get_emin (); + mpfr_exp_t saved_emax = mpfr_get_emax (); + mpfr_prec_t p = mpfr_get_prec (mpc_realref (x)), pp = p + EXTRA; + mpfr_t t; + int inex, ok; + + /* enlarge exponent range to the largest possible */ + mpfr_set_emin (mpfr_get_emin_min ()); + mpfr_set_emax (mpfr_get_emax_max ()); + + mpfr_init2 (t, pp); + MPC_FOO (t, x, MPFR_RNDN); + inex = mpfr_prec_round (t, p, MPFR_RNDN); + /* restore exponent range */ + mpfr_set_emin (saved_emin); + mpfr_set_emax (saved_emax); + inex = mpfr_check_range (t, inex, MPFR_RNDN); + inex = mpfr_subnormalize (t, inex, MPFR_RNDN); + + /* check real parts agree */ + ok = mpfr_agree (z, t, inex); + + if (!ok) + { + mpfr_printf ("Potential bug in mpc_%s for x=(%Re,%Re)\n", BAR, + mpc_realref (x), mpc_imagref (x)); + mpfr_printf (" mpc_%s to precision %lu gives %Re\n", BAR, p, z); + mpfr_printf (" mpc_%s to precision %lu gives %Re\n", BAR, pp, t); + exit (1); + } + + mpfr_clear (t); +} + +static void +FUN (mpfr_prec_t p, unsigned long n) +{ + unsigned long i = 0; + mpc_t x; + mpfr_t z, t; + TYPE complex xx; + TYPE zz; + int inex; + unsigned long errors = 0, max_err = 0; + + if (verbose > 0) + printf ("Testing function %s\n", BAR); + + gmp_randseed_ui (state, seed); + + mpc_init2 (x, p); + mpfr_init2 (z, p); + mpfr_init2 (t, p); + for (i = 0; i < n; i++) + { + mpcheck_random (x); + inex = MPC_FOO (z, x, MPC_RNDNN); + mpfr_subnormalize (z, inex, MPFR_RNDN); + xx = mpc_get_type (x, MPC_RNDNN); + zz = CFOO (xx); + mpfr_set_type (t, zz, MPFR_RNDN); + if (mpcheck_mpfr_cmp (t, z) != 0) + { + unsigned long err = ulp_error (t, z); + if (verbose > 0 && err > max_err) + { + mpfr_printf (" mpc_%s and c%s differ by %lu ulp(s)\n for x=(%Re,%Re)\n", + BAR, BAR, err, + mpc_realref (x), mpc_imagref (x)); + mpfr_printf (" mpc_%s gives %Re\n", BAR, z); + mpfr_printf (" c%s gives %Re\n", BAR, t); + } + if (recheck) + CHECK (x, z); + errors++; + if (err > max_err) + { + max_err = err; + if (max_err > Max_err) + Max_err = max_err; + } + } + } + mpc_clear (x); + mpfr_clear (z); + mpfr_clear (t); + printf ("Errors for %s: %lu/%lu (max %lu) [seed %lu]\n", + BAR, errors, n, max_err, seed); +} + +#undef FOO +#undef CFOO diff --git a/tools/mpcheck/mpcheck-template3.c b/tools/mpcheck/mpcheck-template3.c new file mode 100644 index 0000000..e6718a8 --- /dev/null +++ b/tools/mpcheck/mpcheck-template3.c @@ -0,0 +1,172 @@ +/* template for (complex,complex) -> complex function, for example add + +Copyright (C) 2020 INRIA + +This file is part of GNU MPC. + +GNU MPC 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 3 of the License, or (at your +option) any later version. + +GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . +*/ + +#define CAT1(X,Y) X ## Y +#define CAT2(X,Y) CAT1(X,Y) +#define STRINGIFY(x) #x +#define TOSTRING(x) STRINGIFY(x) + +#define FUN CAT2(test_,FOO) +#define CHECK CAT2(check_,FOO) +#define MPC_FOO CAT2(mpc_,FOO) +#ifndef CFOO +#define CFOO CAT2(CAT2(c,FOO),SUFFIX) +#endif +#define BAR TOSTRING(FOO) + +#define EXTRA 20 + +/* routine to double-check result from MPC with more precision and without + reduced exponent range */ +static void +CHECK (mpc_t x, mpc_t y, mpc_t z) +{ + mpfr_exp_t saved_emin = mpfr_get_emin (); + mpfr_exp_t saved_emax = mpfr_get_emax (); + mpfr_prec_t p = mpfr_get_prec (mpc_realref (x)), pp = p + EXTRA; + mpc_t t; + int ok, inex_re, inex_im; + + /* enlarge exponent range to the largest possible */ + mpfr_set_emin (mpfr_get_emin_min ()); + mpfr_set_emax (mpfr_get_emax_max ()); + + mpc_init2 (t, pp); + MPC_FOO (t, x, y, MPC_RNDNN); + inex_re = mpfr_prec_round (mpc_realref (t), p, MPFR_RNDN); + inex_im = mpfr_prec_round (mpc_imagref (t), p, MPFR_RNDN); + /* restore exponent range */ + mpfr_set_emin (saved_emin); + mpfr_set_emax (saved_emax); + inex_re = mpfr_check_range (mpc_realref (t), inex_re, MPFR_RNDN); + inex_re = mpfr_subnormalize (mpc_realref (t), inex_re, MPFR_RNDN); + inex_im = mpfr_check_range (mpc_imagref (t), inex_im, MPFR_RNDN); + inex_im = mpfr_subnormalize (mpc_imagref (t), inex_im, MPFR_RNDN); + + /* check real parts agree */ + ok = mpfr_agree (mpc_realref (z), mpc_realref (t), inex_re) + && mpfr_agree (mpc_imagref (z), mpc_imagref (t), inex_im); + + if (!ok) + { + mpfr_printf ("Potential bug in mpc_%s for\n", BAR); + mpfr_printf (" x=(%Re,%Re)\n", mpc_realref (x), mpc_imagref (x)); + mpfr_printf (" y=(%Re,%Re)\n", mpc_realref (y), mpc_imagref (y)); + mpfr_printf (" mpc_%s to precision %lu gives (%Re,%Re)\n", BAR, p, + mpc_realref (z), mpc_imagref (z)); + mpfr_printf (" mpc_%s to precision %lu gives (%Re,%Re)\n", BAR, pp, + mpc_realref (t), mpc_imagref (t)); + exit (1); + } + + mpc_clear (t); +} + +static void +FUN (mpfr_prec_t p, unsigned long n) +{ + unsigned long i = 0; + mpc_t x, y, z, t; + TYPE complex xx, yy, zz; + int inex; + unsigned long errors = 0, max_err_re = 0, max_err_im = 0; + unsigned long special_errors = 0; + int cmp; + + if (verbose > 0) + printf ("Testing function %s\n", BAR); + + gmp_randseed_ui (state, seed); + + mpc_init2 (x, p); + mpc_init2 (y, p); + mpc_init2 (z, p); + mpc_init2 (t, p); + for (i = 0; i < n; i++) + { + mpcheck_random (x); + mpcheck_random (y); + inex = MPC_FOO (z, x, y, MPC_RNDNN); + mpfr_subnormalize (mpc_realref (z), MPC_INEX_RE(inex), MPFR_RNDN); + mpfr_subnormalize (mpc_imagref (z), MPC_INEX_IM(inex), MPFR_RNDN); + xx = mpc_get_type (x, MPC_RNDNN); + yy = mpc_get_dc (y, MPC_RNDNN); + zz = CFOO(xx, yy); + mpc_set_type (t, zz, MPFR_RNDN); + cmp = mpcheck_mpc_cmp (t, z); + if (cmp > 1) /* NaN or Inf */ + { + if (verbose > 1 || (verbose == 1 && special_errors == 0)) + { + printf (" mpc_%s and c%s differ\n", BAR, BAR); + mpfr_printf (" for x=(%Re,%Re)\n y=(%Re,%Re)\n", + mpc_realref (x), mpc_imagref (x), + mpc_realref (y), mpc_imagref (y)); + mpfr_printf (" mpc_%s gives (%Re,%Re)\n", BAR, + mpc_realref (z), mpc_imagref (z)); + mpfr_printf (" c%s gives (%Re,%Re)\n", BAR, + mpc_realref (t), mpc_imagref (t)); + special_errors ++; + } + errors ++; + } + else if (cmp != 0) + { + unsigned long err_re = ulp_error (mpc_realref (t), mpc_realref (z)); + unsigned long err_im = ulp_error (mpc_imagref (t), mpc_imagref (z)); + if (verbose > 0 && (err_re > max_err_re || err_im > max_err_im)) + { + printf (" mpc_%s and c%s differ by (%lu,%lu) ulp(s)\n", + BAR, BAR, err_re, err_im); + mpfr_printf (" for x=(%Re,%Re)\n y=(%Re,%Re)\n", + mpc_realref (x), mpc_imagref (x), + mpc_realref (y), mpc_imagref (y)); + mpfr_printf (" mpc_%s gives (%Re,%Re)\n", BAR, + mpc_realref (z), mpc_imagref (z)); + mpfr_printf (" c%s gives (%Re,%Re)\n", BAR, + mpc_realref (t), mpc_imagref (t)); + } + errors ++; + if (err_re > max_err_re) + { + max_err_re = err_re; + if (max_err_re > Max_err_re) + Max_err_re = max_err_re; + } + if (err_im > max_err_im) + { + max_err_im = err_im; + if (max_err_im > Max_err_im) + Max_err_im = max_err_im; + } + } + if (recheck && cmp != 0) + CHECK (x, y, z); + } + mpc_clear (x); + mpc_clear (y); + mpc_clear (z); + mpc_clear (t); + printf ("Errors for %s: %lu/%lu (max %lu,%lu) [seed %lu]\n", BAR, errors, n, + max_err_re, max_err_im, seed); +} + +#undef FOO +#undef CFOO -- cgit v1.2.1