C++ MPFR_ZIV_INIT函数代码示例

int
mpfr_const_euler_internal (mpfr_t x, mpfr_rnd_t rnd)
{
  mpfr_prec_t prec = MPFR_PREC(x), m, log2m;
  mpfr_t y, z;
  unsigned long n;
  int inexact;
  MPFR_ZIV_DECL (loop);

  log2m = MPFR_INT_CEIL_LOG2 (prec);
  m = prec + 2 * log2m + 23;

  mpfr_init2 (y, m);
  mpfr_init2 (z, m);

  MPFR_ZIV_INIT (loop, m);
  for (;;)
    {
      mpfr_exp_t exp_S, err;
      /* since prec >= 1, we have m >= 24 here, which ensures n >= 9 below */
      n = 1 + (unsigned long) ((double) m * LOG2 / 2.0);
      MPFR_ASSERTD (n >= 9);
      mpfr_const_euler_S2 (y, n); /* error <= 3 ulps */
      exp_S = MPFR_EXP(y);
      mpfr_set_ui (z, n, MPFR_RNDN);
      mpfr_log (z, z, MPFR_RNDD); /* error <= 1 ulp */
      mpfr_sub (y, y, z, MPFR_RNDN); /* S'(n) - log(n) */
      /* the error is less than 1/2 + 3*2^(exp_S-EXP(y)) + 2^(EXP(z)-EXP(y))
         <= 1/2 + 2^(exp_S+2-EXP(y)) + 2^(EXP(z)-EXP(y))
         <= 1/2 + 2^(1+MAX(exp_S+2,EXP(z))-EXP(y)) */
      err = 1 + MAX(exp_S + 2, MPFR_EXP(z)) - MPFR_EXP(y);
      err = (err >= -1) ? err + 1 : 0; /* error <= 2^err ulp(y) */
      exp_S = MPFR_EXP(y);
      mpfr_const_euler_R (z, n); /* err <= ulp(1/2) = 2^(-m) */
      mpfr_sub (y, y, z, MPFR_RNDN);
      /* err <= 1/2 ulp(y) + 2^(-m) + 2^(err + exp_S - EXP(y)) ulp(y).
         Since the result is between 0.5 and 1, ulp(y) = 2^(-m).
         So we get 3/2*ulp(y) + 2^(err + exp_S - EXP(y)) ulp(y).
         3/2 + 2^e <= 2^(e+1) for e>=1, and <= 2^2 otherwise */
      err = err + exp_S - MPFR_EXP(y);
      err = (err >= 1) ? err + 1 : 2;
      if (MPFR_LIKELY (MPFR_CAN_ROUND (y, m - err, prec, rnd)))
        break;
      MPFR_ZIV_NEXT (loop, m);
      mpfr_set_prec (y, m);
      mpfr_set_prec (z, m);
    }
  MPFR_ZIV_FREE (loop);

  inexact = mpfr_set (x, y, rnd);

  mpfr_clear (y);
  mpfr_clear (z);

  return inexact; /* always inexact */
}

/* computes tan(x) = sign(x)*sqrt(1/cos(x)^2-1) */
int
mpfr_tan (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
{
  mp_prec_t precy, m;
  int inexact;
  mpfr_t s, c;
  MPFR_ZIV_DECL (loop);
  MPFR_SAVE_EXPO_DECL (expo);
  MPFR_GROUP_DECL (group);

  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_IS_INF(x))
        {
          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);
        }
    }

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

  MPFR_SAVE_EXPO_MARK (expo);

  /* Compute initial precision */
  precy = MPFR_PREC (y);
  m = precy + MPFR_INT_CEIL_LOG2 (precy) + 13;
  MPFR_ASSERTD (m >= 2); /* needed for the error analysis in algorithms.tex */

  MPFR_GROUP_INIT_2 (group, m, s, c);
  MPFR_ZIV_INIT (loop, m);
  for (;;)
    {
      /* The only way to get an overflow is to get ~ Pi/2
         But the result will be ~ 2^Prec(y). */
      mpfr_sin_cos (s, c, x, GMP_RNDN); /* err <= 1/2 ulp on s and c */
      mpfr_div (c, s, c, GMP_RNDN);     /* err <= 4 ulps */
      MPFR_ASSERTD (!MPFR_IS_SINGULAR (c));
      if (MPFR_LIKELY (MPFR_CAN_ROUND (c, m - 2, precy, rnd_mode)))
        break;
      MPFR_ZIV_NEXT (loop, m);
      MPFR_GROUP_REPREC_2 (group, m, s, c);
    }
  MPFR_ZIV_FREE (loop);
  inexact = mpfr_set (y, c, rnd_mode);
  MPFR_GROUP_CLEAR (group);

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (y, inexact, 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);
//.........这里部分代码省略.........

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

  MPFR_LOG_FUNC
    (("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));

  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))
    {
      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)) */

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

