C++ MPFR_LIKELY函数代码示例

int
mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate)
{
  mp_ptr rp;
  mp_prec_t nbits;
  mp_size_t nlimbs;
  mp_size_t k; /* number of high zero limbs */
  mp_exp_t exp;
  int cnt;

  MPFR_CLEAR_FLAGS (rop);

  rp = MPFR_MANT (rop);
  nbits = MPFR_PREC (rop);
  nlimbs = MPFR_LIMB_SIZE (rop);
  MPFR_SET_POS (rop);

  /* Uniform non-normalized significand */
  _gmp_rand (rp, rstate, nlimbs * BITS_PER_MP_LIMB);

  /* If nbits isn't a multiple of BITS_PER_MP_LIMB, mask the low bits */
  cnt = nlimbs * BITS_PER_MP_LIMB - nbits;
  if (MPFR_LIKELY (cnt != 0))
    rp[0] &= ~MPFR_LIMB_MASK (cnt);

  /* Count the null significant limbs and remaining limbs */
  exp = 0;
  k = 0;
  while (nlimbs != 0 && rp[nlimbs - 1] == 0)
    {
      k ++;
      nlimbs --;
      exp -= BITS_PER_MP_LIMB;
    }

  if (MPFR_LIKELY (nlimbs != 0)) /* otherwise value is zero */
    {
      count_leading_zeros (cnt, rp[nlimbs - 1]);
      /* Normalization */
      if (mpfr_set_exp (rop, exp - cnt))
        {
          /* If the exponent is not in the current exponent range, we
             choose to return a NaN as this is probably a user error.
             Indeed this can happen only if the exponent range has been
             reduced to a very small interval and/or the precision is
             huge (very unlikely). */
          MPFR_SET_NAN (rop);
          __gmpfr_flags |= MPFR_FLAGS_NAN; /* Can't use MPFR_RET_NAN */
          return 1;
        }
      if (cnt != 0)
        mpn_lshift (rp + k, rp, nlimbs, cnt);
      if (k != 0)
        MPN_ZERO (rp, k);
    }
  else
    MPFR_SET_ZERO (rop);

  return 0;
}

