C++ MPFR_SAVE_EXPO_FREE函数代码示例

/* Convert an mpz_t to an mpfr_exp_t, restricted to
   the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]. */
mpfr_exp_t
mpfr_ubf_zexp2exp (mpz_ptr ez)
{
  mp_size_t n;
  mpfr_eexp_t e;
  mpfr_t d;
  int inex;
  MPFR_SAVE_EXPO_DECL (expo);

  n = ABSIZ (ez); /* limb size of ez */
  if (n == 0)
    return 0;

  MPFR_SAVE_EXPO_MARK (expo);
  mpfr_init2 (d, n * GMP_NUMB_BITS);
  MPFR_DBGRES (inex = mpfr_set_z (d, ez, MPFR_RNDN));
  MPFR_ASSERTD (inex == 0);
  e = mpfr_get_exp_t (d, MPFR_RNDZ);
  mpfr_clear (d);
  MPFR_SAVE_EXPO_FREE (expo);
  if (MPFR_UNLIKELY (e < MPFR_EXP_MIN))
    return MPFR_EXP_MIN;
  if (MPFR_UNLIKELY (e > MPFR_EXP_MAX))
    return MPFR_EXP_MAX;
  return e;
}

int
mpfr_sqrt_ui (mpfr_ptr r, unsigned long u, mpfr_rnd_t rnd_mode)
{
  if (u)
    {
      mpfr_t uu;
      mp_limb_t up[1];
      unsigned long 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 = (mp_limb_t) u << cnt;

      MPFR_SAVE_EXPO_MARK (expo);
      MPFR_SET_EXP (uu, GMP_NUMB_BITS - cnt);
      inex = mpfr_sqrt(r, uu, rnd_mode);
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_check_range(r, inex, rnd_mode);
    }
  else /* sqrt(0) = 0 */
    {
      MPFR_SET_ZERO(r);
      MPFR_SET_POS(r);
      MPFR_RET(0);
    }
}

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

  if (MPFR_LIKELY(u != 0) )  /* if u=0, do nothing */
    {
      mpfr_t uu;
      mp_limb_t up[1];
      unsigned long cnt;
      int inex;
      MPFR_SAVE_EXPO_DECL (expo);

      MPFR_TMP_INIT1 (up, uu, GMP_NUMB_BITS);
      MPFR_ASSERTD (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_add works even when uu is out-of-range. */
      MPFR_SAVE_EXPO_MARK (expo);
      MPFR_SET_EXP (uu, GMP_NUMB_BITS - cnt);
      inex = mpfr_add(y, x, uu, rnd_mode);
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_check_range(y, inex, rnd_mode);
    }
  else
    /* (unsigned long) 0 is assumed to be a real 0 (unsigned) */
    return mpfr_set (y, x, rnd_mode);
}

/* This function does not change the flags. */
static void
mpfr_get_zexp (mpz_ptr ez, mpfr_srcptr x)
{
  mpz_init (ez);

  if (MPFR_IS_UBF (x))
    mpz_set (ez, MPFR_ZEXP (x));
  else
    {
      mp_limb_t e_limb[MPFR_EXP_LIMB_SIZE];
      mpfr_t e;
      int inex;
      MPFR_SAVE_EXPO_DECL (expo);

      /* TODO: Once this has been tested, optimize based on whether
         _MPFR_EXP_FORMAT <= 3. */
      MPFR_TMP_INIT1 (e_limb, e, sizeof (mpfr_exp_t) * CHAR_BIT);
      MPFR_SAVE_EXPO_MARK (expo);
      MPFR_DBGRES (inex = mpfr_set_exp_t (e, MPFR_GET_EXP (x), MPFR_RNDN));
      MPFR_ASSERTD (inex == 0);
      MPFR_DBGRES (inex = mpfr_get_z (ez, e, MPFR_RNDN));
      MPFR_ASSERTD (inex == 0);
      MPFR_SAVE_EXPO_FREE (expo);
    }
}

int
mpfr_mul_d (mpfr_ptr a, mpfr_srcptr b, double c, mpfr_rnd_t rnd_mode)
{
  int inexact;
  mpfr_t d;
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_LOG_FUNC
    (("b[%Pu]=%.*Rg c=%.20g rnd=%d",
      mpfr_get_prec(b), mpfr_log_prec, b, c, rnd_mode),
     ("a[%Pu]=%.*Rg inexact=%d",
      mpfr_get_prec (a), mpfr_get_prec, a, inexact));

  MPFR_SAVE_EXPO_MARK (expo);

  mpfr_init2 (d, IEEE_DBL_MANT_DIG);
  inexact = mpfr_set_d (d, c, rnd_mode);
  MPFR_ASSERTN (inexact == 0);

  mpfr_clear_flags ();
  inexact = mpfr_mul (a, b, d, rnd_mode);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);

  mpfr_clear(d);
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (a, inexact, rnd_mode);
}

