C++ MPFR_IS_NEG函数代码示例

static void
check_pow_si (void)
{
    mpfr_t x;

    mpfr_init (x);

    mpfr_set_nan (x);
    mpfr_pow_si (x, x, -1, GMP_RNDN);
    MPFR_ASSERTN(mpfr_nan_p (x));

    mpfr_set_inf (x, 1);
    mpfr_pow_si (x, x, -1, GMP_RNDN);
    MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS(x));

    mpfr_set_inf (x, -1);
    mpfr_pow_si (x, x, -1, GMP_RNDN);
    MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_NEG(x));

    mpfr_set_inf (x, -1);
    mpfr_pow_si (x, x, -2, GMP_RNDN);
    MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS(x));

    mpfr_set_ui (x, 0, GMP_RNDN);
    mpfr_pow_si (x, x, -1, GMP_RNDN);
    MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);

    mpfr_set_ui (x, 0, GMP_RNDN);
    mpfr_neg (x, x, GMP_RNDN);
    mpfr_pow_si (x, x, -1, GMP_RNDN);
    MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) < 0);

    mpfr_set_ui (x, 0, GMP_RNDN);
    mpfr_neg (x, x, GMP_RNDN);
    mpfr_pow_si (x, x, -2, GMP_RNDN);
    MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);

    mpfr_set_si (x, 2, GMP_RNDN);
    mpfr_pow_si (x, x, LONG_MAX, GMP_RNDN);  /* 2^LONG_MAX */
    if (LONG_MAX > mpfr_get_emax () - 1)  /* LONG_MAX + 1 > emax */
    {
        MPFR_ASSERTN (mpfr_inf_p (x));
    }
    else
    {
        MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 1, LONG_MAX));
    }

    mpfr_set_si (x, 2, GMP_RNDN);
    mpfr_pow_si (x, x, LONG_MIN, GMP_RNDN);  /* 2^LONG_MIN */
    if (LONG_MIN + 1 < mpfr_get_emin ())
    {
        MPFR_ASSERTN (mpfr_zero_p (x));
    }
    else
    {
        MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 1, LONG_MIN));
    }

    mpfr_set_si (x, 2, GMP_RNDN);
    mpfr_pow_si (x, x, LONG_MIN + 1, GMP_RNDN);  /* 2^(LONG_MIN+1) */
    if (mpfr_nan_p (x))
    {
        printf ("Error in pow_si(2, LONG_MIN+1): got NaN\n");
        exit (1);
    }
    if (LONG_MIN + 2 < mpfr_get_emin ())
    {
        MPFR_ASSERTN (mpfr_zero_p (x));
    }
    else
    {
        MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 1, LONG_MIN + 1));
    }

    mpfr_set_si_2exp (x, 1, -1, GMP_RNDN);  /* 0.5 */
    mpfr_pow_si (x, x, LONG_MIN, GMP_RNDN);  /* 2^(-LONG_MIN) */
    if (LONG_MIN < 1 - mpfr_get_emax ())  /* 1 - LONG_MIN > emax */
    {
        MPFR_ASSERTN (mpfr_inf_p (x));
    }
    else
    {
        MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 2, - (LONG_MIN + 1)));
    }

    mpfr_clear (x);
}

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

  mpfr_init (c);

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

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

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

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

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

/* (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, mp_rnd_t rnd_mode)
{
  mp_prec_t prec, m;
  int neg, reduce;
  mpfr_t c, xr;
  mpfr_srcptr xx;
  mp_exp_t err, expx;
  MPFR_ZIV_DECL (loop);

  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 */
          return mpfr_set_ui (z, 1, rnd_mode);
        }
    }

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

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

  mpfr_init (c);
  mpfr_init (xr);

  MPFR_ZIV_INIT (loop, m);
  for (;;)
    {
      /* the following is copied from sin.c */
      if (expx >= 2) /* reduce the argument */
        {
          reduce = 1;
          mpfr_set_prec (c, expx + m - 1);
          mpfr_set_prec (xr, m);
          mpfr_const_pi (c, GMP_RNDN);
          mpfr_mul_2ui (c, c, 1, GMP_RNDN);
          mpfr_remainder (xr, x, c, GMP_RNDN);
          mpfr_div_2ui (c, c, 1, GMP_RNDN);
          if (MPFR_SIGN (xr) > 0)
            mpfr_sub (c, c, xr, GMP_RNDZ);
          else
            mpfr_add (c, c, xr, GMP_RNDZ);
          if (MPFR_IS_ZERO(xr) || MPFR_EXP(xr) < (mp_exp_t) 3 - (mp_exp_t) m
              || MPFR_EXP(c) < (mp_exp_t) 3 - (mp_exp_t) m)
            goto next_step;
          xx = xr;
        }
      else /* the input argument is already reduced */
        {
          reduce = 0;
          xx = x;
        }

      neg = MPFR_IS_NEG (xx); /* gives sign of sin(x) */
      mpfr_set_prec (c, m);
      mpfr_cos (c, xx, GMP_RNDZ);
      /* If no argument reduction was performed, the error is at most ulp(c),
         otherwise it is at most ulp(c) + 2^(2-m). Since |c| < 1, we have
         ulp(c) <= 2^(-m), thus the error is bounded by 2^(3-m) in that later
         case. */
      if (reduce == 0)
        err = m;
      else
        err = MPFR_GET_EXP (c) + (mp_exp_t) (m - 3);
      if (!mpfr_can_round (c, err, GMP_RNDN, rnd_mode,
                           MPFR_PREC (z) + (rnd_mode == GMP_RNDN)))
        goto next_step;

      mpfr_set (z, c, rnd_mode);
      mpfr_sqr (c, c, GMP_RNDU);
      mpfr_ui_sub (c, 1, c, GMP_RNDN);
      err = 2 + (- MPFR_GET_EXP (c)) / 2;
      mpfr_sqrt (c, c, GMP_RNDN);
      if (neg)
        MPFR_CHANGE_SIGN (c);

      /* the absolute error on c is at most 2^(err-m), which we must put
         in the form 2^(EXP(c)-err). If there was an argument reduction,
         we need to add 2^(2-m); since err >= 2, the error is bounded by
         2^(err+1-m) in that case. */
      err = MPFR_GET_EXP (c) + (mp_exp_t) m - (err + reduce);
      if (mpfr_can_round (c, err, GMP_RNDN, rnd_mode,
                          MPFR_PREC (y) + (rnd_mode == GMP_RNDN)))
        break;
//.........这里部分代码省略.........

static void
special (void)
{
  mpfr_t x, y;
  int inex;

  mpfr_init (x);
  mpfr_init (y);

  mpfr_set_nan (x);
  mpfr_lngamma (y, x, MPFR_RNDN);
  if (!mpfr_nan_p (y))
    {
      printf ("Error for lngamma(NaN)\n");
      exit (1);
    }

  mpfr_set_inf (x, -1);
  mpfr_lngamma (y, x, MPFR_RNDN);
  if (!mpfr_nan_p (y))
    {
      printf ("Error for lngamma(-Inf)\n");
      exit (1);
    }

  mpfr_set_inf (x, 1);
  mpfr_lngamma (y, x, MPFR_RNDN);
  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
    {
      printf ("Error for lngamma(+Inf)\n");
      exit (1);
    }

  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_lngamma (y, x, MPFR_RNDN);
  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
    {
      printf ("Error for lngamma(+0)\n");
      exit (1);
    }

  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_neg (x, x, MPFR_RNDN);
  mpfr_lngamma (y, x, MPFR_RNDN);
  if (!mpfr_nan_p (y))
    {
      printf ("Error for lngamma(-0)\n");
      exit (1);
    }

  mpfr_set_ui (x, 1, MPFR_RNDN);
  mpfr_lngamma (y, x, MPFR_RNDN);
  if (MPFR_IS_NAN (y) || mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y))
    {
      printf ("Error for lngamma(1)\n");
      exit (1);
    }

  mpfr_set_si (x, -1, MPFR_RNDN);
  mpfr_lngamma (y, x, MPFR_RNDN);
  if (!mpfr_nan_p (y))
    {
      printf ("Error for lngamma(-1)\n");
      exit (1);
    }

  mpfr_set_ui (x, 2, MPFR_RNDN);
  mpfr_lngamma (y, x, MPFR_RNDN);
  if (MPFR_IS_NAN (y) || mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y))
    {
      printf ("Error for lngamma(2)\n");
      exit (1);
    }

  mpfr_set_prec (x, 53);
  mpfr_set_prec (y, 53);

