/* Convert an mpz_t to an mpfr_exp_t, restricted to
   the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]. */
mpfr_exp_t
mpfr_ubf_zexp2exp (mpz_ptr ez)
{
  mp_size_t n;
  mpfr_eexp_t e;
  mpfr_t d;
  int inex;
  MPFR_SAVE_EXPO_DECL (expo);

  n = ABSIZ (ez); /* limb size of ez */
  if (n == 0)
    return 0;

  MPFR_SAVE_EXPO_MARK (expo);
  mpfr_init2 (d, n * GMP_NUMB_BITS);
  MPFR_DBGRES (inex = mpfr_set_z (d, ez, MPFR_RNDN));
  MPFR_ASSERTD (inex == 0);
  e = mpfr_get_exp_t (d, MPFR_RNDZ);
  mpfr_clear (d);
  MPFR_SAVE_EXPO_FREE (expo);
  if (MPFR_UNLIKELY (e < MPFR_EXP_MIN))
    return MPFR_EXP_MIN;
  if (MPFR_UNLIKELY (e > MPFR_EXP_MAX))
    return MPFR_EXP_MAX;
  return e;
}

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_add_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mpfr_rnd_t rnd_mode)
{
  MPFR_LOG_FUNC
    (("x[%Pu]=%.*Rg u=%lu rnd=%d",
      mpfr_get_prec(x), mpfr_log_prec, x, u, rnd_mode),
     ("y[%Pu]=%.*Rg", mpfr_get_prec (y), mpfr_log_prec, y));

  if (MPFR_LIKELY(u != 0) )  /* if u=0, do nothing */
    {
      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_ASSERTD (u == (mp_limb_t) u);
      count_leading_zeros(cnt, (mp_limb_t) u);
      up[0] = (mp_limb_t) u << cnt;

      /* Optimization note: Exponent save/restore operations may be
         removed if mpfr_add works even when uu is out-of-range. */
      MPFR_SAVE_EXPO_MARK (expo);
      MPFR_SET_EXP (uu, GMP_NUMB_BITS - cnt);
      inex = mpfr_add(y, x, uu, rnd_mode);
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_check_range(y, inex, rnd_mode);
    }
  else
    /* (unsigned long) 0 is assumed to be a real 0 (unsigned) */
    return mpfr_set (y, x, rnd_mode);
}

/* This function does not change the flags. */
static void
mpfr_get_zexp (mpz_ptr ez, mpfr_srcptr x)
{
  mpz_init (ez);

  if (MPFR_IS_UBF (x))
    mpz_set (ez, MPFR_ZEXP (x));
  else
    {
      mp_limb_t e_limb[MPFR_EXP_LIMB_SIZE];
      mpfr_t e;
      int inex;
      MPFR_SAVE_EXPO_DECL (expo);

      /* TODO: Once this has been tested, optimize based on whether
         _MPFR_EXP_FORMAT <= 3. */
      MPFR_TMP_INIT1 (e_limb, e, sizeof (mpfr_exp_t) * CHAR_BIT);
      MPFR_SAVE_EXPO_MARK (expo);
      MPFR_DBGRES (inex = mpfr_set_exp_t (e, MPFR_GET_EXP (x), MPFR_RNDN));
      MPFR_ASSERTD (inex == 0);
      MPFR_DBGRES (inex = mpfr_get_z (ez, e, MPFR_RNDN));
      MPFR_ASSERTD (inex == 0);
      MPFR_SAVE_EXPO_FREE (expo);
    }
}

int
mpfr_mul_d (mpfr_ptr a, mpfr_srcptr b, double c, mpfr_rnd_t rnd_mode)
{
  int inexact;
  mpfr_t d;
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_LOG_FUNC
    (("b[%Pu]=%.*Rg c=%.20g rnd=%d",
      mpfr_get_prec(b), mpfr_log_prec, b, c, rnd_mode),
     ("a[%Pu]=%.*Rg inexact=%d",
      mpfr_get_prec (a), mpfr_get_prec, a, inexact));

  MPFR_SAVE_EXPO_MARK (expo);

  mpfr_init2 (d, IEEE_DBL_MANT_DIG);
  inexact = mpfr_set_d (d, c, rnd_mode);
  MPFR_ASSERTN (inexact == 0);

  mpfr_clear_flags ();
  inexact = mpfr_mul (a, b, d, rnd_mode);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);

  mpfr_clear(d);
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (a, inexact, rnd_mode);
}

