mpfr-site.patch

   1 diff -Naurd mpfr-2.2.0/lngamma.c mpfr-2.2.0-p1/lngamma.c
   2 --- mpfr-2.2.0/lngamma.c        2005-09-09 15:17:58.000000000 +0000
   3 +++ mpfr-2.2.0-p1/lngamma.c     2005-09-29 11:27:04.000000000 +0000
   4 @@ -167,8 +167,8 @@
   5    compared = mpfr_cmp_ui (z0, 1);
   6
   7  #ifndef IS_GAMMA
   8 -  if (compared == 0) /* lngamma(1) = +0 */
   9 -    return mpfr_set_ui (y, 0, GMP_RNDN);
  10 +  if (compared == 0 || (compared > 0 && mpfr_cmp_ui (z0, 2) == 0))
  11 +    return mpfr_set_ui (y, 0, GMP_RNDN);  /* lngamma(1 or 2) = +0 */
  12  #endif
  13
  14    mpfr_init2 (s, MPFR_PREC_MIN);
  15 diff -Naurd mpfr-2.2.0/tests/tlngamma.c mpfr-2.2.0-p1/tests/tlngamma.c
  16 --- mpfr-2.2.0/tests/tlngamma.c 2005-09-09 15:17:59.000000000 +0000
  17 +++ mpfr-2.2.0-p1/tests/tlngamma.c      2005-09-29 11:20:34.000000000 +0000
  18 @@ -79,7 +79,7 @@
  19
  20    mpfr_set_ui (x, 1, GMP_RNDN);
  21    mpfr_lngamma (y, x, GMP_RNDN);
  22 -  if (mpfr_cmp_ui (y, 0))
  23 +  if (mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y))
  24      {
  25        printf ("Error for lngamma(1)\n");
  26        exit (1);
  27 @@ -93,6 +93,14 @@
  28        exit (1);
  29      }
  30
  31 +  mpfr_set_ui (x, 2, GMP_RNDN);
  32 +  mpfr_lngamma (y, x, GMP_RNDN);
  33 +  if (mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y))
  34 +    {
  35 +      printf ("Error for lngamma(2)\n");
  36 +      exit (1);
  37 +    }
  38 +
  39    mpfr_set_prec (x, 53);
  40    mpfr_set_prec (y, 53);
  41
  42 diff -Naurd mpfr-2.2.0-p1/mpfr.h mpfr-2.2.0-p2/mpfr.h
  43 --- mpfr-2.2.0-p1/mpfr.h        2005-09-06 15:02:12.000000000 +0000
  44 +++ mpfr-2.2.0-p2/mpfr.h        2005-09-29 11:36:36.000000000 +0000
  45 @@ -630,12 +630,17 @@
  46   (__builtin_constant_p (_s) && (_s) >= 0 ? \
  47     mpfr_cmp_ui ((_f), (_s)) :              \
  48     mpfr_cmp_si_2exp ((_f), (_s), 0))
  49 +#if __GNUC__ > 2 || __GNUC_MINOR__ >= 95
  50  #undef mpfr_set_ui
  51  #define mpfr_set_ui(_f,_u,_r)              \
  52   (__builtin_constant_p (_u) && (_u) == 0 ? \
  53 -   ((_f)->_mpfr_sign = 1,                  \
  54 -    (_f)->_mpfr_exp = __MPFR_EXP_ZERO, 0): \
  55 -    mpfr_set_ui (_f,_u,_r))
  56 +   __extension__ ({                        \
  57 +     mpfr_ptr _p = (_f);                   \
  58 +     _p->_mpfr_sign = 1;                   \
  59 +     _p->_mpfr_exp = __MPFR_EXP_ZERO;      \
  60 +     (void) (_r); 0; }) :                  \
  61 +   mpfr_set_ui (_f,_u,_r))
  62 +#endif
  63  #undef mpfr_set_si
  64  #define mpfr_set_si(_f,_s,_r)              \
  65   (__builtin_constant_p (_s) && (_s) >= 0 ? \
  66 diff -Naurd mpfr-2.2.0-p1/tests/tset_si.c mpfr-2.2.0-p2/tests/tset_si.c
  67 --- mpfr-2.2.0-p1/tests/tset_si.c       2005-08-18 17:03:17.000000000 +0000
  68 +++ mpfr-2.2.0-p2/tests/tset_si.c       2005-09-29 09:19:39.000000000 +0000
  69 @@ -72,6 +72,35 @@
  70    mpfr_clear (x);
  71  }
  72
  73 +static void
  74 +test_macros (void)
  75 +{
  76 +  mpfr_t x[3];
  77 +  mpfr_ptr p;
  78 +  mpfr_rnd_t r;
  79 +
  80 +  mpfr_inits (x[0], x[1], x[2], NULL);
  81 +  p = x[0];
  82 +  r = 0;
  83 +  mpfr_set_ui (p++, 0, r++);
  84 +  if (p != x[1] || r != 1)
  85 +    {
  86 +      printf ("Error in mpfr_set_ui macro: p - x[0] = %d (expecting 1), "
  87 +              "r = %d (expecting 1)\n", (int) (p - x[0]), r);
  88 +      exit (1);
  89 +    }
  90 +  p = x[0];
  91 +  r = 0;
  92 +  mpfr_set_si (p++, 0, r++);
  93 +  if (p != x[1] || r != 1)
  94 +    {
  95 +      printf ("Error in mpfr_set_si macro: p - x[0] = %d (expecting 1), "
  96 +              "r = %d (expecting 1)\n", (int) (p - x[0]), r);
  97 +      exit (1);
  98 +    }
  99 +  mpfr_clears (x[0], x[1], x[2], NULL);
 100 +}
 101 +
 102  /* FIXME: Comparing against mpfr_get_si/ui is not ideal, it'd be better to
 103     have all tests examine the bits in mpfr_t for what should come out.  */
 104
 105 @@ -324,6 +353,7 @@
 106    mpfr_clear (x);
 107
 108    test_2exp ();
 109 +  test_macros ();
 110    tests_end_mpfr ();
 111    return 0;
 112  }
 113 diff -Naurd mpfr-2.2.0-p2/configure mpfr-2.2.0-p3/configure
 114 --- mpfr-2.2.0-p2/configure     2005-09-19 13:31:58.000000000 +0000
 115 +++ mpfr-2.2.0-p3/configure     2005-10-02 10:49:55.000000000 +0000
 116 @@ -10519,7 +10519,7 @@
 117        ;;
 118
 119      darwin* | rhapsody*)
 120 -    if test "$GXX" = yes ; then
 121 +    if test "$GCC" = yes ; then
 122        archive_cmds_need_lc=no
 123        case "$host_os" in
 124        rhapsody* | darwin1.[012])
 125 diff -Naurd mpfr-2.2.0-p3/tests/tpow.c mpfr-2.2.0-p4/tests/tpow.c
 126 --- mpfr-2.2.0-p3/tests/tpow.c  2005-06-02 16:12:05.000000000 +0000
 127 +++ mpfr-2.2.0-p4/tests/tpow.c  2005-10-06 09:54:52.000000000 +0000
 128 @@ -509,6 +509,7 @@
 129    for (i = 0; i < 11; i++)
 130      for (j = 0; j < 11; j++)
 131        {
 132 +        double d;
 133          int p;
 134          static int q[11][11] = {
 135            /*          NaN +inf -inf  +0   -0   +1   -1   +2   -2  +0.5 -0.5 */
 136 @@ -527,7 +528,7 @@
 137          test_pow (r, t[i], t[j], GMP_RNDN);
 138          p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 :
 139            mpfr_cmp_ui (r, 0) == 0 ? 2 :
 140 -          (int) (fabs (mpfr_get_d (r, GMP_RNDN)) * 128.0);
 141 +          (d = mpfr_get_d (r, GMP_RNDN), (int) (ABS(d) * 128.0));
 142          if (p != 0 && MPFR_SIGN(r) < 0)
 143            p = -p;
 144          if (p != q[i][j])
 145 diff -Naurd mpfr-2.2.0-p4/cache.c mpfr-2.2.0-p5/cache.c
 146 --- mpfr-2.2.0-p4/cache.c       2005-08-18 16:35:03.000000000 +0000
 147 +++ mpfr-2.2.0-p5/cache.c       2005-11-23 09:04:29.000000000 +0000
 148 @@ -32,9 +32,9 @@
 149  void
 150  mpfr_clear_cache (mpfr_cache_t cache)
 151  {
 152 -  if (MPFR_PREC(cache->x) != 0)
 153 +  if (MPFR_PREC (cache->x) != 0)
 154      mpfr_clear (cache->x);
 155 -  MPFR_PREC(cache->x) = 0;
 156 +  MPFR_PREC (cache->x) = 0;
 157  }
 158
 159  int
 160 @@ -47,45 +47,58 @@
 161
 162    MPFR_SAVE_EXPO_MARK (expo);
 163
 164 -  /* Check if the cache has been already filled */
 165 -  if (MPFR_UNLIKELY(pold == 0))
 166 -    mpfr_init2 (cache->x, prec);
 167 +  if (MPFR_UNLIKELY (prec > pold))
 168 +    {
 169 +      /* No previous result in the cache or the precision of the
 170 +         previous result is not sufficient. */
 171
 172 -  /* Check if we can round with the previous result */
 173 -  else if (MPFR_LIKELY(prec <= pold))
 174 -    goto round;
 175 +      if (MPFR_UNLIKELY (pold == 0))  /* No previous result. */
 176 +        mpfr_init2 (cache->x, prec);
 177
 178 -  /* Update the cache */
 179 -  pold = prec /*MPFR_PREC_MIN + prec + __gmpfr_ceil_exp2 (prec)*/;
 180 -  MPFR_ASSERTD (pold >= prec);
 181 -  mpfr_prec_round (cache->x, pold, GMP_RNDN);
 182 -  cache->inexact = (*cache->func) (cache->x, GMP_RNDN);
 183 +      /* Update the cache. */
 184 +      pold = prec;
 185 +      mpfr_prec_round (cache->x, pold, GMP_RNDN);
 186 +      cache->inexact = (*cache->func) (cache->x, GMP_RNDN);
 187 +    }
 188
 189 - round:
 190 -  /* First check if the cache has the exact value (Unlikely)
 191 -     Else the exact value is between (assuming x=cache->x > 0)
 192 -     x and x+ulp(x) if cache->inexact < 0
 193 -     x-ulp(x) and x if cache->inexact > 0
 194 -     and abs(x-exact) <= ulp(x)/2 */
 195 -  MPFR_ASSERTD (MPFR_IS_POS(cache->x)); /* TODO...*/
 196 -  /* We must use nextbelow instead of sub_one_ulp, since we know
 197 -     that the exact value is < 1/2ulp(x) (We want sub_demi_ulp(x)).
 198 -     Can't use mpfr_set since we need the even flag. */
 199 +  /* First, check if the cache has the exact value (unlikely).
 200 +     Else the exact value is between (assuming x=cache->x > 0):
 201 +       x and x+ulp(x) if cache->inexact < 0,
 202 +       x-ulp(x) and x if cache->inexact > 0,
 203 +     and abs(x-exact) <= ulp(x)/2. */
 204 +  MPFR_ASSERTN (MPFR_IS_POS (cache->x)); /* TODO... */
 205    sign = MPFR_SIGN (cache->x);
 206    MPFR_SET_EXP (dest, MPFR_GET_EXP (cache->x));
 207    MPFR_SET_SIGN (dest, sign);
 208 -  MPFR_RNDRAW_EVEN (inexact, dest,
 209 -                    MPFR_MANT (cache->x), MPFR_PREC (cache->x), rnd, sign,
 210 -                    if (MPFR_UNLIKELY ( ++MPFR_EXP (dest) > __gmpfr_emax))
 211 -                       mpfr_overflow (dest, rnd, sign) );
 212 -  /* inexact = mpfr_set (dest, cache->x, rnd); */
 213 -  if (MPFR_LIKELY(cache->inexact != 0))
 214 +  MPFR_RNDRAW_GEN (inexact, dest,
 215 +                   MPFR_MANT (cache->x), MPFR_PREC (cache->x), rnd, sign,
 216 +                   if (MPFR_UNLIKELY (cache->inexact == 0))
 217 +                     {
 218 +                       if ((sp[0] & ulp) == 0)
 219 +                         {
 220 +                           inexact = -sign;
 221 +                           goto trunc_doit;
 222 +                         }
 223 +                       else
 224 +                         goto addoneulp;
 225 +                     }
 226 +                   else if (cache->inexact < 0)
 227 +                     goto addoneulp;
 228 +                   else
 229 +                     {
 230 +                       inexact = -sign;
 231 +                       goto trunc_doit;
 232 +                     },
 233 +                   if (MPFR_UNLIKELY (++MPFR_EXP (dest) > __gmpfr_emax))
 234 +                     mpfr_overflow (dest, rnd, sign);
 235 +                  );
 236 +  if (MPFR_LIKELY (cache->inexact != 0))
 237      {
 238        switch (rnd)
 239          {
 240          case GMP_RNDZ:
 241          case GMP_RNDD:
 242 -          if (MPFR_UNLIKELY(inexact == 0))
 243 +          if (MPFR_UNLIKELY (inexact == 0))
 244              {
 245                inexact = cache->inexact;
 246                if (inexact > 0)
 247 @@ -93,7 +106,7 @@
 248              }
 249            break;
 250          case GMP_RNDU:
 251 -          if (MPFR_UNLIKELY(inexact == 0))
 252 +          if (MPFR_UNLIKELY (inexact == 0))
 253              {
 254                inexact = cache->inexact;
 255                if (inexact < 0)
 256 @@ -101,16 +114,7 @@
 257              }
 258            break;
 259          default: /* GMP_RNDN */
 260 -          if (MPFR_UNLIKELY(inexact == MPFR_EVEN_INEX ||
 261 -                            inexact == -MPFR_EVEN_INEX))
 262 -            {
 263 -              if (cache->inexact < 0)
 264 -                mpfr_nextabove (dest);
 265 -              else
 266 -                mpfr_nextbelow (dest);
 267 -              inexact = -inexact;
 268 -            }
 269 -          else if (MPFR_UNLIKELY(inexact == 0))
 270 +          if (MPFR_UNLIKELY(inexact == 0))
 271              inexact = cache->inexact;
 272            break;
 273          }
 274 diff -Naurd mpfr-2.2.0-p4/hypot.c mpfr-2.2.0-p5/hypot.c
 275 --- mpfr-2.2.0-p4/hypot.c       2005-06-06 13:39:39.000000000 +0000
 276 +++ mpfr-2.2.0-p5/hypot.c       2005-11-23 09:04:29.000000000 +0000
 277 @@ -73,6 +73,7 @@
 278    Ey = MPFR_GET_EXP (y);
 279    diff_exp = (mp_exp_unsigned_t) Ex - Ey;
 280
 281 +  Nx = MPFR_PREC (x);   /* Precision of input variable */
 282    Nz = MPFR_PREC (z);   /* Precision of output variable */
 283
 284    /* we have x < 2^Ex thus x^2 < 2^(2*Ex),
 285 @@ -87,10 +88,10 @@
 286       if 2*diff_exp > Nx (see above as if Nz = Nx), therefore on Nz bits.
 287       Hence the condition: 2*diff_exp > MAX(Nz,Nx).
 288    */
 289 -  if (diff_exp > MAX (Nz, MPFR_PREC (x)) / 2)
 290 +  if (diff_exp > MAX (Nz, Nx) / 2)
 291      /* result is |x| or |x|+ulp(|x|,Nz) */
 292      {
 293 -      if (rnd_mode == GMP_RNDU)
 294 +      if (MPFR_UNLIKELY (rnd_mode == GMP_RNDU))
 295          {
 296            /* If z > abs(x), then it was already rounded up; otherwise
 297               z = abs(x), and we need to add one ulp due to y. */
 298 @@ -100,14 +101,27 @@
 299          }
 300        else /* GMP_RNDZ, GMP_RNDD, GMP_RNDN */
 301          {
 302 -          inexact = mpfr_abs (z, x, rnd_mode);
 303 -          return (inexact) ? inexact : -1;
 304 +          if (MPFR_LIKELY (Nz >= Nx))
 305 +            {
 306 +              mpfr_abs (z, x, rnd_mode);  /* exact */
 307 +              return -1;
 308 +            }
 309 +          else
 310 +            {
 311 +              MPFR_SET_EXP (z, Ex);
 312 +              MPFR_SET_SIGN (z, 1);
 313 +              MPFR_RNDRAW_GEN (inexact, z, MPFR_MANT (x), Nx, rnd_mode, 1,
 314 +                               goto addoneulp,
 315 +                  if (MPFR_UNLIKELY (++MPFR_EXP (z) > __gmpfr_emax))
 316 +                    return mpfr_overflow (z, rnd_mode, 1);
 317 +                              );
 318 +              return inexact ? inexact : -1;
 319 +            }
 320          }
 321      }
 322
 323    /* General case */
 324
 325 -  Nx = MPFR_PREC(x);   /* Precision of input variable */
 326    Ny = MPFR_PREC(y);   /* Precision of input variable */
 327
 328    /* compute the working precision -- see algorithms.ps */
 329 diff -Naurd mpfr-2.2.0-p4/mpfr-impl.h mpfr-2.2.0-p5/mpfr-impl.h
 330 --- mpfr-2.2.0-p4/mpfr-impl.h   2005-09-11 22:13:24.000000000 +0000
 331 +++ mpfr-2.2.0-p5/mpfr-impl.h   2005-11-23 09:04:29.000000000 +0000
 332 @@ -184,6 +184,14 @@
 333  # define MPFR_THREAD_ATTR
 334  #endif
 335
 336 +/* Cache struct */
 337 +struct __gmpfr_cache_s {
 338 +  mpfr_t x;
 339 +  int inexact;
 340 +  int (*func)(mpfr_ptr, mpfr_rnd_t);
 341 +};
 342 +typedef struct __gmpfr_cache_s mpfr_cache_t[1];
 343 +
 344  #if defined (__cplusplus)
 345  extern "C" {
 346  #endif
 347 @@ -924,113 +932,19 @@
 348   ******************************************************/
 349
 350  /*
 351 - * Round Mantissa (`srcp`, `sprec`) to mpfr_t `dest` using rounding mode `rnd`
 352 - * assuming dest's sign is `sign`.
 353 - * Execute OVERFLOW_HANDLER in case of overflow when rounding (Power 2 case)
 354 + * Note: due to the labels, one cannot use a macro MPFR_RNDRAW* more than
 355 + * once in a function (otherwise these labels would not be unique).
 356   */
 357 -#define MPFR_RNDRAW(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER)\
 358 -  do {                                                                      \
 359 -    mp_size_t dests, srcs;                                                  \
 360 -    mp_limb_t *destp;                                                       \
 361 -    mp_prec_t destprec, srcprec;                                            \
 362 -                                                                            \
 363 -    /* Check Trivial Case when Dest Mantissa has more bits than source */   \
 364 -    srcprec = sprec;                                                        \
 365 -    destprec = MPFR_PREC (dest);                                            \
 366 -    destp = MPFR_MANT (dest);                                               \
 367 -    if (MPFR_UNLIKELY (destprec >= srcprec))                                \
 368 -      {                                                                     \
 369 -        srcs  = (srcprec  + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB;           \
 370 -        dests = (destprec + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB - srcs;    \
 371 -        MPN_COPY (destp + dests, srcp, srcs);                               \
 372 -        MPN_ZERO (destp, dests);                                            \
 373 -        inexact = 0;                                                        \
 374 -      }                                                                     \
 375 -    else                                                                    \
 376 -      {                                                                     \
 377 -        /* Non trivial case: rounding needed */                             \
 378 -        mp_prec_t sh;                                                       \
 379 -        mp_limb_t *sp;                                                      \
 380 -        mp_limb_t rb, sb, ulp;                                              \
 381 -                                                                            \
 382 -        /* Compute Position and shift */                                    \
 383 -        srcs  = (srcprec  + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB;           \
 384 -        dests = (destprec + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB;           \
 385 -        MPFR_UNSIGNED_MINUS_MODULO (sh, destprec);                          \
 386 -        sp = srcp + srcs - dests;                                           \
 387 -                                                                            \
 388 -        /* General case when prec % BITS_PER_MP_LIMB != 0 */                \
 389 -        if (MPFR_LIKELY (sh != 0))                                          \
 390 -          {                                                                 \
 391 -            mp_limb_t mask;                                                 \
 392 -            /* Compute Rounding Bit and Sticky Bit */                       \
 393 -            mask = MPFR_LIMB_ONE << (sh-1);                                 \
 394 -            rb = sp[0] & mask;                                              \
 395 -            sb = sp[0] & (mask-1);                                          \
 396 -            if (MPFR_UNLIKELY (sb == 0))                                    \
 397 -              { /* TODO: Improve it */                                      \
 398 -                mp_limb_t *tmp;                                             \
 399 -                mp_size_t n;                                                \
 400 -                for (tmp = sp, n = srcs - dests ; n != 0 && sb == 0 ; n--)  \
 401 -                  sb = *--tmp;                                              \
 402 -              }                                                             \
 403 -            ulp = 2*mask;                                                   \
 404 -          }                                                                 \
 405 -        else /* sh == 0 */                                                  \
 406 -          {                                                                 \
 407 -            MPFR_ASSERTD (dests < srcs);                                    \
 408 -            /* Compute Rounding Bit and Sticky Bit */                       \
 409 -            rb = sp[-1] & MPFR_LIMB_HIGHBIT;                                \
 410 -            sb = sp[-1] & (MPFR_LIMB_HIGHBIT-1);                            \
 411 -            if (MPFR_UNLIKELY (sb == 0))                                    \
 412 -              {                                                             \
 413 -                mp_limb_t *tmp;                                             \
 414 -                mp_size_t n;                                                \
 415 -                for (tmp = sp-1, n = srcs - dests-1 ; n!=0 && sb==0 ; n--)  \
 416 -                  sb = *--tmp;                                              \
 417 -              }                                                             \
 418 -            ulp = MPFR_LIMB_ONE;                                            \
 419 -          }                                                                 \
 420 -        /* Rounding */                                                      \
 421 -        if (MPFR_LIKELY (rnd == GMP_RNDN))                                  \
 422 -          {                                                                 \
 423 -            if (rb == 0 || MPFR_UNLIKELY (sb == 0 && (sp[0] & ulp) == 0))   \
 424 -              {                                                             \
 425 -              trunc:                                                        \
 426 -                inexact = MPFR_LIKELY ((sb | rb) != 0) ? -sign : 0;         \
 427 -                MPN_COPY (destp, sp, dests);                                \
 428 -                destp[0] &= ~(ulp-1);                                       \
 429 -              }                                                             \
 430 -            else                                                            \
 431 -              {                                                             \
 432 -              addoneulp:                                                    \
 433 -                if (MPFR_UNLIKELY (mpn_add_1 (destp, sp, dests, ulp)))      \
 434 -                  {                                                         \
 435 -                    destp[dests-1] = MPFR_LIMB_HIGHBIT;                     \
 436 -                    OVERFLOW_HANDLER;                                       \
 437 -                  }                                                         \
 438 -                destp[0] &= ~(ulp-1);                                       \
 439 -                inexact = sign;                                             \
 440 -              }                                                             \
 441 -          }                                                                 \
 442 -        else                                                                \
 443 -          { /* Not Rounding to Nearest */                                   \
 444 -            if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN (sign)))\
 445 -                || MPFR_UNLIKELY ((sb | rb) == 0))                          \
 446 -              goto trunc;                                                   \
 447 -             else                                                           \
 448 -              goto addoneulp;                                               \
 449 -          }                                                                 \
 450 -      }                                                                     \
 451 -  } while (0)
 452
 453  /*
 454 - * Round Mantissa (`srcp`, `sprec`) to mpfr_t `dest` using rounding mode `rnd`
 455 - * assuming dest's sign is `sign`.
 456 - * Execute OVERFLOW_HANDLER in case of overflow when rounding (Power 2 case)
 457 - * Return MPFR_EVEN_INEX in case of EVEN rounding
 458 + * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
 459 + * assuming dest's sign is sign.
 460 + * In rounding to nearest mode, execute MIDDLE_HANDLER when the value
 461 + * is the middle of two consecutive numbers in dest precision.
 462 + * Execute OVERFLOW_HANDLER in case of overflow when rounding.
 463   */
 464 -#define MPFR_RNDRAW_EVEN(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER)\
 465 +#define MPFR_RNDRAW_GEN(inexact, dest, srcp, sprec, rnd, sign,              \
 466 +                        MIDDLE_HANDLER, OVERFLOW_HANDLER)                   \
 467    do {                                                                      \
 468      mp_size_t dests, srcs;                                                  \
 469      mp_limb_t *destp;                                                       \
 470 @@ -1105,19 +1019,8 @@
 471                  destp[0] &= ~(ulp-1);                                       \
 472                }                                                             \
 473              else if (MPFR_UNLIKELY (sb == 0))                               \
 474 -              {                                                             \
 475 -                /* EVEN rounding */                                         \
 476 -                if ((sp[0] & ulp) == 0)                                     \
 477 -                 {                                                          \
 478 -                  MPFR_ASSERTD (rb != 0);                                   \
 479 -                  inexact = -MPFR_EVEN_INEX*sign;                           \
 480 -                  goto trunc_doit;                                          \
 481 -                 }                                                          \
 482 -                else                                                        \
 483 -                 {                                                          \
 484 -                  inexact = MPFR_EVEN_INEX*sign;                            \
 485 -                  goto addoneulp_doit;                                      \
 486 -                 }                                                          \
 487 +              { /* Middle of two consecutive representable numbers */       \
 488 +                MIDDLE_HANDLER;                                             \
 489                }                                                             \
 490              else                                                            \
 491                {                                                             \
 492 @@ -1133,16 +1036,58 @@
 493                }                                                             \
 494            }                                                                 \
 495          else                                                                \
 496 -          { /* Not Rounding to Nearest */                                   \
 497 -            if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN (sign)))\
 498 -                || MPFR_UNLIKELY ((sb | rb) == 0))                          \
 499 +          { /* Directed rounding mode */                                    \
 500 +            if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd,                        \
 501 +                                                MPFR_IS_NEG_SIGN (sign))))  \
 502                goto trunc;                                                   \
 503 +             else if (MPFR_UNLIKELY ((sb | rb) == 0))                       \
 504 +               {                                                            \
 505 +                 inexact = 0;                                               \
 506 +                 goto trunc_doit;                                           \
 507 +               }                                                            \
 508               else                                                           \
 509                goto addoneulp;                                               \
 510            }                                                                 \
 511        }                                                                     \
 512    } while (0)
 513
 514 +/*
 515 + * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
 516 + * assuming dest's sign is sign.
 517 + * Execute OVERFLOW_HANDLER in case of overflow when rounding.
 518 + */
 519 +#define MPFR_RNDRAW(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER) \
 520 +   MPFR_RNDRAW_GEN (inexact, dest, srcp, sprec, rnd, sign,                   \
 521 +     if ((sp[0] & ulp) == 0)                                                 \
 522 +       {                                                                     \
 523 +         inexact = -sign;                                                    \
 524 +         goto trunc_doit;                                                    \
 525 +       }                                                                     \
 526 +     else                                                                    \
 527 +       goto addoneulp;                                                       \
 528 +     , OVERFLOW_HANDLER)
 529 +
 530 +/*
 531 + * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
 532 + * assuming dest's sign is sign.
 533 + * Execute OVERFLOW_HANDLER in case of overflow when rounding.
 534 + * Set inexact to +/- MPFR_EVEN_INEX in case of even rounding.
 535 + */
 536 +#define MPFR_RNDRAW_EVEN(inexact, dest, srcp, sprec, rnd, sign, \
 537 +                         OVERFLOW_HANDLER)                      \
 538 +   MPFR_RNDRAW_GEN (inexact, dest, srcp, sprec, rnd, sign,      \
 539 +     if ((sp[0] & ulp) == 0)                                    \
 540 +       {                                                        \
 541 +         inexact = -MPFR_EVEN_INEX * sign;                      \
 542 +         goto trunc_doit;                                       \
 543 +       }                                                        \
 544 +     else                                                       \
 545 +       {                                                        \
 546 +         inexact = MPFR_EVEN_INEX * sign;                       \
 547 +         goto addoneulp_doit;                                   \
 548 +       }                                                        \
 549 +     , OVERFLOW_HANDLER)
 550 +
 551  /* Return TRUE if b is non singular and we can round it to precision 'prec'
 552     with rounding mode 'rnd', and with error at most 'error' */
 553  #define MPFR_CAN_ROUND(b,err,prec,rnd)                                       \
 554 @@ -1498,6 +1443,13 @@
 555  __MPFR_DECLSPEC int mpfr_const_log2_internal _MPFR_PROTO((mpfr_ptr,mp_rnd_t));
 556  __MPFR_DECLSPEC int mpfr_const_euler_internal _MPFR_PROTO((mpfr_ptr, mp_rnd_t));
 557  __MPFR_DECLSPEC int mpfr_const_catalan_internal _MPFR_PROTO((mpfr_ptr, mp_rnd_t));
 558 +
 559 +__MPFR_DECLSPEC void mpfr_init_cache _MPFR_PROTO ((mpfr_cache_t,
 560 +                                           int(*)(mpfr_ptr,mpfr_rnd_t)));
 561 +__MPFR_DECLSPEC void mpfr_clear_cache _MPFR_PROTO ((mpfr_cache_t));
 562 +__MPFR_DECLSPEC int  mpfr_cache _MPFR_PROTO ((mpfr_ptr, mpfr_cache_t,
 563 +                                              mpfr_rnd_t));
 564 +
 565  __MPFR_DECLSPEC void mpfr_mulhigh_n _MPFR_PROTO ((mp_ptr, mp_srcptr,
 566                                                    mp_srcptr, mp_size_t));
 567
 568 diff -Naurd mpfr-2.2.0-p4/mpfr.h mpfr-2.2.0-p5/mpfr.h
 569 --- mpfr-2.2.0-p4/mpfr.h        2005-09-29 11:36:36.000000000 +0000
 570 +++ mpfr-2.2.0-p5/mpfr.h        2005-11-23 09:04:29.000000000 +0000
 571 @@ -118,14 +118,6 @@
 572  /* For those who needs a direct access and fast access to the sign field */
 573  #define MPFR_SIGN(x) ((x)->_mpfr_sign)
 574
 575 -/* Cache struct */
 576 -struct __gmpfr_cache_s {
 577 -  mpfr_t x;
 578 -  int inexact;
 579 -  int (*func)(mpfr_ptr, mpfr_rnd_t);
 580 -};
 581 -typedef struct __gmpfr_cache_s mpfr_cache_t[1];
 582 -
 583  /* Stack interface */
 584  typedef enum {
 585    MPFR_NAN_KIND = 0,
 586 @@ -544,11 +536,6 @@
 587  __MPFR_DECLSPEC int mpfr_sum _MPFR_PROTO ((mpfr_ptr, mpfr_ptr *__gmp_const,
 588                                             unsigned long, mpfr_rnd_t));
 589
 590 -__MPFR_DECLSPEC void mpfr_init_cache _MPFR_PROTO ((mpfr_cache_t,
 591 -                                           int(*)(mpfr_ptr,mpfr_rnd_t)));
 592 -__MPFR_DECLSPEC void mpfr_clear_cache _MPFR_PROTO ((mpfr_cache_t));
 593 -__MPFR_DECLSPEC int  mpfr_cache _MPFR_PROTO ((mpfr_ptr, mpfr_cache_t,
 594 -                                              mpfr_rnd_t));
 595  __MPFR_DECLSPEC void mpfr_free_cache _MPFR_PROTO ((void));
 596
 597  __MPFR_DECLSPEC int  mpfr_subnormalize _MPFR_PROTO ((mpfr_ptr, int,
 598 diff -Naurd mpfr-2.2.0-p4/round_near_x.c mpfr-2.2.0-p5/round_near_x.c
 599 --- mpfr-2.2.0-p4/round_near_x.c        2005-08-18 16:35:12.000000000 +0000
 600 +++ mpfr-2.2.0-p5/round_near_x.c        2005-11-23 09:04:29.000000000 +0000
 601 @@ -21,13 +21,13 @@
 602
 603  #include "mpfr-impl.h"
 604
 605 -/* Uses MPFR_FAST_COMPUTE_IF_SMALL_INPUT instead (a simple wrapper) */
 606 +/* Use MPFR_FAST_COMPUTE_IF_SMALL_INPUT instead (a simple wrapper) */
 607
 608  /* int mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr x, mp_exp_t err, int dir,
 609                            mp_rnd_t rnd)
 610
 611     Assuming y = o(f(x)) = o(x + g(x)) with |g(x)| < 2^(EXP(x)-error)
 612 -   If x is small enought, y ~= x. This function checks and does this.
 613 +   If x is small enough, y ~= x. This function checks and does this.
 614
 615     It assumes that f(x) is not representable exactly as a FP number.
 616     x must not be a singular value (NAN, INF or ZERO).
 617 @@ -42,7 +42,7 @@
 618     Otherwise it returns the ternary flag (It can't return an exact value).
 619  */
 620
 621 -/* What "small enought" means?
 622 +/* What "small enough" means?
 623
 624     We work with the positive values.
 625     Assuming err > Prec (y)+1
 626 @@ -50,7 +50,7 @@
 627     i = [ y = o(x)]   // i = inexact flag
 628     If i == 0
 629         Setting x in y is exact. We have:
 630 -       y = [XXXXXXXXX[...]]0[...] + error where [..] are optionnal zeros
 631 +       y = [XXXXXXXXX[...]]0[...] + error where [..] are optional zeros
 632        if dirError = ToInf,
 633          x < f(x) < x + 2^(EXP(x)-err)
 634          since x=y, and ulp (y)/2 > 2^(EXP(x)-err), we have:
 635 @@ -172,9 +172,17 @@
 636    sign = MPFR_SIGN (x);
 637    MPFR_SET_EXP (y, MPFR_GET_EXP (x));
 638    MPFR_SET_SIGN (y, sign);
 639 -  MPFR_RNDRAW_EVEN (inexact, y, MPFR_MANT (x), MPFR_PREC (x), rnd, sign,
 640 -                    if (MPFR_UNLIKELY ( ++MPFR_EXP (y) > __gmpfr_emax))
 641 -                    mpfr_overflow (y, rnd, sign) );
 642 +  MPFR_RNDRAW_GEN (inexact, y, MPFR_MANT (x), MPFR_PREC (x), rnd, sign,
 643 +                   if (dir == 0)
 644 +                     {
 645 +                       inexact = -sign;
 646 +                       goto trunc_doit;
 647 +                     }
 648 +                   else
 649 +                     goto addoneulp;
 650 +                   , if (MPFR_UNLIKELY (++MPFR_EXP (y) > __gmpfr_emax))
 651 +                       mpfr_overflow (y, rnd, sign)
 652 +                  );
 653
 654    /* Fix it in some cases */
 655    MPFR_ASSERTD (!MPFR_IS_NAN (y) && !MPFR_IS_ZERO (y));
 656 @@ -212,15 +220,6 @@
 657              }
 658          }
 659      }
 660 -  /* The even rule has been used. But due to error term, we should never
 661 -     use this rule. That's why we have to fix some wrong rounding */
 662 -  else if (inexact == MPFR_EVEN_INEX || inexact == -MPFR_EVEN_INEX)
 663 -    {
 664 -      if (inexact*sign > 0 && dir == 0)
 665 -        goto nexttozero;
 666 -      else if (inexact*sign < 0 && dir == 1)
 667 -        goto nexttoinf;
 668 -    }
 669
 670    MPFR_RET (inexact);
 671  }
 672 diff -Naurd mpfr-2.2.0-p4/tests/thypot.c mpfr-2.2.0-p5/tests/thypot.c
 673 --- mpfr-2.2.0-p4/tests/thypot.c        2005-06-02 16:12:04.000000000 +0000
 674 +++ mpfr-2.2.0-p5/tests/thypot.c        2005-11-23 09:04:29.000000000 +0000
 675 @@ -131,7 +131,26 @@
 676    inexact = mpfr_hypot (z, x, y, GMP_RNDN);
 677    if (inexact >= 0 || mpfr_cmp (x, z))
 678      {
 679 -      printf ("Error in test_large_small\n");
 680 +      printf ("Error 1 in test_large_small\n");
 681 +      exit (1);
 682 +    }
 683 +
 684 +  mpfr_mul_ui (x, x, 5, GMP_RNDN);
 685 +  inexact = mpfr_hypot (z, x, y, GMP_RNDN);
 686 +  if (mpfr_cmp (x, z) >= 0)
 687 +    {
 688 +      printf ("Error 2 in test_large_small\n");
 689 +      printf ("x = ");
 690 +      mpfr_out_str (stdout, 2, 0, x, GMP_RNDN);
 691 +      printf ("\n");
 692 +      printf ("y = ");
 693 +      mpfr_out_str (stdout, 2, 0, y, GMP_RNDN);
 694 +      printf ("\n");
 695 +      printf ("z = ");
 696 +      mpfr_out_str (stdout, 2, 0, z, GMP_RNDN);
 697 +      printf (" (in precision 2) instead of\n    ");
 698 +      mpfr_out_str (stdout, 2, 2, x, GMP_RNDU);
 699 +      printf ("\n");
 700        exit (1);
 701      }
 702
 703 diff -Naurd mpfr-2.2.0-p5/div.c mpfr-2.2.0-p6/div.c
 704 --- mpfr-2.2.0-p5/div.c 2005-08-18 17:03:05.000000000 +0000
 705 +++ mpfr-2.2.0-p6/div.c 2005-11-24 21:39:31.000000000 +0000
 706 @@ -298,17 +298,16 @@
 707            MPN_COPY(bp, vp, vsize);
 708          }
 709        sticky_v = sticky_v || mpn_cmpzero (vp, k);
 710 +      k = 0;
 711      }
 712 -  else /* vsize < qsize */
 713 +  else /* vsize < qsize: small divisor case */
 714      {
 715 +      bp = vp;
 716        k = qsize - vsize;
 717 -      bp = (mp_ptr) MPFR_TMP_ALLOC (qsize * sizeof(mp_limb_t));
 718 -      MPN_COPY(bp + k, vp, vsize);
 719 -      MPN_ZERO(bp, k);
 720      }
 721
 722    /* we now can perform the division */
 723 -  qh = mpn_divrem (qp, 0, ap, qqsize, bp, qsize);
 724 +  qh = mpn_divrem (qp, 0, ap + k, qqsize - k, bp, qsize - k);
 725    /* warning: qh may be 1 if u1 == v1, but u < v */
 726  #ifdef DEBUG
 727    printf ("q="); mpn_print (qp, qsize);
 728 diff -Naurd mpfr-2.2.0-p6/sin.c mpfr-2.2.0-p7/sin.c
 729 --- mpfr-2.2.0-p6/sin.c 2005-08-18 17:03:11.000000000 +0000
 730 +++ mpfr-2.2.0-p7/sin.c 2005-12-24 15:17:54.000000000 +0000
 731 @@ -162,10 +162,12 @@
 732          {
 733            /* the absolute error on c is at most 2^(3-m-EXP(c)) */
 734            e = 2 * MPFR_GET_EXP (c) + m - 3;
 735 -          if (mpfr_can_round (c, e, GMP_RNDZ, GMP_RNDZ,
 736 +          if (mpfr_can_round (c, e, GMP_RNDN, GMP_RNDZ,
 737                                precy + (rnd_mode == GMP_RNDN)))
 738 -            /* WARNING: need one more bit for rounding to nearest,
 739 -               to be able to get the inexact flag correct */
 740 +            /* WARNING: even if we know c <= sin(x), don't give GMP_RNDZ
 741 +               as 3rd argument to mpfr_can_round, since if c is exactly
 742 +               representable to the target precision (inexact = 0 below),
 743 +               we would have to add one ulp when rounding away from 0. */
 744              break;
 745
 746            /* check for huge cancellation (Near 0) */
 747 @@ -183,14 +185,8 @@
 748    MPFR_ZIV_FREE (loop);
 749
 750    inexact = mpfr_set (y, c, rnd_mode);
 751 -
 752 -  /* sin(x) is exact only for x = 0, which was treated apart above;
 753 -     nevertheless, we can have inexact = 0 here if the approximation c
 754 -     is exactly representable with PREC(y) bits. Since c is an approximation
 755 -     towards zero, in that case the inexact flag should have the opposite sign
 756 -     as y. */
 757 -  if (MPFR_UNLIKELY (inexact == 0))
 758 -    inexact = -MPFR_INT_SIGN (y);
 759 +  /* inexact cannot be 0, since this would mean that c was representable
 760 +     within the target precision, but in that case mpfr_can_round will fail */
 761
 762    mpfr_clear (c);
 763
 764 diff -Naurd mpfr-2.2.0-p7/get_f.c mpfr-2.2.0-p8/get_f.c
 765 --- mpfr-2.2.0-p7/get_f.c       2005-06-08 09:53:48.000000000 +0000
 766 +++ mpfr-2.2.0-p8/get_f.c       2006-01-13 15:04:34.000000000 +0000
 767 @@ -1,6 +1,6 @@
 768  /* mpfr_get_f -- convert a MPFR number to a GNU MPF number
 769
 770 -Copyright 2005 Free Software Foundation, Inc.
 771 +Copyright 2005, 2006 Free Software Foundation, Inc.
 772
 773  This file is part of the MPFR Library.
 774
 775 @@ -27,8 +27,9 @@
 776  int
 777  mpfr_get_f (mpf_ptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
 778  {
 779 -  unsigned long sx, sy, precx, precy, sh;
 780 -  mp_exp_t ey;
 781 +  mp_size_t sx, sy;
 782 +  mp_prec_t precx, precy;
 783 +  int sh;
 784
 785    if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(y)))
 786      {
 787 @@ -44,20 +45,18 @@
 788    sx = PREC(x); /* number of limbs of the mantissa of x */
 789
 790    precy = MPFR_PREC(y);
 791 -  precx = sx * BITS_PER_MP_LIMB;
 792 -  sy = 1 + (MPFR_PREC(y) - 1) / BITS_PER_MP_LIMB;
 793 +  precx = (mp_prec_t) sx * BITS_PER_MP_LIMB;
 794 +  sy = MPFR_LIMB_SIZE (y);
 795
 796    /* since mpf numbers are represented in base 2^BITS_PER_MP_LIMB,
 797       we loose -EXP(y) % BITS_PER_MP_LIMB bits in the most significant limb */
 798 -  ey = MPFR_GET_EXP(y) % BITS_PER_MP_LIMB;
 799 -  if (ey <= 0)
 800 -    sh = (unsigned long) (-ey);
 801 -  else /* 0 < ey < BITS_PER_MP_LIMB */
 802 -    sh = BITS_PER_MP_LIMB - (unsigned long) ey;
 803 +  sh = MPFR_GET_EXP(y) % BITS_PER_MP_LIMB;
 804 +  sh = sh <= 0 ? - sh : BITS_PER_MP_LIMB - sh;
 805 +  MPFR_ASSERTD (sh >= 0);
 806    if (precy + sh <= precx) /* we can copy directly */
 807      {
 808 -      /* necessarily sy <= sx */
 809 -      if (sh)
 810 +      MPFR_ASSERTN (sx >= sy);
 811 +      if (sh != 0)
 812          mpn_rshift (PTR(x) + sx - sy, MPFR_MANT(y), sy, sh);
 813        else
 814          MPN_COPY (PTR(x) + sx - sy, MPFR_MANT(y), sy);
 815 @@ -68,17 +67,17 @@
 816    else /* we have to round to precx - sh bits */
 817      {
 818        mpfr_t z;
 819 -      unsigned long sz;
 820 +      mp_size_t sz;
 821
 822        mpfr_init2 (z, precx - sh);
 823 -      sz = 1 + (MPFR_PREC(z) - 1) / BITS_PER_MP_LIMB;
 824 +      sz = MPFR_LIMB_SIZE (z);
 825        mpfr_set (z, y, rnd_mode);
 826        /* warning, sh may change due to rounding, but then z is a power of two,
 827           thus we can safely ignore its last bit which is 0 */
 828 -      ey = MPFR_GET_EXP(z) % BITS_PER_MP_LIMB;
 829 -      sh = (ey <= 0) ? (unsigned long) (-ey)
 830 -        : BITS_PER_MP_LIMB - (unsigned long) ey;
 831 -      if (sh)
 832 +      sh = MPFR_GET_EXP(z) % BITS_PER_MP_LIMB;
 833 +      sh = sh <= 0 ? - sh : BITS_PER_MP_LIMB - sh;
 834 +      MPFR_ASSERTD (sh >= 0);
 835 +      if (sh != 0)
 836          mpn_rshift (PTR(x) + sx - sz, MPFR_MANT(z), sz, sh);
 837        else
 838          MPN_COPY (PTR(x) + sx - sz, MPFR_MANT(z), sz);
 839 diff -Naurd mpfr-2.2.0-p7/set_f.c mpfr-2.2.0-p8/set_f.c
 840 --- mpfr-2.2.0-p7/set_f.c       2005-08-18 16:35:13.000000000 +0000
 841 +++ mpfr-2.2.0-p8/set_f.c       2006-01-12 10:34:40.000000000 +0000
 842 @@ -1,6 +1,6 @@
 843  /* mpfr_set_f -- set a MPFR number from a GNU MPF number
 844
 845 -Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005 Free Software Foundation, Inc.
 846 +Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 Free Software Foundation, Inc.
 847
 848  This file is part of the MPFR Library.
 849
 850 @@ -80,7 +80,22 @@
 851        inexact = 0;
 852      }
 853
 854 -  MPFR_SET_EXP(y, EXP(x) * BITS_PER_MP_LIMB - cnt + carry);
 855 +  /* warning: EXP(x) * BITS_PER_MP_LIMB may exceed the maximal exponent */
 856 +  if (EXP(x) > 1 + (__gmpfr_emax - 1) / BITS_PER_MP_LIMB)
 857 +    {
 858 +      /* EXP(x) >= 2 + floor((__gmpfr_emax-1)/BITS_PER_MP_LIMB)
 859 +         EXP(x) >= 2 + (__gmpfr_emax - BITS_PER_MP_LIMB) / BITS_PER_MP_LIMB
 860 +                >= 1 + __gmpfr_emax / BITS_PER_MP_LIMB
 861 +         EXP(x) * BITS_PER_MP_LIMB >= __gmpfr_emax + BITS_PER_MP_LIMB
 862 +         Since 0 <= cnt <= BITS_PER_MP_LIMB-1, and 0 <= carry <= 1,
 863 +         we have then EXP(x) * BITS_PER_MP_LIMB - cnt + carry > __gmpfr_emax */
 864 +      return mpfr_overflow (y, rnd_mode, MPFR_SIGN (y));
 865 +    }
 866 +  else
 867 +    {
 868 +      /* Do not use MPFR_SET_EXP as the exponent may be out of range. */
 869 +      MPFR_EXP (y) = EXP (x) * BITS_PER_MP_LIMB - (mp_exp_t) cnt + carry;
 870 +    }
 871
 872 -  return inexact;
 873 +  return mpfr_check_range (y, inexact, rnd_mode);
 874  }
 875 diff -Naurd mpfr-2.2.0-p7/tests/tget_f.c mpfr-2.2.0-p8/tests/tget_f.c
 876 --- mpfr-2.2.0-p7/tests/tget_f.c        2005-06-02 16:12:04.000000000 +0000
 877 +++ mpfr-2.2.0-p8/tests/tget_f.c        2006-01-13 15:05:14.000000000 +0000
 878 @@ -1,6 +1,6 @@
 879  /* Test file for mpfr_get_f.
 880
 881 -Copyright 2005 Free Software Foundation, Inc.
 882 +Copyright 2005, 2006 Free Software Foundation, Inc.
 883
 884  This file is part of the MPFR Library.
 885
 886 @@ -22,6 +22,7 @@
 887  #include <stdio.h>
 888  #include <stdlib.h>
 889  #include <time.h>
 890 +#include <limits.h>
 891
 892  #include "mpfr-test.h"
 893
 894 @@ -31,6 +32,7 @@
 895    mpf_t x;
 896    mpfr_t y;
 897    unsigned long i;
 898 +  mp_exp_t e;
 899
 900    MPFR_TEST_USE_RANDS ();
 901    tests_start_mpfr ();
 902 @@ -83,11 +85,14 @@
 903            printf ("Error: mpfr_get_f(%lu) fails\n", i);
 904            exit (1);
 905          }
 906 -      mpfr_set_si (y, (signed long) -i, GMP_RNDN);
 907 -      if (mpfr_get_f (x, y, GMP_RNDN) || mpf_cmp_si (x, (signed long) -i))
 908 +      if (i <= - (unsigned long) LONG_MIN)
 909          {
 910 -          printf ("Error: mpfr_get_f(-%lu) fails\n", i);
 911 -          exit (1);
 912 +          mpfr_set_si (y, - (long) i, GMP_RNDN);
 913 +          if (mpfr_get_f (x, y, GMP_RNDN) || mpf_cmp_si (x, - (long) i))
 914 +            {
 915 +              printf ("Error: mpfr_get_f(-%lu) fails\n", i);
 916 +              exit (1);
 917 +            }
 918          }
 919        i *= 2;
 920      }
 921 @@ -113,6 +118,42 @@
 922        i *= 2;
 923      }
 924
 925 +  /* bug reported by Jim White */
 926 +  for (e = 0; e <= 2 * BITS_PER_MP_LIMB; e++)
 927 +    {
 928 +      /* test with 2^(-e) */
 929 +      mpfr_set_ui (y, 1, GMP_RNDN);
 930 +      mpfr_div_2exp (y, y, e, GMP_RNDN);
 931 +      mpfr_get_f (x, y, GMP_RNDN);
 932 +      mpf_mul_2exp (x, x, e);
 933 +      if (mpf_cmp_ui (x, 1) != 0)
 934 +        {
 935 +          printf ("Error: mpfr_get_f(x,y,GMP_RNDN) fails\n");
 936 +          printf ("y=");
 937 +          mpfr_dump (y);
 938 +          printf ("x=");
 939 +          mpf_div_2exp (x, x, e);
 940 +          mpf_dump (x);
 941 +          exit (1);
 942 +        }
 943 +
 944 +      /* test with 2^(e) */
 945 +      mpfr_set_ui (y, 1, GMP_RNDN);
 946 +      mpfr_mul_2exp (y, y, e, GMP_RNDN);
 947 +      mpfr_get_f (x, y, GMP_RNDN);
 948 +      mpf_div_2exp (x, x, e);
 949 +      if (mpf_cmp_ui (x, 1) != 0)
 950 +        {
 951 +          printf ("Error: mpfr_get_f(x,y,GMP_RNDN) fails\n");
 952 +          printf ("y=");
 953 +          mpfr_dump (y);
 954 +          printf ("x=");
 955 +          mpf_mul_2exp (x, x, e);
 956 +          mpf_dump (x);
 957 +          exit (1);
 958 +        }
 959 +    }
 960 +
 961    mpfr_clear (y);
 962    mpf_clear (x);
 963
 964 diff -Naurd mpfr-2.2.0-p7/tests/tset_f.c mpfr-2.2.0-p8/tests/tset_f.c
 965 --- mpfr-2.2.0-p7/tests/tset_f.c        2005-09-09 15:18:00.000000000 +0000
 966 +++ mpfr-2.2.0-p8/tests/tset_f.c        2006-01-12 10:31:42.000000000 +0000
 967 @@ -1,6 +1,6 @@
 968  /* Test file for mpfr_set_f.
 969
 970 -Copyright 1999, 2001, 2002, 2003, 2004, 2005 Free Software Foundation, Inc.
 971 +Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006 Free Software Foundation, Inc.
 972
 973  This file is part of the MPFR Library.
 974
 975 @@ -22,6 +22,7 @@
 976  #include <stdio.h>
 977  #include <stdlib.h>
 978  #include <time.h>
 979 +#include <limits.h> /* for ULONG_MAX */
 980
 981  #include "mpfr-test.h"
 982
 983 @@ -30,6 +31,7 @@
 984  {
 985    mpfr_t x, u;
 986    mpf_t y, z;
 987 +  mp_exp_t emax;
 988    unsigned long k, pr;
 989    int r, inexact;
 990
 991 @@ -87,8 +89,6 @@
 992      }
 993    MPFR_ASSERTN(mpfr_cmp_ui_2exp (x, 1, 901) == 0);
 994
 995 -  mpfr_clear (u);
 996 -
 997    for (k = 1; k <= 100000; k++)
 998      {
 999        pr = 2 + (randlimb () & 255);
1000 @@ -131,7 +131,49 @@
1001        mpf_mul_2exp (y, y, 1);
1002      }
1003
1004 +  mpf_set_ui (y, 1);
1005 +  mpf_mul_2exp (y, y, ULONG_MAX);
1006 +  mpfr_set_f (x, y, GMP_RNDN);
1007 +  mpfr_set_ui (u, 1, GMP_RNDN);
1008 +  mpfr_mul_2ui (u, u, ULONG_MAX, GMP_RNDN);
1009 +  if (!mpfr_equal_p (x, u))
1010 +    {
1011 +      printf ("Error: mpfr_set_f (x, y, GMP_RNDN) for y = 2^ULONG_MAX\n");
1012 +      exit (1);
1013 +    }
1014 +
1015 +  emax = mpfr_get_emax ();
1016 +
1017 +  /* For mpf_mul_2exp, emax must fit in an unsigned long! */
1018 +  if (emax >= 0 && emax <= ULONG_MAX)
1019 +    {
1020 +      mpf_set_ui (y, 1);
1021 +      mpf_mul_2exp (y, y, emax);
1022 +      mpfr_set_f (x, y, GMP_RNDN);
1023 +      mpfr_set_ui_2exp (u, 1, emax, GMP_RNDN);
1024 +      if (!mpfr_equal_p (x, u))
1025 +        {
1026 +          printf ("Error: mpfr_set_f (x, y, GMP_RNDN) for y = 2^emax\n");
1027 +          exit (1);
1028 +        }
1029 +    }
1030 +
1031 +  /* For mpf_mul_2exp, emax - 1 must fit in an unsigned long! */
1032 +  if (emax >= 1 && emax - 1 <= ULONG_MAX)
1033 +    {
1034 +      mpf_set_ui (y, 1);
1035 +      mpf_mul_2exp (y, y, emax - 1);
1036 +      mpfr_set_f (x, y, GMP_RNDN);
1037 +      mpfr_set_ui_2exp (u, 1, emax - 1, GMP_RNDN);
1038 +      if (!mpfr_equal_p (x, u))
1039 +        {
1040 +          printf ("Error: mpfr_set_f (x, y, GMP_RNDN) for y = 2^(emax-1)\n");
1041 +          exit (1);
1042 +        }
1043 +    }
1044 +
1045    mpfr_clear (x);
1046 +  mpfr_clear (u);
1047    mpf_clear (y);
1048    mpf_clear (z);
1049
1050 diff -Naurd mpfr-2.2.0-p8/random2.c mpfr-2.2.0-p9/random2.c
1051 --- mpfr-2.2.0-p8/random2.c     2005-08-18 17:03:10.000000000 +0000
1052 +++ mpfr-2.2.0-p9/random2.c     2006-02-20 10:29:01.000000000 +0000
1053 @@ -2,7 +2,7 @@
1054     long runs of consecutive ones and zeros in the binary representation.
1055     Intended for testing.
1056
1057 -Copyright 1999, 2001, 2002, 2003, 2004 Free Software Foundation, Inc.
1058 +Copyright 1999, 2001, 2002, 2003, 2004, 2006 Free Software Foundation, Inc.
1059
1060  This file is part of the MPFR Library.
1061
1062 @@ -36,7 +36,8 @@
1063
1064    if (MPFR_UNLIKELY(size == 0))
1065      {
1066 -      MPFR_SET_ZERO(x);
1067 +      MPFR_SET_ZERO (x);
1068 +      MPFR_SET_POS (x);
1069        return ;
1070      }
1071    else if (size > 0)
1072 diff -Naurd mpfr-2.2.0-p8/tests/trandom.c mpfr-2.2.0-p9/tests/trandom.c
1073 --- mpfr-2.2.0-p8/tests/trandom.c       2005-08-18 17:03:16.000000000 +0000
1074 +++ mpfr-2.2.0-p9/tests/trandom.c       2006-02-20 10:29:01.000000000 +0000
1075 @@ -1,6 +1,6 @@
1076  /* Test file for the various mpfr_random fonctions.
1077
1078 -Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation.
1079 +Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006 Free Software Foundation, Inc.
1080
1081  This file is part of the MPFR Library.
1082
1083 @@ -120,7 +120,10 @@
1084
1085    /* test size=0 */
1086    mpfr_random2 (x, 0, 0);
1087 -  MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS(x));
1088 +  MPFR_ASSERTN (mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS (x));
1089 +  mpfr_set_si (x, -1, GMP_RNDN); /* x is negative */
1090 +  mpfr_random2 (x, 0, 0);
1091 +  MPFR_ASSERTN (mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS (x));
1092
1093    /* test size < 0 */
1094    mpfr_random2 (x, -1, 0);
1095 diff -Naurd mpfr-2.2.0-p9/asin.c mpfr-2.2.0-p10/asin.c
1096 --- mpfr-2.2.0-p9/asin.c        2005-08-18 17:03:04.000000000 +0000
1097 +++ mpfr-2.2.0-p10/asin.c       2006-03-16 17:47:51.000000000 +0000
1098 @@ -51,7 +51,7 @@
1099      }
1100
1101    /* asin(x) = x + x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
1102 -  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (asin,x, -2*MPFR_GET_EXP (x)+2,1,rnd_mode,);
1103 +  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (asin,x, -2*MPFR_GET_EXP (x)+2,1,rnd_mode,{});
1104
1105    /* Set x_p=|x| (x is a normal number) */
1106    mpfr_init2 (xp, MPFR_PREC (x));
1107 diff -Naurd mpfr-2.2.0-p9/asinh.c mpfr-2.2.0-p10/asinh.c
1108 --- mpfr-2.2.0-p9/asinh.c       2005-08-18 17:03:04.000000000 +0000
1109 +++ mpfr-2.2.0-p10/asinh.c      2006-03-16 17:47:51.000000000 +0000
1110 @@ -62,7 +62,7 @@
1111      }
1112
1113    /* asinh(x) = x - x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
1114 -  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2*MPFR_GET_EXP (x)+2,0,rnd_mode,);
1115 +  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2*MPFR_GET_EXP (x)+2,0,rnd_mode,{});
1116
1117    Ny = MPFR_PREC (y);   /* Precision of output variable */
1118
1119 diff -Naurd mpfr-2.2.0-p9/atan.c mpfr-2.2.0-p10/atan.c
1120 --- mpfr-2.2.0-p9/atan.c        2005-08-18 16:35:03.000000000 +0000
1121 +++ mpfr-2.2.0-p10/atan.c       2006-03-16 17:47:51.000000000 +0000
1122 @@ -185,7 +185,7 @@
1123    /* atan(x) = x - x^3/3 + x^5/5...
1124       so the error is < 2^(3*EXP(x)-1)
1125       so `EXP(x)-(3*EXP(x)-1)` = -2*EXP(x)+1 */
1126 -  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (atan,x, -2*MPFR_GET_EXP (x)+1,0,rnd_mode,);
1127 +  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (atan,x, -2*MPFR_GET_EXP (x)+1,0,rnd_mode,{});
1128
1129    /* Set x_p=|x| */
1130    MPFR_TMP_INIT_ABS (xp, x);
1131 diff -Naurd mpfr-2.2.0-p9/atanh.c mpfr-2.2.0-p10/atanh.c
1132 --- mpfr-2.2.0-p9/atanh.c       2005-08-18 17:03:04.000000000 +0000
1133 +++ mpfr-2.2.0-p10/atanh.c      2006-03-16 17:47:51.000000000 +0000
1134 @@ -71,7 +71,7 @@
1135      }
1136
1137    /* atanh(x) = x + x^3/3 + ... so the error is < 2^(3*EXP(x)-1) */
1138 -  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2*MPFR_GET_EXP (xt)+1,1,rnd_mode,);
1139 +  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2*MPFR_GET_EXP (xt)+1,1,rnd_mode,{});
1140
1141    MPFR_SAVE_EXPO_MARK (expo);
1142
1143 diff -Naurd mpfr-2.2.0-p9/expm1.c mpfr-2.2.0-p10/expm1.c
1144 --- mpfr-2.2.0-p9/expm1.c       2005-08-18 17:03:06.000000000 +0000
1145 +++ mpfr-2.2.0-p10/expm1.c      2006-03-16 17:47:51.000000000 +0000
1146 @@ -61,7 +61,7 @@
1147      }
1148
1149    /* exp(x)-1 = x +x^2/2 + ... so the error is < 2^(2*EXP(x)-1) */
1150 -  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -MPFR_GET_EXP (x)+1,1,rnd_mode,);
1151 +  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -MPFR_GET_EXP (x)+1,1,rnd_mode,{});
1152
1153    MPFR_SAVE_EXPO_MARK (expo);
1154    /* General case */
1155 diff -Naurd mpfr-2.2.0-p9/log1p.c mpfr-2.2.0-p10/log1p.c
1156 --- mpfr-2.2.0-p9/log1p.c       2005-08-18 17:03:08.000000000 +0000
1157 +++ mpfr-2.2.0-p10/log1p.c      2006-03-16 17:47:51.000000000 +0000
1158 @@ -63,7 +63,7 @@
1159      }
1160
1161    /* log(1+x) = x-x^2/2 + ... so the error is < 2^(2*EXP(x)-1) */
1162 -  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -MPFR_GET_EXP (x)+1,0,rnd_mode,);
1163 +  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -MPFR_GET_EXP (x)+1,0,rnd_mode,{});
1164
1165    comp = mpfr_cmp_si (x, -1);
1166    /* log1p(x) is undefined for x < -1 */
1167 diff -Naurd mpfr-2.2.0-p9/sin.c mpfr-2.2.0-p10/sin.c
1168 --- mpfr-2.2.0-p9/sin.c 2005-12-24 15:17:54.000000000 +0000
1169 +++ mpfr-2.2.0-p10/sin.c        2006-03-16 17:47:51.000000000 +0000
1170 @@ -129,7 +129,7 @@
1171      }
1172
1173    /* sin(x) = x - x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
1174 -  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2*MPFR_GET_EXP (x)+2,0,rnd_mode, );
1175 +  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2*MPFR_GET_EXP (x)+2,0,rnd_mode,{});
1176
1177    /* Compute initial precision */
1178    precy = MPFR_PREC (y);
1179 diff -Naurd mpfr-2.2.0-p9/sinh.c mpfr-2.2.0-p10/sinh.c
1180 --- mpfr-2.2.0-p9/sinh.c        2005-08-18 16:35:14.000000000 +0000
1181 +++ mpfr-2.2.0-p10/sinh.c       2006-03-16 17:47:51.000000000 +0000
1182 @@ -57,7 +57,7 @@
1183      }
1184
1185    /* sinh(x) = x + x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
1186 -  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2*MPFR_GET_EXP(xt)+2,1,rnd_mode, );
1187 +  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2*MPFR_GET_EXP(xt)+2,1,rnd_mode,{});
1188
1189    MPFR_TMP_INIT_ABS (x, xt);
1190
1191 diff -Naurd mpfr-2.2.0-p9/tan.c mpfr-2.2.0-p10/tan.c
1192 --- mpfr-2.2.0-p9/tan.c 2005-08-18 16:35:15.000000000 +0000
1193 +++ mpfr-2.2.0-p10/tan.c        2006-03-16 17:47:51.000000000 +0000
1194 @@ -53,7 +53,7 @@
1195      }
1196
1197    /* tan(x) = x + x^3/3 + ... so the error is < 2^(3*EXP(x)-1) */
1198 -  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2*MPFR_GET_EXP (x)+1,1,rnd_mode, );
1199 +  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2*MPFR_GET_EXP (x)+1,1,rnd_mode,{});
1200
1201    MPFR_SAVE_EXPO_MARK (expo);
1202
1203 diff -Naurd mpfr-2.2.0-p9/tanh.c mpfr-2.2.0-p10/tanh.c
1204 --- mpfr-2.2.0-p9/tanh.c        2005-08-18 17:03:12.000000000 +0000
1205 +++ mpfr-2.2.0-p10/tanh.c       2006-03-16 17:47:51.000000000 +0000
1206 @@ -56,7 +56,7 @@
1207      }
1208
1209    /* tanh(x) = x - x^3/3 + ... so the error is < 2^(3*EXP(x)-1) */
1210 -  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2*MPFR_GET_EXP(xt)+1,0,rnd_mode, );
1211 +  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2*MPFR_GET_EXP(xt)+1,0,rnd_mode,{});
1212
1213    MPFR_TMP_INIT_ABS (x, xt);
1214
1215 diff -Naurd mpfr-2.2.0-p9/zeta.c mpfr-2.2.0-p10/zeta.c
1216 --- mpfr-2.2.0-p9/zeta.c        2005-09-16 14:36:54.000000000 +0000
1217 +++ mpfr-2.2.0-p10/zeta.c       2006-03-16 17:47:51.000000000 +0000
1218 @@ -175,7 +175,7 @@
1219        else
1220          err = ((mp_exp_t)1) << err;
1221        err = 1 - (-err+1); /* GET_EXP(one) - (-err+1) = err :) */
1222 -      MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, err, 1, rnd_mode, );
1223 +      MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, err, 1, rnd_mode,{});
1224      }
1225
1226    d = precz + MPFR_INT_CEIL_LOG2(precz) + 10;