int
main (void)
{
  mpfr_t x;
  FILE *f;
  int i;

  tests_start_mpfr ();

  /* Check RND_MODE */
  if (strcmp (mpfr_print_rnd_mode(MPFR_RNDN), "MPFR_RNDN"))
    {
      printf ("Error for printing MPFR_RNDN\n");
      exit (1);
    }
  if (strcmp (mpfr_print_rnd_mode(MPFR_RNDU), "MPFR_RNDU"))
    {
      printf ("Error for printing MPFR_RNDU\n");
      exit (1);
    }
  if (strcmp (mpfr_print_rnd_mode(MPFR_RNDD), "MPFR_RNDD"))
    {
      printf ("Error for printing MPFR_RNDD\n");
      exit (1);
    }
  if (strcmp (mpfr_print_rnd_mode(MPFR_RNDZ), "MPFR_RNDZ"))
    {
      printf ("Error for printing MPFR_RNDZ\n");
      exit (1);
    }
  if (mpfr_print_rnd_mode ((mpfr_rnd_t) -1) != NULL ||
      mpfr_print_rnd_mode (MPFR_RND_MAX) != NULL)
    {
      printf ("Error for illegal rounding mode values.\n");
      exit (1);
    }

  /* Reopen stdout to a file. All errors will be put to stderr
     Can't use tmpname since it is unsecure */
  if (freopen (FILE_NAME, "w", stdout) == NULL)
    {
      printf ("Error can't redirect stdout\n");
      exit (1);
    }
  mpfr_init (x);
  mpfr_set_nan (x);
  mpfr_dump (x);
  mpfr_set_inf (x, -1);
  mpfr_dump (x);
  MPFR_SET_ZERO (x); MPFR_SET_NEG (x);
  mpfr_dump (x);
  mpfr_set_str_binary (x, "0.101010101010111110010001100011000100001E32");
  mpfr_dump (x);
  mpfr_print_mant_binary ("x=",MPFR_MANT(x), MPFR_PREC(x));


  mpfr_clear (x);
  fclose (stdout);
  /* Open it and check for it */
  f = fopen (FILE_NAME, "r");
  if (f == NULL)
    {
      fprintf (stderr, "Can't reopen file!\n");
      exit (1);
    }
  for(i = 0 ; i < sizeof(Buffer)-1 ; i++)
    {
      if (feof (f))
        {
          fprintf (stderr, "Error EOF\n");
          exit (1);
        }
      if (Buffer[i] != fgetc (f))
        {
          fprintf (stderr, "Character mismatch for i=%d / %lu\n",
                  i, (unsigned long) sizeof(Buffer));
          exit (1);
        }
    }
  fclose (f);

  remove (FILE_NAME);
  tests_end_mpfr ();
  return 0;
}