int
mpfr_d_div (mpfr_ptr a, double b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
  int inexact;
  mpfr_t d;
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_LOG_FUNC (
    ("b=%.20g c[%Pu]=%*.Rg rnd=%d", b, mpfr_get_prec (c), mpfr_log_prec, c, rnd_mode),
    ("a[%Pu]=%*.Rg", mpfr_get_prec (a), mpfr_log_prec, a));

  MPFR_SAVE_EXPO_MARK (expo);

  mpfr_init2 (d, IEEE_DBL_MANT_DIG);
  inexact = mpfr_set_d (d, b, rnd_mode);
  MPFR_ASSERTN (inexact == 0);

  mpfr_clear_flags ();
  inexact = mpfr_div (a, d, c, rnd_mode);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);

  mpfr_clear(d);
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (a, inexact, rnd_mode);
}

int
mpfr_mul_d (mpfr_ptr a, mpfr_srcptr b, double c, mpfr_rnd_t rnd_mode)
{
  int inexact;
  mpfr_t d;
  mp_limb_t tmp_man[MPFR_LIMBS_PER_DOUBLE];
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_LOG_FUNC
    (("b[%Pu]=%.*Rg c=%.20g rnd=%d",
      mpfr_get_prec(b), mpfr_log_prec, b, c, rnd_mode),
     ("a[%Pu]=%.*Rg inexact=%d",
      mpfr_get_prec (a), mpfr_log_prec, a, inexact));

  MPFR_SAVE_EXPO_MARK (expo);

  MPFR_TMP_INIT1(tmp_man, d, IEEE_DBL_MANT_DIG);
  inexact = mpfr_set_d (d, c, rnd_mode);
  MPFR_ASSERTD (inexact == 0);

  MPFR_CLEAR_FLAGS ();
  inexact = mpfr_mul (a, b, d, rnd_mode);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (a, inexact, rnd_mode);
}

int
mpfr_sub_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode)
{
  if (MPFR_LIKELY (u != 0))  /* if u=0, do nothing */
    {
      mpfr_t uu;
      mp_limb_t up[1];
      unsigned long cnt;
      int inex;

      MPFR_SAVE_EXPO_DECL (expo);

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

      /* Optimization note: Exponent save/restore operations may be
         removed if mpfr_sub works even when uu is out-of-range. */
      MPFR_SAVE_EXPO_MARK (expo);
      MPFR_SET_EXP (uu, BITS_PER_MP_LIMB - cnt);
      inex = mpfr_sub (y, x, uu, rnd_mode);
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_check_range (y, inex, rnd_mode);
    }
  else
    return mpfr_set (y, x, rnd_mode);
}

int
mpfr_ui_pow (mpfr_ptr y, unsigned long int n, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t t;
  int inexact;
  mp_limb_t tmp_mant[(sizeof (n) - 1) / sizeof (mp_limb_t) + 1];
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_SAVE_EXPO_MARK (expo);
  MPFR_TMP_INIT1(tmp_mant, t, sizeof(n) * CHAR_BIT);
  inexact = mpfr_set_ui (t, n, MPFR_RNDN);
  MPFR_ASSERTD (inexact == 0);
  inexact = mpfr_pow (y, t, x, rnd_mode);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (y, inexact, rnd_mode);
}

int
mpfr_get_z (mpz_ptr z, mpfr_srcptr f, mpfr_rnd_t rnd)
{
  int inex;
  mpfr_t r;
  mpfr_exp_t exp;
  MPFR_SAVE_EXPO_DECL (expo);

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (f)))
    {
      if (MPFR_UNLIKELY (MPFR_NOTZERO (f)))
        MPFR_SET_ERANGEFLAG ();
      mpz_set_ui (z, 0);
      /* The ternary value is 0 even for infinity. Giving the rounding
         direction in this case would not make much sense anyway, and
         the direction would not necessarily match rnd. */
      return 0;
    }

  MPFR_SAVE_EXPO_MARK (expo);

  exp = MPFR_GET_EXP (f);
  /* if exp <= 0, then |f|<1, thus |o(f)|<=1 */
  MPFR_ASSERTN (exp < 0 || exp <= MPFR_PREC_MAX);
  mpfr_init2 (r, (exp < (mpfr_exp_t) MPFR_PREC_MIN ?
                  MPFR_PREC_MIN : (mpfr_prec_t) exp));
  inex = mpfr_rint (r, f, rnd);
  MPFR_ASSERTN (inex != 1 && inex != -1); /* integral part of f is
                                             representable in r */
  MPFR_ASSERTN (MPFR_IS_FP (r));

  /* The flags from mpfr_rint are the wanted ones. In particular,
     it sets the inexact flag when necessary. */
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);

  exp = mpfr_get_z_2exp (z, r);
  if (exp >= 0)
    mpz_mul_2exp (z, z, exp);
  else
    mpz_fdiv_q_2exp (z, z, -exp);
  mpfr_clear (r);

  MPFR_SAVE_EXPO_FREE (expo);

  return inex;
}

