diff options
author | Paul Zimmermann <Paul.Zimmermann@inria.fr> | 2020-02-06 18:35:19 +0100 |
---|---|---|
committer | Paul Zimmermann <Paul.Zimmermann@inria.fr> | 2020-02-06 18:35:19 +0100 |
commit | 0fce03e7e78fe6bf789e3d372fdc6dbd27bf63a6 (patch) | |
tree | 8bd293d0f99c83ca2931616eeef66348b1dea866 | |
parent | fbcde25676931fbdbf707d03fd980c1936a00b02 (diff) | |
download | mpc-git-0fce03e7e78fe6bf789e3d372fdc6dbd27bf63a6.tar.gz |
enlarge internal exponent range for several functions + mpcheck cleanup
-rw-r--r-- | src/acos.c | 12 | ||||
-rw-r--r-- | src/asin.c | 22 | ||||
-rw-r--r-- | src/atan.c | 24 | ||||
-rw-r--r-- | src/pow.c | 15 | ||||
-rw-r--r-- | src/sin_cos.c | 33 | ||||
-rw-r--r-- | src/tan.c | 9 | ||||
-rw-r--r-- | tests/Makefile.am | 2 | ||||
-rw-r--r-- | tests/mpcheck-double.c | 49 | ||||
-rw-r--r-- | tests/mpcheck-float.c | 38 | ||||
-rw-r--r-- | tests/mpcheck-float128.c | 38 | ||||
-rw-r--r-- | tests/mpcheck-longdouble.c | 38 | ||||
-rw-r--r-- | tests/mpcheck-template1.c | 20 | ||||
-rw-r--r-- | tests/mpcheck-template2.c | 9 | ||||
-rw-r--r-- | tests/mpcheck-template3.c | 29 | ||||
-rw-r--r-- | tests/pow_ui.dat | 4 | ||||
-rw-r--r-- | tests/tdiv.c | 1 |
16 files changed, 216 insertions, 127 deletions
@@ -31,6 +31,7 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_exp_t e1, e2; mpfr_rnd_t rnd_im; mpc_rnd_t rnd1; + mpfr_exp_t saved_emin, saved_emax; inex_re = 0; inex_im = 0; @@ -170,6 +171,11 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) return MPC_INEX (inex_re, inex_im); } + saved_emin = mpfr_get_emin (); + saved_emax = mpfr_get_emax (); + mpfr_set_emin (mpfr_get_emin_min ()); + mpfr_set_emax (mpfr_get_emax_max ()); + /* regular complex argument: acos(z) = Pi/2 - asin(z) */ p_re = mpfr_get_prec (mpc_realref(rop)); p_im = mpfr_get_prec (mpc_imagref(rop)); @@ -225,5 +231,11 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpc_clear (z1); mpfr_clear (pi_over_2); + /* restore the exponent range, and check the range of results */ + mpfr_set_emin (saved_emin); + mpfr_set_emax (saved_emax); + inex_re = mpfr_check_range (mpc_realref (rop), inex_re, MPC_RND_RE (rnd)); + inex_im = mpfr_check_range (mpc_imagref (rop), inex_im, MPC_RND_IM (rnd)); + return MPC_INEX(inex_re, inex_im); } @@ -86,7 +86,8 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_prec_t p, p_re, p_im; mpfr_rnd_t rnd_re, rnd_im; mpc_t z1; - int inex, loop = 0; + int inex, inex_re, inex_im, loop = 0; + mpfr_exp_t saved_emin, saved_emax; /* special values */ if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) @@ -112,7 +113,6 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op))) { - int inex_re; if (mpfr_inf_p (mpc_realref (op))) { int inf_im = mpfr_inf_p (mpc_imagref (op)); @@ -137,8 +137,6 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* pure real argument */ if (mpfr_zero_p (mpc_imagref (op))) { - int inex_re; - int inex_im; int s_im; s_im = mpfr_signbit (mpc_imagref (op)); @@ -186,7 +184,6 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* pure imaginary argument */ if (mpfr_zero_p (mpc_realref (op))) { - int inex_im; int s; s = mpfr_signbit (mpc_realref (op)); mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN); @@ -197,6 +194,11 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) return MPC_INEX (0, inex_im); } + saved_emin = mpfr_get_emin (); + saved_emax = mpfr_get_emax (); + mpfr_set_emin (mpfr_get_emin_min ()); + mpfr_set_emax (mpfr_get_emax_max ()); + /* regular complex: asin(z) = -i*log(i*z+sqrt(1-z^2)) */ p_re = mpfr_get_prec (mpc_realref(rop)); p_im = mpfr_get_prec (mpc_imagref(rop)); @@ -286,5 +288,13 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) inex = mpc_set (rop, z1, rnd); mpc_clear (z1); - return inex; + /* restore the exponent range, and check the range of results */ + mpfr_set_emin (saved_emin); + mpfr_set_emax (saved_emax); + inex_re = mpfr_check_range (mpc_realref (rop), MPC_INEX_RE (inex), + MPC_RND_RE (rnd)); + inex_im = mpfr_check_range (mpc_imagref (rop), MPC_INEX_IM (inex), + MPC_RND_IM (rnd)); + + return MPC_INEX (inex_re, inex_im); } @@ -45,11 +45,9 @@ set_pi_over_2 (mpfr_ptr rop, int s, mpfr_rnd_t rnd) int mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { - int s_re; - int s_im; - int inex_re; - int inex_im; - int inex; + int s_re, s_im; + int inex_re, inex_im, inex; + mpfr_exp_t saved_emin, saved_emax; inex_re = 0; inex_im = 0; @@ -216,6 +214,11 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) return MPC_INEX (inex_re, inex_im); } + saved_emin = mpfr_get_emin (); + saved_emax = mpfr_get_emax (); + mpfr_set_emin (mpfr_get_emin_min ()); + mpfr_set_emax (mpfr_get_emax_max ()); + /* regular number argument */ { mpfr_t a, b, x, y; @@ -391,6 +394,15 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) inex = mpc_set_fr_fr (rop, x, y, rnd); mpfr_clears (a, b, x, y, (mpfr_ptr) 0); - return inex; + + /* restore the exponent range, and check the range of results */ + mpfr_set_emin (saved_emin); + mpfr_set_emax (saved_emax); + inex_re = mpfr_check_range (mpc_realref (rop), MPC_INEX_RE (inex), + MPC_RND_RE (rnd)); + inex_im = mpfr_check_range (mpc_imagref (rop), MPC_INEX_IM (inex), + MPC_RND_IM (rnd)); + + return MPC_INEX (inex_re, inex_im); } } @@ -486,6 +486,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) mpfr_prec_t p, pr, pi, maxprec; int saved_underflow, saved_overflow; int inex_re, inex_im; + mpfr_exp_t saved_emin, saved_emax; /* save the underflow or overflow flags from MPFR */ saved_underflow = mpfr_underflow_p (); @@ -645,6 +646,11 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) z_real = 1; } + saved_emin = mpfr_get_emin (); + saved_emax = mpfr_get_emax (); + mpfr_set_emin (mpfr_get_emin_min ()); + mpfr_set_emax (mpfr_get_emax_max ()); + pr = mpfr_get_prec (mpc_realref(z)); pi = mpfr_get_prec (mpc_imagref(z)); p = (pr > pi) ? pr : pi; @@ -837,6 +843,15 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) if (saved_overflow) mpfr_set_overflow (); + /* restore the exponent range, and check the range of results */ + mpfr_set_emin (saved_emin); + mpfr_set_emax (saved_emax); + inex_re = mpfr_check_range (mpc_realref (z), MPC_INEX_RE(ret), + MPC_RND_RE (rnd)); + inex_im = mpfr_check_range (mpc_imagref (z), MPC_INEX_IM(ret), + MPC_RND_IM (rnd)); + ret = MPC_INEX(inex_re, inex_im); + end: return ret; } diff --git a/src/sin_cos.c b/src/sin_cos.c index da0dd17..a1be8a6 100644 --- a/src/sin_cos.c +++ b/src/sin_cos.c @@ -343,6 +343,12 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_prec_t prec; int ok; int inex_re, inex_im, inex_sin, inex_cos, loop = 0; + mpfr_exp_t saved_emin, saved_emax; + + saved_emin = mpfr_get_emin (); + saved_emax = mpfr_get_emax (); + mpfr_set_emin (mpfr_get_emin_min ()); + mpfr_set_emax (mpfr_get_emax_max ()); prec = 2; if (rop_sin != NULL) @@ -376,7 +382,6 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, do { loop ++; - ok = 1; prec += (loop <= 2) ? mpc_ceil_log2 (prec) + 5 : prec / 2; mpfr_set_prec (s, prec); @@ -389,6 +394,8 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_sin_cos (s, c, mpc_realref(op), MPFR_RNDN); mpfr_sinh_cosh (sh, ch, mpc_imagref(op), MPFR_RNDN); + ok = 1; + if (rop_sin != NULL) { /* real part of sine */ mpfr_mul (sch, s, ch, MPFR_RNDN); @@ -458,6 +465,30 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_clear (sch); mpfr_clear (csh); + /* restore the exponent range, and check the range of results */ + mpfr_set_emin (saved_emin); + mpfr_set_emax (saved_emax); + if (rop_sin != NULL) + { + inex_re = mpfr_check_range (mpc_realref (rop_sin), + MPC_INEX_RE (inex_sin), + MPC_RND_RE (rnd_sin)); + inex_im = mpfr_check_range (mpc_imagref (rop_sin), + MPC_INEX_IM (inex_sin), + MPC_RND_IM (rnd_sin)); + inex_sin = MPC_INEX (inex_re, inex_im); + } + if (rop_cos != NULL) + { + inex_re = mpfr_check_range (mpc_realref (rop_cos), + MPC_INEX_RE (inex_cos), + MPC_RND_RE (rnd_cos)); + inex_im = mpfr_check_range (mpc_imagref (rop_cos), + MPC_INEX_IM (inex_cos), + MPC_RND_IM (rnd_cos)); + inex_cos = MPC_INEX (inex_re, inex_im); + } + return (MPC_INEX12 (inex_sin, inex_cos)); } } @@ -29,7 +29,7 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_prec_t prec; mpfr_exp_t err; int ok = 0; - int inex; + int inex, inex_re, inex_im; /* special values */ if (!mpc_fin_p (op)) @@ -80,7 +80,6 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* tan(+Inf +i*Inf) = +/-0 +i */ { const int sign_re = mpfr_signbit (mpc_realref (op)); - int inex_im; mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd)); mpfr_setsign (mpc_realref (rop), mpc_realref (rop), sign_re, MPFR_RNDN); @@ -106,7 +105,6 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { mpfr_t c; mpfr_t s; - int inex_im; mpfr_init (c); mpfr_init (s); @@ -132,8 +130,6 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* tan(-0 -i*y) = -0 +i*tanh(y), when y is finite. */ /* tan(+0 +i*y) = +0 +i*tanh(y), when y is finite. */ { - int inex_im; - mpfr_set (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); inex_im = mpfr_tanh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); @@ -144,8 +140,6 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* tan(x -i*0) = tan(x) -i*0, when x is finite. */ /* tan(x +i*0) = tan(x) +i*0, when x is finite. */ { - int inex_re; - inex_re = mpfr_tan (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); mpfr_set (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); @@ -206,7 +200,6 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) Im(op) was large, in which case the result is sign(tan(Re(op)))*0 + sign(Im(op))*I, where sign(tan(Re(op))) = sign(Re(x))*sign(Re(y)). */ - int inex_re, inex_im; mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN); if (mpfr_sgn (mpc_realref (x)) * mpfr_sgn (mpc_realref (y)) < 0) { diff --git a/tests/Makefile.am b/tests/Makefile.am index 37f1496..6f6ec61 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -67,7 +67,7 @@ DATA_SETS = abs.dat acos.dat acosh.dat add.dat add_fr.dat arg.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-float128.c mpcheck-common.c \ mpcheck-template1.c mpcheck-template2.c mpcheck-template3.c LOG_COMPILER = $(VALGRIND) diff --git a/tests/mpcheck-double.c b/tests/mpcheck-double.c index 7b3e88e..8b93bad 100644 --- a/tests/mpcheck-double.c +++ b/tests/mpcheck-double.c @@ -35,10 +35,13 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include <stdlib.h> #include <string.h> #include <complex.h> -#include "mpc-tests.h" +#include <sys/types.h> +#include <unistd.h> +#include <assert.h> #ifdef __GNUC__ #include <gnu/libc-version.h> #endif +#include "mpc-tests.h" #define PRECISION 53 #define EMAX 1024 @@ -50,28 +53,12 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #define mpfr_set_type mpfr_set_d gmp_randstate_t state; -unsigned long seed = 1; +mpz_t expz; /* global variable used in mpcheck_random */ +unsigned long seed = 0; int verbose = 0; +mpfr_exp_t emin, emax; -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; - - if (mpfr_cmp (x, y) == 0) - return 0; - - 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) */ - mpfr_div_2si (z, z, mpfr_get_exp (y) - p, MPFR_RNDN); - n = mpfr_get_ui (z, MPFR_RNDZ); - mpfr_clear (z); - return n; -} +#include "mpcheck-common.c" #define FOO add #define CFOO(x,y) (x+y) @@ -140,12 +127,19 @@ ulp_error (mpfr_t x, mpfr_t y) #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[]) { @@ -185,16 +179,22 @@ main (int argc, char *argv[]) } } - /* set exponent range for 'double' */ - mpfr_set_emin (-EMAX - PRECISION + 4); /* should be -1073 */ - mpfr_set_emax (EMAX); + /* 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 */ @@ -227,6 +227,7 @@ main (int argc, char *argv[]) test_tanh (p, n); gmp_randclear (state); + mpz_clear (expz); return 0; } diff --git a/tests/mpcheck-float.c b/tests/mpcheck-float.c index e056c86..0adb6b6 100644 --- a/tests/mpcheck-float.c +++ b/tests/mpcheck-float.c @@ -35,6 +35,9 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include <stdlib.h> #include <string.h> #include <complex.h> +#include <sys/types.h> +#include <unistd.h> +#include <assert.h> #include "mpc-tests.h" #ifdef __GNUC__ #include <gnu/libc-version.h> @@ -62,28 +65,12 @@ mpc_set_type (mpc_t x, TYPE complex y, mpc_rnd_t rnd) } gmp_randstate_t state; -unsigned long seed = 1; +mpz_t expz; /* global variable used in mpcheck_random */ +unsigned long seed = 0; int verbose = 0; +mpfr_exp_t emin, emax; -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; - - if (mpfr_cmp (x, y) == 0) - return 0; - - 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) */ - mpfr_div_2si (z, z, mpfr_get_exp (y) - p, MPFR_RNDN); - n = mpfr_get_ui (z, MPFR_RNDZ); - mpfr_clear (z); - return n; -} +#include "mpcheck-common.c" #define FOO add #define CFOO(x,y) (x+y) @@ -198,15 +185,21 @@ main (int argc, char *argv[]) } /* set exponent range */ - mpfr_set_emin (-EMAX - PRECISION + 4); /* should be -148 */ - mpfr_set_emax (EMAX); + 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 */ @@ -239,6 +232,7 @@ main (int argc, char *argv[]) test_tanh (p, n); gmp_randclear (state); + mpz_clear (expz); return 0; } diff --git a/tests/mpcheck-float128.c b/tests/mpcheck-float128.c index 43f5ddf..cc929f3 100644 --- a/tests/mpcheck-float128.c +++ b/tests/mpcheck-float128.c @@ -35,6 +35,9 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #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-tests.h" #ifdef __GNUC__ @@ -67,28 +70,12 @@ mpc_set_type (mpc_t x, TYPE complex y, mpc_rnd_t rnd) } gmp_randstate_t state; -unsigned long seed = 1; +mpz_t expz; /* global variable used in mpcheck_random */ +unsigned long seed = 0; int verbose = 0; +mpfr_exp_t emin, emax; -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; - - if (mpfr_cmp (x, y) == 0) - return 0; - - 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) */ - mpfr_div_2si (z, z, mpfr_get_exp (y) - p, MPFR_RNDN); - n = mpfr_get_ui (z, MPFR_RNDZ); - mpfr_clear (z); - return n; -} +#include "mpcheck-common.c" #define FOO add #define CFOO(x,y) (x+y) @@ -203,15 +190,21 @@ main (int argc, char *argv[]) } /* set exponent range */ - mpfr_set_emin (-EMAX - 64 + 4); /* should be -16444 like for long double */ - mpfr_set_emax (EMAX); + 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 */ @@ -244,6 +237,7 @@ main (int argc, char *argv[]) test_tanh (p, n); gmp_randclear (state); + mpz_clear (expz); return 0; } diff --git a/tests/mpcheck-longdouble.c b/tests/mpcheck-longdouble.c index f991c2e..5f9895e 100644 --- a/tests/mpcheck-longdouble.c +++ b/tests/mpcheck-longdouble.c @@ -35,6 +35,9 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include <stdlib.h> #include <string.h> #include <complex.h> +#include <sys/types.h> +#include <unistd.h> +#include <assert.h> #include "mpc-tests.h" #ifdef __GNUC__ #include <gnu/libc-version.h> @@ -50,28 +53,12 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #define mpfr_set_type mpfr_set_ld gmp_randstate_t state; -unsigned long seed = 1; +mpz_t expz; /* global variable used in mpcheck_random */ +unsigned long seed = 0; int verbose = 0; +mpfr_exp_t emin, emax; -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; - - if (mpfr_cmp (x, y) == 0) - return 0; - - 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) */ - mpfr_div_2si (z, z, mpfr_get_exp (y) - p, MPFR_RNDN); - n = mpfr_get_ui (z, MPFR_RNDZ); - mpfr_clear (z); - return n; -} +#include "mpcheck-common.c" #define FOO add #define CFOO(x,y) (x+y) @@ -186,15 +173,21 @@ main (int argc, char *argv[]) } /* set exponent range */ - mpfr_set_emin (-EMAX - PRECISION + 4); /* should be -16444 */ - mpfr_set_emax (EMAX); + 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 */ @@ -227,6 +220,7 @@ main (int argc, char *argv[]) test_tanh (p, n); gmp_randclear (state); + mpz_clear (expz); return 0; } diff --git a/tests/mpcheck-template1.c b/tests/mpcheck-template1.c index 742f7d3..734f472 100644 --- a/tests/mpcheck-template1.c +++ b/tests/mpcheck-template1.c @@ -36,6 +36,14 @@ FUN (mpfr_prec_t p, unsigned long n) TYPE complex xx, zz; int inex; 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 + + printf ("Testing function %s\n", BAR); gmp_randseed_ui (state, seed); @@ -44,14 +52,14 @@ FUN (mpfr_prec_t p, unsigned long n) mpc_init2 (t, p); for (i = 0; i < n; i++) { - mpc_urandom (x, state); + 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); - if (mpc_cmp (z, t) != 0) + if (mpcheck_mpc_cmp (t, z) != 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)); @@ -75,8 +83,12 @@ FUN (mpfr_prec_t p, unsigned long n) mpc_clear (x); mpc_clear (z); mpc_clear (t); - printf ("Number of errors for %s: %lu/%lu (max %lu,%lu)\n", BAR, - errors, n, max_err_re, max_err_im); + 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 diff --git a/tests/mpcheck-template2.c b/tests/mpcheck-template2.c index d3513bf..2f94bfb 100644 --- a/tests/mpcheck-template2.c +++ b/tests/mpcheck-template2.c @@ -39,6 +39,8 @@ FUN (mpfr_prec_t p, unsigned long n) int inex; unsigned long errors = 0, max_err = 0; + printf ("Testing function %s\n", BAR); + gmp_randseed_ui (state, seed); mpc_init2 (x, p); @@ -46,13 +48,13 @@ FUN (mpfr_prec_t p, unsigned long n) mpfr_init2 (t, p); for (i = 0; i < n; i++) { - mpc_urandom (x, state); + 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 (mpfr_cmp (z, t) != 0) + if (mpcheck_mpfr_cmp (t, z) != 0) { unsigned long err = ulp_error (t, z); if (verbose > 0 && err > max_err) @@ -71,7 +73,8 @@ FUN (mpfr_prec_t p, unsigned long n) mpc_clear (x); mpfr_clear (z); mpfr_clear (t); - printf ("Number of errors for %s: %lu/%lu (max %lu)\n", BAR, errors, n, max_err); + printf ("Errors for %s: %lu/%lu (max %lu) [seed %lu]\n", + BAR, errors, n, max_err, seed); } #undef FOO diff --git a/tests/mpcheck-template3.c b/tests/mpcheck-template3.c index b42cb92..5d754c1 100644 --- a/tests/mpcheck-template3.c +++ b/tests/mpcheck-template3.c @@ -38,6 +38,9 @@ FUN (mpfr_prec_t p, unsigned long n) TYPE complex xx, yy, zz; int inex; unsigned long errors = 0, max_err_re = 0, max_err_im = 0; + int cmp; + + printf ("Testing function %s\n", BAR); gmp_randseed_ui (state, seed); @@ -47,8 +50,8 @@ FUN (mpfr_prec_t p, unsigned long n) mpc_init2 (t, p); for (i = 0; i < n; i++) { - mpc_urandom (x, state); - mpc_urandom (y, state); + 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); @@ -56,7 +59,23 @@ FUN (mpfr_prec_t p, unsigned long n) yy = mpc_get_dc (y, MPC_RNDNN); zz = CFOO(xx, yy); mpc_set_type (t, zz, MPFR_RNDN); - if (mpc_cmp (z, t) != 0) + cmp = mpcheck_mpc_cmp (t, z); + if (cmp > 1) + { + if (verbose > 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)); + } + 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)); @@ -83,8 +102,8 @@ FUN (mpfr_prec_t p, unsigned long n) mpc_clear (y); mpc_clear (z); mpc_clear (t); - printf ("Number of errors for %s: %lu/%lu (max %lu,%lu)\n", BAR, errors, n, - max_err_re, max_err_im); + printf ("Errors for %s: %lu/%lu (max %lu,%lu) [seed %lu]\n", BAR, errors, n, + max_err_re, max_err_im, seed); } #undef FOO diff --git a/tests/pow_ui.dat b/tests/pow_ui.dat index d448a68..b41df7e 100644 --- a/tests/pow_ui.dat +++ b/tests/pow_ui.dat @@ -85,8 +85,8 @@ 0 0 53 1 53 0 53 +1 53 +0 +2 N N 0 0 53 1 53 0 53 +1 53 -0 +2 N N -# overflow -? ? 53 +inf 53 +inf 53 1e100000000 53 1e100000000 1000000000 N N +# overflow: (a+a*I)^n is real for n divisible by 4 +? 0 53 +inf 53 0 53 1e100000000 53 1e100000000 1000000000 N N # underflow ? ? 53 0 53 0 53 1e-100000000 53 1e-100000000 1000000000 N N # cannot round after one loop diff --git a/tests/tdiv.c b/tests/tdiv.c index 66d32ea..6ca3115 100644 --- a/tests/tdiv.c +++ b/tests/tdiv.c @@ -48,7 +48,6 @@ bug20200206 (void) mpfr_set_d (mpc_imagref (x), -4.9188093228194944e-89, MPFR_RNDN); mpfr_set_d (mpc_realref (y), -3.6801500191882962e-14, MPFR_RNDN); mpfr_set_d (mpc_imagref (y), 4.5420247980297260e-145, MPFR_RNDN); - printf ("emin=%ld emax=%ld\n", mpfr_get_emin (), mpfr_get_emax ()); mpc_div (z, x, y, MPC_RNDNN); /* quotient is 1.44844440684571e-202 + 1.33657848108714e-75*I, where both the real and imaginary parts fit in the exponent range */ |