diff options
author | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2021-02-03 04:29:56 +0000 |
---|---|---|
committer | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2021-02-03 04:29:56 +0000 |
commit | 901eb7ea6b98872dd17f84c4120aeff7ff21978f (patch) | |
tree | 5fc7e6ad68b0a97c6d8989c8b4ce1093e02dca1e /src | |
parent | 604fb667ce910a3ce23f54cbc25048e764322dfe (diff) | |
download | mpfr-901eb7ea6b98872dd17f84c4120aeff7ff21978f.tar.gz |
[src/atan2u.c] deal with underflow and overflow in y/x
git-svn-id: https://scm.gforge.inria.fr/anonscm/svn/mpfr/trunk@14329 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'src')
-rw-r--r-- | src/atan2u.c | 100 |
1 files changed, 97 insertions, 3 deletions
diff --git a/src/atan2u.c b/src/atan2u.c index 1839bb7ed..b424f1c68 100644 --- a/src/atan2u.c +++ b/src/atan2u.c @@ -68,6 +68,91 @@ mpfr_atan2u_aux2 (mpfr_ptr z, unsigned long u, int e, int s, return mpfr_check_range (z, inex, rnd_mode); } +/* return round(sign(s)*(u/2-eps),rnd_mode), where eps < 1/2*ulp(u/2) */ +static int +mpfr_atan2u_aux3 (mpfr_ptr z, unsigned long u, int s, mpfr_rnd_t rnd_mode) +{ + mpfr_t t; + mpfr_prec_t prec; + int inex; + + prec = (MPFR_PREC(z) > 64) ? MPFR_PREC(z) + 2 : 66; + /* prec >= PREC(z)+2 and prec >= 64 + 2 */ + mpfr_init2 (t, prec); + mpfr_set_ui_2exp (t, u, -1, MPFR_RNDN); /* exact */ + mpfr_nextbelow (t); + if (s < 0) + mpfr_neg (t, t, MPFR_RNDN); + inex = mpfr_set (z, t, rnd_mode); + mpfr_clear (t); + return inex; +} + +/* return round(sign(y)*(u/4-sign(x)*eps),rnd_mode), + where eps < 1/2*ulp(u/4) */ +static int +mpfr_atan2u_aux4 (mpfr_ptr z, mpfr_srcptr y, mpfr_srcptr x, unsigned long u, + mpfr_rnd_t rnd_mode) +{ + mpfr_t t; + mpfr_prec_t prec; + int inex; + + prec = (MPFR_PREC(z) > 64) ? MPFR_PREC(z) + 2 : 66; + /* prec >= PREC(z)+2 and prec >= 64 + 2 */ + mpfr_init2 (t, prec); + mpfr_set_ui_2exp (t, u, -2, MPFR_RNDN); /* exact */ + if (MPFR_SIGN(x) > 0) + mpfr_nextbelow (t); + else + mpfr_nextabove (t); + if (MPFR_SIGN(y) < 0) + mpfr_neg (t, t, MPFR_RNDN); + inex = mpfr_set (z, t, rnd_mode); + mpfr_clear (t); + return inex; +} + +/* deal with underflow case, i.e., when y/x rounds to zero */ +static int +mpfr_atan2u_underflow (mpfr_ptr z, mpfr_srcptr y, mpfr_srcptr x, + unsigned long u, mpfr_rnd_t rnd_mode) +{ + mpfr_exp_t e = MPFR_GET_EXP(y) - MPFR_GET_EXP(x) + 63; + /* Detect underflow: since |atan(|y/x|)| < |y/x| for |y/x| < 1, + |atan2u(y,x,u)| < |y/x|*u/(2*pi) < 2^62*|y/x| < 2^(EXP(y)-EXP(x)+63). + For x > 0, we have underflow when EXP(y)-EXP(x)+63 < emin. + For x < 0, we have underflow when EXP(y)-EXP(x)+63 < EXP(u/2)-prec. */ + if (MPFR_IS_POS(x)) + { + MPFR_ASSERTN(e < __gmpfr_emin); + return mpfr_underflow (z, + (rnd_mode == MPFR_RNDN) ? MPFR_RNDZ : rnd_mode, MPFR_SIGN(y)); + } + else if (MPFR_IS_NEG(x)) + { + MPFR_ASSERTN(e < 63 - MPFR_PREC(z)); + return mpfr_atan2u_aux3 (z, u, MPFR_SIGN(y), rnd_mode); + } +} + +/* deal with overflow case, i.e., when y/x rounds to infinity */ +static int +mpfr_atan2u_overflow (mpfr_ptr z, mpfr_srcptr y, mpfr_srcptr x, + unsigned long u, mpfr_rnd_t rnd_mode) +{ + /* When t goes to +Inf, pi/2 - 1/t < atan(t) < pi/2, + thus u/4 - u/(2*pi*t) < atanu(t) < u/4. + As soon as u/(2*pi*t) < 1/2*ulp(u/4), the result is either u/4 + or the number just below. + Here t = y/x, thus 1/t <= x/y < 2^(EXP(x)-EXP(y)+1), + and u/(2*pi*t) < 2^(EXP(x)-EXP(y)+62). */ + mpfr_exp_t e = MPFR_GET_EXP(x) - MPFR_GET_EXP(y) + 62; + mpfr_exp_t ulpz = 62 - MPFR_PREC(z); /* ulp(u/4) <= 2^ulpz */ + MPFR_ASSERTN (e < ulpz - 1); + return mpfr_atan2u_aux4 (z, y, x, u, rnd_mode); +} + /* put in z the correctly rounded value of atan2y(y,x,u) */ int mpfr_atan2u (mpfr_ptr z, mpfr_srcptr y, mpfr_srcptr x, unsigned long u, @@ -206,13 +291,22 @@ mpfr_atan2u (mpfr_ptr z, mpfr_srcptr y, mpfr_srcptr x, unsigned long u, Third quadrant: [-1,-1/2] [-u/2,-u/4] Fourth quadrant: [-1/2,0] [-u/4,0] */ inex = mpfr_div (tmp, y, x, MPFR_RNDN); + if (MPFR_IS_ZERO(tmp)) + { + mpfr_clear (tmp); + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_atan2u_underflow (z, y, x, u, rnd_mode); + } + if (MPFR_IS_INF(tmp)) + { + mpfr_clear (tmp); + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_atan2u_overflow (z, y, x, u, rnd_mode); + } MPFR_SET_SIGN (tmp, 1); - /* FIXME: The division may overflow or underflow, in which case - the following MPFR_GET_EXP is incorrect. */ expt = MPFR_GET_EXP(tmp); /* |tmp - y/x| <= e1 := 1/2*ulp(tmp) = 2^(expt-prec-1) */ inex |= mpfr_atanu (tmp, tmp, u, MPFR_RNDN); - MPFR_ASSERTN(inex != 0); /* the derivative of atan(t) is 1/(1+t^2), thus if tmp2 is the new value of tmp, we have |tmp2 - atan(y/x)| <= 1/2*ulp(tmp2) + e1/(1+tmp^2) */ e = (expt < 1) ? 0 : expt - 1; |