C++ MPFR_CAN_ROUND函数代码示例

static void
tgeneric_cci (mpc_function *function, mpc_ptr op1, int op2,
              mpc_ptr rop, mpc_ptr rop4, mpc_ptr rop4rnd, mpc_rnd_t rnd)
{
  known_signs_t ks = {1, 1};

  function->pointer.CCI (rop4, op1, op2, rnd);
  function->pointer.CCI (rop, op1, op2, rnd);
  if (MPFR_CAN_ROUND (mpc_realref (rop4), 1, MPC_PREC_RE (rop),
                      MPC_RND_RE (rnd))
      && MPFR_CAN_ROUND (mpc_imagref (rop4), 1, MPC_PREC_IM (rop),
                         MPC_RND_IM (rnd)))
    mpc_set (rop4rnd, rop4, rnd);
  else
    return;

  if (same_mpc_value (rop, rop4rnd, ks))
    return;

  printf ("Rounding in %s might be incorrect for\n", function->name);
  MPC_OUT (op1);
  printf ("op2=%d\n", op2);
  printf ("with rounding mode (%s, %s)",
          mpfr_print_rnd_mode (MPC_RND_RE (rnd)),
          mpfr_print_rnd_mode (MPC_RND_IM (rnd)));

  printf ("\n%s                     gives ", function->name);
  MPC_OUT (rop);
  printf ("%s quadruple precision gives ", function->name);
  MPC_OUT (rop4);
  printf ("and is rounded to                  ");
  MPC_OUT (rop4rnd);

  exit (1);
}

static void
tgeneric_cc_c (mpc_function *function, mpc_ptr op, mpc_ptr rop1, mpc_ptr rop2,
   mpc_ptr rop14, mpc_ptr rop24, mpc_ptr rop14rnd, mpc_ptr rop24rnd,
   mpc_rnd_t rnd1, mpc_rnd_t rnd2)
{
   /* same as the previous function, but for mpc functions computing two
      results from one argument                                          */
   known_signs_t ks = {1, 1};

   function->pointer.CC_C (rop14, rop24, op, rnd1, rnd2);
   function->pointer.CC_C (rop1,  rop2,  op, rnd1, rnd2);

   if (   MPFR_CAN_ROUND (mpc_realref (rop14), 1, MPC_PREC_RE (rop1),
                          MPC_RND_RE (rnd1))
       && MPFR_CAN_ROUND (mpc_imagref (rop14), 1, MPC_PREC_IM (rop1),
                          MPC_RND_IM (rnd1))
       && MPFR_CAN_ROUND (mpc_realref (rop24), 1, MPC_PREC_RE (rop2),
                          MPC_RND_RE (rnd2))
       && MPFR_CAN_ROUND (mpc_imagref (rop24), 1, MPC_PREC_IM (rop2),
                          MPC_RND_IM (rnd2))) {
     mpc_set (rop14rnd, rop14, rnd1);
     mpc_set (rop24rnd, rop24, rnd2);
   }
   else
     return;

   if (!same_mpc_value (rop1, rop14rnd, ks)) {
      /* rounding failed for first result */
      printf ("Rounding might be incorrect for the first result of %s at\n", function->name);
      MPC_OUT (op);
      printf ("with rounding mode (%s, %s)",
          mpfr_print_rnd_mode (MPC_RND_RE (rnd1)),
          mpfr_print_rnd_mode (MPC_RND_IM (rnd1)));
      printf ("\n%s                     gives ", function->name);
      MPC_OUT (rop1);
      printf ("%s quadruple precision gives ", function->name);
      MPC_OUT (rop14);
      printf ("and is rounded to                  ");
      MPC_OUT (rop14rnd);
      exit (1);
   }
   else if (!same_mpc_value (rop2, rop24rnd, ks)) {
      /* rounding failed for second result */
      printf ("Rounding might be incorrect for the second result of %s at\n", function->name);
      MPC_OUT (op);
      printf ("with rounding mode (%s, %s)",
          mpfr_print_rnd_mode (MPC_RND_RE (rnd2)),
          mpfr_print_rnd_mode (MPC_RND_IM (rnd2)));
      printf ("\n%s                     gives ", function->name);
      MPC_OUT (rop2);
      printf ("%s quadruple precision gives ", function->name);
      MPC_OUT (rop24);
      printf ("and is rounded to                  ");
      MPC_OUT (rop24rnd);
      exit (1);
   }
}

