int
mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
{
  int inexact;
  long xint;
  mpfr_t xfrac;
  MPFR_SAVE_EXPO_DECL (expo);

  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 /* 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))
    {
      mp_rnd_t rnd2 = rnd_mode;
      /* in round to nearest mode, round to zero when x <= __gmpfr_emin-2 */
      if (rnd_mode == GMP_RNDN &&
          mpfr_cmp_si_2exp (x, __gmpfr_emin - 2, 0) <= 0)
        rnd2 = GMP_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_SIGN(x) > 0, rnd_mode, expo, {});

  xint = mpfr_get_si (x, GMP_RNDZ);
  mpfr_init2 (xfrac, MPFR_PREC (x));
  mpfr_sub_si (xfrac, x, xint, GMP_RNDN); /* exact */

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

      /* Declaration of the size variable */
      mp_prec_t Ny = MPFR_PREC(y);              /* target precision */
      mp_prec_t Nt;                             /* working precision */
      mp_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);

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

      /* First computation */
      MPFR_ZIV_INIT (loop, Nt);
      for (;;)
        {
          /* compute exp(x*ln(2))*/
          mpfr_const_log2 (t, GMP_RNDU);       /* ln(2) */
          mpfr_mul (t, xfrac, t, GMP_RNDU);    /* xfrac * ln(2) */
          err = Nt - (MPFR_GET_EXP (t) + 2);   /* Estimate of the error */
          mpfr_exp (t, t, GMP_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);
//.........这里部分代码省略.........

static int
mpfr_mul3 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
    /* Old implementation */
    int sign_product, cc, inexact;
    mpfr_exp_t ax;
    mp_limb_t *tmp;
    mp_limb_t b1;
    mpfr_prec_t bq, cq;
    mp_size_t bn, cn, tn, k;
    MPFR_TMP_DECL(marker);

    /* deal with special cases */
    if (MPFR_ARE_SINGULAR(b,c))
    {
        if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c))
        {
            MPFR_SET_NAN(a);
            MPFR_RET_NAN;
        }
        sign_product = MPFR_MULT_SIGN( MPFR_SIGN(b) , MPFR_SIGN(c) );
        if (MPFR_IS_INF(b))
        {
            if (MPFR_IS_INF(c) || MPFR_NOTZERO(c))
            {
                MPFR_SET_SIGN(a,sign_product);
                MPFR_SET_INF(a);
                MPFR_RET(0); /* exact */
            }
            else
            {
                MPFR_SET_NAN(a);
                MPFR_RET_NAN;
            }
        }
        else if (MPFR_IS_INF(c))
        {
            if (MPFR_NOTZERO(b))
            {
                MPFR_SET_SIGN(a, sign_product);
                MPFR_SET_INF(a);
                MPFR_RET(0); /* exact */
            }
            else
            {
                MPFR_SET_NAN(a);
                MPFR_RET_NAN;
            }
        }
        else
        {
            MPFR_ASSERTD(MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c));
            MPFR_SET_SIGN(a, sign_product);
            MPFR_SET_ZERO(a);
            MPFR_RET(0); /* 0 * 0 is exact */
        }
    }
    sign_product = MPFR_MULT_SIGN( MPFR_SIGN(b) , MPFR_SIGN(c) );

    ax = MPFR_GET_EXP (b) + MPFR_GET_EXP (c);

    bq = MPFR_PREC(b);
    cq = MPFR_PREC(c);

    MPFR_ASSERTD(bq+cq > bq); /* PREC_MAX is /2 so no integer overflow */

    bn = (bq+GMP_NUMB_BITS-1)/GMP_NUMB_BITS; /* number of limbs of b */
    cn = (cq+GMP_NUMB_BITS-1)/GMP_NUMB_BITS; /* number of limbs of c */
    k = bn + cn; /* effective nb of limbs used by b*c (= tn or tn+1) below */
    tn = (bq + cq + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
    /* <= k, thus no int overflow */
    MPFR_ASSERTD(tn <= k);

    /* Check for no size_t overflow*/
    MPFR_ASSERTD((size_t) k <= ((size_t) -1) / BYTES_PER_MP_LIMB);
    MPFR_TMP_MARK(marker);
    tmp = (mp_limb_t *) MPFR_TMP_ALLOC((size_t) k * BYTES_PER_MP_LIMB);

    /* multiplies two mantissa in temporary allocated space */
    b1 = (MPFR_LIKELY(bn >= cn)) ?
         mpn_mul (tmp, MPFR_MANT(b), bn, MPFR_MANT(c), cn)
         : mpn_mul (tmp, MPFR_MANT(c), cn, MPFR_MANT(b), bn);

    /* now tmp[0]..tmp[k-1] contains the product of both mantissa,
       with tmp[k-1]>=2^(GMP_NUMB_BITS-2) */
    b1 >>= GMP_NUMB_BITS - 1; /* msb from the product */

    /* if the mantissas of b and c are uniformly distributed in ]1/2, 1],
       then their product is in ]1/4, 1/2] with probability 2*ln(2)-1 ~ 0.386
       and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */
    tmp += k - tn;
    if (MPFR_UNLIKELY(b1 == 0))
        mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */
    cc = mpfr_round_raw (MPFR_MANT (a), tmp, bq + cq,
                         MPFR_IS_NEG_SIGN(sign_product),
                         MPFR_PREC (a), rnd_mode, &inexact);

    /* cc = 1 ==> result is a power of two */
    if (MPFR_UNLIKELY(cc))
        MPFR_MANT(a)[MPFR_LIMB_SIZE(a)-1] = MPFR_LIMB_HIGHBIT;
//.........这里部分代码省略.........

int
mpfr_sqr (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
{
  int cc, inexact;
  mpfr_exp_t ax;
  mp_limb_t *tmp;
  mp_limb_t b1;
  mpfr_prec_t bq;
  mp_size_t bn, tn;
  MPFR_TMP_DECL(marker);

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

  /* deal with special cases */
  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(b)))
    {
      if (MPFR_IS_NAN(b))
        {
          MPFR_SET_NAN(a);
          MPFR_RET_NAN;
        }
      MPFR_SET_POS (a);
      if (MPFR_IS_INF(b))
        MPFR_SET_INF(a);
      else
        ( MPFR_ASSERTD(MPFR_IS_ZERO(b)), MPFR_SET_ZERO(a) );
      MPFR_RET(0);
    }
  ax = 2 * MPFR_GET_EXP (b);
  bq = MPFR_PREC(b);

  MPFR_ASSERTD (2 * bq > bq); /* PREC_MAX is /2 so no integer overflow */

  bn = MPFR_LIMB_SIZE(b); /* number of limbs of b */
  tn = 1 + (2 * bq - 1) / GMP_NUMB_BITS; /* number of limbs of square,
                                               2*bn or 2*bn-1 */

  MPFR_TMP_MARK(marker);
  tmp = (mp_limb_t *) MPFR_TMP_ALLOC((size_t) 2 * bn * BYTES_PER_MP_LIMB);

  /* Multiplies the mantissa in temporary allocated space */
  mpn_sqr_n (tmp, MPFR_MANT(b), bn);
  b1 = tmp[2 * bn - 1];

  /* now tmp[0]..tmp[2*bn-1] contains the product of both mantissa,
     with tmp[2*bn-1]>=2^(GMP_NUMB_BITS-2) */
  b1 >>= GMP_NUMB_BITS - 1; /* msb from the product */

  /* if the mantissas of b and c are uniformly distributed in ]1/2, 1],
     then their product is in ]1/4, 1/2] with probability 2*ln(2)-1 ~ 0.386
     and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */
  tmp += 2 * bn - tn; /* +0 or +1 */
  if (MPFR_UNLIKELY(b1 == 0))
    mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */

  cc = mpfr_round_raw (MPFR_MANT (a), tmp, 2 * bq, 0,
                       MPFR_PREC (a), rnd_mode, &inexact);
  /* cc = 1 ==> result is a power of two */
  if (MPFR_UNLIKELY(cc))
    MPFR_MANT(a)[MPFR_LIMB_SIZE(a)-1] = MPFR_LIMB_HIGHBIT;

  MPFR_TMP_FREE(marker);
  {
    mpfr_exp_t ax2 = ax + (mpfr_exp_t) (b1 - 1 + cc);
    if (MPFR_UNLIKELY( ax2 > __gmpfr_emax))
      return mpfr_overflow (a, rnd_mode, MPFR_SIGN_POS);
    if (MPFR_UNLIKELY( ax2 < __gmpfr_emin))
      {
        /* In the rounding to the nearest mode, if the exponent of the exact
           result (i.e. before rounding, i.e. without taking cc into account)
           is < __gmpfr_emin - 1 or the exact result is a power of 2 (i.e. if
           both arguments are powers of 2), then round to zero. */
        if (rnd_mode == MPFR_RNDN &&
            (ax + (mpfr_exp_t) b1 < __gmpfr_emin || mpfr_powerof2_raw (b)))
          rnd_mode = MPFR_RNDZ;
        return mpfr_underflow (a, rnd_mode, MPFR_SIGN_POS);
      }
    MPFR_SET_EXP (a, ax2);
    MPFR_SET_POS (a);
  }
  MPFR_RET (inexact);
}