int
mpfr_d_div (mpfr_ptr a, double b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
  int inexact;
  mpfr_t d;
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_LOG_FUNC (
    ("b=%.20g c[%Pu]=%*.Rg rnd=%d", b, mpfr_get_prec (c), mpfr_log_prec, c, rnd_mode),
    ("a[%Pu]=%*.Rg", mpfr_get_prec (a), mpfr_log_prec, a));

  MPFR_SAVE_EXPO_MARK (expo);

  mpfr_init2 (d, IEEE_DBL_MANT_DIG);
  inexact = mpfr_set_d (d, b, rnd_mode);
  MPFR_ASSERTN (inexact == 0);

  mpfr_clear_flags ();
  inexact = mpfr_div (a, d, c, rnd_mode);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);

  mpfr_clear(d);
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (a, inexact, rnd_mode);
}

int
mpfr_mul_d (mpfr_ptr a, mpfr_srcptr b, double c, mpfr_rnd_t rnd_mode)
{
  int inexact;
  mpfr_t d;
  mp_limb_t tmp_man[MPFR_LIMBS_PER_DOUBLE];
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_LOG_FUNC
    (("b[%Pu]=%.*Rg c=%.20g rnd=%d",
      mpfr_get_prec(b), mpfr_log_prec, b, c, rnd_mode),
     ("a[%Pu]=%.*Rg inexact=%d",
      mpfr_get_prec (a), mpfr_log_prec, a, inexact));

  MPFR_SAVE_EXPO_MARK (expo);

  MPFR_TMP_INIT1(tmp_man, d, IEEE_DBL_MANT_DIG);
  inexact = mpfr_set_d (d, c, rnd_mode);
  MPFR_ASSERTD (inexact == 0);

  MPFR_CLEAR_FLAGS ();
  inexact = mpfr_mul (a, b, d, rnd_mode);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);

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

int
mpfr_sub_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode)
{
  if (MPFR_LIKELY (u != 0))  /* if u=0, do nothing */
    {
      mpfr_t uu;
      mp_limb_t up[1];
      unsigned long cnt;
      int inex;

      MPFR_SAVE_EXPO_DECL (expo);

      MPFR_TMP_INIT1 (up, uu, BITS_PER_MP_LIMB);
      MPFR_ASSERTN (u == (mp_limb_t) u);
      count_leading_zeros (cnt, (mp_limb_t) u);
      *up = (mp_limb_t) u << cnt;

      /* Optimization note: Exponent save/restore operations may be
         removed if mpfr_sub works even when uu is out-of-range. */
      MPFR_SAVE_EXPO_MARK (expo);
      MPFR_SET_EXP (uu, BITS_PER_MP_LIMB - cnt);
      inex = mpfr_sub (y, x, uu, rnd_mode);
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_check_range (y, inex, rnd_mode);
    }
  else
    return mpfr_set (y, x, rnd_mode);
}

int
mpfr_ui_pow (mpfr_ptr y, unsigned long int n, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t t;
  int inexact;
  mp_limb_t tmp_mant[(sizeof (n) - 1) / sizeof (mp_limb_t) + 1];
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_SAVE_EXPO_MARK (expo);
  MPFR_TMP_INIT1(tmp_mant, t, sizeof(n) * CHAR_BIT);
  inexact = mpfr_set_ui (t, n, MPFR_RNDN);
  MPFR_ASSERTD (inexact == 0);
  inexact = mpfr_pow (y, t, x, rnd_mode);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (y, inexact, rnd_mode);
}