int
mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate)
{
    mpfr_limb_ptr rp;
    mpfr_prec_t nbits;
    mp_size_t nlimbs;
    mp_size_t k; /* number of high zero limbs */
    mpfr_exp_t exp;
    int cnt;

    rp = MPFR_MANT (rop);
    nbits = MPFR_PREC (rop);
    nlimbs = MPFR_LIMB_SIZE (rop);
    MPFR_SET_POS (rop);
    cnt = nlimbs * GMP_NUMB_BITS - nbits;

    /* Uniform non-normalized significand */
    /* generate exactly nbits so that the random generator stays in the same
       state, independent of the machine word size GMP_NUMB_BITS */
    mpfr_rand_raw (rp, rstate, nbits);
    if (MPFR_LIKELY (cnt != 0)) /* this will put the low bits to zero */
        mpn_lshift (rp, rp, nlimbs, cnt);

    /* Count the null significant limbs and remaining limbs */
    exp = 0;
    k = 0;
    while (nlimbs != 0 && rp[nlimbs - 1] == 0)
    {
        k ++;
        nlimbs --;
        exp -= GMP_NUMB_BITS;
    }

    if (MPFR_LIKELY (nlimbs != 0)) /* otherwise value is zero */
    {
        count_leading_zeros (cnt, rp[nlimbs - 1]);
        /* Normalization */
        if (mpfr_set_exp (rop, exp - cnt))
        {
            /* If the exponent is not in the current exponent range, we
               choose to return a NaN as this is probably a user error.
               Indeed this can happen only if the exponent range has been
               reduced to a very small interval and/or the precision is
               huge (very unlikely). */
            MPFR_SET_NAN (rop);
            __gmpfr_flags |= MPFR_FLAGS_NAN; /* Can't use MPFR_RET_NAN */
            return 1;
        }
        if (cnt != 0)
            mpn_lshift (rp + k, rp, nlimbs, cnt);
        if (k != 0)
            MPN_ZERO (rp, k);
    }
    else
        MPFR_SET_ZERO (rop);

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

MPFR_HOT_FUNCTION_ATTR void
mpfr_mpz_init2 (mpz_t z, mp_bitcnt_t n)
{
  /* The condition on n is used below as the argument n will be ignored if
     the mpz_t is obtained from the MPFR stack of previously used mpz_t.
     Said otherwise, it z is expected to have a large size at the end, then
     it is better to allocate this size directly than to get a mpz_t of
     small size, with possibly several realloc's on it. But if n satisfies
     the condition and is larger than the stacked mpz_t, this may still
     yield useless realloc's. This is not ideal. We might consider to use
     mpz_init2 with the maximum size in mpfr_mpz_init to solve this issue. */
  if (MPFR_LIKELY (n_alloc > 0 && n <= MPFR_POOL_MAX_SIZE * GMP_NUMB_BITS))
    {
      /* Get a mpz_t from the MPFR stack of previously used mpz_t.
         It reduces memory pressure, and it allows to reuse
         a mpz_t that should be sufficiently big. */
      MPFR_ASSERTD (n_alloc <= numberof (mpz_tab));
      memcpy (z, &mpz_tab[--n_alloc], sizeof (mpz_t));
      SIZ(z) = 0;
    }
  else
    {
      /* Call the real GMP function */
      mpz_init2 (z, n);
    }
}

/*
 * Assuming {bp, bn} is an approximation of a non-singular number
 * with error at most equal to 2^(EXP(b)-err0) (`err0' bits of b are known)
 * of direction unknown, check if we can round b toward zero with
 * precision prec.
 */
int
mpfr_round_p (mp_limb_t *bp, mp_size_t bn, mpfr_exp_t err0, mpfr_prec_t prec)
{
  mpfr_prec_t err;
  mp_size_t k, n;
  mp_limb_t tmp, mask;
  int s;

  err = (mpfr_prec_t) bn * GMP_NUMB_BITS;
  if (MPFR_UNLIKELY (err0 <= 0 || (mpfr_uexp_t) err0 <= prec || prec >= err))
    return 0;  /* can't round */
  err = MIN (err, (mpfr_uexp_t) err0);

  k = prec / GMP_NUMB_BITS;
  s = GMP_NUMB_BITS - prec%GMP_NUMB_BITS;
  n = err / GMP_NUMB_BITS - k;

  MPFR_ASSERTD (n >= 0);
  MPFR_ASSERTD (bn > k);

  /* Check first limb */
  bp += bn-1-k;
  tmp = *bp--;
  mask = s == GMP_NUMB_BITS ? MP_LIMB_T_MAX : MPFR_LIMB_MASK (s);
  tmp &= mask;

  if (MPFR_LIKELY (n == 0))
    {
      /* prec and error are in the same limb */
      s = GMP_NUMB_BITS - err % GMP_NUMB_BITS;
      MPFR_ASSERTD (s < GMP_NUMB_BITS);
      tmp  >>= s;
      mask >>= s;
      return tmp != 0 && tmp != mask;
    }

int
mpfr_div_2si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode)
{
  int inexact;

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

  inexact = MPFR_UNLIKELY(y != x) ? mpfr_set (y, x, rnd_mode) : 0;

  if (MPFR_LIKELY( MPFR_IS_PURE_FP(y) ))
    {
      mp_exp_t exp = MPFR_GET_EXP (y);
      if (MPFR_UNLIKELY( n > 0 && (__gmpfr_emin > MPFR_EMAX_MAX - n ||
                                   exp < __gmpfr_emin + n)) )
        {
          if (rnd_mode == GMP_RNDN &&
              (__gmpfr_emin > MPFR_EMAX_MAX - (n - 1) ||
               exp < __gmpfr_emin + (n - 1) ||
               (inexact >= 0 && mpfr_powerof2_raw (y))))
            rnd_mode = GMP_RNDZ;
          return mpfr_underflow (y, rnd_mode, MPFR_SIGN(y));
        }

      if (MPFR_UNLIKELY(n < 0 && (__gmpfr_emax < MPFR_EMIN_MIN - n ||
                                  exp > __gmpfr_emax + n)) )
        return mpfr_overflow (y, rnd_mode, MPFR_SIGN(y));

      MPFR_SET_EXP (y, exp - n);
    }

  return inexact;
}

/*
 * Set f to z, choosing the smallest precision for f
 * so that z = f*(2^BPML)*zs*2^(RetVal)
 */
static int
set_z (mpfr_ptr f, mpz_srcptr z, mp_size_t *zs)
{
  mp_limb_t *p;
  mp_size_t s;
  int c;
  mp_prec_t pf;

  MPFR_ASSERTD (mpz_sgn (z) != 0);

  /* Remove useless ending 0 */
  for (p = PTR (z), s = *zs = ABS (SIZ (z)) ; *p == 0; p++, s--)
    MPFR_ASSERTD (s >= 0);

  /* Get working precision */
  count_leading_zeros (c, p[s-1]);
  pf = s * BITS_PER_MP_LIMB - c;
  if (pf < MPFR_PREC_MIN)
    pf = MPFR_PREC_MIN;
  mpfr_init2 (f, pf);

  /* Copy Mantissa */
  if (MPFR_LIKELY (c))
    mpn_lshift (MPFR_MANT (f), p, s, c);
  else
    MPN_COPY (MPFR_MANT (f), p, s);

  MPFR_SET_SIGN (f, mpz_sgn (z));
  MPFR_SET_EXP (f, 0);

  return -c;
}

/* Put in  rp[n..2n-1] an approximation of the n high limbs
   of {np, n} * {mp, n}. The error is less than n ulps of rp[n] (and the
   approximation is always less or equal to the truncated full product).

   Implements Algorithm ShortMul from [1].
*/
void
mpfr_mulhigh_n (mpfr_limb_ptr rp, mpfr_limb_srcptr np, mpfr_limb_srcptr mp,
                mp_size_t n)
{
  mp_size_t k;

  MPFR_ASSERTN (MPFR_MULHIGH_TAB_SIZE >= 8); /* so that 3*(n/4) > n/2 */
  k = MPFR_LIKELY (n < MPFR_MULHIGH_TAB_SIZE) ? mulhigh_ktab[n] : 3*(n/4);
  /* Algorithm ShortMul from [1] requires k >= (n+3)/2, which translates
     into k >= (n+4)/2 in the C language. */
  MPFR_ASSERTD (k == -1 || k == 0 || (k >= (n+4)/2 && k < n));
  if (k < 0)
    mpn_mul_basecase (rp, np, n, mp, n); /* result is exact, no error */
  else if (k == 0)
    mpfr_mulhigh_n_basecase (rp, np, mp, n); /* basecase error < n ulps */
  else if (n > MUL_FFT_THRESHOLD)
    mpn_mul_n (rp, np, mp, n); /* result is exact, no error */
  else
    {
      mp_size_t l = n - k;
      mp_limb_t cy;

      mpn_mul_n (rp + 2 * l, np + l, mp + l, k); /* fills rp[2l..2n-1] */
      mpfr_mulhigh_n (rp, np + k, mp, l);        /* fills rp[l-1..2l-1] */
      cy = mpn_add_n (rp + n - 1, rp + n - 1, rp + l - 1, l + 1);
      mpfr_mulhigh_n (rp, np, mp + k, l);        /* fills rp[l-1..2l-1] */
      cy += mpn_add_n (rp + n - 1, rp + n - 1, rp + l - 1, l + 1);
      mpn_add_1 (rp + n + l, rp + n + l, k, cy); /* propagate carry */
    }
}

/* Put in  rp[0..n] the n+1 low limbs of {np, n} * {mp, n}.
   Assume 2n limbs are allocated at rp. */
void
mpfr_mullow_n (mpfr_limb_ptr rp, mpfr_limb_srcptr np, mpfr_limb_srcptr mp,
               mp_size_t n)
{
  mp_size_t k;

  MPFR_ASSERTN (MPFR_MULHIGH_TAB_SIZE >= 8); /* so that 3*(n/4) > n/2 */
  k = MPFR_LIKELY (n < MPFR_MULHIGH_TAB_SIZE) ? mulhigh_ktab[n] : 3*(n/4);
  MPFR_ASSERTD (k == -1 || k == 0 || (2 * k >= n && k < n));
  if (k < 0)
    mpn_mul_basecase (rp, np, n, mp, n);
  else if (k == 0)
    mpfr_mullow_n_basecase (rp, np, mp, n);
  else if (n > MUL_FFT_THRESHOLD)
    mpn_mul_n (rp, np, mp, n);
  else
    {
      mp_size_t l = n - k;

      mpn_mul_n (rp, np, mp, k);                      /* fills rp[0..2k] */
      mpfr_mullow_n (rp + n, np + k, mp, l);          /* fills rp[n..n+2l] */
      mpn_add_n (rp + k, rp + k, rp + n, l + 1);
      mpfr_mullow_n (rp + n, np, mp + k, l);          /* fills rp[n..n+2l] */
      mpn_add_n (rp + k, rp + k, rp + n, l + 1);
    }
}

void
mpfr_custom_init_set (mpfr_ptr x, int kind, mp_exp_t exp,
                     mp_prec_t prec, void *mantissa)
{
  mpfr_kind_t t;
  int s;
  mp_exp_t e;

  if (kind >= 0)
    {
      t = (mpfr_kind_t) kind;
      s = MPFR_SIGN_POS;
    }
  else
    {
      t = (mpfr_kind_t) -kind;
      s = MPFR_SIGN_NEG;
    }
  MPFR_ASSERTD (t <= MPFR_REGULAR_KIND);
  e = MPFR_LIKELY (t == MPFR_REGULAR_KIND) ? exp :
    MPFR_UNLIKELY (t == MPFR_NAN_KIND) ? MPFR_EXP_NAN :
    MPFR_UNLIKELY (t == MPFR_INF_KIND) ? MPFR_EXP_INF : MPFR_EXP_ZERO;

  MPFR_PREC (x) = prec;
  MPFR_SET_SIGN (x, s);
  MPFR_EXP (x) = e;
  MPFR_MANT (x) = (mp_limb_t*) mantissa;
  return;
}

/* Put in  rp[n..2n-1] an approximation of the n high limbs
   of {np, n}^2. The error is less than n ulps of rp[n]. */
void
mpfr_sqrhigh_n (mpfr_limb_ptr rp, mpfr_limb_srcptr np, mp_size_t n)
{
  mp_size_t k;

  MPFR_ASSERTN (MPFR_SQRHIGH_TAB_SIZE > 2); /* ensures k < n */
  k = MPFR_LIKELY (n < MPFR_SQRHIGH_TAB_SIZE) ? sqrhigh_ktab[n]
    : (n+4)/2; /* ensures that k >= (n+3)/2 */
  MPFR_ASSERTD (k == -1 || k == 0 || (k >= (n+4)/2 && k < n));
  if (k < 0)
    /* we can't use mpn_sqr_basecase here, since it requires
       n <= SQR_KARATSUBA_THRESHOLD, where SQR_KARATSUBA_THRESHOLD
       is not exported by GMP */
    mpn_sqr_n (rp, np, n);
  else if (k == 0)
    mpfr_mulhigh_n_basecase (rp, np, np, n);
  else
    {
      mp_size_t l = n - k;
      mp_limb_t cy;

      mpn_sqr_n (rp + 2 * l, np + l, k);          /* fills rp[2l..2n-1] */
      mpfr_mulhigh_n (rp, np, np + k, l);         /* fills rp[l-1..2l-1] */
      /* {rp+n-1,l+1} += 2 * {rp+l-1,l+1} */
      cy = mpn_lshift (rp + l - 1, rp + l - 1, l + 1, 1);
      cy += mpn_add_n (rp + n - 1, rp + n - 1, rp + l - 1, l + 1);
      mpn_add_1 (rp + n + l, rp + n + l, k, cy); /* propagate carry */
    }
}

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_div (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t uu;
  mp_limb_t up[1];
  unsigned long cnt;

  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_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;
      MPFR_SET_EXP (uu, GMP_NUMB_BITS - cnt);
      return mpfr_div (y, uu, x, rnd_mode);
    }
  else /* u = 0, and x != 0 */
    {
      MPFR_SET_ZERO(y);         /* if u=0, then set y to 0 */
      MPFR_SET_SAME_SIGN(y, x); /* u considered as +0: sign(+0/x) = sign(x) */
      MPFR_RET(0);
    }
}