int
mpfr_sqrt (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
{
  mp_size_t rsize; /* number of limbs of r (plus 1 if exact limb multiple) */
  mp_size_t rrsize;
  mp_size_t usize; /* number of limbs of u */
  mp_size_t tsize; /* number of limbs of the sqrtrem remainder */
  mp_size_t k;
  mp_size_t l;
  mpfr_limb_ptr rp, rp0;
  mpfr_limb_ptr up;
  mpfr_limb_ptr sp;
  mp_limb_t sticky0; /* truncated part of input */
  mp_limb_t sticky1; /* truncated part of rp[0] */
  mp_limb_t sticky;
  int odd_exp;
  int sh; /* number of extra bits in rp[0] */
  int inexact; /* return ternary flag */
  mpfr_exp_t expr;
  MPFR_TMP_DECL(marker);

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

  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(u)))
    {
      if (MPFR_IS_NAN(u))
        {
          MPFR_SET_NAN(r);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_ZERO(u))
        {
          /* 0+ or 0- */
          MPFR_SET_SAME_SIGN(r, u);
          MPFR_SET_ZERO(r);
          MPFR_RET(0); /* zero is exact */
        }
      else
        {
          MPFR_ASSERTD(MPFR_IS_INF(u));
          /* sqrt(-Inf) = NAN */
          if (MPFR_IS_NEG(u))
            {
              MPFR_SET_NAN(r);
              MPFR_RET_NAN;
            }
          MPFR_SET_POS(r);
          MPFR_SET_INF(r);
          MPFR_RET(0);
        }
    }
  if (MPFR_UNLIKELY(MPFR_IS_NEG(u)))
    {
      MPFR_SET_NAN(r);
      MPFR_RET_NAN;
    }
  MPFR_SET_POS(r);

  MPFR_TMP_MARK (marker);
  MPFR_UNSIGNED_MINUS_MODULO(sh,MPFR_PREC(r));
  if (sh == 0 && rnd_mode == MPFR_RNDN)
    sh = GMP_NUMB_BITS; /* ugly case */
  rsize = MPFR_LIMB_SIZE(r) + (sh == GMP_NUMB_BITS);
  /* rsize is the number of limbs of r + 1 if exact limb multiple and rounding
     to nearest, this is the number of wanted limbs for the square root */
  rrsize = rsize + rsize;
  usize = MPFR_LIMB_SIZE(u); /* number of limbs of u */
  rp0 = MPFR_MANT(r);
  rp = (sh < GMP_NUMB_BITS) ? rp0 : MPFR_TMP_LIMBS_ALLOC (rsize);
  up = MPFR_MANT(u);
  sticky0 = MPFR_LIMB_ZERO; /* truncated part of input */
  sticky1 = MPFR_LIMB_ZERO; /* truncated part of rp[0] */
  odd_exp = (unsigned int) MPFR_GET_EXP (u) & 1;
  inexact = -1; /* return ternary flag */

  sp = MPFR_TMP_LIMBS_ALLOC (rrsize);

  /* copy the most significant limbs of u to {sp, rrsize} */
  if (MPFR_LIKELY(usize <= rrsize)) /* in case r and u have the same precision,
                                       we have indeed rrsize = 2 * usize */
    {
      k = rrsize - usize;
      if (MPFR_LIKELY(k))
        MPN_ZERO (sp, k);
      if (odd_exp)
        {
          if (MPFR_LIKELY(k))
            sp[k - 1] = mpn_rshift (sp + k, up, usize, 1);
          else
            sticky0 = mpn_rshift (sp, up, usize, 1);
        }
      else
        MPN_COPY (sp + rrsize - usize, up, usize);
    }
  else /* usize > rrsize: truncate the input */
    {
      k = usize - rrsize;
//.........这里部分代码省略.........

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);
//.........这里部分代码省略.........