/* try asymptotic expansion when x is large and negative:
   Li2(x) = -log(-x)^2/2 - Pi^2/6 - 1/x + O(1/x^2).
   More precisely for x <= -2 we have for g(x) = -log(-x)^2/2 - Pi^2/6:
   |Li2(x) - g(x)| <= 1/|x|.
   Assumes x <= -7, which ensures |log(-x)^2/2| >= Pi^2/6, and g(x) <= -3.5.
   Return 0 if asymptotic expansion failed (unable to round), otherwise
   returns correct ternary value.
*/
static int
mpfr_li2_asympt_neg (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t g, h;
  mp_prec_t w = MPFR_PREC (y) + 20;
  int inex = 0;

  MPFR_ASSERTN (mpfr_cmp_si (x, -7) <= 0);

  mpfr_init2 (g, w);
  mpfr_init2 (h, w);
  mpfr_neg (g, x, GMP_RNDN);
  mpfr_log (g, g, GMP_RNDN);    /* rel. error <= |(1 + theta) - 1| */
  mpfr_sqr (g, g, GMP_RNDN);    /* rel. error <= |(1 + theta)^3 - 1| <= 2^(2-w) */
  mpfr_div_2ui (g, g, 1, GMP_RNDN);     /* rel. error <= 2^(2-w) */
  mpfr_const_pi (h, GMP_RNDN);  /* error <= 2^(1-w) */
  mpfr_sqr (h, h, GMP_RNDN);    /* rel. error <= 2^(2-w) */
  mpfr_div_ui (h, h, 6, GMP_RNDN);      /* rel. error <= |(1 + theta)^4 - 1|
                                           <= 5 * 2^(-w) */
  MPFR_ASSERTN (MPFR_EXP (g) >= MPFR_EXP (h));
  mpfr_add (g, g, h, GMP_RNDN); /* err <= ulp(g)/2 + g*2^(2-w) + g*5*2^(-w)
                                   <= ulp(g) * (1/2 + 4 + 5) < 10 ulp(g).

                                   If in addition |1/x| <= 4 ulp(g), then the
                                   total error is bounded by 16 ulp(g). */
  if ((MPFR_EXP (x) >= (mp_exp_t) (w - 2) - MPFR_EXP (g)) &&
      MPFR_CAN_ROUND (g, w - 4, MPFR_PREC (y), rnd_mode))
    inex = mpfr_neg (y, g, rnd_mode);

  mpfr_clear (g);
  mpfr_clear (h);

  return inex;
}

static void
tgeneric_fc (mpc_function *function, mpc_ptr op, mpfr_ptr rop,
             mpfr_ptr rop4, mpfr_ptr rop4rnd, mpfr_rnd_t rnd)
{
  function->pointer.FC (rop4, op, rnd);
  function->pointer.FC (rop, op, rnd);
  if (MPFR_CAN_ROUND (rop4, 1, mpfr_get_prec (rop), rnd))
    mpfr_set (rop4rnd, rop4, rnd);
  else
    return;

  if (same_mpfr_value (rop, rop4rnd, 1))
    return;

  printf ("Rounding in %s might be incorrect for\n", function->name);
  MPC_OUT (op);
  printf ("with rounding mode %s", mpfr_print_rnd_mode (rnd));

  printf ("\n%s                     gives ", function->name);
  MPFR_OUT (rop);
  printf ("%s quadruple precision gives ", function->name);
  MPFR_OUT (rop4);
  printf ("and is rounded to                  ");
  MPFR_OUT (rop4rnd);

  exit (1);
}

