C++ MPFR_GET_EXP函数代码示例

int
mpfr_eint (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd)
{
  int inex;
  mpfr_t tmp, ump;
  mp_exp_t err, te;
  mp_prec_t prec;
  MPFR_SAVE_EXPO_DECL (expo);
  MPFR_ZIV_DECL (loop);

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

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      /* exp(NaN) = exp(-Inf) = NaN */
      if (MPFR_IS_NAN (x) || (MPFR_IS_INF (x) && MPFR_IS_NEG(x)))
        {
          MPFR_SET_NAN (y);
          MPFR_RET_NAN;
        }
      /* eint(+inf) = +inf */
      else if (MPFR_IS_INF (x))
        {
          MPFR_SET_INF(y);
          MPFR_SET_POS(y);
          MPFR_RET(0);
        }
      else /* eint(+/-0) = -Inf */
        {
          MPFR_SET_INF(y);
          MPFR_SET_NEG(y);
          MPFR_RET(0);
        }
    }

  /* eint(x) = NaN for x < 0 */
  if (MPFR_IS_NEG(x))
    {
      MPFR_SET_NAN (y);
      MPFR_RET_NAN;
    }

  MPFR_SAVE_EXPO_MARK (expo);

  /* Since eint(x) >= exp(x)/x, we have log2(eint(x)) >= (x-log(x))/log(2).
     Let's compute k <= (x-log(x))/log(2) in a low precision. If k >= emax,
     then log2(eint(x)) >= emax, and eint(x) >= 2^emax, i.e. it overflows. */
  mpfr_init2 (tmp, 64);
  mpfr_init2 (ump, 64);
  mpfr_log (tmp, x, GMP_RNDU);
  mpfr_sub (ump, x, tmp, GMP_RNDD);
  mpfr_const_log2 (tmp, GMP_RNDU);
  mpfr_div (ump, ump, tmp, GMP_RNDD);
  /* FIXME: We really need mpfr_set_exp_t and mpfr_cmp_exp_t functions. */
  MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
  if (mpfr_cmp_ui (ump, __gmpfr_emax) >= 0)
    {
      mpfr_clear (tmp);
      mpfr_clear (ump);
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_overflow (y, rnd, 1);
    }

  /* Init stuff */
  prec = MPFR_PREC (y) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)) + 6;

  /* eint() has a root 0.37250741078136663446..., so if x is near,
     already take more bits */
  if (MPFR_GET_EXP(x) == -1) /* 1/4 <= x < 1/2 */
    {
      double d;
      d = mpfr_get_d (x, GMP_RNDN) - 0.37250741078136663;
      d = (d == 0.0) ? -53 : __gmpfr_ceil_log2 (d);
      prec += -d;
    }

  mpfr_set_prec (tmp, prec);
  mpfr_set_prec (ump, prec);

  MPFR_ZIV_INIT (loop, prec);            /* Initialize the ZivLoop controler */
  for (;;)                               /* Infinite loop */
    {
      /* We need that the smallest value of k!/x^k is smaller than 2^(-p).
         The minimum is obtained for x=k, and it is smaller than e*sqrt(x)/e^x
         for x>=1. */
      if (MPFR_GET_EXP (x) > 0 && mpfr_cmp_d (x, ((double) prec +
                            0.5 * (double) MPFR_GET_EXP (x)) * LOG2 + 1.0) > 0)
        err = mpfr_eint_asympt (tmp, x);
      else
        {
          err = mpfr_eint_aux (tmp, x); /* error <= 2^err ulp(tmp) */
          te = MPFR_GET_EXP(tmp);
          mpfr_const_euler (ump, GMP_RNDN); /* 0.577 -> EXP(ump)=0 */
          mpfr_add (tmp, tmp, ump, GMP_RNDN);
          /* error <= 1/2 + 1/2*2^(EXP(ump)-EXP(tmp)) + 2^(te-EXP(tmp)+err)
             <= 1/2 + 2^(MAX(EXP(ump), te+err+1) - EXP(tmp))
             <= 2^(MAX(0, 1 + MAX(EXP(ump), te+err+1) - EXP(tmp))) */
          err = MAX(1, te + err + 2) - MPFR_GET_EXP(tmp);
          err = MAX(0, err);
//.........这里部分代码省略.........

MPFR_HOT_FUNCTION_ATTR int
mpfr_cmp3 (mpfr_srcptr b, mpfr_srcptr c, int s)
{
  mpfr_exp_t be, ce;
  mp_size_t bn, cn;
  mp_limb_t *bp, *cp;

  s = MPFR_MULT_SIGN( s , MPFR_SIGN(c) );

  if (MPFR_ARE_SINGULAR(b, c))
    {
      if (MPFR_IS_NAN (b) || MPFR_IS_NAN (c))
        {
          MPFR_SET_ERANGEFLAG ();
          return 0;
        }
      else if (MPFR_IS_INF(b))
        {
          if (MPFR_IS_INF(c) && s == MPFR_SIGN(b) )
            return 0;
          else
            return MPFR_SIGN(b);
        }
      else if (MPFR_IS_INF(c))
        return -s;
      else if (MPFR_IS_ZERO(b))
        return MPFR_IS_ZERO(c) ? 0 : -s;
      else /* necessarily c=0 */
        return MPFR_SIGN(b);
    }
  /* b and c are real numbers */
  if (s != MPFR_SIGN(b))
    return MPFR_SIGN(b);

  /* now signs are equal */

  be = MPFR_GET_EXP (b);
  ce = MPFR_GET_EXP (c);
  if (be > ce)
    return s;
  if (be < ce)
    return -s;

  /* both signs and exponents are equal */

  bn = MPFR_LAST_LIMB (b);
  cn = MPFR_LAST_LIMB (c);

  bp = MPFR_MANT(b);
  cp = MPFR_MANT(c);

  for ( ; bn >= 0 && cn >= 0; bn--, cn--)
    {
      if (bp[bn] > cp[cn])
        return s;
      if (bp[bn] < cp[cn])
        return -s;
    }
  for ( ; bn >= 0; bn--)
    if (bp[bn])
      return s;
  for ( ; cn >= 0; cn--)
    if (cp[cn])
      return -s;

   return 0;
}

int
mpfr_cmp_ui_2exp (mpfr_srcptr b, unsigned long int i, mp_exp_t f)
{
  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(b) ))
    {
      if (MPFR_IS_NAN (b))
        {
          MPFR_SET_ERANGE ();
          return 0;
        }
      else if (MPFR_IS_INF(b))
        return MPFR_INT_SIGN (b);
      else /* since b cannot be NaN, b=0 here */
        return i != 0 ? -1 : 0;
    }

  if (MPFR_IS_NEG (b))
    return -1;
  /* now b > 0 */
  else if (MPFR_UNLIKELY(i == 0))
    return 1;
  else /* b > 0, i > 0 */
    {
      mp_exp_t e;
      int k;
      mp_size_t bn;
      mp_limb_t c, *bp;

      /* i must be representable in a mp_limb_t */
      MPFR_ASSERTN(i == (mp_limb_t) i);

      e = MPFR_GET_EXP (b); /* 2^(e-1) <= b < 2^e */
      if (e <= f)
        return -1;
      if (f < MPFR_EMAX_MAX - BITS_PER_MP_LIMB &&
          e > f + BITS_PER_MP_LIMB)
        return 1;

      /* now f < e <= f + BITS_PER_MP_LIMB */
      c = (mp_limb_t) i;
      count_leading_zeros(k, c);
      if ((int) (e - f) > BITS_PER_MP_LIMB - k)
        return 1;
      if ((int) (e - f) < BITS_PER_MP_LIMB - k)
        return -1;

      /* now b and i*2^f have the same exponent */
      c <<= k;
      bn = (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
      bp = MPFR_MANT(b);
      if (bp[bn] > c)
        return 1;
      if (bp[bn] < c)
        return -1;

      /* most significant limbs agree, check remaining limbs from b */
      while (bn > 0)
        if (bp[--bn] != 0)
          return 1;
      return 0;
    }
}

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

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

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

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

/* return non zero iff x^y is exact.
   Assumes x and y are ordinary numbers,
   y is not an integer, x is not a power of 2 and x is positive

   If x^y is exact, it computes it and sets *inexact.
*/
static int
mpfr_pow_is_exact (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
                   mpfr_rnd_t rnd_mode, int *inexact)
{
  mpz_t a, c;
  mpfr_exp_t d, b;
  unsigned long i;
  int res;

  MPFR_ASSERTD (!MPFR_IS_SINGULAR (y));
  MPFR_ASSERTD (!MPFR_IS_SINGULAR (x));
  MPFR_ASSERTD (!mpfr_integer_p (y));
  MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_INT_SIGN (x),
                                  MPFR_GET_EXP (x) - 1) != 0);
  MPFR_ASSERTD (MPFR_IS_POS (x));

  if (MPFR_IS_NEG (y))
    return 0; /* x is not a power of two => x^-y is not exact */

  /* compute d such that y = c*2^d with c odd integer */
  mpz_init (c);
  d = mpfr_get_z_2exp (c, y);
  i = mpz_scan1 (c, 0);
  mpz_fdiv_q_2exp (c, c, i);
  d += i;
  /* now y=c*2^d with c odd */
  /* Since y is not an integer, d is necessarily < 0 */
  MPFR_ASSERTD (d < 0);

  /* Compute a,b such that x=a*2^b */
  mpz_init (a);
  b = mpfr_get_z_2exp (a, x);
  i = mpz_scan1 (a, 0);
  mpz_fdiv_q_2exp (a, a, i);
  b += i;
  /* now x=a*2^b with a is odd */

  for (res = 1 ; d != 0 ; d++)
    {
      /* a*2^b is a square iff
            (i)  a is a square when b is even
            (ii) 2*a is a square when b is odd */
      if (b % 2 != 0)
        {
          mpz_mul_2exp (a, a, 1); /* 2*a */
          b --;
        }
      MPFR_ASSERTD ((b % 2) == 0);
      if (!mpz_perfect_square_p (a))
        {
          res = 0;
          goto end;
        }
      mpz_sqrt (a, a);
      b = b / 2;
    }
  /* Now x = (a'*2^b')^(2^-d) with d < 0
     so x^y = ((a'*2^b')^(2^-d))^(c*2^d)
            = ((a'*2^b')^c with c odd integer */
  {
    mpfr_t tmp;
    mpfr_prec_t p;
    MPFR_MPZ_SIZEINBASE2 (p, a);
    mpfr_init2 (tmp, p); /* prec = 1 should not be possible */
    res = mpfr_set_z (tmp, a, MPFR_RNDN);
    MPFR_ASSERTD (res == 0);
    res = mpfr_mul_2si (tmp, tmp, b, MPFR_RNDN);
    MPFR_ASSERTD (res == 0);
    *inexact = mpfr_pow_z (z, tmp, c, rnd_mode);
    mpfr_clear (tmp);
    res = 1;
  }
 end:
  mpz_clear (a);
  mpz_clear (c);
  return res;
}

