diff options
Diffstat (limited to 'tools')
| -rw-r--r-- | tools/Makefile.am | 3 | ||||
| -rw-r--r-- | tools/mpcheck/Makefile.am | 31 | ||||
| -rw-r--r-- | tools/mpcheck/README | 3 | ||||
| -rw-r--r-- | tools/mpcheck/mpcheck-common.c | 193 | ||||
| -rw-r--r-- | tools/mpcheck/mpcheck-double.c | 242 | ||||
| -rw-r--r-- | tools/mpcheck/mpcheck-float.c | 247 | ||||
| -rw-r--r-- | tools/mpcheck/mpcheck-float128.c | 252 | ||||
| -rw-r--r-- | tools/mpcheck/mpcheck-longdouble.c | 235 | ||||
| -rw-r--r-- | tools/mpcheck/mpcheck-template1.c | 155 | ||||
| -rw-r--r-- | tools/mpcheck/mpcheck-template2.c | 130 | ||||
| -rw-r--r-- | tools/mpcheck/mpcheck-template3.c | 172 |
11 files changed, 1662 insertions, 1 deletions
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 <sys/types.h> +#include <sys/resource.h> + +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 <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <complex.h> +#include <sys/types.h> +#include <unistd.h> +#include <assert.h> +#ifdef __GNUC__ +#include <gnu/libc-version.h> +#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 <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <complex.h> +#include <sys/types.h> +#include <unistd.h> +#include <assert.h> +#include "mpc.h" +#ifdef __GNUC__ +#include <gnu/libc-version.h> +#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 <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <complex.h> +#include <sys/types.h> +#include <unistd.h> +#include <assert.h> +#define MPFR_WANT_FLOAT128 +#include "mpc.h" +#ifdef __GNUC__ +#include <gnu/libc-version.h> +#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 <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <complex.h> +#include <sys/types.h> +#include <unistd.h> +#include <assert.h> +#include "mpc.h" +#ifdef __GNUC__ +#include <gnu/libc-version.h> +#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 |