void
mpz_combit (mpz_ptr d, unsigned long int bit_index)
{
  mp_size_t dsize = ABSIZ(d);
  mp_ptr dp = LIMBS(d);
  
  mp_size_t limb_index = bit_index / GMP_NUMB_BITS;
  mp_limb_t bit = ((mp_limb_t) 1 << (bit_index % GMP_NUMB_BITS));

  if (limb_index >= dsize)
    {
      MPZ_REALLOC(d, limb_index + 1);
      dp = LIMBS(d);
      
      MPN_ZERO(dp + dsize, limb_index + 1 - dsize);
      dsize = limb_index + 1;
    }
    
  if (SIZ(d) >= 0)
    {
      dp[limb_index] ^= bit;
      MPN_NORMALIZE (dp, dsize);
      SIZ(d) = dsize;
    }
  else
    {
      mp_limb_t x = -dp[limb_index];
      mp_size_t i;

      /* non-zero limb below us means ones-complement */
      for (i = limb_index-1; i >= 0; i--)
        if (dp[i] != 0)
          {
            x--;  /* change twos comp to ones comp */
            break;
          }

      if (x & bit)
	{
          mp_limb_t  c;

	  /* Clearing the bit increases the magitude. We might need a carry. */
	  MPZ_REALLOC(d, dsize + 1);
	  dp = LIMBS(d);

          __GMPN_ADD_1 (c, dp+limb_index, dp+limb_index,
                        dsize - limb_index, bit);
          dp[dsize] = c;
          dsize += c;
        }
      else
	/* Setting the bit decreases the magnitude */
	mpn_sub_1(dp+limb_index, dp+limb_index, dsize + limb_index, bit);

      MPN_NORMALIZE (dp, dsize);
      SIZ(d) = -dsize;
    }
}

