static void
special (void)
{
  mpfr_t x;
  int inex;

  mpfr_init (x);

  mpfr_set_nan (x);
  inex = test_log1p (x, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (x) && inex == 0);

  mpfr_set_inf (x, -1);
  inex = test_log1p (x, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (x) && inex == 0);

  mpfr_set_inf (x, 1);
  inex = test_log1p (x, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_inf_p (x) && mpfr_sgn (x) > 0 && inex == 0);

  mpfr_set_ui (x, 0, MPFR_RNDN);
  inex = test_log1p (x, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS (x) && inex == 0);
  mpfr_neg (x, x, MPFR_RNDN);
  inex = test_log1p (x, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_NEG (x) && inex == 0);

  mpfr_set_si (x, -1, MPFR_RNDN);
  inex = test_log1p (x, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_inf_p (x) && mpfr_sgn (x) < 0 && inex == 0);

  mpfr_set_si (x, -2, MPFR_RNDN);
  inex = test_log1p (x, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (x) && inex == 0);

  mpfr_clear (x);
}

static void
special (void)
{
  mpfr_t x, y, z;

  mpfr_init (x);
  mpfr_init (y);
  mpfr_init (z);

  mpfr_set_nan (x);
  mpfr_set_ui (y, 0, MPFR_RNDN);
  mpfr_hypot (z, x, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_nan_p (z));
  mpfr_hypot (z, x, y, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_nan_p (z));
  mpfr_hypot (z, y, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_nan_p (z));

  mpfr_set_inf (x, 1);
  mpfr_set_inf (y, -1);
  mpfr_hypot (z, x, y, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_inf_p (z) && mpfr_sgn (z) > 0);

  mpfr_set_inf (x, -1);
  mpfr_set_nan (y);
  mpfr_hypot (z, x, y, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_inf_p (z) && mpfr_sgn (z) > 0);

  mpfr_set_nan (x);
  mpfr_set_inf (y, -1);
  mpfr_hypot (z, x, y, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_inf_p (z) && mpfr_sgn (z) > 0);

  mpfr_clear (x);
  mpfr_clear (y);
  mpfr_clear (z);
}

static void
check_special (void)
{
  mpfr_t x, y;
  mpfr_exp_t emin;

  mpfr_init (x);
  mpfr_init (y);

  mpfr_set_nan (x);
  mpfr_sqr (y, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (y));

  mpfr_set_inf (x, 1);
  mpfr_sqr (y, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_inf_p (y) && mpfr_sgn (y) > 0);

  mpfr_set_inf (x, -1);
  mpfr_sqr (y, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_inf_p (y) && mpfr_sgn (y) > 0);

  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_sqr (y, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_zero_p (y));

  emin = mpfr_get_emin ();
  mpfr_set_emin (0);
  mpfr_set_ui (x, 1, MPFR_RNDN);
  mpfr_div_2ui (x, x, 1, MPFR_RNDN);
  MPFR_ASSERTN (!mpfr_zero_p (x));
  mpfr_sqr (y, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_zero_p (y));
  mpfr_set_emin (emin);

  mpfr_clear (y);
  mpfr_clear (x);
}

//.........这里部分代码省略.........
                  bbcp1 = (kx>=0);
                  DEBUG (printf ("(Pre) Scan done for %ld\n", (long) kx));
                }
            } /*End of Bcp1 != 0*/
          DEBUG( printf("(Pre) Cp+1=%d C'p+2=%d\n", bbcp!=0, bbcp1!=0) );
        } /* End of "can lose a bit" */

      /* Clean shifted C' */
      mask = ~MPFR_LIMB_MASK (sh);
      cp[0] &= mask;

      /* Subtract the mantissa c from b in a */
      ap = MPFR_MANT(a);
      mpn_sub_n (ap, bp, cp, n);
      DEBUG( mpfr_print_mant_binary("Sub=  ", ap, p) );

     /* Normalize: we lose at max one bit*/
      if (MPFR_UNLIKELY(MPFR_LIMB_MSB(ap[n-1]) == 0))
        {
          /* High bit is not set and we have to fix it! */
          /* Ap >= 010000xxx001 */
          mpn_lshift(ap, ap, n, 1);
          /* Ap >= 100000xxx010 */
          if (MPFR_UNLIKELY(bcp!=0)) /* Check if Cp = -1 */
            /* Since Cp == -1, we have to substract one more */
            {
              mpn_sub_1(ap, ap, n, MPFR_LIMB_ONE<<sh);
              MPFR_ASSERTD(MPFR_LIMB_MSB(ap[n-1]) != 0);
            }
          /* Ap >= 10000xxx001 */
          /* Final exponent -1 since we have shifted the mantissa */
          bx--;
          /* Update bcp and bcp1 */
          MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
          MPFR_ASSERTN(bbcp1 != (mp_limb_t) -1);
          bcp  = bbcp;
          bcp1 = bbcp1;
          /* We dont't have anymore a valid Cp+1!
             But since Ap >= 100000xxx001, the final sub can't unnormalize!*/
        }
      MPFR_ASSERTD( !(ap[0] & ~mask) );

      /* Rounding */
      if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
        {
          if (MPFR_LIKELY(bcp==0))
            goto truncate;
          else if ((bcp1) || ((ap[0] & (MPFR_LIMB_ONE<<sh)) != 0))
            goto sub_one_ulp;
          else
            goto truncate;
        }

      /* Update rounding mode */
      MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
      if (rnd_mode == MPFR_RNDZ && (MPFR_LIKELY(bcp || bcp1)))
        goto sub_one_ulp;
      goto truncate;
    }
  MPFR_RET_NEVER_GO_HERE ();

  /* Sub one ulp to the result */
 sub_one_ulp:
  mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
  /* Result should be smaller than exact value: inexact=-1 */
  inexact = -1;

/* (y, z) <- (sin(x), cos(x)), return value is 0 iff both results are exact
   ie, iff x = 0 */