int
mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd)
{
  MPFR_LOG_FUNC
    (("x[%Pu]=%.*Rg n=%ld rnd=%d",
      mpfr_get_prec (x), mpfr_log_prec, x, n, rnd),
     ("y[%Pu]=%.*Rg", mpfr_get_prec (y), mpfr_log_prec, y));

  if (n >= 0)
    return mpfr_pow_ui (y, x, n, rnd);
  else
    {
      if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
        {
          if (MPFR_IS_NAN (x))
            {
              MPFR_SET_NAN (y);
              MPFR_RET_NAN;
            }
          else
            {
              int positive = MPFR_IS_POS (x) || ((unsigned long) n & 1) == 0;
              if (MPFR_IS_INF (x))
                MPFR_SET_ZERO (y);
              else /* x is zero */
                {
                  MPFR_ASSERTD (MPFR_IS_ZERO (x));
                  MPFR_SET_INF (y);
                  mpfr_set_divby0 ();
                }
              if (positive)
                MPFR_SET_POS (y);
              else
                MPFR_SET_NEG (y);
              MPFR_RET (0);
            }
        }

      /* detect exact powers: x^(-n) is exact iff x is a power of 2 */
      if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0)
        {
          mpfr_exp_t expx = MPFR_EXP (x) - 1, expy;
          MPFR_ASSERTD (n < 0);
          /* Warning: n * expx may overflow!
           *
           * Some systems (apparently alpha-freebsd) abort with
           * LONG_MIN / 1, and LONG_MIN / -1 is undefined.
           * http://www.freebsd.org/cgi/query-pr.cgi?pr=72024
           *
           * Proof of the overflow checking. The expressions below are
           * assumed to be on the rational numbers, but the word "overflow"
           * still has its own meaning in the C context. / still denotes
           * the integer (truncated) division, and // denotes the exact
           * division.
           * - First, (__gmpfr_emin - 1) / n and (__gmpfr_emax - 1) / n
           *   cannot overflow due to the constraints on the exponents of
           *   MPFR numbers.
           * - If n = -1, then n * expx = - expx, which is representable
           *   because of the constraints on the exponents of MPFR numbers.
           * - If expx = 0, then n * expx = 0, which is representable.
           * - If n < -1 and expx > 0:
           *   + If expx > (__gmpfr_emin - 1) / n, then
           *           expx >= (__gmpfr_emin - 1) / n + 1
           *                > (__gmpfr_emin - 1) // n,
           *     and
           *           n * expx < __gmpfr_emin - 1,
           *     i.e.
           *           n * expx <= __gmpfr_emin - 2.
           *     This corresponds to an underflow, with a null result in
           *     the rounding-to-nearest mode.
           *   + If expx <= (__gmpfr_emin - 1) / n, then n * expx cannot
           *     overflow since 0 < expx <= (__gmpfr_emin - 1) / n and
           *           0 > n * expx >= n * ((__gmpfr_emin - 1) / n)
           *                        >= __gmpfr_emin - 1.
           * - If n < -1 and expx < 0:
           *   + If expx < (__gmpfr_emax - 1) / n, then
           *           expx <= (__gmpfr_emax - 1) / n - 1
           *                < (__gmpfr_emax - 1) // n,
           *     and
           *           n * expx > __gmpfr_emax - 1,
           *     i.e.
           *           n * expx >= __gmpfr_emax.
           *     This corresponds to an overflow (2^(n * expx) has an
           *     exponent > __gmpfr_emax).
           *   + If expx >= (__gmpfr_emax - 1) / n, then n * expx cannot
           *     overflow since 0 > expx >= (__gmpfr_emax - 1) / n and
           *           0 < n * expx <= n * ((__gmpfr_emax - 1) / n)
           *                        <= __gmpfr_emax - 1.
           * Note: one could use expx bounds based on MPFR_EXP_MIN and
           * MPFR_EXP_MAX instead of __gmpfr_emin and __gmpfr_emax. The
           * current bounds do not lead to noticeably slower code and
           * allow us to avoid a bug in Sun's compiler for Solaris/x86
           * (when optimizations are enabled); known affected versions:
           *   cc: Sun C 5.8 2005/10/13
           *   cc: Sun C 5.8 Patch 121016-02 2006/03/31
           *   cc: Sun C 5.8 Patch 121016-04 2006/10/18
           */
          expy =
            n != -1 && expx > 0 && expx > (__gmpfr_emin - 1) / n ?
            MPFR_EMIN_MIN - 2 /* Underflow */ :
//.........这里部分代码省略.........

int
mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode)
{
  int inexact;
  MPFR_SAVE_EXPO_DECL (expo);

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
    {
      /* If a is NaN, the result is NaN */
      if (MPFR_IS_NAN (a))
        {
          MPFR_SET_NAN (r);
          MPFR_RET_NAN;
        }
      /* check for infinity before zero */
      else if (MPFR_IS_INF (a))
        {
          if (MPFR_IS_NEG (a))
            /* log(-Inf) = NaN */
            {
              MPFR_SET_NAN (r);
              MPFR_RET_NAN;
            }
          else /* log(+Inf) = +Inf */
            {
              MPFR_SET_INF (r);
              MPFR_SET_POS (r);
              MPFR_RET (0);
            }
        }
      else /* a is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (a));
          MPFR_SET_INF (r);
          MPFR_SET_NEG (r);
          MPFR_RET (0); /* log2(0) is an exact -infinity */
        }
    }

  /* If a is negative, the result is NaN */
  if (MPFR_UNLIKELY (MPFR_IS_NEG (a)))
    {
      MPFR_SET_NAN (r);
      MPFR_RET_NAN;
    }

  /* If a is 1, the result is 0 */
  if (MPFR_UNLIKELY (mpfr_cmp_ui (a, 1) == 0))
    {
      MPFR_SET_ZERO (r);
      MPFR_SET_POS (r);
      MPFR_RET (0); /* only "normal" case where the result is exact */
    }

  /* If a is 2^N, log2(a) is exact*/
  if (MPFR_UNLIKELY (mpfr_cmp_ui_2exp (a, 1, MPFR_GET_EXP (a) - 1) == 0))
    return mpfr_set_si(r, MPFR_GET_EXP (a) - 1, rnd_mode);

  MPFR_SAVE_EXPO_MARK (expo);

  /* General case */
  {
    /* Declaration of the intermediary variable */
    mpfr_t t, tt;
    /* Declaration of the size variable */
    mp_prec_t Ny = MPFR_PREC(r);              /* target precision */
    mp_prec_t Nt;                             /* working precision */
    mp_exp_t err;                             /* error */
    MPFR_ZIV_DECL (loop);

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

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

    /* First computation of log2 */
    MPFR_ZIV_INIT (loop, Nt);
    for (;;)
      {
        /* compute log2 */
        mpfr_const_log2(t,GMP_RNDD); /* log(2) */
        mpfr_log(tt,a,GMP_RNDN);     /* log(a) */
        mpfr_div(t,tt,t,GMP_RNDN); /* log(a)/log(2) */

        /* estimation of the error */
        err = Nt-3;
        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_set_prec (tt, Nt);
      }
    MPFR_ZIV_FREE (loop);

    inexact = mpfr_set (r, t, rnd_mode);
//.........这里部分代码省略.........

int
mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
{
  MPFR_LOG_FUNC (("b[%#R]=%R c[%#R]=%R rnd=%d", b, b, c, c, rnd_mode),
                 ("a[%#R]=%R", a, a));

  if (MPFR_ARE_SINGULAR (b,c))
    {
      if (MPFR_IS_NAN (b) || MPFR_IS_NAN (c))
        {
          MPFR_SET_NAN (a);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF (b))
        {
          if (!MPFR_IS_INF (c) || MPFR_SIGN (b) != MPFR_SIGN(c))
            {
              MPFR_SET_INF (a);
              MPFR_SET_SAME_SIGN (a, b);
              MPFR_RET (0); /* exact */
            }
          else
            {
              MPFR_SET_NAN (a); /* Inf - Inf */
              MPFR_RET_NAN;
            }
        }
      else if (MPFR_IS_INF (c))
        {
          MPFR_SET_INF (a);
          MPFR_SET_OPPOSITE_SIGN (a, c);
          MPFR_RET (0); /* exact */
        }
      else if (MPFR_IS_ZERO (b))
        {
          if (MPFR_IS_ZERO (c))
            {
              int sign = rnd_mode != GMP_RNDD
                ? ((MPFR_IS_NEG(b) && MPFR_IS_POS(c)) ? -1 : 1)
                : ((MPFR_IS_POS(b) && MPFR_IS_NEG(c)) ? 1 : -1);
              MPFR_SET_SIGN (a, sign);
              MPFR_SET_ZERO (a);
              MPFR_RET(0); /* 0 - 0 is exact */
            }
          else
            return mpfr_neg (a, c, rnd_mode);
        }
      else
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (c));
          return mpfr_set (a, b, rnd_mode);
        }
    }
  MPFR_CLEAR_FLAGS (a);
  MPFR_ASSERTD (MPFR_IS_PURE_FP (b) && MPFR_IS_PURE_FP (c));

  if (MPFR_LIKELY (MPFR_SIGN (b) == MPFR_SIGN (c)))
    { /* signs are equal, it's a real subtraction */
      if (MPFR_LIKELY (MPFR_PREC (a) == MPFR_PREC (b)
                       && MPFR_PREC (b) == MPFR_PREC (c)))
        return mpfr_sub1sp (a, b, c, rnd_mode);
      else
        return mpfr_sub1 (a, b, c, rnd_mode);
    }
  else
    { /* signs differ, it's an addition */
      if (MPFR_GET_EXP (b) < MPFR_GET_EXP (c))
         { /* exchange rounding modes toward +/- infinity */
          int inexact;
          rnd_mode = MPFR_INVERT_RND (rnd_mode);
          if (MPFR_LIKELY (MPFR_PREC (a) == MPFR_PREC (b)
                           && MPFR_PREC (b) == MPFR_PREC (c)))
            inexact = mpfr_add1sp (a, c, b, rnd_mode);
          else
            inexact = mpfr_add1 (a, c, b, rnd_mode);
          MPFR_CHANGE_SIGN (a);
          return -inexact;
        }
      else
        {
          if (MPFR_LIKELY (MPFR_PREC (a) == MPFR_PREC (b)
                           && MPFR_PREC (b) == MPFR_PREC (c)))
            return mpfr_add1sp (a, b, c, rnd_mode);
          else
            return mpfr_add1 (a, b, c, rnd_mode);
        }
    }
}

