C++ MPFR_IS_SINGULAR函数代码示例

static void
test_cmp_z (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
{
  mpfr_t x, z;
  mpz_t  y;
  mpfr_prec_t p;
  int res1, res2;
  int n;

  mpfr_init (x);
  mpfr_init2 (z, MPFR_PREC_MIN);
  mpz_init (y);

  /* check the erange flag when x is NaN */
  mpfr_set_nan (x);
  mpz_set_ui (y, 17);
  mpfr_clear_erangeflag ();
  res1 = mpfr_cmp_z (x, y);
  if (res1 != 0 || mpfr_erangeflag_p () == 0)
    {
      printf ("Error for mpfr_cmp_z (NaN, 17)\n");
      printf ("Return value: expected 0, got %d\n", res1);
      printf ("Erange flag: expected set, got %d\n", mpfr_erangeflag_p ());
      exit (1);
    }

  for(p=pmin ; p < pmax ; p++)
    {
      mpfr_set_prec (x, p);
      for ( n = 0; n < nmax ; n++)
        {
          mpfr_urandomb (x, RANDS);
          mpz_urandomb  (y, RANDS, 1024);
          if (!MPFR_IS_SINGULAR (x))
            {
              mpfr_sub_z (z, x, y, MPFR_RNDN);
              res1 = mpfr_sgn (z);
              res2 = mpfr_cmp_z (x, y);
              if (res1 != res2)
                {
                  printf("Error for mpfr_cmp_z: res=%d sub_z gives %d\n",
                         res2, res1);
                  exit (1);
                }
            }
        }
    }
  mpz_clear (y);
  mpfr_clear (x);
  mpfr_clear (z);
}

double
mpfr_get_d_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode)
{
  double ret;
  mpfr_exp_t exp;
  mpfr_t tmp;

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
    {
      int negative;
      *expptr = 0;
      if (MPFR_IS_NAN (src))
        return MPFR_DBL_NAN;
      negative = MPFR_IS_NEG (src);
      if (MPFR_IS_INF (src))
        return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
      MPFR_ASSERTD (MPFR_IS_ZERO(src));
      return negative ? DBL_NEG_ZERO : 0.0;
    }

  tmp[0] = *src;        /* Hack copy mpfr_t */
  MPFR_SET_EXP (tmp, 0);
  ret = mpfr_get_d (tmp, rnd_mode);

  if (MPFR_IS_PURE_FP(src))
    {
      exp = MPFR_GET_EXP (src);

      /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
      if (ret == 1.0)
        {
          ret = 0.5;
          exp++;
        }
      else if (ret == -1.0)
        {
          ret = -0.5;
          exp++;
        }

      MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
                    || (ret <= -0.5 && ret > -1.0));
      MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
    }
  else
    exp = 0;

  *expptr = exp;
  return ret;
}