int
mpfr_sin_cos (mpfr_ptr y, mpfr_ptr z, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_prec_t prec, m;
  int neg, reduce;
  mpfr_t c, xr;
  mpfr_srcptr xx;
  mpfr_exp_t err, expx;
  int inexy, inexz;
  MPFR_ZIV_DECL (loop);
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_ASSERTN (y != z);

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN(x) || MPFR_IS_INF(x))
        {
          MPFR_SET_NAN (y);
          MPFR_SET_NAN (z);
          MPFR_RET_NAN;
        }
      else /* x is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
          MPFR_SET_ZERO (y);
          MPFR_SET_SAME_SIGN (y, x);
          /* y = 0, thus exact, but z is inexact in case of underflow
             or overflow */
          inexy = 0; /* y is exact */
          inexz = mpfr_set_ui (z, 1, rnd_mode);
          return INEX(inexy,inexz);
        }
    }

  MPFR_LOG_FUNC
    (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
     ("sin[%Pu]=%.*Rg cos[%Pu]=%.*Rg", mpfr_get_prec(y), mpfr_log_prec, y,
      mpfr_get_prec (z), mpfr_log_prec, z));

  MPFR_SAVE_EXPO_MARK (expo);

  prec = MAX (MPFR_PREC (y), MPFR_PREC (z));
  m = prec + MPFR_INT_CEIL_LOG2 (prec) + 13;
  expx = MPFR_GET_EXP (x);

  /* When x is close to 0, say 2^(-k), then there is a cancellation of about
     2k bits in 1-cos(x)^2. FIXME: in that case, it would be more efficient
     to compute sin(x) directly. VL: This is partly done by using
     MPFR_FAST_COMPUTE_IF_SMALL_INPUT from the mpfr_sin and mpfr_cos
     functions. Moreover, any overflow on m is avoided. */
  if (expx < 0)
    {
      /* Warning: in case y = x, and the first call to
         MPFR_FAST_COMPUTE_IF_SMALL_INPUT succeeds but the second fails,
         we will have clobbered the original value of x.
         The workaround is to first compute z = cos(x) in that case, since
         y and z are different. */
      if (y != x)
        /* y and x differ, thus we can safely try to compute y first */
        {
          MPFR_FAST_COMPUTE_IF_SMALL_INPUT (
            y, x, -2 * expx, 2, 0, rnd_mode,
            { inexy = _inexact;
              goto small_input; });

static void
special (void)
{
  mpfr_t x;
  int i;

  mpfr_init (x);

  mpfr_set_nan (x);
  mpfr_sinh (x, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_nan_p (x));

  mpfr_set_inf (x, 1);
  mpfr_sinh (x, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);
  mpfr_set_inf (x, -1);
  mpfr_sinh (x, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) < 0);

  mpfr_set_prec (x, 10);
  mpfr_set_str_binary (x, "-0.1001011001");
  mpfr_sinh (x, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_cmp_si_2exp (x, -159, -8) == 0);

  /* corner case */
  mpfr_set_prec (x, 2);
  mpfr_set_str_binary (x, "1E-6");
  mpfr_sinh (x, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_cmp_ui_2exp (x, 1, -6) == 0);

  mpfr_clear_flags ();
  mpfr_set_str_binary (x, "1E1000000000");
  i = mpfr_sinh (x, x, MPFR_RNDN);
  MPFR_ASSERTN (MPFR_IS_INF (x) && MPFR_SIGN (x) > 0);
  MPFR_ASSERTN (mpfr_overflow_p ());
  MPFR_ASSERTN (i == 1);

  mpfr_clear_flags ();
  mpfr_set_str_binary (x, "-1E1000000000");
  i = mpfr_sinh (x, x, MPFR_RNDN);
  MPFR_ASSERTN (MPFR_IS_INF (x) && MPFR_SIGN (x) < 0);
  MPFR_ASSERTN (mpfr_overflow_p () && !mpfr_underflow_p ());
  MPFR_ASSERTN (i == -1);

  mpfr_clear_flags ();
  mpfr_set_str_binary (x, "-1E1000000000");
  i = mpfr_sinh (x, x, MPFR_RNDD);
  MPFR_ASSERTN (MPFR_IS_INF (x) && MPFR_SIGN (x) < 0);
  MPFR_ASSERTN (mpfr_overflow_p () && !mpfr_underflow_p ());
  MPFR_ASSERTN (i == -1);

  mpfr_clear_flags ();
  mpfr_set_str_binary (x, "-1E1000000000");
  i = mpfr_sinh (x, x, MPFR_RNDU);
  MPFR_ASSERTN (!MPFR_IS_INF (x) && MPFR_SIGN (x) < 0);
  MPFR_ASSERTN (mpfr_overflow_p () && !mpfr_underflow_p ());
  MPFR_ASSERTN (i == 1);

  mpfr_clear (x);
}

static void
check_parse (void)
{
  mpfr_t x;
  char *s;
  int res;

  mpfr_init (x);

  /* Invalid data */
  mpfr_set_si (x, -1, MPFR_RNDN);
  res = mpfr_strtofr (x, "  invalid", NULL, 10, MPFR_RNDN);
  if (MPFR_NOTZERO (x) || MPFR_IS_NEG (x))
    {
      printf ("Failed parsing '  invalid' (1)\n X=");
      mpfr_dump (x);
      exit (1);
    }
  MPFR_ASSERTN (res == 0);
  mpfr_set_si (x, -1, MPFR_RNDN);
  res = mpfr_strtofr (x, "  invalid", &s, 0, MPFR_RNDN);
  if (MPFR_NOTZERO (x) || MPFR_IS_NEG (x) || strcmp (s, "  invalid"))
    {
      printf ("Failed parsing '  invalid' (2)\n S=%s\n X=", s);
      mpfr_dump (x);
      exit (1);
    }
  MPFR_ASSERTN (res == 0);
  /* Check if it stops correctly */
  mpfr_strtofr (x, "15*x", &s, 10, MPFR_RNDN);
  if (mpfr_cmp_ui (x, 15) || strcmp (s, "*x"))
    {
      printf ("Failed parsing '15*x'\n S=%s\n X=", s);
      mpfr_dump (x);
      exit (1);
    }
  /* Check for leading spaces */
  mpfr_strtofr (x, "  1.5E-10 *x^2", &s, 10, MPFR_RNDN);
  if (mpfr_cmp_str1 (x, "1.5E-10") || strcmp (s, " *x^2"))
    {
      printf ("Failed parsing '1.5E-10*x^2'\n S=%s\n X=", s);
      mpfr_dump (x);
      exit (1);
    }
  /* Check for leading sign */
  mpfr_strtofr (x, "  +17.5E-42E ", &s, 10, MPFR_RNDN);
  if (mpfr_cmp_str1 (x, "17.5E-42") || strcmp (s, "E "))
    {
      printf ("Failed parsing '+17.5E-42E '\n S=%s\n X=", s);
      mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); putchar ('\n');
      exit (1);
    }
  mpfr_strtofr (x, "-17.5E+42E\n", &s, 10, MPFR_RNDN);
  if (mpfr_cmp_str1 (x, "-17.5E42") || strcmp (s, "E\n"))
    {
      printf ("Failed parsing '-17.5E+42\\n'\n S=%s\n X=", s);
      mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); putchar ('\n');
      exit (1);
    }
  /* P form */
  mpfr_strtofr (x, "0x42P17", &s, 16, MPFR_RNDN);
  if (mpfr_cmp_str (x, "8650752", 10, MPFR_RNDN) || *s != 0)
    {
      printf ("Failed parsing '0x42P17' (base = 16)\n S='%s'\n X=", s);
      mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); putchar ('\n');
      exit (1);
    }
  mpfr_strtofr (x, "-0X42p17", &s, 16, MPFR_RNDN);
  if (mpfr_cmp_str (x, "-8650752", 10, MPFR_RNDN) || *s != 0)
    {
      printf ("Failed parsing '-0x42p17' (base = 16)\n S='%s'\n X=", s);
      mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); putchar ('\n');
      exit (1);
    }
  mpfr_strtofr (x, "42p17", &s, 16, MPFR_RNDN);
  if (mpfr_cmp_str (x, "8650752", 10, MPFR_RNDN) || *s != 0)
    {
      printf ("Failed parsing '42p17' (base = 16)\n S='%s'\n X=", s);
      mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); putchar ('\n');
      exit (1);
    }
  mpfr_strtofr (x, "-42P17", &s, 16, MPFR_RNDN);
  if (mpfr_cmp_str (x, "-8650752", 10, MPFR_RNDN) || *s != 0)
    {
      printf ("Failed parsing '-42P17' (base = 16)\n S='%s'\n X=", s);
      mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); putchar ('\n');
      exit (1);
    }
  mpfr_strtofr (x, "0b1001P17", &s, 2, MPFR_RNDN);
  if (mpfr_cmp_str (x, "1179648", 10, MPFR_RNDN) || *s != 0)
    {
      printf ("Failed parsing '0b1001P17' (base = 2)\n S='%s'\n X=", s);
      mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); putchar ('\n');
      exit (1);
    }
  mpfr_strtofr (x, "-0B1001p17", &s, 2, MPFR_RNDN);
  if (mpfr_cmp_str (x, "-1179648", 10, MPFR_RNDN) || *s != 0)
    {
      printf ("Failed parsing '-0B1001p17' (base = 2)\n S='%s'\n X=", s);
      mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); putchar ('\n');
