diff options
Diffstat (limited to 'src/asin.c')
| -rw-r--r-- | src/asin.c | 75 |
1 files changed, 70 insertions, 5 deletions
@@ -1,6 +1,6 @@ /* mpc_asin -- arcsine of a complex number. -Copyright (C) 2009, 2010, 2011, 2012 INRIA +Copyright (C) 2009, 2010, 2011, 2012, 2013 INRIA This file is part of GNU MPC. @@ -20,13 +20,74 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include "mpc-impl.h" +/* Special case op = 1 + i*y for tiny y (see algorithms.tex). + Return 0 if special formula fails, otherwise put in z1 the approximate + value which needs to be converted to rop. + z1 is a temporary variable with enough precision. + */ +static int +mpc_asin_special (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, mpc_ptr z1) +{ + mpfr_exp_t ey = mpfr_get_exp (mpc_imagref (op)); + mpfr_t abs_y; + mpfr_prec_t p; + int inex; + + /* |Re(asin(1+i*y)) - pi/2| <= y^(1/2) */ + if (ey >= 0 || ((-ey) / 2 < mpfr_get_prec (mpc_realref (z1)))) + return 0; + if (mpfr_get_prec (mpc_imagref (z1)) > (-ey) + 2) /* p > -ey + 2 */ + return 0; + + mpfr_const_pi (mpc_realref (z1), MPFR_RNDN); + mpfr_div_2exp (mpc_realref (z1), mpc_realref (z1), 1, MPFR_RNDN); /* exact */ + p = mpfr_get_prec (mpc_realref (z1)); + /* if z1 has precision p, the error on z1 is 1/2*ulp(z1) = 2^(-p) so far, + and since ey <= -2p, then y^(1/2) <= 1/2*ulp(z1) too, thus the total + error is bounded by ulp(z1) */ + if (!mpfr_can_round (mpc_realref(z1), p, MPFR_RNDN, MPFR_RNDZ, + mpfr_get_prec (mpc_realref(rop)) + + (MPC_RND_RE(rnd) == MPFR_RNDN))) + return 0; + + /* |Im(asin(1+i*y)) - y^(1/2)| <= (1/12) * y^(3/2) for y >= 0 (err >= 0) + |Im(asin(1-i*y)) + y^(1/2)| <= (1/12) * y^(3/2) for y >= 0 (err <= 0) */ + abs_y[0] = mpc_imagref (op)[0]; + if (mpfr_signbit (mpc_imagref (op))) + MPFR_CHANGE_SIGN (abs_y); + inex = mpfr_sqrt (mpc_imagref (z1), abs_y, MPFR_RNDN); /* error <= 1/2 ulp */ + if (mpfr_signbit (mpc_imagref (op))) + MPFR_CHANGE_SIGN (mpc_imagref (z1)); + /* If z1 has precision p, the error on z1 is 1/2*ulp(z1) = 2^(-p) so far, + and (1/12) * y^(3/2) <= (1/8) * y * y^(1/2) <= 2^(ey-3)*2^p*ulp(y^(1/2)) + thus for p+ey-3 <= -1 we have (1/12) * y^(3/2) <= (1/2) * ulp(y^(1/2)), + and the total error is bounded by ulp(z1). + Note: if y^(1/2) is exactly representable, and ends with many zeroes, + then mpfr_can_round below will fail; however in that case we know that + Im(asin(1+i*y)) is away from +/-y^(1/2) wrt zero. */ + if (inex == 0) /* enlarge im(z1) so that the inexact flag is correct */ + { + if (mpfr_signbit (mpc_imagref (op))) + mpfr_nextbelow (mpc_imagref (z1)); + else + mpfr_nextabove (mpc_imagref (z1)); + return 1; + } + p = mpfr_get_prec (mpc_imagref (z1)); + if (!mpfr_can_round (mpc_imagref(z1), p - 1, MPFR_RNDA, MPFR_RNDZ, + mpfr_get_prec (mpc_imagref(rop)) + + (MPC_RND_IM(rnd) == MPFR_RNDN))) + return 0; + return 1; +} + int mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { - mpfr_prec_t p, p_re, p_im, incr_p = 0; + mpfr_prec_t p, p_re, p_im; mpfr_rnd_t rnd_re, rnd_im; mpc_t z1; - int inex; + int inex, loop = 0; /* special values */ if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) @@ -148,11 +209,15 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { mpfr_exp_t ex, ey, err; - p += mpc_ceil_log2 (p) + 3 + incr_p; /* incr_p is zero initially */ - incr_p = p / 2; + loop ++; + p += (loop <= 2) ? mpc_ceil_log2 (p) + 3 : p / 2; mpfr_set_prec (mpc_realref(z1), p); mpfr_set_prec (mpc_imagref(z1), p); + /* try special code for 1+i*y with tiny y */ + if (loop == 1 && mpc_asin_special (rop, op, rnd, z1)) + break; + /* z1 <- z^2 */ mpc_sqr (z1, op, MPC_RNDNN); /* err(x) <= 1/2 ulp(x), err(y) <= 1/2 ulp(y) */ |