C++ MPFR_ZIV_NEXT函数代码示例

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 */
}

//.........这里部分代码省略.........
      MPFR_ASSERTD (inex == 0);
      goto clear_and_return;
    }

  MPFR_GROUP_INIT_4 (group, MPFR_PREC_MIN, b, c, z_pre, f);

  MPFR_ZIV_INIT (loop, d);
  for (;;)
    {
      /* Principal loop: we compute, in z_pre,
         an approximation of Zeta(s), that we send to can_round */
      if (MPFR_GET_EXP (s1) <= -(mp_exp_t) ((mpfr_prec_t) (d-3)/2))
        /* Branch 1: when s-1 is very small, one
           uses the approximation Zeta(s)=1/(s-1)+gamma,
           where gamma is Euler's constant */
        {
          dint = MAX (d + 3, precs);
          MPFR_TRACE (printf ("branch 1\ninternal precision=%d\n", dint));
          MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
          mpfr_div (z_pre, __gmpfr_one, s1, GMP_RNDN);
          mpfr_const_euler (f, GMP_RNDN);
          mpfr_add (z_pre, z_pre, f, GMP_RNDN);
        }
      else /* Branch 2 */
        {
          size_t size;

          MPFR_TRACE (printf ("branch 2\n"));
          /* Computation of parameters n, p and working precision */
          dnep = (double) d * LOG2;
          sd = mpfr_get_d (s, GMP_RNDN);
          /* beta = dnep + 0.61 + sd * log (6.2832 / sd);
             but a larger value is ok */
#define LOG6dot2832 1.83787940484160805532
          beta = dnep + 0.61 + sd * (LOG6dot2832 - LOG2 *
                                     __gmpfr_floor_log2 (sd));
          if (beta <= 0.0)
            {
              p = 0;
              /* n = 1 + (int) (exp ((dnep - LOG2) / sd)); */
              n = 1 + (int) __gmpfr_ceil_exp2 ((d - 1.0) / sd);
            }
          else
            {
              p = 1 + (int) beta / 2;
              n = 1 + (int) ((sd + 2.0 * (double) p - 1.0) / 6.2832);
            }
          MPFR_TRACE (printf ("\nn=%d\np=%d\n",n,p));
          /* add = 4 + floor(1.5 * log(d) / log (2)).
             We should have add >= 10, which is always fulfilled since
             d = precz + 11 >= 12, thus ceil(log2(d)) >= 4 */
          add = 4 + (3 * MPFR_INT_CEIL_LOG2 (d)) / 2;
          MPFR_ASSERTD(add >= 10);
          dint = d + add;
          if (dint < precs)
            dint = precs;

          MPFR_TRACE (printf("internal precision=%d\n",dint));

          size = (p + 1) * sizeof(mpfr_t);
          tc1 = (mpfr_t*) (*__gmp_allocate_func) (size);
          for (l=1; l<=p; l++)
            mpfr_init2 (tc1[l], dint);
          MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);

          MPFR_TRACE (printf ("precision of z =%d\n", precz));

          /* Computation of the coefficients c_k */
          mpfr_zeta_c (p, tc1);
          /* Computation of the 3 parts of the fonction Zeta. */
          mpfr_zeta_part_a (z_pre, s, n);
          mpfr_zeta_part_b (b, s, n, p, tc1);
          /* s1 = s-1 is already computed above */
          mpfr_div (c, __gmpfr_one, s1, GMP_RNDN);
          mpfr_ui_pow (f, n, s1, GMP_RNDN);
          mpfr_div (c, c, f, GMP_RNDN);
          MPFR_TRACE (MPFR_DUMP (c));
          mpfr_add (z_pre, z_pre, c, GMP_RNDN);
          mpfr_add (z_pre, z_pre, b, GMP_RNDN);
          for (l=1; l<=p; l++)
            mpfr_clear (tc1[l]);
          (*__gmp_free_func) (tc1, size);
          /* End branch 2 */
        }

      MPFR_TRACE (MPFR_DUMP (z_pre));
      if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, d-3, precz, rnd_mode)))
        break;
      MPFR_ZIV_NEXT (loop, d);
    }
  MPFR_ZIV_FREE (loop);

  inex = mpfr_set (z, z_pre, rnd_mode);

  MPFR_GROUP_CLEAR (group);
 clear_and_return:
  mpfr_clear (s1);

  return inex;
}

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

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