//.........这里部分代码省略.........

static void
check_nans (void)
{
  mpfr_t  x, y;
  int inexact;

  mpfr_init2 (x, 123);
  mpfr_init2 (y, 123);

  /* 1.0 - nan is nan */
  mpfr_set_nan (x);
  mpfr_clear_flags ();
  inexact = mpfr_d_sub (y, 1.0, x, MPFR_RNDN);
  MPFR_ASSERTN (inexact == 0);
  MPFR_ASSERTN ((__gmpfr_flags ^ MPFR_FLAGS_NAN) == 0);
  MPFR_ASSERTN (mpfr_nan_p (y));

  /* 1.0 - +inf == -inf */
  mpfr_set_inf (x, 1);
  mpfr_clear_flags ();
  inexact = mpfr_d_sub (y, 1.0, x, MPFR_RNDN);
  MPFR_ASSERTN (inexact == 0);
  MPFR_ASSERTN (__gmpfr_flags == 0);
  MPFR_ASSERTN (mpfr_inf_p (y));
  MPFR_ASSERTN (MPFR_IS_NEG (y));

  /* 1.0 - -inf == +inf */
  mpfr_set_inf (x, -1);
  mpfr_clear_flags ();
  inexact = mpfr_d_sub (y, 1.0, x, MPFR_RNDN);
  MPFR_ASSERTN (inexact == 0);
  MPFR_ASSERTN (__gmpfr_flags == 0);
  MPFR_ASSERTN (mpfr_inf_p (y));
  MPFR_ASSERTN (MPFR_IS_POS (y));

  mpfr_clear (x);
  mpfr_clear (y);
}

