fixed problem of tiny input for coth (and new fix again for csc, which

takes into account the sign of the input)


git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@4493 280ebfd0-de03-0410-8827-d642c229c3f4

1 files changed, 14 insertions, 6 deletions

diff --git a/csc.c b/csc.c
index 246a12aa0..d0f438e19 100644
--- a/csc.c
+++ b/csc.c
@@ -35,8 +35,8 @@ MA 02110-1301, USA. */
                               MPFR_RET(0); } while (1)
 /* near x=0, we have csc(x) = 1/x + x/6 + ..., more precisely we have
    |csc(x) - 1/x| <= 0.2 for |x| <= 1. The analysis is similar to that for
-   gamma(x) near x=0 (see gamma.c):
-   then:
+   gamma(x) near x=0 (see gamma.c), except here the error term has the same
+   sign as 1/x, thus |csc(x)| >= |1/x|. Then:
    (i) either x is a power of two, then 1/x is exactly representable, and
        as long as 1/2*ulp(1/x) > 0.2, we can conclude;
    (ii) otherwise assume x has <= n bits, and y has <= n+1 bits, then
@@ -48,16 +48,24 @@ MA 02110-1301, USA. */
 #define ACTION_TINY(y,x,r) \
   if (MPFR_EXP(x) <= -2 * (mp_exp_t) MAX(MPFR_PREC(x), MPFR_PREC(y)))   \
     {                                                                   \
+      int signx = MPFR_SIGN(x);                                         \
       inexact = mpfr_ui_div (y, 1, x, r);                               \
       if (inexact == 0) /* x is a power of two */                       \
-        {                                                               \
+        { /* result always 1/x, except when rounding away from zero */  \
           if (rnd_mode == GMP_RNDU)                                     \
             {                                                           \
-              mpfr_nextabove (y); /* 2^k + epsilon */                   \
+              if (signx > 0)                                            \
+                mpfr_nextabove (y); /* 2^k + epsilon */                 \
               inexact = 1;                                              \
             }                                                           \
-          else /* round to zero, nearest or to -Inf */                  \
-            inexact = -1; /* 2^k */                                     \
+          else if (rnd_mode == GMP_RNDD)                                \
+            {                                                           \
+              if (signx < 0)                                            \
+                mpfr_nextbelow (y); /* -2^k - epsilon */                \
+              inexact = -1;                                             \
+            }                                                           \
+          else /* round to zero, or nearest */                          \
+            inexact = -signx;                                           \
         }                                                               \
       goto end;                                                         \
     }