/* Don't need to save / restore exponent range: the cache does it */
int
mpfr_const_log2_internal (mpfr_ptr x, mpfr_rnd_t rnd_mode)
{
  unsigned long n = MPFR_PREC (x);
  mpfr_prec_t w; /* working precision */
  unsigned long N;
  mpz_t *T, *P, *Q;
  mpfr_t t, q;
  int inexact;
  int ok = 1; /* ensures that the 1st try will give correct rounding */
  unsigned long lgN, i;
  MPFR_GROUP_DECL(group);
  MPFR_TMP_DECL(marker);
  MPFR_ZIV_DECL(loop);

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

  if (n < 1253)
    w = n + 10; /* ensures correct rounding for the four rounding modes,
                   together with N = w / 3 + 1 (see below). */
  else if (n < 2571)
    w = n + 11; /* idem */
  else if (n < 3983)
    w = n + 12;
  else if (n < 4854)
    w = n + 13;
  else if (n < 26248)
    w = n + 14;
  else
    {
      w = n + 15;
      ok = 0;
    }

  MPFR_TMP_MARK(marker);
  MPFR_GROUP_INIT_2(group, w, t, q);

  MPFR_ZIV_INIT (loop, w);
  for (;;)
    {
      N = w / 3 + 1; /* Warning: do not change that (even increasing N!)
                        without checking correct rounding in the above
                        ranges for n. */

      /* the following are needed for error analysis (see algorithms.tex) */
      MPFR_ASSERTD(w >= 3 && N >= 2);

      lgN = MPFR_INT_CEIL_LOG2 (N) + 1;
      T  = (mpz_t *) MPFR_TMP_ALLOC (3 * lgN * sizeof (mpz_t));
      P  = T + lgN;
      Q  = T + 2*lgN;
      for (i = 0; i < lgN; i++)
        {
          mpz_init (T[i]);
          mpz_init (P[i]);
          mpz_init (Q[i]);
        }

      S (T, P, Q, 0, N, 0);

      mpfr_set_z (t, T[0], MPFR_RNDN);
      mpfr_set_z (q, Q[0], MPFR_RNDN);
      mpfr_div (t, t, q, MPFR_RNDN);

      for (i = 0; i < lgN; i++)
        {
          mpz_clear (T[i]);
          mpz_clear (P[i]);
          mpz_clear (Q[i]);
        }

      if (MPFR_LIKELY (ok != 0
                       || mpfr_can_round (t, w - 2, MPFR_RNDN, rnd_mode, n)))
        break;

      MPFR_ZIV_NEXT (loop, w);
      MPFR_GROUP_REPREC_2(group, w, t, q);
    }
  MPFR_ZIV_FREE (loop);

  inexact = mpfr_set (x, t, rnd_mode);

  MPFR_GROUP_CLEAR(group);
  MPFR_TMP_FREE(marker);

  return inexact;
}