int
main (void)
{
  mpfr_t x, y;
  int i, r;

  tests_start_mpfr ();

  mpfr_init2 (x, 256);
  mpfr_init2 (y, 8);

  RND_LOOP (r)
    {

      /* Check NAN */
      mpfr_set_nan (x);
      if (mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
        ERROR1 (1);
      if (mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
        ERROR1 (2);
      if (mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
        ERROR1 (3);
      if (mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
        ERROR1 (4);
      if (mpfr_fits_ushort_p (x, (mpfr_rnd_t) r))
        ERROR1 (5);
      if (mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
        ERROR1 (6);

      /* Check INF */
      mpfr_set_inf (x, 1);
      if (mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
        ERROR1 (7);
      if (mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
        ERROR1 (8);
      if (mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
        ERROR1 (9);
      if (mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
        ERROR1 (10);
      if (mpfr_fits_ushort_p (x, (mpfr_rnd_t) r))
        ERROR1 (11);
      if (mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
        ERROR1 (12);

      /* Check Zero */
      MPFR_SET_ZERO (x);
      if (!mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
        ERROR1 (13);
      if (!mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
        ERROR1 (14);
      if (!mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
        ERROR1 (15);
      if (!mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
        ERROR1 (16);
      if (!mpfr_fits_ushort_p (x, (mpfr_rnd_t) r))
        ERROR1 (17);
      if (!mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
        ERROR1 (18);

      /* Check small positive op */
      mpfr_set_str1 (x, "[email protected]");
      if (!mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
        ERROR1 (19);
      if (!mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
        ERROR1 (20);
      if (!mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
        ERROR1 (21);
      if (!mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
        ERROR1 (22);
      if (!mpfr_fits_ushort_p (x, (mpfr_rnd_t) r))
        ERROR1 (23);
      if (!mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
        ERROR1 (24);

      /* Check 17 */
      mpfr_set_ui (x, 17, MPFR_RNDN);
      if (!mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
        ERROR1 (25);
      if (!mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
        ERROR1 (26);
      if (!mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
        ERROR1 (27);
      if (!mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
        ERROR1 (28);
      if (!mpfr_fits_ushort_p (x, (mpfr_rnd_t) r))
        ERROR1 (29);
      if (!mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
        ERROR1 (30);

      /* Check all other values */
      mpfr_set_ui (x, ULONG_MAX, MPFR_RNDN);
      mpfr_mul_2exp (x, x, 1, MPFR_RNDN);
      if (mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
        ERROR1 (31);
      if (mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
        ERROR1 (32);
      mpfr_mul_2exp (x, x, 40, MPFR_RNDN);
      if (mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
        ERROR1 (33);
      if (mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
//.........这里部分代码省略.........

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

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

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

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

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

  MPFR_ASSERTN ((mpfr_uprec_t) bq + cq <= MPFR_PREC_MAX);

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

  /* Check for no size_t overflow*/
  MPFR_ASSERTD((size_t) k <= ((size_t) -1) / MPFR_BYTES_PER_MP_LIMB);
  MPFR_TMP_MARK(marker);
  tmp = MPFR_TMP_LIMBS_ALLOC (k);

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

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

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

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

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

static void check_intmax (void)
{
#ifdef _MPFR_H_HAVE_INTMAX_T
    mpfr_t x;

    mpfr_init2 (x, sizeof (uintmax_t)*CHAR_BIT);

    /* Check NAN */
    mpfr_set_nan (x);
    if (mpfr_fits_uintmax_p (x, MPFR_RNDN))
        ERROR1;
    if (mpfr_fits_intmax_p (x, MPFR_RNDN))
        ERROR1;

    /* Check INF */
    mpfr_set_inf (x, 1);
    if (mpfr_fits_uintmax_p (x, MPFR_RNDN))
        ERROR1;
    if (mpfr_fits_intmax_p (x, MPFR_RNDN))
        ERROR1;

    /* Check Zero */
    MPFR_SET_ZERO (x);
    if (!mpfr_fits_uintmax_p (x, MPFR_RNDN))
        ERROR2;
    if (!mpfr_fits_intmax_p (x, MPFR_RNDN))
        ERROR2;

    /* Check small op */
    mpfr_set_str1 (x, "[email protected]");
    if (!mpfr_fits_uintmax_p (x, MPFR_RNDN))
        ERROR2;
    if (!mpfr_fits_intmax_p (x, MPFR_RNDN))
        ERROR2;

    /* Check 17 */
    mpfr_set_ui (x, 17, MPFR_RNDN);
    if (!mpfr_fits_uintmax_p (x, MPFR_RNDN))
        ERROR2;
    if (!mpfr_fits_intmax_p (x, MPFR_RNDN))
        ERROR2;

    /* Check hugest */
    mpfr_set_ui_2exp (x, 42, sizeof (uintmax_t) * 32, MPFR_RNDN);
    if (mpfr_fits_uintmax_p (x, MPFR_RNDN))
        ERROR1;
    if (mpfr_fits_intmax_p (x, MPFR_RNDN))
        ERROR1;

    /* Check all other values */
    mpfr_set_uj (x, MPFR_UINTMAX_MAX, MPFR_RNDN);
    mpfr_add_ui (x, x, 1, MPFR_RNDN);
    if (mpfr_fits_uintmax_p (x, MPFR_RNDN))
        ERROR1;
    mpfr_set_uj (x, MPFR_UINTMAX_MAX, MPFR_RNDN);
    if (!mpfr_fits_uintmax_p (x, MPFR_RNDN))
        ERROR2;
    mpfr_set_sj (x, MPFR_INTMAX_MAX, MPFR_RNDN);
    mpfr_add_ui (x, x, 1, MPFR_RNDN);
    if (mpfr_fits_intmax_p (x, MPFR_RNDN))
        ERROR1;
    mpfr_set_sj (x, MPFR_INTMAX_MAX, MPFR_RNDN);
    if (!mpfr_fits_intmax_p (x, MPFR_RNDN))
        ERROR2;
    mpfr_set_sj (x, MPFR_INTMAX_MIN, MPFR_RNDN);
    if (!mpfr_fits_intmax_p (x, MPFR_RNDN))
        ERROR2;
    mpfr_sub_ui (x, x, 1, MPFR_RNDN);
    if (mpfr_fits_intmax_p (x, MPFR_RNDN))
        ERROR1;

    /* Check negative value */
    mpfr_set_si (x, -1, MPFR_RNDN);
    if (!mpfr_fits_intmax_p (x, MPFR_RNDN))
        ERROR2;
    if (mpfr_fits_uintmax_p (x, MPFR_RNDN))
        ERROR1;

    mpfr_clear (x);
#endif
}

int
mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
  mpfr_exp_t bx,cx;
  mpfr_uexp_t d;
  mpfr_prec_t p, sh, cnt;
  mp_size_t n;
  mp_limb_t *ap, *bp, *cp;
  mp_limb_t limb;
  int inexact;
  mp_limb_t bcp,bcp1; /* Cp and C'p+1 */
  mp_limb_t bbcp = (mp_limb_t) -1, bbcp1 = (mp_limb_t) -1; /* Cp+1 and C'p+2,
    gcc claims that they might be used uninitialized. We fill them with invalid
    values, which should produce a failure if so. See README.dev file. */

  MPFR_TMP_DECL(marker);

  MPFR_TMP_MARK(marker);

  MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c));
  MPFR_ASSERTD(MPFR_IS_PURE_FP(b));
  MPFR_ASSERTD(MPFR_IS_PURE_FP(c));

  /* Read prec and num of limbs */
  p = MPFR_PREC (b);
  n = MPFR_PREC2LIMBS (p);

  /* Fast cmp of |b| and |c|*/
  bx = MPFR_GET_EXP (b);
  cx = MPFR_GET_EXP (c);
  if (MPFR_UNLIKELY(bx == cx))
    {
      mp_size_t k = n - 1;
      /* Check mantissa since exponent are equals */
      bp = MPFR_MANT(b);
      cp = MPFR_MANT(c);
      while (k>=0 && MPFR_UNLIKELY(bp[k] == cp[k]))
        k--;
      if (MPFR_UNLIKELY(k < 0))
        /* b == c ! */
        {
          /* Return exact number 0 */
          if (rnd_mode == MPFR_RNDD)
            MPFR_SET_NEG(a);
          else
            MPFR_SET_POS(a);
          MPFR_SET_ZERO(a);
          MPFR_RET(0);
        }
      else if (bp[k] > cp[k])
        goto BGreater;
      else
        {
          MPFR_ASSERTD(bp[k]<cp[k]);
          goto CGreater;
        }
    }
  else if (MPFR_UNLIKELY(bx < cx))
    {
      /* Swap b and c and set sign */
      mpfr_srcptr t;
      mpfr_exp_t tx;
    CGreater:
      MPFR_SET_OPPOSITE_SIGN(a,b);
      t  = b;  b  = c;  c  = t;
      tx = bx; bx = cx; cx = tx;
    }
  else
    {
      /* b > c */
    BGreater:
      MPFR_SET_SAME_SIGN(a,b);
    }

  /* Now b > c */
  MPFR_ASSERTD(bx >= cx);
  d = (mpfr_uexp_t) bx - cx;
  DEBUG (printf ("New with diff=%lu\n", (unsigned long) d));

  if (MPFR_UNLIKELY(d <= 1))
    {
      if (MPFR_LIKELY(d < 1))
        {
          /* <-- b -->
             <-- c --> : exact sub */
          ap = MPFR_MANT(a);
          mpn_sub_n (ap, MPFR_MANT(b), MPFR_MANT(c), n);
          /* Normalize */
        ExactNormalize:
          limb = ap[n-1];
          if (MPFR_LIKELY(limb))
            {
              /* First limb is not zero. */
              count_leading_zeros(cnt, limb);
              /* cnt could be == 0 <= SubD1Lose */
              if (MPFR_LIKELY(cnt))
                {
                  mpn_lshift(ap, ap, n, cnt); /* Normalize number */
                  bx -= cnt; /* Update final expo */
                }
//.........这里部分代码省略.........

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

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

  MPFR_SAVE_EXPO_MARK (expo);

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

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

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

int
main (void)
{
  mpfr_t a;
  mp_limb_t *p, tmp;
  mp_size_t s;
  mpfr_prec_t pr;
  int max;

  tests_start_mpfr ();
  for(pr = MPFR_PREC_MIN ; pr < 500 ; pr++)
    {
      mpfr_init2 (a, pr);
      if (!mpfr_check(a)) ERROR("for init");
      /* Check special cases */
      MPFR_SET_NAN(a);
      if (!mpfr_check(a)) ERROR("for nan");
      MPFR_SET_POS(a);
      MPFR_SET_INF(a);
      if (!mpfr_check(a)) ERROR("for inf");
      MPFR_SET_ZERO(a);
      if (!mpfr_check(a)) ERROR("for zero");
      MPFR_EXP (a) = MPFR_EXP_MIN;
      if (mpfr_check(a))  ERROR("for EXP = MPFR_EXP_MIN");
      /* Check var */
      mpfr_set_ui(a, 2, MPFR_RNDN);
      if (!mpfr_check(a)) ERROR("for set_ui");
      mpfr_clear_overflow();
      max = 1000; /* Allows max 2^1000 bits for the exponent */
      while ((!mpfr_overflow_p()) && (max>0))
        {
          mpfr_mul(a, a, a, MPFR_RNDN);
          if (!mpfr_check(a)) ERROR("for mul");
          max--;
        }
      if (max==0) ERROR("can't reach overflow");
      mpfr_set_ui(a, 2137, MPFR_RNDN);
      /* Corrupt a and check for it */
      MPFR_SIGN(a) = 2;
      if (mpfr_check(a))  ERROR("sgn");
      MPFR_SET_POS(a);
      /* Check prec */
      MPFR_PREC(a) = MPFR_PREC_MIN - 1;
      if (mpfr_check(a))  ERROR("precmin");
#if MPFR_VERSION_MAJOR < 3
      /* Disable the test with MPFR >= 3 since mpfr_prec_t is now signed.
         The "if" below is sufficient, but the MPFR_PREC_MAX+1 generates
         a warning with GCC 4.4.4 even though the test is always false. */
      if ((mpfr_prec_t) 0 - 1 > 0)
        {
          MPFR_PREC(a) = MPFR_PREC_MAX+1;
          if (mpfr_check(a))  ERROR("precmax");
        }
#endif
      MPFR_PREC(a) = pr;
      if (!mpfr_check(a)) ERROR("prec");
      /* Check exponent */
      MPFR_EXP(a) = MPFR_EXP_INVALID;
      if (mpfr_check(a))  ERROR("exp invalid");
      MPFR_EXP(a) = -MPFR_EXP_INVALID;
      if (mpfr_check(a))  ERROR("-exp invalid");
      MPFR_EXP(a) = 0;
      if (!mpfr_check(a)) ERROR("exp 0");
      /* Check Mantissa */
      p = MPFR_MANT(a);
      MPFR_MANT(a) = NULL;
      if (mpfr_check(a))  ERROR("Mantissa Null Ptr");
      MPFR_MANT(a) = p;
      /* Check size */
      s = MPFR_GET_ALLOC_SIZE(a);
      MPFR_SET_ALLOC_SIZE(a, 0);
      if (mpfr_check(a))  ERROR("0 size");
      MPFR_SET_ALLOC_SIZE(a, MP_SIZE_T_MIN);
      if (mpfr_check(a))  ERROR("min size");
      MPFR_SET_ALLOC_SIZE(a, MPFR_LIMB_SIZE(a)-1 );
      if (mpfr_check(a))  ERROR("size < prec");
      MPFR_SET_ALLOC_SIZE(a, s);
      /* Check normal form */
      tmp = MPFR_MANT(a)[0];
      if ((pr % GMP_NUMB_BITS) != 0)
        {
          MPFR_MANT(a)[0] = MPFR_LIMB_MAX;
          if (mpfr_check(a))  ERROR("last bits non 0");
        }
      MPFR_MANT(a)[0] = tmp;
      MPFR_MANT(a)[MPFR_LIMB_SIZE(a)-1] &= MPFR_LIMB_MASK (GMP_NUMB_BITS-1);
      if (mpfr_check(a))  ERROR("last bits non 0");
      /* Final */
      mpfr_set_ui(a, 2137, MPFR_RNDN);
      if (!mpfr_check(a)) ERROR("after last set");
      mpfr_clear (a);
      if (mpfr_check(a))  ERROR("after clear");
    }
  tests_end_mpfr ();
  return 0;
}

int
mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
{
  int sign;
  int rnd_away;
  mp_exp_t exp;

  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ))
    {
      if (MPFR_IS_NAN(u))
        {
          MPFR_SET_NAN(r);
          MPFR_RET_NAN;
        }
      MPFR_SET_SAME_SIGN(r, u);
      if (MPFR_IS_INF(u))
        {
          MPFR_SET_INF(r);
          MPFR_RET(0);  /* infinity is exact */
        }
      else /* now u is zero */
        {
          MPFR_ASSERTD(MPFR_IS_ZERO(u));
          MPFR_SET_ZERO(r);
          MPFR_RET(0);  /* zero is exact */
        }
    }
  MPFR_SET_SAME_SIGN (r, u); /* Does nothing if r==u */

  sign = MPFR_INT_SIGN (u);
  exp = MPFR_GET_EXP (u);

  rnd_away =
    rnd_mode == GMP_RNDD ? sign < 0 :
    rnd_mode == GMP_RNDU ? sign > 0 :
    rnd_mode == GMP_RNDZ ? 0 : -1;

  /* rnd_away:
     1 if round away from zero,
     0 if round to zero,
     -1 if not decided yet.
   */

  if (MPFR_UNLIKELY (exp <= 0))  /* 0 < |u| < 1 ==> round |u| to 0 or 1 */
    {
      /* Note: in the GMP_RNDN mode, 0.5 must be rounded to 0. */
      if (rnd_away != 0 &&
          (rnd_away > 0 ||
           (exp == 0 && (rnd_mode == GMP_RNDNA ||
                         !mpfr_powerof2_raw (u)))))
        {
          mp_limb_t *rp;
          mp_size_t rm;

          rp = MPFR_MANT(r);
          rm = (MPFR_PREC(r) - 1) / BITS_PER_MP_LIMB;
          rp[rm] = MPFR_LIMB_HIGHBIT;
          MPN_ZERO(rp, rm);
          MPFR_SET_EXP (r, 1);  /* |r| = 1 */
          MPFR_RET(sign > 0 ? 2 : -2);
        }
      else
        {
          MPFR_SET_ZERO(r);  /* r = 0 */
          MPFR_RET(sign > 0 ? -2 : 2);
        }
    }
  else  /* exp > 0, |u| >= 1 */
    {
      mp_limb_t *up, *rp;
      mp_size_t un, rn, ui;
      int sh, idiff;
      int uflags;

      /*
       * uflags will contain:
       *   _ 0 if u is an integer representable in r,
       *   _ 1 if u is an integer not representable in r,
       *   _ 2 if u is not an integer.
       */

      up = MPFR_MANT(u);
      rp = MPFR_MANT(r);

      un = MPFR_LIMB_SIZE(u);
      rn = MPFR_LIMB_SIZE(r);
      MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (r));

      MPFR_SET_EXP (r, exp); /* Does nothing if r==u */

      if ((exp - 1) / BITS_PER_MP_LIMB >= un)
        {
          ui = un;
          idiff = 0;
          uflags = 0;  /* u is an integer, representable or not in r */
        }
      else
        {
          mp_size_t uj;

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

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

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

  MPFR_ASSERTN (iop != fop);

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

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

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

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

  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)))
    {
      if (MPFR_IS_NAN(x))
        {
          MPFR_SET_NAN(y);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF(x)) /* u/Inf = 0 */
        {
          MPFR_SET_ZERO(y);
          MPFR_SET_SAME_SIGN(y,x);
          MPFR_RET(0);
        }
      else /* u / 0 */
        {
          MPFR_ASSERTD(MPFR_IS_ZERO(x));
          if (u)
            {
              /* u > 0, so y = sign(x) * Inf */
              MPFR_SET_SAME_SIGN(y, x);
              MPFR_SET_INF(y);
              MPFR_SET_DIVBY0 ();
              MPFR_RET(0);
            }
          else
            {
              /* 0 / 0 */
              MPFR_SET_NAN(y);
              MPFR_RET_NAN;
            }
        }
    }
  else if (MPFR_LIKELY(u != 0))
    {
      mpfr_t uu;
      mp_limb_t up[1];
      int cnt;
      int inex;

      MPFR_SAVE_EXPO_DECL (expo);

      MPFR_TMP_INIT1(up, uu, GMP_NUMB_BITS);
      MPFR_ASSERTN(u == (mp_limb_t) u);
      count_leading_zeros(cnt, (mp_limb_t) u);
      up[0] = (mp_limb_t) u << cnt;

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

int
mpfr_atan2 (mpfr_ptr dest, mpfr_srcptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t tmp, pi;
  int inexact;
  mpfr_prec_t prec;
  mpfr_exp_t e;
  MPFR_SAVE_EXPO_DECL (expo);
  MPFR_ZIV_DECL (loop);

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

  /* Special cases */
  if (MPFR_ARE_SINGULAR (x, y))
    {
      /* atan2(0, 0) does not raise the "invalid" floating-point
         exception, nor does atan2(y, 0) raise the "divide-by-zero"
         floating-point exception.
         -- atan2(±0, -0) returns ±pi.313)
         -- atan2(±0, +0) returns ±0.
         -- atan2(±0, x) returns ±pi, for x < 0.
         -- atan2(±0, x) returns ±0, for x > 0.
         -- atan2(y, ±0) returns -pi/2 for y < 0.
         -- atan2(y, ±0) returns pi/2 for y > 0.
         -- atan2(±oo, -oo) returns ±3pi/4.
         -- atan2(±oo, +oo) returns ±pi/4.
         -- atan2(±oo, x) returns ±pi/2, for finite x.
         -- atan2(±y, -oo) returns ±pi, for finite y > 0.
         -- atan2(±y, +oo) returns ±0, for finite y > 0.
      */
      if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y))
        {
          MPFR_SET_NAN (dest);
          MPFR_RET_NAN;
        }
      if (MPFR_IS_ZERO (y))
        {
          if (MPFR_IS_NEG (x)) /* +/- PI */
            {
            set_pi:
              if (MPFR_IS_NEG (y))
                {
                  inexact =  mpfr_const_pi (dest, MPFR_INVERT_RND (rnd_mode));
                  MPFR_CHANGE_SIGN (dest);
                  return -inexact;
                }
              else
                return mpfr_const_pi (dest, rnd_mode);
            }
          else /* +/- 0 */
            {
            set_zero:
              MPFR_SET_ZERO (dest);
              MPFR_SET_SAME_SIGN (dest, y);
              return 0;
            }
        }
      if (MPFR_IS_ZERO (x))
        {
          return pi_div_2ui (dest, 1, MPFR_IS_NEG (y), rnd_mode);
        }
      if (MPFR_IS_INF (y))
        {
          if (!MPFR_IS_INF (x)) /* +/- PI/2 */
            return pi_div_2ui (dest, 1, MPFR_IS_NEG (y), rnd_mode);
          else if (MPFR_IS_POS (x)) /* +/- PI/4 */
            return pi_div_2ui (dest, 2, MPFR_IS_NEG (y), rnd_mode);
          else /* +/- 3*PI/4: Ugly since we have to round properly */
            {
              mpfr_t tmp2;
              MPFR_ZIV_DECL (loop2);
              mpfr_prec_t prec2 = MPFR_PREC (dest) + 10;

              MPFR_SAVE_EXPO_MARK (expo);
              mpfr_init2 (tmp2, prec2);
              MPFR_ZIV_INIT (loop2, prec2);
              for (;;)
                {
                  mpfr_const_pi (tmp2, MPFR_RNDN);
                  mpfr_mul_ui (tmp2, tmp2, 3, MPFR_RNDN); /* Error <= 2  */
                  mpfr_div_2ui (tmp2, tmp2, 2, MPFR_RNDN);
                  if (mpfr_round_p (MPFR_MANT (tmp2), MPFR_LIMB_SIZE (tmp2),
                                    MPFR_PREC (tmp2) - 2,
                                    MPFR_PREC (dest) + (rnd_mode == MPFR_RNDN)))
                    break;
                  MPFR_ZIV_NEXT (loop2, prec2);
                  mpfr_set_prec (tmp2, prec2);
                }
              MPFR_ZIV_FREE (loop2);
              if (MPFR_IS_NEG (y))
                MPFR_CHANGE_SIGN (tmp2);
              inexact = mpfr_set (dest, tmp2, rnd_mode);
              mpfr_clear (tmp2);
              MPFR_SAVE_EXPO_FREE (expo);
              return mpfr_check_range (dest, inexact, rnd_mode);
//.........这里部分代码省略.........

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

int
mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode)
{
  mpfr_t xp;
  int compared, inexact;
  mp_prec_t prec;
  mp_exp_t xp_exp;
  MPFR_SAVE_EXPO_DECL (expo);
  MPFR_ZIV_DECL (loop);

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

  /* Special cases */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
        {
          MPFR_SET_NAN (asin);
          MPFR_RET_NAN;
        }
      else /* x = 0 */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
          MPFR_SET_ZERO (asin);
          MPFR_SET_SAME_SIGN (asin, x);
          MPFR_RET (0); /* exact result */
        }
    }

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

  /* Set x_p=|x| (x is a normal number) */
  mpfr_init2 (xp, MPFR_PREC (x));
  inexact = mpfr_abs (xp, x, GMP_RNDN);
  MPFR_ASSERTD (inexact == 0);

  compared = mpfr_cmp_ui (xp, 1);

  if (MPFR_UNLIKELY (compared >= 0))
    {
      mpfr_clear (xp);
      if (compared > 0)                  /* asin(x) = NaN for |x| > 1 */
        {
          MPFR_SET_NAN (asin);
          MPFR_RET_NAN;
        }
      else                              /* x = 1 or x = -1 */
        {
          if (MPFR_IS_POS (x)) /* asin(+1) = Pi/2 */
            inexact = mpfr_const_pi (asin, rnd_mode);
          else /* asin(-1) = -Pi/2 */
            {
              inexact = -mpfr_const_pi (asin, MPFR_INVERT_RND(rnd_mode));
              MPFR_CHANGE_SIGN (asin);
            }
          mpfr_div_2ui (asin, asin, 1, rnd_mode); /* May underflow */
          return inexact;
        }
    }

  MPFR_SAVE_EXPO_MARK (expo);

  /* Compute exponent of 1 - ABS(x) */
  mpfr_ui_sub (xp, 1, xp, GMP_RNDD);
  MPFR_ASSERTD (MPFR_GET_EXP (xp) <= 0);
  MPFR_ASSERTD (MPFR_GET_EXP (x) <= 0);
  xp_exp = 2 - MPFR_GET_EXP (xp);

  /* Set up initial prec */
  prec = MPFR_PREC (asin) + 10 + xp_exp;

  /* use asin(x) = atan(x/sqrt(1-x^2)) */
  MPFR_ZIV_INIT (loop, prec);
  for (;;)
    {
      mpfr_set_prec (xp, prec);
      mpfr_sqr (xp, x, GMP_RNDN);
      mpfr_ui_sub (xp, 1, xp, GMP_RNDN);
      mpfr_sqrt (xp, xp, GMP_RNDN);
      mpfr_div (xp, x, xp, GMP_RNDN);
      mpfr_atan (xp, xp, GMP_RNDN);
      if (MPFR_LIKELY (MPFR_CAN_ROUND (xp, prec - xp_exp,
                                       MPFR_PREC (asin), rnd_mode)))
        break;
      MPFR_ZIV_NEXT (loop, prec);
    }
  MPFR_ZIV_FREE (loop);
  inexact = mpfr_set (asin, xp, rnd_mode);

  mpfr_clear (xp);

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

/* (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&q