From ff08a6c7094671336598d8c588dc882dbe61b822 Mon Sep 17 00:00:00 2001 From: thevenyp Date: Wed, 5 Aug 2009 15:33:26 +0000 Subject: src/atan.c: compute special values and pure real/ pure imaginary arguments. src/atanh.c: fix inversions in rounding mode and nonary value. tests/atan.dat: Fix special values, add tests with pure real and pure imaginary argument. tests/atanh.dat: Fix undetermined signs of zero, add tests with pure real and pure imaginary argument. git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/branches/feature-inverse-trigo@642 211d60ee-9f03-0410-a15a-8952a2c7a4e4 --- src/atan.c | 185 +++++++++++++++++++++++++++++++++++++++++++++++++++++++- src/atanh.c | 4 +- tests/atan.dat | 92 ++++++++++++++++++++-------- tests/atanh.dat | 34 ++++++++--- 4 files changed, 277 insertions(+), 38 deletions(-) diff --git a/src/atan.c b/src/atan.c index 60d6c6c..b394b28 100644 --- a/src/atan.c +++ b/src/atan.c @@ -21,10 +21,193 @@ MA 02111-1307, USA. */ #include "mpc-impl.h" +/* set rop to + -pi/2 if s < 0 + +pi/2 else + rounded in the direction rnd +*/ +static int +set_pi_over_2 (mpfr_ptr rop, int s, mpfr_rnd_t rnd) +{ + int inex; + + inex = mpfr_const_pi (rop, s ? INV_RND (rnd) : rnd); + mpfr_div_2ui (rop, rop, 1, GMP_RNDN); + if (s) + { + inex = -inex; + mpfr_neg (rop, rop, GMP_RNDN); + } + + return inex; +} + int mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { + int s_re; + int s_im; + + s_re = mpfr_signbit (MPC_RE (op)); + s_im = mpfr_signbit (MPC_IM (op)); + + /* special values */ + if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op))) + { + int inex_re; + inex_re = 0; + + if (mpfr_nan_p (MPC_RE (op))) + { + mpfr_set_nan (MPC_RE (rop)); + if (mpfr_zero_p (MPC_IM (op)) || mpfr_inf_p (MPC_IM (op))) + { + mpfr_set_ui (MPC_IM (rop), 0, GMP_RNDN); + if (s_im) + mpc_conj (rop, rop, MPC_RNDNN); + } + else + mpfr_set_nan (MPC_IM (rop)); + } + else + { + if (mpfr_inf_p (MPC_RE (op))) + { + inex_re = set_pi_over_2 (MPC_RE (rop), s_re, MPC_RND_RE (rnd)); + mpfr_set_ui (MPC_IM (rop), 0, GMP_RNDN); + } + else + { + mpfr_set_nan (MPC_RE (rop)); + mpfr_set_nan (MPC_IM (rop)); + } + } + return MPC_INEX (inex_re, 0); + } + + if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op))) + { + int inex_re; + inex_re = set_pi_over_2 (MPC_RE (rop), s_re, MPC_RND_RE (rnd)); + + mpfr_set_ui (MPC_IM (rop), 0, GMP_RNDN); + if (s_im) + mpc_conj (rop, rop, GMP_RNDN); + + return MPC_INEX (inex_re, 0); + } + + /* pure real argument */ + if (mpfr_zero_p (MPC_IM (op))) + { + int inex_re; + + inex_re = mpfr_atan (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd)); + + mpfr_set_ui (MPC_IM (rop), 0, GMP_RNDN); + if (s_im) + mpc_conj (rop, rop, GMP_RNDN); + + return MPC_INEX (inex_re, 0); + } + + /* pure imaginary argument */ + if (mpfr_zero_p (MPC_RE (op))) + { + int inex_re; + int inex_im; + int cmp_1; + + inex_re = 0; + inex_im = 0; + + if (s_im) + cmp_1 = -mpfr_cmp_si (MPC_IM (op), -1); + else + cmp_1 = mpfr_cmp_ui (MPC_IM (op), +1); + + if (cmp_1 < 0) + { + /* atan(+0+iy) = +0 +i*atanh(y), if |y| < 1 + atan(-0+iy) = -0 +i*atanh(y), if |y| < 1 */ + + mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN); + if (s_re) + mpfr_neg (MPC_RE (rop), MPC_RE (rop), GMP_RNDN); + + inex_im = mpfr_atanh (MPC_IM (rop), MPC_IM (op), MPC_RND_IM (rnd)); + } + else if (cmp_1 == 0) + { + /* atan(+/-0+i) = NaN +i*inf + atan(+/-0-i) = NaN -i*inf */ + mpfr_set_nan (MPC_RE (rop)); + mpfr_set_inf (MPC_IM (rop), s_im ? -1 : +1); + } + else + { + /* atan(+0+iy) = +pi/2 +i*atanh(1/y), if |y| > 1 + atan(-0+iy) = -pi/2 +i*atanh(1/y), if |y| > 1 */ + mp_rnd_t rnd_im; + mpfr_t y; + mp_prec_t p, p_im; + mp_prec_t err; + int ok; + int cmp_2; + + rnd_im = MPC_RND_IM (rnd); + mpfr_init (y); + p_im = mpfr_get_prec (MPC_IM (op)); + p = p_im; + + /* a = o(1/y) with error(a) < 1 ulp(a) + b = o(atanh(a)) with error(b) < (1+2^{1+Exp(a)-Exp(b)}) ulp(b) + + * if 1 < |y| <= 2 then Exp(a) = 0 and Exp(b) >= 0, so, at most, 2 + bits of precision are lost. + * if |y| > 2, no simplification. */ + if (s_im) + cmp_2 = -mpfr_cmp_si (MPC_IM (op), -2); + else + cmp_2 = mpfr_cmp_ui (MPC_IM (op), +2); + + err = 2; + + do + { + p += mpc_ceil_log2 (p) + err; + mpfr_set_prec (y, p); + inex_im = mpfr_ui_div (y, 1, MPC_IM (op), GMP_RNDZ); + if (mpfr_zero_p (y)) + break; /* underflow. */ + /* FIXME: should we consider the case with unreasonably huge + precision prec(y)>3*exp_min, where atanh(1/Im(op)) could be + representable while 1/Im(op) underflows ? */ + if (cmp_2 > 0) + err = 2 + (mp_prec_t)MPC_MAX (MPFR_EXP (y), 0); + + /* atanh cannot underflow: |atanh(x)| > |x| for |x| < 1 */ + inex_im |= mpfr_atanh (y, y, GMP_RNDZ); + if (cmp_2 > 0) + err += MPC_MAX (SAFE_ABS (mp_prec_t, MPFR_EXP (y)), 0); + + ok = inex_im == 0 + || mpfr_can_round (y, p-err, GMP_RNDZ, rnd_im, + p_im + (rnd_im == GMP_RNDN)); + } while (ok == 0); + + inex_re = set_pi_over_2 (MPC_RE (rop), s_re, MPC_RND_RE (rnd)); + inex_im = mpfr_set (MPC_IM (rop), y, rnd_im); + mpfr_clear (y); + } + return MPC_INEX (inex_re, inex_im); + } + + /* regular number argument */ + + /* TODO */ mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_RE (rop)); + mpfr_set_nan (MPC_IM (rop)); + return 0; } diff --git a/src/atanh.c b/src/atanh.c index 057462a..bb2ad1e 100644 --- a/src/atanh.c +++ b/src/atanh.c @@ -34,7 +34,7 @@ mpc_atanh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) MPFR_CHANGE_SIGN (MPC_RE (z)); inex = mpc_atan (rop, z, - RNDC (MPC_RND_IM (rnd), INV_RND (MPC_RND_RE (rnd)))); + RNDC (INV_RND (MPC_RND_IM (rnd)), MPC_RND_RE (rnd))); /* change rop to -i*rop */ tmp[0] = MPC_RE (rop)[0]; @@ -42,5 +42,5 @@ mpc_atanh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) MPC_IM (rop)[0] = tmp[0]; MPFR_CHANGE_SIGN (MPC_IM (rop)); - return MPC_INEX (-MPC_INEX_IM (inex), MPC_INEX_RE (inex)); + return MPC_INEX (MPC_INEX_IM (inex), -MPC_INEX_RE (inex)); } diff --git a/tests/atan.dat b/tests/atan.dat index c0029b3..bffec72 100644 --- a/tests/atan.dat +++ b/tests/atan.dat @@ -50,47 +50,87 @@ # The sign of the result is checked with "+inf", "-inf", "-0", or "+0". # special values (following ISO C99 standard) -+ 0 53 -0x1921FB54442D18p-52 53 -inf 53 -inf 53 -inf N N -+ 0 53 -0x1921FB54442D18p-52 53 -inf 53 -inf 53 -1 N N -+ 0 53 -0x1921FB54442D18p-52 53 -inf 53 -inf 53 -0 N N -+ 0 53 -0x1921FB54442D18p-52 53 +inf 53 -inf 53 +0 N N -+ 0 53 -0x1921FB54442D18p-52 53 +inf 53 -inf 53 +1 N N -+ 0 53 -0x1921FB54442D18p-52 53 +inf 53 -inf 53 +inf N N -+ 0 53 -0x1921FB54442D18p-52 53 inf 53 -inf 53 nan N N ++ 0 53 -0x1921FB54442D18p-52 53 -0 53 -inf 53 -inf N N ++ 0 53 -0x1921FB54442D18p-52 53 -0 53 -inf 53 -1 N N ++ 0 53 -0x1921FB54442D18p-52 53 -0 53 -inf 53 -0 N N ++ 0 53 -0x1921FB54442D18p-52 53 +0 53 -inf 53 +0 N N ++ 0 53 -0x1921FB54442D18p-52 53 +0 53 -inf 53 +1 N N ++ 0 53 -0x1921FB54442D18p-52 53 +0 53 -inf 53 +inf N N ++ 0 53 -0x1921FB54442D18p-52 53 0 53 -inf 53 nan N N + 0 53 -0x1921FB54442D18p-52 53 -0 53 -6 53 -inf N N + 0 53 -0x1921FB54442D18p-52 53 +0 53 -6 53 +inf N N 0 0 53 nan 53 nan 53 -6 53 nan N N + 0 53 -0x1921FB54442D18p-52 53 -0 53 -0 53 -inf N N -+ 0 53 -0x1921FB54442D18p-52 53 -0 53 -0 53 -0 N N -+ 0 53 -0x1921FB54442D18p-52 53 +0 53 -0 53 +0 N N +0 0 53 -0 53 -0 53 -0 53 -0 N N +0 0 53 -0 53 +0 53 -0 53 +0 N N + 0 53 -0x1921FB54442D18p-52 53 +0 53 -0 53 +inf N N 0 0 53 nan 53 nan 53 -0 53 nan N N -+ 0 53 +0x1921FB54442D18p-52 53 -0 53 +0 53 -inf N N -+ 0 53 +0x1921FB54442D18p-52 53 -0 53 +0 53 -0 N N -+ 0 53 +0x1921FB54442D18p-52 53 +0 53 +0 53 +0 N N -+ 0 53 +0x1921FB54442D18p-52 53 +0 53 +0 53 +inf N N +- 0 53 +0x1921FB54442D18p-52 53 -0 53 +0 53 -inf N N +0 0 53 +0 53 -0 53 +0 53 -0 N N +0 0 53 +0 53 +0 53 +0 53 +0 N N +- 0 53 +0x1921FB54442D18p-52 53 +0 53 +0 53 +inf N N 0 0 53 nan 53 nan 53 +0 53 nan N N -+ 0 53 +0x1921FB54442D18p-52 53 -0 53 +6 53 -inf N N -+ 0 53 +0x1921FB54442D18p-52 53 +0 53 +6 53 +inf N N +- 0 53 +0x1921FB54442D18p-52 53 -0 53 +6 53 -inf N N +- 0 53 +0x1921FB54442D18p-52 53 +0 53 +6 53 +inf N N 0 0 53 nan 53 nan 53 +6 53 nan N N -+ 0 53 +0x1921FB54442D18p-52 53 -0 53 +inf 53 -inf N N -+ 0 53 +0x1921FB54442D18p-52 53 -0 53 +inf 53 -1 N N -+ 0 53 +0x1921FB54442D18p-52 53 -0 53 +inf 53 -0 N N -+ 0 53 +0x1921FB54442D18p-52 53 +0 53 +inf 53 +0 N N -+ 0 53 +0x1921FB54442D18p-52 53 +0 53 +inf 53 +1 N N -+ 0 53 +0x1921FB54442D18p-52 53 +0 53 +inf 53 +inf N N -0 0 53 +0x1921FB54442D18p-52 53 nan 53 +inf 53 nan N N -0 0 53 +inf 53 nan 53 nan 53 -inf N N +- 0 53 +0x1921FB54442D18p-52 53 -0 53 +inf 53 -inf N N +- 0 53 +0x1921FB54442D18p-52 53 -0 53 +inf 53 -1 N N +- 0 53 +0x1921FB54442D18p-52 53 -0 53 +inf 53 -0 N N +- 0 53 +0x1921FB54442D18p-52 53 +0 53 +inf 53 +0 N N +- 0 53 +0x1921FB54442D18p-52 53 +0 53 +inf 53 +1 N N +- 0 53 +0x1921FB54442D18p-52 53 +0 53 +inf 53 +inf N N +- 0 53 +0x1921FB54442D18p-52 53 0 53 +inf 53 nan N N +0 0 53 nan 53 -0 53 nan 53 -inf N N 0 0 53 nan 53 nan 53 nan 53 -1 N N -0 0 53 nan 53 nan 53 nan 53 -0 N N -0 0 53 nan 53 nan 53 nan 53 +0 N N +0 0 53 nan 53 -0 53 nan 53 -0 N N +0 0 53 nan 53 +0 53 nan 53 +0 N N 0 0 53 nan 53 nan 53 nan 53 +1 N N -0 0 53 +inf 53 nan 53 nan 53 +inf N N +0 0 53 nan 53 +0 53 nan 53 +inf N N 0 0 53 nan 53 nan 53 nan 53 nan N N # pure real argument +- 0 53 -0x16DCC57BB565FDp-52 53 -0 53 -7 53 -0 N N +- 0 53 -0x16DCC57BB565FDp-52 53 +0 53 -7 53 +0 N N ++ 0 53 -0x1F730BD281F69Bp-53 53 -0 53 -1.5 53 -0 N N ++ 0 53 -0x1F730BD281F69Bp-53 53 +0 53 -1.5 53 +0 N N ++ 0 53 -0x1921FB54442D18p-53 53 -0 53 -1 53 -0 N N ++ 0 53 -0x1921FB54442D18p-53 53 +0 53 -1 53 +0 N N +- 0 53 -0x1700A7C5784634p-53 53 -0 53 -0.875 53 -0 N N +- 0 53 -0x1700A7C5784634p-53 53 +0 53 -0.875 53 +0 N N +- 0 53 -0x1FD5BA9AAC2F6Ep-56 53 -0 53 -0.125 53 -0 N N +- 0 53 -0x1FD5BA9AAC2F6Ep-56 53 +0 53 -0.125 53 +0 N N ++ 0 53 +0x1FD5BA9AAC2F6Ep-56 53 +0 53 +0.125 53 +0 N N ++ 0 53 +0x1FD5BA9AAC2F6Ep-56 53 -0 53 +0.125 53 -0 N N ++ 0 53 +0x1700A7C5784634p-53 53 +0 53 +0.875 53 +0 N N ++ 0 53 +0x1700A7C5784634p-53 53 -0 53 +0.875 53 -0 N N +- 0 53 +0x1921FB54442D18p-53 53 +0 53 +1 53 +0 N N +- 0 53 +0x1921FB54442D18p-53 53 -0 53 +1 53 -0 N N +- 0 53 +0x1F730BD281F69Bp-53 53 +0 53 +1.5 53 +0 N N +- 0 53 +0x1F730BD281F69Bp-53 53 -0 53 +1.5 53 -0 N N ++ 0 53 +0x16DCC57BB565FDp-52 53 +0 53 +7 53 +0 N N ++ 0 53 +0x16DCC57BB565FDp-52 53 -0 53 +7 53 -0 N N # pure imaginary argument ++ + 53 -0x1921FB54442D18p-52 53 -0x1269621134DB92p-55 53 -0 53 -7 N N +- + 53 +0x1921FB54442D18p-52 53 -0x1269621134DB92p-55 53 +0 53 -7 N N ++ + 53 -0x1921FB54442D18p-52 53 -0x19C041F7ED8D33p-53 53 -0 53 -1.5 N N +- + 53 +0x1921FB54442D18p-52 53 -0x19C041F7ED8D33p-53 53 +0 53 -1.5 N N +0 0 53 nan 53 -inf 53 -0 53 -1 N N +0 0 53 nan 53 -inf 53 +0 53 -1 N N +0 + 53 -0 53 -0x15AA16394D481Fp-52 53 -0 53 -0.875 N N +0 + 53 +0 53 -0x15AA16394D481Fp-52 53 +0 53 -0.875 N N +0 + 53 -0 53 -0x1015891C9EAEF7p-55 53 -0 53 -0.125 N N +0 + 53 +0 53 -0x1015891C9EAEF7p-55 53 +0 53 -0.125 N N +0 - 53 +0 53 +0x1015891C9EAEF7p-55 53 +0 53 +0.125 N N +0 - 53 -0 53 +0x1015891C9EAEF7p-55 53 -0 53 +0.125 N N +0 - 53 +0 53 +0x15AA16394D481Fp-52 53 +0 53 +0.875 N N +0 - 53 -0 53 +0x15AA16394D481Fp-52 53 -0 53 +0.875 N N +0 0 53 nan 53 +inf 53 +0 53 +1 N N +0 0 53 nan 53 +inf 53 -0 53 +1 N N +- - 53 +0x1921FB54442D18p-52 53 +0x19C041F7ED8D33p-53 53 +0 53 +1.5 N N ++ - 53 -0x1921FB54442D18p-52 53 +0x19C041F7ED8D33p-53 53 -0 53 +1.5 N N +- - 53 +0x1921FB54442D18p-52 53 +0x1269621134DB92p-55 53 +0 53 +7 N N ++ - 53 -0x1921FB54442D18p-52 53 +0x1269621134DB92p-55 53 -0 53 +7 N N # IEEE-754 double precision diff --git a/tests/atanh.dat b/tests/atanh.dat index ad075c2..1b4227d 100644 --- a/tests/atanh.dat +++ b/tests/atanh.dat @@ -57,13 +57,13 @@ 0 - 53 -0 53 +0x1921FB54442D18p-52 53 -inf 53 +1 N N 0 - 53 -0 53 +0x1921FB54442D18p-52 53 -inf 53 +inf N N 0 0 53 -0 53 nan 53 -inf 53 nan N N -0 + 53 +0 53 -0x1921FB54442D18p-52 53 -6 53 -inf N N -0 - 53 +0 53 0x1921FB54442D18p-52 53 -6 53 +inf N N +0 + 53 -0 53 -0x1921FB54442D18p-52 53 -6 53 -inf N N +0 - 53 -0 53 +0x1921FB54442D18p-52 53 -6 53 +inf N N 0 0 53 nan 53 nan 53 -6 53 nan N N 0 + 53 -0 53 -0x1921FB54442D18p-52 53 -0 53 -inf N N 0 0 53 -0 53 -0 53 -0 53 -0 N N 0 0 53 -0 53 +0 53 -0 53 +0 N N -0 - 53 -0 53 0x1921FB54442D18p-52 53 -0 53 +inf N N +0 - 53 -0 53 +0x1921FB54442D18p-52 53 -0 53 +inf N N 0 0 53 -0 53 nan 53 -0 53 nan N N 0 + 53 +0 53 -0x1921FB54442D18p-52 53 +0 53 -inf N N 0 0 53 +0 53 -0 53 +0 53 -0 N N @@ -71,26 +71,42 @@ 0 - 53 +0 53 +0x1921FB54442D18p-52 53 +0 53 +inf N N 0 0 53 +0 53 nan 53 +0 53 nan N N 0 + 53 +0 53 -0x1921FB54442D18p-52 53 +6 53 -inf N N -0 - 53 +0 53 0x1921FB54442D18p-52 53 +6 53 +inf N N +0 - 53 +0 53 +0x1921FB54442D18p-52 53 +6 53 +inf N N 0 0 53 nan 53 nan 53 +6 53 nan N N 0 + 53 +0 53 -0x1921FB54442D18p-52 53 +inf 53 -inf N N 0 + 53 +0 53 -0x1921FB54442D18p-52 53 +inf 53 -1 N N 0 + 53 +0 53 -0x1921FB54442D18p-52 53 +inf 53 -0 N N -0 - 53 +0 53 0x1921FB54442D18p-52 53 +inf 53 +0 N N -0 - 53 +0 53 0x1921FB54442D18p-52 53 +inf 53 +1 N N -0 - 53 +0 53 0x1921FB54442D18p-52 53 +inf 53 +inf N N +0 - 53 +0 53 +0x1921FB54442D18p-52 53 +inf 53 +0 N N +0 - 53 +0 53 +0x1921FB54442D18p-52 53 +inf 53 +1 N N +0 - 53 +0 53 +0x1921FB54442D18p-52 53 +inf 53 +inf N N 0 0 53 +0 53 nan 53 +inf 53 nan N N -0 + 53 +0 53 -0x1921FB54442D18p-52 53 nan 53 -inf N N +0 + 53 0 53 -0x1921FB54442D18p-52 53 nan 53 -inf N N 0 0 53 nan 53 nan 53 nan 53 -1 N N 0 0 53 nan 53 nan 53 nan 53 -0 N N 0 0 53 nan 53 nan 53 nan 53 +0 N N 0 0 53 nan 53 nan 53 nan 53 +1 N N -0 - 53 +0 53 0x1921FB54442D18p-52 53 nan 53 +inf N N +0 - 53 0 53 +0x1921FB54442D18p-52 53 nan 53 +inf N N 0 0 53 nan 53 nan 53 nan 53 nan N N # pure real argument +- + 53 -0x1E27076E2AF2E6p-57 53 -0x1921FB54442D18p-52 53 -17 53 -0 N N +- - 53 -0x1E27076E2AF2E6p-57 53 +0x1921FB54442D18p-52 53 -17 53 +0 N N ++ + 53 +0x1E27076E2AF2E6p-57 53 -0x1921FB54442D18p-52 53 +17 53 -0 N N ++ - 53 +0x1E27076E2AF2E6p-57 53 +0x1921FB54442D18p-52 53 +17 53 +0 N N ++ 0 53 -0x1F2272AE325A57p-53 53 -0 53 -.75 53 -0 N N ++ 0 53 -0x1F2272AE325A57p-53 53 +0 53 -.75 53 +0 N N +- 0 53 +0x1F2272AE325A57p-53 53 -0 53 +.75 53 -0 N N +- 0 53 +0x1F2272AE325A57p-53 53 +0 53 +.75 53 +0 N N # pure imaginary argument +0 - 53 -0 53 -0x167D8863BC99BDp-52 53 -0 53 -6 N N +0 - 53 +0 53 -0x167D8863BC99BDp-52 53 +0 53 -6 N N +0 + 53 -0 53 +0x167D8863BC99BDp-52 53 -0 53 +6 N N +0 + 53 +0 53 +0x167D8863BC99BDp-52 53 +0 53 +6 N N +0 + 53 -0 53 -0x1F5B75F92C80DDp-55 53 -0 53 -.25 N N +0 + 53 +0 53 -0x1F5B75F92C80DDp-55 53 +0 53 -.25 N N +0 - 53 -0 53 +0x1F5B75F92C80DDp-55 53 -0 53 +.25 N N +0 - 53 +0 53 +0x1F5B75F92C80DDp-55 53 +0 53 +.25 N N # IEEE-754 double precision -- cgit v1.2.1