#define CHECK_X1 "1.0762904832837976166"
#define CHECK_Y1 "-0.039418362817587634939"

  mpfr_set_str (x, CHECK_X1, 10, MPFR_RNDN);
  mpfr_lngamma (y, x, MPFR_RNDN);
  mpfr_set_str (x, CHECK_Y1, 10, MPFR_RNDN);
  if (MPFR_IS_NAN (y) || mpfr_cmp (y, x))
    {
      printf ("mpfr_lngamma("CHECK_X1") is wrong:\n"
              "expected ");
      mpfr_print_binary (x); putchar ('\n');
      printf ("got      ");
      mpfr_print_binary (y); putchar ('\n');
      exit (1);
    }

#define CHECK_X2 "9.23709516716202383435e-01"
#define CHECK_Y2 "0.049010669407893718563"
  mpfr_set_str (x, CHECK_X2, 10, MPFR_RNDN);
  mpfr_lngamma (y, x, MPFR_RNDN);
  mpfr_set_str (x, CHECK_Y2, 10, MPFR_RNDN);
  if (MPFR_IS_NAN (y) || mpfr_cmp (y, x))
    {
//.........这里部分代码省略.........

static void
check_overflow (void)
{
  mpfr_t x;
  char *s;

  mpfr_init (x);

  /* Huge overflow */
  mpfr_strtofr (x, "123456789E2147483646", &s, 0, MPFR_RNDN);
  if (s[0] != 0 || !MPFR_IS_INF (x) || !MPFR_IS_POS (x) )
    {
      printf ("Check overflow failed (1) with:\n s=%s\n x=", s);
      mpfr_dump (x);
      exit (1);
    }
  mpfr_strtofr (x, "123456789E9223372036854775807", &s, 0, MPFR_RNDN);
  if (s[0] != 0 || !MPFR_IS_INF (x) || !MPFR_IS_POS (x) )
    {
      printf ("Check overflow failed (2) with:\n s='%s'\n x=", s);
      mpfr_dump (x);
      exit (1);
    }
  mpfr_strtofr (x, "123456789E170141183460469231731687303715884105728",
                &s, 0, MPFR_RNDN);
  if (s[0] != 0 || !MPFR_IS_INF (x) || !MPFR_IS_POS (x) )
    {
      printf ("Check overflow failed (3) with:\n s=%s\n x=", s);
      mpfr_dump (x);
      exit (1);
    }

  /* Limit overflow */
  mpfr_strtofr (x, "12E2147483646", &s, 0, MPFR_RNDN);
  if (s[0] != 0 || !MPFR_IS_INF (x) || !MPFR_IS_POS (x) )
    {
      printf ("Check overflow failed (4) with:\n s=%s\n x=", s);
      mpfr_dump (x);
      exit (1);
    }
  mpfr_strtofr (x, "12E2147483645", &s, 0, MPFR_RNDN);
  if (s[0] != 0 || !MPFR_IS_INF (x) || !MPFR_IS_POS (x))
    {
      printf ("Check overflow failed (5) with:\n s=%s\n x=", s);
      mpfr_dump (x);
      exit (1);
    }
  mpfr_strtofr (x, "[email protected]", &s, 16, MPFR_RNDN);
  if (s[0] != 0 || !MPFR_IS_INF (x) || !MPFR_IS_POS (x))
    {
      printf ("Check overflow failed (6) with:\n s=%s\n x=", s);
      mpfr_dump (x);
      exit (1);
    }
  mpfr_strtofr (x, "[email protected]", &s, 16, MPFR_RNDN);
  if (s[0] != 0 || !MPFR_IS_INF (x) || !MPFR_IS_POS (x))
    {
      printf ("Check overflow failed (7) with:\n s=%s\n x=", s);
      mpfr_dump (x);
      exit (1);
    }

  /* Check underflow */
  mpfr_strtofr (x, "123456789E-2147483646", &s, 0, MPFR_RNDN);
  if (s[0] != 0 || !MPFR_IS_ZERO (x) || !MPFR_IS_POS (x) )
    {
      printf ("Check underflow failed (1) with:\n s=%s\n x=", s);
      mpfr_dump (x);
      exit (1);
    }
  mpfr_strtofr (x, "123456789E-9223372036854775807", &s, 0, MPFR_RNDN);
  if (s[0] != 0 || !MPFR_IS_ZERO (x) || !MPFR_IS_POS (x) )
    {
      printf ("Check underflow failed (2) with:\n s='%s'\n x=", s);
      mpfr_dump (x);
      exit (1);
    }
  mpfr_strtofr (x, "-123456789E-170141183460469231731687303715884105728",
                &s, 0, MPFR_RNDN);
  if (s[0] != 0 || !MPFR_IS_ZERO (x) || !MPFR_IS_NEG (x) )
    {
      printf ("Check underflow failed (3) with:\n s=%s\n x=", s);
      mpfr_dump (x);
      exit (1);
    }
  mpfr_strtofr (x, "[email protected]", &s, 16, MPFR_RNDN);
  if (s[0] != 0 || !MPFR_IS_ZERO (x) || !MPFR_IS_POS (x))
    {
      printf ("Check overflow failed (7) with:\n s=%s\n x=", s);
      mpfr_dump (x);
      exit (1);
    }

  mpfr_clear (x);
}

/* Assumes that the exponent range has already been extended and if y is
   an integer, then the result is not exact in unbounded exponent range. */
int
mpfr_pow_general (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
                  mpfr_rnd_t rnd_mode, int y_is_integer, mpfr_save_expo_t *expo)
{
  mpfr_t t, u, k, absx;
  int neg_result = 0;
  int k_non_zero = 0;
  int check_exact_case = 0;
  int inexact;
  /* Declaration of the size variable */
  mpfr_prec_t Nz = MPFR_PREC(z);               /* target precision */
  mpfr_prec_t Nt;                              /* working precision */
  mpfr_exp_t err;                              /* error */
  MPFR_ZIV_DECL (ziv_loop);


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

  /* We put the absolute value of x in absx, pointing to the significand
     of x to avoid allocating memory for the significand of absx. */
  MPFR_ALIAS(absx, x, /*sign=*/ 1, /*EXP=*/ MPFR_EXP(x));

  /* We will compute the absolute value of the result. So, let's
     invert the rounding mode if the result is negative. */
  if (MPFR_IS_NEG (x) && is_odd (y))
    {
      neg_result = 1;
      rnd_mode = MPFR_INVERT_RND (rnd_mode);
    }

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

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

  MPFR_ZIV_INIT (ziv_loop, Nt);
  for (;;)
    {
      MPFR_BLOCK_DECL (flags1);

      /* compute exp(y*ln|x|), using MPFR_RNDU to get an upper bound, so
         that we can detect underflows. */
      mpfr_log (t, absx, MPFR_IS_NEG (y) ? MPFR_RNDD : MPFR_RNDU); /* ln|x| */
      mpfr_mul (t, y, t, MPFR_RNDU);                              /* y*ln|x| */
      if (k_non_zero)
        {
          MPFR_LOG_MSG (("subtract k * ln(2)\n", 0));
          mpfr_const_log2 (u, MPFR_RNDD);
          mpfr_mul (u, u, k, MPFR_RNDD);
          /* Error on u = k * log(2): < k * 2^(-Nt) < 1. */
          mpfr_sub (t, t, u, MPFR_RNDU);
          MPFR_LOG_MSG (("t = y * ln|x| - k * ln(2)\n", 0));
          MPFR_LOG_VAR (t);
        }
      /* estimate of the error -- see pow function in algorithms.tex.
         The error on t is at most 1/2 + 3*2^(EXP(t)+1) ulps, which is
         <= 2^(EXP(t)+3) for EXP(t) >= -1, and <= 2 ulps for EXP(t) <= -2.
         Additional error if k_no_zero: treal = t * errk, with
         1 - |k| * 2^(-Nt) <= exp(-|k| * 2^(-Nt)) <= errk <= 1,
         i.e., additional absolute error <= 2^(EXP(k)+EXP(t)-Nt).
         Total error <= 2^err1 + 2^err2 <= 2^(max(err1,err2)+1). */
      err = MPFR_NOTZERO (t) && MPFR_GET_EXP (t) >= -1 ?
        MPFR_GET_EXP (t) + 3 : 1;
      if (k_non_zero)
        {
          if (MPFR_GET_EXP (k) > err)
            err = MPFR_GET_EXP (k);
          err++;
        }
      MPFR_BLOCK (flags1, mpfr_exp (t, t, MPFR_RNDN));  /* exp(y*ln|x|)*/
      /* We need to test */
      if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (t) || MPFR_UNDERFLOW (flags1)))
        {
          mpfr_prec_t Ntmin;
          MPFR_BLOCK_DECL (flags2);

          MPFR_ASSERTN (!k_non_zero);
          MPFR_ASSERTN (!MPFR_IS_NAN (t));

          /* Real underflow? */
          if (MPFR_IS_ZERO (t))
            {
              /* Underflow. We computed rndn(exp(t)), where t >= y*ln|x|.
                 Therefore rndn(|x|^y) = 0, and we have a real underflow on
                 |x|^y. */
              inexact = mpfr_underflow (z, rnd_mode == MPFR_RNDN ? MPFR_RNDZ
                                        : rnd_mode, MPFR_SIGN_POS);
              if (expo != NULL)
                MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, MPFR_FLAGS_INEXACT
                                             | MPFR_FLAGS_UNDERFLOW);
              break;
//.........这里部分代码省略.........

/* The computation of z = pow(x,y) is done by
   z = exp(y * log(x)) = x^y
   For the special cases, see Section F.9.4.4 of the C standard:
     _ pow(±0, y) = ±inf for y an odd integer < 0.
     _ pow(±0, y) = +inf for y < 0 and not an odd integer.
     _ pow(±0, y) = ±0 for y an odd integer > 0.
     _ pow(±0, y) = +0 for y > 0 and not an odd integer.
     _ pow(-1, ±inf) = 1.
     _ pow(+1, y) = 1 for any y, even a NaN.
     _ pow(x, ±0) = 1 for any x, even a NaN.
     _ pow(x, y) = NaN for finite x < 0 and finite non-integer y.
     _ pow(x, -inf) = +inf for |x| < 1.
     _ pow(x, -inf) = +0 for |x| > 1.
     _ pow(x, +inf) = +0 for |x| < 1.
     _ pow(x, +inf) = +inf for |x| > 1.
     _ pow(-inf, y) = -0 for y an odd integer < 0.
     _ pow(-inf, y) = +0 for y < 0 and not an odd integer.
     _ pow(-inf, y) = -inf for y an odd integer > 0.
     _ pow(-inf, y) = +inf for y > 0 and not an odd integer.
     _ pow(+inf, y) = +0 for y < 0.
     _ pow(+inf, y) = +inf for y > 0. */
int
mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
{
  int inexact;
  int cmp_x_1;
  int y_is_integer;
  MPFR_SAVE_EXPO_DECL (expo);

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

  if (MPFR_ARE_SINGULAR (x, y))
    {
      /* pow(x, 0) returns 1 for any x, even a NaN. */
      if (MPFR_UNLIKELY (MPFR_IS_ZERO (y)))
        return mpfr_set_ui (z, 1, rnd_mode);
      else if (MPFR_IS_NAN (x))
        {
          MPFR_SET_NAN (z);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_NAN (y))
        {
          /* pow(+1, NaN) returns 1. */
          if (mpfr_cmp_ui (x, 1) == 0)
            return mpfr_set_ui (z, 1, rnd_mode);
          MPFR_SET_NAN (z);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF (y))
        {
          if (MPFR_IS_INF (x))
            {
              if (MPFR_IS_POS (y))
                MPFR_SET_INF (z);
              else
                MPFR_SET_ZERO (z);
              MPFR_SET_POS (z);
              MPFR_RET (0);
            }
          else
            {
              int cmp;
              cmp = mpfr_cmpabs (x, __gmpfr_one) * MPFR_INT_SIGN (y);
              MPFR_SET_POS (z);
              if (cmp > 0)
                {
                  /* Return +inf. */
                  MPFR_SET_INF (z);
                  MPFR_RET (0);
                }
              else if (cmp < 0)
                {
                  /* Return +0. */
                  MPFR_SET_ZERO (z);
                  MPFR_RET (0);
                }
              else
                {
                  /* Return 1. */
                  return mpfr_set_ui (z, 1, rnd_mode);
                }
            }
        }
      else if (MPFR_IS_INF (x))
        {
          int negative;
          /* Determine the sign now, in case y and z are the same object */
          negative = MPFR_IS_NEG (x) && is_odd (y);
          if (MPFR_IS_POS (y))
            MPFR_SET_INF (z);
          else
            MPFR_SET_ZERO (z);
          if (negative)
            MPFR_SET_NEG (z);
//.........这里部分代码省略.........

static void
test_2exp (void)
{
  mpfr_t x;
  int res;

  mpfr_init2 (x, 32);

  mpfr_set_ui_2exp (x, 1, 0, GMP_RNDN);
  if (mpfr_cmp_ui(x, 1))
    ERROR("(1U,0)");

  mpfr_set_ui_2exp (x, 1024, -10, GMP_RNDN);
  if (mpfr_cmp_ui(x, 1))
    ERROR("(1024U,-10)");

  mpfr_set_ui_2exp (x, 1024, 10, GMP_RNDN);
  if (mpfr_cmp_ui(x, 1024*1024))
    ERROR("(1024U,+10)");

  mpfr_set_si_2exp (x, -1024L * 1024L, -10, GMP_RNDN);
  if (mpfr_cmp_si(x, -1024))
    ERROR("(1M,-10)");

  mpfr_set_ui_2exp (x, 0x92345678, 16, GMP_RNDN);
  if (mpfr_cmp_str (x, "[email protected]", 16, GMP_RNDN))
    ERROR("(x92345678U,+16)");

  mpfr_set_si_2exp (x, -0x1ABCDEF0, -256, GMP_RNDN);
  if (mpfr_cmp_str (x, "[email protected]", 16, GMP_RNDN))
    ERROR("(-x1ABCDEF0,-256)");

  mpfr_set_prec (x, 2);
  res = mpfr_set_si_2exp (x, 7, 10, GMP_RNDU);
  if (mpfr_cmp_ui (x, 1<<13) || res <= 0)
    ERROR ("Prec 2 + si_2exp");

  res = mpfr_set_ui_2exp (x, 7, 10, GMP_RNDU);
  if (mpfr_cmp_ui (x, 1<<13) || res <= 0)
    ERROR ("Prec 2 + ui_2exp");

  mpfr_clear_flags ();
  mpfr_set_ui_2exp (x, 17, MPFR_EMAX_MAX, GMP_RNDN);
  if (!mpfr_inf_p (x) || MPFR_IS_NEG (x))
    ERROR ("mpfr_set_ui_2exp and overflow (bad result)");
  if (!mpfr_overflow_p ())
    ERROR ("mpfr_set_ui_2exp and overflow (overflow flag not set)");

  mpfr_clear_flags ();
  mpfr_set_si_2exp (x, 17, MPFR_EMAX_MAX, GMP_RNDN);
  if (!mpfr_inf_p (x) || MPFR_IS_NEG (x))
    ERROR ("mpfr_set_si_2exp (pos) and overflow (bad result)");
  if (!mpfr_overflow_p ())
    ERROR ("mpfr_set_si_2exp (pos) and overflow (overflow flag not set)");

  mpfr_clear_flags ();
  mpfr_set_si_2exp (x, -17, MPFR_EMAX_MAX, GMP_RNDN);
  if (!mpfr_inf_p (x) || MPFR_IS_POS (x))
    ERROR ("mpfr_set_si_2exp (neg) and overflow (bad result)");
  if (!mpfr_overflow_p ())
    ERROR ("mpfr_set_si_2exp (neg) and overflow (overflow flag not set)");

  mpfr_clear (x);
}

int
mpfr_digamma (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  int inex;
  MPFR_SAVE_EXPO_DECL (expo);

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


  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)) /* Digamma(+Inf) = +Inf */
            {
              MPFR_SET_SAME_SIGN(y, x);
              MPFR_SET_INF(y);
              MPFR_RET(0);
            }
          else                /* Digamma(-Inf) = NaN */
            {
              MPFR_SET_NAN(y);
              MPFR_RET_NAN;
            }
        }
      else /* Zero case */
        {
          /* the following works also in case of overlap */
          MPFR_SET_INF(y);
          MPFR_SET_OPPOSITE_SIGN(y, x);
          mpfr_set_divby0 ();
          MPFR_RET(0);
        }
    }

  /* Digamma is undefined for negative integers */
  if (MPFR_IS_NEG(x) && mpfr_integer_p (x))
    {
      MPFR_SET_NAN(y);
      MPFR_RET_NAN;
    }

  /* now x is a normal number */

  MPFR_SAVE_EXPO_MARK (expo);
  /* for x very small, we have Digamma(x) = -1/x - gamma + O(x), more precisely
     -1 < Digamma(x) + 1/x < 0 for -0.2 < x < 0.2, thus:
     (i) either x is a power of two, then 1/x is exactly representable, and
         as long as 1/2*ulp(1/x) > 1, we can conclude;
     (ii) otherwise assume x has <= n bits, and y has <= n+1 bits, then
   |y + 1/x| >= 2^(-2n) ufp(y), where ufp means unit in first place.
   Since |Digamma(x) + 1/x| <= 1, if 2^(-2n) ufp(y) >= 2, then
   |y - Digamma(x)| >= 2^(-2n-1)ufp(y), and rounding -1/x gives the correct result.
   If x < 2^E, then y > 2^(-E), thus ufp(y) > 2^(-E-1).
   A sufficient condition is thus EXP(x) <= -2 MAX(PREC(x),PREC(Y)). */
  if (MPFR_EXP(x) < -2)
    {
      if (MPFR_EXP(x) <= -2 * (mpfr_exp_t) MAX(MPFR_PREC(x), MPFR_PREC(y)))
        {
          int signx = MPFR_SIGN(x);
          inex = mpfr_si_div (y, -1, x, rnd_mode);
          if (inex == 0) /* x is a power of two */
            { /* result always -1/x, except when rounding down */
              if (rnd_mode == MPFR_RNDA)
                rnd_mode = (signx > 0) ? MPFR_RNDD : MPFR_RNDU;
              if (rnd_mode == MPFR_RNDZ)
                rnd_mode = (signx > 0) ? MPFR_RNDU : MPFR_RNDD;
              if (rnd_mode == MPFR_RNDU)
                inex = 1;
              else if (rnd_mode == MPFR_RNDD)
                {
                  mpfr_nextbelow (y);
                  inex = -1;
                }
              else /* nearest */
                inex = 1;
            }
          MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
          goto end;
        }
    }

  if (MPFR_IS_NEG(x))
    inex = mpfr_digamma_reflection (y, x, rnd_mode);
  /* if x < 1/2 we use the reflection formula */
  else if (MPFR_EXP(x) < 0)
    inex = mpfr_digamma_reflection (y, x, rnd_mode);
  else
    inex = mpfr_digamma_positive (y, x, rnd_mode);

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

static void
special (void)
{
  mpfr_t x, y;
  mpq_t q;
  mpz_t z;
  int res = 0;

  mpfr_init (x);
  mpfr_init (y);
  mpq_init (q);
  mpz_init (z);

  /* cancellation in mpfr_add_q */
  mpfr_set_prec (x, 60);
  mpfr_set_prec (y, 20);
  mpz_set_str (mpq_numref (q), "-187207494", 10);
  mpz_set_str (mpq_denref (q), "5721", 10);
  mpfr_set_str_binary (x, "11111111101001011011100101100011011110010011100010000100001E-44");
  mpfr_add_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("cancelation in add_q", mpfr_cmp_ui_2exp (y, 256783, -64) == 0);

  mpfr_set_prec (x, 19);
  mpfr_set_str_binary (x, "0.1011110101110011100E0");
  mpz_set_str (mpq_numref (q), "187207494", 10);
  mpz_set_str (mpq_denref (q), "5721", 10);
  mpfr_set_prec (y, 29);
  mpfr_add_q (y, x, q, MPFR_RNDD);
  mpfr_set_prec (x, 29);
  mpfr_set_str_binary (x, "11111111101001110011010001001E-14");
  CHECK_FOR ("cancelation in add_q", mpfr_cmp (x,y) == 0);

  /* Inf */
  mpfr_set_inf (x, 1);
  mpz_set_str (mpq_numref (q), "395877315", 10);
  mpz_set_str (mpq_denref (q), "3508975966", 10);
  mpfr_set_prec (y, 118);
  mpfr_add_q (y, x, q, MPFR_RNDU);
  CHECK_FOR ("inf", mpfr_inf_p (y) && mpfr_sgn (y) > 0);
  mpfr_sub_q (y, x, q, MPFR_RNDU);
  CHECK_FOR ("inf", mpfr_inf_p (y) && mpfr_sgn (y) > 0);

  /* Nan */
  MPFR_SET_NAN (x);
  mpfr_add_q (y, x, q, MPFR_RNDU);
  CHECK_FOR ("nan", mpfr_nan_p (y));
  mpfr_sub_q (y, x, q, MPFR_RNDU);
  CHECK_FOR ("nan", mpfr_nan_p (y));

  /* Exact value */
  mpfr_set_prec (x, 60);
  mpfr_set_prec (y, 60);
  mpfr_set_str1 (x, "0.5");
  mpz_set_str (mpq_numref (q), "3", 10);
  mpz_set_str (mpq_denref (q), "2", 10);
  res = mpfr_add_q (y, x, q, MPFR_RNDU);
  CHECK_FOR ("0.5+3/2", mpfr_cmp_ui(y, 2)==0 && res==0);
  res = mpfr_sub_q (y, x, q, MPFR_RNDU);
  CHECK_FOR ("0.5-3/2", mpfr_cmp_si(y, -1)==0 && res==0);

  /* Inf Rationnal */
  mpq_set_ui (q, 1, 0);
  mpfr_set_str1 (x, "0.5");
  res = mpfr_add_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("0.5+1/0", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0);
  res = mpfr_sub_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("0.5-1/0", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0);
  mpq_set_si (q, -1, 0);
  res = mpfr_add_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("0.5+ -1/0", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0);
  res = mpfr_sub_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("0.5- -1/0", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0);
  res = mpfr_div_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("0.5 / (-1/0)", mpfr_zero_p (y) && MPFR_IS_NEG (y) && res == 0);
  mpq_set_ui (q, 1, 0);
  mpfr_set_inf (x, 1);
  res = mpfr_add_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("+Inf + +Inf", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0);
  res = mpfr_sub_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("+Inf - +Inf", MPFR_IS_NAN (y) && res == 0);
  mpfr_set_inf (x, -1);
  res = mpfr_add_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("-Inf + +Inf", MPFR_IS_NAN (y) && res == 0);
  res = mpfr_sub_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("-Inf - +Inf", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0);
  mpq_set_si (q, -1, 0);
  mpfr_set_inf (x, 1);
  res = mpfr_add_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("+Inf + -Inf", MPFR_IS_NAN (y) && res == 0);
  res = mpfr_sub_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("+Inf - -Inf", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0);
  mpfr_set_inf (x, -1);
  res = mpfr_add_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("-Inf + -Inf", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0);
  res = mpfr_sub_q (y, x, q, MPFR_RNDN);
  CHECK_FOR ("-Inf - -Inf", MPFR_IS_NAN (y) && res == 0);

  /* 0 */
  mpq_set_ui (q, 0, 1);
  mpfr_set_ui (x, 42, MPFR_RNDN);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      puts ("");
      exit (1);
    }
  inex = mpfr_set_ui (x, 77617, GMP_RNDD); /* should be 65536 */
  if (MPFR_MANT(x)[0] != ((mp_limb_t)1 << (mp_bits_per_limb-1))
      || inex >= 0)
    {
      printf ("Error in mpfr_set_ui(x:3, 77617, GMP_RNDD)\n");
      mpfr_print_binary (x);
      puts ("");
      exit (1);
    }

  mpfr_set_prec (x, 2);
  inex = mpfr_set_si (x, 33096, GMP_RNDU);
  if (mpfr_get_si (x, GMP_RNDZ) != 49152 || inex <= 0)
    {
      printf ("Error in mpfr_set_si, exp. 49152, got %ld, inex %d\n",
              mpfr_get_si (x, GMP_RNDZ), inex);
      exit (1);
    }
  inex = mpfr_set_ui (x, 33096, GMP_RNDU);
  if (mpfr_get_si (x, GMP_RNDZ) != 49152)
    {
      printf ("Error in mpfr_set_ui, exp. 49152, got %ld, inex %d\n",
              mpfr_get_si (x, GMP_RNDZ), inex);
      exit (1);
    }

  for (r = 0 ; r < GMP_RND_MAX ; r++)
    {
      mpfr_set_si (x, -1, (mp_rnd_t) r);
      mpfr_set_ui (x, 0, (mp_rnd_t) r);
      if (MPFR_IS_NEG (x) || mpfr_get_ui (x, (mp_rnd_t) r) != 0)
        {
          printf ("mpfr_set_ui (x, 0) gives -0 for %s\n",
                  mpfr_print_rnd_mode ((mp_rnd_t) r));
          exit (1);
        }

      mpfr_set_si (x, -1, (mp_rnd_t) r);
      mpfr_set_si (x, 0, (mp_rnd_t) r);
      if (MPFR_IS_NEG (x) || mpfr_get_si (x, (mp_rnd_t) r) != 0)
        {
          printf ("mpfr_set_si (x, 0) gives -0 for %s\n",
                  mpfr_print_rnd_mode ((mp_rnd_t) r));
          exit (1);
        }
    }

  /* check potential bug in case mp_limb_t is unsigned */
  emax = mpfr_get_emax ();
  set_emax (0);
  mpfr_set_si (x, -1, GMP_RNDN);
  if (mpfr_sgn (x) >= 0)
    {
      printf ("mpfr_set_si (x, -1) fails\n");
      exit (1);
    }
  set_emax (emax);

  emax = mpfr_get_emax ();
  set_emax (5);
  mpfr_set_prec (x, 2);
  mpfr_set_si (x, -31, GMP_RNDN);
  if (mpfr_sgn (x) >= 0)

static void
overflowed_sin_cos0 (void)
{
  mpfr_t x, y, z;
  int emax, inex, rnd, err = 0;
  mpfr_exp_t old_emax;

  old_emax = mpfr_get_emax ();

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

  for (emax = -1; emax <= 0; emax++)
    {
      mpfr_set_ui_2exp (z, 1, emax, MPFR_RNDN);
      mpfr_nextbelow (z);
      set_emax (emax);  /* 1 is not representable. */
      /* and if emax < 0, 1 - eps is not representable either. */
      RND_LOOP (rnd)
        {
          mpfr_set_si (x, 0, MPFR_RNDN);
          mpfr_neg (x, x, MPFR_RNDN);
          mpfr_clear_flags ();
          inex = mpfr_sin_cos (x, y, x, (mpfr_rnd_t) rnd);
          if (! mpfr_overflow_p ())
            {
              printf ("Error in overflowed_sin_cos0 (rnd = %s):\n"
                      "  The overflow flag is not set.\n",
                      mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
              err = 1;
            }
          if (! (mpfr_zero_p (x) && MPFR_IS_NEG (x)))
            {
              printf ("Error in overflowed_sin_cos0 (rnd = %s):\n"
                      "  Got sin = ", mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
              mpfr_print_binary (x);
              printf (" instead of -0.\n");
              err = 1;
            }
          if (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)
            {
              if (inex == 0)
                {
                  printf ("Error in overflowed_sin_cos0 (rnd = %s):\n"
                          "  The inexact value must be non-zero.\n",
                          mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
                  err = 1;
                }
              if (! mpfr_equal_p (y, z))
                {
                  printf ("Error in overflowed_sin_cos0 (rnd = %s):\n"
                          "  Got cos = ",
                          mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
                  mpfr_print_binary (y);
                  printf (" instead of 0.11111111E%d.\n", emax);
                  err = 1;
                }
            }
          else
            {
              if (inex == 0)
                {
                  printf ("Error in overflowed_sin_cos0 (rnd = %s):\n"
                          "  The inexact value must be non-zero.\n",
                          mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
                  err = 1;
                }
              if (! (mpfr_inf_p (y) && MPFR_IS_POS (y)))
                {
                  printf ("Error in overflowed_sin_cos0 (rnd = %s):\n"
                          "  Got cos = ",
                          mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
                  mpfr_print_binary (y);
                  printf (" instead of +Inf.\n");
                  err = 1;
                }
            }
        }
      set_emax (old_emax);
    }

  if (err)
    exit (1);
  mpfr_clear (x);
  mpfr_clear (y);
  mpfr_clear (z);
}

static void
check_pow_ui (void)
{
    mpfr_t a, b;
    int res;

    mpfr_init2 (a, 53);
    mpfr_init2 (b, 53);

    /* check in-place operations */
    mpfr_set_str (b, "0.6926773", 10, GMP_RNDN);
    mpfr_pow_ui (a, b, 10, GMP_RNDN);
    mpfr_pow_ui (b, b, 10, GMP_RNDN);
    if (mpfr_cmp (a, b))
    {
        printf ("Error for mpfr_pow_ui (b, b, ...)\n");
        exit (1);
    }

    /* check large exponents */
    mpfr_set_ui (b, 1, GMP_RNDN);
    mpfr_pow_ui (a, b, 4294967295UL, GMP_RNDN);

    mpfr_set_inf (a, -1);
    mpfr_pow_ui (a, a, 4049053855UL, GMP_RNDN);
    if (!mpfr_inf_p (a) || (mpfr_sgn (a) >= 0))
    {
        printf ("Error for (-Inf)^4049053855\n");
        exit (1);
    }

    mpfr_set_inf (a, -1);
    mpfr_pow_ui (a, a, (unsigned long) 30002752, GMP_RNDN);
    if (!mpfr_inf_p (a) || (mpfr_sgn (a) <= 0))
    {
        printf ("Error for (-Inf)^30002752\n");
        exit (1);
    }

    /* Check underflow */
    mpfr_set_str_binary (a, "1E-1");
    res = mpfr_pow_ui (a, a, -mpfr_get_emin (), GMP_RNDN);
    if (MPFR_GET_EXP (a) != mpfr_get_emin () + 1)
    {
        printf ("Error for (1e-1)^MPFR_EMAX_MAX\n");
        mpfr_dump (a);
        exit (1);
    }

    mpfr_set_str_binary (a, "1E-10");
    res = mpfr_pow_ui (a, a, -mpfr_get_emin (), GMP_RNDZ);
    if (!MPFR_IS_ZERO (a))
    {
        printf ("Error for (1e-10)^MPFR_EMAX_MAX\n");
        mpfr_dump (a);
        exit (1);
    }

    /* Check overflow */
    mpfr_set_str_binary (a, "1E10");
    res = mpfr_pow_ui (a, a, ULONG_MAX, GMP_RNDN);
    if (!MPFR_IS_INF (a) || MPFR_SIGN (a) < 0)
    {
        printf ("Error for (1e10)^ULONG_MAX\n");
        exit (1);
    }

    /* Check 0 */
    MPFR_SET_ZERO (a);
    MPFR_SET_POS (a);
    mpfr_set_si (b, -1, GMP_RNDN);
    res = mpfr_pow_ui (b, a, 1, GMP_RNDN);
    if (res != 0 || MPFR_IS_NEG (b))
    {
        printf ("Error for (0+)^1\n");
        exit (1);
    }
    MPFR_SET_ZERO (a);
    MPFR_SET_NEG (a);
    mpfr_set_ui (b, 1, GMP_RNDN);
    res = mpfr_pow_ui (b, a, 5, GMP_RNDN);
    if (res != 0 || MPFR_IS_POS (b))
    {
        printf ("Error for (0-)^5\n");
        exit (1);
    }
    MPFR_SET_ZERO (a);
    MPFR_SET_NEG (a);
    mpfr_set_si (b, -1, GMP_RNDN);
    res = mpfr_pow_ui (b, a, 6, GMP_RNDN);
    if (res != 0 || MPFR_IS_NEG (b))
    {
        printf ("Error for (0-)^6\n");
        exit (1);
    }

    mpfr_set_prec (a, 122);
    mpfr_set_prec (b, 122);
    mpfr_set_str_binary (a, "0.10000010010000111101001110100101101010011110011100001111000001001101000110011001001001001011001011010110110110101000111011E1");
    mpfr_set_str_binary (b, "0.11111111100101001001000001000001100011100000001110111111100011111000111011100111111111110100011000111011000100100011001011E51290375");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    mpfr_set_prec (y, 21);
    mpfr_set_prec (z, 21);
    mpfr_set_str (x, "1.11111100100001100101", 2, GMP_RNDN);
    test_pow (z, x, t, GMP_RNDZ);
    mpfr_set_str (y, "1.01101011010001100000e-1", 2, GMP_RNDN);
    if (mpfr_cmp (z, y))
    {
        printf ("Error in mpfr_pow for prec=21, rnd=GMP_RNDZ\n");
        printf ("Expected ");
        mpfr_out_str (stdout, 2, 0, y, GMP_RNDN);
        printf ("\nGot      ");
        mpfr_out_str (stdout, 2, 0, z, GMP_RNDN);
        printf ("\n");
        exit (1);
    }

    /* From http://www.terra.es/personal9/ismaeljc/hall.htm */
    mpfr_set_prec (x, 113);
    mpfr_set_prec (y, 2);
    mpfr_set_prec (z, 169);
    mpfr_set_str1 (x, "6078673043126084065007902175846955");
    mpfr_set_ui_2exp (y, 3, -1, GMP_RNDN);
    test_pow (z, x, y, GMP_RNDZ);
    if (mpfr_cmp_str1 (z, "473928882491000966028828671876527456070714790264144"))
    {
        printf ("Error in mpfr_pow for 6078673043126084065007902175846955");
        printf ("^(3/2), GMP_RNDZ\nExpected ");
        printf ("4.73928882491000966028828671876527456070714790264144e50");
        printf ("\nGot      ");
        mpfr_out_str (stdout, 10, 0, z, GMP_RNDN);
        printf ("\n");
        exit (1);
    }
    test_pow (z, x, y, GMP_RNDU);
    if (mpfr_cmp_str1 (z, "473928882491000966028828671876527456070714790264145"))
    {
        printf ("Error in mpfr_pow for 6078673043126084065007902175846955");
        printf ("^(3/2), GMP_RNDU\nExpected ");
        printf ("4.73928882491000966028828671876527456070714790264145e50");
        printf ("\nGot      ");
        mpfr_out_str (stdout, 10, 0, z, GMP_RNDN);
        printf ("\n");
        exit (1);
    }

    mpfr_set_inf (x, 1);
    mpfr_set_prec (y, 2);
    mpfr_set_str_binary (y, "1E10");
    test_pow (z, x, y, GMP_RNDN);
    MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
    mpfr_set_inf (x, -1);
    test_pow (z, x, y, GMP_RNDN);
    MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
    mpfr_set_prec (y, 10);
    mpfr_set_str_binary (y, "1.000000001E9");
    test_pow (z, x, y, GMP_RNDN);
    MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_NEG(z));
    mpfr_set_str_binary (y, "1.000000001E8");
    test_pow (z, x, y, GMP_RNDN);
    MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));

    mpfr_set_inf (x, -1);
    mpfr_set_prec (y, 2 * mp_bits_per_limb);
    mpfr_set_ui (y, 1, GMP_RNDN);
    mpfr_mul_2exp (y, y, mp_bits_per_limb - 1, GMP_RNDN);
    /* y = 2^(mp_bits_per_limb - 1) */
    test_pow (z, x, y, GMP_RNDN);
    MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
    mpfr_nextabove (y);
    test_pow (z, x, y, GMP_RNDN);
    /* y = 2^(mp_bits_per_limb - 1) + epsilon */
    MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
    mpfr_nextbelow (y);
    mpfr_div_2exp (y, y, 1, GMP_RNDN);
    mpfr_nextabove (y);
    test_pow (z, x, y, GMP_RNDN);
    /* y = 2^(mp_bits_per_limb - 2) + epsilon */
    MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));

    mpfr_set_si (x, -1, GMP_RNDN);
    mpfr_set_prec (y, 2);
    mpfr_set_str_binary (y, "1E10");
    test_pow (z, x, y, GMP_RNDN);
    MPFR_ASSERTN(mpfr_cmp_ui (z, 1) == 0);

    /* Check (-0)^(17.0001) */
    mpfr_set_prec (x, 6);
    mpfr_set_prec (y, 640);
    MPFR_SET_ZERO (x);
    MPFR_SET_NEG (x);
    mpfr_set_ui (y, 17, GMP_RNDN);
    mpfr_nextabove (y);
    test_pow (z, x, y, GMP_RNDN);
    MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_POS (z));

    mpfr_clear (x);
    mpfr_clear (y);
    mpfr_clear (z);
    mpfr_clear (t);
}

int
main (int argc, char *argv[])
{
  mpfr_t x, y;

  tests_start_mpfr ();

  mpfr_init (x);
  mpfr_init (y);

  /* special values */
  mpfr_set_nan (x);
  mpfr_y0 (y, x, GMP_RNDN);
  MPFR_ASSERTN(mpfr_nan_p (y));

  mpfr_set_inf (x, 1); /* +Inf */
  mpfr_y0 (y, x, GMP_RNDN);
  MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_POS (y));

  mpfr_set_inf (x, -1); /* -Inf */
  mpfr_y0 (y, x, GMP_RNDN);
  MPFR_ASSERTN(mpfr_nan_p (y));

  mpfr_set_ui (x, 0, GMP_RNDN); /* +0 */
  mpfr_y0 (y, x, GMP_RNDN);
  MPFR_ASSERTN(mpfr_inf_p (y) && MPFR_IS_NEG (y)); /* y0(+0)=-Inf */

  mpfr_set_ui (x, 0, GMP_RNDN);
  mpfr_neg (x, x, GMP_RNDN); /* -0 */
  mpfr_y0 (y, x, GMP_RNDN);
  MPFR_ASSERTN(mpfr_inf_p (y) && MPFR_IS_NEG (y)); /* y0(-0)=-Inf */

  mpfr_set_prec (x, 53);
  mpfr_set_prec (y, 53);

  mpfr_set_ui (x, 1, GMP_RNDN);
  mpfr_y0 (y, x, GMP_RNDN);
  mpfr_set_str_binary (x, "0.00010110100110000000001000100110111100110101100011011111");
  if (mpfr_cmp (x, y))
    {
      printf ("Error in mpfr_y0 for x=1, rnd=GMP_RNDN\n");
      printf ("Expected "); mpfr_dump (x);
      printf ("Got      "); mpfr_dump (y);
      exit (1);
    }

  mpfr_set_si (x, -1, GMP_RNDN);
  mpfr_y0 (y, x, GMP_RNDN);
  if (!mpfr_nan_p (y))
    {
      printf ("Error in mpfr_y0 for x=-1, rnd=GMP_RNDN\n");
      printf ("Expected NaN\n");
      printf ("Got      "); mpfr_dump (y);
      exit (1);
    }

  mpfr_clear (x);
  mpfr_clear (y);

  test_generic (2, 100, 1);

  data_check ("data/y0", mpfr_y0, "mpfr_y0");

  tests_end_mpfr ();

  return 0;
}