C++ MPN_ZERO函数代码示例

void
mpz_set_f (mpz_ptr w, mpf_srcptr u)
{
  mp_ptr    wp, up;
  mp_size_t size;
  mp_exp_t  exp;

  /* abs(u)<1 truncates to zero */
  exp = EXP (u);
  if (exp <= 0)
    {
      SIZ(w) = 0;
      return;
    }

  wp = MPZ_REALLOC (w, exp);
  up = PTR(u);

  size = SIZ (u);
  SIZ(w) = (size >= 0 ? exp : -exp);
  size = ABS (size);

  if (exp > size)
    {
      /* pad with low zeros to get a total "exp" many limbs */
      mp_size_t  zeros = exp - size;
      MPN_ZERO (wp, zeros);
      wp += zeros;
    }
  else
    {
      /* exp<=size, trucate to the high "exp" many limbs */
      up += (size - exp);
      size = exp;
    }

  MPN_COPY (wp, up, size);
}

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

void tc4_addlsh1_unsigned(mp_ptr rp, mp_size_t * rn, mp_srcptr xp, mp_size_t xn)
{
	if (xn)
	{
		if (xn >= *rn)
		{
            mp_limb_t cy;
			if (xn > *rn) MPN_ZERO(rp + *rn, xn - *rn);
#if HAVE_NATIVE_mpn_addlsh1_n
            cy = mpn_addlsh1_n(rp, rp, xp, xn);
#else
            cy = mpn_add_n(rp, rp, xp, xn);
            cy += mpn_add_n(rp, rp, xp, xn);
#endif
			if (cy) 
			{
				rp[xn] = cy;
				*rn = xn + 1;
			} else *rn = xn;
		} else
	   {
		   mp_limb_t cy;
#if HAVE_NATIVE_mpn_addlsh1_n
            cy = mpn_addlsh1_n(rp, rp, xp, xn);
#else
            cy = mpn_add_n(rp, rp, xp, xn);
            cy += mpn_add_n(rp, rp, xp, xn);
#endif
	      if (cy) cy = mpn_add_1(rp + xn, rp + xn, *rn - xn, cy);
		   if (cy) 
		   {
			   rp[*rn] = cy;
			   (*rn)++;
		   }
		}
	}
}

