void
check_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.  */
void
check_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
}

int
mpz_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 void
cfdiv_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_t
DO_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 void
gmp_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 */
    }
}

void
check_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 ();
}

void
mpn_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.  */
void
mpz_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. */
void
check_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);
}

void
gmp_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;

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

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

mp_limb_t
r_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 void
cfdiv_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 void
mpf_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);
}

void
check_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();
        }
    }
}