static void
set_special (mpfr_ptr x, unsigned int select)
{
  MPFR_ASSERTN (select < SPECIAL_MAX);
  switch (select)
    {
    case 0:
      MPFR_SET_NAN (x);
      break;
    case 1:
      MPFR_SET_INF (x);
      MPFR_SET_POS (x);
      break;
    case 2:
      MPFR_SET_INF (x);
      MPFR_SET_NEG (x);
      break;
    case 3:
      MPFR_SET_ZERO (x);
      MPFR_SET_POS  (x);
      break;
    case 4:
      MPFR_SET_ZERO (x);
      MPFR_SET_NEG  (x);
      break;
    case 5:
      mpfr_set_str_binary (x, "1");
      break;
    case 6:
      mpfr_set_str_binary (x, "-1");
      break;
    case 7:
      mpfr_set_str_binary (x, "1e-1");
      break;
    case 8:
      mpfr_set_str_binary (x, "1e+1");
      break;
    case 9:
      mpfr_const_pi (x, MPFR_RNDN);
      break;
    case 10:
      mpfr_const_pi (x, MPFR_RNDN);
      MPFR_SET_EXP (x, MPFR_GET_EXP (x)-1);
      break;
    default:
      mpfr_urandomb (x, RANDS);
      if (randlimb () & 1)
        mpfr_neg (x, x, MPFR_RNDN);
      break;
    }
}

void
mpfr_set_str_binary (mpfr_ptr x, const char *str)
{
  int has_sign;
  int res;

  if (*str == 'N')
    {
      MPFR_SET_NAN(x);
      __gmpfr_flags |= MPFR_FLAGS_NAN;
      return;
    }

  has_sign = *str == '-' || *str == '+';
  if (str[has_sign] == 'I')
    {
      MPFR_SET_INF(x);
      if (*str == '-')
        MPFR_SET_NEG(x);
      else
        MPFR_SET_POS(x);
      return;
    }

  res = mpfr_strtofr (x, str, 0, 2, MPFR_RNDZ);
  MPFR_ASSERTN (res == 0);
}

int
mpfr_sqrt_ui (mpfr_ptr r, unsigned long u, mpfr_rnd_t rnd_mode)
{
  if (u)
    {
      mpfr_t uu;
      mp_limb_t up[1];
      unsigned long cnt;
      int inex;
      MPFR_SAVE_EXPO_DECL (expo);

      MPFR_TMP_INIT1 (up, uu, GMP_NUMB_BITS);
      MPFR_ASSERTN (u == (mp_limb_t) u);
      count_leading_zeros (cnt, (mp_limb_t) u);
      *up = (mp_limb_t) u << cnt;

      MPFR_SAVE_EXPO_MARK (expo);
      MPFR_SET_EXP (uu, GMP_NUMB_BITS - cnt);
      inex = mpfr_sqrt(r, uu, rnd_mode);
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_check_range(r, inex, rnd_mode);
    }
  else /* sqrt(0) = 0 */
    {
      MPFR_SET_ZERO(r);
      MPFR_SET_POS(r);
      MPFR_RET(0);
    }
}

int
mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate)
{
  mp_ptr rp;
  mp_prec_t nbits;
  mp_size_t nlimbs;
  mp_size_t k; /* number of high zero limbs */
  mp_exp_t exp;
  int cnt;

  MPFR_CLEAR_FLAGS (rop);

  rp = MPFR_MANT (rop);
  nbits = MPFR_PREC (rop);
  nlimbs = MPFR_LIMB_SIZE (rop);
  MPFR_SET_POS (rop);

  /* Uniform non-normalized significand */
  _gmp_rand (rp, rstate, nlimbs * BITS_PER_MP_LIMB);

  /* If nbits isn't a multiple of BITS_PER_MP_LIMB, mask the low bits */
  cnt = nlimbs * BITS_PER_MP_LIMB - nbits;
  if (MPFR_LIKELY (cnt != 0))
    rp[0] &= ~MPFR_LIMB_MASK (cnt);

  /* Count the null significant limbs and remaining limbs */
  exp = 0;
  k = 0;
  while (nlimbs != 0 && rp[nlimbs - 1] == 0)
    {
      k ++;
      nlimbs --;
      exp -= BITS_PER_MP_LIMB;
    }

  if (MPFR_LIKELY (nlimbs != 0)) /* otherwise value is zero */
    {
      count_leading_zeros (cnt, rp[nlimbs - 1]);
      /* Normalization */
      if (mpfr_set_exp (rop, exp - cnt))
        {
          /* If the exponent is not in the current exponent range, we
             choose to return a NaN as this is probably a user error.
             Indeed this can happen only if the exponent range has been
             reduced to a very small interval and/or the precision is
             huge (very unlikely). */
          MPFR_SET_NAN (rop);
          __gmpfr_flags |= MPFR_FLAGS_NAN; /* Can't use MPFR_RET_NAN */
          return 1;
        }
      if (cnt != 0)
        mpn_lshift (rp + k, rp, nlimbs, cnt);
      if (k != 0)
        MPN_ZERO (rp, k);
    }
  else
    MPFR_SET_ZERO (rop);

  return 0;
}

