void
mpz_tdiv_q (mpz_ptr quot, mpz_srcptr num, mpz_srcptr den)
{
  mp_size_t ql;
  mp_size_t ns, ds, nl, dl;
  mp_ptr np, dp, qp;
  TMP_DECL;

  ns = SIZ (num);
  ds = SIZ (den);
  nl = ABS (ns);
  dl = ABS (ds);
  ql = nl - dl + 1;

  if (dl == 0)
    DIVIDE_BY_ZERO;

  if (ql <= 0)
    {
      SIZ (quot) = 0;
      return;
    }

  MPZ_REALLOC (quot, ql);

  TMP_MARK;
  qp = PTR (quot);
  np = PTR (num);
  dp = PTR (den);

  /* Copy denominator to temporary space if it overlaps with the quotient.  */
  if (dp == qp)
    {
      mp_ptr tp;
      tp = TMP_ALLOC_LIMBS (dl);
      MPN_COPY (tp, dp, dl);
      dp = tp;
    }
  /* Copy numerator to temporary space if it overlaps with the quotient.  */
  if (np == qp)
    {
      mp_ptr tp;
      tp = TMP_ALLOC_LIMBS (nl + 1);
      MPN_COPY (tp, np, nl);
      /* Overlap dividend and scratch.  */
      mpn_div_q (qp, tp, nl, dp, dl, tp);
    }
  else
    {
      mp_ptr tp;
      tp = TMP_ALLOC_LIMBS (nl + 1);
      mpn_div_q (qp, np, nl, dp, dl, tp);
    }

  ql -=  qp[ql - 1] == 0;

  SIZ (quot) = (ns ^ ds) >= 0 ? ql : -ql;
  TMP_FREE;
}

int
mpz_root (mpz_ptr root, mpz_srcptr u, mpir_ui nth)
{
  mp_ptr rootp, up;
  mp_size_t us, un, rootn, remn;
  TMP_DECL;

  us = SIZ(u);

  /* even roots of negatives provoke an exception */
  if (us < 0 && (nth & 1) == 0)
    SQRT_OF_NEGATIVE;

  /* root extraction interpreted as c^(1/nth) means a zeroth root should
     provoke a divide by zero, do this even if c==0 */
  if (nth == 0)
    DIVIDE_BY_ZERO;

  if (us == 0)
    {
      if (root != NULL)
	SIZ(root) = 0;
      return 1;			/* exact result */
    }

  un = ABS (us);
  rootn = (un - 1) / nth + 1;

  TMP_MARK;

  /* FIXME: Perhaps disallow root == NULL */
  if (root != NULL && u != root)
    rootp = MPZ_REALLOC (root, rootn);
  else
    rootp = TMP_ALLOC_LIMBS (rootn);

  up = PTR(u);

  if (nth == 1)
    {
      MPN_COPY (rootp, up, un);
      remn = 0;
    }
  else
    {
      remn = mpn_rootrem (rootp, NULL, up, un, (mp_limb_t) nth);
    }

  if (root != NULL)
    {
      SIZ(root) = us >= 0 ? rootn : -rootn;
      if (u == root)
	MPN_COPY (up, rootp, rootn);
    }

  TMP_FREE;
  return remn == 0;
}

static void
matrix_copy (struct matrix *R, const struct matrix *M)
{
  R->n = M->n;
  MPN_COPY (R->e00, M->e00, M->n);
  MPN_COPY (R->e01, M->e01, M->n);
  MPN_COPY (R->e10, M->e10, M->n);
  MPN_COPY (R->e11, M->e11, M->n);
}