int
mpfr_get_z (mpz_ptr z, mpfr_srcptr f, mpfr_rnd_t rnd)
{
  int inex;
  mpfr_t r;
  mpfr_exp_t exp;
  MPFR_SAVE_EXPO_DECL (expo);

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (f)))
    {
      if (MPFR_UNLIKELY (MPFR_NOTZERO (f)))
        MPFR_SET_ERANGEFLAG ();
      mpz_set_ui (z, 0);
      /* The ternary value is 0 even for infinity. Giving the rounding
         direction in this case would not make much sense anyway, and
         the direction would not necessarily match rnd. */
      return 0;
    }

  MPFR_SAVE_EXPO_MARK (expo);

  exp = MPFR_GET_EXP (f);
  /* if exp <= 0, then |f|<1, thus |o(f)|<=1 */
  MPFR_ASSERTN (exp < 0 || exp <= MPFR_PREC_MAX);
  mpfr_init2 (r, (exp < (mpfr_exp_t) MPFR_PREC_MIN ?
                  MPFR_PREC_MIN : (mpfr_prec_t) exp));
  inex = mpfr_rint (r, f, rnd);
  MPFR_ASSERTN (inex != 1 && inex != -1); /* integral part of f is
                                             representable in r */
  MPFR_ASSERTN (MPFR_IS_FP (r));

  /* The flags from mpfr_rint are the wanted ones. In particular,
     it sets the inexact flag when necessary. */
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);

  exp = mpfr_get_z_2exp (z, r);
  if (exp >= 0)
    mpz_mul_2exp (z, z, exp);
  else
    mpz_fdiv_q_2exp (z, z, -exp);
  mpfr_clear (r);

  MPFR_SAVE_EXPO_FREE (expo);

  return inex;
}

static int
pi_div_2ui (mpfr_ptr dest, int i, int neg, mpfr_rnd_t rnd_mode)
{
  int inexact;
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_SAVE_EXPO_MARK (expo);
  if (neg) /* -PI/2^i */
    {
      inexact = - mpfr_const_pi (dest, MPFR_INVERT_RND (rnd_mode));
      MPFR_CHANGE_SIGN (dest);
    }
  else /* PI/2^i */
    {
      inexact = mpfr_const_pi (dest, rnd_mode);
    }
  mpfr_div_2ui (dest, dest, i, rnd_mode);  /* exact */
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (dest, inexact, rnd_mode);
}

int
mpfr_cmp_d (mpfr_srcptr b, double d)
{
  mpfr_t tmp;
  int res;
  mp_limb_t tmp_man[MPFR_LIMBS_PER_DOUBLE];
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_SAVE_EXPO_MARK (expo);

  MPFR_TMP_INIT1(tmp_man, tmp, IEEE_DBL_MANT_DIG);
  res = mpfr_set_d (tmp, d, MPFR_RNDN);
  MPFR_ASSERTD (res == 0);

  MPFR_CLEAR_FLAGS ();
  res = mpfr_cmp (b, tmp);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);

  MPFR_SAVE_EXPO_FREE (expo);
  return res;
}

int
mpfr_rint_trunc (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
{
  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
    return mpfr_set (r, u, rnd_mode);
  else
    {
      mpfr_t tmp;
      int inex;
      MPFR_SAVE_EXPO_DECL (expo);

      MPFR_SAVE_EXPO_MARK (expo);
      mpfr_init2 (tmp, MPFR_PREC (u));
      /* trunc(u) is always representable in tmp */
      mpfr_trunc (tmp, u);
      inex = mpfr_set (r, tmp, rnd_mode);
      mpfr_clear (tmp);
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_check_range (r, inex, rnd_mode);
    }
}

int
mpfr_div_d (mpfr_ptr a, mpfr_srcptr b, double c, mp_rnd_t rnd_mode)
{
  int inexact;
  mpfr_t d;
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_LOG_FUNC (("b[%#R]=%R c%.20g rnd=%d", b, b, c, rnd_mode),
                 ("a[%#R]=%R", a, a));

  MPFR_SAVE_EXPO_MARK (expo);

  mpfr_init2 (d, IEEE_DBL_MANT_DIG);
  inexact = mpfr_set_d (d, c, rnd_mode);
  MPFR_ASSERTN (inexact == 0);

  mpfr_clear_flags ();
  inexact = mpfr_div (a, b, d, rnd_mode);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);

  mpfr_clear(d);
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (a, inexact, rnd_mode);
}

int
mpfr_rint_floor (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
{
  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
    return mpfr_set (r, u, rnd_mode);
  else
    {
      mpfr_t tmp;
      int inex;
      MPFR_SAVE_EXPO_DECL (expo);
      MPFR_BLOCK_DECL (flags);

      MPFR_SAVE_EXPO_MARK (expo);
      mpfr_init2 (tmp, MPFR_PREC (u));
      /* floor(u) is representable in tmp unless an overflow occurs */
      MPFR_BLOCK (flags, mpfr_floor (tmp, u));
      inex = (MPFR_OVERFLOW (flags)
              ? mpfr_overflow (r, rnd_mode, MPFR_SIGN_NEG)
              : mpfr_set (r, tmp, rnd_mode));
      mpfr_clear (tmp);
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_check_range (r, inex, rnd_mode);
    }
}