/* (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, mp_rnd_t rnd_mode)
{
  mp_prec_t prec, m;
  int neg, reduce;
  mpfr_t c, xr;
  mpfr_srcptr xx;
  mp_exp_t err, expx;
  MPFR_ZIV_DECL (loop);

  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 */
          return mpfr_set_ui (z, 1, rnd_mode);
        }
    }

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

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

  mpfr_init (c);
  mpfr_init (xr);

  MPFR_ZIV_INIT (loop, m);
  for (;;)
    {
      /* the following is copied from sin.c */
      if (expx >= 2) /* reduce the argument */
        {
          reduce = 1;
          mpfr_set_prec (c, expx + m - 1);
          mpfr_set_prec (xr, m);
          mpfr_const_pi (c, GMP_RNDN);
          mpfr_mul_2ui (c, c, 1, GMP_RNDN);
          mpfr_remainder (xr, x, c, GMP_RNDN);
          mpfr_div_2ui (c, c, 1, GMP_RNDN);
          if (MPFR_SIGN (xr) > 0)
            mpfr_sub (c, c, xr, GMP_RNDZ);
          else
            mpfr_add (c, c, xr, GMP_RNDZ);
          if (MPFR_IS_ZERO(xr) || MPFR_EXP(xr) < (mp_exp_t) 3 - (mp_exp_t) m
              || MPFR_EXP(c) < (mp_exp_t) 3 - (mp_exp_t) m)
            goto next_step;
          xx = xr;
        }
      else /* the input argument is already reduced */
        {
          reduce = 0;
          xx = x;
        }

      neg = MPFR_IS_NEG (xx); /* gives sign of sin(x) */
      mpfr_set_prec (c, m);
      mpfr_cos (c, xx, GMP_RNDZ);
      /* If no argument reduction was performed, the error is at most ulp(c),
         otherwise it is at most ulp(c) + 2^(2-m). Since |c| < 1, we have
         ulp(c) <= 2^(-m), thus the error is bounded by 2^(3-m) in that later
         case. */
      if (reduce == 0)
        err = m;
      else
        err = MPFR_GET_EXP (c) + (mp_exp_t) (m - 3);
      if (!mpfr_can_round (c, err, GMP_RNDN, rnd_mode,
                           MPFR_PREC (z) + (rnd_mode == GMP_RNDN)))
        goto next_step;

      mpfr_set (z, c, rnd_mode);
      mpfr_sqr (c, c, GMP_RNDU);
      mpfr_ui_sub (c, 1, c, GMP_RNDN);
      err = 2 + (- MPFR_GET_EXP (c)) / 2;
      mpfr_sqrt (c, c, GMP_RNDN);
      if (neg)
        MPFR_CHANGE_SIGN (c);

      /* the absolute error on c is at most 2^(err-m), which we must put
         in the form 2^(EXP(c)-err). If there was an argument reduction,
         we need to add 2^(2-m); since err >= 2, the error is bounded by
         2^(err+1-m) in that case. */
      err = MPFR_GET_EXP (c) + (mp_exp_t) m - (err + reduce);
      if (mpfr_can_round (c, err, GMP_RNDN, rnd_mode,
                          MPFR_PREC (y) + (rnd_mode == GMP_RNDN)))
        break;
//.........这里部分代码省略.........

int
mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
{
  int comp, inexact;
  mp_exp_t ex;
  MPFR_SAVE_EXPO_DECL (expo);

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN (x))
        {
          MPFR_SET_NAN (y);
          MPFR_RET_NAN;
        }
      /* check for inf or -inf (result is not defined) */
      else if (MPFR_IS_INF (x))
        {
          if (MPFR_IS_POS (x))
            {
              MPFR_SET_INF (y);
              MPFR_SET_POS (y);
              MPFR_RET (0);
            }
          else
            {
              MPFR_SET_NAN (y);
              MPFR_RET_NAN;
            }
        }
      else /* x is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
          MPFR_SET_ZERO (y);   /* log1p(+/- 0) = +/- 0 */
          MPFR_SET_SAME_SIGN (y, x);
          MPFR_RET (0);
        }
    }

  ex = MPFR_GET_EXP (x);
  if (ex < 0)  /* -0.5 < x < 0.5 */
    {
      /* For x > 0,    abs(log(1+x)-x) < x^2/2.
         For x > -0.5, abs(log(1+x)-x) < x^2. */
      if (MPFR_IS_POS (x))
        MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex - 1, 0, 0, rnd_mode, {});
      else
        MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 0, 1, rnd_mode, {});
    }

  comp = mpfr_cmp_si (x, -1);
  /* log1p(x) is undefined for x < -1 */
  if (MPFR_UNLIKELY(comp <= 0))
    {
      if (comp == 0)
        /* x=0: log1p(-1)=-inf (division by zero) */
        {
          MPFR_SET_INF (y);
          MPFR_SET_NEG (y);
          MPFR_RET (0);
        }
      MPFR_SET_NAN (y);
      MPFR_RET_NAN;
    }

  MPFR_SAVE_EXPO_MARK (expo);

  /* General case */
  {
    /* 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 + MPFR_INT_CEIL_LOG2 (Ny) + 6;

    /* if |x| is smaller than 2^(-e), we will loose about e bits
       in log(1+x) */
    if (MPFR_EXP(x) < 0)
      Nt += -MPFR_EXP(x);

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

    /* First computation of log1p */
    MPFR_ZIV_INIT (loop, Nt);
    for (;;)
      {
        /* compute log1p */
        inexact = mpfr_add_ui (t, x, 1, GMP_RNDN);      /* 1+x */
        /* if inexact = 0, then t = x+1, and the result is simply log(t) */
        if (inexact == 0)
          {
            inexact = mpfr_log (y, t, rnd_mode);
            goto end;
          }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
           *     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).
           */
          expy =
            n != -1 && expx > 0 && expx > (__gmpfr_emin - 1) / n ?
            MPFR_EMIN_MIN - 2 /* Underflow */ :
            n != -1 && expx < 0 && expx < (__gmpfr_emax - 1) / n ?
            MPFR_EMAX_MAX /* Overflow */ : n * expx;
          return mpfr_set_si_2exp (y, n % 2 ? MPFR_INT_SIGN (x) : 1,
                                   expy, rnd);
        }

      /* General case */
      {
        /* 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 */
        int inexact;
        unsigned long abs_n;
        MPFR_SAVE_EXPO_DECL (expo);
        MPFR_ZIV_DECL (loop);

        abs_n = - (unsigned long) n;

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

        MPFR_SAVE_EXPO_MARK (expo);

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

        MPFR_ZIV_INIT (loop, Nt);
        for (;;)
          {
            /* compute 1/(x^n), with n > 0 */
            mpfr_pow_ui (t, x, abs_n, GMP_RNDN);
            mpfr_ui_div (t, 1, t, GMP_RNDN);
            /* FIXME: old code improved, but I think this is still incorrect. */
            if (MPFR_UNLIKELY (MPFR_IS_ZERO (t)))
              {
                MPFR_ZIV_FREE (loop);
                mpfr_clear (t);
                MPFR_SAVE_EXPO_FREE (expo);
                return mpfr_underflow (y, rnd == GMP_RNDN ? GMP_RNDZ : rnd,
                                       abs_n & 1 ? MPFR_SIGN (x) :
                                       MPFR_SIGN_POS);
              }
            if (MPFR_UNLIKELY (MPFR_IS_INF (t)))
              {
                MPFR_ZIV_FREE (loop);
                mpfr_clear (t);
                MPFR_SAVE_EXPO_FREE (expo);
                return mpfr_overflow (y, rnd, abs_n & 1 ? MPFR_SIGN (x) :
                                      MPFR_SIGN_POS);
              }
            /* error estimate -- see pow function in algorithms.ps */
            err = Nt - 3;
            if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd)))
              break;

            /* actualisation of the precision */
            Nt += BITS_PER_MP_LIMB;
            mpfr_set_prec (t, Nt);
          }
        MPFR_ZIV_FREE (loop);

        inexact = mpfr_set (y, t, rnd);
        mpfr_clear (t);
        MPFR_SAVE_EXPO_FREE (expo);
        return mpfr_check_range (y, inexact, rnd);
      }
    }
}