static void
test_overflow1 (void)
{
  mpfr_t x, y, z, r;
  int inex;

  mpfr_inits2 (8, x, y, z, r, (void *) 0);
  MPFR_SET_POS (x);
  mpfr_setmax (x, mpfr_get_emax ());  /* x = 2^emax - ulp */
  mpfr_set_ui (y, 2, GMP_RNDN);       /* y = 2 */
  mpfr_neg (z, x, GMP_RNDN);          /* z = -x = -(2^emax - ulp) */
  mpfr_clear_flags ();
  /* The intermediate multiplication x * y overflows, but x * y + z = x
     is representable. */
  inex = mpfr_fma (r, x, y, z, GMP_RNDN);
  if (inex || ! mpfr_equal_p (r, x))
    {
      printf ("Error in test_overflow1\nexpected ");
      mpfr_out_str (stdout, 2, 0, x, GMP_RNDN);
      printf (" with inex = 0\n     got ");
      mpfr_out_str (stdout, 2, 0, r, GMP_RNDN);
      printf (" with inex = %d\n", inex);
      exit (1);
    }
  if (mpfr_overflow_p ())
    {
      printf ("Error in test_overflow1: overflow flag set\n");
      exit (1);
    }
  mpfr_clears (x, y, z, r, (void *) 0);
}

int
mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate)
{
    mpfr_limb_ptr rp;
    mpfr_prec_t nbits;
    mp_size_t nlimbs;
    mp_size_t k; /* number of high zero limbs */
    mpfr_exp_t exp;
    int cnt;

    rp = MPFR_MANT (rop);
    nbits = MPFR_PREC (rop);
    nlimbs = MPFR_LIMB_SIZE (rop);
    MPFR_SET_POS (rop);
    cnt = nlimbs * GMP_NUMB_BITS - nbits;

    /* Uniform non-normalized significand */
    /* generate exactly nbits so that the random generator stays in the same
       state, independent of the machine word size GMP_NUMB_BITS */
    mpfr_rand_raw (rp, rstate, nbits);
    if (MPFR_LIKELY (cnt != 0)) /* this will put the low bits to zero */
        mpn_lshift (rp, rp, nlimbs, cnt);

    /* Count the null significant limbs and remaining limbs */
    exp = 0;
    k = 0;
    while (nlimbs != 0 && rp[nlimbs - 1] == 0)
    {
        k ++;
        nlimbs --;
        exp -= GMP_NUMB_BITS;
    }

    if (MPFR_LIKELY (nlimbs != 0)) /* otherwise value is zero */
    {
        count_leading_zeros (cnt, rp[nlimbs - 1]);
        /* Normalization */
        if (mpfr_set_exp (rop, exp - cnt))
        {
            /* If the exponent is not in the current exponent range, we
               choose to return a NaN as this is probably a user error.
               Indeed this can happen only if the exponent range has been
               reduced to a very small interval and/or the precision is
               huge (very unlikely). */
            MPFR_SET_NAN (rop);
            __gmpfr_flags |= MPFR_FLAGS_NAN; /* Can't use MPFR_RET_NAN */
            return 1;
        }
        if (cnt != 0)
            mpn_lshift (rp + k, rp, nlimbs, cnt);
        if (k != 0)
            MPN_ZERO (rp, k);
    }
    else
        MPFR_SET_ZERO (rop);

    return 0;
}