static int
pi_div_2ui (mpfr_ptr dest, int i, int neg, mpfr_rnd_t rnd_mode)
{
  int inexact;
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_SAVE_EXPO_MARK (expo);
  if (neg) /* -PI/2^i */
    {
      inexact = - mpfr_const_pi (dest, MPFR_INVERT_RND (rnd_mode));
      MPFR_CHANGE_SIGN (dest);
    }
  else /* PI/2^i */
    {
      inexact = mpfr_const_pi (dest, rnd_mode);
    }
  mpfr_div_2ui (dest, dest, i, rnd_mode);  /* exact */
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (dest, inexact, rnd_mode);
}

int
mpfr_cmp_d (mpfr_srcptr b, double d)
{
  mpfr_t tmp;
  int res;
  mp_limb_t tmp_man[MPFR_LIMBS_PER_DOUBLE];
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_SAVE_EXPO_MARK (expo);

  MPFR_TMP_INIT1(tmp_man, tmp, IEEE_DBL_MANT_DIG);
  res = mpfr_set_d (tmp, d, MPFR_RNDN);
  MPFR_ASSERTD (res == 0);

  MPFR_CLEAR_FLAGS ();
  res = mpfr_cmp (b, tmp);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);

  MPFR_SAVE_EXPO_FREE (expo);
  return res;
}

int
mpfr_rint_trunc (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
{
  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
    return mpfr_set (r, u, rnd_mode);
  else
    {
      mpfr_t tmp;
      int inex;
      MPFR_SAVE_EXPO_DECL (expo);

      MPFR_SAVE_EXPO_MARK (expo);
      mpfr_init2 (tmp, MPFR_PREC (u));
      /* trunc(u) is always representable in tmp */
      mpfr_trunc (tmp, u);
      inex = mpfr_set (r, tmp, rnd_mode);
      mpfr_clear (tmp);
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_check_range (r, inex, rnd_mode);
    }
}

int
mpfr_div_d (mpfr_ptr a, mpfr_srcptr b, double c, mp_rnd_t rnd_mode)
{
  int inexact;
  mpfr_t d;
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_LOG_FUNC (("b[%#R]=%R c%.20g rnd=%d", b, b, c, rnd_mode),
                 ("a[%#R]=%R", a, a));

  MPFR_SAVE_EXPO_MARK (expo);

  mpfr_init2 (d, IEEE_DBL_MANT_DIG);
  inexact = mpfr_set_d (d, c, rnd_mode);
  MPFR_ASSERTN (inexact == 0);

  mpfr_clear_flags ();
  inexact = mpfr_div (a, b, d, rnd_mode);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);

  mpfr_clear(d);
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (a, inexact, rnd_mode);
}