/* special code for IEEE 754 little-endian extended format */
long double
mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
    mpfr_long_double_t ld;
    mpfr_t tmp;
    int inex;
    MPFR_SAVE_EXPO_DECL (expo);

    MPFR_SAVE_EXPO_MARK (expo);

    mpfr_init2 (tmp, MPFR_LDBL_MANT_DIG);
    inex = mpfr_set (tmp, x, rnd_mode);

    mpfr_set_emin (-16382-63);
    mpfr_set_emax (16384);
    mpfr_subnormalize (tmp, mpfr_check_range (tmp, inex, rnd_mode), rnd_mode);
    mpfr_prec_round (tmp, 64, MPFR_RNDZ); /* exact */
    if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tmp)))
        ld.ld = (long double) mpfr_get_d (tmp, rnd_mode);
    else
    {
        mp_limb_t *tmpmant;
        mpfr_exp_t e, denorm;

        tmpmant = MPFR_MANT (tmp);
        e = MPFR_GET_EXP (tmp);
        /* The smallest positive normal number is 2^(-16382), which is
           0.5*2^(-16381) in MPFR, thus any exponent <= -16382 corresponds to a
           subnormal number. The smallest positive subnormal number is 2^(-16445)
           which is 0.5*2^(-16444) in MPFR thus 0 <= denorm <= 63. */
        denorm = MPFR_UNLIKELY (e <= -16382) ? - e - 16382 + 1 : 0;
        MPFR_ASSERTD (0 <= denorm && denorm < 64);
#if GMP_NUMB_BITS >= 64
        ld.s.manl = (tmpmant[0] >> denorm);
        ld.s.manh = (tmpmant[0] >> denorm) >> 32;
#elif GMP_NUMB_BITS == 32
        if (MPFR_LIKELY (denorm == 0))
        {
            ld.s.manl = tmpmant[0];
            ld.s.manh = tmpmant[1];
        }
        else if (denorm < 32)
        {
            ld.s.manl = (tmpmant[0] >> denorm) | (tmpmant[1] << (32 - denorm));
            ld.s.manh = tmpmant[1] >> denorm;
        }
        else /* 32 <= denorm < 64 */
        {

int
(mpfr_sgn) (mpfr_srcptr a)
{
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
    {
      if (MPFR_LIKELY (MPFR_IS_ZERO (a)))
        return 0;
      if (MPFR_UNLIKELY (MPFR_IS_NAN (a)))
        {
          MPFR_SET_ERANGEFLAG ();
          return 0;
        }
      /* Remains infinity, handled by the return below. */
    }
  return MPFR_INT_SIGN (a);
}

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

/* Return 1 if y is an odd integer, 0 otherwise.
   Assumes y is not singular. */
int
mpfr_odd_p (mpfr_srcptr y)
{
  mpfr_exp_t expo;
  mpfr_prec_t prec;
  mp_size_t yn;
  mp_limb_t *yp;

  /* NAN, INF or ZERO are not allowed */
  MPFR_ASSERTD (!MPFR_IS_SINGULAR (y));

  expo = MPFR_GET_EXP (y);
  if (expo <= 0)
    return 0;  /* |y| < 1 and not 0 */

  prec = MPFR_PREC(y);
  if ((mpfr_prec_t) expo > prec)
    return 0;  /* y is a multiple of 2^(expo-prec), thus not odd */

  /* 0 < expo <= prec:
     y = 1xxxxxxxxxt.zzzzzzzzzzzzzzzzzz[000]
          expo bits   (prec-expo) bits

     We have to check that:
     (a) the bit 't' is set
     (b) all the 'z' bits are zero
  */

  prec = MPFR_PREC2LIMBS (prec) * GMP_NUMB_BITS - expo;
  /* number of z+0 bits */

  yn = prec / GMP_NUMB_BITS;
  MPFR_ASSERTN(yn >= 0);
  /* yn is the index of limb containing the 't' bit */

  yp = MPFR_MANT(y);
  /* if expo is a multiple of GMP_NUMB_BITS, t is bit 0 */
  if (expo % GMP_NUMB_BITS == 0 ? (yp[yn] & 1) == 0
      : yp[yn] << ((expo % GMP_NUMB_BITS) - 1) != MPFR_LIMB_HIGHBIT)
    return 0;
  while (--yn >= 0)
    if (yp[yn] != 0)
      return 0;
  return 1;
}

int
mpfr_erfc (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
{
  int inex;
  mpfr_t tmp;
  mpfr_exp_t te, err;
  mpfr_prec_t prec;
  mpfr_exp_t emin = mpfr_get_emin ();
  MPFR_SAVE_EXPO_DECL (expo);
  MPFR_ZIV_DECL (loop);

  MPFR_LOG_FUNC
    (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
     ("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;
        }
      /* erfc(+inf) = 0+, erfc(-inf) = 2 erfc (0) = 1 */
      else if (MPFR_IS_INF (x))
        return mpfr_set_ui (y, MPFR_IS_POS (x) ? 0 : 2, rnd);
      else
        return mpfr_set_ui (y, 1, rnd);
    }

  if (MPFR_SIGN (x) > 0)
    {
      /* by default, emin = 1-2^30, thus the smallest representable
         number is 1/2*2^emin = 2^(-2^30):
         for x >= 27282, erfc(x) < 2^(-2^30-1), and
         for x >= 1787897414, erfc(x) < 2^(-2^62-1).
      */
      if ((emin >= -1073741823 && mpfr_cmp_ui (x, 27282) >= 0) ||
          mpfr_cmp_ui (x, 1787897414) >= 0)
        {
          /* May be incorrect if MPFR_EMAX_MAX >= 2^62. */
          MPFR_ASSERTN ((MPFR_EMAX_MAX >> 31) >> 31 == 0);
          return mpfr_underflow (y, (rnd == MPFR_RNDN) ? MPFR_RNDZ : rnd, 1);
        }
    }

int
mpfr_ui_sub (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 (u == 0))
    return mpfr_neg (y, x, rnd_mode);

  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 = -Inf and u - -Inf = +Inf  */
          MPFR_SET_INF(y);
          MPFR_SET_OPPOSITE_SIGN(y,x);
          MPFR_RET(0); /* +/-infinity is exact */
        }
      else /* x is zero */
        /* u - 0 = u */
        return mpfr_set_ui(y, u, rnd_mode);
    }
  else
    {
      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_SET_EXP (uu, GMP_NUMB_BITS - cnt);
      return mpfr_sub (y, uu, x, rnd_mode);
    }
}

/* We can't use fits_u.h <= mpfr_cmp_ui */
int
mpfr_fits_uintmax_p (mpfr_srcptr f, mpfr_rnd_t rnd)
{
  mpfr_exp_t e;
  int prec;
  uintmax_t s;
  mpfr_t x;
  int res;

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (f)))
    /* Zero always fit */
    return MPFR_IS_ZERO (f) ? 1 : 0;
  else if (MPFR_IS_NEG (f))
    /* Negative numbers don't fit */
    return 0;
  /* now it fits if
     (a) f <= MAXIMUM
     (b) round(f, prec(slong), rnd) <= MAXIMUM */

  e = MPFR_GET_EXP (f);

  /* first compute prec(MAXIMUM); fits in an int */
  for (s = MPFR_UINTMAX_MAX, prec = 0; s != 0; s /= 2, prec ++);

  /* MAXIMUM needs prec bits, i.e. MAXIMUM = 2^prec - 1 */

  /* if e <= prec - 1, then f < 2^(prec-1) < MAXIMUM */
  if (e <= prec - 1)
    return 1;

  /* if e >= prec + 1, then f >= 2^prec > MAXIMUM */
  if (e >= prec + 1)
    return 0;

  MPFR_ASSERTD (e == prec);

  /* hard case: first round to prec bits, then check */
  mpfr_init2 (x, prec);
  mpfr_set (x, f, rnd);
  res = MPFR_GET_EXP (x) == e;
  mpfr_clear (x);
  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;
      unsigned int saved_flags = __gmpfr_flags;

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

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