static void
ref_mpn_mul (mp_ptr wp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
{
  mp_ptr tp;
  mp_size_t tn;
  mp_limb_t cy;

  if (vn < TOOM3_THRESHOLD)
    {
      /* In the mpn_mul_basecase and mpn_kara_mul_n range, use our own
	 mul_basecase.  */
      if (vn != 0)
	mul_basecase (wp, up, un, vp, vn);
      else
	MPN_ZERO (wp, un);
      return;
    }

  if (vn < FFT_THRESHOLD)
    {
      /* In the mpn_toom3_mul_n and mpn_toom4_mul_n range, use mpn_kara_mul_n.  */
      tn = 2 * vn + MPN_KARA_MUL_N_TSIZE (vn);
      tp = __GMP_ALLOCATE_FUNC_LIMBS (tn);
      mpn_kara_mul_n (tp, up, vp, vn, tp + 2 * vn);
    }
  else
    {
      /* Finally, for the largest operands, use mpn_toom3_mul_n.  */
      /* The "- 63 + 255" tweaks the allocation to allow for huge operands.
	 See the definition of this macro in gmp-impl.h to understand this.  */
      tn = 2 * vn + MPN_TOOM3_MUL_N_TSIZE (vn) - 63 + 255;
      tp = __GMP_ALLOCATE_FUNC_LIMBS (tn);
      mpn_toom3_mul_n (tp, up, vp, vn, tp + 2 * vn);
    }

  if (un != vn)
    {
      if (un - vn < vn)
	ref_mpn_mul (wp + vn, vp, vn, up + vn, un - vn);
      else
	ref_mpn_mul (wp + vn, up + vn, un - vn, vp, vn);

      MPN_COPY (wp, tp, vn);
      cy = mpn_add_n (wp + vn, wp + vn, tp + vn, vn);
      mpn_incr_u (wp + 2 * vn, cy);
    }
  else
    {
      MPN_COPY (wp, tp, 2 * vn);
    }

  __GMP_FREE_FUNC_LIMBS (tp, tn);
}

void
mpf_sqrt (mpf_ptr r, mpf_srcptr u)
{
  mp_size_t usize;
  mp_ptr up, tp;
  mp_size_t prec, tsize;
  mp_exp_t uexp, expodd;
  TMP_DECL;

  usize = u->_mp_size;
  if (usize <= 0)
    {
      if (usize < 0)
        SQRT_OF_NEGATIVE;
      r->_mp_size = 0;
      r->_mp_exp = 0;
      return;
    }

  TMP_MARK;

  uexp = u->_mp_exp;
  prec = r->_mp_prec;
  up = u->_mp_d;

  expodd = (uexp & 1);
  tsize = 2 * prec - expodd;
  r->_mp_size = prec;
  r->_mp_exp = (uexp + expodd) / 2;    /* ceil(uexp/2) */

  /* root size is ceil(tsize/2), this will be our desired "prec" limbs */
  ASSERT ((tsize + 1) / 2 == prec);

  tp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB);

  if (usize > tsize)
    {
      up += usize - tsize;
      usize = tsize;
      MPN_COPY (tp, up, tsize);
    }
  else
    {
      MPN_ZERO (tp, tsize - usize);
      MPN_COPY (tp + (tsize - usize), up, usize);
    }

  mpn_sqrtrem (r->_mp_d, NULL, tp, tsize);

  TMP_FREE;
}

/* Compute t = a mod m, a is defined by (ap,an), m is defined by (mp,mn), and
   t is defined by (tp,mn).  */
static void
reduce (mp_ptr tp, mp_srcptr ap, mp_size_t an, mp_srcptr mp, mp_size_t mn, gmp_pi1_t *dinv)
{
  mp_ptr rp, scratch;
  TMP_DECL;
  TMP_MARK;

  rp = TMP_ALLOC_LIMBS (an);
  scratch = TMP_ALLOC_LIMBS (an - mn + 1);
  MPN_COPY (rp, ap, an);
  mod (rp, an, mp, mn, dinv, scratch);
  MPN_COPY (tp, rp, mn);

  TMP_FREE;
}

void
mpi_set_secure( MPI a )
{
    mpi_ptr_t ap, bp;

    if( (a->flags & 1) )
	return;
    a->flags |= 1;
    ap = a->d;
    if( !a->nlimbs ) {
	assert(!ap);
	return;
    }
#ifdef M_DEBUG
    bp = mpi_debug_alloc_limb_space( a->nlimbs, 1, "set_secure" );
#else
    bp = mpi_alloc_limb_space( a->nlimbs, 1 );
#endif
    MPN_COPY( bp, ap, a->nlimbs );
    a->d = bp;
#ifdef M_DEBUG
    mpi_debug_free_limb_space(ap, "set_secure");
#else
    mpi_free_gpg_limb_space(ap);
#endif
}

