C++ MPFR_IS_INF函数代码示例

int
main (int argc, char *argv[])
{
#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
  unsigned int prec, yprec;
  int rnd;
  mpfr_t x, y, z, t;
  unsigned long n;
  int inex;

  tests_start_mpfr ();

  mpfr_init (x);
  mpfr_init (y);
  mpfr_init (z);
  mpfr_init (t);

  if (argc >= 3) /* tzeta_ui n prec [rnd] */
    {
      mpfr_set_prec (x, atoi (argv[2]));
      mpfr_zeta_ui (x, atoi (argv[1]),
                    argc > 3 ? (mpfr_rnd_t) atoi (argv[3]) : MPFR_RNDN);
      mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
      printf ("\n");
      goto clear_and_exit;
    }

  mpfr_set_prec (x, 33);
  mpfr_set_prec (y, 33);
  mpfr_zeta_ui (x, 3, MPFR_RNDZ);
  mpfr_set_str_binary (y, "0.100110011101110100000000001001111E1");
  if (mpfr_cmp (x, y))
    {
      printf ("Error for zeta(3), prec=33, MPFR_RNDZ\n");
      printf ("expected "); mpfr_dump (y);
      printf ("got      "); mpfr_dump (x);
      exit (1);
    }

  mpfr_clear_divby0 ();
  inex = mpfr_zeta_ui (x, 0, MPFR_RNDN);
  MPFR_ASSERTN (inex == 0 && mpfr_cmp_si_2exp (x, -1, -1) == 0
                && !mpfr_divby0_p ());

  mpfr_clear_divby0 ();
  inex = mpfr_zeta_ui (x, 1, MPFR_RNDN);
  MPFR_ASSERTN (inex == 0 && MPFR_IS_INF (x) && MPFR_IS_POS (x)
                && mpfr_divby0_p ());

  for (prec = MPFR_PREC_MIN; prec <= 100; prec++)
    {
      mpfr_set_prec (x, prec);
      mpfr_set_prec (z, prec);
      mpfr_set_prec (t, prec);
      yprec = prec + 10;
      mpfr_set_prec (y, yprec);

      for (n = 0; n < 50; n++)
        for (rnd = 0; rnd < MPFR_RND_MAX; rnd++)
          {
            mpfr_zeta_ui (y, n, MPFR_RNDN);
            if (mpfr_can_round (y, yprec, MPFR_RNDN, MPFR_RNDZ, prec
                                + (rnd == MPFR_RNDN)))
              {
                mpfr_set (t, y, (mpfr_rnd_t) rnd);
                mpfr_zeta_ui (z, n, (mpfr_rnd_t) rnd);
                if (mpfr_cmp (t, z))
                  {
                    printf ("results differ for n=%lu", n);
                    printf (" prec=%u rnd_mode=%s\n", prec,
                            mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
                    printf ("   got      ");
                    mpfr_dump (z);
                    printf ("   expected ");
                    mpfr_dump (t);
                    printf ("   approx   ");
                    mpfr_dump (y);
                    exit (1);
                  }
              }
          }
    }

 clear_and_exit:
  mpfr_clear (x);
  mpfr_clear (y);
  mpfr_clear (z);
  mpfr_clear (t);

  tests_end_mpfr ();
#endif
  return 0;
}

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

  tests_start_mpfr ();

  special_overflow ();
  emax_m_eps ();
  exp_range ();

  mpfr_init (x);
  mpfr_init (y);

  mpfr_set_ui (x, 4, MPFR_RNDN);
  mpfr_exp2 (y, x, MPFR_RNDN);
  if (mpfr_cmp_ui (y, 16) != 0)
    {
      printf ("Error for 2^4, MPFR_RNDN\n");
      exit (1);
    }
  mpfr_exp2 (y, x, MPFR_RNDD);
  if (mpfr_cmp_ui (y, 16) != 0)
    {
      printf ("Error for 2^4, MPFR_RNDD\n");
      exit (1);
    }
  mpfr_exp2 (y, x, MPFR_RNDU);
  if (mpfr_cmp_ui (y, 16) != 0)
    {
      printf ("Error for 2^4, MPFR_RNDU\n");
      exit (1);
    }

  mpfr_set_si (x, -4, MPFR_RNDN);
  mpfr_exp2 (y, x, MPFR_RNDN);
  if (mpfr_cmp_ui_2exp (y, 1, -4) != 0)
    {
      printf ("Error for 2^(-4), MPFR_RNDN\n");
      exit (1);
    }
  mpfr_exp2 (y, x, MPFR_RNDD);
  if (mpfr_cmp_ui_2exp (y, 1, -4) != 0)
    {
      printf ("Error for 2^(-4), MPFR_RNDD\n");
      exit (1);
    }
  mpfr_exp2 (y, x, MPFR_RNDU);
  if (mpfr_cmp_ui_2exp (y, 1, -4) != 0)
    {
      printf ("Error for 2^(-4), MPFR_RNDU\n");
      exit (1);
    }

  mpfr_set_prec (x, 53);
  mpfr_set_prec (y, 53);
  mpfr_set_str (x, /*-1683977482443233.0 / 2199023255552.0*/
                "-7.6578429909351734750089235603809357e2", 10, MPFR_RNDN);
  mpfr_exp2 (y, x, MPFR_RNDN);
  if (mpfr_cmp_str1 (y, "2.991959870867646566478e-231"))
    {
      printf ("Error for x=-1683977482443233/2^41\n");
      exit (1);
    }

  mpfr_set_prec (x, 10);
  mpfr_set_prec (y, 10);
  /* save emin */
  emin = mpfr_get_emin ();
  set_emin (-10);
  mpfr_set_si (x, -12, MPFR_RNDN);
  mpfr_exp2 (y, x, MPFR_RNDN);
  if (mpfr_cmp_ui (y, 0) || mpfr_sgn (y) < 0)
    {
      printf ("Error for x=emin-2, RNDN\n");
      printf ("Expected +0\n");
      printf ("Got      "); mpfr_print_binary (y); puts ("");
      exit (1);
    }
  /* restore emin */
  set_emin (emin);

  /* save emax */
  emax = mpfr_get_emax ();
  set_emax (10);
  mpfr_set_ui (x, 11, MPFR_RNDN);
  mpfr_exp2 (y, x, MPFR_RNDN);
  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
    {
      printf ("Error for x=emax+1, RNDN\n");
      exit (1);
    }
  /* restore emax */
  set_emax (emax);

  MPFR_SET_INF(x);
  MPFR_SET_POS(x);
  mpfr_exp2 (y, x, MPFR_RNDN);
  if(!MPFR_IS_INF(y))
//.........这里部分代码省略.........

static void
check_max(void)
{
  mpfr_t xx, yy, zz;
  mpfr_exp_t emin;

  mpfr_init2(xx, 4);
  mpfr_init2(yy, 4);
  mpfr_init2(zz, 4);
  mpfr_set_str1 (xx, "0.68750");
  mpfr_mul_2si(xx, xx, MPFR_EMAX_DEFAULT/2, MPFR_RNDN);
  mpfr_set_str1 (yy, "0.68750");
  mpfr_mul_2si(yy, yy, MPFR_EMAX_DEFAULT - MPFR_EMAX_DEFAULT/2 + 1, MPFR_RNDN);
  mpfr_clear_flags();
  test_mul(zz, xx, yy, MPFR_RNDU);
  if (!(mpfr_overflow_p() && MPFR_IS_INF(zz)))
    {
      printf("check_max failed (should be an overflow)\n");
      exit(1);
    }

  mpfr_clear_flags();
  test_mul(zz, xx, yy, MPFR_RNDD);
  if (mpfr_overflow_p() || MPFR_IS_INF(zz))
    {
      printf("check_max failed (should NOT be an overflow)\n");
      exit(1);
    }
  mpfr_set_str1 (xx, "0.93750");
  mpfr_mul_2si(xx, xx, MPFR_EMAX_DEFAULT, MPFR_RNDN);
  if (!(MPFR_IS_FP(xx) && MPFR_IS_FP(zz)))
    {
      printf("check_max failed (internal error)\n");
      exit(1);
    }
  if (mpfr_cmp(xx, zz) != 0)
    {
      printf("check_max failed: got ");
      mpfr_out_str(stdout, 2, 0, zz, MPFR_RNDZ);
      printf(" instead of ");
      mpfr_out_str(stdout, 2, 0, xx, MPFR_RNDZ);
      printf("\n");
      exit(1);
    }

  /* check underflow */
  emin = mpfr_get_emin ();
  set_emin (0);
  mpfr_set_str_binary (xx, "0.1E0");
  mpfr_set_str_binary (yy, "0.1E0");
  test_mul (zz, xx, yy, MPFR_RNDN);
  /* exact result is 0.1E-1, which should round to 0 */
  MPFR_ASSERTN(mpfr_cmp_ui (zz, 0) == 0 && MPFR_IS_POS(zz));
  set_emin (emin);

  /* coverage test for mpfr_powerof2_raw */
  emin = mpfr_get_emin ();
  set_emin (0);
  mpfr_set_prec (xx, mp_bits_per_limb + 1);
  mpfr_set_str_binary (xx, "0.1E0");
  mpfr_nextabove (xx);
  mpfr_set_str_binary (yy, "0.1E0");
  test_mul (zz, xx, yy, MPFR_RNDN);
  /* exact result is just above 0.1E-1, which should round to minfloat */
  MPFR_ASSERTN(mpfr_cmp (zz, yy) == 0);
  set_emin (emin);

  mpfr_clear(xx);
  mpfr_clear(yy);
  mpfr_clear(zz);
}

int
mpfr_cmp_si_2exp (mpfr_srcptr b, long int i, mpfr_exp_t f)
{
  int si;

  si = i < 0 ? -1 : 1; /* sign of i */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (b)))
    {
      if (MPFR_IS_INF(b))
        return MPFR_INT_SIGN(b);
      else if (MPFR_IS_ZERO(b))
        return i != 0 ? -si : 0;
      /* NAN */
      MPFR_SET_ERANGEFLAG ();
      return 0;
    }
  else if (MPFR_SIGN(b) != si || i == 0)
    return MPFR_INT_SIGN (b);
  else /* b and i are of same sign si */
    {
      mpfr_exp_t e;
      unsigned long ai;
      int k;
      mp_size_t bn;
      mp_limb_t c, *bp;

      ai = SAFE_ABS(unsigned long, i);

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

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

      /* now f < e <= f + GMP_NUMB_BITS */
      c = (mp_limb_t) ai;
      count_leading_zeros(k, c);
      if ((int) (e - f) > GMP_NUMB_BITS - k)
        return si;
      if ((int) (e - f) < GMP_NUMB_BITS - k)
        return -si;

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

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

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

/* computes tan(x) = sign(x)*sqrt(1/cos(x)^2-1) */
int
mpfr_tan (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
{
  mp_prec_t precy, m;
  int inexact;
  mpfr_t s, c;
  MPFR_ZIV_DECL (loop);
  MPFR_SAVE_EXPO_DECL (expo);
  MPFR_GROUP_DECL (group);

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

  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)))
    {
      if (MPFR_IS_NAN(x) || MPFR_IS_INF(x))
        {
          MPFR_SET_NAN(y);
          MPFR_RET_NAN;
        }
      else /* x is zero */
        {
          MPFR_ASSERTD(MPFR_IS_ZERO(x));
          MPFR_SET_ZERO(y);
          MPFR_SET_SAME_SIGN(y, x);
          MPFR_RET(0);
        }
    }

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

  MPFR_SAVE_EXPO_MARK (expo);

  /* Compute initial precision */
  precy = MPFR_PREC (y);
  m = precy + MPFR_INT_CEIL_LOG2 (precy) + 13;
  MPFR_ASSERTD (m >= 2); /* needed for the error analysis in algorithms.tex */

  MPFR_GROUP_INIT_2 (group, m, s, c);
  MPFR_ZIV_INIT (loop, m);
  for (;;)
    {
      /* The only way to get an overflow is to get ~ Pi/2
         But the result will be ~ 2^Prec(y). */
      mpfr_sin_cos (s, c, x, GMP_RNDN); /* err <= 1/2 ulp on s and c */
      mpfr_div (c, s, c, GMP_RNDN);     /* err <= 4 ulps */
      MPFR_ASSERTD (!MPFR_IS_SINGULAR (c));
      if (MPFR_LIKELY (MPFR_CAN_ROUND (c, m - 2, precy, rnd_mode)))
        break;
      MPFR_ZIV_NEXT (loop, m);
      MPFR_GROUP_REPREC_2 (group, m, s, c);
    }
  MPFR_ZIV_FREE (loop);
  inexact = mpfr_set (y, c, rnd_mode);
  MPFR_GROUP_CLEAR (group);

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

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

  mpfr_init (x);
  mpfr_init (y);

  MPFR_SET_INF(x);
  mpfr_set_ui (y, 1, GMP_RNDN);
  mpfr_acosh (x, y, GMP_RNDN);
  if (MPFR_IS_INF(x) || MPFR_IS_NAN(x) )
    {
      printf ("Inf flag not clears in acosh!\n");
      exit (1);
    }
  if (mpfr_cmp_ui (x, 0))
    {
      printf ("Error: mpfr_acosh(1) <> 0\n");
      exit (1);
    }

  MPFR_SET_NAN(x);
  mpfr_acosh (x, y, GMP_RNDN);
  if (MPFR_IS_NAN(x) || MPFR_IS_INF(x) )
    {
      printf ("NAN flag not clears in acosh!\n");
      exit (1);
    }

  mpfr_set_ui (x, 0, GMP_RNDN);
  mpfr_acosh (y, x, GMP_RNDN);
  if (!mpfr_nan_p (y))
    {
      printf ("Error: mpfr_acosh(0) <> NaN\n");
      exit (1);
    }

  mpfr_set_si (x, -1, GMP_RNDN);
  mpfr_acosh (y, x, GMP_RNDN);
  if (!mpfr_nan_p (y))
    {
      printf ("Error: mpfr_acosh(-1) <> NaN\n");
      exit (1);
    }

  MPFR_SET_NAN(x);
  mpfr_acosh (y, x, GMP_RNDN);
  if (!mpfr_nan_p (y))
    {
      printf ("Error: mpfr_acosh(NaN) <> NaN\n");
      exit (1);
    }

  mpfr_set_inf (x, 1);
  mpfr_acosh (y, x, GMP_RNDN);
  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
    {
      printf ("Error: mpfr_acosh(+Inf) <> +Inf\n");
      exit (1);
    }

  mpfr_set_inf (x, -1);
  mpfr_acosh (y, x, GMP_RNDN);
  if (!mpfr_nan_p (y))
    {
      printf ("Error: mpfr_acosh(-Inf) <> NaN\n");
      exit (1);
    }

  mpfr_set_ui (x, 1, GMP_RNDN);
  mpfr_div_2exp (x, x, 1, GMP_RNDN);
  mpfr_acosh (y, x, GMP_RNDN);
  if (!mpfr_nan_p (y))
    {
      printf ("Error: mpfr_acosh(1/2) <> NaN\n");
      exit (1);
    }

  mpfr_set_prec (x, 32);
  mpfr_set_prec (y, 32);
  mpfr_set_str_binary (x, "1.000001101011101111001011");
  mpfr_acosh (y, x, GMP_RNDN);
  mpfr_set_str_binary (x, "0.111010100101101001010001101001E-2");
  if (mpfr_cmp (x, y))
    {
      printf ("Error: mpfr_acosh (1)\n");
      exit (1);
    }

  mpfr_clear (x);
  mpfr_clear (y);
}

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
main (int argc, char *argv[])
{
  mpfr_t w,z;
  unsigned long k;
  int i;

  tests_start_mpfr ();

  mpfr_inits2 (53, w, z, (mpfr_ptr) 0);

  for (i = 0; i < 3; i++)
    {
      mpfr_set_inf (w, 1);
      test_mul (i, 0, w, w, 10, GMP_RNDZ);
      if (!MPFR_IS_INF(w))
        {
          printf ("Result is not Inf (i = %d)\n", i);
          exit (1);
        }

      mpfr_set_nan (w);
      test_mul (i, 0, w, w, 10, GMP_RNDZ);
      if (!MPFR_IS_NAN(w))
        {
          printf ("Result is not NaN (i = %d)\n", i);
          exit (1);
        }

      for (k = 0 ; k < numberof(val) ; k+=3)
        {
          mpfr_set_str (w, val[k], 16, GMP_RNDN);
          test_mul (i, 0, z, w, 10, GMP_RNDZ);
          if (mpfr_cmp_str (z, val[k+1], 16, GMP_RNDN))
            {
              printf ("ERROR for x * 2^n (i = %d) for %s\n", i, val[k]);
              printf ("Expected: %s\n"
                      "Got     : ", val[k+1]);
              mpfr_out_str (stdout, 16, 0, z, GMP_RNDN);
              putchar ('\n');
              exit (1);
            }
          test_mul (i, 1, z, w, 10, GMP_RNDZ);
          if (mpfr_cmp_str (z, val[k+2], 16, GMP_RNDN))
            {
              printf ("ERROR for x / 2^n (i = %d) for %s\n", i, val[k]);
              printf ("Expected: %s\n"
                      "Got     : ", val[k+2]);
              mpfr_out_str (stdout, 16, 0, z, GMP_RNDN);
              putchar ('\n');
              exit (1);
            }
        }

      mpfr_set_inf (w, 1);
      mpfr_nextbelow (w);
      test_mul (i, 0, w, w, 1, GMP_RNDN);
      if (!mpfr_inf_p (w))
        {
          printf ("Overflow error (i = %d)!\n", i);
          exit (1);
        }
      mpfr_set_ui (w, 0, GMP_RNDN);
      mpfr_nextabove (w);
      test_mul (i, 1, w, w, 1, GMP_RNDN);
      if (mpfr_cmp_ui (w, 0))
        {
          printf ("Underflow error (i = %d)!\n", i);
          exit (1);
        }
    }

  if (MPFR_EXP_MAX >= LONG_MAX/2 && MPFR_EXP_MIN <= LONG_MAX/2-LONG_MAX-1)
    {
      unsigned long lmp1 = (unsigned long) LONG_MAX + 1;

      mpfr_set_ui (w, 1, GMP_RNDN);
      mpfr_mul_2ui (w, w, LONG_MAX/2, GMP_RNDZ);
      mpfr_div_2ui (w, w, lmp1, GMP_RNDZ);
      mpfr_mul_2ui (w, w, lmp1 - LONG_MAX/2, GMP_RNDZ);
      if (!mpfr_cmp_ui (w, 1))
        {
          printf ("Underflow LONG_MAX error!\n");
          exit (1);
        }
    }

  mpfr_clears (w, z, (mpfr_ptr) 0);

  underflow0 ();

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

//.........这里部分代码省略.........
                       flags == oldflags))
                  {
                    printf ("tgeneric: error for " MAKE_STR(TEST_FUNCTION)
                            ", reduced exponent range [%"
                            MPFR_EXP_FSPEC "d,%" MPFR_EXP_FSPEC "d] on:\n",
                            (mpfr_eexp_t) emin, (mpfr_eexp_t) emax);
                    printf ("x = ");
                    mpfr_dump (x);
#if defined(TWO_ARGS_ALL)
                    printf ("u = ");
                    mpfr_dump (u);
#endif
                    printf ("yprec = %u, rnd_mode = %s\n",
                            (unsigned int) yprec, mpfr_print_rnd_mode (rnd));
                    printf ("Expected:\n  y = ");
                    mpfr_dump (y);
                    printf ("  inex = %d, flags =", compare);
                    flags_out (oldflags);
                    printf ("Got:\n  w = ");
                    mpfr_dump (w);
                    printf ("  inex = %d, flags =", inexact);
                    flags_out (flags);
                    exit (1);
                  }
              }
            __gmpfr_flags = oldflags;  /* restore the flags */
          }

          if (MPFR_IS_SINGULAR (y))
            {
              if (MPFR_IS_NAN (y) || mpfr_nanflag_p ())
                TGENERIC_CHECK ("bad NaN flag",
                                MPFR_IS_NAN (y) && mpfr_nanflag_p ());
              else if (MPFR_IS_INF (y))
                {
                  TGENERIC_CHECK ("bad overflow flag",
                                  (compare != 0) ^ (mpfr_overflow_p () == 0));
                  TGENERIC_CHECK ("bad divide-by-zero flag",
                                  (compare == 0 && !infinite_input) ^
                                  (mpfr_divby0_p () == 0));
                }
              else if (MPFR_IS_ZERO (y))
                TGENERIC_CHECK ("bad underflow flag",
                                (compare != 0) ^ (mpfr_underflow_p () == 0));
            }
          else if (mpfr_divby0_p ())
            {
              TGENERIC_CHECK ("both overflow and divide-by-zero",
                              ! mpfr_overflow_p ());
              TGENERIC_CHECK ("both underflow and divide-by-zero",
                              ! mpfr_underflow_p ());
              TGENERIC_CHECK ("bad compare value (divide-by-zero)",
                              compare == 0);
            }
          else if (mpfr_overflow_p ())
            {
              TGENERIC_CHECK ("both underflow and overflow",
                              ! mpfr_underflow_p ());
              TGENERIC_CHECK ("bad compare value (overflow)", compare != 0);
              mpfr_nexttoinf (y);
              TGENERIC_CHECK ("should have been max MPFR number (overflow)",
                              MPFR_IS_INF (y));
            }
          else if (mpfr_underflow_p ())
            {
              TGENERIC_CHECK ("bad compare value (underflow)", compare != 0);

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

  if (MPFR_ARE_SINGULAR (b, c))
    {
      if (MPFR_IS_NAN (b) || MPFR_IS_NAN (c))
        {
          MPFR_SET_ERANGE ();
          return 0;
        }
      else if (MPFR_IS_INF (b))
        return ! MPFR_IS_INF (c);
      else if (MPFR_IS_INF (c))
        return -1;
      else if (MPFR_IS_ZERO (c))
        return ! MPFR_IS_ZERO (b);
      else /* b == 0 */
        return -1;
    }

  MPFR_ASSERTD (MPFR_IS_PURE_FP (b));
  MPFR_ASSERTD (MPFR_IS_PURE_FP (c));

  /* Now that we know that b and c are pure FP numbers (i.e. they have
     a meaningful exponent), we use MPFR_EXP instead of MPFR_GET_EXP to
     allow exponents outside the current exponent range. For instance,
     this is useful for mpfr_pow, which compares values to __gmpfr_one.
     This is for internal use only! For compatibility with other MPFR
     versions, the user must still provide values that are representable
     in the current exponent range. */
  be = MPFR_EXP (b);
  ce = MPFR_EXP (c);
  if (be > ce)
    return 1;
  if (be < ce)
    return -1;

  /* exponents are equal */

  bn = MPFR_LIMB_SIZE(b)-1;
  cn = MPFR_LIMB_SIZE(c)-1;

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

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

  for ( ; bn >= 0; bn--)
    if (bp[bn])
      return 1;

  for ( ; cn >= 0; cn--)
    if (cp[cn])
      return -1;

   return 0;
}

