自学教程：C++ CNST_LIMB函数代码示例

voidcheck_limb (void){  int        i;  mp_limb_t  limb;  mpz_t      z;  char       *s;  check_one ("0", "%Md", CNST_LIMB(0));  check_one ("1", "%Md", CNST_LIMB(1));  /* "i" many 1 bits, tested against mpz_get_str in decimal and hex */  limb = 1;  mpz_init_set_ui (z, 1L);  for (i = 1; i <= GMP_LIMB_BITS; i++)    {      s = mpz_get_str (NULL, 10, z);      check_one (s, "%Mu", limb);      (*__gmp_free_func) (s, strlen (s) + 1);      s = mpz_get_str (NULL, 16, z);      check_one (s, "%Mx", limb);      (*__gmp_free_func) (s, strlen (s) + 1);      s = mpz_get_str (NULL, -16, z);      check_one (s, "%MX", limb);      (*__gmp_free_func) (s, strlen (s) + 1);      limb = 2*limb + 1;      mpz_mul_2exp (z, z, 1L);      mpz_add_ui (z, z, 1L);    }  mpz_clear (z);}

void tc4_rshift_inplace(mp_ptr rp, mp_size_t * rn, mp_size_t bits){   if (*rn)	{		if ((*rn) > 0) 		{			mpn_rshift(rp, rp, *rn, bits);	      if (rp[(*rn) - 1] == CNST_LIMB(0)) (*rn)--;		} else 		{			mpn_rshift(rp, rp, -(*rn), bits);	      if (rp[-(*rn) - 1] == CNST_LIMB(0)) (*rn)++;		}	}}