int
mpfr_tanh (mpfr_ptr y, mpfr_srcptr xt , mpfr_rnd_t rnd_mode)
{
  /****** Declaration ******/
  mpfr_t x;
  int inexact;
  MPFR_SAVE_EXPO_DECL (expo);

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

  /* Special value checking */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
    {
      if (MPFR_IS_NAN (xt))
        {
          MPFR_SET_NAN (y);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF (xt))
        {
          /* tanh(inf) = 1 && tanh(-inf) = -1 */
          return mpfr_set_si (y, MPFR_INT_SIGN (xt), rnd_mode);
        }
      else /* tanh (0) = 0 and xt is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO(xt));
          MPFR_SET_ZERO (y);
          MPFR_SET_SAME_SIGN (y, xt);
          MPFR_RET (0);
        }
    }

  /* tanh(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, 0,
                                    rnd_mode, {});

  MPFR_TMP_INIT_ABS (x, xt);

  MPFR_SAVE_EXPO_MARK (expo);

  /* General case */
  {
    /* Declaration of the intermediary variable */
    mpfr_t t, te;
    mpfr_exp_t d;

    /* Declaration of the size variable */
    mpfr_prec_t Ny = MPFR_PREC(y);   /* target precision */
    mpfr_prec_t Nt;                  /* working precision */
    long int err;                  /* error */
    int sign = MPFR_SIGN (xt);
    MPFR_ZIV_DECL (loop);
    MPFR_GROUP_DECL (group);

    /* First check for BIG overflow of exp(2*x):
       For x > 0, exp(2*x) > 2^(2*x)
       If 2 ^(2*x) > 2^emax or x>emax/2, there is an overflow */
    if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emax/2) >= 0)) {
      /* initialise of intermediary variables
         since 'set_one' label assumes the variables have been
         initialize */
      MPFR_GROUP_INIT_2 (group, MPFR_PREC_MIN, t, te);
      goto set_one;
    }

    /* Compute the precision of intermediary variable */
    /* The optimal number of bits: see algorithms.tex */
    Nt = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 4;
    /* if x is small, there will be a cancellation in exp(2x)-1 */
    if (MPFR_GET_EXP (x) < 0)
      Nt += -MPFR_GET_EXP (x);

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

    MPFR_ZIV_INIT (loop, Nt);
    for (;;) {
      /* tanh = (exp(2x)-1)/(exp(2x)+1) */
      mpfr_mul_2ui (te, x, 1, MPFR_RNDN);  /* 2x */
      /* since x > 0, we can only have an overflow */
      mpfr_exp (te, te, MPFR_RNDN);        /* exp(2x) */
      if (MPFR_UNLIKELY (MPFR_IS_INF (te))) {
      set_one:
        inexact = MPFR_FROM_SIGN_TO_INT (sign);
        mpfr_set4 (y, __gmpfr_one, MPFR_RNDN, sign);
        if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG_SIGN (sign)))
          {
            inexact = -inexact;
            mpfr_nexttozero (y);
          }
        break;
      }
      d = MPFR_GET_EXP (te);              /* For Error calculation */
      mpfr_add_ui (t, te, 1, MPFR_RNDD);   /* exp(2x) + 1*/
      mpfr_sub_ui (te, te, 1, MPFR_RNDU);  /* exp(2x) - 1*/
      d = d - MPFR_GET_EXP (te);
      mpfr_div (t, te, t, MPFR_RNDN);      /* (exp(2x)-1)/(exp(2x)+1)*/

      /* Calculation of the error */
