C++ MPFR_ZIV_DECL函数代码示例

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

int
mpfr_erfc (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
{
  int inex;
  mpfr_t tmp;
  mpfr_exp_t te, err;
  mpfr_prec_t prec;
  mpfr_exp_t emin = mpfr_get_emin ();
  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),
     ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inex));

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN (x))
        {
          MPFR_SET_NAN (y);
          MPFR_RET_NAN;
        }
      /* erfc(+inf) = 0+, erfc(-inf) = 2 erfc (0) = 1 */
      else if (MPFR_IS_INF (x))
        return mpfr_set_ui (y, MPFR_IS_POS (x) ? 0 : 2, rnd);
      else
        return mpfr_set_ui (y, 1, rnd);
    }

  if (MPFR_SIGN (x) > 0)
    {
      /* by default, emin = 1-2^30, thus the smallest representable
         number is 1/2*2^emin = 2^(-2^30):
         for x >= 27282, erfc(x) < 2^(-2^30-1), and
         for x >= 1787897414, erfc(x) < 2^(-2^62-1).
      */
      if ((emin >= -1073741823 && mpfr_cmp_ui (x, 27282) >= 0) ||
          mpfr_cmp_ui (x, 1787897414) >= 0)
        {
          /* May be incorrect if MPFR_EMAX_MAX >= 2^62. */
          MPFR_ASSERTN ((MPFR_EMAX_MAX >> 31) >> 31 == 0);
          return mpfr_underflow (y, (rnd == MPFR_RNDN) ? MPFR_RNDZ : rnd, 1);
        }
    }

int
mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode)
{
  mpfr_t x;
  int inexact;

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

  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))
        {
          MPFR_SET_INF (y);
          MPFR_SET_SAME_SIGN (y, xt);
          MPFR_RET (0);
        }
      else /* xt is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (xt));
          MPFR_SET_ZERO (y);   /* sinh(0) = 0 */
          MPFR_SET_SAME_SIGN (y, xt);
          MPFR_RET (0);
        }
    }

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

  MPFR_TMP_INIT_ABS (x, xt);

  {
    mpfr_t t, ti;
    mp_exp_t d;
    mp_prec_t Nt;    /* Precision of the intermediary variable */
    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 */
    Nt = MAX (MPFR_PREC (x), MPFR_PREC (y));
    /* the optimal number of bits : see algorithms.ps */
    Nt = Nt + MPFR_INT_CEIL_LOG2 (Nt) + 4;
    /* If x is near 0, exp(x) - 1/exp(x) = 2*x+x^3/3+O(x^5) */
    if (MPFR_GET_EXP (x) < 0)
      Nt -= 2*MPFR_GET_EXP (x);

    /* initialise of intermediary variables */
    MPFR_GROUP_INIT_2 (group, Nt, t, ti);

    /* First computation of sinh */
    MPFR_ZIV_INIT (loop, Nt);
    for (;;) {
      /* compute sinh */
      mpfr_clear_flags ();
      mpfr_exp (t, x, GMP_RNDD);        /* exp(x) */
      /* exp(x) can overflow! */
      /* BUG/TODO/FIXME: exp can overflow but sinh may be representable! */
      if (MPFR_UNLIKELY (mpfr_overflow_p ())) {
        inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN (xt));
        MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
        break;
      }
      d = MPFR_GET_EXP (t);
      mpfr_ui_div (ti, 1, t, GMP_RNDU); /* 1/exp(x) */
      mpfr_sub (t, t, ti, GMP_RNDN);    /* exp(x) - 1/exp(x) */
      mpfr_div_2ui (t, t, 1, GMP_RNDN);  /* 1/2(exp(x) - 1/exp(x)) */

      /* it may be that t is zero (in fact, it can only occur when te=1,
         and thus ti=1 too) */
      if (MPFR_IS_ZERO (t))
        err = Nt; /* double the precision */
      else
        {
          /* calculation of the error */
          d = d - MPFR_GET_EXP (t) + 2;
          /* error estimate: err = Nt-(__gmpfr_ceil_log2(1+pow(2,d)));*/
          err = Nt - (MAX (d, 0) + 1);
          if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, MPFR_PREC (y), rnd_mode)))
            {
              inexact = mpfr_set4 (y, t, rnd_mode, MPFR_SIGN (xt));
              break;
            }
        }
      /* actualisation of the precision */
      Nt += err;
      MPFR_ZIV_NEXT (loop, Nt);
      MPFR_GROUP_REPREC_2 (group, Nt, t, ti);
    }
    MPFR_ZIV_FREE (loop);
    MPFR_GROUP_CLEAR (group);
