C++ MPFR_MANT函数代码示例

/*
 * x : OUT : MPFR number extracted from the binary buffer, should have the same
 *           precision than the number in the binary format
 * buffer : IN : limb of the MPFR number x in a binary format,
 * buffer_size : IN/OUT : size of the buffer => size used in the buffer
 * return 0 if successful
 */
static int
mpfr_fpif_read_limbs (mpfr_t x, unsigned char *buffer, size_t *buffer_size)
{
  mpfr_prec_t precision;
  size_t nb_byte;
  size_t mp_bytes_per_limb;
  size_t nb_partial_byte;
  size_t i, j;

  precision = mpfr_get_prec (x);
  nb_byte = (precision + 7) >> 3;
  mp_bytes_per_limb = mp_bits_per_limb >> 3;
  nb_partial_byte = nb_byte % mp_bytes_per_limb;

  if ((buffer == NULL) || (*buffer_size < nb_byte))
    {
      *buffer_size = 0;
      return 1;
    }
  *buffer_size = nb_byte;

  if (nb_partial_byte > 0)
    {
      memset (MPFR_MANT(x), 0, sizeof(mp_limb_t));
      getBigEndianData ((unsigned char*) MPFR_MANT(x), buffer,
                        sizeof(mp_limb_t), nb_partial_byte);
    }
  for (i = nb_partial_byte, j = (nb_partial_byte == 0) ? 0 : 1; i < nb_byte;
       i += mp_bytes_per_limb, j++)
    getLittleEndianData ((unsigned char*) (MPFR_MANT(x) + j), buffer + i,
                         sizeof(mp_limb_t), sizeof(mp_limb_t));

  return 0;
}

void
mpfr_swap (mpfr_ptr u, mpfr_ptr v)
{
  mpfr_prec_t p1, p2;
  mpfr_sign_t s1, s2;
  mpfr_exp_t e1, e2;
  mp_limb_t *m1, *m2;

  p1 = MPFR_PREC(u);
  p2 = MPFR_PREC(v);
  MPFR_PREC(v) = p1;
  MPFR_PREC(u) = p2;

  s1 = MPFR_SIGN(u);
  s2 = MPFR_SIGN(v);
  MPFR_SIGN(v) = s1;
  MPFR_SIGN(u) = s2;

  e1 = MPFR_EXP(u);
  e2 = MPFR_EXP(v);
  MPFR_EXP(v) = e1;
  MPFR_EXP(u) = e2;

  m1 = MPFR_MANT(u);
  m2 = MPFR_MANT(v);
  MPFR_MANT(v) = m1;
  MPFR_MANT(u) = m2;
}

/*
 * buffer : OUT : store the limb of the MPFR number x in a binary format,
 *                can be null (may be reallocated if too small)
 * buffer_size : IN/OUT : size of the buffer => size used in the buffer
 * x : IN : MPFR number
 * return pointer to a buffer storing the limb of the MPFR number x in a binary
 *        format
 */
static unsigned char*
mpfr_fpif_store_limbs (unsigned char *buffer, size_t *buffer_size, mpfr_t x)
{
  unsigned char *result;
  mpfr_prec_t precision;
  size_t nb_byte;
  size_t nb_limb, mp_bytes_per_limb;
  size_t nb_partial_byte;
  size_t i, j;

  precision = mpfr_get_prec (x);
  nb_byte = (precision + 7) >> 3;
  mp_bytes_per_limb = mp_bits_per_limb >> 3;
  nb_partial_byte = nb_byte % mp_bytes_per_limb;
  nb_limb = (nb_byte + mp_bytes_per_limb - 1) / mp_bytes_per_limb;

  result = buffer;
  ALLOC_RESULT(result, buffer_size, nb_byte);

  putBigEndianData (result, (unsigned char*) MPFR_MANT(x),
                    sizeof(mp_limb_t), nb_partial_byte);
  for (i = nb_partial_byte, j = (nb_partial_byte == 0) ? 0 : 1; j < nb_limb;
       i += mp_bytes_per_limb, j++)
    putLittleEndianData (result + i, (unsigned char*) (MPFR_MANT(x) + j),
                         sizeof(mp_limb_t), sizeof(mp_limb_t));

  return result;
}

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

