diff options
| -rw-r--r-- | coth.c | 14 | ||||
| -rw-r--r-- | csch.c | 10 | ||||
| -rw-r--r-- | gen_inverse.h | 19 | ||||
| -rw-r--r-- | tests/tcoth.c | 28 | ||||
| -rw-r--r-- | tests/tcsch.c | 16 | ||||
| -rw-r--r-- | tests/tsech.c | 16 |
6 files changed, 97 insertions, 6 deletions
@@ -34,4 +34,18 @@ MA 02110-1301, USA. */ #define ACTION_ZERO(y,x) do { MPFR_SET_SAME_SIGN(y,x); MPFR_SET_ZERO(y); \ MPFR_RET(0); } while (1) +/* We know |coth(x)| > 1, thus if the approximation z is such that + 1 <= z <= 1 + 2^(-p) where p is the target precision, then the + result is either 1 or nextabove(1) = 1 + 2^(1-p). */ +#define ACTION_SPECIAL \ + if (MPFR_GET_EXP(z) == 1) /* 1 <= |z| < 2 */ \ + { \ + mpfr_sub_ui (z, z, MPFR_SIGN(z) > 0 ? 1 : -1, GMP_RNDN); \ + if (MPFR_IS_ZERO(z) || MPFR_GET_EXP(z) <= - (mp_exp_t) precy) \ + { \ + mpfr_add_ui (z, z, MPFR_SIGN(z) > 0 ? 1 : -1, GMP_RNDN); \ + break; \ + } \ + } + #include "gen_inverse.h" @@ -20,11 +20,11 @@ the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ /* the hyperbolic cosecant is defined by csch(x) = 1/sinh(x). - csc (NaN) = NaN. - csc (+Inf) = 0+. - csc (-Inf) = 0-. - csc (+0) = +Inf. - csc (-0) = -Inf. + csch (NaN) = NaN. + csch (+Inf) = +0. + csch (-Inf) = -0. + csch (+0) = +Inf. + csch (-0) = -Inf. */ #define FUNCTION mpfr_csch diff --git a/gen_inverse.h b/gen_inverse.h index d2a69d8e7..b4f23a305 100644 --- a/gen_inverse.h +++ b/gen_inverse.h @@ -22,6 +22,10 @@ MA 02110-1301, USA. */ #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" +#ifndef ACTION_SPECIAL +#define ACTION_SPECIAL +#endif + /* example of use: #define FUNCTION mpfr_sec #define INVERSE mpfr_cos @@ -60,13 +64,26 @@ FUNCTION (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_ZIV_INIT (loop, m); for(;;) { - INVERSE (z, x, GMP_RNDZ); /* error k_u < 1 ulp */ + INVERSE (z, x, GMP_RNDZ); /* error k_u < 1 ulp */ + /* the following assumes that if an overflow happens with + MPFR_EMAX_MAX, then necessarily an underflow happens with + __gmpfr_emin */ + if (mpfr_overflow_p ()) + { + int s = MPFR_SIGN(z); + MPFR_ZIV_FREE (loop); + mpfr_clear (z); + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_underflow (y, (rnd_mode == GMP_RNDN) ? + GMP_RNDZ : rnd_mode, s); + } mpfr_ui_div (z, 1, z, GMP_RNDN); /* the error is less than c_w + 2*c_u*k_u (see algorithms.tex), where c_w = 1/2, c_u = 1 since z was rounded towards zero, thus 1/2 + 2 < 4 */ if (MPFR_LIKELY (MPFR_CAN_ROUND (z, m - 2, precy, rnd_mode))) break; + ACTION_SPECIAL; MPFR_ZIV_NEXT (loop, m); mpfr_set_prec (z, m); } diff --git a/tests/tcoth.c b/tests/tcoth.c index 94605256d..1fef81189 100644 --- a/tests/tcoth.c +++ b/tests/tcoth.c @@ -98,6 +98,34 @@ check_bugs (void) exit (1); } + mpfr_set_prec (x, 53); + mpfr_set_prec (y, 53); + + mpfr_set_str (x, "18.368400284838550", 10, GMP_RNDN); + mpfr_set_str (y, "1.0000000000000002", 10, GMP_RNDN); + mpfr_coth (x, x, GMP_RNDN); + if (mpfr_cmp (x, y) != 0) + { + printf ("Error for coth(18.368400284838550)\n"); + exit (1); + } + + mpfr_set_str (x, "18.714973875118520", 10, GMP_RNDN); + mpfr_coth (x, x, GMP_RNDN); + if (mpfr_cmp (x, y) != 0) + { + printf ("Error for coth(18.714973875118520)\n"); + exit (1); + } + + mpfr_set_str (x, "18.714973875118524", 10, GMP_RNDN); + mpfr_coth (x, x, GMP_RNDN); + if (mpfr_cmp_ui (x, 1) != 0) + { + printf ("Error for coth(18.714973875118524)\n"); + exit (1); + } + mpfr_clear (x); mpfr_clear (y); } diff --git a/tests/tcsch.c b/tests/tcsch.c index 172625488..271666518 100644 --- a/tests/tcsch.c +++ b/tests/tcsch.c @@ -75,6 +75,22 @@ check_specials (void) exit (1); } + /* check huge x */ + mpfr_set_str (x, "8e8", 10, GMP_RNDN); + mpfr_csch (y, x, GMP_RNDN); + if (! (mpfr_zero_p (y) && MPFR_SIGN (y) > 0)) + { + printf ("Error: csch(8e8) != +0\n"); + exit (1); + } + mpfr_set_str (x, "-8e8", 10, GMP_RNDN); + mpfr_csch (y, x, GMP_RNDN); + if (! (mpfr_zero_p (y) && MPFR_SIGN (y) < 0)) + { + printf ("Error: csch(-8e8) != -0\n"); + exit (1); + } + mpfr_clear (x); mpfr_clear (y); } diff --git a/tests/tsech.c b/tests/tsech.c index 3a585763e..623def916 100644 --- a/tests/tsech.c +++ b/tests/tsech.c @@ -75,6 +75,22 @@ check_specials (void) exit (1); } + /* check huge x */ + mpfr_set_str (x, "8e8", 10, GMP_RNDN); + mpfr_sech (y, x, GMP_RNDN); + if (! (mpfr_zero_p (y) && MPFR_SIGN (y) > 0)) + { + printf ("Error: sech(8e8) != +0\n"); + exit (1); + } + mpfr_set_str (x, "-8e8", 10, GMP_RNDN); + mpfr_sech (y, x, GMP_RNDN); + if (! (mpfr_zero_p (y) && MPFR_SIGN (y) > 0)) + { + printf ("Error: sech(-8e8) != +0\n"); + exit (1); + } + mpfr_clear (x); mpfr_clear (y); } |