int
mpfr_rint_floor (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
{
  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
    return mpfr_set (r, u, rnd_mode);
  else
    {
      mpfr_t tmp;
      int inex;
      MPFR_SAVE_EXPO_DECL (expo);
      MPFR_BLOCK_DECL (flags);

      MPFR_SAVE_EXPO_MARK (expo);
      mpfr_init2 (tmp, MPFR_PREC (u));
      /* floor(u) is representable in tmp unless an overflow occurs */
      MPFR_BLOCK (flags, mpfr_floor (tmp, u));
      inex = (MPFR_OVERFLOW (flags)
              ? mpfr_overflow (r, rnd_mode, MPFR_SIGN_NEG)
              : mpfr_set (r, tmp, rnd_mode));
      mpfr_clear (tmp);
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_check_range (r, inex, rnd_mode);
    }
}

int
mpfr_atanh (mpfr_ptr y, mpfr_srcptr xt , mpfr_rnd_t rnd_mode)
{
  int inexact;
  mpfr_t x, t, te;
  mpfr_prec_t Nx, Ny, Nt;
  mpfr_exp_t err;
  MPFR_ZIV_DECL (loop);
  MPFR_SAVE_EXPO_DECL (expo);

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

  /* Special cases */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
    {
      /* atanh(NaN) = NaN, and atanh(+/-Inf) = NaN since tanh gives a result
         between -1 and 1 */
      if (MPFR_IS_NAN (xt) || MPFR_IS_INF (xt))
        {
          MPFR_SET_NAN (y);
          MPFR_RET_NAN;
        }
      else /* necessarily xt is 0 */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (xt));
          MPFR_SET_ZERO (y);   /* atanh(0) = 0 */
          MPFR_SET_SAME_SIGN (y,xt);
          MPFR_RET (0);
        }
    }

  /* atanh (x) = NaN as soon as |x| > 1, and arctanh(+/-1) = +/-Inf */
  if (MPFR_UNLIKELY (MPFR_GET_EXP (xt) > 0))
    {
      if (MPFR_GET_EXP (xt) == 1 && mpfr_powerof2_raw (xt))
        {
          MPFR_SET_INF (y);
          MPFR_SET_SAME_SIGN (y, xt);
          mpfr_set_divby0 ();
          MPFR_RET (0);
        }
      MPFR_SET_NAN (y);
      MPFR_RET_NAN;
    }

  /* atanh(x) = x + x^3/3 + ... so the error is < 2^(3*EXP(x)-1) */
  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2 * MPFR_GET_EXP (xt), 1, 1,
                                    rnd_mode, {});

  MPFR_SAVE_EXPO_MARK (expo);

  /* Compute initial precision */
  Nx = MPFR_PREC (xt);
  MPFR_TMP_INIT_ABS (x, xt);
  Ny = MPFR_PREC (y);
  Nt = MAX (Nx, Ny);
  /* the optimal number of bits : see algorithms.ps */
  Nt = Nt + MPFR_INT_CEIL_LOG2 (Nt) + 4;

  /* initialise of intermediary variable */
  mpfr_init2 (t, Nt);
  mpfr_init2 (te, Nt);

  /* First computation of cosh */
  MPFR_ZIV_INIT (loop, Nt);
  for (;;)
    {
      /* compute atanh */
      mpfr_ui_sub (te, 1, x, MPFR_RNDU);   /* (1-xt)*/
      mpfr_add_ui (t,  x, 1, MPFR_RNDD);   /* (xt+1)*/
      mpfr_div (t, t, te, MPFR_RNDN);      /* (1+xt)/(1-xt)*/
      mpfr_log (t, t, MPFR_RNDN);          /* ln((1+xt)/(1-xt))*/
      mpfr_div_2ui (t, t, 1, MPFR_RNDN);   /* (1/2)*ln((1+xt)/(1-xt))*/

      /* error estimate: see algorithms.tex */
      /* FIXME: this does not correspond to the value in algorithms.tex!!! */
      /* err=Nt-__gmpfr_ceil_log2(1+5*pow(2,1-MPFR_EXP(t)));*/
      err = Nt - (MAX (4 - MPFR_GET_EXP (t), 0) + 1);

      if (MPFR_LIKELY (MPFR_IS_ZERO (t)
                       || MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
        break;

      /* reactualisation of the precision */
      MPFR_ZIV_NEXT (loop, Nt);
      mpfr_set_prec (t, Nt);
      mpfr_set_prec (te, Nt);
    }
  MPFR_ZIV_FREE (loop);

  inexact = mpfr_set4 (y, t, rnd_mode, MPFR_SIGN (xt));

  mpfr_clear(t);
  mpfr_clear(te);

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (y, inexact, rnd_mode);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    ("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
    ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y,
     inexact));

  /* special values */
  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))
        {
          MPFR_SET_INF (y);
          MPFR_SET_SAME_SIGN (y, x);
          MPFR_RET (0);
        }
      /* case 0: cbrt(+/- 0) = +/- 0 */
      else /* x is necessarily 0 */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
          MPFR_SET_ZERO (y);
          MPFR_SET_SAME_SIGN (y, x);
          MPFR_RET (0);
        }
    }

  /* General case */
  MPFR_SAVE_EXPO_MARK (expo);
  mpz_init (m);

  e = mpfr_get_z_2exp (m, x);                /* x = m * 2^e */
  if ((negative = MPFR_IS_NEG(x)))
    mpz_neg (m, m);
  r = e % 3;
  if (r < 0)
    r += 3;
  /* x = (m*2^r) * 2^(e-r) = (m*2^r) * 2^(3*q) */

  MPFR_MPZ_SIZEINBASE2 (size_m, m);
  n = MPFR_PREC (y) + (rnd_mode == MPFR_RNDN);

  /* we want 3*n-2 <= size_m + 3*sh + r <= 3*n
     i.e. 3*sh + size_m + r <= 3*n */
  sh = (3 * (mpfr_exp_t) n - (mpfr_exp_t) size_m - r) / 3;
  sh = 3 * sh + r;
  if (sh >= 0)
    {
      mpz_mul_2exp (m, m, sh);
      e = e - sh;
    }
  else if (r > 0)
    {
      mpz_mul_2exp (m, m, r);
      e = e - r;
    }

  /* invariant: x = m*2^e, with e divisible by 3 */

  /* we reuse the variable m to store the cube root, since it is not needed
     any more: we just need to know if the root is exact */
  inexact = mpz_root (m, m, 3) == 0;

  MPFR_MPZ_SIZEINBASE2 (tmp, m);
  sh = tmp - n;
  if (sh > 0) /* we have to flush to 0 the last sh bits from m */
    {
      inexact = inexact || ((mpfr_exp_t) mpz_scan1 (m, 0) < sh);
      mpz_fdiv_q_2exp (m, m, sh);
      e += 3 * sh;
    }

  if (inexact)
    {
      if (negative)
        rnd_mode = MPFR_INVERT_RND (rnd_mode);
      if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA
          || (rnd_mode == MPFR_RNDN && mpz_tstbit (m, 0)))
        inexact = 1, mpz_add_ui (m, m, 1);
      else
        inexact = -1;
    }

  /* either inexact is not zero, and the conversion is exact, i.e. inexact
     is not changed; or inexact=0, and inexact is set only when
     rnd_mode=MPFR_RNDN and bit (n+1) from m is 1 */
  inexact += mpfr_set_z (y, m, MPFR_RNDN);
  MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e / 3);

  if (negative)
    {
      MPFR_CHANGE_SIGN (y);
      inexact = -inexact;
    }

  mpz_clear (m);
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (y, inexact, rnd_mode);
}

