C++ CNST_LIMB函数代码示例

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;

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