diff options
| -rw-r--r-- | NEWS | 3 | ||||
| -rw-r--r-- | src/pow_ui.c | 15 | ||||
| -rw-r--r-- | src/sqr.c | 1 | ||||
| -rw-r--r-- | tests/Makefile.am | 4 | ||||
| -rw-r--r-- | tests/mpc-tests.h | 2 | ||||
| -rw-r--r-- | tests/pow_ui.dat | 93 | ||||
| -rw-r--r-- | tests/read_data.c | 78 | ||||
| -rw-r--r-- | tests/tpow_ui.c | 6 |
8 files changed, 188 insertions, 14 deletions
@@ -1,3 +1,6 @@ +Changes in version 0.8.2: + - Speed-up of mpc_pow_ui through binary exponentiation + Changes in version 0.8.1: - Bug fixes: - acosh, asinh, atanh: swap of precisions between real and imaginary parts diff --git a/src/pow_ui.c b/src/pow_ui.c index ddcb9df..876614a 100644 --- a/src/pow_ui.c +++ b/src/pow_ui.c @@ -47,13 +47,20 @@ mpc_pow_ui (mpc_ptr z, mpc_srcptr x, unsigned long y, mpc_rnd_t rnd) mp_exp_t diff; int has3; /* non-zero if y has '11' in its binary representation */ - if (y == 1) + mp_prec_t exp_r = mpfr_get_exp (MPC_RE (x)), + exp_i = mpfr_get_exp (MPC_IM (x)); + if (!mpc_fin_p (x) || mpfr_zero_p (MPC_IM(x)) || y == 0 + || MPC_MAX (exp_r, exp_i) > mpfr_get_emax () / y + /* heuristic for overflow */ + || MPC_MAX (-exp_r, -exp_i) > (-mpfr_get_emin ()) / y + /* heuristic for underflow */ + ) + /* let mpc_pow deal with special cases */ + return mpc_pow_ui_naive (z, x, y, rnd); + else if (y == 1) return mpc_set (z, x, rnd); else if (y == 2) return mpc_sqr (z, x, rnd); - else if (!mpc_fin_p (x) || mpfr_zero_p (MPC_IM(x)) || y == 0) - /* let mpc_pow deal with special cases */ - return mpc_pow_ui_naive (z, x, y, rnd); for (l = 0, u = y, has3 = u&3; u > 3; l ++, u >>= 1, has3 |= u&3); /* l>0 is the number of bits of y, minus 2, thus y has bits: @@ -208,7 +208,6 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { mpfr_set_ui (MPC_RE (rop), 0, rnd_re); inex_re = 1; - } else /* round down or away from zero */ inex_re = -1; diff --git a/tests/Makefile.am b/tests/Makefile.am index 6907381..29f81a5 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -4,10 +4,10 @@ AM_CPPFLAGS = -I$(top_srcdir)/src LDADD = libmpc-tests.la $(top_builddir)/src/libmpc.la LOADLIBES=$(DEFS) -I$(top_srcdir)/src $(CPPFLAGS) $(CFLAGS) $(top_builddir)/src/.libs/libmpc.a $(LIBS) -check_PROGRAMS = tabs tacos tacosh tadd tadd_fr tadd_ui targ tasin tasinh \ +check_PROGRAMS = tpow_ui tabs tacos tacosh tadd tadd_fr tadd_ui targ tasin tasinh \ tatan tatanh tconj tcos tcosh tdiv tdiv_2exp tdiv_fr tdiv_ui texp tfr_div \ tfr_sub timag tio_str tlog tmul tmul_2exp tmul_fr tmul_i tmul_si \ -tmul_ui tneg tnorm tpow tpow_ld tpow_d tpow_fr tpow_si tpow_ui tpow_z tprec \ +tmul_ui tneg tnorm tpow tpow_ld tpow_d tpow_fr tpow_si tpow_z tprec \ tproj treal treimref tset tsin tsinh tsqr tsqrt tstrtoc tsub tsub_fr tsub_ui \ ttan ttanh tui_div tui_ui_sub tget_version diff --git a/tests/mpc-tests.h b/tests/mpc-tests.h index 24eef36..5c89db5 100644 --- a/tests/mpc-tests.h +++ b/tests/mpc-tests.h @@ -152,7 +152,7 @@ typedef struct with random numbers: - with precision ranging from prec_min to prec_max with an increment of step, - - with exponent between -prec_max and prec_max. + - with exponent between -exp_max and exp_max. It also checks parameter reuse (it is assumed here that either two mpc_t variables are equal or they are different, in the sense that the real part of diff --git a/tests/pow_ui.dat b/tests/pow_ui.dat new file mode 100644 index 0000000..be8c369 --- /dev/null +++ b/tests/pow_ui.dat @@ -0,0 +1,93 @@ +# Data file for mpc_pow_ui. +# +# Copyright (C) 2010 Andreas Enge +# +# This file is part of the MPC Library. +# +# The MPC Library is free software; you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation; either version 2.1 of the License, or (at your +# option) any later version. +# +# The MPC Library is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +# License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with the MPC Library; see the file COPYING.LIB. If not, write to +# the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +# MA 02111-1307, USA. +# +# The line format respects the parameter order in function prototype as +# follow: +# +# INEX_RE INEX_IM PREC_ROP_RE ROP_RE PREC_ROP_IM ROP_IM PREC_OP1_RE OP1_RE PREC_OP1_IM OP1_IM OP2 RND_RE RND_IM +# +# For further details, see add_fr.dat. + +# special cases, copied from pow.dat +0 0 53 +1 53 0 53 nan 53 +0 +0 N N +0 0 53 nan 53 nan 53 nan 53 +0 +1 N N +0 0 53 inf 53 nan 53 +inf 53 +0 +1 N N +0 0 53 +inf 53 nan 53 +inf 53 +1 +1 N N +0 0 53 +inf 53 nan 53 +inf 53 -1 +1 N N +0 0 53 +inf 53 nan 53 -inf 53 +0 +1 N N +0 0 53 +inf 53 nan 53 -inf 53 +1 +1 N N +0 0 53 +inf 53 nan 53 -inf 53 -1 +1 N N + +0 0 53 nan 53 nan 53 +0 53 +0 +0 N N +0 0 53 0 53 0 53 +0 53 +0 +1 N N + +0 0 53 +1 53 +0 53 +0 53 +1 +0 N N +0 0 53 +1 53 -0 53 +0 53 +1 +0 N D +0 0 53 +1 53 +0 53 -0 53 +1 +0 N N +0 0 53 +1 53 +0 53 -1 53 +0 +0 N N +0 0 53 +1 53 -0 53 -1 53 -0 +0 N N +0 0 53 +1 53 -0 53 -0 53 -1 +0 N N +0 0 53 +1 53 -0 53 +0 53 -1 +0 N N + +0 0 53 +1 53 +0 53 +inf 53 +2 +0 N N +0 0 53 +1 53 +0 53 +inf 53 -0 +0 N N +0 0 53 +1 53 +0 53 +2 53 +inf +0 N N +0 0 53 +1 53 +0 53 +2 53 +0 +0 N N +0 0 53 +1 53 +0 53 +0 53 +2 +0 N N +0 0 53 +1 53 +0 53 +0 53 +inf +0 N N +0 0 53 +1 53 +0 53 -0 53 +2 +0 N N +0 0 53 +1 53 +0 53 -0 53 +inf +0 N N +0 0 53 +1 53 +0 53 -5 53 +inf +0 N N +0 0 53 +1 53 +0 53 -2 53 +0 +0 N N +0 0 53 +1 53 +0 53 -inf 53 +0 +0 N N +0 0 53 +1 53 +0 53 -inf 53 +3 +0 N N +0 0 53 +1 53 -0 53 +0.5 53 -0.5 +0 N N +0 0 53 +1 53 -0 53 +0.5 53 +0 +0 N N +0 0 53 +1 53 -0 53 +0.5 53 -0 +0 N N +0 0 53 +1 53 -0 53 -0.5 53 -0 +0 N N +0 0 53 +1 53 -0 53 +0 53 -0.5 +0 N N +0 0 53 +1 53 -0 53 -0 53 -0.5 +0 N N +0 0 53 +9 53 +0 53 +3 53 +0 +2 N N +0 0 53 +1 53 +0 53 +1 53 +0 +4 N N +0 0 53 +1 53 -0 53 +1 53 -0 +4 N N +0 0 53 0.25 53 -0 53 +0.5 53 -0 +2 N N + +0 0 53 1 53 0 53 +2 53 -1 +0 N N +0 0 53 1 53 0 53 -2 53 -1 +0 N N +0 0 53 1 53 0 53 -2 53 -0 +0 N N +0 0 53 1 53 0 53 +0.5 53 +0.5 +0 N N +0 0 53 1 53 0 53 -0.5 53 +0.5 +0 N N +0 0 53 1 53 0 53 -0.5 53 +0 +0 N N +0 0 53 1 53 0 53 +0 53 +0.5 +0 N N +0 0 53 1 53 0 53 -0 53 +0.5 +0 N N +0 0 53 1 53 0 53 -0 53 -4 +0 N N +0 0 53 1 53 0 53 +0 53 -4 +0 N N +0 0 53 1 53 0 53 -1 53 -0 +0 N N +0 0 53 1 53 0 53 -1 53 +0 +0 N N + +0 0 53 4 53 0 53 +2 53 -0 +2 N N +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 +# underflow +? ? 53 0 53 0 53 1e-100000000 53 1e-100000000 1000000000 N N diff --git a/tests/read_data.c b/tests/read_data.c index cca7dfc..4129027 100644 --- a/tests/read_data.c +++ b/tests/read_data.c @@ -1,6 +1,6 @@ /* Read data file and check function. -Copyright (C) 2008, 2009 Andreas Enge, Philippe Th\'eveny +Copyright (C) 2008, 2009, 2010 Andreas Enge, Philippe Th\'eveny This file is part of the MPC Library. @@ -282,7 +282,7 @@ read_int (FILE *fp, int *nread, const char *name) if (nextchar == EOF) { - printf ("Error: Unexpected EOF when reading mpfr precision " + printf ("Error: Unexpected EOF when reading int " "in file '%s' line %lu\n", pathname, line_number); exit (1); @@ -299,6 +299,30 @@ read_int (FILE *fp, int *nread, const char *name) skip_whitespace_comments (fp); } +void +read_uint (FILE *fp, unsigned long int *ui) +{ + int n = 0; + + if (nextchar == EOF) + { + printf ("Error: Unexpected EOF when reading uint " + "in file '%s' line %lu\n", + pathname, line_number); + exit (1); + } + ungetc (nextchar, fp); + n = fscanf (fp, "%lu", ui); + if (ferror (fp) || n == 0 || n == EOF) + { + printf ("Error: Cannot read uint in file '%s' line %lu\n", + pathname, line_number); + exit (1); + } + nextchar = getc (fp); + skip_whitespace_comments (fp); +} + mpfr_prec_t read_mpfr_prec (FILE *fp) { @@ -458,6 +482,21 @@ read_ccf (FILE *fp, int *inex_re, int *inex_im, mpc_ptr expected, check_compatible (*inex_im, MPC_IM(expected), MPC_RND_IM(*rnd), "imag"); } +static void +read_ccu (FILE *fp, int *inex_re, int *inex_im, mpc_ptr expected, + known_signs_t *signs, mpc_ptr op1, unsigned long int *op2, mpc_rnd_t *rnd) +{ + test_line_number = line_number; + read_ternary (fp, inex_re); + read_ternary (fp, inex_im); + read_mpc (fp, expected, signs); + read_mpc (fp, op1, NULL); + read_uint (fp, op2); + read_mpc_rounding_mode (fp, rnd); + check_compatible (*inex_re, MPC_RE(expected), MPC_RND_RE(*rnd), "real"); + check_compatible (*inex_im, MPC_IM(expected), MPC_RND_IM(*rnd), "imag"); +} + /* data_check (function, data_file_name) checks function results against precomputed data in a file.*/ void @@ -473,6 +512,9 @@ data_check (mpc_function function, const char *file_name) int inex_im; mpc_t z1, z2, z3, z4; mpc_rnd_t rnd = MPC_RNDNN; + + unsigned long int ui; + known_signs_t signs; int inex = 0; @@ -486,7 +528,7 @@ data_check (mpc_function function, const char *file_name) mpfr_init (x1); mpfr_init (x2); break; - case CC: + case CC: case CCU: mpc_init2 (z2, 2); mpc_init2 (z3, 2); break; @@ -689,6 +731,34 @@ data_check (mpc_function function, const char *file_name) } break; + case CCU: /* example mpc_pow_ui */ + read_ccu (fp, &inex_re, &inex_im, z1, &signs, z2, &ui, &rnd); + mpfr_set_prec (MPC_RE(z3), MPC_PREC_RE (z1)); + mpfr_set_prec (MPC_IM(z3), MPC_PREC_IM (z1)); + inex = function.pointer.CCU (z3, z2, ui, rnd); + if (!MPC_INEX_CMP (inex_re, inex_im, inex) + || !same_mpc_value (z3, z1, signs)) + { + /* display sensible variable names */ + mpc_t op1, got, expected; + op1[0] = z2[0]; + expected[0]= z1[0]; + got[0] = z3[0]; + printf ("%s(op) failed (line %lu)\nwith rounding mode %s\n", + function.name, test_line_number, rnd_mode[rnd]); + if (!MPC_INEX_CMP (inex_re, inex_im, inex)) + printf("ternary value: got %s, expected (%s, %s)\n", + MPC_INEX_STR (inex), + MPFR_INEX_STR (inex_re), MPFR_INEX_STR (inex_im)); + OUT (op1); + printf ("op2 %lu\n ", ui); + OUT (got); + OUT (expected); + + exit (1); + } + break; + default: ; } @@ -702,7 +772,7 @@ data_check (mpc_function function, const char *file_name) mpfr_clear (x1); mpfr_clear (x2); break; - case CC: + case CC: case CCU: mpc_clear (z2); mpc_clear (z3); break; diff --git a/tests/tpow_ui.c b/tests/tpow_ui.c index feda442..4795849 100644 --- a/tests/tpow_ui.c +++ b/tests/tpow_ui.c @@ -1,6 +1,6 @@ /* test file for mpc_pow_ui. -Copyright (C) 2009 Paul Zimmermann +Copyright (C) 2009, 2010 Paul Zimmermann, Andreas Enge This file is part of the MPC Library. @@ -31,7 +31,7 @@ compare_mpc_pow (mp_prec_t pmax, int iter, unsigned long nbits) int i, inex_pow, inex_pow_ui; gmp_randstate_t state; mpc_rnd_t rnd; - + gmp_randinit_default (state); mpc_init3 (y, sizeof (unsigned long) * CHAR_BIT, MPFR_PREC_MIN); for (p = MPFR_PREC_MIN; p <= pmax; p++) @@ -107,7 +107,9 @@ main (int argc, char *argv[]) return 0; } + DECL_FUNC (CCU, f, mpc_pow_ui); test_start (); + data_check (f, "pow_ui.dat"); mpc_init2 (z, 5); mpc_set_ui_ui (z, 3, 2, MPC_RNDNN); |