void
mpz_sqrt (mpz_ptr root, mpz_srcptr op)
{
  mp_size_t op_size, root_size;
  mp_ptr root_ptr, op_ptr;
  mp_ptr free_me = NULL;
  mp_size_t free_me_size;
  TMP_DECL (marker);

  TMP_MARK (marker);
  op_size = op->_mp_size;
  if (op_size <= 0)
    {
      if (op_size < 0)
        SQRT_OF_NEGATIVE;
      SIZ(root) = 0;
      return;
    }

  /* The size of the root is accurate after this simple calculation.  */
  root_size = (op_size + 1) / 2;

  root_ptr = root->_mp_d;
  op_ptr = op->_mp_d;

  if (root->_mp_alloc < root_size)
    {
      if (root_ptr == op_ptr)
	{
	  free_me = root_ptr;
	  free_me_size = root->_mp_alloc;
	}
      else
	(*__gmp_free_func) (root_ptr, root->_mp_alloc * BYTES_PER_MP_LIMB);

      root->_mp_alloc = root_size;
      root_ptr = (mp_ptr) (*__gmp_allocate_func) (root_size * BYTES_PER_MP_LIMB);
      root->_mp_d = root_ptr;
    }
  else
    {
      /* Make OP not overlap with ROOT.  */
      if (root_ptr == op_ptr)
	{
	  /* ROOT and OP are identical.  Allocate temporary space for OP.  */
	  op_ptr = (mp_ptr) TMP_ALLOC (op_size * BYTES_PER_MP_LIMB);
	  /* Copy to the temporary space.  Hack: Avoid temporary variable
	     by using ROOT_PTR.  */
	  MPN_COPY (op_ptr, root_ptr, op_size);
	}
    }

  mpn_sqrtrem (root_ptr, NULL, op_ptr, op_size);

  root->_mp_size = root_size;

  if (free_me != NULL)
    (*__gmp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);
  TMP_FREE (marker);
}

void
mpf_init_set (mpf_ptr r, mpf_srcptr s)
{
    mp_ptr rp, sp;
    mp_size_t ssize, size;
    mp_size_t prec;

    prec = __gmp_default_fp_limb_precision;
    r->_mp_d = __GMP_ALLOCATE_FUNC_LIMBS (prec + 1);
    r->_mp_prec = prec;

    prec++;		/* lie not to lose precision in assignment */
    ssize = s->_mp_size;
    size = ABS (ssize);

    rp = r->_mp_d;
    sp = s->_mp_d;

    if (size > prec)
    {
        sp += size - prec;
        size = prec;
    }

    r->_mp_exp = s->_mp_exp;
    r->_mp_size = ssize >= 0 ? size : -size;

    MPN_COPY (rp, sp, size);
}

int
test_invert (mp_ptr xp, mp_srcptr ap, mp_size_t n)
{
  int res = 1;
  mp_size_t i;
  mp_ptr tp, up;
  mp_limb_t cy;
  TMP_DECL;

  TMP_MARK;
  tp = TMP_ALLOC_LIMBS (2 * n);
  up = TMP_ALLOC_LIMBS (2 * n);

  /* first check X*A < B^(2*n) */
  mpn_mul_n (tp, xp, ap, n);
  cy = mpn_add_n (tp + n, tp + n, ap, n); /* A * msb(X) */
  if (cy != 0)
    return 0;

  /* now check B^(2n) - X*A <= A */
  mpn_com_n (tp, tp, 2 * n);
  mpn_add_1 (tp, tp, 2 * n, 1); /* B^(2n) - X*A */
  MPN_ZERO (up, 2 * n);
  MPN_COPY (up, ap, n);
  res = mpn_cmp (tp, up, 2 * n) <= 0;
  TMP_FREE;
  return res;
}

mp_limb_t	mpn_sumdiff_n(mp_ptr s,mp_ptr d,mp_srcptr x,mp_srcptr y,mp_size_t n)
{mp_limb_t ret;mp_ptr t;

ASSERT(n>0);
ASSERT_MPN(x,n);ASSERT_MPN(y,n);//ASSERT_SPACE(s,n);ASSERT_SPACE(d,n);
ASSERT(MPN_SAME_OR_SEPARATE_P(s,x,n));
ASSERT(MPN_SAME_OR_SEPARATE_P(s,y,n));
ASSERT(MPN_SAME_OR_SEPARATE_P(d,x,n));
ASSERT(MPN_SAME_OR_SEPARATE_P(d,y,n));
ASSERT(!MPN_OVERLAP_P(s,n,d,n));

if( (s==x && d==y)||(s==y && d==x) )
  {t=__GMP_ALLOCATE_FUNC_LIMBS(n);
   ret=mpn_sub_n(t,x,y,n);
   ret+=2*mpn_add_n(s,x,y,n);
   MPN_COPY(d,t,n);
   __GMP_FREE_FUNC_LIMBS(t,n);
   return ret;}
if(s==x || s==y)
  {ret=mpn_sub_n(d,x,y,n);
   ret+=2*mpn_add_n(s,x,y,n);
   return ret;}
ret=2*mpn_add_n(s,x,y,n);
ret+=mpn_sub_n(d,x,y,n);
return ret;}