/* 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;
}

//.........这里部分代码省略.........
        {
          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;
      /* check for huge cancellation */
      if (err < (mp_exp_t) MPFR_PREC (y))
        m += MPFR_PREC (y) - err;
      /* Check if near 1 */
      if (MPFR_GET_EXP (c) == 1
          && MPFR_MANT (c)[MPFR_LIMB_SIZE (c)-1] == MPFR_LIMB_HIGHBIT)
        m += m;

    next_step:
      MPFR_ZIV_NEXT (loop, m);
      mpfr_set_prec (c, m);
    }
  MPFR_ZIV_FREE (loop);

  mpfr_set (y, c, rnd_mode);

  mpfr_clear (c);
  mpfr_clear (xr);

  MPFR_RET (1); /* Always inexact */
}

//.........这里部分代码省略.........
            {
              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;
          }
        mpfr_log (t, t, GMP_RNDN);        /* log(1+x) */

        /* the error is bounded by (1/2+2^(1-EXP(t))*ulp(t) (cf algorithms.tex)
           if EXP(t)>=2, then error <= ulp(t)
           if EXP(t)<=1, then error <= 2^(2-EXP(t))*ulp(t) */
        err = Nt - MAX (0, 2 - MPFR_GET_EXP (t));

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

        /* increase the precision */
        MPFR_ZIV_NEXT (loop, Nt);
        mpfr_set_prec (t, Nt);
      }
    inexact = mpfr_set (y, t, rnd_mode);

  end:
    MPFR_ZIV_FREE (loop);
    mpfr_clear (t);
  }

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

//.........这里部分代码省略.........
        {
          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 */
      d = MAX(3, d + 1);
      err = Nt - (d + 1);

      if (MPFR_LIKELY ((d <= Nt / 2) && MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
        {
          inexact = mpfr_set4 (y, t, rnd_mode, sign);
          break;
        }

      /* if t=1, we still can round since |sinh(x)| < 1 */
      if (MPFR_GET_EXP (t) == 1)
        goto set_one;

      /* Actualisation of the precision */
      MPFR_ZIV_NEXT (loop, Nt);
      MPFR_GROUP_REPREC_2 (group, Nt, t, te);
    }
    MPFR_ZIV_FREE (loop);
    MPFR_GROUP_CLEAR (group);
  }
  MPFR_SAVE_EXPO_FREE (expo);
  inexact = mpfr_check_range (y, inexact, rnd_mode);

  return inexact;
}

//.........这里部分代码省略.........
        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));
                if (rnd == MPFR_RNDN)
                  {
                    mpfr_t y2, nn;

                    /* We cannot decide now whether the result should be
                       rounded toward zero or away from zero. So, like
                       in mpfr_pow_pos_z, let's use the general case of
                       mpfr_pow in precision 2. */
                    MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_SIGN (x),
                                                    MPFR_EXP (x) - 1) != 0);
                    mpfr_init2 (y2, 2);
                    mpfr_init2 (nn, sizeof (long) * CHAR_BIT);
                    inexact = mpfr_set_si (nn, n, MPFR_RNDN);
                    MPFR_ASSERTN (inexact == 0);
                    inexact = mpfr_pow_general (y2, x, nn, rnd, 1,
                                                (mpfr_save_expo_t *) NULL);
                    mpfr_clear (nn);
                    mpfr_set (y, y2, MPFR_RNDN);
                    mpfr_clear (y2);
                    MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_UNDERFLOW);
                    goto end;
                  }
                else
                  {
                    MPFR_SAVE_EXPO_FREE (expo);
                    return mpfr_underflow (y, rnd, abs_n & 1 ?
                                           MPFR_SIGN (x) : MPFR_SIGN_POS);
                  }
              }
            /* error estimate -- see pow function in algorithms.ps */
            if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - size_n - 2, Ny, rnd)))
              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);
        mpfr_clear (t);

      end:
        MPFR_SAVE_EXPO_FREE (expo);
        return mpfr_check_range (y, inexact, rnd);
      }
    }
}

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

  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);
      mpfr_sqrt (c, c, MPFR_RNDZ);
      if (MPFR_IS_NEG_SIGN(sign))
        MPFR_CHANGE_SIGN(c);

      /* Warning: c may be 0! */
      if (MPFR_UNLIKELY (MPFR_IS_ZERO (c)))
        {
          /* Huge cancellation: increase prec a lot! */
          m = MAX (m, MPFR_PREC (x));
          m = 2 * m;
        }
      else
        {
          /* the absolute error on c is at most 2^(3-m-EXP(c)),
             plus 2^(2-m) if there was an argument reduction.
             Since EXP(c) <= 1, 3-m-EXP(c) >= 2-m, thus the error
             is at most 2^(3-m-EXP(c)) in case of argument reduction. */
          err = 2 * MPFR_GET_EXP (c) + (mpfr_exp_t) m - 3 - (reduce != 0);
          if (MPFR_CAN_ROUND (c, err, precy, rnd_mode))
            break;

          /* check for huge cancellation (Near 0) */
          if (err < (mpfr_exp_t) MPFR_PREC (y))
            m += MPFR_PREC (y) - err;
          /* Check if near 1 */
          if (MPFR_GET_EXP (c) == 1)
            m += m;
        }

    ziv_next:
      /* Else generic increase */
      MPFR_ZIV_NEXT (loop, m);
    }
  MPFR_ZIV_FREE (loop);

  inexact = mpfr_set (y, c, rnd_mode);
  /* inexact cannot be 0, since this would mean that c was representable
     within the target precision, but in that case mpfr_can_round will fail */

  mpfr_clear (c);
  mpfr_clear (xr);

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

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

          /* we have 7 ulps of error from the above roundings,
             4 ulps from the 4/s^2 second order term,
             plus the canceled bits */
          if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp1, p-cancel-4, q, rnd_mode)))
            break;

          /* VL: I think it is better to have an increment that it isn't
             too low; in particular, the increment must be positive even
             if cancel = 0 (can this occur?). */
          p += cancel >= 8 ? cancel : 8;
        }
      else
        {
          /* TODO: find why this case can occur and what is best to do
             with it. */
          p += 32;
        }

      MPFR_ZIV_NEXT (loop, p);
      MPFR_GROUP_REPREC_2 (group, p, tmp1, tmp2);
    }
  MPFR_ZIV_FREE (loop);
  inexact = mpfr_set (r, tmp1, rnd_mode);
  /* We clean */
  MPFR_GROUP_CLEAR (group);

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

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

