From 1b66d90c9dc94acbad0427571ed08d1be287149f Mon Sep 17 00:00:00 2001 From: Paul Zimmermann Date: Wed, 19 Feb 2020 17:57:15 +0100 Subject: improved asin (try series for |z| < 1, and special case with tiny Im(z)) --- src/asin.c | 232 +++++++++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 180 insertions(+), 52 deletions(-) diff --git a/src/asin.c b/src/asin.c index e69e2b0..9b09634 100644 --- a/src/asin.c +++ b/src/asin.c @@ -1,6 +1,6 @@ /* mpc_asin -- arcsine of a complex number. -Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014 INRIA +Copyright (C) 2009-2020 INRIA This file is part of GNU MPC. @@ -80,62 +80,182 @@ mpc_asin_special (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, mpc_ptr z1) return 1; } -/* Special case for |z| <= 1/2 and tiny imaginary part (see algorithms.tex). - Return 0 if the special case fails (this cannot be confused with an exact - ternary value since this cannot happen). - Otherwise puts an approximation in z1, so that rounding it in rop will get - the correct result. */ +/* Put in s an approximation of asin(z) using: + asin z = z + 1/2*z^3/3 + (1*3)/(2*4)*z^5/5 + ... + Assume |Re(z)|, |Im(z)| < 1/2. + Return non-zero if we can get the correct result by rounding s: + mpc_set (rop, s, ...) */ static int -mpc_asin_special2 (mpc_ptr rop, mpc_srcptr op, mpc_ptr z1) +mpc_asin_series (mpc_srcptr rop, mpc_ptr s, mpc_srcptr z, mpc_rnd_t rnd) { - /* for |z| <= 1/2, we know that: - |Re(asin z) - Re(z)| < 0.198 |z|^3 - |Im(asin z) - Im(z)| < 0.667 |z|^2 |Im(z)| - thus if 0.667 |z|^2 < 1/2 ulp(1) we can round to Re(z) and Im(z) */ - mpfr_exp_t ex = mpfr_get_exp (mpc_realref (op)); - mpfr_exp_t ey = mpfr_get_exp (mpc_imagref (op)); - mpfr_prec_t px = mpfr_get_prec (mpc_realref (rop)); - mpfr_prec_t py = mpfr_get_prec (mpc_imagref (rop)); - int sign_re, sign_im, sign_re_t, sign_im_t; - mpc_t t; - - /* |x| < 2^ex and |y| < 2^ey thus |z^2| < 2^(2*max(ex,ey)+1) */ - ex = (ex > ey) ? 2 * ex + 1 : 2 * ey + 1; - /* 0.667 |z|^2 < 1/2 ulp(1) if 2^ex <= 2^-p thus ex <= -p */ - if (ex + px > 0 || ex + py > 0) - return 0; /* special case fails */ - - /* if the precision of the input is larger than that of the output, - we fail since we should deal with a double-rounding, which is not - implemented below */ - if (mpfr_get_prec (mpc_realref (op)) > px || - mpfr_get_prec (mpc_imagref (op)) > py) + mpc_t w, t; + unsigned long k, err, kx, ky; + mpfr_prec_t p; + mpfr_exp_t ex, ey, e; + + /* assume z = (x,y) with |x|,|y| < 2^(-e) with e >= 1, see the error + analysis in algorithms.tex */ + ex = mpfr_get_exp (mpc_realref (z)); + MPC_ASSERT(ex <= -1); + ey = mpfr_get_exp (mpc_imagref (z)); + MPC_ASSERT(ey <= -1); + e = (ex >= ey) ? ex : ey; + e = -e; + /* now e >= 1 */ + + p = mpfr_get_prec (mpc_realref (s)); /* working precision */ + MPC_ASSERT(mpfr_get_prec (mpc_imagref (s)) == p); + + mpc_init2 (w, p); + mpc_init2 (t, p); + mpc_set (s, z, MPC_RNDNN); + mpc_sqr (w, z, MPC_RNDNN); + mpc_set (t, z, MPC_RNDNN); + for (k = 1; ;k++) + { + mpfr_exp_t exp_x, exp_y; + mpc_mul (t, t, w, MPC_RNDNN); + mpc_mul_ui (t, t, (2 * k - 1) * (2 * k - 1), MPC_RNDNN); + mpc_div_ui (t, t, (2 * k) * (2 * k + 1), MPC_RNDNN); + exp_x = mpfr_get_exp (mpc_realref (s)); + exp_y = mpfr_get_exp (mpc_imagref (s)); + if (mpfr_get_exp (mpc_realref (t)) < exp_x - p && + mpfr_get_exp (mpc_imagref (t)) < exp_y - p) + /* Re(t) < 1/2 ulp(Re(s)) and Im(t) < 1/2 ulp(Im(s)), + thus adding t to s will not change s */ + break; + mpc_add (s, s, t, MPC_RNDNN); + } + mpc_clear (w); + mpc_clear (t); + /* check (2k-1)^2 is exactly representable */ + MPC_ASSERT(2 * k - 1 <= ULONG_MAX / (2 * k - 1)); + /* maximal absolute error on Re(s),Im(s) is: + (5k-3)k/2*2^(-1-p) for e=1 + 5k/2*2^(-e-p) for e >= 2 */ + if (e == 1) + { + MPC_ASSERT(5 * k - 3 <= ULONG_MAX / k); + kx = (5 * k - 3) * k; + } + else + kx = 5 * k; + kx = (kx + 1) / 2; /* takes into account the 1/2 factor in both cases */ + /* now (5k-3)k/2 <= kx for e=1, and 5k/2 <= kx for e >= 2, thus + the maximal absolute error on Re(s),Im(s) is bounded by kx*2^(-e-p) */ + + e = -e; + ky = kx; + + /* for the real part, convert the maximal absolute error kx*2^(e-p) into + relative error */ + ex = mpfr_get_exp (mpc_realref (s)); + /* ulp(Re(s)) = 2^(ex+1-p) */ + if (ex+1 > e) /* divide kx by 2^(ex+1-e) */ + while (ex+1 > e) + { + kx = (kx + 1) / 2; + ex --; + } + else /* multiply kx by 2^(e-(ex+1)) */ + kx <<= e - (ex+1); + /* now the rounding error is bounded by kx*ulp(Re(s)), add the + mathematical error which is bounded by ulp(Re(s)): the first neglected + term is less than 1/2*ulp(Re(s)), and each term decreases by at least + a factor 2, since |z^2| <= 1/2. */ + kx ++; + for (err = 0; kx > 2; err ++, kx = (kx + 1) / 2); + /* can we round Re(s) with error less than 2^(EXP(Re(s))-err) ? */ + if (!mpfr_can_round (mpc_realref (s), p - err, MPFR_RNDN, MPFR_RNDZ, + mpfr_get_prec (mpc_realref (rop)) + + (MPC_RND_RE(rnd) == MPFR_RNDN))) return 0; - mpc_init2 (t, 32); - mpc_pow_ui (t, op, 3, MPC_RNDNN); - sign_re_t = MPFR_SIGN (mpc_realref (t)); - sign_im_t = MPFR_SIGN (mpc_imagref (t)); - mpc_clear (t); + /* same for the imaginary part */ + ey = mpfr_get_exp (mpc_imagref (s)); + /* we take for e the exponent of Im(z), which amounts to divide the error by + 2^delta where delta is the exponent difference between Re(z) and Im(z) + (see algorithms.tex) */ + e = mpfr_get_exp (mpc_imagref (z)); + /* ulp(Im(s)) = 2^(ey+1-p) */ + if (ey+1 > e) /* divide ky by 2^(ey+1-e) */ + while (ey+1 > e) + { + ky = (ky + 1) / 2; + ey --; + } + else /* multiply ky by 2^(e-(ey+1)) */ + ky <<= e - (ey+1); + /* now the rounding error is bounded by ky*ulp(Im(s)), add the + mathematical error which is bounded by ulp(Im(s)): the first neglected + term is less than 1/2*ulp(Im(s)), and each term decreases by at least + a factor 2, since |z^2| <= 1/2. */ + ky ++; + for (err = 0; ky > 2; err ++, ky = (ky + 1) / 2); + /* can we round Im(s) with error less than 2^(EXP(Im(s))-err) ? */ + return mpfr_can_round (mpc_imagref (s), p - err, MPFR_RNDN, MPFR_RNDZ, + mpfr_get_prec (mpc_imagref (rop)) + + (MPC_RND_IM(rnd) == MPFR_RNDN)); +} - /* approximate the real part: since asin z = z + z^3/6 + ..., the real part - has to be rounded according to the sign of Re(t) */ - mpfr_set (mpc_realref (z1), mpc_realref (op), MPFR_RNDN); - sign_re = MPFR_SIGN (mpc_realref (op)); - if (sign_re * sign_re_t > 0) - MPFR_ADD_ONE_ULP (mpc_realref (z1)); - else - MPFR_SUB_ONE_ULP (mpc_realref (z1)); +/* Put in s an approximation of asin(z) for |Re(z)| < 1 and tiny Im(y) + (see algorithms.tex). Assume z = x + i*y: + |Re(asin(z)) - asin(x)| <= Pi*beta^2/(1-beta) + |Im(asin(z)) - y/sqrt(1-x^2)| <= Pi/2*beta^2/(1-beta) + where beta = |y|/(1-|x|). + We assume |x| >= 1/2 > |y| here, thus beta < 1, and the bounds + simplify to 16y^2. + Assume Re(s) and Im(s) have the same precision. + Return non-zero if we can get the correct result by rounding s: + mpc_set (rop, s, ...) */ +static int +mpc_asin_tiny (mpc_srcptr rop, mpc_ptr s, mpc_srcptr z, mpc_rnd_t rnd) +{ + mpfr_exp_t ey, e1, e2, es; + mpfr_prec_t p = mpfr_get_prec (mpc_realref (s)), err; - /* approximate the imaginary part */ - mpfr_set (mpc_imagref (z1), mpc_imagref (op), MPFR_RNDN); - sign_im = MPFR_SIGN (mpc_imagref (op)); - if (sign_im * sign_im_t > 0) - MPFR_ADD_ONE_ULP (mpc_imagref (z1)); - else - MPFR_SUB_ONE_ULP (mpc_imagref (z1)); + ey = mpfr_get_exp (mpc_imagref (z)); - return 1; + MPC_ASSERT(mpfr_get_exp (mpc_realref (z)) >= 0); /* |x| >= 1/2 */ + MPC_ASSERT(ey < 0); /* |y| < 1/2 */ + + e2 = 2 * ey + 4; + /* for |x| >= 0.5, |asin x| >= 0.5, thus no need to compute asin(x) + if 16y^2 >= 1/2 ulp(1) for the target variable */ + if (e2 >= - (mpfr_exp_t) mpfr_get_prec (mpc_realref (rop))) + return 0; + + /* real part */ + mpfr_asin (mpc_realref (s), mpc_realref (z), MPFR_RNDN); + /* check that we can round with error < 16y^2 */ + e1 = mpfr_get_exp (mpc_realref (s)) - p; + err = (e1 >= e2) ? 0 : e2 - e1; + if (!mpfr_can_round (mpc_realref (s), p - err, MPFR_RNDN, MPFR_RNDZ, + mpfr_get_prec (mpc_realref (rop)) + + (MPC_RND_RE(rnd) == MPFR_RNDN))) + return 0; + + /* now compute the approximate imaginary part y/sqrt(1-x^2) */ + mpfr_sqr (mpc_imagref (s), mpc_realref (z), MPFR_RNDN); + /* now Im(s) approximates x^2, with 1/4 <= Im(s) <= 1 and absolute error + less than 2^-p */ + mpfr_ui_sub (mpc_imagref (s), 1, mpc_imagref (s), MPFR_RNDN); + /* now Im(s) approximates 1-x^2, with 0 <= Im(s) <= 3/4, assuming p >= 2, + and absolute error less than 2^(1-p) */ + mpfr_sqrt (mpc_imagref (s), mpc_imagref (s), MPFR_RNDN); /* sqrt(1-x^2) */ + es = mpfr_get_exp (mpc_imagref (s)); + /* now Im(s) approximates sqrt(1-x^2), with 0 <= Im(s) <= 7/8, + assuming p >= 3, and absolute error less than 2^-p + 2^(1-p)/s + <= 3*2^-p/s, thus relative error less than 3*2^-p/s^2. + Since |s| >= 2^(es-1), the relative error is less than 3*2^(-p+2-2*es) */ + mpfr_div (mpc_imagref (s), mpc_imagref (z), mpc_imagref (s), MPFR_RNDN); + /* now Im(s) approximates y/sqrt(1-x^2), with relative error less than + (3*2^(2-2*es)+2)*2^-p. Since es <= 0, 2-2*es >= 2, thus the relative + error is less than 2^(4-2*es-p) */ + err = 4 - 2 * es; + return mpfr_can_round (mpc_imagref (s), p - err, MPFR_RNDN, MPFR_RNDZ, + mpfr_get_prec (mpc_imagref (rop)) + + (MPC_RND_IM(rnd) == MPFR_RNDN)); } int @@ -281,8 +401,16 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpc_asin_special (rop, op, rnd, z1)) break; - /* try special code for 1+i*y with tiny y */ - if (loop == 1 && mpc_asin_special2 (rop, op, z1)) + /* try special code for small z */ + if (mpfr_get_exp (mpc_realref (op)) <= -1 && + mpfr_get_exp (mpc_imagref (op)) <= -1 && + mpc_asin_series (rop, z1, op, rnd)) + break; + + /* try special code for 1/2 <= |x| < 1 and |y| < 1/2 */ + if (mpfr_get_exp (mpc_realref (op)) == 0 && + mpfr_get_exp (mpc_imagref (op)) <= -1 && + mpc_asin_tiny (rop, z1, op, rnd)) break; /* z1 <- z^2 */ -- cgit v1.2.1