/* Multiply M by M1 from the right. Since the M1 elements fit in
   GMP_NUMB_BITS - 1 bits, M grows by at most one limb. Needs
   temporary space M->n */
static void
ngcd_matrix_mul_1 (struct ngcd_matrix *M, const struct ngcd_matrix1 *M1)
{
  unsigned row;
  mp_limb_t grow;
  for (row = 0, grow = 0; row < 2; row++)
    {
      mp_limb_t c0, c1;

      /* Compute (u, u') <-- (r00 u + r10 u', r01 u + r11 u') as

	  t   = u
	  u  *= r00
	  u  += r10 * u'
	  u' *= r11
	  u' += r01 * t
      */

      MPN_COPY (M->tp, M->p[row][0], M->n);
      c0 =     mpn_mul_1 (M->p[row][0], M->p[row][0], M->n, M1->u[0][0]);
      c0 += mpn_addmul_1 (M->p[row][0], M->p[row][1], M->n, M1->u[1][0]);
      M->p[row][0][M->n] = c0;

      c1 =     mpn_mul_1 (M->p[row][1], M->p[row][1], M->n, M1->u[1][1]);
      c1 += mpn_addmul_1 (M->p[row][1], M->tp,        M->n, M1->u[0][1]);
      M->p[row][1][M->n] = c1;

      grow |= (c0 | c1);
    }
  M->n += (grow != 0);
  ASSERT (M->n < M->alloc);
}

void mpn_invert_trunc(mp_ptr x_new, mp_size_t m, mp_srcptr xp, mp_size_t n, mp_srcptr ap)
{
  mp_ptr tp;
  mp_limb_t cy;
  TMP_DECL;

  TMP_MARK;
  tp = TMP_ALLOC_LIMBS (2 * m);
  
  MPN_COPY(x_new, xp + n - m, m);
  ap += (n - m);

  mpn_mul_n (tp, x_new, ap, m);
  mpn_add_n (tp + m, tp + m, ap, m); /* A * msb(X) */
  
  /* now check B^(2n) - X*A <= A */
  mpn_not (tp, 2 * m);
  mpn_add_1 (tp, tp, 2 * m, 1); /* B^(2m) - X*A */
  
  while (tp[m] || mpn_cmp (tp, ap, m) > 0)
  {
     mpn_add_1(x_new, x_new, m, 1);
     tp[m] -= mpn_sub_n(tp, tp, ap, m);
  }
  TMP_FREE;
}

void
mpf_neg (mpf_ptr r, mpf_srcptr u)
{
  mp_size_t size;

  size = -u->_mp_size;
  if (r != u)
    {
      mp_size_t prec;
      mp_size_t asize;
      mp_ptr rp, up;

      prec = r->_mp_prec + 1;	/* lie not to lose precision in assignment */
      asize = ABS (size);
      rp = r->_mp_d;
      up = u->_mp_d;

      if (asize > prec)
	{
	  up += asize - prec;
	  asize = prec;
	}

      MPN_COPY (rp, up, asize);
      r->_mp_exp = u->_mp_exp;
      size = size >= 0 ? asize : -asize;
    }
  r->_mp_size = size;
}

void
_gst_mpz_div_2exp (gst_mpz *w, const gst_mpz *u, unsigned cnt)
{
    mp_size_t usize = u->size;
    mp_size_t wsize;
    mp_size_t abs_usize = ABS (usize);
    mp_size_t limb_cnt;

    limb_cnt = cnt / BITS_PER_MP_LIMB;
    wsize = abs_usize - limb_cnt;
    if (wsize <= 0)
        wsize = 0;
    else
    {
        if (w->alloc < wsize)
            gst_mpz_realloc (w, wsize);

        if (cnt % BITS_PER_MP_LIMB)
            mpn_rshift (w->d, u->d + limb_cnt, abs_usize - limb_cnt,
                        cnt % BITS_PER_MP_LIMB);
        else
            MPN_COPY (w->d, u->d + limb_cnt, abs_usize - limb_cnt);

        wsize -= w->d[wsize - 1] == 0;
    }

    w->size = (usize >= 0) ? wsize : -wsize;
}

/*
 * Set f to z, choosing the smallest precision for f
 * so that z = f*(2^BPML)*zs*2^(RetVal)
 */