/* 
   c is the top bits of the inputs, (fully reduced)
   c & 2 is the top bit of y
   c & 1 is the top bit of z
*/
int
mpn_mulmod_2expp1_basecase (mp_ptr xp, mp_srcptr yp, mp_srcptr zp, 
                                           int c, mpir_ui b, mp_ptr tp)
{
  int cy, cz;
  mp_size_t n, k;

  cy = c & 2;
  cz = c & 1;

  n = BITS_TO_LIMBS (b);
  k = GMP_NUMB_BITS * n - b;

  ASSERT(b > 0);
  ASSERT(n > 0);
  ASSERT_MPN(yp, n);
  ASSERT_MPN(zp, n);
  ASSERT(!MPN_OVERLAP_P (tp, 2 * n, yp, n));
  ASSERT(!MPN_OVERLAP_P (tp, 2 * n, zp, n));
  ASSERT(MPN_SAME_OR_SEPARATE_P (xp, tp, n));
  ASSERT(MPN_SAME_OR_SEPARATE_P (xp, tp + n, n));
  ASSERT(k == 0 || yp[n - 1] >> (GMP_NUMB_BITS - k) == 0);
  ASSERT(k == 0 || zp[n - 1] >> (GMP_NUMB_BITS - k) == 0);

#if WANT_ASSERT
  {
     mp_size_t t = n;

     MPN_NORMALIZE(yp, t);
     ASSERT(cy == 0 || t == 0);
     
     t = n; 
     MPN_NORMALIZE(zp, t);
     ASSERT(cz == 0 || t == 0);
  }
#endif

  if (LIKELY (cy == 0))
    {
      if (LIKELY (cz == 0))
	{
	  c = mpn_mulmod_2expp1_internal (xp, yp, zp, b, tp);
	}
      else
	{
	  c = mpn_neg_n (xp, yp, n);
	  c = mpn_add_1 (xp, xp, n, c);
	  xp[n - 1] &= GMP_NUMB_MASK >> k;
	}
    }
  else
    {
      if (LIKELY (cz == 0))

void
mpn_trace_file (const char *filename, mp_srcptr ptr, mp_size_t size)
{
  FILE   *fp;
  mpz_t  z;

  fp = fopen (filename, "w");
  if (fp == NULL)
    {
      perror ("fopen");
      abort();
    }

  MPN_NORMALIZE (ptr, size);
  PTR(z) = (mp_ptr) ptr;
  SIZ(z) = (int) size;

  mpz_out_str (fp, mp_trace_base, z);
  fprintf (fp, "\n");

  if (ferror (fp) || fclose (fp) != 0)
    {
      printf ("error writing %s\n", filename);
      abort();
    }
}

void _tc4_add(mp_ptr rp, mp_size_t * rn, mp_srcptr r1, mp_size_t r1n, 
                                                 mp_srcptr r2, mp_size_t r2n)
{
   mp_limb_t cy;
   mp_size_t s1 = ABS(r1n);
   mp_size_t s2 = ABS(r2n);
   
   if (!s1)
   {
      *rn = 0;
   } else if (!s2)
   {
      if (rp != r1) MPN_COPY(rp, r1, s1);
		*rn = r1n;
   } else if ((r1n ^ r2n) >= 0)
   {
      *rn = r1n;
      cy = mpn_add(rp, r1, s1, r2, s2);
      if (cy) 
      {
         rp[s1] = cy;
         if ((*rn) < 0) (*rn)--;
         else (*rn)++;
      }
   } else
   {
      mp_size_t ct;
		if (s1 != s2) ct = 1;
		else MPN_CMP(ct, r1, r2, s1); 
		    
      if (!ct) *rn = 0;
      else if (ct > 0) 
      {
         mpn_sub(rp, r1, s1, r2, s2);
         *rn = s1;
         MPN_NORMALIZE(rp, (*rn));
			if (r1n < 0) *rn = -(*rn);
      }
      else
      {
         mpn_sub_n(rp, r2, r1, s1);
         *rn = s1;
         MPN_NORMALIZE(rp, (*rn));
			if (r1n > 0) *rn = -(*rn);
      }
   }
}

void
mpz_limbs_finish (mpz_ptr x, mp_size_t n)
{
  assert (n >= 0);
  MPN_NORMALIZE (PTR(x), n);

  SIZ (x) = n;
}

static int
mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize)
{
  mp_srcptr ap = a->_mp_d;
  mp_size_t asize = a->_mp_size;

  MPN_NORMALIZE (bp, bsize);
  return asize == bsize && mpn_cmp (ap, bp, asize) == 0;
}

void
mpz_set_n (mpz_ptr z, mp_srcptr p, mp_size_t size)
{
  ASSERT (size >= 0);
  MPN_NORMALIZE (p, size);
  MPZ_REALLOC (z, size);
  MPN_COPY (PTR(z), p, size);
  SIZ(z) = size;
}

void
mpz_init_set_n (mpz_ptr z, mp_srcptr p, mp_size_t size)
{
  ASSERT (size >= 0);

  MPN_NORMALIZE (p, size);
  ALLOC(z) = MAX (size, 1);
  PTR(z) = __GMP_ALLOCATE_FUNC_LIMBS (ALLOC(z));
  SIZ(z) = size;
  MPN_COPY (PTR(z), p, size);
}

mpz_srcptr
mpz_roinit_n (mpz_ptr x, mp_srcptr xp, mp_size_t xs)
{
  mp_size_t xn = ABS(xs);
  MPN_NORMALIZE (xp, xn);

  ALLOC (x) = 0;
  SIZ (x) = xs < 0 ? -xn : xn;
  PTR (x) = (mp_ptr) xp;
  return x;
}

/* Needs some ugly casts. */
mpz_srcptr
mpz_roinit_n (mpz_ptr x, const mp_limb_t *xp, mp_size_t xs)
{
  mp_size_t xn = ABS (xs);
  
  MPN_NORMALIZE (xp, xn);

  x->_mp_size = xs < 0 ? -xn : xn;
  x->_mp_alloc = 0;
  x->_mp_d = (mp_limb_t *) xp;
  return x;
}

void
mpz_divexact (mpz_ptr quot, mpz_srcptr num, mpz_srcptr den)
{
  mp_ptr qp;
  mp_size_t qn;
  mp_srcptr np, dp;
  mp_size_t nn, dn;
  TMP_DECL;

#if WANT_ASSERT
  {
    mpz_t  rem;
    mpz_init (rem);
    mpz_tdiv_r (rem, num, den);
    ASSERT (SIZ(rem) == 0);
    mpz_clear (rem);
  }
#endif

  nn = ABSIZ (num);
  dn = ABSIZ (den);

  if (nn < dn)
    {
      /* This special case avoids segfaults below when the function is
	 incorrectly called with |N| < |D|, N != 0.  It also handles the
	 well-defined case N = 0.  */
      SIZ(quot) = 0;
      return;
    }

  qn = nn - dn + 1;

  TMP_MARK;

  if (quot == num || quot == den)
    qp = TMP_ALLOC_LIMBS (qn);
  else
    qp = MPZ_REALLOC (quot, qn);

  np = PTR(num);
  dp = PTR(den);

  mpn_divexact (qp, np, nn, dp, dn);
  MPN_NORMALIZE (qp, qn);

  if (qp != PTR(quot))
    MPN_COPY (MPZ_REALLOC (quot, qn), qp, qn);

  SIZ(quot) = (SIZ(num) ^ SIZ(den)) >= 0 ? qn : -qn;

  TMP_FREE;
}

/* this function is useful in debug mode to print non-normalized residues */
static void
mpresn_print (mpres_t x, mpmod_t n)
{
  mp_size_t m, xn;

  xn = SIZ(x);
  m = ABSIZ(x);
  MPN_NORMALIZE(PTR(x), m);
  SIZ(x) = xn >= 0 ? m : -m;
  gmp_printf ("%Zd\n", x);
  SIZ(x) = xn;
}

void
mpz_urandomb (mpz_ptr rop, gmp_randstate_t rstate, unsigned long int nbits)
{
  mp_ptr rp;
  mp_size_t size;

  size = BITS_TO_LIMBS (nbits);
  rp = MPZ_REALLOC (rop, size);

  _gmp_rand (rp, rstate, nbits);
  MPN_NORMALIZE (rp, size);
  SIZ (rop) = size;
}

/* Print "name=value\n" to stdout for an mpn style ptr,size. */
void
mpn_trace (const char *name, mp_srcptr ptr, mp_size_t size)
{
  mpz_t  z;
  if (ptr == NULL)
    {
      mpz_trace (name, NULL);
      return;
    }
  MPN_NORMALIZE (ptr, size);
  PTR(z) = (mp_ptr) ptr;
  SIZ(z) = size;
  ALLOC(z) = size;
  mpz_trace (name, z);
}

void
mpz_tdiv_r_2exp (mpz_ptr res, mpz_srcptr in, mp_bitcnt_t cnt)
{
  mp_size_t in_size = ABS (in->_mp_size);
  mp_size_t res_size;
  mp_size_t limb_cnt = cnt / GMP_NUMB_BITS;
  mp_srcptr in_ptr = in->_mp_d;

  if (in_size > limb_cnt)
    {
      /* The input operand is (probably) greater than 2**CNT.  */
      mp_limb_t x;

      x = in_ptr[limb_cnt] & (((mp_limb_t) 1 << cnt % GMP_NUMB_BITS) - 1);
      if (x != 0)
	{
	  res_size = limb_cnt + 1;
	  if (res->_mp_alloc < res_size)
	    _mpz_realloc (res, res_size);

	  res->_mp_d[limb_cnt] = x;
	}
      else
	{
	  res_size = limb_cnt;
	  MPN_NORMALIZE (in_ptr, res_size);

	  if (res->_mp_alloc < res_size)
	    _mpz_realloc (res, res_size);

	  limb_cnt = res_size;
	}
    }
  else
    {
      /* The input operand is smaller than 2**CNT.  We perform a no-op,
	 apart from that we might need to copy IN to RES.  */
      res_size = in_size;
      if (res->_mp_alloc < res_size)
	_mpz_realloc (res, res_size);

      limb_cnt = res_size;
    }

  if (res != in)
    MPN_COPY (res->_mp_d, in->_mp_d, limb_cnt);
  res->_mp_size = in->_mp_size >= 0 ? res_size : -res_size;
}

/****************
 * Shift A by N bits to the right
 * FIXME: should use alloc_limb if X and A are same.
 */
void
mpi_rshift( MPI x, MPI a, unsigned n )
{
    mpi_ptr_t xp;
    mpi_size_t xsize;

    xsize = a->nlimbs;
    x->sign = a->sign;
    RESIZE_IF_NEEDED(x, xsize);
    xp = x->d;

    if( xsize ) {
	mpihelp_rshift( xp, a->d, xsize, n);
	MPN_NORMALIZE( xp, xsize);
    }
    x->nlimbs = xsize;
}

/* this function is useful in debug mode to print residues */
static void
mpres_print (mpres_t x, char* name, mpmod_t n)
{
  mp_size_t m, xn;
  mpres_t t;
  mpres_init(t, n);
  mpz_set_ui(t, 1);
  mpres_mul (t, x, t, n);

  xn = SIZ(t);
  m = ABSIZ(t);
  MPN_NORMALIZE(PTR(t), m);
  SIZ(t) = xn >= 0 ? m : -m;
  gmp_printf ("%s=%Zd\n", name, t);
  SIZ(t) = xn;
  mpres_clear (t, n);
}

void
mpz_tdiv_r_2exp (mpz_ptr res, mpz_srcptr in, mp_bitcnt_t cnt)
{
  mp_size_t in_size = ABSIZ (in);
  mp_size_t res_size;
  mp_size_t limb_cnt = cnt / GMP_NUMB_BITS;
  mp_srcptr in_ptr = PTR (in);

  if (in_size > limb_cnt)
    {
      /* The input operand is (probably) greater than 2**CNT.  */
      mp_limb_t x;

      x = in_ptr[limb_cnt] & (((mp_limb_t) 1 << cnt % GMP_NUMB_BITS) - 1);
      if (x != 0)
	{
	  res_size = limb_cnt + 1;
	  MPZ_REALLOC (res, res_size);

	  PTR (res)[limb_cnt] = x;
	}
      else
	{
	  res_size = limb_cnt;
	  MPN_NORMALIZE (in_ptr, res_size);

	  MPZ_REALLOC (res, res_size);

	  limb_cnt = res_size;
	}
    }
  else
    {
      /* The input operand is smaller than 2**CNT.  We perform a no-op,
	 apart from that we might need to copy IN to RES.  */
      res_size = in_size;
      MPZ_REALLOC (res, res_size);

      limb_cnt = res_size;
    }

  if (res != in)
    MPN_COPY (PTR (res), PTR (in), limb_cnt);
  SIZ (res) = SIZ (in) >= 0 ? res_size : -res_size;
}

static void
hgcd_hook (void *p, mp_srcptr gp, mp_size_t gn,
	   mp_srcptr qp, mp_size_t qn, int d)
{
  ASSERT (!gp);
  ASSERT (d >= 0);
  ASSERT (d <= 1);

  MPN_NORMALIZE (qp, qn);
  if (qn > 0)
    {
      struct hgcd_matrix *M = (struct hgcd_matrix *) p;
      /* NOTES: This is a bit ugly. A tp area is passed to
	 gcd_subdiv_step, which stores q at the start of that area. We
	 now use the rest. */
      mp_ptr tp = (mp_ptr) qp + qn;
      mpn_hgcd_matrix_update_q (M, qp, qn, d, tp);
    }
}

void
gmpmee_redc_decode (mpz_ptr r, mpz_srcptr u, gmpmee_redc_t mrt)
{
  mp_ptr rp, mp;
  mp_size_t n;
  mpz_srcptr m;

  gmpmee_redc_mod (r, u, mrt);

  m = MPZ(mrt);
  n = SIZ(m);

  rp = PTR(r);
  mp = PTR(m);

  if (mpn_cmp (rp, mp, n) >= 0)
    mpn_sub_n (rp, rp, mp, n);

  MPN_NORMALIZE (rp, n);
  SIZ(r) = n;
}