/* functions with one input, one output */
static void
tgeneric_cc (mpc_function *function, mpc_ptr op, mpc_ptr rop,
             mpc_ptr rop4, mpc_ptr rop4rnd, mpc_rnd_t rnd)
{
  known_signs_t ks = {1, 1};

  /* We compute the result with four times the precision and check whether the
     rounding is correct. Error reports in this part of the algorithm might
     still be wrong, though, since there are two consecutive roundings (but we
     try to avoid them).  */
  function->pointer.CC (rop4, op, rnd);
  function->pointer.CC (rop, op, rnd);

  /* can't use the mpfr_can_round function when argument is singular,
     use a custom macro instead. */
  if (MPFR_CAN_ROUND (mpc_realref (rop4), 1, MPC_PREC_RE (rop),
                      MPC_RND_RE (rnd))
      && MPFR_CAN_ROUND (mpc_imagref (rop4), 1, MPC_PREC_IM (rop),
                         MPC_RND_IM (rnd)))
    mpc_set (rop4rnd, rop4, rnd);
  else
    /* avoid double rounding error */
    return;

  if (same_mpc_value (rop, rop4rnd, ks))
    return;

  /* rounding failed */
  printf ("Rounding in %s might be incorrect for\n", function->name);
  MPC_OUT (op);

  printf ("with rounding mode (%s, %s)",
          mpfr_print_rnd_mode (MPC_RND_RE (rnd)),
          mpfr_print_rnd_mode (MPC_RND_IM (rnd)));

  printf ("\n%s                     gives ", function->name);
  MPC_OUT (rop);
  printf ("%s quadruple precision gives ", function->name);
  MPC_OUT (rop4);
  printf ("and is rounded to                  ");
  MPC_OUT (rop4rnd);

  exit (1);
}

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_atanh (mpfr_ptr y, mpfr_srcptr xt , mpfr_rnd_t rnd_mode)
{
  int inexact;
  mpfr_t x, t, te;
  mpfr_prec_t Nx, Ny, Nt;
  mpfr_exp_t err;
  MPFR_ZIV_DECL (loop);
  MPFR_SAVE_EXPO_DECL (expo);

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

  /* Special cases */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
    {
      /* atanh(NaN) = NaN, and atanh(+/-Inf) = NaN since tanh gives a result
         between -1 and 1 */
      if (MPFR_IS_NAN (xt) || MPFR_IS_INF (xt))
        {
          MPFR_SET_NAN (y);
          MPFR_RET_NAN;
        }
      else /* necessarily xt is 0 */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (xt));
          MPFR_SET_ZERO (y);   /* atanh(0) = 0 */
          MPFR_SET_SAME_SIGN (y,xt);
          MPFR_RET (0);
        }
    }

  /* atanh (x) = NaN as soon as |x| > 1, and arctanh(+/-1) = +/-Inf */
  if (MPFR_UNLIKELY (MPFR_GET_EXP (xt) > 0))
    {
      if (MPFR_GET_EXP (xt) == 1 && mpfr_powerof2_raw (xt))
        {
          MPFR_SET_INF (y);
          MPFR_SET_SAME_SIGN (y, xt);
          mpfr_set_divby0 ();
          MPFR_RET (0);
        }
      MPFR_SET_NAN (y);
      MPFR_RET_NAN;
    }

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

  MPFR_SAVE_EXPO_MARK (expo);

  /* Compute initial precision */
  Nx = MPFR_PREC (xt);
  MPFR_TMP_INIT_ABS (x, xt);
  Ny = MPFR_PREC (y);
  Nt = MAX (Nx, Ny);
  /* the optimal number of bits : see algorithms.ps */
  Nt = Nt + MPFR_INT_CEIL_LOG2 (Nt) + 4;

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

  /* First computation of cosh */
  MPFR_ZIV_INIT (loop, Nt);
  for (;;)
    {
      /* compute atanh */
      mpfr_ui_sub (te, 1, x, MPFR_RNDU);   /* (1-xt)*/
      mpfr_add_ui (t,  x, 1, MPFR_RNDD);   /* (xt+1)*/
      mpfr_div (t, t, te, MPFR_RNDN);      /* (1+xt)/(1-xt)*/
      mpfr_log (t, t, MPFR_RNDN);          /* ln((1+xt)/(1-xt))*/
      mpfr_div_2ui (t, t, 1, MPFR_RNDN);   /* (1/2)*ln((1+xt)/(1-xt))*/

      /* error estimate: see algorithms.tex */
      /* FIXME: this does not correspond to the value in algorithms.tex!!! */
      /* err=Nt-__gmpfr_ceil_log2(1+5*pow(2,1-MPFR_EXP(t)));*/
      err = Nt - (MAX (4 - MPFR_GET_EXP (t), 0) + 1);

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

      /* reactualisation of the precision */
      MPFR_ZIV_NEXT (loop, Nt);
      mpfr_set_prec (t, Nt);
      mpfr_set_prec (te, Nt);
    }
  MPFR_ZIV_FREE (loop);

  inexact = mpfr_set4 (y, t, rnd_mode, MPFR_SIGN (xt));

  mpfr_clear(t);
  mpfr_clear(te);

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (y, inexact, rnd_mode);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            MPFR_SET_INF (y);
          else
            MPFR_SET_ZERO (y);
          MPFR_SET_POS (y);
          MPFR_RET (0);
        }
      else /* 2^0 = 1 */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO(x));
          return mpfr_set_ui (y, 1, rnd_mode);
        }
    }

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

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

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

  MPFR_SAVE_EXPO_MARK (expo);

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

  xint = mpfr_get_si (x, MPFR_RNDZ);
  mpfr_init2 (xfrac, MPFR_PREC (x));
  mpfr_sub_si (xfrac, x, xint, MPFR_RNDN); /* exact */

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

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

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

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

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

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

          /* Actualisation of the precision */
          MPFR_ZIV_NEXT (loop, Nt);
          mpfr_set_prec (t, Nt);
        }
      MPFR_ZIV_FREE (loop);

      inexact = mpfr_set (y, t, rnd_mode);

      mpfr_clear (t);
    }

  mpfr_clear (xfrac);
  MPFR_CLEAR_FLAGS ();
  mpfr_mul_2si (y, y, xint, MPFR_RNDN); /* exact or overflow */
  /* Note: We can have an overflow only when t was rounded up to 2. */
  MPFR_ASSERTD (MPFR_IS_PURE_FP (y) || inexact > 0);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (y, inexact, rnd_mode);
}