//.........这里部分代码省略.........

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

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

/* agm(x,y) is between x and y, so we don't need to save exponent range */
int
mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode)
{
  int compare, inexact;
  mp_size_t s;
  mp_prec_t p, q;
  mp_limb_t *up, *vp, *tmpp;
  mpfr_t u, v, tmp;
  unsigned long n; /* number of iterations */
  unsigned long err = 0;
  MPFR_ZIV_DECL (loop);
  MPFR_TMP_DECL(marker);

  MPFR_LOG_FUNC (("op2[%#R]=%R op1[%#R]=%R rnd=%d", op2,op2,op1,op1,rnd_mode),
                 ("r[%#R]=%R inexact=%d", r, r, inexact));

  /* Deal with special values */
  if (MPFR_ARE_SINGULAR (op1, op2))
    {
      /* If a or b is NaN, the result is NaN */
      if (MPFR_IS_NAN(op1) || MPFR_IS_NAN(op2))
        {
          MPFR_SET_NAN(r);
          MPFR_RET_NAN;
        }
      /* now one of a or b is Inf or 0 */
      /* If a and b is +Inf, the result is +Inf.
         Otherwise if a or b is -Inf or 0, the result is NaN */
      else if (MPFR_IS_INF(op1) || MPFR_IS_INF(op2))
        {
          if (MPFR_IS_STRICTPOS(op1) && MPFR_IS_STRICTPOS(op2))
            {
              MPFR_SET_INF(r);
              MPFR_SET_SAME_SIGN(r, op1);
              MPFR_RET(0); /* exact */
            }
          else
            {
              MPFR_SET_NAN(r);
              MPFR_RET_NAN;
            }
        }
      else /* a and b are neither NaN nor Inf, and one is zero */
        {  /* If a or b is 0, the result is +0 since a sqrt is positive */
          MPFR_ASSERTD (MPFR_IS_ZERO (op1) || MPFR_IS_ZERO (op2));
          MPFR_SET_POS (r);
          MPFR_SET_ZERO (r);
          MPFR_RET (0); /* exact */
        }
    }
  MPFR_CLEAR_FLAGS (r);

  /* If a or b is negative (excluding -Infinity), the result is NaN */
  if (MPFR_UNLIKELY(MPFR_IS_NEG(op1) || MPFR_IS_NEG(op2)))
    {
      MPFR_SET_NAN(r);
      MPFR_RET_NAN;
    }

  /* Precision of the following calculus */
  q = MPFR_PREC(r);
  p = q + MPFR_INT_CEIL_LOG2(q) + 15;
  MPFR_ASSERTD (p >= 7); /* see algorithms.tex */
  s = (p - 1) / BITS_PER_MP_LIMB + 1;

  /* b (op2) and a (op1) are the 2 operands but we want b >= a */
  compare = mpfr_cmp (op1, op2);
  if (MPFR_UNLIKELY( compare == 0 ))
    {
      mpfr_set (r, op1, rnd_mode);
      MPFR_RET (0); /* exact */
    }
  else if (compare > 0)
    {
      mpfr_srcptr t = op1;
      op1 = op2;
      op2 = t;
    }
  /* Now b(=op2) >= a (=op1) */

  MPFR_TMP_MARK(marker);

  /* Main loop */
  MPFR_ZIV_INIT (loop, p);
  for (;;)
    {
      mp_prec_t eq;

      /* Init temporary vars */
      MPFR_TMP_INIT (up, u, p, s);
      MPFR_TMP_INIT (vp, v, p, s);
      MPFR_TMP_INIT (tmpp, tmp, p, s);

      /* Calculus of un and vn */
      mpfr_mul (u, op1, op2, GMP_RNDN); /* Faster since PREC(op) < PREC(u) */
      mpfr_sqrt (u, u, GMP_RNDN);
      mpfr_add (v, op1, op2, GMP_RNDN); /* add with !=prec is still good*/
      mpfr_div_2ui (v, v, 1, GMP_RNDN);
      n = 1;
//.........这里部分代码省略.........

/* We use the reflection formula
  Gamma(1+t) Gamma(1-t) = - Pi t / sin(Pi (1 + t))
  in order to treat the case x <= 1,
  i.e. with x = 1-t, then Gamma(x) = -Pi*(1-x)/sin(Pi*(2-x))/GAMMA(2-x)
*/
int
mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t xp, GammaTrial, tmp, tmp2;
  mpz_t fact;
  mpfr_prec_t realprec;
  int compared, is_integer;
  int inex = 0;  /* 0 means: result gamma not set yet */
  MPFR_GROUP_DECL (group);
  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),
     ("gamma[%Pu]=%.*Rg inexact=%d",
      mpfr_get_prec (gamma), mpfr_log_prec, gamma, inex));

  /* Trivial cases */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN (x))
        {
          MPFR_SET_NAN (gamma);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF (x))
        {
          if (MPFR_IS_NEG (x))
            {
              MPFR_SET_NAN (gamma);
              MPFR_RET_NAN;
            }
          else
            {
              MPFR_SET_INF (gamma);
              MPFR_SET_POS (gamma);
              MPFR_RET (0);  /* exact */
            }
        }
      else /* x is zero */
        {
          MPFR_ASSERTD(MPFR_IS_ZERO(x));
          MPFR_SET_INF(gamma);
          MPFR_SET_SAME_SIGN(gamma, x);
          MPFR_SET_DIVBY0 ();
          MPFR_RET (0);  /* exact */
        }
    }

  /* Check for tiny arguments, where gamma(x) ~ 1/x - euler + ....
     We know from "Bound on Runs of Zeros and Ones for Algebraic Functions",
     Proceedings of Arith15, T. Lang and J.-M. Muller, 2001, that the maximal
     number of consecutive zeroes or ones after the round bit is n-1 for an
     input of n bits. But we need a more precise lower bound. Assume x has
     n bits, and 1/x is near a floating-point number y of n+1 bits. We can
     write x = X*2^e, y = Y/2^f with X, Y integers of n and n+1 bits.
     Thus X*Y^2^(e-f) is near from 1, i.e., X*Y is near from 2^(f-e).
     Two cases can happen:
     (i) either X*Y is exactly 2^(f-e), but this can happen only if X and Y
         are themselves powers of two, i.e., x is a power of two;
     (ii) or X*Y is at distance at least one from 2^(f-e), thus
          |xy-1| >= 2^(e-f), or |y-1/x| >= 2^(e-f)/x = 2^(-f)/X >= 2^(-f-n).
          Since ufp(y) = 2^(n-f) [ufp = unit in first place], this means
          that the distance |y-1/x| >= 2^(-2n) ufp(y).
          Now assuming |gamma(x)-1/x| <= 1, which is true for x <= 1,
          if 2^(-2n) ufp(y) >= 2, the error is at most 2^(-2n-1) ufp(y),
          and round(1/x) with precision >= 2n+2 gives the correct result.
          If x < 2^E, then y > 2^(-E), thus ufp(y) > 2^(-E-1).
          A sufficient condition is thus EXP(x) + 2 <= -2 MAX(PREC(x),PREC(Y)).
  */
  if (MPFR_GET_EXP (x) + 2
      <= -2 * (mpfr_exp_t) MAX(MPFR_PREC(x), MPFR_PREC(gamma)))
    {
      int sign = MPFR_SIGN (x); /* retrieve sign before possible override */
      int special;
      MPFR_BLOCK_DECL (flags);

      MPFR_SAVE_EXPO_MARK (expo);

      /* for overflow cases, see below; this needs to be done
         before x possibly gets overridden. */
      special =
        MPFR_GET_EXP (x) == 1 - MPFR_EMAX_MAX &&
        MPFR_IS_POS_SIGN (sign) &&
        MPFR_IS_LIKE_RNDD (rnd_mode, sign) &&
        mpfr_powerof2_raw (x);

      MPFR_BLOCK (flags, inex = mpfr_ui_div (gamma, 1, x, rnd_mode));
      if (inex == 0) /* x is a power of two */
        {
          /* return RND(1/x - euler) = RND(+/- 2^k - eps) with eps > 0 */
          if (rnd_mode == MPFR_RNDN || MPFR_IS_LIKE_RNDU (rnd_mode, sign))
            inex = 1;
          else
            {
//.........这里部分代码省略.........

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

/* Don't need to save/restore exponent range: the cache does it */
int
mpfr_const_pi_internal (mpfr_ptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t a, A, B, D, S;
  mpfr_prec_t px, p, cancel, k, kmax;
  MPFR_ZIV_DECL (loop);
  int inex;

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

  px = MPFR_PREC (x);

  /* we need 9*2^kmax - 4 >= px+2*kmax+8 */
  for (kmax = 2; ((px + 2 * kmax + 12) / 9) >> kmax; kmax ++);

  p = px + 3 * kmax + 14; /* guarantees no recomputation for px <= 10000 */

  mpfr_init2 (a, p);
  mpfr_init2 (A, p);
  mpfr_init2 (B, p);
  mpfr_init2 (D, p);
  mpfr_init2 (S, p);

  MPFR_ZIV_INIT (loop, p);
  for (;;) {
    mpfr_set_ui (a, 1, MPFR_RNDN);          /* a = 1 */
    mpfr_set_ui (A, 1, MPFR_RNDN);          /* A = a^2 = 1 */
    mpfr_set_ui_2exp (B, 1, -1, MPFR_RNDN); /* B = b^2 = 1/2 */
    mpfr_set_ui_2exp (D, 1, -2, MPFR_RNDN); /* D = 1/4 */

#define b B
#define ap a
#define Ap A
#define Bp B
    for (k = 0; ; k++)
      {
        /* invariant: 1/2 <= B <= A <= a < 1 */
        mpfr_add (S, A, B, MPFR_RNDN); /* 1 <= S <= 2 */
        mpfr_div_2ui (S, S, 2, MPFR_RNDN); /* exact, 1/4 <= S <= 1/2 */
        mpfr_sqrt (b, B, MPFR_RNDN); /* 1/2 <= b <= 1 */
        mpfr_add (ap, a, b, MPFR_RNDN); /* 1 <= ap <= 2 */
        mpfr_div_2ui (ap, ap, 1, MPFR_RNDN); /* exact, 1/2 <= ap <= 1 */
        mpfr_mul (Ap, ap, ap, MPFR_RNDN); /* 1/4 <= Ap <= 1 */
        mpfr_sub (Bp, Ap, S, MPFR_RNDN); /* -1/4 <= Bp <= 3/4 */
        mpfr_mul_2ui (Bp, Bp, 1, MPFR_RNDN); /* -1/2 <= Bp <= 3/2 */
        mpfr_sub (S, Ap, Bp, MPFR_RNDN);
        MPFR_ASSERTN (mpfr_cmp_ui (S, 1) < 0);
        cancel = mpfr_cmp_ui (S, 0) ? (mpfr_uexp_t) -mpfr_get_exp(S) : p;
        /* MPFR_ASSERTN (cancel >= px || cancel >= 9 * (1 << k) - 4); */
        mpfr_mul_2ui (S, S, k, MPFR_RNDN);
        mpfr_sub (D, D, S, MPFR_RNDN);
        /* stop when |A_k - B_k| <= 2^(k-p) i.e. cancel >= p-k */
        if (cancel + k >= p)
          break;
      }
#undef b
#undef ap
#undef Ap
#undef Bp

      mpfr_div (A, B, D, MPFR_RNDN);

      /* MPFR_ASSERTN(p >= 2 * k + 8); */
      if (MPFR_LIKELY (MPFR_CAN_ROUND (A, p - 2 * k - 8, px, rnd_mode)))
        break;

      p += kmax;
      MPFR_ZIV_NEXT (loop, p);
      mpfr_set_prec (a, p);
      mpfr_set_prec (A, p);
      mpfr_set_prec (B, p);
      mpfr_set_prec (D, p);
      mpfr_set_prec (S, p);
  }
  MPFR_ZIV_FREE (loop);
  inex = mpfr_set (x, A, rnd_mode);

  mpfr_clear (a);
  mpfr_clear (A);
  mpfr_clear (B);
  mpfr_clear (D);
  mpfr_clear (S);

  return inex;
}

int
mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mpfr_rnd_t rnd_mode)
{
  mpfr_t x;
  int inexact;

  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));

  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))
        {
          MPFR_SET_INF (y);
          MPFR_SET_SAME_SIGN (y, xt);
          MPFR_RET (0);
        }
      else /* xt is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (xt));
          MPFR_SET_ZERO (y);   /* sinh(0) = 0 */
          MPFR_SET_SAME_SIGN (y, xt);
          MPFR_RET (0);
        }
    }

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

  MPFR_TMP_INIT_ABS (x, xt);

  {
    mpfr_t t, ti;
    mpfr_exp_t d;
    mpfr_prec_t Nt;    /* Precision of the intermediary variable */
    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 */
    Nt = MAX (MPFR_PREC (x), MPFR_PREC (y));
    /* the optimal number of bits : see algorithms.ps */
    Nt = Nt + MPFR_INT_CEIL_LOG2 (Nt) + 4;
    /* If x is near 0, exp(x) - 1/exp(x) = 2*x+x^3/3+O(x^5) */
    if (MPFR_GET_EXP (x) < 0)
      Nt -= 2*MPFR_GET_EXP (x);

    /* initialise of intermediary variables */
    MPFR_GROUP_INIT_2 (group, Nt, t, ti);

    /* First computation of sinh */
    MPFR_ZIV_INIT (loop, Nt);
    for (;;)
      {
        MPFR_BLOCK_DECL (flags);

        /* compute sinh */
        MPFR_BLOCK (flags, mpfr_exp (t, x, MPFR_RNDD));
        if (MPFR_OVERFLOW (flags))
          /* exp(x) does overflow */
          {
            /* sinh(x) = 2 * sinh(x/2) * cosh(x/2) */
            mpfr_div_2ui (ti, x, 1, MPFR_RNDD); /* exact */

            /* t <- cosh(x/2): error(t) <= 1 ulp(t) */
            MPFR_BLOCK (flags, mpfr_cosh (t, ti, MPFR_RNDD));
            if (MPFR_OVERFLOW (flags))
              /* when x>1 we have |sinh(x)| >= cosh(x/2), so sinh(x)
                 overflows too */
              {
                inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN (xt));
                MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
                break;
              }

            /* ti <- sinh(x/2): , error(ti) <= 1 ulp(ti)
               cannot overflow because 0 < sinh(x) < cosh(x) when x > 0 */
            mpfr_sinh (ti, ti, MPFR_RNDD);

            /* multiplication below, error(t) <= 5 ulp(t) */
            MPFR_BLOCK (flags, mpfr_mul (t, t, ti, MPFR_RNDD));
            if (MPFR_OVERFLOW (flags))
              {
                inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN (xt));
                MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
                break;
              }

            /* doubling below, exact */
//.........这里部分代码省略.........

/* Implements asymptotic expansion for jn or yn (formulae 9.2.5 and 9.2.6
   from Abramowitz & Stegun).
   Assumes |z| > p log(2)/2, where p is the target precision
   (z can be negative only for jn).
   Return 0 if the expansion does not converge enough (the value 0 as inexact
   flag should not happen for normal input).
*/
static int
FUNCTION (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r)
{
  mpfr_t s, c, P, Q, t, iz, err_t, err_s, err_u;
  mpfr_prec_t w;
  long k;
  int inex, stop, diverge = 0;
  mpfr_exp_t err2, err;
  MPFR_ZIV_DECL (loop);

  mpfr_init (c);

  w = MPFR_PREC(res) + MPFR_INT_CEIL_LOG2(MPFR_PREC(res)) + 4;

  MPFR_ZIV_INIT (loop, w);
  for (;;)
    {
      mpfr_set_prec (c, w);
      mpfr_init2 (s, w);
      mpfr_init2 (P, w);
      mpfr_init2 (Q, w);
      mpfr_init2 (t, w);
      mpfr_init2 (iz, w);
      mpfr_init2 (err_t, 31);
      mpfr_init2 (err_s, 31);
      mpfr_init2 (err_u, 31);

      /* Approximate sin(z) and cos(z). In the following, err <= k means that
         the approximate value y and the true value x are related by
         y = x * (1 + u)^k with |u| <= 2^(-w), following Higham's method. */
      mpfr_sin_cos (s, c, z, MPFR_RNDN);
      if (MPFR_IS_NEG(z))
        mpfr_neg (s, s, MPFR_RNDN); /* compute jn/yn(|z|), fix sign later */
      /* The absolute error on s/c is bounded by 1/2 ulp(1/2) <= 2^(-w-1). */
      mpfr_add (t, s, c, MPFR_RNDN);
      mpfr_sub (c, s, c, MPFR_RNDN);
      mpfr_swap (s, t);
      /* now s approximates sin(z)+cos(z), and c approximates sin(z)-cos(z),
         with total absolute error bounded by 2^(1-w). */

      /* precompute 1/(8|z|) */
      mpfr_si_div (iz, MPFR_IS_POS(z) ? 1 : -1, z, MPFR_RNDN);   /* err <= 1 */
      mpfr_div_2ui (iz, iz, 3, MPFR_RNDN);

      /* compute P and Q */
      mpfr_set_ui (P, 1, MPFR_RNDN);
      mpfr_set_ui (Q, 0, MPFR_RNDN);
      mpfr_set_ui (t, 1, MPFR_RNDN); /* current term */
      mpfr_set_ui (err_t, 0, MPFR_RNDN); /* error on t */
      mpfr_set_ui (err_s, 0, MPFR_RNDN); /* error on P and Q (sum of errors) */
      for (k = 1, stop = 0; stop < 4; k++)
        {
          /* compute next term: t(k)/t(k-1) = (2n+2k-1)(2n-2k+1)/(8kz) */
          mpfr_mul_si (t, t, 2 * (n + k) - 1, MPFR_RNDN); /* err <= err_k + 1 */
          mpfr_mul_si (t, t, 2 * (n - k) + 1, MPFR_RNDN); /* err <= err_k + 2 */
          mpfr_div_ui (t, t, k, MPFR_RNDN);               /* err <= err_k + 3 */
          mpfr_mul (t, t, iz, MPFR_RNDN);                 /* err <= err_k + 5 */
          /* the relative error on t is bounded by (1+u)^(5k)-1, which is
             bounded by 6ku for 6ku <= 0.02: first |5 log(1+u)| <= |5.5u|
             for |u| <= 0.15, then |exp(5.5u)-1| <= 6u for |u| <= 0.02. */
          mpfr_mul_ui (err_t, t, 6 * k, MPFR_IS_POS(t) ? MPFR_RNDU : MPFR_RNDD);
          mpfr_abs (err_t, err_t, MPFR_RNDN); /* exact */
          /* the absolute error on t is bounded by err_t * 2^(-w) */
          mpfr_abs (err_u, t, MPFR_RNDU);
          mpfr_mul_2ui (err_u, err_u, w, MPFR_RNDU); /* t * 2^w */
          mpfr_add (err_u, err_u, err_t, MPFR_RNDU); /* max|t| * 2^w */
          if (stop >= 2)
            {
              /* take into account the neglected terms: t * 2^w */
              mpfr_div_2ui (err_s, err_s, w, MPFR_RNDU);
              if (MPFR_IS_POS(t))
                mpfr_add (err_s, err_s, t, MPFR_RNDU);
              else
                mpfr_sub (err_s, err_s, t, MPFR_RNDU);
              mpfr_mul_2ui (err_s, err_s, w, MPFR_RNDU);
              stop ++;
            }
          /* if k is odd, add to Q, otherwise to P */
          else if (k & 1)
            {
              /* if k = 1 mod 4, add, otherwise subtract */
              if ((k & 2) == 0)
                mpfr_add (Q, Q, t, MPFR_RNDN);
              else
                mpfr_sub (Q, Q, t, MPFR_RNDN);
              /* check if the next term is smaller than ulp(Q): if EXP(err_u)
                 <= EXP(Q), since the current term is bounded by
                 err_u * 2^(-w), it is bounded by ulp(Q) */
              if (MPFR_EXP(err_u) <= MPFR_EXP(Q))
                stop ++;
              else
                stop = 0;
            }
//.........这里部分代码省略.........

int
mpfr_acos (mpfr_ptr acos, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t xp, arcc, tmp;
  mpfr_exp_t supplement;
  mpfr_prec_t prec;
  int sign, compared, inexact;
  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),
     ("acos[%Pu]=%.*Rg inexact=%d",
      mpfr_get_prec(acos), mpfr_log_prec, acos, inexact));

  /* Singular cases */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
        {
          MPFR_SET_NAN (acos);
          MPFR_RET_NAN;
        }
      else /* necessarily x=0 */
        {
          MPFR_ASSERTD(MPFR_IS_ZERO(x));
          /* acos(0)=Pi/2 */
          MPFR_SAVE_EXPO_MARK (expo);
          inexact = mpfr_const_pi (acos, rnd_mode);
          mpfr_div_2ui (acos, acos, 1, rnd_mode); /* exact */
          MPFR_SAVE_EXPO_FREE (expo);
          return mpfr_check_range (acos, inexact, rnd_mode);
        }
    }

  /* Set x_p=|x| */
  sign = MPFR_SIGN (x);
  mpfr_init2 (xp, MPFR_PREC (x));
  mpfr_abs (xp, x, MPFR_RNDN); /* Exact */

  compared = mpfr_cmp_ui (xp, 1);

  if (MPFR_UNLIKELY (compared >= 0))
    {
      mpfr_clear (xp);
      if (compared > 0) /* acos(x) = NaN for x > 1 */
        {
          MPFR_SET_NAN(acos);
          MPFR_RET_NAN;
        }
      else
        {
          if (MPFR_IS_POS_SIGN (sign)) /* acos(+1) = +0 */
            return mpfr_set_ui (acos, 0, rnd_mode);
          else /* acos(-1) = Pi */
            return mpfr_const_pi (acos, rnd_mode);
        }
    }

  MPFR_SAVE_EXPO_MARK (expo);

  /* Compute the supplement */
  mpfr_ui_sub (xp, 1, xp, MPFR_RNDD);
  if (MPFR_IS_POS_SIGN (sign))
    supplement = 2 - 2 * MPFR_GET_EXP (xp);
  else
    supplement = 2 - MPFR_GET_EXP (xp);
  mpfr_clear (xp);

  prec = MPFR_PREC (acos);
  prec += MPFR_INT_CEIL_LOG2(prec) + 10 + supplement;

  /* VL: The following change concerning prec comes from r3145
     "Optimize mpfr_acos by choosing a better initial precision."
     but it doesn't seem to be correct and leads to problems (assertion
     failure or very important inefficiency) with tiny arguments.
     Therefore, I've disabled it. */
  /* If x ~ 2^-N, acos(x) ~ PI/2 - x - x^3/6
     If Prec < 2*N, we can't round since x^3/6 won't be counted. */
#if 0
  if (MPFR_PREC (acos) >= MPFR_PREC (x) && MPFR_GET_EXP (x) < 0)
    {
      mpfr_uexp_t pmin = (mpfr_uexp_t) (-2 * MPFR_GET_EXP (x)) + 5;
      MPFR_ASSERTN (pmin <= MPFR_PREC_MAX);
      if (prec < pmin)
        prec = pmin;
    }
#endif

  mpfr_init2 (tmp, prec);
  mpfr_init2 (arcc, prec);

  MPFR_ZIV_INIT (loop, prec);
  for (;;)
    {
      /* acos(x) = Pi/2 - asin(x) = Pi/2 - atan(x/sqrt(1-x^2)) */
      mpfr_sqr (tmp, x, MPFR_RNDN);
      mpfr_ui_sub (tmp, 1, tmp, MPFR_RNDN);
      mpfr_sqrt (tmp, tmp, MPFR_RNDN);
      mpfr_div (tmp, x, tmp, MPFR_RNDN);
//.........这里部分代码省略.........

int
mpfr_exp_3 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
{
  mpfr_t t, x_copy, tmp;
  mpz_t uk;
  mp_exp_t ttt, shift_x;
  unsigned long twopoweri;
  mpz_t *P;
  mp_prec_t *mult;
  int i, k, loop;
  int prec_x;
  mp_prec_t realprec, Prec;
  int iter;
  int inexact = 0;
  MPFR_SAVE_EXPO_DECL (expo);
  MPFR_ZIV_DECL (ziv_loop);

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

  MPFR_SAVE_EXPO_MARK (expo);

  /* decompose x */
  /* we first write x = 1.xxxxxxxxxxxxx
     ----- k bits -- */
  prec_x = MPFR_INT_CEIL_LOG2 (MPFR_PREC (x)) - MPFR_LOG2_BITS_PER_MP_LIMB;
  if (prec_x < 0)
    prec_x = 0;

  ttt = MPFR_GET_EXP (x);
  mpfr_init2 (x_copy, MPFR_PREC(x));
  mpfr_set (x_copy, x, GMP_RNDD);

  /* we shift to get a number less than 1 */
  if (ttt > 0)
    {
      shift_x = ttt;
      mpfr_div_2ui (x_copy, x, ttt, GMP_RNDN);
      ttt = MPFR_GET_EXP (x_copy);
    }
  else
    shift_x = 0;
  MPFR_ASSERTD (ttt <= 0);

  /* Init prec and vars */
  realprec = MPFR_PREC (y) + MPFR_INT_CEIL_LOG2 (prec_x + MPFR_PREC (y));
  Prec = realprec + shift + 2 + shift_x;
  mpfr_init2 (t, Prec);
  mpfr_init2 (tmp, Prec);
  mpz_init (uk);

  /* Main loop */
  MPFR_ZIV_INIT (ziv_loop, realprec);
  for (;;)
    {
      int scaled = 0;
      MPFR_BLOCK_DECL (flags);

      k = MPFR_INT_CEIL_LOG2 (Prec) - MPFR_LOG2_BITS_PER_MP_LIMB;

      /* now we have to extract */
      twopoweri = BITS_PER_MP_LIMB;

      /* Allocate tables */
      P    = (mpz_t*) (*__gmp_allocate_func) (3*(k+2)*sizeof(mpz_t));
      for (i = 0; i < 3*(k+2); i++)
        mpz_init (P[i]);
      mult = (mp_prec_t*) (*__gmp_allocate_func) (2*(k+2)*sizeof(mp_prec_t));

      /* Particular case for i==0 */
      mpfr_extract (uk, x_copy, 0);
      MPFR_ASSERTD (mpz_cmp_ui (uk, 0) != 0);
      mpfr_exp_rational (tmp, uk, shift + twopoweri - ttt, k + 1, P, mult);
      for (loop = 0; loop < shift; loop++)
        mpfr_sqr (tmp, tmp, GMP_RNDD);
      twopoweri *= 2;

      /* General case */
      iter = (k <= prec_x) ? k : prec_x;
      for (i = 1; i <= iter; i++)
        {
          mpfr_extract (uk, x_copy, i);
          if (MPFR_LIKELY (mpz_cmp_ui (uk, 0) != 0))
            {
              mpfr_exp_rational (t, uk, twopoweri - ttt, k  - i + 1, P, mult);
              mpfr_mul (tmp, tmp, t, GMP_RNDD);
            }
          MPFR_ASSERTN (twopoweri <= LONG_MAX/2);
          twopoweri *=2;
        }

      /* Clear tables */
      for (i = 0; i < 3*(k+2); i++)
        mpz_clear (P[i]);
      (*__gmp_free_func) (P, 3*(k+2)*sizeof(mpz_t));
      (*__gmp_free_func) (mult, 2*(k+2)*sizeof(mp_prec_t));

      if (shift_x > 0)
        {
          MPFR_BLOCK (flags, {
//.........这里部分代码省略.........

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