int
mpfr_check_range (mpfr_ptr x, int t, mp_rnd_t rnd_mode)
{
  if (MPFR_LIKELY( MPFR_IS_PURE_FP(x)) )
    { /* x is a non-zero FP */
      mp_exp_t exp = MPFR_EXP (x);  /* Do not use MPFR_GET_EXP */
      if (MPFR_UNLIKELY( exp < __gmpfr_emin) )
        {
          /* The following test is necessary because in the rounding to the
           * nearest mode, mpfr_underflow always rounds away from 0. In
           * this rounding mode, we need to round to 0 if:
           *   _ |x| < 2^(emin-2), or
           *   _ |x| = 2^(emin-2) and the absolute value of the exact
           *     result is <= 2^(emin-2).
           */
          if (rnd_mode == GMP_RNDN &&
              (exp + 1 < __gmpfr_emin ||
               (mpfr_powerof2_raw(x) &&
                (MPFR_IS_NEG(x) ? t <= 0 : t >= 0))))
            rnd_mode = GMP_RNDZ;
          return mpfr_underflow(x, rnd_mode, MPFR_SIGN(x));
        }
      if (MPFR_UNLIKELY( exp > __gmpfr_emax) )
        return mpfr_overflow(x, rnd_mode, MPFR_SIGN(x));
    }
  MPFR_RET (t);  /* propagate inexact ternary value, unlike most functions */
}

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 */
}