/* Input: s - a floating-point number >= 1/2.
          rnd_mode - a rounding mode.
          Assumes s is neither NaN nor Infinite.
   Output: z - Zeta(s) rounded to the precision of z with direction rnd_mode
*/
static int
mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
{
  mpfr_t b, c, z_pre, f, s1;
  double beta, sd, dnep;
  mpfr_t *tc1;
  mp_prec_t precz, precs, d, dint;
  int p, n, l, add;
  int inex;
  MPFR_GROUP_DECL (group);
  MPFR_ZIV_DECL (loop);

  MPFR_ASSERTD (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0);

  precz = MPFR_PREC (z);
  precs = MPFR_PREC (s);

  /* Zeta(x) = 1+1/2^x+1/3^x+1/4^x+1/5^x+O(1/6^x)
     so with 2^(EXP(x)-1) <= x < 2^EXP(x)
     So for x > 2^3, k^x > k^8, so 2/k^x < 2/k^8
     Zeta(x) = 1 + 1/2^x*(1+(2/3)^x+(2/4)^x+...)
             = 1 + 1/2^x*(1+sum((2/k)^x,k=3..infinity))
            <= 1 + 1/2^x*(1+sum((2/k)^8,k=3..infinity))
     And sum((2/k)^8,k=3..infinity) = -257+128*Pi^8/4725 ~= 0.0438035
     So Zeta(x) <= 1 + 1/2^x*2 for x >= 8
     The error is < 2^(-x+1) <= 2^(-2^(EXP(x)-1)+1) */
  if (MPFR_GET_EXP (s) > 3)
    {
      mp_exp_t err;
      err = MPFR_GET_EXP (s) - 1;
      if (err > (mp_exp_t) (sizeof (mp_exp_t)*CHAR_BIT-2))
        err = MPFR_EMAX_MAX;
      else
        err = ((mp_exp_t)1) << err;
      err = 1 - (-err+1); /* GET_EXP(one) - (-err+1) = err :) */
      MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, err, 0, 1,
                                        rnd_mode, {});
    }

  d = precz + MPFR_INT_CEIL_LOG2(precz) + 10;

  /* we want that s1 = s-1 is exact, i.e. we should have PREC(s1) >= EXP(s) */
  dint = (mpfr_uexp_t) MPFR_GET_EXP (s);
  mpfr_init2 (s1, MAX (precs, dint));
  inex = mpfr_sub (s1, s, __gmpfr_one, GMP_RNDN);
  MPFR_ASSERTD (inex == 0);

  /* case s=1 */
  if (MPFR_IS_ZERO (s1))
    {
      MPFR_SET_INF (z);
      MPFR_SET_POS (z);
      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
//.........这里部分代码省略.........

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

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN (x))
        {
          MPFR_SET_NAN (y);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF (x))
        {
          if (MPFR_IS_POS (x))
            MPFR_SET_INF (y);
          else
            MPFR_SET_ZERO (y);
          MPFR_SET_POS (y);
          MPFR_RET (0);
        }
      else /* 2^0 = 1 */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO(x));
          return mpfr_set_ui (y, 1, rnd_mode);
        }
    }

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

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

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

  MPFR_SAVE_EXPO_MARK (expo);

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

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

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

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

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

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

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

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

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