//.........这里部分代码省略.........
#ifdef MPFR_JN
          mpfr_mul (c, c, Q, MPFR_RNDN); /* Q * (sin - cos) */
          mpfr_mul (s, s, P, MPFR_RNDN); /* P * (sin + cos) */
#else
          mpfr_mul (c, c, P, MPFR_RNDN); /* P * (sin - cos) */
          mpfr_mul (s, s, Q, MPFR_RNDN); /* Q * (sin + cos) */
#endif
          err = MPFR_EXP(c);
          if (MPFR_EXP(s) > err)
            err = MPFR_EXP(s);
#ifdef MPFR_JN
          mpfr_sub (s, s, c, MPFR_RNDN);
#else
          mpfr_add (s, s, c, MPFR_RNDN);
#endif
        }
      else /* n odd: P * (sin - cos) + Q (cos + sin) for jn,
                     Q * (sin - cos) - P (cos + sin) for yn */
        {
#ifdef MPFR_JN
          mpfr_mul (c, c, P, MPFR_RNDN); /* P * (sin - cos) */
          mpfr_mul (s, s, Q, MPFR_RNDN); /* Q * (sin + cos) */
#else
          mpfr_mul (c, c, Q, MPFR_RNDN); /* Q * (sin - cos) */
          mpfr_mul (s, s, P, MPFR_RNDN); /* P * (sin + cos) */
#endif
          err = MPFR_EXP(c);
          if (MPFR_EXP(s) > err)
            err = MPFR_EXP(s);
#ifdef MPFR_JN
          mpfr_add (s, s, c, MPFR_RNDN);
#else
          mpfr_sub (s, c, s, MPFR_RNDN);
#endif
        }
      if ((n & 2) != 0)
        mpfr_neg (s, s, MPFR_RNDN);
      if (MPFR_EXP(s) > err)
        err = MPFR_EXP(s);
      /* the absolute error on s is bounded by P*err(s/c) + Q*err(s/c)
         + err(P)*(s/c) + err(Q)*(s/c) + 3 * 2^(err - w - 1)
         <= (|P|+|Q|) * 2^(1-w) + err_s * 2^(1-w) + 2^err * 2^(1-w),
         since |c|, |old_s| <= 2. */
      err2 = (MPFR_EXP(P) >= MPFR_EXP(Q)) ? MPFR_EXP(P) + 2 : MPFR_EXP(Q) + 2;
      /* (|P| + |Q|) * 2^(1 - w) <= 2^(err2 - w) */
      err = MPFR_EXP(err_s) >= err ? MPFR_EXP(err_s) + 2 : err + 2;
      /* err_s * 2^(1-w) + 2^old_err * 2^(1-w) <= 2^err * 2^(-w) */
      err2 = (err >= err2) ? err + 1 : err2 + 1;
      /* now the absolute error on s is bounded by 2^(err2 - w) */

      /* multiply by sqrt(1/(Pi*z)) */
      mpfr_const_pi (c, MPFR_RNDN);     /* Pi, err <= 1 */
      mpfr_mul (c, c, z, MPFR_RNDN);    /* err <= 2 */
      mpfr_si_div (c, MPFR_IS_POS(z) ? 1 : -1, c, MPFR_RNDN); /* err <= 3 */
      mpfr_sqrt (c, c, MPFR_RNDN);      /* err<=5/2, thus the absolute error is
                                          bounded by 3*u*|c| for |u| <= 0.25 */
      mpfr_mul (err_t, c, s, MPFR_SIGN(c)==MPFR_SIGN(s) ? MPFR_RNDU : MPFR_RNDD);
      mpfr_abs (err_t, err_t, MPFR_RNDU);
      mpfr_mul_ui (err_t, err_t, 3, MPFR_RNDU);
      /* 3*2^(-w)*|old_c|*|s| [see below] is bounded by err_t * 2^(-w) */
      err2 += MPFR_EXP(c);
      /* |old_c| * 2^(err2 - w) [see below] is bounded by 2^(err2-w) */
      mpfr_mul (c, c, s, MPFR_RNDN);    /* the absolute error on c is bounded by
                                          1/2 ulp(c) + 3*2^(-w)*|old_c|*|s|
                                          + |old_c| * 2^(err2 - w) */
      /* compute err_t * 2^(-w) + 1/2 ulp(c) = (err_t + 2^EXP(c)) * 2^(-w) */
      err = (MPFR_EXP(err_t) > MPFR_EXP(c)) ? MPFR_EXP(err_t) + 1 : MPFR_EXP(c) + 1;
      /* err_t * 2^(-w) + 1/2 ulp(c) <= 2^(err - w) */
      /* now err_t * 2^(-w) bounds 1/2 ulp(c) + 3*2^(-w)*|old_c|*|s| */
      err = (err >= err2) ? err + 1 : err2 + 1;
      /* the absolute error on c is bounded by 2^(err - w) */

      mpfr_clear (s);
      mpfr_clear (P);
      mpfr_clear (Q);
      mpfr_clear (t);
      mpfr_clear (iz);
      mpfr_clear (err_t);
      mpfr_clear (err_s);
      mpfr_clear (err_u);

      err -= MPFR_EXP(c);
      if (MPFR_LIKELY (MPFR_CAN_ROUND (c, w - err, MPFR_PREC(res), r)))
        break;
      if (diverge != 0)
        {
          mpfr_set (c, z, r); /* will force inex=0 below, which means the
                               asymptotic expansion failed */
          break;
        }
      MPFR_ZIV_NEXT (loop, w);
    }
  MPFR_ZIV_FREE (loop);

  inex = (MPFR_IS_POS(z) || ((n & 1) == 0)) ? mpfr_set (res, c, r)
    : mpfr_neg (res, c, r);
  mpfr_clear (c);

  return inex;
}