static void
check_nan (void)
{
  mpfr_t  d, q;

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

  /* 1/+inf == 0 */
  MPFR_CLEAR_FLAGS (d);
  MPFR_SET_INF (d);
  MPFR_SET_POS (d);
  MPFR_ASSERTN (mpfr_ui_div (q, 1L, d, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_number_p (q));
  MPFR_ASSERTN (mpfr_sgn (q) == 0);

  /* 1/-inf == -0 */
  MPFR_CLEAR_FLAGS (d);
  MPFR_SET_INF (d);
  MPFR_SET_NEG (d);
  MPFR_ASSERTN (mpfr_ui_div (q, 1L, d, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_number_p (q));
  MPFR_ASSERTN (mpfr_sgn (q) == 0);

  /* 1/nan == nan */
  MPFR_SET_NAN (d);
  MPFR_ASSERTN (mpfr_ui_div (q, 1L, d, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_nan_p (q));

  /* 0/0 == nan */
  mpfr_set_ui (d, 0L, GMP_RNDN);
  MPFR_ASSERTN (mpfr_ui_div (q, 0L, d, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_nan_p (q));

  /* 1/+0 = +inf */
  mpfr_set_ui (d, 0L, GMP_RNDN);
  MPFR_ASSERTN (mpfr_ui_div (q, 1L, d, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0);

  /* 1/-0 = -inf */
  mpfr_set_ui (d, 0L, GMP_RNDN);
  mpfr_neg (d, d, GMP_RNDN);
  MPFR_ASSERTN (mpfr_ui_div (q, 1L, d, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0);

  /* 0/1 = +0 */
  mpfr_set_ui (d, 1L, GMP_RNDN);
  MPFR_ASSERTN (mpfr_ui_div (q, 0L, d, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_cmp_ui (q, 0) == 0 && MPFR_IS_POS (q));

  /* 0/-1 = -0 */
  mpfr_set_si (d, -1, GMP_RNDN);
  MPFR_ASSERTN (mpfr_ui_div (q, 0L, d, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_cmp_ui (q, 0) == 0 && MPFR_IS_NEG (q));

  mpfr_clear (d);
  mpfr_clear (q);
}

void
mpfr_set_inf (mpfr_ptr x, int sign)
{
  MPFR_SET_INF(x);
  if (sign >= 0)
    MPFR_SET_POS(x);
  else
    MPFR_SET_NEG(x);
}

int
mpfr_set_si (mpfr_ptr x, long i, mp_rnd_t rnd_mode)
{
  int inex;
  mp_size_t xn;
  unsigned int cnt, nbits;
  mp_limb_t ai, *xp;

  MPFR_CLEAR_FLAGS(x);
  if (i == 0)
    {
      MPFR_SET_ZERO(x);
      MPFR_SET_POS(x);
      MPFR_RET(0);
    }

  xn = (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB;
  ai = SAFE_ABS(long, i);
  count_leading_zeros(cnt, ai);

  xp = MPFR_MANT(x);
  xp[xn] = ai << cnt;
  /* don't forget to put zero in lower limbs */
  MPN_ZERO(xp, xn);
  /* set sign */
  if ((i < 0) ^ (MPFR_SIGN(x) < 0))
    MPFR_CHANGE_SIGN(x);

  MPFR_EXP(x) = nbits = BITS_PER_MP_LIMB - cnt;
  inex = mpfr_check_range(x, rnd_mode);
  if (inex)
    return inex; /* underflow or overflow */

  /* round if MPFR_PREC(x) smaller than length of i */
  if (MPFR_PREC(x) < nbits)
    {
      int carry;

      carry = mpfr_round_raw(xp+xn, xp+xn, nbits, (i < 0), MPFR_PREC(x),
                             rnd_mode, &inex);
      if (carry)
        {
          mp_exp_t exp = MPFR_EXP(x);

          if (exp == __mpfr_emax)
            return mpfr_set_overflow(x, rnd_mode, (i < 0 ? -1 : 1));

          MPFR_EXP(x)++;
          xp[xn] = GMP_LIMB_HIGHBIT;
        }
    }

  MPFR_RET(inex);
}

static int
my_setstr (mpfr_ptr t, const char *s)
{
  if (strcmp (s, "min") == 0)
    {
      mpfr_setmin (t, mpfr_get_emin ());
      MPFR_SET_POS (t);
      return 0;
    }
  if (strcmp (s, "min+") == 0)
    {
      mpfr_setmin (t, mpfr_get_emin ());
      MPFR_SET_POS (t);
      mpfr_nextabove (t);
      return 0;
    }
  if (strcmp (s, "max") == 0)
    {
      mpfr_setmax (t, mpfr_get_emax ());
      MPFR_SET_POS (t);
      return 0;
    }
  return mpfr_set_str (t, s, 10, MPFR_RNDN);
}

/* set f to the integer z */
int 
mpfr_set_z (mpfr_ptr f, mpz_srcptr z, mp_rnd_t rnd_mode)
{
  mp_size_t fn, zn, dif;
  int k, sign_z, inex;
  mp_limb_t *fp, *zp;
  mp_exp_t exp;

  MPFR_CLEAR_FLAGS (f); /* z cannot be NaN nor Inf */

  sign_z = mpz_cmp_ui (z, 0);

  if (sign_z == 0)
    {
      MPFR_SET_ZERO(f);
      MPFR_SET_POS(f);
      MPFR_RET(0);
    }

  fp = MPFR_MANT(f);
  fn = 1 + (MPFR_PREC(f) - 1) / BITS_PER_MP_LIMB;
  zn = ABS(SIZ(z));
  dif = zn - fn;
  zp = PTR(z);
  count_leading_zeros(k, zp[zn-1]);

  exp = (mp_prec_t) zn * BITS_PER_MP_LIMB - k;
  /* The exponent will be exp or exp + 1 (due to rounding) */
  if (exp > __mpfr_emax)
    return mpfr_set_overflow(f, rnd_mode, sign_z);
  if (exp + 1 < __mpfr_emin)
    return mpfr_set_underflow(f, rnd_mode, sign_z);

  if (MPFR_SIGN(f) * sign_z < 0)
    MPFR_CHANGE_SIGN(f);

  if (dif >= 0)
    {
      mp_limb_t cc;
      int sh;

      /* number has to be truncated */
      if (k != 0)
        {
          mpn_lshift(fp, zp + dif, fn, k);
          if (dif != 0)
            fp[0] += zp[dif - 1] >> (BITS_PER_MP_LIMB - k);
        }

static void
check_singular (void)
{
  mpfr_t  x, got;

  mpfr_init2 (x, 100L);
  mpfr_init2 (got, 100L);

  /* sqrt(NaN) == NaN */
  MPFR_SET_NAN (x);
  MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_nan_p (got));

  /* sqrt(-1) == NaN */
  mpfr_set_si (x, -1L, MPFR_RNDZ);
  MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_nan_p (got));

  /* sqrt(+inf) == +inf */
  MPFR_SET_INF (x);
  MPFR_SET_POS (x);
  MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_inf_p (got));

  /* sqrt(-inf) == NaN */
  MPFR_SET_INF (x);
  MPFR_SET_NEG (x);
  MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_nan_p (got));

  /* sqrt(+0) == +0 */
  mpfr_set_si (x, 0L, MPFR_RNDZ);
  MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_number_p (got));
  MPFR_ASSERTN (mpfr_cmp_ui (got, 0L) == 0);
  MPFR_ASSERTN (MPFR_IS_POS (got));

  /* sqrt(-0) == -0 */
  mpfr_set_si (x, 0L, MPFR_RNDZ);
  MPFR_SET_NEG (x);
  MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_number_p (got));
  MPFR_ASSERTN (mpfr_cmp_ui (got, 0L) == 0);
  MPFR_ASSERTN (MPFR_IS_NEG (got));

  mpfr_clear (x);
  mpfr_clear (got);
}

int
mpfr_set_si_2exp (mpfr_ptr x, long i, mp_exp_t e, mp_rnd_t rnd_mode)
{
  if (i == 0)
    {
      MPFR_SET_ZERO (x);
      MPFR_SET_POS (x);
      MPFR_RET (0);
    }
  else
    {
      mp_size_t xn;
      unsigned int cnt, nbits;
      mp_limb_t ai, *xp;
      int inex = 0;

      /* FIXME: support int limbs (e.g. 16-bit limbs on 16-bit proc) */
      ai = SAFE_ABS (unsigned long, i);
      MPFR_ASSERTN (SAFE_ABS (unsigned long, i) == ai);

      /* Position of the highest limb */
      xn = (MPFR_PREC (x) - 1) / BITS_PER_MP_LIMB;
      count_leading_zeros (cnt, ai);
      MPFR_ASSERTD (cnt < BITS_PER_MP_LIMB);  /* OK since i != 0 */

      xp = MPFR_MANT(x);
      xp[xn] = ai << cnt;
      /* Zero the xn lower limbs. */
      MPN_ZERO(xp, xn);
      MPFR_SET_SIGN (x, i < 0 ? MPFR_SIGN_NEG : MPFR_SIGN_POS);

      nbits = BITS_PER_MP_LIMB - cnt;
      e += nbits;  /* exponent _before_ the rounding */

      /* round if MPFR_PREC(x) smaller than length of i */
      if (MPFR_UNLIKELY (MPFR_PREC (x) < nbits) &&
          MPFR_UNLIKELY (mpfr_round_raw (xp + xn, xp + xn, nbits, i < 0,
                                         MPFR_PREC (x), rnd_mode, &inex)))
        {
          e++;
          xp[xn] = MPFR_LIMB_HIGHBIT;
        }

      MPFR_CLEAR_FLAGS (x);
      MPFR_EXP (x) = e;
      return mpfr_check_range (x, inex, rnd_mode);
    }
}

/* set f to the rational q */
int
mpfr_set_q (mpfr_ptr f, mpq_srcptr q, mp_rnd_t rnd)
{
  mpz_srcptr num, den;
  mpfr_t n, d;
  int inexact;
  mp_prec_t prec;

  MPFR_CLEAR_FLAGS (f);
  num = mpq_numref (q);
  if (mpz_cmp_ui (num, 0) == 0)
    {
      MPFR_SET_ZERO (f);
      MPFR_SET_POS (f);
      MPFR_RET (0);
    }

  den = mpq_denref (q);
  mpfr_save_emin_emax ();
  prec = mpz_sizeinbase (num, 2);
  if (prec < MPFR_PREC_MIN)
    prec = MPFR_PREC_MIN;
  mpfr_init2 (n, prec);
  if (mpfr_set_z (n, num, GMP_RNDZ)) /* result is exact unless overflow */
    {
      mpfr_clear (n);
      mpfr_restore_emin_emax ();
      MPFR_SET_NAN (f);
      MPFR_RET_NAN;
    }
  prec = mpz_sizeinbase(den, 2);
  if (prec < MPFR_PREC_MIN)
    prec = MPFR_PREC_MIN;
  mpfr_init2 (d, prec);
  if (mpfr_set_z (d, den, GMP_RNDZ)) /* result is exact unless overflow */
    {
      mpfr_clear (d);
      mpfr_clear (n);
      mpfr_restore_emin_emax ();
      MPFR_SET_NAN (f);
      MPFR_RET_NAN;
    }
  inexact = mpfr_div (f, n, d, rnd);
  mpfr_clear (n);
  mpfr_clear (d);
  MPFR_RESTORE_RET (inexact, f, rnd);
}

static void
check_sgn(void)
{
  mpfr_t x;
  int i, s1, s2;

  mpfr_init(x);
  for(i = 0 ; i < 100 ; i++)
    {
      mpfr_urandomb (x, RANDS);
      if (i&1)
        {
          MPFR_SET_POS(x);
          s2 = 1;
        }
      else
        {
          MPFR_SET_NEG(x);
          s2 = -1;
        }
      s1 = mpfr_sgn(x);
      if (s1 < -1 || s1 > 1)
        {
          printf("Error for sgn: out of range.\n");
          goto lexit;
        }
      else if (MPFR_IS_NAN(x) || MPFR_IS_ZERO(x))
        {
          if (s1 != 0)
            {
              printf("Error for sgn: Nan or Zero should return 0.\n");
              goto lexit;
            }
        }
      else if (s1 != s2)
        {
          printf("Error for sgn. Return %d instead of %d.\n", s1, s2);
          goto lexit;
        }
    }
  mpfr_clear(x);
  return;

 lexit:
  mpfr_clear(x);
  exit(1);
}

int
mpfr_set_ui_2exp (mpfr_ptr x, unsigned long i, mpfr_exp_t e, mpfr_rnd_t rnd_mode)
{
  MPFR_SET_POS (x);

  if (i == 0)
    {
      MPFR_SET_ZERO (x);
      MPFR_RET (0);
    }
  else
    {
      mp_size_t xn;
      unsigned int cnt, nbits;
      mp_limb_t *xp;
      int inex = 0;

      /* FIXME: support int limbs (e.g. 16-bit limbs on 16-bit proc) */
      MPFR_ASSERTD (i == (mp_limb_t) i);

      /* Position of the highest limb */
      xn = (MPFR_PREC (x) - 1) / GMP_NUMB_BITS;
      count_leading_zeros (cnt, (mp_limb_t) i);
      MPFR_ASSERTD (cnt < GMP_NUMB_BITS);  /* OK since i != 0 */

      xp = MPFR_MANT(x);
      xp[xn] = ((mp_limb_t) i) << cnt;
      /* Zero the xn lower limbs. */
      MPN_ZERO(xp, xn);

      nbits = GMP_NUMB_BITS - cnt;
      e += nbits;  /* exponent _before_ the rounding */

      /* round if MPFR_PREC(x) smaller than length of i */
      if (MPFR_UNLIKELY (MPFR_PREC (x) < nbits) &&
          MPFR_UNLIKELY (mpfr_round_raw (xp + xn, xp + xn, nbits, 0,
                                         MPFR_PREC (x), rnd_mode, &inex)))
        {
          e++;
          xp[xn] = MPFR_LIMB_HIGHBIT;
        }

      MPFR_EXP (x) = e;
      return mpfr_check_range (x, inex, rnd_mode);
    }
}

int
mpfr_dim (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
{
  if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y))
    {
      MPFR_SET_NAN(z);
      MPFR_RET_NAN;
    }

  if (mpfr_cmp (x,y) > 0)
    return mpfr_sub (z, x, y, rnd_mode);
  else
    {
      MPFR_SET_ZERO(z);
      MPFR_SET_POS(z);
      MPFR_RET(0);
    }
}

/*
   Parameters:
   s - the input floating-point number
   n, p - parameters from the algorithm
   tc - an array of p floating-point numbers tc[1]..tc[p]
   Output:
   b is the result, i.e.
   sum(tc[i]*product((s+2j)*(s+2j-1)/n^2,j=1..i-1), i=1..p)*s*n^(-s-1)
*/
static void
mpfr_zeta_part_b (mpfr_t b, mpfr_srcptr s, int n, int p, mpfr_t *tc)
{
  mpfr_t s1, d, u;
  unsigned long n2;
  int l, t;
  MPFR_GROUP_DECL (group);

  if (p == 0)
    {
      MPFR_SET_ZERO (b);
      MPFR_SET_POS (b);
      return;
    }

  n2 = n * n;
  MPFR_GROUP_INIT_3 (group, MPFR_PREC (b), s1, d, u);

  /* t equals 2p-2, 2p-3, ... ; s1 equals s+t */
  t = 2 * p - 2;
  mpfr_set (d, tc[p], GMP_RNDN);
  for (l = 1; l < p; l++)
    {
      mpfr_add_ui (s1, s, t, GMP_RNDN); /* s + (2p-2l) */
      mpfr_mul (d, d, s1, GMP_RNDN);
      t = t - 1;
      mpfr_add_ui (s1, s, t, GMP_RNDN); /* s + (2p-2l-1) */
      mpfr_mul (d, d, s1, GMP_RNDN);
      t = t - 1;
      mpfr_div_ui (d, d, n2, GMP_RNDN);
      mpfr_add (d, d, tc[p-l], GMP_RNDN);
      /* since s is positive and the tc[i] have alternate signs,
         the following is unlikely */
      if (MPFR_UNLIKELY (mpfr_cmpabs (d, tc[p-l]) > 0))
        mpfr_set (d, tc[p-l], GMP_RNDN);
    }
  mpfr_mul (d, d, s, GMP_RNDN);
  mpfr_add (s1, s, __gmpfr_one, GMP_RNDN);
  mpfr_neg (s1, s1, GMP_RNDN);
  mpfr_ui_pow (u, n, s1, GMP_RNDN);
  mpfr_mul (b, d, u, GMP_RNDN);

  MPFR_GROUP_CLEAR (group);
}

static void
check_special (void)
{
  mpfr_t x;
  int ret = 0;

  mpfr_init (x);
  MPFR_SET_ZERO (x);
  if ((mpfr_sgn) (x) != 0 || mpfr_sgn (x) != 0)
    {
      printf("Sgn error for 0.\n");
      ret = 1;
    }
  MPFR_SET_INF (x);
  MPFR_SET_POS (x);
  if ((mpfr_sgn) (x) != 1 || mpfr_sgn (x) != 1)
    {
      printf("Sgn error for +Inf.\n");
      ret = 1;
    }
  MPFR_SET_INF (x);
  MPFR_SET_NEG (x);
  if ((mpfr_sgn) (x) != -1 || mpfr_sgn (x) != -1)
    {
      printf("Sgn error for -Inf.\n");
      ret = 1;
    }
  MPFR_SET_NAN (x);
  mpfr_clear_flags ();
  if ((mpfr_sgn) (x) != 0 || !mpfr_erangeflag_p ())
    {
      printf("Sgn error for NaN.\n");
      ret = 1;
    }
  mpfr_clear_flags ();
  if (mpfr_sgn (x) != 0 || !mpfr_erangeflag_p ())
    {
      printf("Sgn error for NaN.\n");
      ret = 1;
    }
  mpfr_clear (x);
  if (ret)
    exit (ret);
}

int
mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_exp_t expx;
  mpfr_prec_t precy;
  int inexact;
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd_mode),
                 ("y[%#R]=%R inexact=%d", y, y, inexact));

  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 (MPFR_IS_POS(x))
            MPFR_SET_INF(y);
          else
            MPFR_SET_ZERO(y);
          MPFR_SET_POS(y);
          MPFR_RET(0);
        }
      else
        {
          MPFR_ASSERTD(MPFR_IS_ZERO(x));
          return mpfr_set_ui (y, 1, rnd_mode);
        }
    }

  /* First, let's detect most overflow and underflow cases. */
  {
    mpfr_t e, bound;

    /* We must extended the exponent range and save the flags now. */
    MPFR_SAVE_EXPO_MARK (expo);

    mpfr_init2 (e, sizeof (mpfr_exp_t) * CHAR_BIT);
    mpfr_init2 (bound, 32);

    inexact = mpfr_set_exp_t (e, expo.saved_emax, MPFR_RNDN);
    MPFR_ASSERTD (inexact == 0);
    mpfr_const_log2 (bound, expo.saved_emax < 0 ? MPFR_RNDD : MPFR_RNDU);
    mpfr_mul (bound, bound, e, MPFR_RNDU);
    if (MPFR_UNLIKELY (mpfr_cmp (x, bound) >= 0))
      {
        /* x > log(2^emax), thus exp(x) > 2^emax */
        mpfr_clears (e, bound, (mpfr_ptr) 0);
        MPFR_SAVE_EXPO_FREE (expo);
        return mpfr_overflow (y, rnd_mode, 1);
      }

    inexact = mpfr_set_exp_t (e, expo.saved_emin, MPFR_RNDN);
    MPFR_ASSERTD (inexact == 0);
    inexact = mpfr_sub_ui (e, e, 2, MPFR_RNDN);
    MPFR_ASSERTD (inexact == 0);
    mpfr_const_log2 (bound, expo.saved_emin < 0 ? MPFR_RNDU : MPFR_RNDD);
    mpfr_mul (bound, bound, e, MPFR_RNDD);
    if (MPFR_UNLIKELY (mpfr_cmp (x, bound) <= 0))
      {
        /* x < log(2^(emin - 2)), thus exp(x) < 2^(emin - 2) */
        mpfr_clears (e, bound, (mpfr_ptr) 0);
        MPFR_SAVE_EXPO_FREE (expo);
        return mpfr_underflow (y, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode,
                               1);
      }

    /* Other overflow/underflow cases must be detected
       by the generic routines. */
    mpfr_clears (e, bound, (mpfr_ptr) 0);
    MPFR_SAVE_EXPO_FREE (expo);
  }

  expx  = MPFR_GET_EXP (x);
  precy = MPFR_PREC (y);

  /* if x < 2^(-precy), then exp(x) i.e. gives 1 +/- 1 ulp(1) */
  if (MPFR_UNLIKELY (expx < 0 && (mpfr_uexp_t) (-expx) > precy))
    {
      mpfr_exp_t emin = __gmpfr_emin;
      mpfr_exp_t emax = __gmpfr_emax;
      int signx = MPFR_SIGN (x);

      MPFR_SET_POS (y);
      if (MPFR_IS_NEG_SIGN (signx) && (rnd_mode == MPFR_RNDD ||
                                       rnd_mode == MPFR_RNDZ))
        {
          __gmpfr_emin = 0;
          __gmpfr_emax = 0;
          mpfr_setmax (y, 0);  /* y = 1 - epsilon */
          inexact = -1;
        }
      else
        {
          __gmpfr_emin = 1;
          __gmpfr_emax = 1;
//.........这里部分代码省略.........