int
mpfr_mul_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mpfr_rnd_t rnd_mode)
{
  mp_limb_t *yp;
  mp_size_t xn;
  int cnt, inexact;
  MPFR_TMP_DECL (marker);

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN (x))
        {
          MPFR_SET_NAN (y);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF (x))
        {
          if (u != 0)
            {
              MPFR_SET_INF (y);
              MPFR_SET_SAME_SIGN (y, x);
              MPFR_RET (0); /* infinity is exact */
            }
          else /* 0 * infinity */
            {
              MPFR_SET_NAN (y);
              MPFR_RET_NAN;
            }
        }
      else /* x is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
          MPFR_SET_ZERO (y);
          MPFR_SET_SAME_SIGN (y, x);
          MPFR_RET (0); /* zero is exact */
        }
    }
  else if (MPFR_UNLIKELY (u <= 1))
    {
      if (u < 1)
        {
          MPFR_SET_ZERO (y);
          MPFR_SET_SAME_SIGN (y, x);
          MPFR_RET (0); /* zero is exact */
        }
      else
        return mpfr_set (y, x, rnd_mode);
    }
  else if (MPFR_UNLIKELY (IS_POW2 (u)))
    return mpfr_mul_2si (y, x, MPFR_INT_CEIL_LOG2 (u), rnd_mode);

  yp = MPFR_MANT (y);
  xn = MPFR_LIMB_SIZE (x);

  MPFR_ASSERTD (xn < MP_SIZE_T_MAX);
  MPFR_TMP_MARK(marker);
  yp = MPFR_TMP_LIMBS_ALLOC (xn + 1);

  MPFR_ASSERTN (u == (mp_limb_t) u);
  yp[xn] = mpn_mul_1 (yp, MPFR_MANT (x), xn, u);

  /* x * u is stored in yp[xn], ..., yp[0] */

  /* since the case u=1 was treated above, we have u >= 2, thus
     yp[xn] >= 1 since x was msb-normalized */
  MPFR_ASSERTD (yp[xn] != 0);
  if (MPFR_LIKELY (MPFR_LIMB_MSB (yp[xn]) == 0))
    {
      count_leading_zeros (cnt, yp[xn]);
      mpn_lshift (yp, yp, xn + 1, cnt);
    }
  else
    {
      cnt = 0;
    }

  /* now yp[xn], ..., yp[0] is msb-normalized too, and has at most
     PREC(x) + (GMP_NUMB_BITS - cnt) non-zero bits */
  MPFR_RNDRAW (inexact, y, yp, (mpfr_prec_t) (xn + 1) * GMP_NUMB_BITS,
               rnd_mode, MPFR_SIGN (x), cnt -- );

  MPFR_TMP_FREE (marker);

  cnt = GMP_NUMB_BITS - cnt;
  if (MPFR_UNLIKELY (__gmpfr_emax < MPFR_EMAX_MIN + cnt
                     || MPFR_GET_EXP (x) > __gmpfr_emax - cnt))
    return mpfr_overflow (y, rnd_mode, MPFR_SIGN(x));

  MPFR_SET_EXP (y, MPFR_GET_EXP (x) + cnt);
  MPFR_SET_SAME_SIGN (y, x);

  return inexact;
}

int
main (void)
{
   mpfr_t x;
   mpfr_exp_t emax;

   tests_start_mpfr ();

   mpfr_init (x);

   mpfr_set_nan (x);
   mpfr_prec_round (x, 2, MPFR_RNDN);
   MPFR_ASSERTN(mpfr_nan_p (x));

   mpfr_set_inf (x, 1);
   mpfr_prec_round (x, 2, MPFR_RNDN);
   MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);

   mpfr_set_inf (x, -1);
   mpfr_prec_round (x, 2, MPFR_RNDN);
   MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) < 0);

   mpfr_set_ui (x, 0, MPFR_RNDN);
   mpfr_prec_round (x, 2, MPFR_RNDN);
   MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS(x));

   mpfr_set_ui (x, 0, MPFR_RNDN);
   mpfr_neg (x, x, MPFR_RNDN);
   mpfr_prec_round (x, 2, MPFR_RNDN);
   MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_NEG(x));

   emax = mpfr_get_emax ();
   set_emax (0);
   mpfr_set_prec (x, 3);
   mpfr_set_str_binary (x, "0.111");
   mpfr_prec_round (x, 2, MPFR_RNDN);
   MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);
   set_emax (emax);

   mpfr_set_prec (x, mp_bits_per_limb + 2);
   mpfr_set_ui (x, 1, MPFR_RNDN);
   mpfr_nextbelow (x);
   mpfr_prec_round (x, mp_bits_per_limb + 1, MPFR_RNDN);
   MPFR_ASSERTN(mpfr_cmp_ui (x, 1) == 0);

   mpfr_set_prec (x, 3);
   mpfr_set_ui (x, 5, MPFR_RNDN);
   mpfr_prec_round (x, 2, MPFR_RNDN);
   if (mpfr_cmp_ui(x, 4))
     {
       printf ("Error in tround: got ");
       mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
       printf (" instead of 4\n");
       exit (1);
     }

   /* check case when reallocation is needed */
   mpfr_set_prec (x, 3);
   mpfr_set_ui (x, 5, MPFR_RNDN); /* exact */
   mpfr_prec_round (x, mp_bits_per_limb + 1, MPFR_RNDN);
   if (mpfr_cmp_ui(x, 5))
     {
       printf ("Error in tround: got ");
       mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
       printf (" instead of 5\n");
       exit (1);
     }

   mpfr_clear(x);
   mpfr_init2 (x, 3);
   mpfr_set_si (x, -5, MPFR_RNDN); /* exact */
   mpfr_prec_round (x, mp_bits_per_limb + 1, MPFR_RNDN);
   if (mpfr_cmp_si(x, -5))
     {
       printf ("Error in tround: got ");
       mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
       printf (" instead of -5\n");
       exit (1);
     }

   /* check case when new precision needs less limbs */
   mpfr_set_prec (x, mp_bits_per_limb + 1);
   mpfr_set_ui (x, 5, MPFR_RNDN); /* exact */
   mpfr_prec_round (x, 3, MPFR_RNDN); /* exact */
   if (mpfr_cmp_ui(x, 5))
     {
       printf ("Error in tround: got ");
       mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
       printf (" instead of 5\n");
       exit (1);
     }

   mpfr_clear(x);

   tests_end_mpfr ();
   return 0;
}