static void *
start_routine (void *arg)
{
  mpfr_prec_t p;
  mpfr_t x;
  mpfr_prec_t inc = *(int *) arg;
  mp_limb_t *m;

  for (p = 100; p < 20000; p += 64 + 100 * (inc % 10))
    {
      mpfr_init2 (x, p);
      m = MPFR_MANT (x);
      mpfr_const_pi (x, MPFR_RNDD);
      mpfr_prec_round (x, 53, MPFR_RNDD);
      if (mpfr_cmp_str1 (x, "3.141592653589793116"))
        {
          printf ("mpfr_const_pi failed with threading\n");
          mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); putchar('\n');
          exit (1);
        }
      /* Check that no reallocation has been performed */
      MPFR_ASSERTN (m == MPFR_MANT (x));
      mpfr_clear (x);
    }

  pthread_exit (NULL);
}

int
mpfr_cmpabs (mpfr_srcptr b, mpfr_srcptr c)
{
  mp_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;
    }

  be = MPFR_GET_EXP (b);
  ce = MPFR_GET_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;
}

/* set bit n of x to b, where bit 0 is the most significant one */
static void
set_bit (mpfr_t x, unsigned int n, int b)
{
  unsigned l;
  mp_size_t xn;

  xn = (MPFR_PREC(x) - 1) / mp_bits_per_limb;
  l = n / mp_bits_per_limb;
  n %= mp_bits_per_limb;
  n = mp_bits_per_limb - 1 - n;
  if (b)
    MPFR_MANT(x)[xn - l] |= (mp_limb_t) 1 << n;
  else
    MPFR_MANT(x)[xn - l] &= ~((mp_limb_t) 1 << n);
}

void
mpfr_random2 (mpfr_ptr x, mp_size_t size, mp_exp_t exp)
{
  mp_size_t xn;
  unsigned long cnt;
  mp_ptr xp = MPFR_MANT(x), yp[1];
  mp_size_t prec = (MPFR_PREC(x) - 1)/BITS_PER_MP_LIMB; 

  MPFR_CLEAR_FLAGS(x);
  xn = ABS (size);
  if (xn != 0)
    {
      if (xn > prec + 1)
	xn = prec + 1;

      mpn_random2 (xp, xn);
    }

  if (exp != 0) {
    /* use mpn_random instead of random since that function is not
       available on all platforms (for example HPUX, DEC OSF, ...) */
    mpn_random ((mp_limb_t*) yp, 1);
    exp = (mp_exp_t) yp[0] % (2 * exp) - exp;
  }

  count_leading_zeros(cnt, xp[xn - 1]); 
  if (cnt) mpn_lshift(xp, xp, xn, cnt); 
  MPFR_EXP(x) = exp-cnt; 
  cnt = xn*BITS_PER_MP_LIMB - prec; 
  /* cnt is the number of non significant bits in the low limb */
  xp[0] &= ~((MP_LIMB_T_ONE << cnt) - MP_LIMB_T_ONE);
}

MPFR_HOT_FUNCTION_ATTR void
mpfr_clear (mpfr_ptr m)
{
  (*__gmp_free_func) (MPFR_GET_REAL_PTR (m),
                      MPFR_MALLOC_SIZE (MPFR_GET_ALLOC_SIZE (m)));
  MPFR_MANT (m) = (mp_limb_t *) 0;
}