//.........这里部分代码省略.........
                  mpfr_set_si (x, n == 0 ? 1 : -1, GMP_RNDN);
                  mpfr_setmax (x, REDUCE_EMAX);
#if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
                  mpfr_set_si (u, randlimb () % 2 == 0 ? 1 : -1, GMP_RNDN);
                  mpfr_setmax (u, mpfr_get_emax ());
#endif
                }
            }

          rnd = RND_RAND ();
          mpfr_clear_flags ();
#ifdef DEBUG_TGENERIC
          TGENERIC_INFO (TEST_FUNCTION, MPFR_PREC (y));
#endif
#if defined(TWO_ARGS)
          compare = TEST_FUNCTION (y, x, u, rnd);
#elif defined(DOUBLE_ARG1)
          d = mpfr_get_d (u, rnd);
          compare = TEST_FUNCTION (y, d, x, rnd);
#elif defined(DOUBLE_ARG2)
          d = mpfr_get_d (u, rnd);
          compare = TEST_FUNCTION (y, x, d, rnd);
#else
          compare = TEST_FUNCTION (y, x, rnd);
#endif
          TGENERIC_CHECK ("Bad inexact flag",
                          (compare != 0) ^ (mpfr_inexflag_p () == 0));
          ctrt++;
          if (MPFR_IS_SINGULAR (y))
            {
              if (MPFR_IS_NAN (y) || mpfr_nanflag_p ())
                TGENERIC_CHECK ("Bad NaN flag",
                                MPFR_IS_NAN (y) && mpfr_nanflag_p ());
              else if (MPFR_IS_INF (y))
                TGENERIC_CHECK ("Bad overflow flag",
                                (compare != 0) ^ (mpfr_overflow_p () == 0));
              else if (MPFR_IS_ZERO (y))
                TGENERIC_CHECK ("Bad underflow flag",
                                (compare != 0) ^ (mpfr_underflow_p () == 0));
            }
          else if (mpfr_overflow_p ())
            {
              TGENERIC_CHECK ("Bad compare value (overflow)", compare != 0);
              mpfr_nexttoinf (y);
              TGENERIC_CHECK ("Should have been max MPFR number",
                              MPFR_IS_INF (y));
            }
          else if (mpfr_underflow_p ())
            {
              TGENERIC_CHECK ("Bad compare value (underflow)", compare != 0);
              mpfr_nexttozero (y);
              TGENERIC_CHECK ("Should have been min MPFR number",
                              MPFR_IS_ZERO (y));
            }
          else if (mpfr_can_round (y, yprec, rnd, rnd, prec))
            {
              ctrn++;
              mpfr_set (t, y, rnd);
              /* Risk of failures are known when some flags are already set
                 before the function call. Do not set the erange flag, as
                 it will remain set after the function call and no checks
                 are performed in such a case (see the mpfr_erangeflag_p
                 test below). */
              if (randlimb () & 1)
                __gmpfr_flags = MPFR_FLAGS_ALL ^ MPFR_FLAGS_ERANGE;
#ifdef DEBUG_TGENERIC

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_tanh (mpfr_ptr y, mpfr_srcptr xt , mpfr_rnd_t rnd_mode)
{
  /****** Declaration ******/
  mpfr_t x;
  int inexact;
  MPFR_SAVE_EXPO_DECL (expo);

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

  /* Special value checking */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
    {
      if (MPFR_IS_NAN (xt))
        {
          MPFR_SET_NAN (y);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF (xt))
        {
          /* tanh(inf) = 1 && tanh(-inf) = -1 */
          return mpfr_set_si (y, MPFR_INT_SIGN (xt), rnd_mode);
        }
      else /* tanh (0) = 0 and xt is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO(xt));
          MPFR_SET_ZERO (y);
          MPFR_SET_SAME_SIGN (y, xt);
          MPFR_RET (0);
        }
    }

  /* tanh(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, 0,
                                    rnd_mode, {});

  MPFR_TMP_INIT_ABS (x, xt);

  MPFR_SAVE_EXPO_MARK (expo);

  /* General case */
  {
    /* Declaration of the intermediary variable */
    mpfr_t t, te;
    mpfr_exp_t d;

    /* Declaration of the size variable */
    mpfr_prec_t Ny = MPFR_PREC(y);   /* target precision */
    mpfr_prec_t Nt;                  /* working precision */
    long int err;                  /* error */
    int sign = MPFR_SIGN (xt);
    MPFR_ZIV_DECL (loop);
    MPFR_GROUP_DECL (group);

    /* First check for BIG overflow of exp(2*x):
       For x > 0, exp(2*x) > 2^(2*x)
       If 2 ^(2*x) > 2^emax or x>emax/2, there is an overflow */
    if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emax/2) >= 0)) {
      /* initialise of intermediary variables
         since 'set_one' label assumes the variables have been
         initialize */
      MPFR_GROUP_INIT_2 (group, MPFR_PREC_MIN, t, te);
      goto set_one;
    }

    /* Compute the precision of intermediary variable */
    /* The optimal number of bits: see algorithms.tex */
    Nt = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 4;
    /* if x is small, there will be a cancellation in exp(2x)-1 */
    if (MPFR_GET_EXP (x) < 0)
      Nt += -MPFR_GET_EXP (x);

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

    MPFR_ZIV_INIT (loop, Nt);
    for (;;) {
      /* tanh = (exp(2x)-1)/(exp(2x)+1) */
      mpfr_mul_2ui (te, x, 1, MPFR_RNDN);  /* 2x */
      /* since x > 0, we can only have an overflow */
      mpfr_exp (te, te, MPFR_RNDN);        /* exp(2x) */
      if (MPFR_UNLIKELY (MPFR_IS_INF (te))) {
      set_one:
        inexact = MPFR_FROM_SIGN_TO_INT (sign);
        mpfr_set4 (y, __gmpfr_one, MPFR_RNDN, sign);
        if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG_SIGN (sign)))
          {
            inexact = -inexact;
            mpfr_nexttozero (y);
          }
        break;
      }
      d = MPFR_GET_EXP (te);              /* For Error calculation */
      mpfr_add_ui (t, te, 1, MPFR_RNDD);   /* exp(2x) + 1*/
      mpfr_sub_ui (te, te, 1, MPFR_RNDU);  /* exp(2x) - 1*/
      d = d - MPFR_GET_EXP (te);
      mpfr_div (t, te, t, MPFR_RNDN);      /* (exp(2x)-1)/(exp(2x)+1)*/

      /* Calculation of the error */
//.........这里部分代码省略.........