static void
check_special (void)
{
  mpfr_t  x, y, m;

  mpfr_init2 (x, 123L);
  mpfr_init2 (y, 123L);
  mpfr_init2 (m, 123L);

  /* agm(1,nan) is NaN */
  mpfr_set_ui (x, 1L, MPFR_RNDN);
  mpfr_set_nan (y);
  mpfr_agm (m, x, y, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (m));
  mpfr_agm (m, y, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (m));

  /* agm(1,+inf) == +inf */
  mpfr_set_ui (x, 1L, MPFR_RNDN);
  mpfr_set_inf (y, 1);
  mpfr_agm (m, x, y, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_inf_p (m));
  MPFR_ASSERTN (mpfr_sgn (m) > 0);
  mpfr_agm (m, y, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_inf_p (m));
  MPFR_ASSERTN (mpfr_sgn (m) > 0);

  /* agm(+inf,+inf) == +inf */
  mpfr_set_inf (x, 1);
  mpfr_set_inf (y, 1);
  mpfr_agm (m, x, y, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_inf_p (m));
  MPFR_ASSERTN (mpfr_sgn (m) > 0);

  /* agm(-inf,+inf) is NaN */
  mpfr_set_inf (x, -1);
  mpfr_set_inf (y, 1);
  mpfr_agm (m, x, y, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (m));
  mpfr_agm (m, y, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (m));

  /* agm(+0,+inf) is NaN */
  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_set_inf (y, 1);
  mpfr_agm (m, x, y, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (m));
  mpfr_agm (m, y, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (m));

  /* agm(-0,+inf) is NaN */
  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_neg (x, x, MPFR_RNDN);
  mpfr_set_inf (y, 1);
  mpfr_agm (m, x, y, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (m));
  mpfr_agm (m, y, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (m));

  /* agm(+0,-inf) is NaN */
  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_set_inf (y, -1);
  mpfr_agm (m, x, y, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (m));
  mpfr_agm (m, y, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (m));

  /* agm(-0,-inf) is NaN */
  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_neg (x, x, MPFR_RNDN);
  mpfr_set_inf (y, -1);
  mpfr_agm (m, x, y, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (m));
  mpfr_agm (m, y, x, MPFR_RNDN);
  MPFR_ASSERTN (mpfr_nan_p (m));

  /* agm(+0,1) == +0 */
  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_set_ui (y, 1, MPFR_RNDN);
  mpfr_agm (m, x, y, MPFR_RNDN);
  MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
  mpfr_agm (m, y, x, MPFR_RNDN);
  MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));

  /* agm(-0,1) == +0 */
  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_neg (x, x, MPFR_RNDN);
  mpfr_set_ui (y, 1, MPFR_RNDN);
  mpfr_agm (m, x, y, MPFR_RNDN);
  MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
  mpfr_agm (m, y, x, MPFR_RNDN);
  MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));

  /* agm(+0,-1) == +0 */
  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_set_si (y, -1, MPFR_RNDN);
  mpfr_agm (m, x, y, MPFR_RNDN);
  MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
  mpfr_agm (m, y, x, MPFR_RNDN);
  MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
//.........这里部分代码省略.........

static void
check_special (void)
{
  mpfr_t  a, d, q;
  mpfr_exp_t emax, emin;
  int i;

  mpfr_init2 (a, 100L);
  mpfr_init2 (d, 100L);
  mpfr_init2 (q, 100L);

  /* 1/nan == nan */
  mpfr_set_ui (a, 1L, MPFR_RNDN);
  MPFR_SET_NAN (d);
  mpfr_clear_flags ();
  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN);

  /* nan/1 == nan */
  MPFR_SET_NAN (a);
  mpfr_set_ui (d, 1L, MPFR_RNDN);
  mpfr_clear_flags ();
  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN);

  /* +inf/1 == +inf */
  MPFR_SET_INF (a);
  MPFR_SET_POS (a);
  mpfr_set_ui (d, 1L, MPFR_RNDN);
  mpfr_clear_flags ();
  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_inf_p (q));
  MPFR_ASSERTN (mpfr_sgn (q) > 0);
  MPFR_ASSERTN (__gmpfr_flags == 0);

  /* +inf/-1 == -inf */
  MPFR_SET_INF (a);
  MPFR_SET_POS (a);
  mpfr_set_si (d, -1, MPFR_RNDN);
  mpfr_clear_flags ();
  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_inf_p (q));
  MPFR_ASSERTN (mpfr_sgn (q) < 0);
  MPFR_ASSERTN (__gmpfr_flags == 0);

  /* -inf/1 == -inf */
  MPFR_SET_INF (a);
  MPFR_SET_NEG (a);
  mpfr_set_ui (d, 1L, MPFR_RNDN);
  mpfr_clear_flags ();
  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_inf_p (q));
  MPFR_ASSERTN (mpfr_sgn (q) < 0);
  MPFR_ASSERTN (__gmpfr_flags == 0);

  /* -inf/-1 == +inf */
  MPFR_SET_INF (a);
  MPFR_SET_NEG (a);
  mpfr_set_si (d, -1, MPFR_RNDN);
  mpfr_clear_flags ();
  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_inf_p (q));
  MPFR_ASSERTN (mpfr_sgn (q) > 0);
  MPFR_ASSERTN (__gmpfr_flags == 0);

  /* 1/+inf == +0 */
  mpfr_set_ui (a, 1L, MPFR_RNDN);
  MPFR_SET_INF (d);
  MPFR_SET_POS (d);
  mpfr_clear_flags ();
  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_number_p (q));
  MPFR_ASSERTN (mpfr_sgn (q) == 0);
  MPFR_ASSERTN (MPFR_IS_POS (q));
  MPFR_ASSERTN (__gmpfr_flags == 0);

  /* 1/-inf == -0 */
  mpfr_set_ui (a, 1L, MPFR_RNDN);
  MPFR_SET_INF (d);
  MPFR_SET_NEG (d);
  mpfr_clear_flags ();
  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_number_p (q));
  MPFR_ASSERTN (mpfr_sgn (q) == 0);
  MPFR_ASSERTN (MPFR_IS_NEG (q));
  MPFR_ASSERTN (__gmpfr_flags == 0);

  /* -1/+inf == -0 */
  mpfr_set_si (a, -1, MPFR_RNDN);
  MPFR_SET_INF (d);
  MPFR_SET_POS (d);
  mpfr_clear_flags ();
  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_number_p (q));
  MPFR_ASSERTN (mpfr_sgn (q) == 0);
  MPFR_ASSERTN (MPFR_IS_NEG (q));
  MPFR_ASSERTN (__gmpfr_flags == 0);

  /* -1/-inf == +0 */
  mpfr_set_si (a, -1, MPFR_RNDN);