//.........这里部分代码省略.........
  k = 0;
  for (i = 1 ; i < n ; i++) {
    k++;

#ifdef A
#  ifdef B
    mpz_set_ui (T[k], (A1 + A2*i)*(B1+B2*i));
#  else
    mpz_set_ui (T[k], A1 + A2*i);
#  endif
#endif

#ifdef C
#  ifdef NO_FACTORIAL
    mpz_set_ui (P[k], (C1 + C2 * (i-1)));
    mpz_set_ui (S[k], 1);
#  else
    mpz_set_ui (P[k], (i+1) * (C1 + C2 * (i-1)));
    mpz_set_ui (S[k], i+1);
#  endif
#else
#  ifdef NO_FACTORIAL
    mpz_set_ui (P[k], 1);
#  else
    mpz_set_ui (P[k], i+1);
#  endif
    mpz_set (S[k], P[k]);
#endif

    for (j = i+1, l = 0 ; (j & 1) == 0 ; l++, j>>=1, k--) {
      if (!is_p_one)
        mpz_mul (S[k], S[k], ptoj[l]);
#ifdef A
#  ifdef B
#    if (A2*B2) != 1
      mpz_mul_ui (P[k], P[k], A2*B2);
#    endif
#  else
#    if A2 != 1
      mpz_mul_ui (P[k], P[k], A2);
#  endif
#endif
      mpz_mul (S[k], S[k], T[k-1]);
#endif
      mpz_mul (S[k-1], S[k-1], P[k]);
#ifdef R_IS_RATIONAL
      mpz_mul (S[k-1], S[k-1], qtoj[l]);
#else
      mpz_mul_2exp (S[k-1], S[k-1], r*(1<<l));
#endif
      mpz_add (S[k-1], S[k-1], S[k]);
      mpz_mul (P[k-1], P[k-1], P[k]);
#ifdef A
      mpz_mul (T[k-1], T[k-1], T[k]);
#endif
    }
  }

  diff = mpz_sizeinbase(S[0],2) - 2*precy;
  expo = diff;
  if (diff >= 0)
    mpz_div_2exp(S[0],S[0],diff);
  else
    mpz_mul_2exp(S[0],S[0],-diff);
  diff = mpz_sizeinbase(P[0],2) - precy;
  expo -= diff;
  if (diff >=0)
    mpz_div_2exp(P[0],P[0],diff);
  else
    mpz_mul_2exp(P[0],P[0],-diff);

  mpz_tdiv_q(S[0], S[0], P[0]);
  mpfr_set_z(y, S[0], GMP_RNDD);
  MPFR_SET_EXP (y, MPFR_GET_EXP (y) + expo);