//.........这里部分代码省略.........

/* The computation of z = pow(x,y) is done by
   z = exp(y * log(x)) = x^y
   For the special cases, see Section F.9.4.4 of the C standard:
     _ pow(±0, y) = ±inf for y an odd integer < 0.
     _ pow(±0, y) = +inf for y < 0 and not an odd integer.
     _ pow(±0, y) = ±0 for y an odd integer > 0.
     _ pow(±0, y) = +0 for y > 0 and not an odd integer.
     _ pow(-1, ±inf) = 1.
     _ pow(+1, y) = 1 for any y, even a NaN.
     _ pow(x, ±0) = 1 for any x, even a NaN.
     _ pow(x, y) = NaN for finite x < 0 and finite non-integer y.
     _ pow(x, -inf) = +inf for |x| < 1.
     _ pow(x, -inf) = +0 for |x| > 1.
     _ pow(x, +inf) = +0 for |x| < 1.
     _ pow(x, +inf) = +inf for |x| > 1.
     _ pow(-inf, y) = -0 for y an odd integer < 0.
     _ pow(-inf, y) = +0 for y < 0 and not an odd integer.
     _ pow(-inf, y) = -inf for y an odd integer > 0.
     _ pow(-inf, y) = +inf for y > 0 and not an odd integer.
     _ pow(+inf, y) = +0 for y < 0.
     _ pow(+inf, y) = +inf for y > 0. */