//.........这里部分代码省略.........

static int
mpfr_all_div (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t r)
{
  mpfr_t a2;
  unsigned int oldflags, newflags;
  int inex, inex2;

  oldflags = __gmpfr_flags;
  inex = mpfr_div (a, b, c, r);

  if (a == b || a == c)
    return inex;

  newflags = __gmpfr_flags;

  mpfr_init2 (a2, MPFR_PREC (a));

  if (mpfr_integer_p (b) && ! (MPFR_IS_ZERO (b) && MPFR_IS_NEG (b)))
    {
      /* b is an integer, but not -0 (-0 is rejected as
         it becomes +0 when converted to an integer). */
      if (mpfr_fits_ulong_p (b, MPFR_RNDA))
        {
          __gmpfr_flags = oldflags;
          inex2 = mpfr_ui_div (a2, mpfr_get_ui (b, MPFR_RNDN), c, r);
          MPFR_ASSERTN (SAME_SIGN (inex2, inex));
          MPFR_ASSERTN (__gmpfr_flags == newflags);
          check_equal (a, a2, "mpfr_ui_div", b, c, r);
        }
      if (mpfr_fits_slong_p (b, MPFR_RNDA))
        {
          __gmpfr_flags = oldflags;
          inex2 = mpfr_si_div (a2, mpfr_get_si (b, MPFR_RNDN), c, r);
          MPFR_ASSERTN (SAME_SIGN (inex2, inex));
          MPFR_ASSERTN (__gmpfr_flags == newflags);
          check_equal (a, a2, "mpfr_si_div", b, c, r);
        }
    }

  if (mpfr_integer_p (c) && ! (MPFR_IS_ZERO (c) && MPFR_IS_NEG (c)))
    {
      /* c is an integer, but not -0 (-0 is rejected as
         it becomes +0 when converted to an integer). */
      if (mpfr_fits_ulong_p (c, MPFR_RNDA))
        {
          __gmpfr_flags = oldflags;
          inex2 = mpfr_div_ui (a2, b, mpfr_get_ui (c, MPFR_RNDN), r);
          MPFR_ASSERTN (SAME_SIGN (inex2, inex));
          MPFR_ASSERTN (__gmpfr_flags == newflags);
          check_equal (a, a2, "mpfr_div_ui", b, c, r);
        }
      if (mpfr_fits_slong_p (c, MPFR_RNDA))
        {
          __gmpfr_flags = oldflags;
          inex2 = mpfr_div_si (a2, b, mpfr_get_si (c, MPFR_RNDN), r);
          MPFR_ASSERTN (SAME_SIGN (inex2, inex));
          MPFR_ASSERTN (__gmpfr_flags == newflags);
          check_equal (a, a2, "mpfr_div_si", b, c, r);
        }
    }

  mpfr_clear (a2);

  return inex;
}

/* check concordance between mpfr_asprintf result with a regular mpfr float
   and with a regular double float */