static int
set_z (mpfr_ptr f, mpz_srcptr z, mp_size_t *zs)
{
  mp_limb_t *p;
  mp_size_t s;
  int c;
  mp_prec_t pf;

  MPFR_ASSERTD (mpz_sgn (z) != 0);

  /* Remove useless ending 0 */
  for (p = PTR (z), s = *zs = ABS (SIZ (z)) ; *p == 0; p++, s--)
    MPFR_ASSERTD (s >= 0);

  /* Get working precision */
  count_leading_zeros (c, p[s-1]);
  pf = s * BITS_PER_MP_LIMB - c;
  if (pf < MPFR_PREC_MIN)
    pf = MPFR_PREC_MIN;
  mpfr_init2 (f, pf);

  /* Copy Mantissa */
  if (MPFR_LIKELY (c))
    mpn_lshift (MPFR_MANT (f), p, s, c);
  else
    MPN_COPY (MPFR_MANT (f), p, s);

  MPFR_SET_SIGN (f, mpz_sgn (z));
  MPFR_SET_EXP (f, 0);

  return -c;
}

static void
randiset_lc (gmp_randstate_ptr dst, gmp_randstate_srcptr src)
{
  gmp_rand_lc_struct *dstp, *srcp;

  srcp = (gmp_rand_lc_struct *) RNG_STATE (src);
  dstp = (*__gmp_allocate_func) (sizeof (gmp_rand_lc_struct));

  RNG_STATE (dst) = (void *) dstp;
  RNG_FNPTR (dst) = (void *) &Linear_Congruential_Generator;

  /* _mp_seed and _mp_a might be unnormalized (high zero limbs), but
     mpz_init_set won't worry about that */
  mpz_init_set (dstp->_mp_seed, srcp->_mp_seed);
  mpz_init_set (dstp->_mp_a,    srcp->_mp_a);

  dstp->_cn = srcp->_cn;

  dstp->_cp[0] = srcp->_cp[0];
  if (LIMBS_PER_ULONG > 1)
    dstp->_cp[1] = srcp->_cp[1];
  if (LIMBS_PER_ULONG > 2)  /* usually there's only 1 or 2 */
    MPN_COPY (dstp->_cp + 2, srcp->_cp + 2, LIMBS_PER_ULONG - 2);

  dstp->_mp_m2exp = srcp->_mp_m2exp;
}

void _tc4_add_unsigned(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 = r1n;
   mp_size_t s2 = r2n;
   
   if (!s2)
   {
      if (!s1) *rn = 0;
      else
      {
         if (rp != r1) MPN_COPY(rp, r1, s1);
		   *rn = r1n;
		}
   } else
   {
      *rn = r1n;
      cy = mpn_add(rp, r1, s1, r2, s2);
      if (cy) 
      {
         rp[s1] = cy;
         if ((*rn) < 0) (*rn)--;
         else (*rn)++;
      }
   } 
}

/* Multiply M by M1 from the right. Since the M1 elements fit in
   GMP_NUMB_BITS - 1 bits, M grows by at most one limb. Needs
   temporary space M->n */
static void
hgcd_matrix_mul_1 (struct hgcd_matrix *M, const struct hgcd_matrix1 *M1,
           mp_ptr tp)
{
  mp_size_t n0, n1;

  /* Could avoid copy by some swapping of pointers. */
  MPN_COPY (tp, M->p[0][0], M->n);
  n0 = mpn_hgcd_mul_matrix1_vector (M1, M->p[0][0], tp, M->p[0][1], M->n);
  MPN_COPY (tp, M->p[1][0], M->n);
  n1 = mpn_hgcd_mul_matrix1_vector (M1, M->p[1][0], tp, M->p[1][1], M->n);

  /* Depends on zero initialization */
  M->n = MAX(n0, n1);
  ASSERT (M->n < M->alloc);
}

void
mpf_init_set (mpf_ptr r, mpf_srcptr s)
{
    mp_ptr rp, sp;
    mp_size_t ssize, size;
    mp_size_t prec;

    prec = __gmp_default_fp_limb_precision;
    r->_mp_d = (mp_ptr) (*__gmp_allocate_func) ((prec + 1) * BYTES_PER_MP_LIMB);
    r->_mp_prec = prec;

    prec++;		/* lie not to lose precision in assignment */
    ssize = s->_mp_size;
    size = ABS (ssize);

    rp = r->_mp_d;
    sp = s->_mp_d;

    if (size > prec)
    {
        sp += size - prec;
        size = prec;
    }

    r->_mp_exp = s->_mp_exp;
    r->_mp_size = ssize >= 0 ? size : -size;

    MPN_COPY (rp, sp, size);
}