/* Since MPFR-3.0, return the usual inexact value.
   The erange flag is set if an error occurred in the conversion
   (y is NaN, +Inf, or -Inf that have no equivalent in mpf)
*/
int
mpfr_get_f (mpf_ptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
{
  int inex;
  mp_size_t sx, sy;
  mpfr_prec_t precx, precy;
  mp_limb_t *xp;
  int sh;

  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(y)))
    {
      if (MPFR_IS_ZERO(y))
        {
          mpf_set_ui (x, 0);
          return 0;
        }
      else if (MPFR_IS_NAN (y))
        {
          MPFR_SET_ERANGEFLAG ();
          return 0;
        }
      else /* y is plus infinity (resp. minus infinity), set x to the maximum
              value (resp. the minimum value) in precision PREC(x) */
        {
          int i;
          mp_limb_t *xp;

          MPFR_SET_ERANGEFLAG ();

          /* To this day, [mp_exp_t] and mp_size_t are #defined as the same
             type */
          EXP (x) = MP_SIZE_T_MAX;

          sx = PREC (x);
          SIZ (x) = sx;
          xp = PTR (x);
          for (i = 0; i < sx; i++)
            xp[i] = MPFR_LIMB_MAX;

          if (MPFR_IS_POS (y))
            return -1;
          else
            {
              mpf_neg (x, x);
              return +1;
            }
        }
    }

  sx = PREC(x); /* number of limbs of the mantissa of x */

  precy = MPFR_PREC(y);
  precx = (mpfr_prec_t) sx * GMP_NUMB_BITS;
  sy = MPFR_LIMB_SIZE (y);

  xp = PTR (x);

  /* since mpf numbers are represented in base 2^GMP_NUMB_BITS,
     we loose -EXP(y) % GMP_NUMB_BITS bits in the most significant limb */
  sh = MPFR_GET_EXP(y) % GMP_NUMB_BITS;
  sh = sh <= 0 ? - sh : GMP_NUMB_BITS - sh;
  MPFR_ASSERTD (sh >= 0);
  if (precy + sh <= precx) /* we can copy directly */
    {
      mp_size_t ds;

      MPFR_ASSERTN (sx >= sy);
      ds = sx - sy;

      if (sh != 0)
        {
          mp_limb_t out;
          out = mpn_rshift (xp + ds, MPFR_MANT(y), sy, sh);
          MPFR_ASSERTN (ds > 0 || out == 0);
          if (ds > 0)
            xp[--ds] = out;
        }
      else
        MPN_COPY (xp + ds, MPFR_MANT (y), sy);
      if (ds > 0)
        MPN_ZERO (xp, ds);
      EXP(x) = (MPFR_GET_EXP(y) + sh) / GMP_NUMB_BITS;
      inex = 0;
    }
  else /* we have to round to precx - sh bits */
    {
      mpfr_t z;
      mp_size_t sz;

      /* Recall that precx = (mpfr_prec_t) sx * GMP_NUMB_BITS, thus removing
         sh bits (sh < GMP_NUMB_BITSS) won't reduce the number of limbs. */
      mpfr_init2 (z, precx - sh);
      sz = MPFR_LIMB_SIZE (z);
      MPFR_ASSERTN (sx == sz);

      inex = mpfr_set (z, y, rnd_mode);
//.........这里部分代码省略.........

void
mpn_toom3_sqr_n (mp_ptr c, mp_srcptr a, mp_size_t n, mp_ptr t)
{
  mp_size_t k, k1, kk1, r, twok, twor;
  mp_limb_t cy, saved, vinf0, cinf0;
  mp_ptr trec;
  int sa;
  mp_ptr c1, c2, c3, c4;

  ASSERT(GMP_NUMB_BITS >= 6);
  ASSERT(n >= 17); /* so that r <> 0 and 5k+3 <= 2n */

  /* the algorithm is the same as mpn_mul_n_tc3, with b=a */

  k = (n + 2) / 3; /* ceil(n/3) */
  twok = 2 * k;
  k1 = k + 1;
  kk1 = k + k1;
  r = n - twok;   /* last chunk */
  twor = 2 * r;

  c1 = c + k;
  c2 = c1 + k;
  c3 = c2 + k;
  c4 = c3 + k;

  trec = t + 4 * k + 3; /* trec = v2 + (2k+2) */

  cy = mpn_add_n (c, a, a + twok, r);
  if (r < k)
    __GMPN_ADD_1 (cy, c + r, a + r, k - r, cy);
  c3[2] = (c1[0] = cy) + mpn_add_n (c2 + 2, c, a + k, k);

#define v2 (t+2*k+1)
#define vinf (t+4*k+2)

  TOOM3_SQR_REC (t, c2 + 2, k1, trec);

  sa = (c[k] != 0) ? 1 : mpn_cmp (c, a + k, k);
  c[k] = (sa >= 0) ? c[k] - mpn_sub_n (c, c, a + k, k)
    : mpn_sub_n (c, a + k, c, k);

  TOOM3_SQR_REC (c2, c, k1, trec);

#ifdef HAVE_NATIVE_mpn_addlsh1_n
  c1[0] = mpn_addlsh1_n (c, a + k, a + twok, r);
  if (r < k)
    __GMPN_ADD_1 (c1[0], c + r, a + k + r, k - r, c1[0]);
  c1[0] = 2 * c1[0] + mpn_addlsh1_n (c, a, c, k);
#else
  c[r] = mpn_lshift (c, a + twok, r, 1);
  if (r < k)
    MPN_ZERO(c + r + 1, k - r);
  c1[0] += mpn_add_n (c, c, a + k, k);
  mpn_lshift (c, c, k1, 1);
  c1[0] += mpn_add_n (c, c, a, k);
#endif

  TOOM3_SQR_REC (v2, c, k1, trec);

  TOOM3_SQR_REC (c, a, k, trec);

#ifdef HAVE_NATIVE_mpn_addlsh1_n
  mpn_addlsh1_n (v2, v2, c2, kk1);
#else
  mpn_lshift (t + 4 * k + 2, c2, kk1, 1);
  mpn_add_n (v2, v2, t + 4 * k + 2, kk1);
#endif

  saved = c4[0];
  TOOM3_SQR_REC (c4, a + twok, r, trec);
  cinf0 = mpn_add_n (vinf, c4, c, twor);
  vinf0 = c4[0];
  c4[0] = saved;

  toom3_interpolate (c, t, v2, c2, vinf, k, r, 1, vinf0, cinf0, vinf + twor);

#undef v2
#undef vinf
}

  choke me
#endif

void
mpq_set_d (mpq_ptr dest, double d)
{
  int negative;
  mp_exp_t exp;
  mp_limb_t tp[LIMBS_PER_DOUBLE];
  mp_ptr np, dp;
  mp_size_t nn, dn;
  int c;

  negative = d < 0;
  d = ABS (d);

  exp = __gmp_extract_double (tp, d);

  /* There are two main version of the conversion.  The `then' arm handles
     things that have a fractional part, while the `else' part handles
     only integers.  */
#if BITS_PER_MP_LIMB == 32
  if (exp <= 1 || (exp == 2 && tp[0] != 0))
#else
  if (exp <= 1)
#endif
    {
      if (d == 0.0)
	{
	  SIZ(&(dest->_mp_num)) = 0;
	  SIZ(&(dest->_mp_den)) = 1;
	  PTR(&(dest->_mp_den))[0] = 1;
	  return;
	}

      dn = -exp;
      if (dest->_mp_num._mp_alloc < 3)
	_mpz_realloc (&(dest->_mp_num), 3);
      np = PTR(&(dest->_mp_num));
#if BITS_PER_MP_LIMB == 32
      if ((tp[0] | tp[1]) == 0)
	np[0] = tp[2], nn = 1;
      else if (tp[0] == 0)
	np[1] = tp[2], np[0] = tp[1], nn = 2;
      else
	np[2] = tp[2], np[1] = tp[1], np[0] = tp[0], nn = 3;
#else
      if (tp[0] == 0)
	np[0] = tp[1], nn = 1;
      else
	np[1] = tp[1], np[0] = tp[0], nn = 2;
#endif
      dn += nn + 1;
      if (dest->_mp_den._mp_alloc < dn)
	_mpz_realloc (&(dest->_mp_den), dn);
      dp = PTR(&(dest->_mp_den));
      MPN_ZERO (dp, dn - 1);
      dp[dn - 1] = 1;
      count_trailing_zeros (c, np[0] | dp[0]);
      if (c != 0)
	{
	  mpn_rshift (np, np, nn, c);
	  nn -= np[nn - 1] == 0;
	  mpn_rshift (dp, dp, dn, c);
	  dn -= dp[dn - 1] == 0;
	}
      SIZ(&(dest->_mp_den)) = dn;
      SIZ(&(dest->_mp_num)) = negative ? -nn : nn;
    }
  else
    {
      nn = exp;
      if (dest->_mp_num._mp_alloc < nn)
	_mpz_realloc (&(dest->_mp_num), nn);
      np = PTR(&(dest->_mp_num));
#if BITS_PER_MP_LIMB == 32
      switch (nn)
        {
	default:
          MPN_ZERO (np, nn - 3);
          np += nn - 3;
	  /* fall through */
	case 3:
	  np[2] = tp[2], np[1] = tp[1], np[0] = tp[0];
	  break;
	case 2:
	  np[1] = tp[2], np[0] = tp[1];
	  break;
	}
#else
      switch (nn)
        {
	default:
	  MPN_ZERO (np, nn - 2);
	  np += nn - 2;
	  /* fall through */
	case 2:
	  np[1] = tp[1], np[0] = tp[0];
	  break;
	}
//.........这里部分代码省略.........

/* FIXME:
    x Take scratch parameter, and figure out scratch need.

    x Use some fallback for small M->n?
*/
static mp_size_t
hgcd_matrix_apply (const struct hgcd_matrix *M,
		   mp_ptr ap, mp_ptr bp,
		   mp_size_t n)
{
  mp_size_t an, bn, un, vn, nn;
  mp_size_t mn[2][2];
  mp_size_t modn;
  mp_ptr tp, sp, scratch;
  mp_limb_t cy;
  unsigned i, j;

  TMP_DECL;

  ASSERT ( (ap[n-1] | bp[n-1]) > 0);

  an = n;
  MPN_NORMALIZE (ap, an);
  bn = n;
  MPN_NORMALIZE (bp, bn);

  for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
      {
	mp_size_t k;
	k = M->n;
	MPN_NORMALIZE (M->p[i][j], k);
	mn[i][j] = k;
      }

  ASSERT (mn[0][0] > 0);
  ASSERT (mn[1][1] > 0);
  ASSERT ( (mn[0][1] | mn[1][0]) > 0);

  TMP_MARK;

  if (mn[0][1] == 0)
    {
      /* A unchanged, M = (1, 0; q, 1) */
      ASSERT (mn[0][0] == 1);
      ASSERT (M->p[0][0][0] == 1);
      ASSERT (mn[1][1] == 1);
      ASSERT (M->p[1][1][0] == 1);

      /* Put B <-- B - q A */
      nn = submul (bp, bn, ap, an, M->p[1][0], mn[1][0]);
    }
  else if (mn[1][0] == 0)
    {
      /* B unchanged, M = (1, q; 0, 1) */
      ASSERT (mn[0][0] == 1);
      ASSERT (M->p[0][0][0] == 1);
      ASSERT (mn[1][1] == 1);
      ASSERT (M->p[1][1][0] == 1);

      /* Put A  <-- A - q * B */
      nn = submul (ap, an, bp, bn, M->p[0][1], mn[0][1]);
    }
  else
    {
      /* A = m00 a + m01 b  ==> a <= A / m00, b <= A / m01.
	 B = m10 a + m11 b  ==> a <= B / m10, b <= B / m11. */
      un = MIN (an - mn[0][0], bn - mn[1][0]) + 1;
      vn = MIN (an - mn[0][1], bn - mn[1][1]) + 1;

      nn = MAX (un, vn);
      /* In the range of interest, mulmod_bnm1 should always beat mullo. */
      modn = mpn_mulmod_bnm1_next_size (nn + 1);

      scratch = TMP_ALLOC_LIMBS (mpn_mulmod_bnm1_itch (modn, modn, M->n));
      tp = TMP_ALLOC_LIMBS (modn);
      sp = TMP_ALLOC_LIMBS (modn);

      ASSERT (n <= 2*modn);

      if (n > modn)
	{
	  cy = mpn_add (ap, ap, modn, ap + modn, n - modn);
	  MPN_INCR_U (ap, modn, cy);

	  cy = mpn_add (bp, bp, modn, bp + modn, n - modn);
	  MPN_INCR_U (bp, modn, cy);

	  n = modn;
	}

      mpn_mulmod_bnm1 (tp, modn, ap, n, M->p[1][1], mn[1][1], scratch);
      mpn_mulmod_bnm1 (sp, modn, bp, n, M->p[0][1], mn[0][1], scratch);

      /* FIXME: Handle the small n case in some better way. */
      if (n + mn[1][1] < modn)
	MPN_ZERO (tp + n + mn[1][1], modn - n - mn[1][1]);
      if (n + mn[0][1] < modn)
	MPN_ZERO (sp + n + mn[0][1], modn - n - mn[0][1]);

//.........这里部分代码省略.........

int
mpfr_sqrt (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
{
  mp_size_t rsize; /* number of limbs of r (plus 1 if exact limb multiple) */
  mp_size_t rrsize;
  mp_size_t usize; /* number of limbs of u */
  mp_size_t tsize; /* number of limbs of the sqrtrem remainder */
  mp_size_t k;
  mp_size_t l;
  mpfr_limb_ptr rp, rp0;
  mpfr_limb_ptr up;
  mpfr_limb_ptr sp;
  mp_limb_t sticky0; /* truncated part of input */
  mp_limb_t sticky1; /* truncated part of rp[0] */
  mp_limb_t sticky;
  int odd_exp;
  int sh; /* number of extra bits in rp[0] */
  int inexact; /* return ternary flag */
  mpfr_exp_t expr;
  MPFR_TMP_DECL(marker);

  MPFR_LOG_FUNC
    (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (u), mpfr_log_prec, u, rnd_mode),
     ("y[%Pu]=%.*Rg inexact=%d",
      mpfr_get_prec (r), mpfr_log_prec, r, inexact));

  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(u)))
    {
      if (MPFR_IS_NAN(u))
        {
          MPFR_SET_NAN(r);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_ZERO(u))
        {
          /* 0+ or 0- */
          MPFR_SET_SAME_SIGN(r, u);
          MPFR_SET_ZERO(r);
          MPFR_RET(0); /* zero is exact */
        }
      else
        {
          MPFR_ASSERTD(MPFR_IS_INF(u));
          /* sqrt(-Inf) = NAN */
          if (MPFR_IS_NEG(u))
            {
              MPFR_SET_NAN(r);
              MPFR_RET_NAN;
            }
          MPFR_SET_POS(r);
          MPFR_SET_INF(r);
          MPFR_RET(0);
        }
    }
  if (MPFR_UNLIKELY(MPFR_IS_NEG(u)))
    {
      MPFR_SET_NAN(r);
      MPFR_RET_NAN;
    }
  MPFR_SET_POS(r);

  MPFR_TMP_MARK (marker);
  MPFR_UNSIGNED_MINUS_MODULO(sh,MPFR_PREC(r));
  if (sh == 0 && rnd_mode == MPFR_RNDN)
    sh = GMP_NUMB_BITS; /* ugly case */
  rsize = MPFR_LIMB_SIZE(r) + (sh == GMP_NUMB_BITS);
  /* rsize is the number of limbs of r + 1 if exact limb multiple and rounding
     to nearest, this is the number of wanted limbs for the square root */
  rrsize = rsize + rsize;
  usize = MPFR_LIMB_SIZE(u); /* number of limbs of u */
  rp0 = MPFR_MANT(r);
  rp = (sh < GMP_NUMB_BITS) ? rp0 : MPFR_TMP_LIMBS_ALLOC (rsize);
  up = MPFR_MANT(u);
  sticky0 = MPFR_LIMB_ZERO; /* truncated part of input */
  sticky1 = MPFR_LIMB_ZERO; /* truncated part of rp[0] */
  odd_exp = (unsigned int) MPFR_GET_EXP (u) & 1;
  inexact = -1; /* return ternary flag */

  sp = MPFR_TMP_LIMBS_ALLOC (rrsize);

  /* copy the most significant limbs of u to {sp, rrsize} */
  if (MPFR_LIKELY(usize <= rrsize)) /* in case r and u have the same precision,
                                       we have indeed rrsize = 2 * usize */
    {
      k = rrsize - usize;
      if (MPFR_LIKELY(k))
        MPN_ZERO (sp, k);
      if (odd_exp)
        {
          if (MPFR_LIKELY(k))
            sp[k - 1] = mpn_rshift (sp + k, up, usize, 1);
          else
            sticky0 = mpn_rshift (sp, up, usize, 1);
        }
      else
        MPN_COPY (sp + rrsize - usize, up, usize);
    }
  else /* usize > rrsize: truncate the input */
    {
      k = usize - rrsize;
//.........这里部分代码省略.........

/* rp[n-1..0] = bp[bn-1..0] ^ ep[en-1..0] mod mp[n-1..0]
   Requires that mp[n-1..0] is odd.  FIXME: is this true?
   Requires that ep[en-1..0] is > 1.
   Uses scratch space at tp of 3n+1 limbs.  */
void
mpn_powm_sec (mp_ptr rp, mp_srcptr bp, mp_size_t bn,
	      mp_srcptr ep, mp_size_t en,
	      mp_srcptr mp, mp_size_t n, mp_ptr tp)
{
  mp_limb_t minv;
  int cnt;
  mp_bitcnt_t ebi;
  int windowsize, this_windowsize;
  mp_limb_t expbits;
  mp_ptr pp, this_pp;
  long i;
  int cnd;

  ASSERT (en > 1 || (en == 1 && ep[0] > 0));
  ASSERT (n >= 1 && ((mp[0] & 1) != 0));

  count_leading_zeros (cnt, ep[en - 1]);
  ebi = (mp_bitcnt_t) en * GMP_LIMB_BITS - cnt;

  windowsize = win_size (ebi);

  binvert_limb (minv, mp[0]);
  minv = -minv;

  pp = tp + 4 * n;

  this_pp = pp;
  this_pp[n] = 1;
  redcify (this_pp, this_pp + n, 1, mp, n, tp + 6 * n);
  this_pp += n;
  redcify (this_pp, bp, bn, mp, n, tp + 6 * n);

  /* Precompute powers of b and put them in the temporary area at pp.  */
  for (i = (1 << windowsize) - 2; i > 0; i--)
    {
      mpn_mul_basecase (tp, this_pp, n, pp + n, n);
      this_pp += n;
      mpn_redc_1_sec (this_pp, tp, mp, n, minv);
    }

  expbits = getbits (ep, ebi, windowsize);
  if (ebi < windowsize)
    ebi = 0;
  else
    ebi -= windowsize;

#if WANT_CACHE_SECURITY
  mpn_tabselect (rp, pp, n, 1 << windowsize, expbits);
#else
  MPN_COPY (rp, pp + n * expbits, n);
#endif

  while (ebi != 0)
    {
      expbits = getbits (ep, ebi, windowsize);
      this_windowsize = windowsize;
      if (ebi < windowsize)
	{
	  this_windowsize -= windowsize - ebi;
	  ebi = 0;
	}
      else
	ebi -= windowsize;

      do
	{
	  mpn_local_sqr (tp, rp, n, tp + 2 * n);
	  mpn_redc_1_sec (rp, tp, mp, n, minv);
	  this_windowsize--;
	}
      while (this_windowsize != 0);

#if WANT_CACHE_SECURITY
      mpn_tabselect (tp + 2*n, pp, n, 1 << windowsize, expbits);
      mpn_mul_basecase (tp, rp, n, tp + 2*n, n);
#else
      mpn_mul_basecase (tp, rp, n, pp + n * expbits, n);
#endif
      mpn_redc_1_sec (rp, tp, mp, n, minv);
    }

  MPN_COPY (tp, rp, n);
  MPN_ZERO (tp + n, n);
  mpn_redc_1_sec (rp, tp, mp, n, minv);
  cnd = mpn_sub_n (tp, rp, mp, n);	/* we need just retval */
  mpn_subcnd_n (rp, rp, mp, n, !cnd);
}

void
mpq_set_f (mpq_ptr q, mpf_srcptr f)
{
  mp_size_t  fexp = EXP(f);
  mp_ptr     fptr = PTR(f);
  mp_size_t  fsize = SIZ(f);
  mp_size_t  abs_fsize = ABS(fsize);
  mp_limb_t  flow;

  if (fsize == 0)
    {
      /* set q=0 */
      SIZ(NUM(q)) = 0;
      SIZ(DEN(q)) = 1;
      PTR(DEN(q))[0] = 1;
      return;
    }

  /* strip low zero limbs from f */
  flow = *fptr;
  MPN_STRIP_LOW_ZEROS_NOT_ZERO (fptr, abs_fsize, flow);

  if (fexp >= abs_fsize)
    {
      /* radix point is to the right of the limbs, no denominator */
      mp_ptr  num_ptr;

      num_ptr = MPZ_NEWALLOC (mpq_numref (q), fexp);
      MPN_ZERO (num_ptr, fexp - abs_fsize);
      MPN_COPY (num_ptr + fexp - abs_fsize, fptr, abs_fsize);

      SIZ(NUM(q)) = fsize >= 0 ? fexp : -fexp;
      SIZ(DEN(q)) = 1;
      PTR(DEN(q))[0] = 1;
    }
  else
    {
      /* radix point is within or to the left of the limbs, use denominator */
      mp_ptr     num_ptr, den_ptr;
      mp_size_t  den_size;

      den_size = abs_fsize - fexp;
      num_ptr = MPZ_NEWALLOC (mpq_numref (q), abs_fsize);
      den_ptr = MPZ_NEWALLOC (mpq_denref (q), den_size+1);

      if (flow & 1)
        {
          /* no powers of two to strip from numerator */

          MPN_COPY (num_ptr, fptr, abs_fsize);
          MPN_ZERO (den_ptr, den_size);
          den_ptr[den_size] = 1;
        }
      else
        {
          /* right shift numerator, adjust denominator accordingly */
          int  shift;

          den_size--;
          count_trailing_zeros (shift, flow);

          mpn_rshift (num_ptr, fptr, abs_fsize, shift);
          abs_fsize -= (num_ptr[abs_fsize-1] == 0);

          MPN_ZERO (den_ptr, den_size);
          den_ptr[den_size] = GMP_LIMB_HIGHBIT >> (shift-1);
        }

      SIZ(NUM(q)) = fsize >= 0 ? abs_fsize : -abs_fsize;
      SIZ(DEN(q)) = den_size + 1;
    }
}

//.........这里部分代码省略.........
	      mp_srcptr mp, mp_size_t n, mp_ptr tp)
{
  mp_limb_t mip[2];
  int cnt;
  long ebi;
  int windowsize, this_windowsize;
  mp_limb_t expbits;
  mp_ptr pp, this_pp, last_pp;
  long i;
  int redc_x;
  TMP_DECL;

  ASSERT (en > 1 || (en == 1 && ep[0] > 1));
  ASSERT (n >= 1 && ((mp[0] & 1) != 0));

  TMP_MARK;

  count_leading_zeros (cnt, ep[en - 1]);
  ebi = en * GMP_LIMB_BITS - cnt;

  windowsize = win_size (ebi);

  if (BELOW_THRESHOLD (n, REDC_2_THRESHOLD))
    {
      binvert_limb (mip[0], mp[0]);
      mip[0] = -mip[0];
      redc_x = 1;
    }
#if defined (HAVE_NATIVE_mpn_addmul_2)
  else
    {
      mpn_binvert (mip, mp, 2, tp);
      mip[0] = -mip[0]; mip[1] = ~mip[1];
      redc_x = 2;
    }
#endif
#if 0
  mpn_binvert (mip, mp, n, tp);
  redc_x = 0;
#endif

  pp = TMP_ALLOC_LIMBS (n << windowsize);

  this_pp = pp;
  this_pp[n] = 1;
  redcify (this_pp, this_pp + n, 1, mp, n);
  this_pp += n;
  redcify (this_pp, bp, bn, mp, n);

  /* Precompute powers of b and put them in the temporary area at pp.  */
  for (i = (1 << windowsize) - 2; i > 0; i--)
    {
      last_pp = this_pp;
      this_pp += n;
      mpn_mul_n (tp, last_pp, pp + n, n);
      MPN_REDC_X (this_pp, tp, mp, n, mip);
    }

  expbits = getbits (ep, ebi, windowsize);
  ebi -= windowsize;
  if (ebi < 0)
    ebi = 0;

  MPN_COPY (rp, pp + n * expbits, n);

  while (ebi != 0)
    {
      expbits = getbits (ep, ebi, windowsize);
      ebi -= windowsize;
      this_windowsize = windowsize;
      if (ebi < 0)
	{
	  this_windowsize += ebi;
	  ebi = 0;
	}

      do
	{
	  mpn_sqr_n (tp, rp, n);
	  MPN_REDC_X (rp, tp, mp, n, mip);
	  this_windowsize--;
	}
      while (this_windowsize != 0);

#if WANT_CACHE_SECURITY
      mpn_tabselect (tp + 2*n, pp, n, 1 << windowsize, expbits);
      mpn_mul_n (tp, rp, tp + 2*n, n);
#else
      mpn_mul_n (tp, rp, pp + n * expbits, n);
#endif
      MPN_REDC_X (rp, tp, mp, n, mip);
    }

  MPN_COPY (tp, rp, n);
  MPN_ZERO (tp + n, n);
  MPN_REDC_X (rp, tp, mp, n, mip);
  if (mpn_cmp (rp, mp, n) >= 0)
    mpn_sub_n (rp, rp, mp, n);
  TMP_FREE;
}

void
mpn_toom22_mul (mp_ptr pp,
                mp_srcptr ap, mp_size_t an,
                mp_srcptr bp, mp_size_t bn,
                mp_ptr scratch)
{
    const int __gmpn_cpuvec_initialized = 1;
    mp_size_t n, s, t;
    int vm1_neg;
    mp_limb_t cy, cy2;
    mp_ptr asm1;
    mp_ptr bsm1;

#define a0  ap
#define a1  (ap + n)
#define b0  bp
#define b1  (bp + n)

    s = an >> 1;
    n = an - s;
    t = bn - n;

    ASSERT (an >= bn);

    ASSERT (0 < s && s <= n && s >= n - 1);
    ASSERT (0 < t && t <= s);

    asm1 = pp;
    bsm1 = pp + n;

    vm1_neg = 0;

    /* Compute asm1.  */
    if (s == n)
    {
        if (mpn_cmp (a0, a1, n) < 0)
        {
            mpn_sub_n (asm1, a1, a0, n);
            vm1_neg = 1;
        }
        else
        {
            mpn_sub_n (asm1, a0, a1, n);
        }
    }
    else /* n - s == 1 */
    {
        if (a0[s] == 0 && mpn_cmp (a0, a1, s) < 0)
        {
            mpn_sub_n (asm1, a1, a0, s);
            asm1[s] = 0;
            vm1_neg = 1;
        }
        else
        {
            asm1[s] = a0[s] - mpn_sub_n (asm1, a0, a1, s);
        }
    }

    /* Compute bsm1.  */
    if (t == n)
    {
        if (mpn_cmp (b0, b1, n) < 0)
        {
            mpn_sub_n (bsm1, b1, b0, n);
            vm1_neg ^= 1;
        }
        else
        {
            mpn_sub_n (bsm1, b0, b1, n);
        }
    }
    else
    {
        if (mpn_zero_p (b0 + t, n - t) && mpn_cmp (b0, b1, t) < 0)
        {
            mpn_sub_n (bsm1, b1, b0, t);
            MPN_ZERO (bsm1 + t, n - t);
            vm1_neg ^= 1;
        }
        else
        {
            mpn_sub (bsm1, b0, n, b1, t);
        }
    }

#define v0	pp				/* 2n */
#define vinf	(pp + 2 * n)			/* s+t */
#define vm1	scratch				/* 2n */
#define scratch_out	scratch + 2 * n

    /* vm1, 2n limbs */
    TOOM22_MUL_N_REC (vm1, asm1, bsm1, n, scratch_out);

    if (s > t)  TOOM22_MUL_REC (vinf, a1, s, b1, t, scratch_out);
    else        TOOM22_MUL_N_REC (vinf, a1, b1, s, scratch_out);

    /* v0, 2n limbs */
    TOOM22_MUL_N_REC (v0, ap, bp, n, scratch_out);

//.........这里部分代码省略.........

//.........这里部分代码省略.........
    if (usize == 0 || (vsize == 1 && vp[0] == 1))
    {
        gst_mpz_copy_abs (g, v);
        return;
    }

    /* GCD(U, 0) == GCD (1, V) == U.  */
    if (vsize == 0 || (usize == 1 && up[0] == 1))
    {
        gst_mpz_copy_abs (g, u);
        return;
    }

    if (usize == 1)
    {
        gst_mpz_realloc (g, 1);
        g->size = 1;
        g->d[0] = mpn_gcd_1 (vp, vsize, up[0]);
        return;
    }

    if (vsize == 1)
    {
        gst_mpz_realloc (g, 1);
        g->size = 1;
        g->d[0] = mpn_gcd_1 (up, usize, vp[0]);
        return;
    }

    /*  Eliminate low zero bits from U and V and move to temporary storage.  */
    u_zero_bits = mpn_scan1 (up, 0);
    u_zero_limbs = u_zero_bits / BITS_PER_MP_LIMB;
    u_zero_bits &= BITS_PER_MP_LIMB - 1;
    up += u_zero_limbs;
    usize -= u_zero_limbs;

    /* Operands could be destroyed for big-endian case, but let's be tidy.  */
    tp = up;
    up = (mp_ptr) alloca (usize * SIZEOF_MP_LIMB_T);
    if (u_zero_bits != 0)
    {
        mpn_rshift (up, tp, usize, u_zero_bits);
        usize -= up[usize - 1] == 0;
    }
    else
        MPN_COPY (up, tp, usize);

    v_zero_bits = mpn_scan1 (vp, 0);
    v_zero_limbs = v_zero_bits / BITS_PER_MP_LIMB;
    v_zero_bits &= BITS_PER_MP_LIMB - 1;
    vp += v_zero_limbs;
    vsize -= v_zero_limbs;

    /* Operands could be destroyed for big-endian case, but let's be tidy.  */
    tp = vp;
    vp = (mp_ptr) alloca (vsize * SIZEOF_MP_LIMB_T);
    if (v_zero_bits != 0)
    {
        mpn_rshift (vp, tp, vsize, v_zero_bits);
        vsize -= vp[vsize - 1] == 0;
    }
    else
        MPN_COPY (vp, tp, vsize);

    if (u_zero_limbs > v_zero_limbs)
    {
        g_zero_limbs = v_zero_limbs;
        g_zero_bits = v_zero_bits;
    }
    else if (u_zero_limbs < v_zero_limbs)
    {
        g_zero_limbs = u_zero_limbs;
        g_zero_bits = u_zero_bits;
    }
    else  /*  Equal.  */
    {
        g_zero_limbs = u_zero_limbs;
        g_zero_bits = MIN (u_zero_bits, v_zero_bits);
    }

    /*  Call mpn_gcd.  The 2nd argument must not have more bits than the 1st.  */
    vsize = (usize < vsize || (usize == vsize && up[usize-1] < vp[vsize-1]))
            ? mpn_gcd (vp, vp, vsize, up, usize)
            : mpn_gcd (vp, up, usize, vp, vsize);

    /*  Here G <-- V << (g_zero_limbs*BITS_PER_MP_LIMB + g_zero_bits).  */
    gsize = vsize + g_zero_limbs;
    if (g_zero_bits != 0)
    {
        mp_limb_t cy_limb;
        gsize += (vp[vsize - 1] >> (BITS_PER_MP_LIMB - g_zero_bits)) != 0;
        if (g->alloc < gsize)
            gst_mpz_realloc (g, gsize);
        MPN_ZERO (g->d, g_zero_limbs);

        tp = g->d + g_zero_limbs;
        cy_limb = mpn_lshift (tp, vp, vsize, g_zero_bits);
        if (cy_limb != 0)
            tp[vsize] = cy_limb;
    }

static unsigned long int
lc (mp_ptr rp, gmp_randstate_t rstate)
{
  mp_ptr tp, seedp, ap;
  mp_size_t ta;
  mp_size_t tn, seedn, an;
  unsigned long int m2exp;
  unsigned long int bits;
  int cy;
  mp_size_t xn;
  gmp_rand_lc_struct *p;
  TMP_DECL;

  p = (gmp_rand_lc_struct *) RNG_STATE (rstate);

  m2exp = p->_mp_m2exp;

  seedp = PTR (p->_mp_seed);
  seedn = SIZ (p->_mp_seed);

  ap = PTR (p->_mp_a);
  an = SIZ (p->_mp_a);

  /* Allocate temporary storage.  Let there be room for calculation of
     (A * seed + C) % M, or M if bigger than that.  */

  TMP_MARK;

  ta = an + seedn + 1;
  tn = BITS_TO_LIMBS (m2exp);
  if (ta <= tn) /* that is, if (ta < tn + 1) */
    {
      mp_size_t tmp = an + seedn;
      ta = tn + 1;
      tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
      MPN_ZERO (&tp[tmp], ta - tmp); /* mpn_mul won't zero it out.  */
    }
  else
    tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);

  /* t = a * seed.  NOTE: an is always > 0; see initialization.  */
  ASSERT (seedn >= an && an > 0);
  mpn_mul (tp, seedp, seedn, ap, an);

  /* t = t + c.  NOTE: tn is always >= p->_cn (precondition for __GMPN_ADD);
     see initialization.  */
  ASSERT (tn >= p->_cn);
  __GMPN_ADD (cy, tp, tp, tn, p->_cp, p->_cn);

  /* t = t % m */
  tp[m2exp / GMP_NUMB_BITS] &= (CNST_LIMB (1) << m2exp % GMP_NUMB_BITS) - 1;

  /* Save result as next seed.  */
  MPN_COPY (PTR (p->_mp_seed), tp, tn);

  /* Discard the lower m2exp/2 of the result.  */
  bits = m2exp / 2;
  xn = bits / GMP_NUMB_BITS;

  tn -= xn;
  if (tn > 0)
    {
      unsigned int cnt = bits % GMP_NUMB_BITS;
      if (cnt != 0)
	{
	  mpn_rshift (tp, tp + xn, tn, cnt);
	  MPN_COPY_INCR (rp, tp, xn + 1);
	}
      else			/* Even limb boundary.  */
	MPN_COPY_INCR (rp, tp + xn, tn);
    }

  TMP_FREE;

  /* Return number of valid bits in the result.  */
  return (m2exp + 1) / 2;
}

void
mpz_combit (mpz_ptr d, mp_bitcnt_t bit_index)
{
  mp_size_t dsize = SIZ(d);
  mp_ptr dp = PTR(d);

  mp_size_t limb_index = bit_index / GMP_NUMB_BITS;
  mp_limb_t bit = (CNST_LIMB (1) << (bit_index % GMP_NUMB_BITS));

  /* Check for the most common case: Positive input, no realloc or
     normalization needed. */
  if (limb_index + 1 < dsize)
    dp[limb_index] ^= bit;

  /* Check for the hairy case. d < 0, and we have all zero bits to the
     right of the bit to toggle. */
  else if (limb_index < -dsize
	   && (limb_index == 0 || mpn_zero_p (dp, limb_index))
	   && (dp[limb_index] & (bit - 1)) == 0)
    {
      ASSERT (dsize < 0);
      dsize = -dsize;

      if (dp[limb_index] & bit)
	{
	  /* We toggle the least significant one bit. Corresponds to
	     an add, with potential carry propagation, on the absolute
	     value. */
	  dp = MPZ_REALLOC (d, 1 + dsize);
	  dp[dsize] = 0;
	  MPN_INCR_U (dp + limb_index, 1 + dsize - limb_index, bit);
	  SIZ(d) = - dsize - dp[dsize];
	}
      else
	{
	  /* We toggle a zero bit, subtract from the absolute value. */
	  MPN_DECR_U (dp + limb_index, dsize - limb_index, bit);
	  /* The absolute value shrinked by at most one bit. */
	  dsize -= dp[dsize - 1] == 0;
	  ASSERT (dsize > 0 && dp[dsize - 1] != 0);
	  SIZ (d) = -dsize;
	}
    }
  else
    {
      /* Simple case: Toggle the bit in the absolute value. */
      dsize = ABS(dsize);
      if (limb_index < dsize)
	{
	  mp_limb_t	 dlimb;
	  dlimb = dp[limb_index] ^ bit;
	  dp[limb_index] = dlimb;

	  /* Can happen only when limb_index = dsize - 1. Avoid SIZ(d)
	     bookkeeping in the common case. */
	  if (UNLIKELY ((dlimb == 0) + limb_index == dsize)) /* dsize == limb_index + 1 */
	    {
	      /* high limb became zero, must normalize */
	      MPN_NORMALIZE (dp, limb_index);
	      SIZ (d) = SIZ (d) >= 0 ? limb_index : -limb_index;
	    }
	}
      else
	{
	  dp = MPZ_REALLOC (d, limb_index + 1);
	  MPN_ZERO(dp + dsize, limb_index - dsize);
	  dp[limb_index++] = bit;
	  SIZ(d) = SIZ(d) >= 0 ? limb_index : -limb_index;
	}
    }
}