//.........这里部分代码省略.........
            MPFR_SET_INF (y);
          else
            MPFR_SET_ZERO (y);
          MPFR_SET_POS (y);
          MPFR_RET (0);
        }
      else /* 2^0 = 1 */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO(x));
          return mpfr_set_ui (y, 1, rnd_mode);
        }
    }

  /* since the smallest representable non-zero float is 1/2*2^__gmpfr_emin,
     if x < __gmpfr_emin - 1, the result is either 1/2*2^__gmpfr_emin or 0 */
  MPFR_ASSERTN (MPFR_EMIN_MIN >= LONG_MIN + 2);
  if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emin - 1) < 0))
    {
      mpfr_rnd_t rnd2 = rnd_mode;
      /* in round to nearest mode, round to zero when x <= __gmpfr_emin-2 */
      if (rnd_mode == MPFR_RNDN &&
          mpfr_cmp_si_2exp (x, __gmpfr_emin - 2, 0) <= 0)
        rnd2 = MPFR_RNDZ;
      return mpfr_underflow (y, rnd2, 1);
    }

  MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
  if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emax) >= 0))
    return mpfr_overflow (y, rnd_mode, 1);

  /* We now know that emin - 1 <= x < emax. */

  MPFR_SAVE_EXPO_MARK (expo);

  /* 2^x = 1 + x*log(2) + O(x^2) for x near zero, and for |x| <= 1 we have
     |2^x - 1| <= x < 2^EXP(x). If x > 0 we must round away from 0 (dir=1);
     if x < 0 we must round toward 0 (dir=0). */
  MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, __gmpfr_one, - MPFR_GET_EXP (x), 0,
                                    MPFR_IS_POS (x), rnd_mode, expo, {});

  xint = mpfr_get_si (x, MPFR_RNDZ);
  mpfr_init2 (xfrac, MPFR_PREC (x));
  mpfr_sub_si (xfrac, x, xint, MPFR_RNDN); /* exact */

  if (MPFR_IS_ZERO (xfrac))
    {
      mpfr_set_ui (y, 1, MPFR_RNDN);
      inexact = 0;
    }
  else
    {
      /* Declaration of the intermediary variable */
      mpfr_t t;

      /* Declaration of the size variable */
      mpfr_prec_t Ny = MPFR_PREC(y);              /* target precision */
      mpfr_prec_t Nt;                             /* working precision */
      mpfr_exp_t err;                             /* error */
      MPFR_ZIV_DECL (loop);

      /* compute the precision of intermediary variable */
      /* the optimal number of bits : see algorithms.tex */
      Nt = Ny + 5 + MPFR_INT_CEIL_LOG2 (Ny);

      /* initialize of intermediary variable */
      mpfr_init2 (t, Nt);

      /* First computation */
      MPFR_ZIV_INIT (loop, Nt);
      for (;;)
        {
          /* compute exp(x*ln(2))*/
          mpfr_const_log2 (t, MPFR_RNDU);       /* ln(2) */
          mpfr_mul (t, xfrac, t, MPFR_RNDU);    /* xfrac * ln(2) */
          err = Nt - (MPFR_GET_EXP (t) + 2);   /* Estimate of the error */
          mpfr_exp (t, t, MPFR_RNDN);           /* exp(xfrac * ln(2)) */

          if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
            break;

          /* Actualisation of the precision */
          MPFR_ZIV_NEXT (loop, Nt);
          mpfr_set_prec (t, Nt);
        }
      MPFR_ZIV_FREE (loop);

      inexact = mpfr_set (y, t, rnd_mode);

      mpfr_clear (t);
    }

  mpfr_clear (xfrac);
  MPFR_CLEAR_FLAGS ();
  mpfr_mul_2si (y, y, xint, MPFR_RNDN); /* exact or overflow */
  /* Note: We can have an overflow only when t was rounded up to 2. */
  MPFR_ASSERTD (MPFR_IS_PURE_FP (y) || inexact > 0);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (y, inexact, rnd_mode);
}

