diff options
-rw-r--r-- | BUGS | 4 | ||||
-rw-r--r-- | src/acos.c | 31 | ||||
-rw-r--r-- | src/atan.c | 35 | ||||
-rw-r--r-- | src/tan.c | 44 | ||||
-rw-r--r-- | tests/tan.dat | 8 |
5 files changed, 80 insertions, 42 deletions
@@ -9,7 +9,3 @@ sin(x)*sinh(y) is representable. If furthermore an underflow occurred in sin(x) (which has not been observed in practice), then the return value would be NaN*(+-inf)=NaN. - As another example, tan is computed as sin/cos; if there is an overflow - in both sin and cos, then NaN+i*NaN is returned even if the result - may be representable. - @@ -1,6 +1,6 @@ /* mpc_acos -- arccosine of a complex number. -Copyright (C) 2009 Philippe Th\'eveny, Paul Zimmermann +Copyright (C) 2009, 2010 Philippe Th\'eveny, Paul Zimmermann This file is part of the MPC Library. @@ -206,19 +206,22 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) inex_im = MPC_INEX_IM(inex); /* inex_im is in {-1, 0, 1} */ e2 = mpfr_get_exp (MPC_RE(z1)); mpfr_sub (MPC_RE(z1), pi_over_2, MPC_RE(z1), GMP_RNDN); - /* the error on x=Re(z1) is bounded by 1/2 ulp(x) + 2^(e1-p-1) + - 2^(e2-p-1) */ - e1 = e1 >= e2 ? e1 + 1 : e2 + 1; - /* the error on x is bounded by 1/2 ulp(x) + 2^(e1-p-1) */ - e1 -= mpfr_get_exp (MPC_RE(z1)); - /* the error on x is bounded by 1/2 ulp(x) [1 + 2^e1] */ - e1 = e1 <= 0 ? 0 : e1; - /* the error on x is bounded by 2^e1 * ulp(x) */ - mpfr_neg (MPC_IM(z1), MPC_IM(z1), GMP_RNDN); /* exact */ - inex_im = -inex_im; - if (mpfr_can_round (MPC_RE(z1), p - e1, GMP_RNDN, GMP_RNDZ, - p_re + (MPC_RND_RE(rnd) == GMP_RNDN))) - break; + if (!mpfr_zero_p (MPC_RE(z1))) + { + /* the error on x=Re(z1) is bounded by 1/2 ulp(x) + 2^(e1-p-1) + + 2^(e2-p-1) */ + e1 = e1 >= e2 ? e1 + 1 : e2 + 1; + /* the error on x is bounded by 1/2 ulp(x) + 2^(e1-p-1) */ + e1 -= mpfr_get_exp (MPC_RE(z1)); + /* the error on x is bounded by 1/2 ulp(x) [1 + 2^e1] */ + e1 = e1 <= 0 ? 0 : e1; + /* the error on x is bounded by 2^e1 * ulp(x) */ + mpfr_neg (MPC_IM(z1), MPC_IM(z1), GMP_RNDN); /* exact */ + inex_im = -inex_im; + if (mpfr_can_round (MPC_RE(z1), p - e1, GMP_RNDN, GMP_RNDZ, + p_re + (MPC_RND_RE(rnd) == GMP_RNDN))) + break; + } } inex = mpc_set (rop, z1, rnd); inex_re = MPC_INEX_RE(inex); @@ -1,6 +1,6 @@ /* mpc_atan -- arctangent of a complex number. -Copyright (C) 2009 Philippe Th\'eveny, Paul Zimmermann +Copyright (C) 2009, 2010 Philippe Th\'eveny, Paul Zimmermann This file is part of the MPC Library. @@ -331,21 +331,26 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_sub (y, a, b, GMP_RNDU); - expo = MPC_MAX (mpfr_get_exp (a), mpfr_get_exp (b)); - expo = expo - mpfr_get_exp (y) + 1; - err = 3 - mpfr_get_exp (y); - /* error(j) <= [1 + 2^expo + 7*2^err] ulp(j) */ - if (expo <= err) /* error(j) <= [1 + 2^{err+1}] ulp(j) */ - err = (err < 0) ? 1 : err + 2; - else - err = (expo < 0) ? 1 : expo + 2; - - mpfr_div_2ui (y, y, 2, GMP_RNDN); if (mpfr_zero_p (y)) - err = 2; /* underflow */ - - ok = mpfr_can_round (y, p - err, GMP_RNDU, GMP_RNDD, - prec + (MPC_RND_IM (rnd) == GMP_RNDN)); + ok = 0; + else + { + expo = MPC_MAX (mpfr_get_exp (a), mpfr_get_exp (b)); + expo = expo - mpfr_get_exp (y) + 1; + err = 3 - mpfr_get_exp (y); + /* error(j) <= [1 + 2^expo + 7*2^err] ulp(j) */ + if (expo <= err) /* error(j) <= [1 + 2^{err+1}] ulp(j) */ + err = (err < 0) ? 1 : err + 2; + else + err = (expo < 0) ? 1 : expo + 2; + + mpfr_div_2ui (y, y, 2, GMP_RNDN); + if (mpfr_zero_p (y)) + err = 2; /* underflow */ + + ok = mpfr_can_round (y, p - err, GMP_RNDU, GMP_RNDD, + prec + (MPC_RND_IM (rnd) == GMP_RNDN)); + } } while (ok == 0); inex = mpc_set_fr_fr (rop, x, y, rnd); @@ -1,6 +1,6 @@ /* mpc_tan -- tangent of a complex number. -Copyright (C) 2008, 2009, 2010 Philippe Th\'eveny, Andreas Enge +Copyright (C) 2008, 2009, 2010 Philippe Th\'eveny, Andreas Enge, Paul Zimmermann This file is part of the MPC Library. @@ -19,6 +19,8 @@ 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. */ +#include <stdio.h> /* for MPC_ASSERT */ +#include <limits.h> #include "mpc-impl.h" int @@ -181,6 +183,7 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) do { mpfr_exp_t k, exr, eyr, eyi, ezr; + int isinf = 0; ok = 0; @@ -193,16 +196,49 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) above), sin x and cos y are never exact, so rounding away from 0 is rounding towards 0 and adding one ulp to the absolute value */ mpc_sin (x, op, MPC_RNDZZ); + isinf = mpfr_inf_p (MPC_RE (x)) || mpfr_inf_p (MPC_IM (x)); + /* if the real or imaginary part of x is infinite, it means that 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)) */ mpfr_signbit (MPC_RE (x)) ? mpfr_nextbelow (MPC_RE (x)) : mpfr_nextabove (MPC_RE (x)); mpfr_signbit (MPC_IM (x)) ? mpfr_nextbelow (MPC_IM (x)) : mpfr_nextabove (MPC_IM (x)); + MPC_ASSERT (mpfr_zero_p (MPC_RE (x)) == 0); exr = mpfr_get_exp (MPC_RE (x)); mpc_cos (y, op, MPC_RNDZZ); + isinf = isinf || mpfr_inf_p (MPC_RE (y)) || mpfr_inf_p (MPC_IM (y)); mpfr_signbit (MPC_RE (y)) ? mpfr_nextbelow (MPC_RE (y)) : mpfr_nextabove (MPC_RE (y)); mpfr_signbit (MPC_IM (y)) ? mpfr_nextbelow (MPC_IM (y)) : mpfr_nextabove (MPC_IM (y)); + + if (isinf) + { + int inex_re, inex_im; + mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN); + if (mpfr_sgn (MPC_RE (x)) * mpfr_sgn (MPC_RE (y)) < 0) + { + mpfr_neg (MPC_RE (rop), MPC_RE (rop), GMP_RNDN); + inex_re = 1; + } + else + inex_re = -1; /* +0 is rounded down */ + if (mpfr_sgn (MPC_IM (op)) > 0) + { + mpfr_set_ui (MPC_IM (rop), 1, GMP_RNDN); + inex_im = 1; + } + else + { + mpfr_set_si (MPC_IM (rop), -1, GMP_RNDN); + inex_im = -1; + } + inex = MPC_INEX(inex_re, inex_im); + goto end; + } + eyr = mpfr_get_exp (MPC_RE (y)); eyi = mpfr_get_exp (MPC_IM (y)); @@ -225,6 +261,7 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) if (MPC_INEX_IM (inex)) mpfr_signbit (MPC_IM (x)) ? mpfr_nextbelow (MPC_IM (x)) : mpfr_nextabove (MPC_IM (x)); + MPC_ASSERT (mpfr_zero_p (MPC_RE (x)) == 0); ezr = mpfr_get_exp (MPC_RE (x)); /* FIXME: compute @@ -233,14 +270,14 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) k = exr - ezr + MPC_MAX(-eyr, eyr - 2 * eyi); err = k < 2 ? 7 : (k == 2 ? 8 : (5 + k)); - /* Can the real part be rounded ? */ + /* Can the real part be rounded? */ ok = (!mpfr_number_p (MPC_RE (x))) || mpfr_can_round (MPC_RE(x), prec - err, GMP_RNDN, GMP_RNDZ, MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == GMP_RNDN)); if (ok) { - /* Can the imaginary part be rounded ? */ + /* Can the imaginary part be rounded? */ ok = (!mpfr_number_p (MPC_IM (x))) || mpfr_can_round (MPC_IM(x), prec - 6, GMP_RNDN, GMP_RNDZ, MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == GMP_RNDN)); @@ -250,6 +287,7 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) inex = mpc_set (rop, x, rnd); + end: mpc_clear (x); mpc_clear (y); diff --git a/tests/tan.dat b/tests/tan.dat index 974c882..0a97e48 100644 --- a/tests/tan.dat +++ b/tests/tan.dat @@ -126,9 +126,5 @@ - + 9 -0x9bp-51 9 -1 9 -0x16dp-8 9 -0x77p-3 N N # huge values -# should yield finite result as in commented lines, in particular simple imaginary part; -# in current implementation, yields NaN due to intermediate overflow in sin and cos -#- + 53 -0xE72F36A66DDF7p-2765858364927957 53 -1 53 0x4580CBF242683p-3 53 -0x1B3E8A3660D279p-3 N N -#- - 53 -0xD10502370373Fp-441000018035230 53 1 53 -0x1B3E8A3660D279p-3 53 0x4580CBF242683p-3 N N -0 0 53 nan 53 nan 53 0x4580CBF242683p-3 53 -0x1B3E8A3660D279p-3 N N -0 0 53 nan 53 nan 53 -0x1B3E8A3660D279p-3 53 0x4580CBF242683p-3 N N ++ - 53 -0 53 -1 53 0x4580CBF242683p-3 53 -0x1B3E8A3660D279p-3 N N ++ + 53 -0 53 +1 53 -0x1B3E8A3660D279p-3 53 0x4580CBF242683p-3 N N |