From cbed615f69ba2e239106fdaf3dd0f073ea50d109 Mon Sep 17 00:00:00 2001 From: zimmerma Date: Wed, 30 Sep 2009 09:04:39 +0000 Subject: [acos.c] completed case of regular argument [acosh.dat,asin.dat,atanh.dat,acos.dat,asinh.dat] added test cases git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/branches/feature-inverse-trigo@684 211d60ee-9f03-0410-a15a-8952a2c7a4e4 --- src/acos.c | 73 ++++++++++++++++++++++++++++++++++++++++++++++++++------- tests/acos.dat | 9 ++++--- tests/acosh.dat | 4 ++-- tests/asin.dat | 2 +- tests/asinh.dat | 6 +++-- tests/atanh.dat | 4 ++-- 6 files changed, 80 insertions(+), 18 deletions(-) diff --git a/src/acos.c b/src/acos.c index e7eee2e..a32eaa9 100644 --- a/src/acos.c +++ b/src/acos.c @@ -1,6 +1,6 @@ /* mpc_acos -- arccosine of a complex number. -Copyright (C) 2009 Philippe Th\'eveny +Copyright (C) 2009 Philippe Th\'eveny, Paul Zimmermann This file is part of the MPC Library. @@ -19,6 +19,7 @@ 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 /* for MPC_ASSERT */ #include "mpc-impl.h" extern int set_pi_over_2 (mpfr_ptr rop, int s, mpfr_rnd_t rnd); @@ -26,8 +27,13 @@ extern int set_pi_over_2 (mpfr_ptr rop, int s, mpfr_rnd_t rnd); int mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { - int inex_re; - int inex_im; + int inex_re, inex_im, inex; + mp_prec_t p_re, p_im, p; + mpc_t z1; + mpfr_t pi_over_2; + mp_exp_t e1, e2; + mp_rnd_t rnd_im; + mpc_rnd_t rnd1; inex_re = 0; inex_im = 0; @@ -167,9 +173,60 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) return MPC_INEX (inex_re, inex_im); } - /* regular complex argument */ - /* TODO */ - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_RE (rop)); - return 0; + /* regular complex argument: acos(z) = Pi/2 - asin(z) */ + p_re = mpfr_get_prec (MPC_RE(rop)); + p_im = mpfr_get_prec (MPC_IM(rop)); + p = p_re; + mpc_init3 (z1, p, p_im); /* we round directly the imaginary part to p_im, + with rounding mode opposite to rnd_im */ + rnd_im = MPC_RND_IM(rnd); + /* the imaginary part of asin(z) has the same sign as Im(z), thus if + Im(z) > 0 and rnd_im = RNDZ, we want to round the Im(asin(z)) to -Inf + so that -Im(asin(z)) is rounded to zero */ + if (rnd_im == GMP_RNDZ) + rnd_im = mpfr_sgn (MPC_IM(op)) > 0 ? GMP_RNDD : GMP_RNDU; + else + rnd_im = rnd_im == GMP_RNDU ? GMP_RNDD + : rnd_im == GMP_RNDD ? GMP_RNDU +#if MPFR_VERSION_MAJOR >= 3 + : rnd_im == GMP_RNDA ? GMP_RNDZ +#endif + : rnd_im; + rnd1 = RNDC(GMP_RNDN, rnd_im); + mpfr_init2 (pi_over_2, p); + for (;;) + { + p += mpc_ceil_log2 (p) + 3; + + mpfr_set_prec (MPC_RE(z1), p); + mpfr_set_prec (pi_over_2, p); + + mpfr_const_pi (pi_over_2, GMP_RNDN); + mpfr_div_2exp (pi_over_2, pi_over_2, 1, GMP_RNDN); /* Pi/2 */ + e1 = 1; /* Exp(pi_over_2) */ + inex = mpc_asin (z1, op, rnd1); /* asin(z) */ + MPC_ASSERT (mpfr_sgn (MPC_IM(z1)) * mpfr_sgn (MPC_IM(op)) > 0); + 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; + } + inex = mpc_set (rop, z1, rnd); + inex_re = MPC_INEX_RE(inex); + mpc_clear (z1); + mpfr_clear (pi_over_2); + + return MPC_INEX(inex_re, inex_im); } diff --git a/tests/acos.dat b/tests/acos.dat index 55508e4..3a51faa 100644 --- a/tests/acos.dat +++ b/tests/acos.dat @@ -1,6 +1,6 @@ # Data file for mpc_acos. # -# Copyright (C) 2009 Philippe Th\'eveny +# Copyright (C) 2009 Philippe Th\'eveny, Paul Zimmermann # # This file is part of the MPC Library. # @@ -122,5 +122,8 @@ - - 53 0x1921FB54442D18p-52 53 -0x1D185B507EDC0Ep-52 53 -0 53 +3 N N - - 53 0x1921FB54442D18p-52 53 -0x1D185B507EDC0Ep-52 53 +0 53 +3 N N -# IEEE-754 double precision - +# regular argument for various precisions +- + 2 0.5 2 -1 2 2 2 1 N Z ++ + 9 0x5Dp-6 9 0x9Fp-5 9 8.5 9 -71 N U ++ + 2 0x3p-9 2 1.5 2 2 2 -0x1p-7 U N ++ - 53 0x74C141310E695p-53 53 -0x1D6D2CFA9F3F11p-52 53 0x3243F6A8885A3p-48 53 0x162E42FEFA39EFp-53 N N diff --git a/tests/acosh.dat b/tests/acosh.dat index 68fea87..c9b2ead 100644 --- a/tests/acosh.dat +++ b/tests/acosh.dat @@ -1,6 +1,6 @@ # Data file for mpc_acosh. # -# Copyright (C) 2009 Philippe Th\'eveny +# Copyright (C) 2009 Philippe Th\'eveny, Paul Zimmermann # # This file is part of the MPC Library. # @@ -120,4 +120,4 @@ - - 53 0x1C37C174A83DEDp-51 53 0x1921FB54442D18p-52 53 +0 53 +17 N N # IEEE-754 double precision - ++ + 53 0x1D6D2CFA9F3F11p-52 53 0x74C141310E695p-53 53 0x3243F6A8885A3p-48 53 0x162E42FEFA39EFp-53 N N diff --git a/tests/asin.dat b/tests/asin.dat index cd7efcc..d41d3f3 100644 --- a/tests/asin.dat +++ b/tests/asin.dat @@ -1,6 +1,6 @@ # Data file for mpc_asin. # -# Copyright (C) 2009 Philippe Th\'eveny +# Copyright (C) 2009 Philippe Th\'eveny, Paul Zimmermann # # This file is part of the MPC Library. # diff --git a/tests/asinh.dat b/tests/asinh.dat index 8aa1395..e8d2f53 100644 --- a/tests/asinh.dat +++ b/tests/asinh.dat @@ -1,6 +1,6 @@ # Data file for mpc_asinh. # -# Copyright (C) 2009 Philippe Th\'eveny +# Copyright (C) 2009 Philippe Th\'eveny, Paul Zimmermann # # This file is part of the MPC Library. # @@ -116,5 +116,7 @@ - - 53 -0x1ECC2CAEC5160Ap-53 53 0x1921FB54442D18p-52 53 -0 53 +1.5 N N + - 53 0x1ECC2CAEC5160Ap-53 53 0x1921FB54442D18p-52 53 +0 53 +1.5 N N -# IEEE-754 double precision +# regular arguments ++ + 53 0x1E20C7792ECE6Bp-52 53 0x3526776219EEBp-52 53 0x3243F6A8885A3p-48 53 0x162E42FEFA39EFp-53 N N + diff --git a/tests/atanh.dat b/tests/atanh.dat index 1b4227d..15f873e 100644 --- a/tests/atanh.dat +++ b/tests/atanh.dat @@ -1,6 +1,6 @@ # Data file for mpc_atanh. # -# Copyright (C) 2009 Philippe Th\'eveny +# Copyright (C) 2009 Philippe Th\'eveny, Paul Zimmermann # # This file is part of the MPC Library. # @@ -109,4 +109,4 @@ 0 - 53 +0 53 +0x1F5B75F92C80DDp-55 53 +0 53 +.25 N N # IEEE-754 double precision - +- + 53 0x13F3F785301CE9p-54 53 0xBFA43C2A868B3p-51 53 0x3243F6A8885A3p-48 53 0x162E42FEFA39EFp-53 N N -- cgit v1.2.1