int
mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t xp;
  int compared, inexact;
  mpfr_prec_t prec;
  mpfr_exp_t xp_exp;
  MPFR_SAVE_EXPO_DECL (expo);
  MPFR_ZIV_DECL (loop);

  MPFR_LOG_FUNC (
    ("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
    ("asin[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (asin), mpfr_log_prec, asin,
     inexact));

  /* Special cases */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
        {
          MPFR_SET_NAN (asin);
          MPFR_RET_NAN;
        }
      else /* x = 0 */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
          MPFR_SET_ZERO (asin);
          MPFR_SET_SAME_SIGN (asin, x);
          MPFR_RET (0); /* exact result */
        }
    }

  /* asin(x) = x + x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (asin, x, -2 * MPFR_GET_EXP (x), 2, 1,
                                    rnd_mode, {});

  /* Set x_p=|x| (x is a normal number) */
  mpfr_init2 (xp, MPFR_PREC (x));
  inexact = mpfr_abs (xp, x, MPFR_RNDN);
  MPFR_ASSERTD (inexact == 0);

  compared = mpfr_cmp_ui (xp, 1);

  MPFR_SAVE_EXPO_MARK (expo);

  if (MPFR_UNLIKELY (compared >= 0))
    {
      mpfr_clear (xp);
      if (compared > 0)                  /* asin(x) = NaN for |x| > 1 */
        {
          MPFR_SAVE_EXPO_FREE (expo);
          MPFR_SET_NAN (asin);
          MPFR_RET_NAN;
        }
      else                              /* x = 1 or x = -1 */
        {
          if (MPFR_IS_POS (x)) /* asin(+1) = Pi/2 */
            inexact = mpfr_const_pi (asin, rnd_mode);
          else /* asin(-1) = -Pi/2 */
            {
              inexact = -mpfr_const_pi (asin, MPFR_INVERT_RND(rnd_mode));
              MPFR_CHANGE_SIGN (asin);
            }
          mpfr_div_2ui (asin, asin, 1, rnd_mode);
        }
    }
  else
    {
      /* Compute exponent of 1 - ABS(x) */
      mpfr_ui_sub (xp, 1, xp, MPFR_RNDD);
      MPFR_ASSERTD (MPFR_GET_EXP (xp) <= 0);
      MPFR_ASSERTD (MPFR_GET_EXP (x) <= 0);
      xp_exp = 2 - MPFR_GET_EXP (xp);

      /* Set up initial prec */
      prec = MPFR_PREC (asin) + 10 + xp_exp;

      /* use asin(x) = atan(x/sqrt(1-x^2)) */
      MPFR_ZIV_INIT (loop, prec);
      for (;;)
        {
          mpfr_set_prec (xp, prec);
          mpfr_sqr (xp, x, MPFR_RNDN);
          mpfr_ui_sub (xp, 1, xp, MPFR_RNDN);
          mpfr_sqrt (xp, xp, MPFR_RNDN);
          mpfr_div (xp, x, xp, MPFR_RNDN);
          mpfr_atan (xp, xp, MPFR_RNDN);
          if (MPFR_LIKELY (MPFR_CAN_ROUND (xp, prec - xp_exp,
                                           MPFR_PREC (asin), rnd_mode)))
            break;
          MPFR_ZIV_NEXT (loop, prec);
        }
      MPFR_ZIV_FREE (loop);
      inexact = mpfr_set (asin, xp, rnd_mode);

      mpfr_clear (xp);
    }

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (asin, inexact, rnd_mode);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    mp_prec_t prec;
    mpfr_t t, slong;
    mp_rnd_t rnd;
    int inexact, compare;
    unsigned int n;

    mp_prec_t p0=2, p1=200;
    unsigned int N=200;

    mpfr_init (t);
    mpfr_init (slong);

    /* generic test */
    for (prec = p0; prec <= p1; prec++)
    {
      mpfr_set_prec (x, prec);
      mpfr_set_prec (y, prec);
      mpfr_set_prec (z, prec);
      mpfr_set_prec (s, prec);
      mpfr_set_prec (t, prec);

      for (n=0; n<N; n++)
        {
          mpfr_urandomb (x, RANDS);
          mpfr_urandomb (y, RANDS);
          mpfr_urandomb (z, RANDS);

          if (randlimb () % 2)
            mpfr_neg (x, x, GMP_RNDN);
          if (randlimb () % 2)
            mpfr_neg (y, y, GMP_RNDN);
          if (randlimb () % 2)
            mpfr_neg (z, z, GMP_RNDN);

          rnd = RND_RAND ();
          mpfr_set_prec (slong, 2 * prec);
          if (mpfr_mul (slong, x, y, rnd))
            {
              printf ("x*y should be exact\n");
              exit (1);
            }
          compare = mpfr_sub (t, slong, z, rnd);
          inexact = mpfr_fms (s, x, y, z, rnd);
          if (mpfr_cmp (s, t))
            {
              printf ("results differ for x=");
              mpfr_out_str (stdout, 2, prec, x, GMP_RNDN);
              printf ("  y=");
              mpfr_out_str (stdout, 2, prec, y, GMP_RNDN);
              printf ("  z=");
              mpfr_out_str (stdout, 2, prec, z, GMP_RNDN);
              printf (" prec=%u rnd_mode=%s\n", (unsigned int) prec,
                      mpfr_print_rnd_mode (rnd));
              printf ("got      ");
              mpfr_out_str (stdout, 2, prec, s, GMP_RNDN);
              puts ("");
              printf ("expected ");
              mpfr_out_str (stdout, 2, prec, t, GMP_RNDN);
              puts ("");
              printf ("approx  ");
              mpfr_print_binary (slong);
              puts ("");
              exit (1);
            }
          if (((inexact == 0) && (compare != 0)) ||
              ((inexact < 0) && (compare >= 0)) ||
              ((inexact > 0) && (compare <= 0)))
            {
              printf ("Wrong inexact flag for rnd=%s: expected %d, got %d\n",
                      mpfr_print_rnd_mode (rnd), compare, inexact);
              printf (" x="); mpfr_out_str (stdout, 2, 0, x, GMP_RNDN);
              printf (" y="); mpfr_out_str (stdout, 2, 0, y, GMP_RNDN);
              printf (" z="); mpfr_out_str (stdout, 2, 0, z, GMP_RNDN);
              printf (" s="); mpfr_out_str (stdout, 2, 0, s, GMP_RNDN);
              printf ("\n");
              exit (1);
            }
        }
    }
  mpfr_clear (t);
  mpfr_clear (slong);

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

  test_exact ();

  MPFR_SAVE_EXPO_MARK (expo);
  test_overflow1 ();
  test_overflow2 ();
  test_underflow1 ();
  test_underflow2 ();
  MPFR_SAVE_EXPO_FREE (expo);

  tests_end_mpfr ();
  return 0;
}