diff options
| author | enge <enge@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2010-07-30 15:31:30 +0000 |
|---|---|---|
| committer | enge <enge@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2010-07-30 15:31:30 +0000 |
| commit | f5e44f7b814899e0c8a01cf4803e99c37de9fecb (patch) | |
| tree | 0c7e502aa06d2b5e7edda453b797e5d84673d9fc | |
| parent | 669f79ebbfb19e5548a8798a575e4cf2406161b1 (diff) | |
| download | mpc-f5e44f7b814899e0c8a01cf4803e99c37de9fecb.tar.gz | |
exp: closed bug http://lists.gforge.inria.fr/pipermail/mpc-discuss/2010-July/000757.html
high precision is needed for argument close to 0
git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@813 211d60ee-9f03-0410-a15a-8952a2c7a4e4
| -rw-r--r-- | NEWS | 1 | ||||
| -rw-r--r-- | src/exp.c | 16 | ||||
| -rw-r--r-- | tests/exp.dat | 3 |
3 files changed, 13 insertions, 7 deletions
@@ -3,6 +3,7 @@ Changes in version 0.8.3: mpc_pow_z when the exponent fits in a long - Bug fixes: - trigonometric functions: infinite loop due to overflow for large arguments + - exp: close to infinite loop for argument close to 0 Changes in version 0.8.2: - Speed-up of mpc_pow_ui through binary exponentiation @@ -148,8 +148,13 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* from now on, both parts of op are regular numbers */ - prec = MPC_MAX_PREC(rop); - + prec = MPC_MAX_PREC(rop) + + MPC_MAX (MPC_MAX (-mpfr_get_exp (MPC_RE (op)), 0), + -mpfr_get_exp (MPC_IM (op))); + /* When op is close to 0, then exp is close to 1+Re(op), while + cos is close to 1-Im(op); to decide on the ternary value of exp*cos, + we need a high enough precision so that none of exp or cos is + computed as 1. */ mpfr_init2 (x, 2); mpfr_init2 (y, 2); mpfr_init2 (z, 2); @@ -166,10 +171,9 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) could be represented in the precision of rop. */ mpfr_clear_overflow (); mpfr_clear_underflow (); - mpfr_exp (x, MPC_RE(op), GMP_RNDN); - mpfr_sin_cos (z, y, MPC_IM(op), GMP_RNDN); - mpfr_mul (y, y, x, GMP_RNDN); - + mpfr_exp (x, MPC_RE(op), GMP_RNDN); /* error <= 0.5ulp */ + mpfr_sin_cos (z, y, MPC_IM(op), GMP_RNDN); /* errors <= 0.5ulp */ + mpfr_mul (y, y, x, GMP_RNDN); /* error <= 2ulp */ ok = mpfr_overflow_p () || mpfr_zero_p (x) || mpfr_can_round (y, prec - 2, GMP_RNDN, GMP_RNDZ, MPFR_PREC(MPC_RE(rop)) + (MPC_RND_RE(rnd) == GMP_RNDN)); diff --git a/tests/exp.dat b/tests/exp.dat index f9c5fa3..5c9c978 100644 --- a/tests/exp.dat +++ b/tests/exp.dat @@ -1,6 +1,6 @@ # Data file for mpc_exp. # -# Copyright (C) 2008, 2010 Philippe Th\'eveny, Andreas Enge +# Copyright (C) 2008, 2010 Philippe Th\'eveny, Andreas Enge, Paul Zimmermann # # This file is part of the MPC Library. # @@ -116,4 +116,5 @@ # overflow - - 2 -inf 2 -inf 53 0x3312ae437f94441ec@-9 53 0xe45f7bab0595dd700@-10 N N +# input close to 0 ? ? 53 1 53 0x5D7A2148071Fp-7213522 53 0x1E02AE0D0F6Fp-7213521 53 0x5D7A2148071Fp-7213522 N N |