static int
random_double (void)
{
    mpfr_t x; /* random regular mpfr float */
    double y; /* regular double float (equal to x) */

    char flag[] =
    {
        '-',
        '+',
        ' ',
        '#',
        '0', /* no ambiguity: first zeros are flag zero*/
        '\''
    };
    /* no 'a': mpfr and glibc do not have the same semantic */
    char specifier[] =
    {
        'e',
        'f',
        'g',
        'E',
        'f', /* SUSv2 doesn't accept %F, but %F and %f are the same for
              regular numbers */
        'G',
    };
    int spec; /* random index in specifier[] */
    int prec; /* random value for precision field */

    /* in the format string for mpfr_t variable, the maximum length is
       reached by something like "%-+ #0'.*Rf", that is 12 characters. */
#define FMT_MPFR_SIZE 12
    char fmt_mpfr[FMT_MPFR_SIZE];
    char *ptr_mpfr;

    /* in the format string for double variable, the maximum length is
       reached by something like "%-+ #0'.*f", that is 11 characters. */
#define FMT_SIZE 11
    char fmt[FMT_SIZE];
    char *ptr;

    int xi;
    char *xs;
    int yi;
    char *ys;

    int i, j, jmax;

    mpfr_init2 (x, MPFR_LDBL_MANT_DIG);

    for (i = 0; i < 1000; ++i)
    {
        /* 1. random double */
        do
        {
            y = DBL_RAND ();
        }
#ifdef HAVE_DENORMS
        while (0);
#else
        while (ABS(y) < DBL_MIN);
#endif

        if (randlimb () % 2 == 0)
            y = -y;

        mpfr_set_d (x, y, MPFR_RNDN);
        if (y != mpfr_get_d (x, MPFR_RNDN))
            /* conversion error: skip this one */
            continue;

        /* 2. build random format strings fmt_mpfr and fmt */
        ptr_mpfr = fmt_mpfr;
        ptr = fmt;
        *ptr_mpfr++ = *ptr++ = '%';
        /* random specifier 'e', 'f', 'g', 'E', 'F', or 'G' */
        spec = (int) (randlimb() % 6);
        /* random flags, but no ' flag with %e */
        jmax = (spec == 0 || spec == 3) ? 5 : 6;
        for (j = 0; j < jmax; j++)
        {
            if (randlimb() % 3 == 0)
                *ptr_mpfr++ = *ptr++ = flag[j];
        }
        *ptr_mpfr++ = *ptr++ = '.';
        *ptr_mpfr++ = *ptr++ = '*';
        *ptr_mpfr++ = 'R';
        *ptr_mpfr++ = *ptr++ = specifier[spec];
        *ptr_mpfr = *ptr = '\0';
        MPFR_ASSERTN (ptr - fmt < FMT_SIZE);
        MPFR_ASSERTN (ptr_mpfr - fmt_mpfr < FMT_MPFR_SIZE);

        /* advantage small precision */
        if (randlimb() % 2 == 0)
            prec = (int) (randlimb() % 10);
        else
            prec = (int) (randlimb() % prec_max_printf);

//.........这里部分代码省略.........

static void
special_atan2 (void)
{
    mpfr_t x, y, z;

    mpfr_inits2 (4, x, y, z, (mpfr_ptr) 0);

    /* Anything with NAN should be set to NAN */
    mpfr_set_ui (y, 0, MPFR_RNDN);
    mpfr_set_nan (x);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (MPFR_IS_NAN (z));
    mpfr_swap (x, y);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (MPFR_IS_NAN (z));

    /* 0+ 0+ --> 0+ */
    mpfr_set_ui (y, 0, MPFR_RNDN);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_POS (z));
    /* 0- 0+ --> 0- */
    MPFR_CHANGE_SIGN (y);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_NEG (z));
    /* 0- 0- --> -PI */
    MPFR_CHANGE_SIGN (x);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "-3.1415", 10, MPFR_RNDN) == 0);
    /* 0+ 0- --> +PI */
    MPFR_CHANGE_SIGN (y);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "3.1415", 10, MPFR_RNDN) == 0);
    /* 0+ -1 --> PI */
    mpfr_set_si (x, -1, MPFR_RNDN);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "3.1415", 10, MPFR_RNDN) == 0);
    /* 0- -1 --> -PI */
    MPFR_CHANGE_SIGN (y);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "-3.1415", 10, MPFR_RNDN) == 0);
    /* 0- +1 --> 0- */
    mpfr_set_ui (x, 1, MPFR_RNDN);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_NEG (z));
    /* 0+ +1 --> 0+ */
    MPFR_CHANGE_SIGN (y);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_POS (z));
    /* +1 0+ --> PI/2 */
    mpfr_swap (x, y);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "1.57075", 10, MPFR_RNDN) == 0);
    /* +1 0- --> PI/2 */
    MPFR_CHANGE_SIGN (x);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "1.57075", 10, MPFR_RNDN) == 0);
    /* -1 0- --> -PI/2 */
    MPFR_CHANGE_SIGN (y);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "-1.57075", 10, MPFR_RNDN) == 0);
    /* -1 0+ --> -PI/2 */
    MPFR_CHANGE_SIGN (x);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "-1.57075", 10, MPFR_RNDN) == 0);

    /* -1 +INF --> -0 */
    MPFR_SET_INF (x);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_NEG (z));
    /* +1 +INF --> +0 */
    MPFR_CHANGE_SIGN (y);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_POS (z));
    /* +1 -INF --> +PI */
    MPFR_CHANGE_SIGN (x);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "3.1415", 10, MPFR_RNDN) == 0);
    /* -1 -INF --> -PI */
    MPFR_CHANGE_SIGN (y);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "-3.1415", 10, MPFR_RNDN) == 0);
    /* -INF -1 --> -PI/2 */
    mpfr_swap (x, y);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "-1.57075", 10, MPFR_RNDN) == 0);
    /* +INF -1  --> PI/2 */
    MPFR_CHANGE_SIGN (y);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "1.57075", 10, MPFR_RNDN) == 0);
    /* +INF -INF --> 3*PI/4 */
    MPFR_SET_INF (x);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "2.356194490192344928", 10, MPFR_RNDN) == 0);
    /* +INF +INF --> PI/4 */
    MPFR_CHANGE_SIGN (x);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
    MPFR_ASSERTN (mpfr_cmp_str (z, "0.785375", 10, MPFR_RNDN) == 0);
    /* -INF +INF --> -PI/4 */
    MPFR_CHANGE_SIGN (y);
    mpfr_atan2 (z, y, x, MPFR_RNDN);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  mpfr_set_prec (y, 32);
  mpfr_set_si_2exp (x, -1, 200, GMP_RNDN);
  mpfr_sub_ui (x, x, 1, GMP_RNDN);
  mpfr_div_2ui (x, x, 1, GMP_RNDN);
  sign = -17;
  mpfr_lgamma (y, &sign, x, GMP_RNDN);
  mpfr_set_prec (x, 32);
  mpfr_set_str_binary (x, "-0.10001000111011111011000010100010E207");
  if (mpfr_equal_p (x, y) == 0 || sign != -1)
    {
      printf ("Error for lgamma(-2^199-0.5)\n");
      printf ("Got        ");
      mpfr_dump (y);
      printf ("with sign %d instead of ", sign);
      mpfr_dump (x);
      printf ("with sign -1.\n");
      exit (1);
    }

  mpfr_set_prec (x, 10);
  mpfr_set_prec (y, 10);
  mpfr_set_str_binary (x, "-0.1101111000E-3");
  inex = mpfr_lgamma (y, &sign, x, GMP_RNDN);
  mpfr_set_str_binary (x, "10.01001011");
  if (mpfr_equal_p (x, y) == 0 || sign != -1 || inex >= 0)
    {
      printf ("Error for lgamma(-0.1101111000E-3)\n");
      printf ("Got        ");
      mpfr_dump (y);
      printf ("instead of ");
      mpfr_dump (x);
      printf ("with sign %d instead of -1 (inex=%d).\n", sign, inex);
      exit (1);
    }

  mpfr_set_prec (x, 18);
  mpfr_set_prec (y, 28);
  mpfr_set_str_binary (x, "-1.10001101010001101e-196");
  inex = mpfr_lgamma (y, &sign, x, GMP_RNDN);
  mpfr_set_prec (x, 28);
  mpfr_set_str_binary (x, "0.100001110110101011011010011E8");
  MPFR_ASSERTN (mpfr_equal_p (x, y) && inex < 0);

  /* values reported by Kaveh Ghazi on 14 Jul 2007, where mpfr_lgamma()
     takes forever */