/* 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);
}

//.........这里部分代码省略.........
      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 */
      precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s));
      sd = mpfr_get_d (s, GMP_RNDN) - 1.0;
      if (sd < 0.0)
        sd = -sd; /* now sd = abs(s-1.0) */
      /* Precision prec1 is the precision on elementary computations;
         it ensures a final precision prec1 - add for zeta(s) */
      /* eps = pow (2.0, - (double) precz - 14.0); */
      eps = __gmpfr_ceil_exp2 (- (double) precz - 14.0);
      m1 = 1.0 + MAX(1.0 / eps,  2.0 * sd) * (1.0 + eps);
      c = (1.0 + eps) * (1.0 + eps * MAX(8.0, m1));
      /* add = 1 + floor(log(c*c*c*(13 + m1))/log(2)); */
      add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1));
      prec1 = precz + add;
      prec1 = MAX (prec1, precs1) + 10;

      MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p);
      MPFR_ZIV_INIT (loop, prec1);
      for (;;)
        {
          mpfr_sub (s1, __gmpfr_one, s, GMP_RNDN);/* s1 = 1-s */
          mpfr_zeta_pos (z_pre, s1, GMP_RNDN);   /* zeta(1-s)  */
          mpfr_gamma (y, s1, GMP_RNDN);          /* gamma(1-s) */
          if (MPFR_IS_INF (y)) /* Zeta(s) < 0 for -4k-2 < s < -4k,
                                  Zeta(s) > 0 for -4k < s < -4k+2 */
            {
              MPFR_SET_INF (z_pre);
              mpfr_div_2ui (s1, s, 2, GMP_RNDN); /* s/4, exact */
              mpfr_frac (s1, s1, GMP_RNDN); /* exact, -1 < s1 < 0 */
              if (mpfr_cmp_si_2exp (s1, -1, -1) > 0)
                MPFR_SET_NEG (z_pre);
              else
                MPFR_SET_POS (z_pre);
              break;
            }
          mpfr_mul (z_pre, z_pre, y, GMP_RNDN);  /* gamma(1-s)*zeta(1-s) */
          mpfr_const_pi (p, GMP_RNDD);
          mpfr_mul (y, s, p, GMP_RNDN);
          mpfr_div_2ui (y, y, 1, GMP_RNDN);      /* s*Pi/2 */
          mpfr_sin (y, y, GMP_RNDN);             /* sin(Pi*s/2) */
          mpfr_mul (z_pre, z_pre, y, GMP_RNDN);
          mpfr_mul_2ui (y, p, 1, GMP_RNDN);      /* 2*Pi */
          mpfr_neg (s1, s1, GMP_RNDN);           /* s-1 */
          mpfr_pow (y, y, s1, GMP_RNDN);         /* (2*Pi)^(s-1) */
          mpfr_mul (z_pre, z_pre, y, GMP_RNDN);
          mpfr_mul_2ui (z_pre, z_pre, 1, GMP_RNDN);

          if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz,
                                           rnd_mode)))
            break;

          MPFR_ZIV_NEXT (loop, prec1);
          MPFR_GROUP_REPREC_4 (group, prec1, z_pre, s1, y, p);
        }
      MPFR_ZIV_FREE (loop);
      inex = mpfr_set (z, z_pre, rnd_mode);
      MPFR_GROUP_CLEAR (group);
    }

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (z, inex, rnd_mode);
}