int
mpfr_sinh_cosh (mpfr_ptr sh, mpfr_ptr ch, mpfr_srcptr xt, mpfr_rnd_t rnd_mode)
{
    mpfr_t x;
    int inexact_sh, inexact_ch;

    MPFR_ASSERTN (sh != ch);

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

    if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
    {
        if (MPFR_IS_NAN (xt))
        {
            MPFR_SET_NAN (ch);
            MPFR_SET_NAN (sh);
            MPFR_RET_NAN;
        }
        else if (MPFR_IS_INF (xt))
        {
            MPFR_SET_INF (sh);
            MPFR_SET_SAME_SIGN (sh, xt);
            MPFR_SET_INF (ch);
            MPFR_SET_POS (ch);
            MPFR_RET (0);
        }
        else /* xt is zero */
        {
            MPFR_ASSERTD (MPFR_IS_ZERO (xt));
            MPFR_SET_ZERO (sh);                   /* sinh(0) = 0 */
            MPFR_SET_SAME_SIGN (sh, xt);
            inexact_sh = 0;
            inexact_ch = mpfr_set_ui (ch, 1, rnd_mode); /* cosh(0) = 1 */
            return INEX(inexact_sh,inexact_ch);
        }
    }

    /* Warning: if we use MPFR_FAST_COMPUTE_IF_SMALL_INPUT here, make sure
       that the code also works in case of overlap (see sin_cos.c) */

    MPFR_TMP_INIT_ABS (x, xt);

    {
        mpfr_t s, c, ti;
        mpfr_exp_t d;
        mpfr_prec_t N;    /* Precision of the intermediary variables */
        long int err;    /* Precision of error */
        MPFR_ZIV_DECL (loop);
        MPFR_SAVE_EXPO_DECL (expo);
        MPFR_GROUP_DECL (group);

        MPFR_SAVE_EXPO_MARK (expo);

        /* compute the precision of intermediary variable */
        N = MPFR_PREC (ch);
        N = MAX (N, MPFR_PREC (sh));
        /* the optimal number of bits : see algorithms.ps */
        N = N + MPFR_INT_CEIL_LOG2 (N) + 4;

        /* initialise of intermediary variables */
        MPFR_GROUP_INIT_3 (group, N, s, c, ti);

        /* First computation of sinh_cosh */
        MPFR_ZIV_INIT (loop, N);
        for (;;)
        {
            MPFR_BLOCK_DECL (flags);

            /* compute sinh_cosh */
            MPFR_BLOCK (flags, mpfr_exp (s, x, MPFR_RNDD));
            if (MPFR_OVERFLOW (flags))
                /* exp(x) does overflow */
            {
                /* since cosh(x) >= exp(x), cosh(x) overflows too */
                inexact_ch = mpfr_overflow (ch, rnd_mode, MPFR_SIGN_POS);
                /* sinh(x) may be representable */
                inexact_sh = mpfr_sinh (sh, xt, rnd_mode);
                MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
                break;
            }
            d = MPFR_GET_EXP (s);
            mpfr_ui_div (ti, 1, s, MPFR_RNDU);  /* 1/exp(x) */
            mpfr_add (c, s, ti, MPFR_RNDU);     /* exp(x) + 1/exp(x) */
            mpfr_sub (s, s, ti, MPFR_RNDN);     /* exp(x) - 1/exp(x) */
            mpfr_div_2ui (c, c, 1, MPFR_RNDN);  /* 1/2(exp(x) + 1/exp(x)) */
            mpfr_div_2ui (s, s, 1, MPFR_RNDN);  /* 1/2(exp(x) - 1/exp(x)) */

            /* it may be that s is zero (in fact, it can only occur when exp(x)=1,
               and thus ti=1 too) */
            if (MPFR_IS_ZERO (s))
                err = N; /* double the precision */
            else
            {
                /* calculation of the error */
                d = d - MPFR_GET_EXP (s) + 2;
//.........这里部分代码省略.........

int
mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_prec_t K0, K, precy, m, k, l;
  int inexact, reduce = 0;
  mpfr_t r, s, xr, c;
  mpfr_exp_t exps, cancel = 0, expx;
  MPFR_ZIV_DECL (loop);
  MPFR_SAVE_EXPO_DECL (expo);
  MPFR_GROUP_DECL (group);

  MPFR_LOG_FUNC (
    ("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));

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
        {
          MPFR_SET_NAN (y);
          MPFR_RET_NAN;
        }
      else
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
          return mpfr_set_ui (y, 1, rnd_mode);
        }
    }

  MPFR_SAVE_EXPO_MARK (expo);

  /* cos(x) = 1-x^2/2 + ..., so error < 2^(2*EXP(x)-1) */
  expx = MPFR_GET_EXP (x);
  MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, __gmpfr_one, -2 * expx,
                                    1, 0, rnd_mode, expo, {});

  /* Compute initial precision */
  precy = MPFR_PREC (y);

  if (precy >= MPFR_SINCOS_THRESHOLD)
    {
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_cos_fast (y, x, rnd_mode);
    }

  K0 = __gmpfr_isqrt (precy / 3);
  m = precy + 2 * MPFR_INT_CEIL_LOG2 (precy) + 2 * K0;

  if (expx >= 3)
    {
      reduce = 1;
      /* As expx + m - 1 will silently be converted into mpfr_prec_t
         in the mpfr_init2 call, the assert below may be useful to
         avoid undefined behavior. */
      MPFR_ASSERTN (expx + m - 1 <= MPFR_PREC_MAX);
      mpfr_init2 (c, expx + m - 1);
      mpfr_init2 (xr, m);
    }

  MPFR_GROUP_INIT_2 (group, m, r, s);
  MPFR_ZIV_INIT (loop, m);
  for (;;)
    {
      /* If |x| >= 4, first reduce x cmod (2*Pi) into xr, using mpfr_remainder:
         let e = EXP(x) >= 3, and m the target precision:
         (1) c <- 2*Pi              [precision e+m-1, nearest]
         (2) xr <- remainder (x, c) [precision m, nearest]
         We have |c - 2*Pi| <= 1/2ulp(c) = 2^(3-e-m)
                 |xr - x - k c| <= 1/2ulp(xr) <= 2^(1-m)
                 |k| <= |x|/(2*Pi) <= 2^(e-2)
         Thus |xr - x - 2kPi| <= |k| |c - 2Pi| + 2^(1-m) <= 2^(2-m).
         It follows |cos(xr) - cos(x)| <= 2^(2-m). */
      if (reduce)
        {
          mpfr_const_pi (c, MPFR_RNDN);
          mpfr_mul_2ui (c, c, 1, MPFR_RNDN); /* 2Pi */
          mpfr_remainder (xr, x, c, MPFR_RNDN);
          if (MPFR_IS_ZERO(xr))
            goto ziv_next;
          /* now |xr| <= 4, thus r <= 16 below */
          mpfr_mul (r, xr, xr, MPFR_RNDU); /* err <= 1 ulp */
        }
      else
        mpfr_mul (r, x, x, MPFR_RNDU); /* err <= 1 ulp */

      /* now |x| < 4 (or xr if reduce = 1), thus |r| <= 16 */

      /* we need |r| < 1/2 for mpfr_cos2_aux, i.e., EXP(r) - 2K <= -1 */
      K = K0 + 1 + MAX(0, MPFR_GET_EXP(r)) / 2;
      /* since K0 >= 0, if EXP(r) < 0, then K >= 1, thus EXP(r) - 2K <= -3;
         otherwise if EXP(r) >= 0, then K >= 1/2 + EXP(r)/2, thus
         EXP(r) - 2K <= -1 */

      MPFR_SET_EXP (r, MPFR_GET_EXP (r) - 2 * K); /* Can't overflow! */

      /* s <- 1 - r/2! + ... + (-1)^l r^l/(2l)! */
      l = mpfr_cos2_aux (s, r);
      /* l is the error bound in ulps on s */
      MPFR_SET_ONE (r);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      /* General case */
      {
        /* Declaration of the intermediary variable */
        mpfr_t t;
        /* Declaration of the size variable */
        mpfr_prec_t Ny;                              /* target precision */
        mpfr_prec_t Nt;                              /* working precision */
        mpfr_rnd_t rnd1;
        int size_n;
        int inexact;
        unsigned long abs_n;
        MPFR_SAVE_EXPO_DECL (expo);
        MPFR_ZIV_DECL (loop);

        abs_n = - (unsigned long) n;
        count_leading_zeros (size_n, (mp_limb_t) abs_n);
        size_n = GMP_NUMB_BITS - size_n;

        /* initial working precision */
        Ny = MPFR_PREC (y);
        Nt = Ny + size_n + 3 + MPFR_INT_CEIL_LOG2 (Ny);

        MPFR_SAVE_EXPO_MARK (expo);

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

        /* We will compute rnd(rnd1(1/x) ^ |n|), where rnd1 is the rounding
           toward sign(x), to avoid spurious overflow or underflow, as in
           mpfr_pow_z. */
        rnd1 = MPFR_EXP (x) < 1 ? MPFR_RNDZ :
          (MPFR_SIGN (x) > 0 ? MPFR_RNDU : MPFR_RNDD);

        MPFR_ZIV_INIT (loop, Nt);
        for (;;)
          {
            MPFR_BLOCK_DECL (flags);

            /* compute (1/x)^|n| */
            MPFR_BLOCK (flags, mpfr_ui_div (t, 1, x, rnd1));
            MPFR_ASSERTD (! MPFR_UNDERFLOW (flags));
            /* t = (1/x)*(1+theta) where |theta| <= 2^(-Nt) */
            if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
              goto overflow;
            MPFR_BLOCK (flags, mpfr_pow_ui (t, t, abs_n, rnd));
            /* t = (1/x)^|n|*(1+theta')^(|n|+1) where |theta'| <= 2^(-Nt).
               If (|n|+1)*2^(-Nt) <= 1/2, which is satisfied as soon as
               Nt >= bits(n)+2, then we can use Lemma \ref{lemma_graillat}
               from algorithms.tex, which yields x^n*(1+theta) with
               |theta| <= 2(|n|+1)*2^(-Nt), thus the error is bounded by
               2(|n|+1) ulps <= 2^(bits(n)+2) ulps. */
            if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
              {
              overflow:
                MPFR_ZIV_FREE (loop);
                mpfr_clear (t);
                MPFR_SAVE_EXPO_FREE (expo);
                MPFR_LOG_MSG (("overflow\n", 0));
                return mpfr_overflow (y, rnd, abs_n & 1 ?
                                      MPFR_SIGN (x) : MPFR_SIGN_POS);
              }
            if (MPFR_UNLIKELY (MPFR_UNDERFLOW (flags)))
              {
                MPFR_ZIV_FREE (loop);
                mpfr_clear (t);
                MPFR_LOG_MSG (("underflow\n", 0));

int
mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t c, xr;
  mpfr_srcptr xx;
  mpfr_exp_t expx, err;
  mpfr_prec_t precy, m;
  int inexact, sign, reduce;
  MPFR_ZIV_DECL (loop);
  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_IS_INF (x))
        {
          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);
        }
    }

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

  MPFR_SAVE_EXPO_MARK (expo);

  /* Compute initial precision */
  precy = MPFR_PREC (y);

  if (precy >= MPFR_SINCOS_THRESHOLD)
    return mpfr_sin_fast (y, x, rnd_mode);

  m = precy + MPFR_INT_CEIL_LOG2 (precy) + 13;
  expx = MPFR_GET_EXP (x);

  mpfr_init (c);
  mpfr_init (xr);

  MPFR_ZIV_INIT (loop, m);
  for (;;)
    {
      /* first perform argument reduction modulo 2*Pi (if needed),
         also helps to determine the sign of sin(x) */
      if (expx >= 2) /* If Pi < x < 4, we need to reduce too, to determine
                        the sign of sin(x). For 2 <= |x| < Pi, we could avoid
                        the reduction. */
        {
          reduce = 1;
          /* As expx + m - 1 will silently be converted into mpfr_prec_t
             in the mpfr_set_prec call, the assert below may be useful to
             avoid undefined behavior. */
          MPFR_ASSERTN (expx + m - 1 <= MPFR_PREC_MAX);
          mpfr_set_prec (c, expx + m - 1);
          mpfr_set_prec (xr, m);
          mpfr_const_pi (c, MPFR_RNDN);
          mpfr_mul_2ui (c, c, 1, MPFR_RNDN);
          mpfr_remainder (xr, x, c, MPFR_RNDN);
          /* The analysis is similar to that of cos.c:
             |xr - x - 2kPi| <= 2^(2-m). Thus we can decide the sign
             of sin(x) if xr is at distance at least 2^(2-m) of both
             0 and +/-Pi. */
          mpfr_div_2ui (c, c, 1, MPFR_RNDN);
          /* Since c approximates Pi with an error <= 2^(2-expx-m) <= 2^(-m),
             it suffices to check that c - |xr| >= 2^(2-m). */
          if (MPFR_SIGN (xr) > 0)
            mpfr_sub (c, c, xr, MPFR_RNDZ);
          else
            mpfr_add (c, c, xr, MPFR_RNDZ);
          if (MPFR_IS_ZERO(xr)
              || MPFR_EXP(xr) < (mpfr_exp_t) 3 - (mpfr_exp_t) m
              || MPFR_EXP(c) < (mpfr_exp_t) 3 - (mpfr_exp_t) m)
            goto ziv_next;

          /* |xr - x - 2kPi| <= 2^(2-m), thus |sin(xr) - sin(x)| <= 2^(2-m) */
          xx = xr;
        }
      else /* the input argument is already reduced */
        {
          reduce = 0;
          xx = x;
        }

      sign = MPFR_SIGN(xx);
      /* now that the argument is reduced, precision m is enough */
      mpfr_set_prec (c, m);
      mpfr_cos (c, xx, MPFR_RNDZ);    /* can't be exact */
      mpfr_nexttoinf (c);           /* now c = cos(x) rounded away */
      mpfr_mul (c, c, c, MPFR_RNDU); /* away */
      mpfr_ui_sub (c, 1, c, MPFR_RNDZ);
//.........这里部分代码省略.........

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, tmp2, MPFR_RNDN);  /* pi/2*AG(1,4/s), err<=5ulps  */
      mpfr_const_log2 (tmp1, MPFR_RNDN);      /* compute log(2),  err<=1ulp   */
      mpfr_mul_si (tmp1, tmp1, m, MPFR_RNDN); /* compute m*log(2),err<=2ulps  */
      mpfr_sub (tmp1, tmp2, tmp1, MPFR_RNDN); /* log(a),    err<=7ulps+cancel */

      if (MPFR_LIKELY (MPFR_IS_PURE_FP (tmp1) && MPFR_IS_PURE_FP (tmp2)))
        {
          cancel = MPFR_GET_EXP (tmp2) - MPFR_GET_EXP (tmp1);
          MPFR_LOG_MSG (("canceled bits=%ld\n", (long) cancel));
          MPFR_LOG_VAR (tmp1);
          if (MPFR_UNLIKELY (cancel < 0))
//.........这里部分代码省略.........

int
mpfr_ui_pow_ui (mpfr_ptr x, unsigned long int y, unsigned long int n,
                mpfr_rnd_t rnd)
{
  mpfr_exp_t err;
  unsigned long m;
  mpfr_t res;
  mpfr_prec_t prec;
  int size_n;
  int inexact;
  MPFR_ZIV_DECL (loop);
  MPFR_SAVE_EXPO_DECL (expo);

  if (MPFR_UNLIKELY (n <= 1))
    {
      if (n == 1)
        return mpfr_set_ui (x, y, rnd);     /* y^1 = y */
      else
        return mpfr_set_ui (x, 1, rnd);     /* y^0 = 1 for any y */
    }
  else if (MPFR_UNLIKELY (y <= 1))
    {
      if (y == 1)
        return mpfr_set_ui (x, 1, rnd);     /* 1^n = 1 for any n > 0 */
      else
        return mpfr_set_ui (x, 0, rnd);     /* 0^n = 0 for any n > 0 */
    }

  for (size_n = 0, m = n; m; size_n++, m >>= 1);

  MPFR_SAVE_EXPO_MARK (expo);
  prec = MPFR_PREC (x) + 3 + size_n;
  mpfr_init2 (res, prec);

  MPFR_ZIV_INIT (loop, prec);
  for (;;)
    {
      int i = size_n;

      inexact = mpfr_set_ui (res, y, MPFR_RNDU);
      err = 1;
      /* now 2^(i-1) <= n < 2^i: i=1+floor(log2(n)) */
      for (i -= 2; i >= 0; i--)
        {
          inexact |= mpfr_mul (res, res, res, MPFR_RNDU);
          err++;
          if (n & (1UL << i))
            inexact |= mpfr_mul_ui (res, res, y, MPFR_RNDU);
        }
      /* since the loop is executed floor(log2(n)) times,
         we have err = 1+floor(log2(n)).
         Since prec >= MPFR_PREC(x) + 4 + floor(log2(n)), prec > err */
      err = prec - err;

      if (MPFR_LIKELY (inexact == 0
                       || MPFR_CAN_ROUND (res, err, MPFR_PREC (x), rnd)))
        break;

      /* Actualisation of the precision */
      MPFR_ZIV_NEXT (loop, prec);
      mpfr_set_prec (res, prec);
    }
  MPFR_ZIV_FREE (loop);

  inexact = mpfr_set (x, res, rnd);

  mpfr_clear (res);

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (x, inexact, rnd);
}

int
mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mpfr_rnd_t rnd_mode)
{
  int inexact;
  MPFR_SAVE_EXPO_DECL (expo);

  /* If a is NaN, the result is NaN */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
    {
      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))
            /* log10(-Inf) = NaN */
            {
              MPFR_SET_NAN (r);
              MPFR_RET_NAN;
            }
          else /* log10(+Inf) = +Inf */
            {
              MPFR_SET_INF (r);
              MPFR_SET_POS (r);
              MPFR_RET (0); /* exact */
            }
        }
      else /* a = 0 */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (a));
          MPFR_SET_INF (r);
          MPFR_SET_NEG (r);
          MPFR_RET (0); /* log10(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_cmp_ui (a, 1) == 0)
    {
      MPFR_SET_ZERO (r);
      MPFR_SET_POS (r);
      MPFR_RET (0); /* result is exact */
    }

  MPFR_SAVE_EXPO_MARK (expo);

  /* General case */
  {
    /* Declaration of the intermediary variable */
    mpfr_t t, tt;
    MPFR_ZIV_DECL (loop);
    /* Declaration of the size variable */
    mpfr_prec_t Ny = MPFR_PREC(r);   /* Precision of output variable */
    mpfr_prec_t Nt;        /* Precision of the intermediary variable */
    mpfr_exp_t  err;                           /* Precision of error */

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

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

    /* First computation of log10 */
    MPFR_ZIV_INIT (loop, Nt);
    for (;;)
      {
        /* compute log10 */
        mpfr_set_ui (t, 10, MPFR_RNDN);   /* 10 */
        mpfr_log (t, t, MPFR_RNDD);       /* log(10) */
        mpfr_log (tt, a, MPFR_RNDN);      /* log(a) */
        mpfr_div (t, tt, t, MPFR_RNDN);   /* log(a)/log(10) */

        /* estimation of the error */
        err = Nt - 4;
        if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
          break;

        /* log10(10^n) is exact:
           FIXME: Can we have 10^n exactly representable as a mpfr_t
           but n can't fit an unsigned long? */
        if (MPFR_IS_POS (t)
            && mpfr_integer_p (t) && mpfr_fits_ulong_p (t, MPFR_RNDN)
            && !mpfr_ui_pow_ui (tt, 10, mpfr_get_ui (t, MPFR_RNDN), MPFR_RNDN)
            && mpfr_cmp (a, tt) == 0)
          break;

        /* actualisation of the precision */
        MPFR_ZIV_NEXT (loop, Nt);
//.........这里部分代码省略.........