#define VAL1 "-0.11100001001010110111001010001001001011110100110000110E-55"
#define OUT1 "100110.01000000010111001110110101110101001001100110111"
#define VAL2 "-0.11100001001010110111001010001001001011110011111111100E-55"
#define OUT2 "100110.0100000001011100111011010111010100100110011111"
#define VAL3 "-0.11100001001010110111001010001001001001110101101010100E-55"
#define OUT3 "100110.01000000010111001110110101110101001011110111011"
#define VAL4 "-0.10001111110110110100100100000000001111110001001001011E-57"
#define OUT4 "101000.0001010111110011101101000101111111010001100011"
#define VAL5 "-0.10001111110110110100100100000000001111011111100001000E-57"
#define OUT5 "101000.00010101111100111011010001011111110100111000001"
#define VAL6 "-0.10001111110110110100100100000000001111011101100011001E-57"
#define OUT6 "101000.0001010111110011101101000101111111010011101111"

  mpfr_set_prec (x, 53);
  mpfr_set_prec (y, 53);

  mpfr_set_str_binary (x, VAL1);
  mpfr_lgamma (y, &sign, x, GMP_RNDN);
  mpfr_set_str_binary (x, OUT1);
  MPFR_ASSERTN(sign == -1 && mpfr_equal_p(x, y));

  mpfr_set_str_binary (x, VAL2);
  mpfr_lgamma (y, &sign, x, GMP_RNDN);
  mpfr_set_str_binary (x, OUT2);
  MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y));

  mpfr_set_str_binary (x, VAL3);
  mpfr_lgamma (y, &sign, x, GMP_RNDN);
  mpfr_set_str_binary (x, OUT3);
  MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y));

  mpfr_set_str_binary (x, VAL4);
  mpfr_lgamma (y, &sign, x, GMP_RNDN);
  mpfr_set_str_binary (x, OUT4);
  MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y));

  mpfr_set_str_binary (x, VAL5);
  mpfr_lgamma (y, &sign, x, GMP_RNDN);
  mpfr_set_str_binary (x, OUT5);
  MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y));

  mpfr_set_str_binary (x, VAL6);
  mpfr_lgamma (y, &sign, x, GMP_RNDN);
  mpfr_set_str_binary (x, OUT6);
  MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y));

  /* further test from Kaveh Ghazi */
  mpfr_set_str_binary (x, "-0.10011010101001010010001110010111010111011101010111001E-53");
  mpfr_lgamma (y, &sign, x, GMP_RNDN);
  mpfr_set_str_binary (x, "100101.00111101101010000000101010111010001111001101111");
  MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y));

  mpfr_clear (x);
  mpfr_clear (y);
}

static void
special (void)
{
  mpfr_t x, y;
  unsigned xprec, yprec;

  mpfr_init (x);
  mpfr_init (y);

  mpfr_set_prec (x, 32);
  mpfr_set_prec (y, 32);
  mpfr_set_ui (x, 1, MPFR_RNDN);
  mpfr_div_ui (y, x, 3, MPFR_RNDN);

  mpfr_set_prec (x, 100);
  mpfr_set_prec (y, 100);
  mpfr_urandomb (x, RANDS);
  mpfr_div_ui (y, x, 123456, MPFR_RNDN);
  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_div_ui (y, x, 123456789, MPFR_RNDN);
  if (mpfr_cmp_ui (y, 0))
    {
      printf ("mpfr_div_ui gives non-zero for 0/ui\n");
      exit (1);
    }

  /* bug found by Norbert Mueller, 21 Aug 2001 */
  mpfr_set_prec (x, 110);
  mpfr_set_prec (y, 60);
  mpfr_set_str_binary (x, "0.110101110011111110011111001110011001110111000000111110001000111011000011E-44");
  mpfr_div_ui (y, x, 17, MPFR_RNDN);
  mpfr_set_str_binary (x, "0.11001010100101100011101110000001100001010110101001010011011E-48");
  if (mpfr_cmp (x, y))
    {
      printf ("Error in x/17 for x=1/16!\n");
      printf ("Expected ");
      mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
      printf ("\nGot      ");
      mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
      printf ("\n");
      exit (1);
    }

  /* corner case */
  mpfr_set_prec (x, 2 * mp_bits_per_limb);
  mpfr_set_prec (y, 2);
  mpfr_set_ui (x, 4, MPFR_RNDN);
  mpfr_nextabove (x);
  mpfr_div_ui (y, x, 2, MPFR_RNDN); /* exactly in the middle */
  MPFR_ASSERTN(mpfr_cmp_ui (y, 2) == 0);

  mpfr_set_prec (x, 3 * mp_bits_per_limb);
  mpfr_set_prec (y, 2);
  mpfr_set_ui (x, 2, MPFR_RNDN);
  mpfr_nextabove (x);
  mpfr_div_ui (y, x, 2, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0);

  mpfr_set_prec (x, 3 * mp_bits_per_limb);
  mpfr_set_prec (y, 2);
  mpfr_set_si (x, -4, MPFR_RNDN);
  mpfr_nextbelow (x);
  mpfr_div_ui (y, x, 2, MPFR_RNDD);
  MPFR_ASSERTN(mpfr_cmp_si (y, -3) == 0);

  for (xprec = 53; xprec <= 128; xprec++)
    {
      mpfr_set_prec (x, xprec);
      mpfr_set_str_binary (x, "0.1100100100001111110011111000000011011100001100110111E2");
      for (yprec = 53; yprec <= 128; yprec++)
        {
          mpfr_set_prec (y, yprec);
          mpfr_div_ui (y, x, 1, MPFR_RNDN);
          if (mpfr_cmp(x,y))
            {
              printf ("division by 1.0 fails for xprec=%u, yprec=%u\n", xprec, yprec);
              printf ("expected "); mpfr_print_binary (x); puts ("");
              printf ("got      "); mpfr_print_binary (y); puts ("");
              exit (1);
            }
        }
    }

  /* Bug reported by Mark Dickinson, 6 Nov 2007 */
  mpfr_set_si (x, 0, MPFR_RNDN);
  mpfr_set_si (y, -1, MPFR_RNDN);
  mpfr_div_ui (y, x, 4, MPFR_RNDN);
  MPFR_ASSERTN(MPFR_IS_ZERO(y) && MPFR_IS_POS(y));

  mpfr_clear (x);
  mpfr_clear (y);
}