long double
mpfr_get_ld (mpfr_srcptr x, mp_rnd_t rnd_mode)
{
  mpfr_long_double_t ld;
  mpfr_t tmp;
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_SAVE_EXPO_MARK (expo);
  mpfr_set_emin (-16382-63);
  mpfr_set_emax (16383);

  mpfr_init2 (tmp, MPFR_LDBL_MANT_DIG);
  mpfr_subnormalize(tmp, mpfr_set (tmp, x, rnd_mode), rnd_mode);
  mpfr_prec_round (tmp, 64, GMP_RNDZ); /* exact */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tmp)))
    ld.ld = (long double) mpfr_get_d (tmp, rnd_mode);
  else
    {
      mp_limb_t *tmpmant;
      mp_exp_t e, denorm;

      tmpmant = MPFR_MANT (tmp);
      e = MPFR_GET_EXP (tmp);
      denorm = MPFR_UNLIKELY (e < -16382) ? - e - 16382 + 1 : 0;
#if BITS_PER_MP_LIMB >= 64
      ld.s.manl = (tmpmant[0] >> denorm);
      ld.s.manh = (tmpmant[0] >> denorm) >> 32;
#elif BITS_PER_MP_LIMB == 32
      if (MPFR_LIKELY (denorm == 0))
        {
          ld.s.manl = tmpmant[0];
          ld.s.manh = tmpmant[1];
        }
      else if (denorm < 32)
        {
          ld.s.manl = (tmpmant[0] >> denorm) | (tmpmant[1] << (32 - denorm));
          ld.s.manh = tmpmant[1] >> denorm;
        }
      else /* 32 <= denorm <= 64 */
        {

static void
test_cmp_f (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
{
  mpfr_t x, z;
  mpf_t  y;
  mpfr_prec_t p;
  int res1, res2;
  int n;

  mpfr_init (x);
  mpfr_init2 (z, pmax+GMP_NUMB_BITS);
  mpf_init2 (y, MPFR_PREC_MIN);
  for(p=pmin ; p < pmax ; p+=3)
    {
      mpfr_set_prec (x, p);
      mpf_set_prec (y, p);
      for ( n = 0; n < nmax ; n++)
        {
          mpfr_urandomb (x, RANDS);
          mpf_urandomb  (y, RANDS, p);
          if (!MPFR_IS_SINGULAR (x))
            {
              mpfr_set_f (z, y, MPFR_RNDN);
              mpfr_sub   (z, x, z, MPFR_RNDN);
              res1 = mpfr_sgn (z);
              res2 = mpfr_cmp_f (x, y);
              if (res1 != res2)
                {
                  printf("Error for mpfr_cmp_f: res=%d sub gives %d\n",
                         res2, res1);
                  exit (1);
                }
            }
        }
    }
  mpf_clear (y);
  mpfr_clear (x);
  mpfr_clear (z);
}

int
mpfr_ui_sub (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mp_rnd_t rnd_mode)
{
  mpfr_t uu;
  mp_limb_t up[1];
  unsigned long cnt;

  if (MPFR_UNLIKELY (u == 0))
    return mpfr_neg (y, x, rnd_mode);

  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 = -Inf and u - -Inf = +Inf  */
          MPFR_SET_INF(y);
          MPFR_SET_OPPOSITE_SIGN(y,x);
          MPFR_RET(0); /* +/-infinity is exact */
        }
      else /* x is zero */
        /* u - 0 = u */
        return mpfr_set_ui(y, u, rnd_mode);
    }
  else
    {
      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;
      MPFR_SET_EXP (uu, BITS_PER_MP_LIMB - cnt);
      return mpfr_sub (y, uu, x, rnd_mode);
    }
}

void
mpfr_extract (mpz_ptr y, mpfr_srcptr p, unsigned int i)
{
  unsigned long two_i = 1UL << i;
  unsigned long two_i_2 = i ? two_i / 2 : 1;
  mp_size_t size_p = MPFR_LIMB_SIZE (p);

  /* as 0 <= |p| < 1, we don't have to care with infinities, NaN, ... */
  MPFR_ASSERTD (!MPFR_IS_SINGULAR (p));

  _mpz_realloc (y, two_i_2);
  if ((mpfr_uexp_t) size_p < two_i)
    {
      MPN_ZERO (PTR(y), two_i_2);
      if ((mpfr_uexp_t) size_p >= two_i_2)
        MPN_COPY (PTR(y) + two_i - size_p, MPFR_MANT(p), size_p - two_i_2);
    }
  else
    MPN_COPY (PTR(y), MPFR_MANT(p) + size_p - two_i, two_i_2);

  MPN_NORMALIZE (PTR(y), two_i_2);
  SIZ(y) = (MPFR_IS_NEG (p)) ? -two_i_2 : two_i_2;
}