/* HGCD definitions */

/* Limited by 2 + twice the bitsize of mp_size_t */
#define QSTACK_MAX_QUOTIENTS 82

/* Name mangling */
#define qstack_itch __gmpn_qstack_itch
#define qstack_init __gmpn_qstack_init
#define qstack_reset __gmpn_qstack_reset
#define qstack_rotate __gmpn_qstack_rotate

#define mpn_hgcd2 __gmpn_hgcd2
#define mpn_hgcd2_fix __gmpn_hgcd2_fix
#define mpn_hgcd2_lehmer_step __gmpn_hgcd2_lehmer_step
#define mpn_hgcd_max_recursion __gmpn_hgcd_max_recursion
#define mpn_hgcd_init_itch __gmpn_hgcd_init_itch
#define mpn_hgcd_init __gmpn_hgcd_init
#define mpn_hgcd_lehmer_itch __gmpn_hgcd_lehmer_itch
#define mpn_hgcd_lehmer __gmpn_hgcd_lehmer
#define mpn_hgcd_itch __gmpn_hgcd_itch
#define mpn_hgcd __gmpn_hgcd
#define mpn_hgcd_equal __gmpn_hgcd_equal
#define mpn_hgcd_fix __gmpn_hgcd_fix

struct qstack
{
  /* Throughout the code we represent q = 1 with qsize = 0. */
  mp_size_t size[QSTACK_MAX_QUOTIENTS];
  mp_ptr limb;
  mp_size_t limb_alloc;

  /* Number of quotients to keep when we discard old quotients */
  unsigned nkeep;

  /* Top quotient is of size size[size_next-1], and starts at
     limb+limb_next - size[size_next-1]. We use size_next == 0 for an
     empty stack.*/
  unsigned size_next;
  mp_size_t limb_next;
};

mp_size_t
qstack_itch __GMP_PROTO ((mp_size_t size));

void
qstack_init __GMP_PROTO ((struct qstack *stack,
			  mp_size_t asize,
			  mp_limb_t *limbs, mp_size_t alloc));

void
qstack_reset __GMP_PROTO ((struct qstack *stack,
			   mp_size_t asize));

void
qstack_rotate __GMP_PROTO ((struct qstack *stack,
			    mp_size_t size));

#if WANT_ASSERT
void
__gmpn_qstack_sanity __GMP_PROTO ((struct qstack *stack));
#define ASSERT_QSTACK __gmpn_qstack_sanity
#else
#define ASSERT_QSTACK(stack)
#endif

struct hgcd2_row
{
  /* r = (-)u a + (-)v b */
  mp_limb_t u;
  mp_limb_t v;
};

struct hgcd2
{
  /* Sign of the first row, sign >= 0 implies that u >= 0 and v <= 0,
     sign < 0 implies u <= 0, v >= 0 */
  int sign;
  struct hgcd2_row row[4];
};

int
mpn_hgcd2 __GMP_PROTO ((struct hgcd2 *hgcd,
			mp_limb_t ah, mp_limb_t al,
			mp_limb_t bh, mp_limb_t bl,
			struct qstack *quotients));

mp_size_t
mpn_hgcd2_fix __GMP_PROTO ((mp_ptr rp, mp_size_t ralloc,
			    int sign,
			    mp_limb_t u, mp_srcptr ap, mp_size_t asize,
			    mp_limb_t v, mp_srcptr bp, mp_size_t bsize));

int
mpn_hgcd2_lehmer_step __GMP_PROTO ((struct hgcd2 *hgcd,
				    mp_srcptr ap, mp_size_t asize,
				    mp_srcptr bp, mp_size_t bsize,
				    struct qstack *quotients));

unsigned
mpn_hgcd_max_recursion __GMP_PROTO ((mp_size_t n));

struct hgcd_row
{
  /* [rp, rsize] should always be normalized. */
  mp_ptr rp; mp_size_t rsize;
  mp_ptr uvp[2];
};

struct hgcd
{
  int sign;
  /* Space allocated for the uv entries, for sanity checking */
  mp_size_t alloc;
  /* Size of the largest u,v entry, usually row[3].uvp[1]. This
     element should be normalized. Smaller elements must be zero
     padded, and all unused limbs (i.e. between size and alloc) must
     be zero. */
  mp_size_t size;
  struct hgcd_row row[4];
};

mp_size_t
mpn_hgcd_init_itch __GMP_PROTO ((mp_size_t size));

void
mpn_hgcd_init __GMP_PROTO ((struct hgcd *hgcd,
			    mp_size_t asize,
			    mp_limb_t *limbs));

mp_size_t
mpn_hgcd_lehmer_itch __GMP_PROTO ((mp_size_t asize));

int
mpn_hgcd_lehmer __GMP_PROTO ((struct hgcd *hgcd,
			      mp_srcptr ap, mp_size_t asize,
			      mp_srcptr bp, mp_size_t bsize,
			      struct qstack *quotients,
			      mp_ptr tp, mp_size_t talloc));

mp_size_t
mpn_hgcd_itch __GMP_PROTO ((mp_size_t size));

int
mpn_hgcd __GMP_PROTO ((struct hgcd *hgcd,
		       mp_srcptr ap, mp_size_t asize,
		       mp_srcptr bp, mp_size_t bsize,
		       struct qstack *quotients,
		       mp_ptr tp, mp_size_t talloc));

#if WANT_ASSERT
void
__gmpn_hgcd_sanity __GMP_PROTO ((const struct hgcd *hgcd,
				 mp_srcptr ap, mp_size_t asize,
				 mp_srcptr bp, mp_size_t bsize,
				 unsigned start, unsigned end));
#define ASSERT_HGCD __gmpn_hgcd_sanity
#else
#define ASSERT_HGCD(hgcd, ap, asize, bp, bsize, start, end)
#endif

int
mpn_hgcd_equal __GMP_PROTO ((const struct hgcd *A, const struct hgcd *B));

mp_size_t
mpn_hgcd_fix __GMP_PROTO ((mp_size_t k,
			   mp_ptr rp, mp_size_t ralloc,
			   int sign, mp_size_t uvsize,
			   const struct hgcd_row *s,
			   mp_srcptr ap, mp_srcptr bp,
			   mp_ptr tp, mp_size_t talloc));

#ifndef HGCD_SCHOENHAGE_THRESHOLD
#define HGCD_SCHOENHAGE_THRESHOLD 150
#endif

#if 0
#ifndef GCD_LEHMER_THRESHOLD
#define GCD_LEHMER_THRESHOLD 200
#endif
#endif

#ifndef GCD_SCHOENHAGE_THRESHOLD
#define GCD_SCHOENHAGE_THRESHOLD 1000
#endif

#ifndef GCDEXT_SCHOENHAGE_THRESHOLD
#define GCDEXT_SCHOENHAGE_THRESHOLD 600
#endif