int
mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
{
  int inexact;
  int cmp_x_1;
  int y_is_integer;
  MPFR_SAVE_EXPO_DECL (expo);

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

  if (MPFR_ARE_SINGULAR (x, y))
    {
      /* pow(x, 0) returns 1 for any x, even a NaN. */
      if (MPFR_UNLIKELY (MPFR_IS_ZERO (y)))
        return mpfr_set_ui (z, 1, rnd_mode);
      else if (MPFR_IS_NAN (x))
        {
          MPFR_SET_NAN (z);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_NAN (y))
        {
          /* pow(+1, NaN) returns 1. */
          if (mpfr_cmp_ui (x, 1) == 0)
            return mpfr_set_ui (z, 1, rnd_mode);
          MPFR_SET_NAN (z);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF (y))
        {
          if (MPFR_IS_INF (x))
            {
              if (MPFR_IS_POS (y))
                MPFR_SET_INF (z);
              else
                MPFR_SET_ZERO (z);
              MPFR_SET_POS (z);
              MPFR_RET (0);
            }
          else
            {
              int cmp;
              cmp = mpfr_cmpabs (x, __gmpfr_one) * MPFR_INT_SIGN (y);
              MPFR_SET_POS (z);
              if (cmp > 0)
                {
                  /* Return +inf. */
                  MPFR_SET_INF (z);
                  MPFR_RET (0);
                }
              else if (cmp < 0)
                {
                  /* Return +0. */
                  MPFR_SET_ZERO (z);
                  MPFR_RET (0);
                }
              else
                {
                  /* Return 1. */
                  return mpfr_set_ui (z, 1, rnd_mode);
                }
            }
        }
      else if (MPFR_IS_INF (x))
        {
          int negative;
          /* Determine the sign now, in case y and z are the same object */
          negative = MPFR_IS_NEG (x) && is_odd (y);
          if (MPFR_IS_POS (y))
            MPFR_SET_INF (z);
          else
            MPFR_SET_ZERO (z);
          if (negative)
            MPFR_SET_NEG (z);
//.........这里部分代码省略.........

int
mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
{
  mpfr_t z_pre, s1, y, p;
  double sd, eps, m1, c;
  long add;
  mp_prec_t precz, prec1, precs, precs1;
  int inex;
  MPFR_GROUP_DECL (group);
  MPFR_ZIV_DECL (loop);
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_LOG_FUNC (("s[%#R]=%R rnd=%d", s, s, rnd_mode),
                 ("z[%#R]=%R inexact=%d", z, z, inex));

  /* Zero, Nan or Inf ? */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (s)))
    {
      if (MPFR_IS_NAN (s))
        {
          MPFR_SET_NAN (z);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF (s))
        {
          if (MPFR_IS_POS (s))
            return mpfr_set_ui (z, 1, GMP_RNDN); /* Zeta(+Inf) = 1 */
          MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */
          MPFR_RET_NAN;
        }
      else /* s iz zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (s));
          mpfr_set_ui (z, 1, rnd_mode);
          mpfr_div_2ui (z, z, 1, rnd_mode);
          MPFR_CHANGE_SIGN (z);
          MPFR_RET (0);
        }
    }

  /* s is neither Nan, nor Inf, nor Zero */

  /* check tiny s: we have zeta(s) = -1/2 - 1/2 log(2 Pi) s + ... around s=0,
     and for |s| <= 0.074, we have |zeta(s) + 1/2| <= |s|.
     Thus if |s| <= 1/4*ulp(1/2), we can deduce the correct rounding
     (the 1/4 covers the case where |zeta(s)| < 1/2 and rounding to nearest).
     A sufficient condition is that EXP(s) + 1 < -PREC(z). */
  if (MPFR_EXP(s) + 1 < - (mp_exp_t) MPFR_PREC(z))
    {
      int signs = MPFR_SIGN(s);
      mpfr_set_si_2exp (z, -1, -1, rnd_mode); /* -1/2 */
      if ((rnd_mode == GMP_RNDU || rnd_mode == GMP_RNDZ) && signs < 0)
        {
          mpfr_nextabove (z); /* z = -1/2 + epsilon */
          inex = 1;
        }
      else if (rnd_mode == GMP_RNDD && signs > 0)
        {
          mpfr_nextbelow (z); /* z = -1/2 - epsilon */
          inex = -1;
        }
      else
        {
          if (rnd_mode == GMP_RNDU) /* s > 0: z = -1/2 */
            inex = 1;
          else if (rnd_mode == GMP_RNDD)
            inex = -1;              /* s < 0: z = -1/2 */
          else /* (GMP_RNDZ and s > 0) or GMP_RNDN: z = -1/2 */
            inex = (signs > 0) ? 1 : -1;
        }
      return mpfr_check_range (z, inex, rnd_mode);
    }

  /* Check for case s= -2n */
  if (MPFR_IS_NEG (s))
    {
      mpfr_t tmp;
      tmp[0] = *s;
      MPFR_EXP (tmp) = MPFR_EXP (s) - 1;
      if (mpfr_integer_p (tmp))
        {
          MPFR_SET_ZERO (z);
          MPFR_SET_POS (z);
          MPFR_RET (0);
        }
    }

  MPFR_SAVE_EXPO_MARK (expo);

  /* Compute Zeta */
  if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */
    inex = mpfr_zeta_pos (z, s, rnd_mode);
  else /* use reflection formula
          zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */
    {
      precz = MPFR_PREC (z);
      precs = MPFR_PREC (s);

      /* Precision precs1 needed to represent 1 - s, and s + 2,
         without any truncation */
//.........这里部分代码省略.........

/* set f to the rational q */
int
mpfr_set_q (mpfr_ptr f, mpq_srcptr q, mpfr_rnd_t rnd)
{
  mpz_srcptr num, den;
  mpfr_t n, d;
  int inexact;
  int cn, cd;
  long shift;
  mp_size_t sn, sd;
  MPFR_SAVE_EXPO_DECL (expo);

  num = mpq_numref (q);
  den = mpq_denref (q);
  /* NAN and INF for mpq are not really documented, but could be found */
  if (MPFR_UNLIKELY (mpz_sgn (num) == 0))
    {
      if (MPFR_UNLIKELY (mpz_sgn (den) == 0))
        {
          MPFR_SET_NAN (f);
          MPFR_RET_NAN;
        }
      else
        {
          MPFR_SET_ZERO (f);
          MPFR_SET_POS (f);
          MPFR_RET (0);
        }
    }
  if (MPFR_UNLIKELY (mpz_sgn (den) == 0))
    {
      MPFR_SET_INF (f);
      MPFR_SET_SIGN (f, mpz_sgn (num));
      MPFR_RET (0);
    }

  MPFR_SAVE_EXPO_MARK (expo);

  cn = set_z (n, num, &sn);
  cd = set_z (d, den, &sd);

  sn -= sd;
  if (MPFR_UNLIKELY (sn > MPFR_EMAX_MAX / GMP_NUMB_BITS))
    {
      MPFR_SAVE_EXPO_FREE (expo);
      inexact = mpfr_overflow (f, rnd, MPFR_SIGN (f));
      goto end;
    }
  if (MPFR_UNLIKELY (sn < MPFR_EMIN_MIN / GMP_NUMB_BITS -1))
    {
      MPFR_SAVE_EXPO_FREE (expo);
      if (rnd == MPFR_RNDN)
        rnd = MPFR_RNDZ;
      inexact = mpfr_underflow (f, rnd, MPFR_SIGN (f));
      goto end;
    }

  inexact = mpfr_div (f, n, d, rnd);
  shift = GMP_NUMB_BITS*sn+cn-cd;
  MPFR_ASSERTD (shift == GMP_NUMB_BITS*sn+cn-cd);
  cd = mpfr_mul_2si (f, f, shift, rnd);
  MPFR_SAVE_EXPO_FREE (expo);
  if (MPFR_UNLIKELY (cd != 0))
    inexact = cd;
  else
    inexact = mpfr_check_range (f, inexact, rnd);
 end:
  mpfr_clear (d);
  mpfr_clear (n);
  MPFR_RET (inexact);
}

/* Set iop to the integral part of op and fop to its fractional part */
int
mpfr_modf (mpfr_ptr iop, mpfr_ptr fop, mpfr_srcptr op, mpfr_rnd_t rnd_mode)
{
  mpfr_exp_t ope;
  mpfr_prec_t opq;
  int inexi, inexf;

  MPFR_LOG_FUNC
    (("op[%Pu]=%.*Rg rnd=%d",
      mpfr_get_prec (op), mpfr_log_prec, op, rnd_mode),
     ("iop[%Pu]=%.*Rg fop[%Pu]=%.*Rg",
      mpfr_get_prec (iop), mpfr_log_prec, iop,
      mpfr_get_prec (fop), mpfr_log_prec, fop));

  MPFR_ASSERTN (iop != fop);

  if ( MPFR_UNLIKELY (MPFR_IS_SINGULAR (op)) )
    {
      if (MPFR_IS_NAN (op))
        {
          MPFR_SET_NAN (iop);
          MPFR_SET_NAN (fop);
          MPFR_RET_NAN;
        }
      MPFR_SET_SAME_SIGN (iop, op);
      MPFR_SET_SAME_SIGN (fop, op);
      if (MPFR_IS_INF (op))
        {
          MPFR_SET_INF (iop);
          MPFR_SET_ZERO (fop);
          MPFR_RET (0);
        }
      else /* op is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (op));
          MPFR_SET_ZERO (iop);
          MPFR_SET_ZERO (fop);
          MPFR_RET (0);
        }
    }

  ope = MPFR_GET_EXP (op);
  opq = MPFR_PREC (op);

  if (ope <= 0)   /* 0 < |op| < 1 */
    {
      inexf = (fop != op) ? mpfr_set (fop, op, rnd_mode) : 0;
      MPFR_SET_SAME_SIGN (iop, op);
      MPFR_SET_ZERO (iop);
      MPFR_RET (INEX(0, inexf));
    }
  else if (ope >= opq) /* op has no fractional part */
    {
      inexi = (iop != op) ? mpfr_set (iop, op, rnd_mode) : 0;
      MPFR_SET_SAME_SIGN (fop, op);
      MPFR_SET_ZERO (fop);
      MPFR_RET (INEX(inexi, 0));
    }
  else /* op has both integral and fractional parts */
    {
      if (iop != op)
        {
          inexi = mpfr_rint_trunc (iop, op, rnd_mode);
          inexf = mpfr_frac (fop, op, rnd_mode);
        }
      else
        {
          MPFR_ASSERTN (fop != op);
          inexf = mpfr_frac (fop, op, rnd_mode);
          inexi = mpfr_rint_trunc (iop, op, rnd_mode);
        }
      MPFR_RET (INEX(inexi, inexf));
    }
}

int
mpfr_ui_div (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  MPFR_LOG_FUNC
    (("u=%lu x[%Pu]=%.*Rg rnd=%d",
      u, mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode),
     ("y[%Pu]=%.*Rg", mpfr_get_prec(y), mpfr_log_prec, y));

  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)) /* u/Inf = 0 */
        {
          MPFR_SET_ZERO(y);
          MPFR_SET_SAME_SIGN(y,x);
          MPFR_RET(0);
        }
      else /* u / 0 */
        {
          MPFR_ASSERTD(MPFR_IS_ZERO(x));
          if (u)
            {
              /* u > 0, so y = sign(x) * Inf */
              MPFR_SET_SAME_SIGN(y, x);
              MPFR_SET_INF(y);
              MPFR_SET_DIVBY0 ();
              MPFR_RET(0);
            }
          else
            {
              /* 0 / 0 */
              MPFR_SET_NAN(y);
              MPFR_RET_NAN;
            }
        }
    }
  else if (MPFR_LIKELY(u != 0))
    {
      mpfr_t uu;
      mp_limb_t up[1];
      int 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[0] = (mp_limb_t) u << cnt;

      /* Optimization note: Exponent save/restore operations may be
         removed if mpfr_div works even when uu is out-of-range. */
      MPFR_SAVE_EXPO_MARK (expo);
      MPFR_SET_EXP (uu, GMP_NUMB_BITS - cnt);
      inex = mpfr_div (y, uu, x, rnd_mode);
      MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_check_range (y, inex, rnd_mode);
    }
  else /* u = 0, and x != 0 */
    {
      MPFR_SET_ZERO(y);         /* if u=0, then set y to 0 */
      MPFR_SET_SAME_SIGN(y, x); /* u considered as +0: sign(+0/x) = sign(x) */
      MPFR_RET(0);
    }
}

int
mpfr_log (mpfr_ptr r, mpfr_srcptr a, mpfr_rnd_t rnd_mode)
{
  int inexact;
  mpfr_prec_t p, q;
  mpfr_t tmp1, tmp2;
  MPFR_SAVE_EXPO_DECL (expo);
  MPFR_ZIV_DECL (loop);
  MPFR_GROUP_DECL(group);

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

  /* Special cases */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
    {
      /* If a is NaN, the result is NaN */
      if (MPFR_IS_NAN (a))
        {
          MPFR_SET_NAN (r);
          MPFR_RET_NAN;
        }
      /* check for infinity before zero */
      else if (MPFR_IS_INF (a))
        {
          if (MPFR_IS_NEG (a))
            /* log(-Inf) = NaN */
            {
              MPFR_SET_NAN (r);
              MPFR_RET_NAN;
            }
          else /* log(+Inf) = +Inf */
            {
              MPFR_SET_INF (r);
              MPFR_SET_POS (r);
              MPFR_RET (0);
            }
        }
      else /* a is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (a));
          MPFR_SET_INF (r);
          MPFR_SET_NEG (r);
          mpfr_set_divby0 ();
          MPFR_RET (0); /* log(0) is an exact -infinity */
        }
    }
  /* If a is negative, the result is NaN */
  else if (MPFR_UNLIKELY (MPFR_IS_NEG (a)))
    {
      MPFR_SET_NAN (r);
      MPFR_RET_NAN;
    }
  /* If a is 1, the result is 0 */
  else if (MPFR_UNLIKELY (MPFR_GET_EXP (a) == 1 && mpfr_cmp_ui (a, 1) == 0))
    {
      MPFR_SET_ZERO (r);
      MPFR_SET_POS (r);
      MPFR_RET (0); /* only "normal" case where the result is exact */
    }

  q = MPFR_PREC (r);

  /* use initial precision about q+lg(q)+5 */
  p = q + 5 + 2 * MPFR_INT_CEIL_LOG2 (q);
  /* % ~(mpfr_prec_t)GMP_NUMB_BITS  ;
     m=q; while (m) { p++; m >>= 1; }  */
  /* if (MPFR_LIKELY(p % GMP_NUMB_BITS != 0))
      p += GMP_NUMB_BITS - (p%GMP_NUMB_BITS); */

  MPFR_SAVE_EXPO_MARK (expo);
  MPFR_GROUP_INIT_2 (group, p, tmp1, tmp2);

  MPFR_ZIV_INIT (loop, p);
  for (;;)
    {
      long m;
      mpfr_exp_t cancel;

      /* Calculus of m (depends on p) */
      m = (p + 1) / 2 - MPFR_GET_EXP (a) + 1;

      mpfr_mul_2si (tmp2, a, m, MPFR_RNDN);    /* s=a*2^m,        err<=1 ulp  */
      mpfr_div (tmp1, __gmpfr_four, tmp2, MPFR_RNDN);/* 4/s,      err<=2 ulps */
      mpfr_agm (tmp2, __gmpfr_one, tmp1, MPFR_RNDN); /* AG(1,4/s),err<=3 ulps */
      mpfr_mul_2ui (tmp2, tmp2, 1, MPFR_RNDN); /* 2*AG(1,4/s),    err<=3 ulps */
      mpfr_const_pi (tmp1, MPFR_RNDN);         /* compute pi,     err<=1ulp   */
      mpfr_div (tmp2, tmp1, tmp