#ifdef R_IS_RATIONAL
  /* exact division */
  mpz_div_ui (qtoj[m], qtoj[m], r);
  mpfr_init2 (tmp, MPFR_PREC(y));
  mpfr_set_z (tmp, qtoj[m] , GMP_RNDD);
  mpfr_div (y, y, tmp, GMP_RNDD);
  mpfr_clear (tmp);
#else
  mpfr_div_2ui(y, y, r*(i-1), GMP_RNDN);
#endif
  for (i = 0 ; i <= m ; i++)
    {
      mpz_clear (P[i]);
      mpz_clear (S[i]);
      mpz_clear (ptoj[i]);
#ifdef R_IS_RATIONAL
      mpz_clear (qtoj[i]);
#endif
#ifdef A
      mpz_clear (T[i]);
#endif
    }
  MPFR_TMP_FREE (marker);
  return 0;
}

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

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

  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); /* log2(0) is an exact -infinity */
        }
    }

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

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

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

  MPFR_SAVE_EXPO_MARK (expo);

  /* General case */
  {
    /* Declaration of the intermediary variable */
    mpfr_t t, tt;
    /* Declaration of the size variable */
    mpfr_prec_t Ny = MPFR_PREC(r);              /* 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 + 3 + MPFR_INT_CEIL_LOG2 (Ny);

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

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

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

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

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

double
mpfr_get_d (mpfr_srcptr src, mpfr_rnd_t rnd_mode)
{
    double d;
    int negative;
    mpfr_exp_t e;

    if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
    {
        if (MPFR_IS_NAN (src))
            return MPFR_DBL_NAN;

        negative = MPFR_IS_NEG (src);

        if (MPFR_IS_INF (src))
            return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;

        MPFR_ASSERTD (MPFR_IS_ZERO(src));
        return negative ? DBL_NEG_ZERO : 0.0;
    }

    e = MPFR_GET_EXP (src);
    negative = MPFR_IS_NEG (src);

    if (MPFR_UNLIKELY(rnd_mode == MPFR_RNDA))
        rnd_mode = negative ? MPFR_RNDD : MPFR_RNDU;

    /* the smallest normalized number is 2^(-1022)=0.1e-1021, and the smallest
       subnormal is 2^(-1074)=0.1e-1073 */
    if (MPFR_UNLIKELY (e < -1073))
    {
        /* Note: Avoid using a constant expression DBL_MIN * DBL_EPSILON
           as this gives 0 instead of the correct result with gcc on some
           Alpha machines. */
        d = negative ?
            (rnd_mode == MPFR_RNDD ||
             (rnd_mode == MPFR_RNDN && mpfr_cmp_si_2exp(src, -1, -1075) < 0)
             ? -DBL_MIN : DBL_NEG_ZERO) :
            (rnd_mode == MPFR_RNDU ||
             (rnd_mode == MPFR_RNDN && mpfr_cmp_si_2exp(src, 1, -1075) > 0)
             ? DBL_MIN : 0.0);
        if (d != 0.0) /* we multiply DBL_MIN = 2^(-1022) by DBL_EPSILON = 2^(-52)
                       to get +-2^(-1074) */
            d *= DBL_EPSILON;
    }
    /* the largest normalized number is 2^1024*(1-2^(-53))=0.111...111e1024 */
    else if (MPFR_UNLIKELY (e > 1024))
    {
        d = negative ?
            (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDU ?
             -DBL_MAX : MPFR_DBL_INFM) :
            (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDD ?
             DBL_MAX : MPFR_DBL_INFP);
    }
    else
    {
        int nbits;
        mp_size_t np, i;
        mp_limb_t tp[ MPFR_LIMBS_PER_DOUBLE ];
        int carry;

        nbits = IEEE_DBL_MANT_DIG; /* 53 */
        if (MPFR_UNLIKELY (e < -1021))
            /*In the subnormal case, compute the exact number of significant bits*/
        {
            nbits += (1021 + e);
            MPFR_ASSERTD (nbits >= 1);
        }
        np = MPFR_PREC2LIMBS (nbits);
        MPFR_ASSERTD ( np <= MPFR_LIMBS_PER_DOUBLE );
        carry = mpfr_round_raw_4 (tp, MPFR_MANT(src), MPFR_PREC(src), negative,
                                  nbits, rnd_mode);
        if (MPFR_UNLIKELY(carry))
            d = 1.0;
        else
        {
            /* The following computations are exact thanks to the previous
               mpfr_round_raw. */
            d = (double) tp[0] / MP_BASE_AS_DOUBLE;
            for (i = 1 ; i < np ; i++)
                d = (d + tp[i]) / MP_BASE_AS_DOUBLE;
            /* d is the mantissa (between 1/2 and 1) of the argument rounded
               to 53 bits */
        }
        d = mpfr_scale2 (d, e);
        if (negative)
            d = -d;
    }

    return d;
}

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[%#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
        {
          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_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);
      for (k = 0; k < K; k++)
        {
//.........这里部分代码省略.........

/* compute in y an approximation of sum(x^k/k/k!, k=1..infinity),
   and return e such that the absolute error is bound by 2^e ulp(y) */
static mp_exp_t
mpfr_eint_aux (mpfr_t y, mpfr_srcptr x)
{
  mpfr_t eps; /* dynamic (absolute) error bound on t */
  mpfr_t erru, errs;
  mpz_t m, s, t, u;
  mp_exp_t e, sizeinbase;
  mp_prec_t w = MPFR_PREC(y);
  unsigned long k;
  MPFR_GROUP_DECL (group);

  /* for |x| <= 1, we have S := sum(x^k/k/k!, k=1..infinity) = x + R(x)
     where |R(x)| <= (x/2)^2/(1-x/2) <= 2*(x/2)^2
     thus |R(x)/x| <= |x|/2
     thus if |x| <= 2^(-PREC(y)) we have |S - o(x)| <= ulp(y) */

  if (MPFR_GET_EXP(x) <= - (mp_exp_t) w)
    {
      mpfr_set (y, x, GMP_RNDN);
      return 0;
    }

  mpz_init (s); /* initializes to 0 */
  mpz_init (t);
  mpz_init (u);
  mpz_init (m);
  MPFR_GROUP_INIT_3 (group, 31, eps, erru, errs);
  e = mpfr_get_z_exp (m, x); /* x = m * 2^e */
  MPFR_ASSERTD (mpz_sizeinbase (m, 2) == MPFR_PREC (x));
  if (MPFR_PREC (x) > w)
    {
      e += MPFR_PREC (x) - w;
      mpz_tdiv_q_2exp (m, m, MPFR_PREC (x) - w);
    }
  /* remove trailing zeroes from m: this will speed up much cases where
     x is a small integer divided by a power of 2 */
  k = mpz_scan1 (m, 0);
  mpz_tdiv_q_2exp (m, m, k);
  e += k;
  /* initialize t to 2^w */
  mpz_set_ui (t, 1);
  mpz_mul_2exp (t, t, w);
  mpfr_set_ui (eps, 0, GMP_RNDN); /* eps[0] = 0 */
  mpfr_set_ui (errs, 0, GMP_RNDN);
  for (k = 1;; k++)
    {
      /* let eps[k] be the absolute error on t[k]:
         since t[k] = trunc(t[k-1]*m*2^e/k), we have
         eps[k+1] <= 1 + eps[k-1]*m*2^e/k + t[k-1]*m*2^(1-w)*2^e/k
                  =  1 + (eps[k-1] + t[k-1]*2^(1-w))*m*2^e/k
                  = 1 + (eps[k-1]*2^(w-1) + t[k-1])*2^(1-w)*m*2^e/k */
      mpfr_mul_2ui (eps, eps, w - 1, GMP_RNDU);
      mpfr_add_z (eps, eps, t, GMP_RNDU);
      MPFR_MPZ_SIZEINBASE2 (sizeinbase, m);
      mpfr_mul_2si (eps, eps, sizeinbase - (w - 1) + e, GMP_RNDU);
      mpfr_div_ui (eps, eps, k, GMP_RNDU);
      mpfr_add_ui (eps, eps, 1, GMP_RNDU);
      mpz_mul (t, t, m);
      if (e < 0)
        mpz_tdiv_q_2exp (t, t, -e);
      else
        mpz_mul_2exp (t, t, e);
      mpz_tdiv_q_ui (t, t, k);
      mpz_tdiv_q_ui (u, t, k);
      mpz_add (s, s, u);
      /* the absolute error on u is <= 1 + eps[k]/k */
      mpfr_div_ui (erru, eps, k, GMP_RNDU);
      mpfr_add_ui (erru, erru, 1, GMP_RNDU);
      /* and that on s is the sum of all errors on u */
      mpfr_add (errs, errs, erru, GMP_RNDU);
      /* we are done when t is smaller than errs */
      if (mpz_sgn (t) == 0)
        sizeinbase = 0;
      else
        MPFR_MPZ_SIZEINBASE2 (sizeinbase, t);
      if (sizeinbase < MPFR_GET_EXP (errs))
        break;
    }
  /* the truncation error is bounded by (|t|+eps)/k*(|x|/k + |x|^2/k^2 + ...)
     <= (|t|+eps)/k*|x|/(k-|x|) */
  mpz_abs (t, t);
  mpfr_add_z (eps, eps, t, GMP_RNDU);
  mpfr_div_ui (eps, eps, k, GMP_RNDU);
  mpfr_abs (erru, x, GMP_RNDU); /* |x| */
  mpfr_mul (eps, eps, erru, GMP_RNDU);
  mpfr_ui_sub (erru, k, erru, GMP_RNDD);
  if (MPFR_IS_NEG (erru))
    {
      /* the truncated series does not converge, return fail */
      e = w;
    }
  else
    {
      mpfr_div (eps, eps, erru, GMP_RNDU);
      mpfr_add (errs, errs, eps, GMP_RNDU);
      mpfr_set_z (y, s, GMP_RNDN);
      mpfr_div_2ui (y, y, w, GMP_RNDN);
      /* errs was an absolute error bound on s. We must convert it to an error
//.........这里部分代码省略.........

/* (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, mpfr_rnd_t rnd_mode)
{
  mpfr_prec_t prec, m;
  int neg, reduce;
  mpfr_t c, xr;
  mpfr_srcptr xx;
  mpfr_exp_t err, expx;
  int inexy, inexz;
  MPFR_ZIV_DECL (loop);
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_ASSERTN (y != z);

  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 */
          inexy = 0; /* y is exact */
          inexz = mpfr_set_ui (z, 1, rnd_mode);
          return INEX(inexy,inexz);
        }
    }

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

  MPFR_SAVE_EXPO_MARK (expo);

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

  /* When x is close to 0, say 2^(-k), then there is a cancellation of about
     2k bits in 1-cos(x)^2. FIXME: in that case, it would be more efficient
     to compute sin(x) directly. VL: This is partly done by using
     MPFR_FAST_COMPUTE_IF_SMALL_INPUT from the mpfr_sin and mpfr_cos
     functions. Moreover, any overflow on m is avoided. */
  if (expx < 0)
    {
      /* Warning: in case y = x, and the first call to
         MPFR_FAST_COMPUTE_IF_SMALL_INPUT succeeds but the second fails,
         we will have clobbered the original value of x.
         The workaround is to first compute z = cos(x) in that case, since
         y and z are different. */
      if (y != x)
        /* y and x differ, thus we can safely try to compute y first */
        {
          MPFR_FAST_COMPUTE_IF_SMALL_INPUT (
            y, x, -2 * expx, 2, 0, rnd_mode,
            { inexy = _inexact;
              goto small_input; });