//.........这里部分代码省略.........
            {
              MPFR_SET_NAN (y);
              MPFR_RET_NAN;
            }
        }
      else /* x is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
          MPFR_SET_ZERO (y);   /* log1p(+/- 0) = +/- 0 */
          MPFR_SET_SAME_SIGN (y, x);
          MPFR_RET (0);
        }
    }

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

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

  MPFR_SAVE_EXPO_MARK (expo);

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

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

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

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

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

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

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

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

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

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

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

//.........这里部分代码省略.........
              mpfr_log (t, absx, MPFR_IS_NEG (y) ? MPFR_RNDU : MPFR_RNDD);
              mpfr_mul (t, y, t, MPFR_RNDD);            /* y * ln|x| */
              MPFR_BLOCK (flags2, mpfr_exp (t, t, MPFR_RNDD));
              /* t = lower bound on exp(y * ln|x|) */
              if (MPFR_OVERFLOW (flags2))
                {
                  /* We have computed a lower bound on |x|^y, and it
                     overflowed. Therefore we have a real overflow
                     on |x|^y. */
                  inexact = mpfr_overflow (z, rnd_mode, MPFR_SIGN_POS);
                  if (expo != NULL)
                    MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, MPFR_FLAGS_INEXACT
                                                 | MPFR_FLAGS_OVERFLOW);
                  break;
                }
            }

          k_non_zero = 1;
          Ntmin = sizeof(mpfr_exp_t) * CHAR_BIT;
          if (Ntmin > Nt)
            {
              Nt = Ntmin;
              mpfr_set_prec (t, Nt);
            }
          mpfr_init2 (u, Nt);
          mpfr_init2 (k, Ntmin);
          mpfr_log2 (k, absx, MPFR_RNDN);
          mpfr_mul (k, y, k, MPFR_RNDN);
          mpfr_round (k, k);
          MPFR_LOG_VAR (k);
          /* |y| < 2^Ntmin, therefore |k| < 2^Nt. */
          continue;
        }
      if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - err, Nz, rnd_mode)))
        {
          inexact = mpfr_set (z, t, rnd_mode);
          break;
        }

      /* check exact power, except when y is an integer (since the
         exact cases for y integer have already been filtered out) */
      if (check_exact_case == 0 && ! y_is_integer)
        {
          if (mpfr_pow_is_exact (z, absx, y, rnd_mode, &inexact))
            break;
          check_exact_case = 1;
        }

      /* reactualisation of the precision */
      MPFR_ZIV_NEXT (ziv_loop, Nt);
      mpfr_set_prec (t, Nt);
      if (k_non_zero)
        mpfr_set_prec (u, Nt);
    }
  MPFR_ZIV_FREE (ziv_loop);

  if (k_non_zero)
    {
      int inex2;
      long lk;

      /* The rounded result in an unbounded exponent range is z * 2^k. As
       * MPFR chooses underflow after rounding, the mpfr_mul_2si below will
       * correctly detect underflows and overflows. However, in rounding to
       * nearest, if z * 2^k = 2^(emin - 2), then the double rounding may
       * affect the result. We need to cope with that before overwriting z.

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

  mpfr_init (c);
  mpfr_init (xr);

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

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

      sign = MPFR_SIGN(xx);
      /* now that the argument is reduced, precision m is enough */
      mpfr_set_prec (c, m);
      mpfr_cos (c, xx, MPFR_RNDZ);    /* can't be exact */
      mpfr_nexttoinf (c);           /* now c = cos(x) rounded away */
      mpfr_mul (c, c, c, MPFR_RNDU); /* away */
      mpfr_ui_sub (c, 1, c, MPFR_RNDZ);
      mpfr_sqrt (c, c, MPFR_RNDZ);
      if (MPFR_IS_NEG_SIGN(sign))
        MPFR_CHANGE_SIGN(c);

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

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

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

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

  mpfr_clear (c);
  mpfr_clear (xr);

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

//.........这里部分代码省略.........
              MPFR_RET_NAN;
            }
          else /* log(+Inf) = +Inf */
            {
              MPFR_SET_INF (r);
              MPFR_SET_POS (r);
              MPFR_RET (0);
            }
        }
      else /* a is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (a));
          MPFR_SET_INF (r);
          MPFR_SET_NEG (r);
          mpfr_set_divby0 ();
          MPFR_RET (0); /* log(0) is an exact -infinity */
        }
    }
  /* If a is negative, the result is NaN */
  else if (MPFR_UNLIKELY (MPFR_IS_NEG (a)))
    {
      MPFR_SET_NAN (r);
      MPFR_RET_NAN;
    }
  /* If a is 1, the result is 0 */
  else if (MPFR_UNLIKELY (MPFR_GET_EXP (a) == 1 && mpfr_cmp_ui (a, 1) == 0))
    {
      MPFR_SET_ZERO (r);
      MPFR_SET_POS (r);
      MPFR_RET (0); /* only "normal" case where the result is exact */
    }

  q = MPFR_PREC (r);

  /* use initial precision about q+lg(q)+5 */
  p = q + 5 + 2 * MPFR_INT_CEIL_LOG2 (q);
  /* % ~(mpfr_prec_t)GMP_NUMB_BITS  ;
     m=q; while (m) { p++; m >>= 1; }  */
  /* if (MPFR_LIKELY(p % GMP_NUMB_BITS != 0))
      p += GMP_NUMB_BITS - (p%GMP_NUMB_BITS); */

  MPFR_SAVE_EXPO_MARK (expo);
  MPFR_GROUP_INIT_2 (group, p, tmp1, tmp2);

  MPFR_ZIV_INIT (loop, p);
  for (;;)
    {
      long m;
      mpfr_exp_t cancel;

      /* Calculus of m (depends on p) */
      m = (p + 1) / 2 - MPFR_GET_EXP (a) + 1;

      mpfr_mul_2si (tmp2, a, m, MPFR_RNDN);    /* s=a*2^m,        err<=1 ulp  */
      mpfr_div (tmp1, __gmpfr_four, tmp2, MPFR_RNDN);/* 4/s,      err<=2 ulps */
      mpfr_agm (tmp2, __gmpfr_one, tmp1, MPFR_RNDN); /* AG(1,4/s),err<=3 ulps */
      mpfr_mul_2ui (tmp2, tmp2, 1, MPFR_RNDN); /* 2*AG(1,4/s),    err<=3 ulps */
      mpfr_const_pi (tmp1, MPFR_RNDN);         /* compute pi,     err<=1ulp   */
      mpfr_div (tmp2, tmp1, tmp2, MPFR_RNDN);  /* pi/2*AG(1,4/s), err<=5ulps  */
      mpfr_const_log2 (tmp1, MPFR_RNDN);      /* compute log(2),  err<=1ulp   */
      mpfr_mul_si (tmp1, tmp1, m, MPFR_RNDN); /* compute m*log(2),err<=2ulps  */
      mpfr_sub (tmp1, tmp2, tmp1, MPFR_RNDN); /* log(a),    err<=7ulps+cancel */

      if (MPFR_LIKELY (MPFR_IS_PURE_FP (tmp1) && MPFR_IS_PURE_FP (tmp2)))
        {
          cancel = MPFR_GET_EXP (tmp2) - MPFR_GET_EXP (tmp1);
          MPFR_LOG_MSG (("canceled bits=%ld\n", (long) cancel));
          MPFR_LOG_VAR (tmp1);
          if (MPFR_UNLIKELY (cancel < 0))
            cancel = 0;

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

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

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

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

int
mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t xp;
  int compared, inexact;
  mpfr_prec_t prec;
  mpfr_exp_t xp_exp;
  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),
    ("asin[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (asin), mpfr_log_prec, 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, MPFR_RNDN);
  MPFR_ASSERTD (inexact == 0);

  compared = mpfr_cmp_ui (xp, 1);

  MPFR_SAVE_EXPO_MARK (expo);

  if (MPFR_UNLIKELY (compared >= 0))
    {
      mpfr_clear (xp);
      if (compared > 0)                  /* asin(x) = NaN for |x| > 1 */
        {
          MPFR_SAVE_EXPO_FREE (expo);
          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);
        }
    }
  else
    {
      /* Compute exponent of 1 - ABS(x) */
      mpfr_ui_sub (xp, 1, xp, MPFR_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, MPFR_RNDN);
          mpfr_ui_sub (xp, 1, xp, MPFR_RNDN);
          mpfr_sqrt (xp, xp, MPFR_RNDN);
          mpfr_div (xp, x, xp, MPFR_RNDN);
          mpfr_atan (xp, xp, MPFR_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);
//.........这里部分代码省略.........