/* sets x to x+sign(x)*ulp(x) */
int
mpfr_add_one_ulp (mpfr_ptr x, mp_rnd_t rnd_mode)
{
  mp_size_t xn;
  int sh;
  mp_limb_t *xp;

  if (MPFR_IS_NAN(x))
    MPFR_RET_NAN;

  if (MPFR_IS_INF(x) || MPFR_IS_ZERO(x))
    return 0;

  xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB;
  sh = xn * BITS_PER_MP_LIMB - MPFR_PREC(x);
  xp = MPFR_MANT(x);
  if (mpn_add_1 (xp, xp, xn, MP_LIMB_T_ONE << sh)) /* got 1.0000... */
    {
      mp_exp_t exp = MPFR_EXP(x);
      if (exp == __mpfr_emax)
        return mpfr_set_overflow(x, rnd_mode, MPFR_SIGN(x));
      else
        {
          MPFR_EXP(x)++;
          xp[xn-1] = GMP_LIMB_HIGHBIT;
        }
    }
  return 0;
}

int
mpfr_set_overflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign)
{
  int inex;

  MPFR_CLEAR_FLAGS(x);
  if ((rnd_mode == GMP_RNDU && sign < 0)
   || (rnd_mode == GMP_RNDD && sign > 0))
    {
      mp_size_t xn, i;
      int sh;
      mp_limb_t *xp;

      MPFR_EXP(x) = __mpfr_emax;
      xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB;
      sh = xn * BITS_PER_MP_LIMB - MPFR_PREC(x);
      xp = MPFR_MANT(x);
      xp[0] = MP_LIMB_T_MAX << sh;
      for (i = 1; i < xn; i++)
        xp[i] = MP_LIMB_T_MAX;
      inex = -1;
    }
  else
    {
      MPFR_SET_INF(x);
      inex = 1;
    }
  if (MPFR_SIGN(x) != sign)
    MPFR_CHANGE_SIGN(x);
  __mpfr_flags |= MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
  return sign > 0 ? inex : -inex;
}

int
mpfr_set_underflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign)
{
  int inex;

  MPFR_CLEAR_FLAGS(x);
  if ((rnd_mode == GMP_RNDU && sign > 0)
   || (rnd_mode == GMP_RNDD && sign < 0))
    {
      mp_size_t xn;
      mp_limb_t *xp;

      MPFR_EXP(x) = __mpfr_emin;
      xn = (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB;
      xp = MPFR_MANT(x);
      xp[xn] = GMP_LIMB_HIGHBIT;
      MPN_ZERO(xp, xn);
      inex = 1;
    }
  else
    {
      MPFR_SET_ZERO(x);
      inex = -1;
    }
  if (MPFR_SIGN(x) != sign)
    MPFR_CHANGE_SIGN(x);
  __mpfr_flags |= MPFR_FLAGS_INEXACT | MPFR_FLAGS_UNDERFLOW;
  return sign > 0 ? inex : -inex;
}

int
mpfr_integer_p (mpfr_srcptr x)
{
  mpfr_exp_t expo;
  mpfr_prec_t prec;
  mp_size_t xn;
  mp_limb_t *xp;

  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x)))
    return (MPFR_IS_ZERO(x));

  expo = MPFR_GET_EXP (x);
  if (expo <= 0)
    return 0;

  prec = MPFR_PREC(x);
  if ((mpfr_uexp_t) expo >= (mpfr_uexp_t) prec)
    return 1;

  /* 0 < expo < prec */

  xn =  (mp_size_t) ((prec - 1) / GMP_NUMB_BITS);  /* index of last limb */
  xn -= (mp_size_t) (expo / GMP_NUMB_BITS);
  /* now the index of the last limb containing bits of the fractional part */

  xp = MPFR_MANT(x);
  MPFR_ASSERTN(xn >= 0);
  if (xp[xn] << (expo % GMP_NUMB_BITS) != 0)
    return 0;
  while (--xn >= 0)
    if (xp[xn] != 0)
      return 0;
  return 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;
}

// remove trailing zeros (only by limbs)
void mpfr_remove_trailing_zeros(mpfr_t x) {
  unsigned int xn = MPFR_LIMB_SIZE(x); 
  mp_limb_t* xp = MPFR_MANT(x);
  unsigned int i = 0;
  while (i < xn && xp[i] == 0) i++;
  if (i > 0 && i < xn) mpfr_round_prec(x, GMP_RNDN, (xn-i)*BITS_PER_MP_LIMB); 
}