/* ARM gcc 2.95.4 was seen generating bad code for ulong->double   conversions, resulting in for instance 0x81c25113 incorrectly converted.   This test exercises that value, to see mpn_get_d has avoided the   problem.  */voidcheck_0x81c25113 (void){#if GMP_NUMB_BITS >= 32  double     want = 2176995603.0;  double     got;  mp_limb_t  np[4];  mp_size_t  nsize;  long       exp;  if (tests_dbl_mant_bits() < 32)    return;  for (nsize = 1; nsize <= numberof (np); nsize++)    {      refmpn_zero (np, nsize-1);      np[nsize-1] = CNST_LIMB(0x81c25113);      exp = - (nsize-1) * GMP_NUMB_BITS;      got = mpn_get_d (np, nsize, (mp_size_t) 0, exp);      if (got != want)        {          printf  ("mpn_get_d wrong on 2176995603 (0x81c25113)/n");          printf  ("  nsize  %ld/n", (long) nsize);          printf  ("  exp    %ld/n", exp);          d_trace ("  got    ", got);          d_trace ("  want   ", want);          abort ();        }    }#endif}

intmpz_divisible_2exp_p (mpz_srcptr a, unsigned long d){  unsigned long  i, dlimbs, dbits;  mp_ptr         ap;  mp_limb_t      dmask;  mp_size_t      asize;  asize = ABSIZ(a);  dlimbs = d / GMP_NUMB_BITS;  /* if d covers the whole of a, then only a==0 is divisible */  if (asize <= dlimbs)    return asize == 0;  /* whole limbs must be zero */  ap = PTR(a);  for (i = 0; i < dlimbs; i++)    if (ap[i] != 0)      return 0;  /* left over bits must be zero */  dbits = d % GMP_NUMB_BITS;  dmask = (CNST_LIMB(1) << dbits) - 1;  return (ap[dlimbs] & dmask) == 0;}

REGPARM_ATTR (1) static voidcfdiv_q_2exp (mpz_ptr w, mpz_srcptr u, mp_bitcnt_t cnt, int dir){  mp_size_t  wsize, usize, abs_usize, limb_cnt, i;  mp_srcptr  up;  mp_ptr     wp;  mp_limb_t  round, rmask;  usize = SIZ (u);  abs_usize = ABS (usize);  limb_cnt = cnt / GMP_NUMB_BITS;  wsize = abs_usize - limb_cnt;  if (wsize <= 0)    {      /* u < 2**cnt, so result 1, 0 or -1 according to rounding */      PTR(w)[0] = 1;      SIZ(w) = (usize == 0 || (usize ^ dir) < 0 ? 0 : dir);      return;    }  /* +1 limb to allow for mpn_add_1 below */  MPZ_REALLOC (w, wsize+1);  /* Check for rounding if direction matches u sign.     Set round if we're skipping non-zero limbs.  */  up = PTR(u);  round = 0;  rmask = ((usize ^ dir) >= 0 ? MP_LIMB_T_MAX : 0);  if (rmask != 0)    for (i = 0; i < limb_cnt && round == 0; i++)      round = up[i];  wp = PTR(w);  cnt %= GMP_NUMB_BITS;  if (cnt != 0)    {      round |= rmask & mpn_rshift (wp, up + limb_cnt, wsize, cnt);      wsize -= (wp[wsize - 1] == 0);    }  else    MPN_COPY_INCR (wp, up + limb_cnt, wsize);  if (round != 0)    {      if (wsize != 0)	{	  mp_limb_t cy;	  cy = mpn_add_1 (wp, wp, wsize, CNST_LIMB(1));	  wp[wsize] = cy;	  wsize += cy;	}      else	{	  /* We shifted something to zero.  */	  wp[0] = 1;	  wsize = 1;	}    }  SIZ(w) = (usize >= 0 ? wsize : -wsize);}

void tc4_divexact_by3(mp_ptr rp, mp_size_t * rn, mp_ptr x, mp_size_t xn){	if (xn)	{		mp_size_t xu = ABS(xn);		mpn_divexact_by3(rp, x, xu);		if (xn > 0)		{			if (rp[xu - 1] == CNST_LIMB(0)) *rn = xn - 1;			else *rn = xn;		} else		{			if (rp[xu - 1] == CNST_LIMB(0)) *rn = xn + 1;			else *rn = xn;		}		} else *rn = 0;}

void tc4_divexact_by15(mp_ptr rp, mp_size_t * rn, mp_ptr x, mp_size_t xn){	if (xn)	{		mp_size_t xu = ABS(xn);		mpn_divexact_byfobm1(rp, x, xu, CNST_LIMB(15), CNST_LIMB((~0)/15)); /* works for 32 and 64 bits */		if (xn > 0)		{			if (rp[xu - 1] == CNST_LIMB(0)) *rn = xn - 1;			else *rn = xn;		} else		{			if (rp[xu - 1] == CNST_LIMB(0)) *rn = xn + 1;			else *rn = xn;		}		} else *rn = 0;}

static mp_limb_tDO_mpn_addlsh_n(mp_ptr dst, mp_srcptr src, mp_size_t n, unsigned int s, mp_ptr ws){#if USE_MUL_1 && 0  return mpn_addmul_1(dst,src,n,CNST_LIMB(1) <<(s));#else  mp_limb_t __cy;  __cy = mpn_lshift(ws,src,n,s);  return    __cy + mpn_add_n(dst,dst,ws,n);#endif}

static voidgmp_rrandomb (mp_ptr rp, gmp_randstate_t rstate, mpir_ui nbits){  mpir_ui bi;  mp_limb_t ranm;		/* buffer for random bits */  unsigned cap_chunksize, chunksize;  mp_size_t i;  /* Set entire result to 111..1  */  i = (nbits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS - 1;  rp[i] = GMP_NUMB_MAX >> (GMP_NUMB_BITS - (nbits % GMP_NUMB_BITS)) % GMP_NUMB_BITS;  for (i = i - 1; i >= 0; i--)    rp[i] = GMP_NUMB_MAX;  _gmp_rand (&ranm, rstate, BITS_PER_RANDCALL);  cap_chunksize = nbits / (ranm % 4 + 1);  cap_chunksize += cap_chunksize == 0; /* make it at least 1 */  bi = nbits;  for (;;)    {      _gmp_rand (&ranm, rstate, BITS_PER_RANDCALL);      chunksize = 1 + ranm % cap_chunksize;      bi = (bi < chunksize) ? 0 : bi - chunksize;      if (bi == 0)	break;			/* low chunk is ...1 */      rp[bi / GMP_NUMB_BITS] ^= CNST_LIMB (1) << bi % GMP_NUMB_BITS;      _gmp_rand (&ranm, rstate, BITS_PER_RANDCALL);      chunksize = 1 + ranm % cap_chunksize;      bi = (bi < chunksize) ? 0 : bi - chunksize;      mpn_incr_u (rp + bi / GMP_NUMB_BITS, CNST_LIMB (1) << bi % GMP_NUMB_BITS);      if (bi == 0)	break;			/* low chunk is ...0 */    }}

voidcheck_twobits (void){#define TWOBITS(a, b) /  ((CNST_LIMB(1) << (a)) | (CNST_LIMB(1) << (b)))  refmpn_zero (x, SIZE);  x[0] = TWOBITS (1, 0);  check ();  refmpn_zero (x, SIZE);  x[0] = TWOBITS (GMP_NUMB_BITS-1, 1);  check ();  refmpn_zero (x, SIZE);  x[0] = CNST_LIMB(1);  x[1] = CNST_LIMB(1);  check ();  refmpn_zero (x, SIZE);  x[0] = CNST_LIMB(1) << (GMP_NUMB_BITS-1);  x[1] = CNST_LIMB(1);  check ();  refmpn_zero (x, SIZE);  x[1] = TWOBITS (1, 0);  check ();  refmpn_zero (x, SIZE);  x[1] = CNST_LIMB(1);  x[2] = CNST_LIMB(1);  check ();}

voidmpn_invert (mp_ptr ip, mp_srcptr dp, mp_size_t n, mp_ptr scratch){  ASSERT (n > 0);  ASSERT (dp[n-1] & GMP_NUMB_HIGHBIT);  ASSERT (! MPN_OVERLAP_P (ip, n, dp, n));  ASSERT (! MPN_OVERLAP_P (ip, n, scratch, mpn_invertappr_itch(n)));  ASSERT (! MPN_OVERLAP_P (dp, n, scratch, mpn_invertappr_itch(n)));  if (n == 1)    invert_limb (*ip, *dp);  else if (BELOW_THRESHOLD (n, INV_APPR_THRESHOLD))    {	/* Maximum scratch needed by this branch: 2*n */	mp_size_t i;	mp_ptr xp;	xp = scratch;				/* 2 * n limbs */	/* n > 1 here */	i = n;	do	  xp[--i] = GMP_NUMB_MAX;	while (i);	mpn_com (xp + n, dp, n);	if (n == 2) {	  mpn_divrem_2 (ip, 0, xp, 4, dp);	} else {	  gmp_pi1_t inv;	  invert_pi1 (inv, dp[n-1], dp[n-2]);	  /* FIXME: should we use dcpi1_div_q, for big sizes? */	  mpn_sbpi1_div_q (ip, xp, 2 * n, dp, n, inv.inv32);	}    }  else { /* Use approximated inverse; correct the result if needed. */      mp_limb_t e; /* The possible error in the approximate inverse */      ASSERT ( mpn_invert_itch (n) >= mpn_invertappr_itch (n) );      e = mpn_ni_invertappr (ip, dp, n, scratch);      if (UNLIKELY (e)) { /* Assume the error can only be "0" (no error) or "1". */	/* Code to detect and correct the "off by one" approximation. */	mpn_mul_n (scratch, ip, dp, n);	e = mpn_add_n (scratch, scratch, dp, n); /* FIXME: we only need e.*/	if (LIKELY(e)) /* The high part can not give a carry by itself. */	  e = mpn_add_nc (scratch + n, scratch + n, dp, n, e); /* FIXME:e */	/* If the value was wrong (no carry), correct it (increment). */	e ^= CNST_LIMB (1);	MPN_INCR_U (ip, n, e);      }  }}

/* Create a fake mpz consisting of just a single 1 bit, with totbits being   the total number of bits, inclusive of that 1 bit.  */voidmpz_fake_bits (mpz_ptr z, unsigned long totbits){  static mp_limb_t  n;  unsigned long     zero_bits, zero_limbs;  zero_bits = totbits - 1;  zero_limbs = zero_bits / GMP_NUMB_BITS;  zero_bits %= GMP_NUMB_BITS;  SIZ(z) = zero_limbs + 1;  PTR(z) = (&n) - (SIZ(z) - 1);  n = CNST_LIMB(1) << zero_bits;  ASSERT_ALWAYS (mpz_sizeinbase (z, 2) == totbits);}

/* Exercise values 2^n+1, while such a value fits the mantissa of a double. */voidcheck_twobit (void){  int        i, mant_bits;  double     got, want;  mp_size_t  nsize, sign;  mp_ptr     np;  mant_bits = tests_dbl_mant_bits ();  if (mant_bits == 0)    return;  np = refmpn_malloc_limbs (BITS_TO_LIMBS (mant_bits));  want = 3.0;  for (i = 1; i < mant_bits; i++)    {      nsize = BITS_TO_LIMBS (i+1);      refmpn_zero (np, nsize);      np[i/GMP_NUMB_BITS] = CNST_LIMB(1) << (i % GMP_NUMB_BITS);      np[0] |= 1;      for (sign = 0; sign >= -1; sign--)        {          got = mpn_get_d (np, nsize, sign, 0);          if (got != want)            {              printf    ("mpn_get_d wrong on 2^%d + 1/n", i);              printf    ("   sign     %ld/n", (long) sign);              mpn_trace ("   n        ", np, nsize);              printf    ("   nsize    %ld/n", (long) nsize);              d_trace   ("   want     ", want);              d_trace   ("   got      ", got);              abort();            }          want = -want;        }      want = 2.0 * want - 1.0;    }  free (np);}

voidgmp_randinit_lc_2exp (gmp_randstate_t rstate,		      mpz_srcptr a,		      unsigned long int c,		      mp_bitcnt_t m2exp){  gmp_rand_lc_struct *p;  mp_size_t seedn = BITS_TO_LIMBS (m2exp);  ASSERT_ALWAYS (m2exp != 0);  p = __GMP_ALLOCATE_FUNC_TYPE (1, gmp_rand_lc_struct);  RNG_STATE (rstate) = (void *) p;  RNG_FNPTR (rstate) = (void *) &Linear_Congruential_Generator;  /* allocate m2exp bits of space for p->_mp_seed, and initial seed "1" */  mpz_init2 (p->_mp_seed, m2exp);  MPN_ZERO (PTR (p->_mp_seed), seedn);  SIZ (p->_mp_seed) = seedn;  PTR (p->_mp_seed)[0] = 1;  /* "a", forced to 0 to 2^m2exp-1 */  mpz_init (p->_mp_a);  mpz_fdiv_r_2exp (p->_mp_a, a, m2exp);  /* Avoid SIZ(a) == 0 to avoid checking for special case in lc().  */  if (SIZ (p->_mp_a) == 0)    {      SIZ (p->_mp_a) = 1;      PTR (p->_mp_a)[0] = CNST_LIMB (0);    }  MPN_SET_UI (p->_cp, p->_cn, c);  /* Internally we may discard any bits of c above m2exp.  The following     code ensures that __GMPN_ADD in lc() will always work.  */  if (seedn < p->_cn)    p->_cn = (p->_cp[0] != 0);  p->_mp_m2exp = m2exp;}

void tc4_addmul_1(mp_ptr wp, mp_size_t * wn, mp_srcptr xp, mp_size_t xn, mp_limb_t y){  mp_size_t  sign, wu, xu, ws, new_wn, min_size, dsize;  mp_limb_t  cy;  /* w unaffected if x==0 or y==0 */  if (xn == 0 || y == 0)    return;  sign = xn;  xu = ABS (xn);  ws = *wn;  if (*wn == 0)  {      /* nothing to add to, just set x*y, "sign" gives the sign */      cy = mpn_mul_1 (wp, xp, xu, y);      if (cy)		{			wp[xu] = cy;         xu = xu + 1;		}       *wn = (sign >= 0 ? xu : -xu);      return;  }    sign ^= *wn;  wu = ABS (*wn);  new_wn = MAX (wu, xu);  min_size = MIN (wu, xu);  if (sign >= 0)  {      /* addmul of absolute values */      cy = mpn_addmul_1 (wp, xp, min_size, y);            dsize = xu - wu;#if HAVE_NATIVE_mpn_mul_1c      if (dsize > 0)        cy = mpn_mul_1c (wp + min_size, xp + min_size, dsize, y, cy);      else if (dsize < 0)      {          dsize = -dsize;          cy = mpn_add_1 (wp + min_size, wp + min_size, dsize, cy);      }#else      if (dsize != 0)      {          mp_limb_t cy2;          if (dsize > 0)            cy2 = mpn_mul_1 (wp + min_size, xp + min_size, dsize, y);          else          {              dsize = -dsize;              cy2 = 0;          }          cy = cy2 + mpn_add_1 (wp + min_size, wp + min_size, dsize, cy);      }#endif      if (cy)		{			wp[dsize + min_size] = cy;         new_wn ++;		}   } else   {      /* submul of absolute values */      cy = mpn_submul_1 (wp, xp, min_size, y);      if (wu >= xu)      {          /* if w bigger than x, then propagate borrow through it */          if (wu != xu)            cy = mpn_sub_1 (wp + xu, wp + xu, wu - xu, cy);          if (cy != 0)          {              /* Borrow out of w, take twos complement negative to get                 absolute value, flip sign of w.  */              wp[new_wn] = ~-cy;  /* extra limb is 0-cy */              mpn_not (wp, new_wn);              new_wn++;              MPN_INCR_U (wp, new_wn, CNST_LIMB(1));              ws = -*wn;          }      } else /* wu < xu */      {          /* x bigger than w, so want x*y-w.  Submul has given w-x*y, so             take twos complement and use an mpn_mul_1 for the rest.  */          mp_limb_t  cy2;          /* -(-cy*b^n + w-x*y) = (cy-1)*b^n + ~(w-x*y) + 1 */          mpn_not (wp, wu);          cy += mpn_add_1 (wp, wp, wu, CNST_LIMB(1));          cy -= 1;//.........这里部分代码省略.........

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

mp_limb_tr_string (const char *s){  const char  *s_orig = s;  long        n;  if (strcmp (s, "aas") == 0)    return GMP_NUMB_0xAA;  {    mpz_t      z;    mp_limb_t  l;    int        set, siz;    mpz_init (z);    set = mpz_set_str (z, s, 0);    siz = SIZ(z);    l = (siz == 0 ? 0 : siz > 0 ? PTR(z)[0] : -PTR(z)[0]);    mpz_clear (z);    if (set == 0)      {        if (siz > 1 || siz < -1)          printf ("Warning, r parameter %s truncated to %d bits/n",                  s_orig, BITS_PER_MP_LIMB);        return l;      }  }  if (s[0] == '0' && (s[1] == 'x' || s[1] == 'X'))    n = strtoul (s+2, (char **) &s, 16);  else    n = strtol (s, (char **) &s, 10);  if (strcmp (s, "bits") == 0)    {      mp_limb_t  l;      if (n > BITS_PER_MP_LIMB)        {          fprintf (stderr, "%ld bit parameter invalid (max %d bits)/n",                    n, BITS_PER_MP_LIMB);          exit (1);        }      mpn_random (&l, 1);      return (l | (CNST_LIMB(1) << (n-1))) & LIMB_ONES(n);    }  else  if (strcmp (s, "ones") == 0)    {      if (n > BITS_PER_MP_LIMB)        {          fprintf (stderr, "%ld bit parameter invalid (max %d bits)/n",                    n, BITS_PER_MP_LIMB);          exit (1);        }      return LIMB_ONES (n);    }  else if (*s != '/0')    {      fprintf (stderr, "invalid r parameter: %s/n", s_orig);      exit (1);    }  return n;}

REGPARM_ATTR (1) static voidcfdiv_r_2exp (mpz_ptr w, mpz_srcptr u, unsigned long cnt, int dir){  mp_size_t  usize, abs_usize, limb_cnt, i;  mp_srcptr  up;  mp_ptr     wp;  mp_limb_t  high;  usize = SIZ(u);  if (usize == 0)    {      SIZ(w) = 0;      return;    }  limb_cnt = cnt / GMP_NUMB_BITS;  cnt %= GMP_NUMB_BITS;  abs_usize = ABS (usize);  /* MPZ_REALLOC(w) below is only when w!=u, so we can fetch PTR(u) here     nice and early */  up = PTR(u);  if ((usize ^ dir) < 0)    {      /* Round towards zero, means just truncate */      if (w == u)        {          /* if already smaller than limb_cnt then do nothing */          if (abs_usize <= limb_cnt)            return;          wp = PTR(w);        }      else        {          i = MIN (abs_usize, limb_cnt+1);          MPZ_REALLOC (w, i);          wp = PTR(w);          MPN_COPY (wp, up, i);          /* if smaller than limb_cnt then only the copy is needed */          if (abs_usize <= limb_cnt)            {              SIZ(w) = usize;              return;            }        }    }  else    {      /* Round away from zero, means twos complement if non-zero */      /* if u!=0 and smaller than divisor, then must negate */      if (abs_usize <= limb_cnt)        goto negate;      /* if non-zero low limb, then must negate */      for (i = 0; i < limb_cnt; i++)        if (up[i] != 0)          goto negate;      /* if non-zero partial limb, then must negate */      if ((up[limb_cnt] & LOW_MASK (cnt)) != 0)        goto negate;      /* otherwise low bits of u are zero, so that's the result */      SIZ(w) = 0;      return;    negate:      /* twos complement negation to get 2**cnt-u */      MPZ_REALLOC (w, limb_cnt+1);      up = PTR(u);      wp = PTR(w);      /* Ones complement */      i = MIN (abs_usize, limb_cnt+1);      mpn_com_n (wp, up, i);      for ( ; i <= limb_cnt; i++)        wp[i] = GMP_NUMB_MAX;      /* Twos complement.  Since u!=0 in the relevant part, the twos         complement never gives 0 and a carry, so can use MPN_INCR_U. */      MPN_INCR_U (wp, limb_cnt+1, CNST_LIMB(1));      usize = -usize;    }  /* Mask the high limb */  high = wp[limb_cnt];  high &= LOW_MASK (cnt);  wp[limb_cnt] = high;  /* Strip any consequent high zeros */  while (high == 0)    {      limb_cnt--;      if (limb_cnt < 0)//.........这里部分代码省略.........

static voidmpf_ceil_or_floor (mpf_ptr r, mpf_srcptr u, int dir){  mp_ptr     rp, up, p;  mp_size_t  size, asize, prec;  mp_exp_t   exp;  size = SIZ(u);  if (size == 0)    {    zero:      SIZ(r) = 0;      EXP(r) = 0;      return;    }  rp = PTR(r);  exp = EXP(u);  if (exp <= 0)    {      /* u is only a fraction */      if ((size ^ dir) < 0)        goto zero;      rp[0] = 1;      EXP(r) = 1;      SIZ(r) = dir;      return;    }  EXP(r) = exp;  up = PTR(u);  asize = ABS (size);  up += asize;  /* skip fraction part of u */  asize = MIN (asize, exp);  /* don't lose precision in the copy */  prec = PREC (r) + 1;  /* skip excess over target precision */  asize = MIN (asize, prec);  up -= asize;  if ((size ^ dir) >= 0)    {      /* rounding direction matches sign, must increment if ignored part is         non-zero */      for (p = PTR(u); p != up; p++)        {          if (*p != 0)            {              if (mpn_add_1 (rp, up, asize, CNST_LIMB(1)))                {                  /* was all 0xFF..FFs, which have become zeros, giving just                     a carry */                  rp[0] = 1;                  asize = 1;                  EXP(r)++;                }              SIZ(r) = (size >= 0 ? asize : -asize);              return;            }        }    }  SIZ(r) = (size >= 0 ? asize : -asize);  if (rp != up)    MPN_COPY_INCR (rp, up, asize);}

voidcheck_limbdata (void){#define M  GMP_NUMB_MAX    static const struct {    mp_exp_t       exp;    mp_size_t      size;    mp_limb_t      d[10];    unsigned long  want;  } data[] = {    /* in the comments here, a "_" indicates a digit (ie. limb) position not       included in the d data, and therefore zero */    { 0, 0, { 0 }, 0L },    /* 0 */    { 1,  1, { 1 }, 1L },   /* 1 */    { 1, -1, { 1 }, -1L },  /* -1 */    { 0,  1, { 1 }, 0L },   /* .1 */    { 0, -1, { 1 }, 0L },   /* -.1 */    { -1,  1, { 1 }, 0L },  /* ._1 */    { -1, -1, { 1 }, 0L },  /* -._1 */    { -999,          1, { 1 }, 0L },   /* .___1 small */    { MP_EXP_T_MIN,  1, { 1 }, 0L },   /* .____1 very small */    { 999,          1, { 1 }, 0L },    /* 1____. big */    { MP_EXP_T_MAX, 1, { 1 }, 0L },    /* 1_____. very big */    { 1, 2, { 999, 2 }, 2L },                  /* 2.9 */    { 5, 8, { 7, 8, 9, 3, 0, 0, 0, 1 }, 3L },  /* 10003.987 */    { 2, 2, { M, M },    LONG_MAX }, /* FF. */    { 2, 2, { M, M, M }, LONG_MAX }, /* FF.F */    { 3, 3, { M, M, M }, LONG_MAX }, /* FFF. */#if GMP_NUMB_BITS >= BITS_PER_ULONG    /* normal case, numb bigger than long */    { 2,  1, { 1 },    0L },      /* 1_. */    { 2,  2, { 0, 1 }, 0L },      /* 10. */    { 2,  2, { 999, 1 }, 999L },  /* 19. */    { 3,  2, { 999, 1 }, 0L },    /* 19_. */#else    /* nails case, numb smaller than long */    { 2,  1, { 1 }, 1L << GMP_NUMB_BITS },  /* 1_. */    { 3,  1, { 1 }, 0L },                   /* 1__. */    { 2,  2, { 99, 1 },    99L + (1L << GMP_NUMB_BITS) },  /* 19. */    { 3,  2, { 1, 99 },    1L << GMP_NUMB_BITS },          /* 91_. */    { 3,  3, { 0, 1, 99 }, 1L << GMP_NUMB_BITS },          /* 910. */#endif  };  mpf_t          f;  unsigned long  got;  int            i;  mp_limb_t      buf[20 + numberof(data[i].d)];  for (i = 0; i < numberof (data); i++)    {      refmpn_fill (buf, 10, CNST_LIMB(0xDEADBEEF));      refmpn_copy (buf+10, data[i].d, ABS(data[i].size));      refmpn_fill (buf+10+ABS(data[i].size), 10, CNST_LIMB(0xDEADBEEF));      PTR(f) = buf+10;      EXP(f) = data[i].exp;      SIZ(f) = data[i].size;      PREC(f) = numberof (data[i].d);      MPF_CHECK_FORMAT (f);      got = mpf_get_si (f);      if (got != data[i].want)        {          printf    ("mpf_get_si wrong at limb data[%d]/n", i);          mpf_trace ("  f", f);          mpn_trace ("  d", data[i].d, data[i].size);          printf    ("  size %ld/n", (long) data[i].size);          printf    ("  exp %ld/n", (long) data[i].exp);          printf    ("  got   %lu (0x%lX)/n", got, got);          printf    ("  want  %lu (0x%lX)/n", data[i].want, data[i].want);          abort();        }    }}

//.........这里部分代码省略.........	  xp[n] = mpn_mul_fft (xp, n, ap1, anp, ap1, anp, k);	else if (UNLIKELY (ap1 == a0))	  {	    ASSERT (anp <= n);	    ASSERT (2*anp > n);	    mpn_sqr (xp, a0, an);	    anp = 2*an - n;	    cy = mpn_sub (xp, xp, n, xp + n, anp);	    xp[n] = 0;	    MPN_INCR_U (xp, n+1, cy);	  }	else	  mpn_bc_sqrmod_bnp1 (xp, ap1, n, xp);      }      /* Here the CRT recomposition begins.	 xm <- (xp + xm)/2 = (xp + xm)B^n/2 mod (B^n-1)	 Division by 2 is a bitwise rotation.	 Assumes xp normalised mod (B^n+1).	 The residue class [0] is represented by [B^n-1]; except when	 both input are ZERO.      */#if HAVE_NATIVE_mpn_rsh1add_n || HAVE_NATIVE_mpn_rsh1add_nc#if HAVE_NATIVE_mpn_rsh1add_nc      cy = mpn_rsh1add_nc(rp, rp, xp, n, xp[n]); /* B^n = 1 */      hi = cy << (GMP_NUMB_BITS - 1);      cy = 0;      /* next update of rp[n-1] will set cy = 1 only if rp[n-1]+=hi	 overflows, i.e. a further increment will not overflow again. */#else /* ! _nc */      cy = xp[n] + mpn_rsh1add_n(rp, rp, xp, n); /* B^n = 1 */      hi = (cy<<(GMP_NUMB_BITS-1))&GMP_NUMB_MASK; /* (cy&1) << ... */      cy >>= 1;      /* cy = 1 only if xp[n] = 1 i.e. {xp,n} = ZERO, this implies that	 the rsh1add was a simple rshift: the top bit is 0. cy=1 => hi=0. */#endif#if GMP_NAIL_BITS == 0      add_ssaaaa(cy, rp[n-1], cy, rp[n-1], CNST_LIMB(0), hi);#else      cy += (hi & rp[n-1]) >> (GMP_NUMB_BITS-1);      rp[n-1] ^= hi;#endif#else /* ! HAVE_NATIVE_mpn_rsh1add_n */#if HAVE_NATIVE_mpn_add_nc      cy = mpn_add_nc(rp, rp, xp, n, xp[n]);#else /* ! _nc */      cy = xp[n] + mpn_add_n(rp, rp, xp, n); /* xp[n] == 1 implies {xp,n} == ZERO */#endif      cy += (rp[0]&1);      mpn_rshift(rp, rp, n, 1);      ASSERT (cy <= 2);      hi = (cy<<(GMP_NUMB_BITS-1))&GMP_NUMB_MASK; /* (cy&1) << ... */      cy >>= 1;      /* We can have cy != 0 only if hi = 0... */      ASSERT ((rp[n-1] & GMP_NUMB_HIGHBIT) == 0);      rp[n-1] |= hi;      /* ... rp[n-1] + cy can not overflow, the following INCR is correct. */#endif      ASSERT (cy <= 1);      /* Next increment can not overflow, read the previous comments about cy. */      ASSERT ((cy == 0) || ((rp[n-1] & GMP_NUMB_HIGHBIT) == 0));      MPN_INCR_U(rp, n, cy);      /* Compute the highest half:	 ([(xp + xm)/2 mod (B^n-1)] - xp ) * B^n       */      if (UNLIKELY (2*an < rn))	{	  /* Note that in this case, the only way the result can equal	     zero mod B^{rn} - 1 is if the input is zero, and	     then the output of both the recursive calls and this CRT	     reconstruction is zero, not B^{rn} - 1. */	  cy = mpn_sub_n (rp + n, rp, xp, 2*an - n);	  /* FIXME: This subtraction of the high parts is not really	     necessary, we do it to get the carry out, and for sanity	     checking. */	  cy = xp[n] + mpn_sub_nc (xp + 2*an - n, rp + 2*an - n,				   xp + 2*an - n, rn - 2*an, cy);	  ASSERT (mpn_zero_p (xp + 2*an - n+1, rn - 1 - 2*an));	  cy = mpn_sub_1 (rp, rp, 2*an, cy);	  ASSERT (cy == (xp + 2*an - n)[0]);	}      else	{	  cy = xp[n] + mpn_sub_n (rp + n, rp, xp, n);	  /* cy = 1 only if {xp,n+1} is not ZERO, i.e. {rp,n} is not ZERO.	     DECR will affect _at most_ the lowest n limbs. */	  MPN_DECR_U (rp, 2*n, cy);	}#undef a0#undef a1#undef xp#undef sp1    }}

//.........这里部分代码省略.........	      MPN_COPY (rp, vp, vsize);	      rsize = vsize;	      negate ^= 1;	      goto done;	    }	  if (up[0] != 0)	    break;	  up++, usize--;	}      ASSERT (usize > 0 && vsize > 0);      TMP_MARK;      tp = TMP_ALLOC_LIMBS (prec);      /* uuuu     |  uuuu     |  uuuu     |  uuuu     |  uuuu    */      /* vvvvvvv  |  vv       |    vvvvv  |    v      |       vv */      if (usize > ediff)	{	  /* U and V partially overlaps.  */	  if (ediff == 0)	    {	      ASSERT (usize == 1 && vsize >= 1 && ulimb == *up); /* usize is 1>ediff, vsize >= 1 */	      if (1 < vsize)		{		  /* u        */		  /* vvvvvvv  */		  rsize = vsize;		  vsize -= 1;		  /* mpn_cmp (up, vp + vsize - usize, usize) > 0 */		  if (ulimb > vp[vsize])		    {		      tp[vsize] = ulimb - vp[vsize] - 1;		      ASSERT_CARRY (mpn_neg (tp, vp, vsize));		    }		  else		    {		      /* vvvvvvv  */  /* Swap U and V. */		      /* u        */		      MPN_COPY (tp, vp, vsize);		      tp[vsize] = vp[vsize] - ulimb;		      negate = 1;		    }		}	      else /* vsize == usize == 1 */		{		  /* u     */		  /* v     */		  rsize = 1;		  negate = ulimb < vp[0];		  tp[0] = negate ? vp[0] - ulimb: ulimb - vp[0];		}	    }	  else	    {	      ASSERT (vsize + ediff <= usize);	      ASSERT (vsize == 1 && usize >= 2 && ulimb == *vp);		{		  /* uuuu     */		  /*   v      */		  mp_size_t size;		  size = usize - ediff - 1;		  MPN_COPY (tp, up, size);		  ASSERT_NOCARRY (mpn_sub_1 (tp + size, up + size, usize - size, ulimb));		  rsize = usize;		}		/* Other cases are not possible */		/* uuuu     */		/*   vvvvv  */	    }	}      else	{	  /* uuuu     */	  /*      vv  */	  mp_size_t size, i;	  ASSERT_CARRY (mpn_neg (tp, vp, vsize));	  rsize = vsize + ediff;	  size = rsize - usize;	  for (i = vsize; i < size; i++)	    tp[i] = GMP_NUMB_MAX;	  ASSERT_NOCARRY (mpn_sub_1 (tp + size, up, usize, CNST_LIMB (1)));	}      /* Full normalize.  Optimize later.  */      while (rsize != 0 && tp[rsize - 1] == 0)	{	  rsize--;	  uexp--;	}      MPN_COPY (rp, tp, rsize);      TMP_FREE;    } done:  r->_mp_size = negate ? -rsize : rsize;  r->_mp_exp = uexp;#endif}

static unsigned long intlc (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;}

mp_limb_tmpn_div_qr_1n_pi1 (mp_ptr qp, mp_srcptr up, mp_size_t n, mp_limb_t u1,		   mp_limb_t d, mp_limb_t dinv){  mp_limb_t B2;  mp_limb_t u0, u2;  mp_limb_t q0, q1;  mp_limb_t p0, p1;  mp_limb_t t;  mp_size_t j;  ASSERT (d & GMP_LIMB_HIGHBIT);  ASSERT (n > 0);  ASSERT (u1 < d);  if (n == 1)    {      udiv_qrnnd_preinv (qp[0], u1, u1, up[0], d, dinv);      return u1;    }  /* FIXME: Could be precomputed */  B2 = -d*dinv;  umul_ppmm (q1, q0, dinv, u1);  umul_ppmm (p1, p0, B2, u1);  q1 += u1;  ASSERT (q1 >= u1);  u0 = up[n-1];	/* Early read, to allow qp == up. */  qp[n-1] = q1;  add_mssaaaa (u2, u1, u0, u0, up[n-2], p1, p0);  /* FIXME: Keep q1 in a variable between iterations, to reduce number     of memory accesses. */  for (j = n-2; j-- > 0; )    {      mp_limb_t q2, cy;      /* Additions for the q update:       *	+-------+       *        |u1 * v |       *        +---+---+       *        | u1|       *    +---+---+       *    | 1 | v |  (conditional on u2)       *    +---+---+       *        | 1 |  (conditional on u0 + u2 B2 carry)       *        +---+       * +      | q0|       *   -+---+---+---+       *    | q2| q1| q0|       *    +---+---+---+      */      umul_ppmm (p1, t, u1, dinv);      add_ssaaaa (q2, q1, -u2, u2 & dinv, CNST_LIMB(0), u1);      add_ssaaaa (q2, q1, q2, q1, CNST_LIMB(0), p1);      add_ssaaaa (q2, q1, q2, q1, CNST_LIMB(0), q0);      q0 = t;      umul_ppmm (p1, p0, u1, B2);      ADDC_LIMB (cy, u0, u0, u2 & B2);      u0 -= (-cy) & d;      /* Final q update */      add_ssaaaa (q2, q1, q2, q1, CNST_LIMB(0), cy);      qp[j+1] = q1;      MPN_INCR_U (qp+j+2, n-j-2, q2);      add_mssaaaa (u2, u1, u0, u0, up[j], p1, p0);    }  q1 = (u2 > 0);  u1 -= (-q1) & d;  t = (u1 >= d);  q1 += t;  u1 -= (-t) & d;  udiv_qrnnd_preinv (t, u0, u1, u0, d, dinv);  add_ssaaaa (q1, q0, q1, q0, CNST_LIMB(0), t);  MPN_INCR_U (qp+1, n-1, q1);  qp[0] = q0;  return u0;}

voidmpn_toom4_sqr_n (mp_ptr rp, mp_srcptr up, mp_size_t n){  mp_size_t len1, ind;  mp_limb_t cy, r30, r31;  mp_ptr tp;  mp_size_t a0n, a1n, a2n, a3n, sn, n1, n2, n3, n4, n5, n6, n7, n8, n9, rpn, t4;  len1 = n;  ASSERT (n >= 1);  MPN_NORMALIZE(up, len1);    sn = (n - 1) / 4 + 1;  /* a0 - a3 are defined in mpn_toom4_mul_n above */     TC4_NORM(a0, a0n, sn);	TC4_NORM(a1, a1n, sn);	TC4_NORM(a2, a2n, sn);	TC4_NORM(a3, a3n, n - 3*sn);    t4 = 2*sn+2; // allows mult of 2 integers of sn + 1 limbs   tp = __GMP_ALLOCATE_FUNC_LIMBS(4*t4 + 4*(sn + 1));   tc4_add_unsigned(u5, &n5, a3, a3n, a1, a1n);    tc4_add_unsigned(u4, &n4, a2, a2n, a0, a0n); 	tc4_add_unsigned(u2, &n2, u4, n4, u5, n5);    tc4_sub(u3, &n3, u4, n4, u5, n5);	SQR_TC4(r4, n4, u3, n3);   SQR_TC4_UNSIGNED(r3, n3, u2, n2);		tc4_lshift(r1, &n1, a0, a0n, 3);	tc4_addlsh1_unsigned(r1, &n1, a2, a2n); 	tc4_lshift(r2, &n8, a1, a1n, 2);   tc4_add(r2, &n8, r2, n8, a3, a3n);   tc4_add(u4, &n9, r1, n1, r2, n8);   tc4_sub(u5, &n5, r1, n1, r2, n8);   	r30 = r3[0];	if (!n3) r30 = CNST_LIMB(0);   r31 = r3[1];	SQR_TC4(r6, n6, u5, n5);   SQR_TC4_UNSIGNED(r5, n5, u4, n9);   r3[1] = r31;   tc4_lshift(u2, &n8, a3, a3n, 3);   tc4_addmul_1(u2, &n8, a2, a2n, 4);	tc4_addlsh1_unsigned(u2, &n8, a1, a1n);	tc4_add(u2, &n8, u2, n8, a0, a0n);   	SQR_TC4_UNSIGNED(r2, n2, u2, n8);   SQR_TC4_UNSIGNED(r1, n1, a3, a3n);   SQR_TC4_UNSIGNED(r7, n7, a0, a0n);	TC4_DENORM(r1, n1,  t4 - 1);   TC4_DENORM(r2, n2,  t4 - 1);   if (n3)     TC4_DENORM(r3, n3,  t4 - 1);   else {     /* MPN_ZERO defeats gcc 4.1.2 here, hence the explicit for loop */     for (ind = 1 ; ind < t4 - 1; ind++)        (r3)[ind] = CNST_LIMB(0);   }   TC4_DENORM(r4, n4,  t4 - 1);   TC4_DENORM(r5, n5,  t4 - 1);   TC4_DENORM(r6, n6,  t4 - 1);   TC4_DENORM(r7, n7,  t4 - 2); // we treat r7 differently (it cannot exceed t4-2 in length)/*	rp        rp1          rp2           rp3          rp4           rp5         rp6           rp7<----------- r7-----------><------------r5-------------->                                                                   <-------------r3------------->              <-------------r6------------->                        < -----------r2------------>{           }                                         <-------------r4-------------->         <--------------r1---->*/	mpn_toom4_interpolate(rp, &rpn, sn, tp, t4 - 1, n4, n6, r30);	if (rpn != 2*n) 	{		MPN_ZERO((rp + rpn), 2*n - rpn);	}   __GMP_FREE_FUNC_LIMBS (tp, 4*t4 + 4*(sn+1));}

#include "mpir.h"#include "gmp-impl.h"#if GMP_NUMB_BITS == 32const struct bases mp_bases[257] ={  /*   0 */ { 0, 0.0, 0 },  /*   1 */ { 0, 1e37, 0 },  /*   2 */ { 32, 1.0000000000000000, 0x1 },  /*   3 */ { 20, 0.6309297535714574, CNST_LIMB(0xcfd41b91), CNST_LIMB(0x3b563c24) },  /*   4 */ { 16, 0.5000000000000000, 0x2 },  /*   5 */ { 13, 0.4306765580733931, CNST_LIMB(0x48c27395), CNST_LIMB(0xc25c2684) },  /*   6 */ { 12, 0.3868528072345416, CNST_LIMB(0x81bf1000), CNST_LIMB(0xf91bd1b6) },  /*   7 */ { 11, 0.3562071871080222, CNST_LIMB(0x75db9c97), CNST_LIMB(0x1607a2cb) },  /*   8 */ { 10, 0.3333333333333334, 0x3 },  /*   9 */ { 10, 0.3154648767857287, CNST_LIMB(0xcfd41b91), CNST_LIMB(0x3b563c24) },  /*  10 */ { 9, 0.3010299956639811, CNST_LIMB(0x3b9aca00), CNST_LIMB(0x12e0be82) },  /*  11 */ { 9, 0.2890648263178878, CNST_LIMB(0x8c8b6d2b), CNST_LIMB(0xd24cde04) },  /*  12 */ { 8, 0.2789429456511298, CNST_LIMB(0x19a10000), CNST_LIMB(0x3fa39ab5) },  /*  13 */ { 8, 0.2702381544273197, CNST_LIMB(0x309f1021), CNST_LIMB(0x50f8ac5f) },  /*  14 */ { 8, 0.2626495350371936, CNST_LIMB(0x57f6c100), CNST_LIMB(0x74843b1e) },  /*  15 */ { 8, 0.2559580248098155, CNST_LIMB(0x98c29b81), CNST_LIMB(0xad0326c2) },  /*  16 */ { 8, 0.2500000000000000, 0x4 },  /*  17 */ { 7, 0.2446505421182260, CNST_LIMB(0x18754571), CNST_LIMB(0x4ef0b6bd) },  /*  18 */ { 7, 0.2398124665681315, CNST_LIMB(0x247dbc80), CNST_LIMB(0xc0fc48a1) },  /*  19 */ { 7, 0.2354089133666382, CNST_LIMB(0x3547667b), CNST_LIMB(0x33838942) },  /*  20 */ { 7, 0.2313782131597592, CNST_LIMB(0x4c4b4000), CNST_LIMB(0xad7f29ab) },  /*  21 */ { 7, 0.2276702486969530, CNST_LIMB(0x6b5a6e1d), CNST_LIMB(0x313c3d15) },  /*  22 */ { 7, 0.2242438242175754, CNST_LIMB(0x94ace180), CNST_LIMB(0xb8cca9e0) },

/*   Toom 4 interpolation. Interpolates the value at 2^(sn*B) of a 	polynomial p(x) with 7 coefficients given the values 	p(oo), p(2), p(1), p(-1), 2^6*p(1/2), 2^6*p(-1/2), p(0).	The output is placed in rp and the final number of limbs of the	output is given in rpn.	The 4th and 6th values may be negative, and if so, n4 and n6 	should be set to a negative value respectively.   To save space we pass r3, r5, r7 in place in the output rp.	The other r's are stored separately in space tp.	The low limb of r3 is stored in r30, as it will be overwritten	by the high limb of r5.rp          rp1          rp2           rp3          rp4           rp5         rp6           rp7<----------- r7-----------><------------r5-------------->                                                                  <-------------r3------------->   We assume that r1 is stored at tp, r2 at (tp + t4), r4 at (tp + 2*t4) 	and r6 (tp + 3*t4). Each of these r's has t4 = s4 + 1 limbs allocated.*/void mpn_toom4_interpolate(mp_ptr rp, mp_size_t * rpn, mp_size_t sn,  		       mp_ptr tp, mp_size_t s4, mp_size_t n4, mp_size_t n6, mp_limb_t r30){	mp_size_t n1, n2, n3, n5, n7, t4;	mp_limb_t saved, saved2, cy;   t4 = s4 + 1;    	mpn_add_n(r2, r2, r5, s4);   if (n6 < 0) 		mpn_add_n(r6, r5, r6, s4);	else      mpn_sub_n(r6, r5, r6, s4);	/* r6 is now in twos complement format */	saved = r3[0];	r3[0] = r30;	if (n4 < 0) 		mpn_add_n(r4, r3, r4, s4);	else      mpn_sub_n(r4, r3, r4, s4);	r3[0] = saved;	/* r4 is now in twos complement format */		mpn_sub_n(r5, r5, r1, s4);#if HAVE_NATIVE_mpn_sublsh_n	r5[s4-1] -= mpn_sublsh_n(r5, r5, r7, s4-1, 6);#else	r5[s4-1] -= mpn_submul_1(r5, r7, s4-1, 64);#endif      TC4_RSHIFT1(r4, s4); 		saved = r3[0];	r3[0] = r30;	mpn_sub_n(r3, r3, r4, s4);	r30 = r3[0];	r3[0] = saved;	mpn_double(r5, s4); 	mpn_sub_n(r5, r5, r6, s4);   saved = r3[0];	r3[0] = r30;	mpn_submul_1(r2, r3, s4, 65);   r3[0] = saved;		saved2 = r7[s4-1];	r7[s4-1] = CNST_LIMB(0); // r7 is always positive so no sign extend needed	saved = r3[0];	r3[0] = r30;#if HAVE_NATIVE_mpn_subadd_n	mpn_subadd_n(r3, r3, r7, r1, s4);#else    mpn_sub_n(r3, r3, r7, s4);    mpn_sub_n(r3, r3, r1, s4);#endif	r7[s4-1] = saved2;   r30 = r3[0];	   mpn_addmul_1(r2, r3, s4, 45);#if HAVE_NATIVE_mpn_sublsh_n   cy = mpn_sublsh_n(r5, r5, r3, s4 - 1, 3);#else   cy = mpn_submul_1(r5, r3, s4 - 1, 8);#endif   r3[0] = saved;	r3[0] -= (cy + 8*r3[s4-1]);   	mpn_rshift(r5, r5, s4, 3); 	mpn_divexact_by3(r5, r5, s4);    	mpn_sub_n(r6, r6, r2, s4);#if HAVE_NATIVE_mpn_sublsh_n//.........这里部分代码省略.........

intmpz_congruent_2exp_p (mpz_srcptr a, mpz_srcptr c, mp_bitcnt_t d){  mp_size_t      i, dlimbs;  unsigned       dbits;  mp_ptr         ap, cp;  mp_limb_t      dmask, alimb, climb, sum;  mp_size_t      asize_signed, csize_signed, asize, csize;  if (ABSIZ(a) < ABSIZ(c))    MPZ_SRCPTR_SWAP (a, c);  dlimbs = d / GMP_NUMB_BITS;  dbits = d % GMP_NUMB_BITS;  dmask = (CNST_LIMB(1) << dbits) - 1;  ap = PTR(a);  cp = PTR(c);  asize_signed = SIZ(a);  asize = ABS(asize_signed);  csize_signed = SIZ(c);  csize = ABS(csize_signed);  if (csize_signed == 0)    goto a_zeros;  if ((asize_signed ^ csize_signed) >= 0)    {      /* same signs, direct comparison */      /* a==c for limbs in common */      if (mpn_cmp (ap, cp, MIN (csize, dlimbs)) != 0)        return 0;      /* if that's all of dlimbs, then a==c for remaining bits */      if (csize > dlimbs)        return ((ap[dlimbs]-cp[dlimbs]) & dmask) == 0;    a_zeros:      /* a remains, need all zero bits */      /* if d covers all of a and c, then must be exactly equal */      if (asize <= dlimbs)        return asize == csize;      /* whole limbs zero */      for (i = csize; i < dlimbs; i++)        if (ap[i] != 0)          return 0;      /* partial limb zero */      return (ap[dlimbs] & dmask) == 0;    }  else    {      /* different signs, negated comparison */      /* common low zero limbs, stopping at first non-zeros, which must         match twos complement */      i = 0;      for (;;)        {          ASSERT (i < csize);  /* always have a non-zero limb on c */          alimb = ap[i];          climb = cp[i];          sum = (alimb + climb) & GMP_NUMB_MASK;          if (i >= dlimbs)            return (sum & dmask) == 0;          i++;          /* require both zero, or first non-zeros as twos-complements */          if (sum != 0)            return 0;          if (alimb != 0)            break;        }      /* further limbs matching as ones-complement */      for (;;)        {          if (i >= csize)            break;          alimb = ap[i];          climb = cp[i];          sum = (alimb + climb + 1) & GMP_NUMB_MASK;          if (i >= dlimbs)            return (sum & dmask) == 0;          if (sum != 0)            return 0;          i++;        }//.........这里部分代码省略.........

/* Evaluates a polynomial of degree k > 2, in the points +2^shift and -2^shift. */intmpn_toom_eval_pm2exp (mp_ptr xp2, mp_ptr xm2, unsigned k,		      mp_srcptr xp, mp_size_t n, mp_size_t hn, unsigned shift,		      mp_ptr tp){  unsigned i;  int neg;#ifdef HAVE_NATIVE_mpn_addlsh_n  mp_limb_t cy;#endif  ASSERT (k >= 3);  ASSERT (shift*k < GMP_NUMB_BITS);  ASSERT (hn > 0);  ASSERT (hn <= n);  /* The degree k is also the number of full-size coefficients, so   * that last coefficient, of size hn, starts at xp + k*n. */#ifdef HAVE_NATIVE_mpn_addlsh_n  xp2[n] = mpn_addlsh_n (xp2, xp, xp + 2*n, n, 2*shift);  for (i = 4; i < k; i += 2)    xp2[n] += mpn_addlsh_n (xp2, xp2, xp + i*n, n, i*shift);  tp[n] = mpn_lshift (tp, xp+n, n, shift);  for (i = 3; i < k; i+= 2)    tp[n] += mpn_addlsh_n (tp, tp, xp+i*n, n, i*shift);  if (k & 1)    {      cy = mpn_addlsh_n (tp, tp, xp+k*n, hn, k*shift);      MPN_INCR_U (tp + hn, n+1 - hn, cy);    }  else    {      cy = mpn_addlsh_n (xp2, xp2, xp+k*n, hn, k*shift);      MPN_INCR_U (xp2 + hn, n+1 - hn, cy);    }#else /* !HAVE_NATIVE_mpn_addlsh_n */  xp2[n] = mpn_lshift (tp, xp+2*n, n, 2*shift);  xp2[n] += mpn_add_n (xp2, xp, tp, n);  for (i = 4; i < k; i += 2)    {      xp2[n] += mpn_lshift (tp, xp + ((mp_size_t) i)*n, n, i*shift);      xp2[n] += mpn_add_n (xp2, xp2, tp, n);    }  tp[n] = mpn_lshift (tp, xp+n, n, shift);  for (i = 3; i < k; i+= 2)    {      tp[n] += mpn_lshift (xm2, xp + ((mp_size_t) i)*n, n, i*shift);      tp[n] += mpn_add_n (tp, tp, xm2, n);    }  xm2[hn] = mpn_lshift (xm2, xp + ((mp_size_t) k)*n, hn, k*shift);  if (k & 1)    mpn_add (tp, tp, n+1, xm2, hn+1);  else    mpn_add (xp2, xp2, n+1, xm2, hn+1);#endif /* !HAVE_NATIVE_mpn_addlsh_n */  neg = (mpn_cmp (xp2, tp, n + 1) < 0) ? ~0 : 0;#ifdef HAVE_NATIVE_mpn_add_n_sub_n  if (neg)    mpn_add_n_sub_n (xp2, xm2, tp, xp2, n + 1);  else    mpn_add_n_sub_n (xp2, xm2, xp2, tp, n + 1);#else /* !HAVE_NATIVE_mpn_add_n_sub_n */  if (neg)    mpn_sub_n (xm2, tp, xp2, n + 1);  else    mpn_sub_n (xm2, xp2, tp, n + 1);  mpn_add_n (xp2, xp2, tp, n + 1);#endif /* !HAVE_NATIVE_mpn_add_n_sub_n */  /* FIXME: the following asserts are useless if (k+1)*shift >= GMP_LIMB_BITS */  ASSERT ((k+1)*shift >= GMP_LIMB_BITS ||	  xp2[n] < ((CNST_LIMB(1)<<((k+1)*shift))-1)/((CNST_LIMB(1)<<shift)-1));  ASSERT ((k+2)*shift >= GMP_LIMB_BITS ||	  xm2[n] < ((CNST_LIMB(1)<<((k+2)*shift))-((k&1)?(CNST_LIMB(1)<<shift):1))/((CNST_LIMB(1)<<(2*shift))-1));  return neg;}

voidcheck_functions (void){  mp_limb_t  wp[2], wp2[2], xp[2], yp[2], r;  int  i;  memcpy (&__gmpn_cpuvec, &initial_cpuvec, sizeof (__gmpn_cpuvec));  for (i = 0; i < 2; i++)    {      xp[0] = 123;      yp[0] = 456;      mpn_add_n (wp, xp, yp, (mp_size_t) 1);      ASSERT_ALWAYS (wp[0] == 579);    }  memcpy (&__gmpn_cpuvec, &initial_cpuvec, sizeof (__gmpn_cpuvec));  for (i = 0; i < 2; i++)    {      xp[0] = 123;      wp[0] = 456;      r = mpn_addmul_1 (wp, xp, (mp_size_t) 1, CNST_LIMB(2));      ASSERT_ALWAYS (wp[0] == 702);      ASSERT_ALWAYS (r == 0);    }#if HAVE_NATIVE_mpn_copyd  memcpy (&__gmpn_cpuvec, &initial_cpuvec, sizeof (__gmpn_cpuvec));  for (i = 0; i < 2; i++)    {      xp[0] = 123;      xp[1] = 456;      mpn_copyd (xp+1, xp, (mp_size_t) 1);      ASSERT_ALWAYS (xp[1] == 123);    }#endif#if HAVE_NATIVE_mpn_copyi  memcpy (&__gmpn_cpuvec, &initial_cpuvec, sizeof (__gmpn_cpuvec));  for (i = 0; i < 2; i++)    {      xp[0] = 123;      xp[1] = 456;      mpn_copyi (xp, xp+1, (mp_size_t) 1);      ASSERT_ALWAYS (xp[0] == 456);    }#endif  memcpy (&__gmpn_cpuvec, &initial_cpuvec, sizeof (__gmpn_cpuvec));  for (i = 0; i < 2; i++)    {      xp[0] = 1605;      mpn_divexact_1 (wp, xp, (mp_size_t) 1, CNST_LIMB(5));      ASSERT_ALWAYS (wp[0] == 321);    }  memcpy (&__gmpn_cpuvec, &initial_cpuvec, sizeof (__gmpn_cpuvec));  for (i = 0; i < 2; i++)    {      xp[0] = 1296;      r = mpn_divexact_by3c (wp, xp, (mp_size_t) 1, CNST_LIMB(0));      ASSERT_ALWAYS (wp[0] == 432);      ASSERT_ALWAYS (r == 0);    }  memcpy (&__gmpn_cpuvec, &initial_cpuvec, sizeof (__gmpn_cpuvec));  for (i = 0; i < 2; i++)    {      xp[0] = 578;      r = mpn_divexact_byfobm1 (wp, xp, (mp_size_t) 1, CNST_LIMB(17),CNST_LIMB(-1)/CNST_LIMB(17));      ASSERT_ALWAYS (wp[0] == 34);      ASSERT_ALWAYS (r == 0);    }  memcpy (&__gmpn_cpuvec, &initial_cpuvec, sizeof (__gmpn_cpuvec));  for (i = 0; i < 2; i++)    {      xp[0] = 287;      r = mpn_divrem_1 (wp, (mp_size_t) 1, xp, (mp_size_t) 1, CNST_LIMB(7));      ASSERT_ALWAYS (wp[1] == 41);      ASSERT_ALWAYS (wp[0] == 0);      ASSERT_ALWAYS (r == 0);    }  memcpy (&__gmpn_cpuvec, &initial_cpuvec, sizeof (__gmpn_cpuvec));  for (i = 0; i < 2; i++)    {      xp[0] = 290;      r = mpn_divrem_euclidean_qr_1 (wp, 0, xp, (mp_size_t) 1, CNST_LIMB(7));      ASSERT_ALWAYS (wp[0] == 41);      ASSERT_ALWAYS (r == 3);    }  memcpy (&__gmpn_cpuvec, &initial_cpuvec, sizeof (__gmpn_cpuvec));  for (i = 0; i < 2; i++)    {      xp[0] = 12;      r = mpn_gcd_1 (xp, (mp_size_t) 1, CNST_LIMB(9));      ASSERT_ALWAYS (r == 3);    }//.........这里部分代码省略.........