textures/base/pack/rangeview_btn.png

/* mini-gmp, a minimalistic implementation of a GNU GMP subset.

   Contributed to the GNU project by Niels Möller

Copyright 1991-1997, 1999-2019 Free Software Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */

/* NOTE: All functions in this file which are not declared in
   mini-gmp.h are internal, and are not intended to be compatible
   neither with GMP nor with future versions of mini-gmp. */

/* Much of the material copied from GMP files, including: gmp-impl.h,
   longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c,
   mpn/generic/lshift.c, mpn/generic/mul_1.c,
   mpn/generic/mul_basecase.c, mpn/generic/rshift.c,
   mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c,
   mpn/generic/submul_1.c. */

#include <assert.h>
#include <ctype.h>
#include <limits.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "mini-gmp.h"

#if !defined(MINI_GMP_DONT_USE_FLOAT_H)
#include <float.h>
#endif


/* Macros */
#define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)

#define GMP_LIMB_MAX ((mp_limb_t) ~ (mp_limb_t) 0)
#define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1))

#define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2))
#define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1)

#define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT)
#define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1))

#define GMP_ABS(x) ((x) >= 0 ? (x) : -(x))
#define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1))

#define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
#define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))

#define GMP_CMP(a,b) (((a) > (b)) - ((a) < (b)))

#if defined(DBL_MANT_DIG) && FLT_RADIX == 2
#define GMP_DBL_MANT_BITS DBL_MANT_DIG
#else
#define GMP_DBL_MANT_BITS (53)
#endif

/* Return non-zero if xp,xsize and yp,ysize overlap.
   If xp+xsize<=yp there's no overlap, or if yp+ysize<=xp there's no
   overlap.  If both these are false, there's an overlap. */
#define GMP_MPN_OVERLAP_P(xp, xsize, yp, ysize)				\
  ((xp) + (xsize) > (yp) && (yp) + (ysize) > (xp))

#define gmp_assert_nocarry(x) do { \
    mp_limb_t __cy = (x);	   \
    assert (__cy == 0);		   \
  } while (0)

#define gmp_clz(count, x) do {						\
    mp_limb_t __clz_x = (x);						\
    unsigned __clz_c = 0;						\
    int LOCAL_SHIFT_BITS = 8;						\
    if (GMP_LIMB_BITS > LOCAL_SHIFT_BITS)				\
      for (;								\
	   (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0;	\
	   __clz_c += 8)						\
	{ __clz_x <<= LOCAL_SHIFT_BITS;	}				\
    for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++)		\
      __clz_x <<= 1;							\
    (count) = __clz_c;							\
  } while (0)

#define gmp_ctz(count, x) do {						\
    mp_limb_t __ctz_x = (x);						\
    unsigned __ctz_c = 0;						\
    gmp_clz (__ctz_c, __ctz_x & - __ctz_x);				\
    (count) = GMP_LIMB_BITS - 1 - __ctz_c;				\
  } while (0)

#define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \
  do {									\
    mp_limb_t __x;							\
    __x = (al) + (bl);							\
    (sh) = (ah) + (bh) + (__x < (al));					\
    (sl) = __x;								\
  } while (0)

#define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \
  do {									\
    mp_limb_t __x;							\
    __x = (al) - (bl);							\
    (sh) = (ah) - (bh) - ((al) < (bl));					\
    (sl) = __x;								\
  } while (0)

#define gmp_umul_ppmm(w1, w0, u, v)					\
  do {									\
    int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS;				\
    if (sizeof(unsigned int) * CHAR_BIT >= 2 * GMP_LIMB_BITS)		\
      {									\
	unsigned int __ww = (unsigned int) (u) * (v);			\
	w0 = (mp_limb_t) __ww;						\
	w1 = (mp_limb_t) (__ww >> LOCAL_GMP_LIMB_BITS);			\
      }									\
    else if (GMP_ULONG_BITS >= 2 * GMP_LIMB_BITS)			\
      {									\
	unsigned long int __ww = (unsigned long int) (u) * (v);		\
	w0 = (mp_limb_t) __ww;						\
	w1 = (mp_limb_t) (__ww >> LOCAL_GMP_LIMB_BITS);			\
      }									\
    else {								\
      mp_limb_t __x0, __x1, __x2, __x3;					\
      unsigned __ul, __vl, __uh, __vh;					\
      mp_limb_t __u = (u), __v = (v);					\
									\
      __ul = __u & GMP_LLIMB_MASK;					\
      __uh = __u >> (GMP_LIMB_BITS / 2);				\
      __vl = __v & GMP_LLIMB_MASK;					\
      __vh = __v >> (GMP_LIMB_BITS / 2);				\
									\
      __x0 = (mp_limb_t) __ul * __vl;					\
      __x1 = (mp_limb_t) __ul * __vh;					\
      __x2 = (mp_limb_t) __uh * __vl;					\
      __x3 = (mp_limb_t) __uh * __vh;					\
									\
      __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */	\
      __x1 += __x2;		/* but this indeed can */		\
      if (__x1 < __x2)		/* did we get it? */			\
	__x3 += GMP_HLIMB_BIT;	/* yes, add it in the proper pos. */	\
									\
      (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2));			\
      (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK);	\
    }									\
  } while (0)

#define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di)			\
  do {									\
    mp_limb_t _qh, _ql, _r, _mask;					\
    gmp_umul_ppmm (_qh, _ql, (nh), (di));				\
    gmp_add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl));		\
    _r = (nl) - _qh * (d);						\
    _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */		\
    _qh += _mask;							\
    _r += _mask & (d);							\
    if (_r >= (d))							\
      {									\
	_r -= (d);							\
	_qh++;								\
      }									\
									\
    (r) = _r;								\
    (q) = _qh;								\
  } while (0)

#define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv)		\
  do {									\
    mp_limb_t _q0, _t1, _t0, _mask;					\
    gmp_umul_ppmm ((q), _q0, (n2), (dinv));				\
    gmp_add_ssaaaa ((q), _q0, (q), _q0, (n2), (n1));			\
									\
    /* Compute the two most significant limbs of n - q'd */		\
    (r1) = (n1) - (d1) * (q);						\
    gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0));		\
    gmp_umul_ppmm (_t1, _t0, (d0), (q));				\
    gmp_sub_ddmmss ((r1), (r0), (r1), (r0), _t1, _t0);			\
    (q)++;								\
									\
    /* Conditionally adjust q and the remainders */			\
    _mask = - (mp_limb_t) ((r1) >= _q0);				\
    (q) += _mask;							\
    gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \
    if ((r1) >= (d1))							\
      {									\
	if ((r1) > (d1) || (r0) >= (d0))				\
	  {								\
	    (q)++;							\
	    gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0));	\
	  }								\
      }									\
  } while (0)

/* Swap macros. */
#define MP_LIMB_T_SWAP(x, y)						\
  do {									\
    mp_limb_t __mp_limb_t_swap__tmp = (x);				\
    (x) = (y);								\
    (y) = __mp_limb_t_swap__tmp;					\
  } while (0)
#define MP_SIZE_T_SWAP(x, y)						\
  do {									\
    mp_size_t __mp_size_t_swap__tmp = (x);				\
    (x) = (y);								\
    (y) = __mp_size_t_swap__tmp;					\
  } while (0)
#define MP_BITCNT_T_SWAP(x,y)			\
  do {						\
    mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x);	\
    (x) = (y);					\
    (y) = __mp_bitcnt_t_swap__tmp;		\
  } while (0)
#define MP_PTR_SWAP(x, y)						\
  do {									\
    mp_ptr __mp_ptr_swap__tmp = (x);					\
    (x) = (y);								\
    (y) = __mp_ptr_swap__tmp;						\
  } while (0)
#define MP_SRCPTR_SWAP(x, y)						\
  do {									\
    mp_srcptr __mp_srcptr_swap__tmp = (x);				\
    (x) = (y);								\
    (y) = __mp_srcptr_swap__tmp;					\
  } while (0)

#define MPN_PTR_SWAP(xp,xs, yp,ys)					\
  do {									\
    MP_PTR_SWAP (xp, yp);						\
    MP_SIZE_T_SWAP (xs, ys);						\
  } while(0)
#define MPN_SRCPTR_SWAP(xp,xs, yp,ys)					\
  do {									\
    MP_SRCPTR_SWAP (xp, yp);						\
    MP_SIZE_T_SWAP (xs, ys);						\
  } while(0)

#define MPZ_PTR_SWAP(x, y)						\
  do {									\
    mpz_ptr __mpz_ptr_swap__tmp = (x);					\
    (x) = (y);								\
    (y) = __mpz_ptr_swap__tmp;						\
  } while (0)
#define MPZ_SRCPTR_SWAP(x, y)						\
  do {									\
    mpz_srcptr __mpz_srcptr_swap__tmp = (x);				\
    (x) = (y);								\
    (y) = __mpz_srcptr_swap__tmp;					\
  } while (0)

const int mp_bits_per_limb = GMP_LIMB_BITS;


/* Memory allocation and other helper functions. */
static void
gmp_die (const char *msg)
{
  fprintf (stderr, "%s\n", msg);
  abort();
}

static void *
gmp_default_alloc (size_t size)
{
  void *p;

  assert (size > 0);

  p = malloc (size);
  if (!p)
    gmp_die("gmp_default_alloc: Virtual memory exhausted.");

  return p;
}

static void *
gmp_default_realloc (void *old, size_t unused_old_size, size_t new_size)
{
  void * p;

  p = realloc (old, new_size);

  if (!p)
    gmp_die("gmp_default_realloc: Virtual memory exhausted.");

  return p;
}

static void
gmp_default_free (void *p, size_t unused_size)
{
  free (p);
}

static void * (*gmp_allocate_func) (size_t) = gmp_default_alloc;
static void * (*gmp_reallocate_func) (void *, size_t, size_t) = gmp_default_realloc;
static void (*gmp_free_func) (void *, size_t) = gmp_default_free;

void
mp_get_memory_functions (void *(**alloc_func) (size_t),
			 void *(**realloc_func) (void *, size_t, size_t),
			 void (**free_func) (void *, size_t))
{
  if (alloc_func)
    *alloc_func = gmp_allocate_func;

  if (realloc_func)
    *realloc_func = gmp_reallocate_func;

  if (free_func)
    *free_func = gmp_free_func;
}

void
mp_set_memory_functions (void *(*alloc_func) (size_t),
			 void *(*realloc_func) (void *, size_t, size_t),
			 void (*free_func) (void *, size_t))
{
  if (!alloc_func)
    alloc_func = gmp_default_alloc;
  if (!realloc_func)
    realloc_func = gmp_default_realloc;
  if (!free_func)
    free_func = gmp_default_free;

  gmp_allocate_func = alloc_func;
  gmp_reallocate_func = realloc_func;
  gmp_free_func = free_func;
}

#define gmp_xalloc(size) ((*gmp_allocate_func)((size)))
#define gmp_free(p) ((*gmp_free_func) ((p), 0))

static mp_ptr
gmp_xalloc_limbs (mp_size_t size)
{
  return (mp_ptr) gmp_xalloc (size * sizeof (mp_limb_t));
}

static mp_ptr
gmp_xrealloc_limbs (mp_ptr old, mp_size_t size)
{
  assert (size > 0);
  return (mp_ptr) (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t));
}


/* MPN interface */

void
mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n)
{
  mp_size_t i;
  for (i = 0; i < n; i++)
    d[i] = s[i];
}

void
mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
{
  while (--n >= 0)
    d[n] = s[n];
}

int
mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
{
  while (--n >= 0)
    {
      if (ap[n] != bp[n])
	return ap[n] > bp[n] ? 1 : -1;
    }
  return 0;
}

static int
mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
{
  if (an != bn)
    return an < bn ? -1 : 1;
  else
    return mpn_cmp (ap, bp, an);
}

static mp_size_t
mpn_normalized_size (mp_srcptr xp, mp_size_t n)
{
  while (n > 0 && xp[n-1] == 0)
    --n;
  return n;
}

int
mpn_zero_p(mp_srcptr rp, mp_size_t n)
{
  return mpn_normalized_size (rp, n) == 0;
}

void
mpn_zero (mp_ptr rp, mp_size_t n)
{
  while (--n >= 0)
    rp[n] = 0;
}

mp_limb_t
mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
{
  mp_size_t i;

  assert (n > 0);
  i = 0;
  do
    {
      mp_limb_t r = ap[i] + b;
      /* Carry out */
      b = (r < b);
      rp[i] = r;
    }
  while (++i < n);

  return b;
}

mp_limb_t
mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
{
  mp_size_t i;
  mp_limb_t cy;

  for (i = 0, cy = 0; i < n; i++)
    {
      mp_limb_t a, b, r;
      a = ap[i]; b = bp[i];
      r = a + cy;
      cy = (r < cy);
      r += b;
      cy += (r < b);
      rp[i] = r;
    }
  return cy;
}

mp_limb_t
mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
{
  mp_limb_t cy;

  assert (an >= bn);

  cy = mpn_add_n (rp, ap, bp, bn);
  if (an > bn)
    cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy);
  return cy;
}

mp_limb_t
mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
{
  mp_size_t i;

  assert (n > 0);

  i = 0;
  do
    {
      mp_limb_t a = ap[i];
      /* Carry out */
      mp_limb_t cy = a < b;
      rp[i] = a - b;
      b = cy;
    }
  while (++i < n);

  return b;
}

mp_limb_t
mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
{
  mp_size_t i;
  mp_limb_t cy;

  for (i = 0, cy = 0; i < n; i++)
    {
      mp_limb_t a, b;
      a = ap[i]; b = bp[i];
      b += cy;
      cy = (b < cy);
      cy += (a < b);
      rp[i] = a - b;
    }
  return cy;
}

mp_limb_t
mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
{
  mp_limb_t cy;

  assert (an >= bn);

  cy = mpn_sub_n (rp, ap, bp, bn);
  if (an > bn)
    cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy);
  return cy;
}

mp_limb_t
mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
{
  mp_limb_t ul, cl, hpl, lpl;

  assert (n >= 1);

  cl = 0;
  do
    {
      ul = *up++;
      gmp_umul_ppmm (hpl, lpl, ul, vl);

      lpl += cl;
      cl = (lpl < cl) + hpl;

      *rp++ = lpl;
    }
  while (--n != 0);

  return cl;
}

mp_limb_t
mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
{
  mp_limb_t ul, cl, hpl, lpl, rl;

  assert (n >= 1);

  cl = 0;
  do
    {
      ul = *up++;
      gmp_umul_ppmm (hpl, lpl, ul, vl);

      lpl += cl;
      cl = (lpl < cl) + hpl;

      rl = *rp;
      lpl = rl + lpl;
      cl += lpl < rl;
      *rp++ = lpl;
    }
  while (--n != 0);

  return cl;
}

mp_limb_t
mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
{
  mp_limb_t ul, cl, hpl, lpl, rl;

  assert (n >= 1);

  cl = 0;
  do
    {
      ul = *up++;
      gmp_umul_ppmm (hpl, lpl, ul, vl);

      lpl += cl;
      cl = (lpl < cl) + hpl;

      rl = *rp;
      lpl = rl - lpl;
      cl += lpl > rl;
      *rp++ = lpl;
    }
  while (--n != 0);

  return cl;
}

mp_limb_t
mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
{
  assert (un >= vn);
  assert (vn >= 1);
  assert (!GMP_MPN_OVERLAP_P(rp, un + vn, up, un));
  assert (!GMP_MPN_OVERLAP_P(rp, un + vn, vp, vn));

  /* We first multiply by the low order limb. This result can be
     stored, not added, to rp. We also avoid a loop for zeroing this
     way. */

  rp[un] = mpn_mul_1 (rp, up, un, vp[0]);

  /* Now accumulate the product of up[] and the next higher limb from
     vp[]. */

  while (--vn >= 1)
    {
      rp += 1, vp += 1;
      rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
    }
  return rp[un];
}

void
mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
{
  mpn_mul (rp, ap, n, bp, n);
}

void
mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n)
{
  mpn_mul (rp, ap, n, ap, n);
}

mp_limb_t
mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
{
  mp_limb_t high_limb, low_limb;
  unsigned int tnc;
  mp_limb_t retval;

  assert (n >= 1);
  assert (cnt >= 1);
  assert (cnt < GMP_LIMB_BITS);

  up += n;
  rp += n;

  tnc = GMP_LIMB_BITS - cnt;
  low_limb = *--up;
  retval = low_limb >> tnc;
  high_limb = (low_limb << cnt);

  while (--n != 0)
    {
      low_limb = *--up;
      *--rp = high_limb | (low_limb >> tnc);
      high_limb = (low_limb << cnt);
    }
  *--rp = high_limb;

  return retval;
}

mp_limb_t
mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
{
  mp_limb_t high_limb, low_limb;
  unsigned int tnc;
  mp_limb_t retval;

  assert (n >= 1);
  assert (cnt >= 1);
  assert (cnt < GMP_LIMB_BITS);

  tnc = GMP_LIMB_BITS - cnt;
  high_limb = *up++;
  retval = (high_limb << tnc);
  low_limb = high_limb >> cnt;

  while (--n != 0)
    {
      high_limb = *up++;
      *rp++ = low_limb | (high_limb << tnc);
      low_limb = high_limb >> cnt;
    }
  *rp = low_limb;

  return retval;
}

static mp_bitcnt_t
mpn_common_scan (mp_limb_t limb, mp_size_t i, mp_srcptr up, mp_size_t un,
		 mp_limb_t ux)
{
  unsigned cnt;

  assert (ux == 0 || ux == GMP_LIMB_MAX);
  assert (0 <= i && i <= un );

  while (limb == 0)
    {
      i++;
      if (i == un)
	return (ux == 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
      limb = ux ^ up[i];
    }
  gmp_ctz (cnt, limb);
  return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
}

mp_bitcnt_t
mpn_scan1 (mp_srcptr ptr, mp_bitcnt_t bit)
{
  mp_size_t i;
  i = bit / GMP_LIMB_BITS;

  return mpn_common_scan ( ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
			  i, ptr, i, 0);
}

mp_bitcnt_t
mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit)
{
  mp_size_t i;
  i = bit / GMP_LIMB_BITS;

  return mpn_common_scan (~ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
			  i, ptr, i, GMP_LIMB_MAX);
}

void
mpn_com (mp_ptr rp, mp_srcptr up, mp_size_t n)
{
  while (--n >= 0)
    *rp++ = ~ *up++;
}

mp_limb_t
mpn_neg (mp_ptr rp, mp_srcptr up, mp_size_t n)
{
  while (*up == 0)
    {
      *rp = 0;
      if (!--n)
	return 0;
      ++up; ++rp;
    }
  *rp = - *up;
  mpn_com (++rp, ++up, --n);
  return 1;
}


/* MPN division interface. */

/* The 3/2 inverse is defined as

     m = floor( (B^3-1) / (B u1 + u0)) - B
*/
mp_limb_t
mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
{
  mp_limb_t r, m;

  {
    mp_limb_t p, ql;
    unsigned ul, uh, qh;

    /* For notation, let b denote the half-limb base, so that B = b^2.
       Split u1 = b uh + ul. */
    ul = u1 & GMP_LLIMB_MASK;
    uh = u1 >> (GMP_LIMB_BITS / 2);

    /* Approximation of the high half of quotient. Differs from the 2/1
       inverse of the half limb uh, since we have already subtracted
       u0. */
    qh = (u1 ^ GMP_LIMB_MAX) / uh;

    /* Adjust to get a half-limb 3/2 inverse, i.e., we want

       qh' = floor( (b^3 - 1) / u) - b = floor ((b^3 - b u - 1) / u
	   = floor( (b (~u) + b-1) / u),

       and the remainder

       r = b (~u) + b-1 - qh (b uh + ul)
       = b (~u - qh uh) + b-1 - qh ul

       Subtraction of qh ul may underflow, which implies adjustments.
       But by normalization, 2 u >= B > qh ul, so we need to adjust by
       at most 2.
    */

    r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK;

    p = (mp_limb_t) qh * ul;
    /* Adjustment steps taken from udiv_qrnnd_c */
    if (r < p)
      {
	qh--;
	r += u1;
	if (r >= u1) /* i.e. we didn't get carry when adding to r */
	  if (r < p)
	    {
	      qh--;
	      r += u1;
	    }
      }
    r -= p;

    /* Low half of the quotient is

       ql = floor ( (b r + b-1) / u1).

       This is a 3/2 division (on half-limbs), for which qh is a
       suitable inverse. */

    p = (r >> (GMP_LIMB_BITS / 2)) * qh + r;
    /* Unlike full-limb 3/2, we can add 1 without overflow. For this to
       work, it is essential that ql is a full mp_limb_t. */
    ql = (p >> (GMP_LIMB_BITS / 2)) + 1;

    /* By the 3/2 trick, we don't need the high half limb. */
    r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1;

    if (r >= (GMP_LIMB_MAX & (p << (GMP_LIMB_BITS / 2))))
      {
	ql--;
	r += u1;
      }
    m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql;
    if (r >= u1)
      {
	m++;
	r -= u1;
      }
  }

  /* Now m is the 2/1 inverse of u1. If u0 > 0, adjust it to become a
     3/2 inverse. */
  if (u0 > 0)
    {
      mp_limb_t th, tl;
      r = ~r;
      r += u0;
      if (r < u0)
	{
	  m--;
	  if (r >= u1)
	    {
	      m--;
	      r -= u1;
	    }
	  r -= u1;
	}
      gmp_umul_ppmm (th, tl, u0, m);
      r += th;
      if (r < th)
	{
	  m--;
	  m -= ((r > u1) | ((r == u1) & (tl > u0)));
	}
    }

  return m;
}

struct gmp_div_inverse
{
  /* Normalization shift count. */
  unsigned shift;
  /* Normalized divisor (d0 unused for mpn_div_qr_1) */
  mp_limb_t d1, d0;
  /* Inverse, for 2/1 or 3/2. */
  mp_limb_t di;
};

static void
mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d)
{
  unsigned shift;

  assert (d > 0);
  gmp_clz (shift, d);
  inv->shift = shift;
  inv->d1 = d << shift;
  inv->di = mpn_invert_limb (inv->d1);
}

static void
mpn_div_qr_2_invert (struct gmp_div_inverse *inv,
		     mp_limb_t d1, mp_limb_t d0)
{
  unsigned shift;

  assert (d1 > 0);
  gmp_clz (shift, d1);
  inv->shift = shift;
  if (shift > 0)
    {
      d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
      d0 <<= shift;
    }
  inv->d1 = d1;
  inv->d0 = d0;
  inv->di = mpn_invert_3by2 (d1, d0);
}

static void
mpn_div_qr_invert (struct gmp_div_inverse *inv,
		   mp_srcptr dp, mp_size_t dn)
{
  assert (dn > 0);

  if (dn == 1)
    mpn_div_qr_1_invert (inv, dp[0]);
  else if (dn == 2)
    mpn_div_qr_2_invert (inv, dp[1], dp[0]);
  else
    {
      unsigned shift;
      mp_limb_t d1, d0;

      d1 = dp[dn-1];
      d0 = dp[dn-2];
      assert (d1 > 0);
      gmp_clz (shift, d1);
      inv->shift = shift;
      if (shift > 0)
	{
	  d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
	  d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift));
	}
      inv->d1 = d1;
      inv->d0 = d0;
      inv->di = mpn_invert_3by2 (d1, d0);
    }
}

/* Not matching current public gmp interface, rather corresponding to
   the sbpi1_div_* functions. */
static mp_limb_t
mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
		     const struct gmp_div_inverse *inv)
{
  mp_limb_t d, di;
  mp_limb_t r;
  mp_ptr tp = NULL;

  if (inv->shift > 0)
    {
      /* Shift, reusing qp area if possible. In-place shift if qp == np. */
      tp = qp ? qp : gmp_xalloc_limbs (nn);
      r = mpn_lshift (tp, np, nn, inv->shift);
      np = tp;
    }
  else
    r = 0;

  d = inv->d1;
  di = inv->di;
  while (--nn >= 0)
    {
      mp_limb_t q;

      gmp_udiv_qrnnd_preinv (q, r, r, np[nn], d, di);
      if (qp)
	qp[nn] = q;
    }
  if ((inv->shift > 0) && (tp != qp))
    gmp_free (tp);

  return r >> inv->shift;
}

static void
mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
		     const struct gmp_div_inverse *inv)
{
  unsigned shift;
  mp_size_t i;
  mp_limb_t d1, d0, di, r1, r0;

  assert (nn >= 2);
  shift = inv->shift;
  d1 = inv->d1;
  d0 = inv->d0;
  di = inv->di;

  if (shift > 0)
    r1 = mpn_lshift (np, np, nn, shift);
  else
    r1 = 0;

  r0 = np[nn - 1];

  i = nn - 2;
  do
    {
      mp_limb_t n0, q;
      n0 = np[i];
      gmp_udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di);

      if (qp)
	qp[i] = q;
    }
  while (--i >= 0);

  if (shift > 0)
    {
      assert ((r0 & (GMP_LIMB_MAX >> (GMP_LIMB_BITS - shift))) == 0);
      r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift));
      r1 >>= shift;
    }

  np[1] = r1;
  np[0] = r0;
}

static void
mpn_div_qr_pi1 (mp_ptr qp,
		mp_ptr np, mp_size_t nn, mp_limb_t n1,
		mp_srcptr dp, mp_size_t dn,
		mp_limb_t dinv)
{
  mp_size_t i;

  mp_limb_t d1, d0;
  mp_limb_t cy, cy1;
  mp_limb_t q;

  assert (dn > 2);
  assert (nn >= dn);

  d1 = dp[dn - 1];
  d0 = dp[dn - 2];

  assert ((d1 & GMP_LIMB_HIGHBIT) != 0);
  /* Iteration variable is the index of the q limb.
   *
   * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]>
   * by            <d1,          d0,        dp[dn-3],  ..., dp[0] >
   */

  i = nn - dn;
  do
    {
      mp_limb_t n0 = np[dn-1+i];

      if (n1 == d1 && n0 == d0)
	{
	  q = GMP_LIMB_MAX;
	  mpn_submul_1 (np+i, dp, dn, q);
	  n1 = np[dn-1+i];	/* update n1, last loop's value will now be invalid */
	}
      else
	{
	  gmp_udiv_qr_3by2 (q, n1, n0, n1, n0, np[dn-2+i], d1, d0, dinv);

	  cy = mpn_submul_1 (np + i, dp, dn-2, q);

	  cy1 = n0 < cy;
	  n0 = n0 - cy;
	  cy = n1 < cy1;
	  n1 = n1 - cy1;
	  np[dn-2+i] = n0;

	  if (cy != 0)
	    {
	      n1 += d1 + mpn_add_n (np + i, np + i, dp, dn - 1);
	      q--;
	    }
	}

      if (qp)
	qp[i] = q;
    }
  while (--i >= 0);

  np[dn - 1] = n1;
}

static void
mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
		   mp_srcptr dp, mp_size_t dn,
		   const struct gmp_div_inverse *inv)
{
  assert (dn > 0);
  assert (nn >= dn);

  if (dn == 1)
    np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv);
  else if (dn == 2)
    mpn_div_qr_2_preinv (qp, np, nn, inv);
  else
    {
      mp_limb_t nh;
      unsigned shift;

      assert (inv->d1 == dp[dn-1]);
      assert (inv->d0 == dp[dn-2]);
      assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0);

      shift = inv->shift;
      if (shift > 0)
	nh = mpn_lshift (np, np, nn, shift);
      else
	nh = 0;

      mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);

      if (shift > 0)
	gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
    }
}

static void
mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
{
  struct gmp_div_inverse inv;
  mp_ptr tp = NULL;

  assert (dn > 0);
  assert (nn >= dn);

  mpn_div_qr_invert (&inv, dp, dn);
  if (dn > 2 && inv.shift > 0)
    {
      tp = gmp_xalloc_limbs (dn);
      gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift));
      dp = tp;
    }
  mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv);
  if (tp)
    gmp_free (tp);
}


/* MPN base conversion. */
static unsigned
mpn_base_power_of_two_p (unsigned b)
{
  switch (b)
    {
    case 2: return 1;
    case 4: return 2;
    case 8: return 3;
    case 16: return 4;
    case 32: return 5;
    case 64: return 6;
    case 128: return 7;
    case 256: return 8;
    default: return 0;
    }
}

struct mpn_base_info
{
  /* bb is the largest power of the base which fits in one limb, and
     exp is the corresponding exponent. */
  unsigned exp;
  mp_limb_t bb;
};

static void
mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b)
{
  mp_limb_t m;
  mp_limb_t p;
  unsigned exp;

  m = GMP_LIMB_MAX / b;
  for (exp = 1, p = b; p <= m; exp++)
    p *= b;

  info->exp = exp;
  info->bb = p;
}

static mp_bitcnt_t
mpn_limb_size_in_base_2 (mp_limb_t u)
{
  unsigned shift;

  assert (u > 0);
  gmp_clz (shift, u);
  return GMP_LIMB_BITS - shift;
}

static size_t
mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un)
{
  unsigned char mask;
  size_t sn, j;
  mp_size_t i;
  unsigned shift;

  sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1])
	+ bits - 1) / bits;

  mask = (1U << bits) - 1;

  for (i = 0, j = sn, shift = 0; j-- > 0;)
    {
      unsigned char digit = up[i] >> shift;

      shift += bits;

      if (shift >= GMP_LIMB_BITS && ++i < un)
	{
	  shift -= GMP_LIMB_BITS;
	  digit |= up[i] << (bits - shift);
	}
      sp[j] = digit & mask;
    }
  return sn;
}

/* We generate digits from the least significant end, and reverse at
   the end. */
static size_t
mpn_limb_get_str (unsigned char *sp, mp_limb_t w,
		  const struct gmp_div_inverse *binv)
{
  mp_size_t i;
  for (i = 0; w > 0; i++)
    {
      mp_limb_t h, l, r;

      h = w >> (GMP_LIMB_BITS - binv->shift);
      l = w << binv->shift;

      gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di);
      assert ((r & (GMP_LIMB_MAX >> (GMP_LIMB_BITS - binv->shift))) == 0);
      r >>= binv->shift;

      sp[i] = r;
    }
  return i;
}

static size_t
mpn_get_str_other (unsigned char *sp,
		   int base, const struct mpn_base_info *info,
		   mp_ptr up, mp_size_t un)
{
  struct gmp_div_inverse binv;
  size_t sn;
  size_t i;

  mpn_div_qr_1_invert (&binv, base);

  sn = 0;

  if (un > 1)
    {
      struct gmp_div_inverse bbinv;
      mpn_div_qr_1_invert (&bbinv, info->bb);

      do
	{
	  mp_limb_t w;
	  size_t done;
	  w = mpn_div_qr_1_preinv (up, up, un, &bbinv);
	  un -= (up[un-1] == 0);
	  done = mpn_limb_get_str (sp + sn, w, &binv);

	  for (sn += done; done < info->exp; done++)
	    sp[sn++] = 0;
	}
      while (un > 1);
    }
  sn += mpn_limb_get_str (sp + sn, up[0], &binv);

  /* Reverse order */
  for (i = 0; 2*i + 1 < sn; i++)
    {
      unsigned char t = sp[i];
      sp[i] = sp[sn - i - 1];
      sp[sn - i - 1] = t;
    }

  return sn;
}

size_t
mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un)
{
  unsigned bits;

  assert (un > 0);
  assert (up[un-1] > 0);

  bits = mpn_base_power_of_two_p (base);
  if (bits)
    return mpn_get_str_bits (sp, bits, up, un);
  else
    {
      struct mpn_base_info info;

      mpn_get_base_info (&info, base);
      return mpn_get_str_other (sp, base, &info, up, un);
    }
}

static mp_size_t
mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn,
		  unsigned bits)
{
  mp_size_t rn;
  size_t j;
  unsigned shift;

  for (j = sn, rn = 0, shift = 0; j-- > 0; )
    {
      if (shift == 0)
	{
	  rp[rn++] = sp[j];
	  shift += bits;
	}
      else
	{
	  rp[rn-1] |= (mp_limb_t) sp[j] << shift;
	  shift += bits;
	  if (shift >= GMP_LIMB_BITS)
	    {
	      shift -= GMP_LIMB_BITS;
	      if (shift > 0)
		rp[rn++] = (mp_limb_t) sp[j] >> (bits - shift);
	    }
	}
    }
  rn = mpn_normalized_size (rp, rn);
  return rn;
}

/* Result is usually normalized, except for all-zero input, in which
   case a single zero limb is written at *RP, and 1 is returned. */
static mp_size_t
mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
		   mp_limb_t b, const struct mpn_base_info *info)
{
  mp_size_t rn;
  mp_limb_t w;
  unsigned k;
  size_t j;

  assert (sn > 0);

  k = 1 + (sn - 1) % info->exp;

  j = 0;
  w = sp[j++];
  while (--k != 0)
    w = w * b + sp[j++];

  rp[0] = w;

  for (rn = 1; j < sn;)
    {
      mp_limb_t cy;

      w = sp[j++];
      for (k = 1; k < info->exp; k++)
	w = w * b + sp[j++];

      cy = mpn_mul_1 (rp, rp, rn, info->bb);
      cy += mpn_add_1 (rp, rp, rn, w);
      if (cy > 0)
	rp[rn++] = cy;
    }
  assert (j == sn);

  return rn;
}

mp_size_t
mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base)
{
  unsigned bits;

  if (sn == 0)
    return 0;

  bits = mpn_base_power_of_two_p (base);
  if (bits)
    return mpn_set_str_bits (rp, sp, sn, bits);
  else
    {
      struct mpn_base_info info;

      mpn_get_base_info (&info, base);
      return mpn_set_str_other (rp, sp, sn, base, &info);
    }
}


/* MPZ interface */
void
mpz_init (mpz_t r)
{
  static const mp_limb_t dummy_limb = GMP_LIMB_MAX & 0xc1a0;

  r->_mp_alloc = 0;
  r->_mp_size = 0;
  r->_mp_d = (mp_ptr) &dummy_limb;
}

/* The utility of this function is a bit limited, since many functions
   assigns the result variable using mpz_swap. */
void
mpz_init2 (mpz_t r, mp_bitcnt_t bits)
{
  mp_size_t rn;

  bits -= (bits != 0);		/* Round down, except if 0 */
  rn = 1 + bits / GMP_LIMB_BITS;

  r->_mp_alloc = rn;
  r->_mp_size = 0;
  r->_mp_d = gmp_xalloc_limbs (rn);
}

void
mpz_clear (mpz_t r)
{
  if (r->_mp_alloc)
    gmp_free (r->_mp_d);
}

static mp_ptr
mpz_realloc (mpz_t r, mp_size_t size)
{
  size = GMP_MAX (size, 1);

  if (r->_mp_alloc)
    r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size);
  else
    r->_mp_d = gmp_xalloc_limbs (size);
  r->_mp_alloc = size;

  if (GMP_ABS (r->_mp_size) > size)
    r->_mp_size = 0;

  return r->_mp_d;
}

/* Realloc for an mpz_t WHAT if it has less than NEEDED limbs.  */
#define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc			\
			  ? mpz_realloc(z,n)			\
			  : (z)->_mp_d)

/* MPZ assignment and basic conversions. */
void
mpz_set_si (mpz_t r, signed long int x)
{
  if (x >= 0)
    mpz_set_ui (r, x);
  else /* (x < 0) */
    if (GMP_LIMB_BITS < GMP_ULONG_BITS)
      {
	mpz_set_ui (r, GMP_NEG_CAST (unsigned long int, x));
	mpz_neg (r, r);
      }
  else
    {
      r->_mp_size = -1;
      MPZ_REALLOC (r, 1)[0] = GMP_NEG_CAST (unsigned long int, x);
    }
}

void
mpz_set_ui (mpz_t r, unsigned long int x)
{
  if (x > 0)
    {
      r->_mp_size = 1;
      MPZ_REALLOC (r, 1)[0] = x;
      if (GMP_LIMB_BITS < GMP_ULONG_BITS)
	{
	  int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS;
	  while (x >>= LOCAL_GMP_LIMB_BITS)
	    {
	      ++ r->_mp_size;
	      MPZ_REALLOC (r, r->_mp_size)[r->_mp_size - 1] = x;
	    }
	}
    }
  else
    r->_mp_size = 0;
}

void
mpz_set (mpz_t r, const mpz_t x)
{
  /* Allow the NOP r == x */
  if (r != x)
    {
      mp_size_t n;
      mp_ptr rp;

      n = GMP_ABS (x->_mp_size);
      rp = MPZ_REALLOC (r, n);

      mpn_copyi (rp, x->_mp_d, n);
      r->_mp_size = x->_mp_size;
    }
}

void
mpz_init_set_si (mpz_t r, signed long int x)
{
  mpz_init (r);
  mpz_set_si (r, x);
}

void
mpz_init_set_ui (mpz_t r, unsigned long int x)
{
  mpz_init (r);
  mpz_set_ui (r, x);
}

void
mpz_init_set (mpz_t r, const mpz_t x)
{
  mpz_init (r);
  mpz_set (r, x);
}

int
mpz_fits_slong_p (const mpz_t u)
{
  return (LONG_MAX + LONG_MIN == 0 || mpz_cmp_ui (u, LONG_MAX) <= 0) &&
    mpz_cmpabs_ui (u, GMP_NEG_CAST (unsigned long int, LONG_MIN)) <= 0;
}

static int
mpn_absfits_ulong_p (mp_srcptr up, mp_size_t un)
{
  int ulongsize = GMP_ULONG_BITS / GMP_LIMB_BITS;
  mp_limb_t ulongrem = 0;

  if (GMP_ULONG_BITS % GMP_LIMB_BITS != 0)
    ulongrem = (mp_limb_t) (ULONG_MAX >> GMP_LIMB_BITS * ulongsize) + 1;

  return un <= ulongsize || (up[ulongsize] < ulongrem && un == ulongsize + 1);
}

int
mpz_fits_ulong_p (const mpz_t u)
{
  mp_size_t us = u->_mp_size;

  return us >= 0 && mpn_absfits_ulong_p (u->_mp_d, us);
}

long int
mpz_get_si (const mpz_t u)
{
  unsigned long r = mpz_get_ui (u);
  unsigned long c = -LONG_MAX - LONG_MIN;

  if (u->_mp_size < 0)
    /* This expression is necessary to properly handle -LONG_MIN */
    return -(long) c - (long) ((r - c) & LONG_MAX);
  else
    return (long) (r & LONG_MAX);
}

unsigned long int
mpz_get_ui (const mpz_t u)
{
  if (GMP_LIMB_BITS < GMP_ULONG_BITS)
    {
      int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS;
      unsigned long r = 0;
      mp_size_t n = GMP_ABS (u->_mp_size);
      n = GMP_MIN (n, 1 + (mp_size_t) (GMP_ULONG_BITS - 1) / GMP_LIMB_BITS);
      while (--n >= 0)
	r = (r << LOCAL_GMP_LIMB_BITS) + u->_mp_d[n];
      return r;
    }

  return u->_mp_size == 0 ? 0 : u->_mp_d[0];
}

size_t
mpz_size (const mpz_t u)
{
  return GMP_ABS (u->_mp_size);
}

mp_limb_t
mpz_getlimbn (const mpz_t u, mp_size_t n)
{
  if (n >= 0 && n < GMP_ABS (u->_mp_size))
    return u->_mp_d[n];
  else
    return 0;
}

void
mpz_realloc2 (mpz_t x, mp_bitcnt_t n)
{
  mpz_realloc (x, 1 + (n - (n != 0)) / GMP_LIMB_BITS);
}

mp_srcptr
mpz_limbs_read (mpz_srcptr x)
{
  return x->_mp_d;
}

mp_ptr
mpz_limbs_modify (mpz_t x, mp_size_t n)
{
  assert (n > 0);
  return MPZ_REALLOC (x, n);
}

mp_ptr
mpz_limbs_write (mpz_t x, mp_size_t n)
{
  return mpz_limbs_modify (x, n);
}

void
mpz_limbs_finish (mpz_t x, mp_size_t xs)
{
  mp_size_t xn;
  xn = mpn_normalized_size (x->_mp_d, GMP_ABS (xs));
  x->_mp_size = xs < 0 ? -xn : xn;
}

static mpz_srcptr
mpz_roinit_normal_n (mpz_t x, mp_srcptr xp, mp_size_t xs)
{
  x->_mp_alloc = 0;
  x->_mp_d = (mp_ptr) xp;
  x->_mp_size = xs;
  return x;
}

mpz_srcptr
mpz_roinit_n (mpz_t x, mp_srcptr xp, mp_size_t xs)
{
  mpz_roinit_normal_n (x, xp, xs);
  mpz_limbs_finish (x, xs);
  return x;
}


/* Conversions and comparison to double. */
void
mpz_set_d (mpz_t r, double x)
{
  int sign;
  mp_ptr rp;
  mp_size_t rn, i;
  double B;
  double Bi;
  mp_limb_t f;

  /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
     zero or infinity. */
  if (x != x || x == x * 0.5)
    {
      r->_mp_size = 0;
      return;
    }

  sign = x < 0.0 ;
  if (sign)
    x = - x;

  if (x < 1.0)
    {
      r->_mp_size = 0;
      return;
    }
  B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1);
  Bi = 1.0 / B;
  for (rn = 1; x >= B; rn++)
    x *= Bi;

  rp = MPZ_REALLOC (r, rn);

  f = (mp_limb_t) x;
  x -= f;
  assert (x < 1.0);
  i = rn-1;
  rp[i] = f;
  while (--i >= 0)
    {
      x = B * x;
      f = (mp_limb_t) x;
      x -= f;
      assert (x < 1.0);
      rp[i] = f;
    }

  r->_mp_size = sign ? - rn : rn;
}

void
mpz_init_set_d (mpz_t r, double x)
{
  mpz_init (r);
  mpz_set_d (r, x);
}

double
mpz_get_d (const mpz_t u)
{
  int m;
  mp_limb_t l;
  mp_size_t un;
  double x;
  double B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1);

  un = GMP_ABS (u->_mp_size);

  if (un == 0)
    return 0.0;

  l = u->_mp_d[--un];
  gmp_clz (m, l);
  m = m + GMP_DBL_MANT_BITS - GMP_LIMB_BITS;
  if (m < 0)
    l &= GMP_LIMB_MAX << -m;

  for (x = l; --un >= 0;)
    {
      x = B*x;
      if (m > 0) {
	l = u->_mp_d[un];
	m -= GMP_LIMB_BITS;
	if (m < 0)
	  l &= GMP_LIMB_MAX << -m;
	x += l;
      }
    }

  if (u->_mp_size < 0)
    x = -x;

  return x;
}

int
mpz_cmpabs_d (const mpz_t x, double d)
{
  mp_size_t xn;
  double B, Bi;
  mp_size_t i;

  xn = x->_mp_size;
  d = GMP_ABS (d);

  if (xn != 0)
    {
      xn = GMP_ABS (xn);

      B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1);
      Bi = 1.0 / B;

      /* Scale d so it can be compared with the top limb. */
      for (i = 1; i < xn; i++)
	d *= Bi;

      if (d >= B)
	return -1;

      /* Compare floor(d) to top limb, subtract and cancel when equal. */
      for (i = xn; i-- > 0;)
	{
	  mp_limb_t f, xl;

	  f = (mp_limb_t) d;
	  xl = x->_mp_d[i];
	  if (xl > f)
	    return 1;
	  else if (xl < f)
	    return -1;
	  d = B * (d - f);
	}
    }
  return - (d > 0.0);
}

int
mpz_cmp_d (const mpz_t x, double d)
{
  if (x->_mp_size < 0)
    {
      if (d >= 0.0)
	return -1;
      else
	return -mpz_cmpabs_d (x, d);
    }
  else
    {
      if (d < 0.0)
	return 1;
      else
	return mpz_cmpabs_d (x, d);
    }
}


/* MPZ comparisons and the like. */
int
mpz_sgn (const mpz_t u)
{
  return GMP_CMP (u->_mp_size, 0);
}

int
mpz_cmp_si (const mpz_t u, long v)
{
  mp_size_t usize = u->_mp_size;

  if (v >= 0)
    return mpz_cmp_ui (u, v);
  else if (usize >= 0)
    return 1;
  else
    return - mpz_cmpabs_ui (u, GMP_NEG_CAST (unsigned long int, v));
}

int
mpz_cmp_ui (const mpz_t u, unsigned long v)
{
  mp_size_t usize = u->_mp_size;

  if (usize < 0)
    return -1;
  else
    return mpz_cmpabs_ui (u, v);
}

int
mpz_cmp (const mpz_t a, const mpz_t b)
{
  mp_size_t asize = a->_mp_size;
  mp_size_t bsize = b->_mp_size;

  if (asize != bsize)
    return (asize < bsize) ? -1 : 1;
  else if (asize >= 0)
    return mpn_cmp (a->_mp_d, b->_mp_d, asize);
  else
    return mpn_cmp (b->_mp_d, a->_mp_d, -asize);
}

int
mpz_cmpabs_ui (const mpz_t u, unsigned long v)
{
  mp_size_t un = GMP_ABS (u->_mp_size);

  if (! mpn_absfits_ulong_p (u->_mp_d, un))
    return 1;
  else
    {
      unsigned long uu = mpz_get_ui (u);
      return GMP_CMP(uu, v);
    }
}

int
mpz_cmpabs (const mpz_t u, const mpz_t v)
{
  return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size),
		   v->_mp_d, GMP_ABS (v->_mp_size));
}

void
mpz_abs (mpz_t r, const mpz_t u)
{
  mpz_set (r, u);
  r->_mp_size = GMP_ABS (r->_mp_size);
}

void
mpz_neg (mpz_t r, const mpz_t u)
{
  mpz_set (r, u);
  r->_mp_size = -r->_mp_size;
}

void
mpz_swap (mpz_t u, mpz_t v)
{
  MP_SIZE_T_SWAP (u->_mp_size, v->_mp_size);
  MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc);
  MP_PTR_SWAP (u->_mp_d, v->_mp_d);
}


/* MPZ addition and subtraction */


void
mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b)
{
  mpz_t bb;
  mpz_init_set_ui (bb, b);
  mpz_add (r, a, bb);
  mpz_clear (bb);
}

void
mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
{
  mpz_ui_sub (r, b, a);
  mpz_neg (r, r);
}

void
mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b)
{
  mpz_neg (r, b);
  mpz_add_ui (r, r, a);
}

static mp_size_t
mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b)
{
  mp_size_t an = GMP_ABS (a->_mp_size);
  mp_size_t bn = GMP_ABS (b->_mp_size);
  mp_ptr rp;
  mp_limb_t cy;

  if (an < bn)
    {
      MPZ_SRCPTR_SWAP (a, b);
      MP_SIZE_T_SWAP (an, bn);
    }

  rp = MPZ_REALLOC (r, an + 1);
  cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);

  rp[an] = cy;

  return an + cy;
}

static mp_size_t
mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b)
{
  mp_size_t an = GMP_ABS (a->_mp_size);
  mp_size_t bn = GMP_ABS (b->_mp_size);
  int cmp;
  mp_ptr rp;

  cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn);
  if (cmp > 0)
    {
      rp = MPZ_REALLOC (r, an);
      gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn));
      return mpn_normalized_size (rp, an);
    }
  else if (cmp < 0)
    {
      rp = MPZ_REALLOC (r, bn);
      gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an));
      return -mpn_normalized_size (rp, bn);
    }
  else
    return 0;
}

void
mpz_add (mpz_t r, const mpz_t a, const mpz_t b)
{
  mp_size_t rn;

  if ( (a->_mp_size ^ b->_mp_size) >= 0)
    rn = mpz_abs_add (r, a, b);
  else
    rn = mpz_abs_sub (r, a, b);

  r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
}

void
mpz_sub (mpz_t r, const mpz_t a, const mpz_t b)
{
  mp_size_t rn;

  if ( (a->_mp_size ^ b->_mp_size) >= 0)
    rn = mpz_abs_sub (r, a, b);
  else
    rn = mpz_abs_add (r, a, b);

  r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
}


/* MPZ multiplication */
void
mpz_mul_si (mpz_t r, const mpz_t u, long int v)
{
  if (v < 0)
    {
      mpz_mul_ui (r, u, GMP_NEG_CAST (unsigned long int, v));
      mpz_neg (r, r);
    }
  else
    mpz_mul_ui (r, u, v);
}

void
mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v)
{
  mpz_t vv;
  mpz_init_set_ui (vv, v);
  mpz_mul (r, u, vv);
  mpz_clear (vv);
  return;
}

void
mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
{
  int sign;
  mp_size_t un, vn, rn;
  mpz_t t;
  mp_ptr tp;

  un = u->_mp_size;
  vn = v->_mp_size;

  if (un == 0 || vn == 0)
    {
      r->_mp_size = 0;
      return;
    }

  sign = (un ^ vn) < 0;

  un = GMP_ABS (un);
  vn = GMP_ABS (vn);

  mpz_init2 (t, (un + vn) * GMP_LIMB_BITS);

  tp = t->_mp_d;
  if (un >= vn)
    mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn);
  else
    mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un);

  rn = un + vn;
  rn -= tp[rn-1] == 0;

  t->_mp_size = sign ? - rn : rn;
  mpz_swap (r, t);
  mpz_clear (t);
}

void
mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
{
  mp_size_t un, rn;
  mp_size_t limbs;
  unsigned shift;
  mp_ptr rp;

  un = GMP_ABS (u->_mp_size);
  if (un == 0)
    {
      r->_mp_size = 0;
      return;
    }

  limbs = bits / GMP_LIMB_BITS;
  shift = bits % GMP_LIMB_BITS;

  rn = un + limbs + (shift > 0);
  rp = MPZ_REALLOC (r, rn);
  if (shift > 0)
    {
      mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift);
      rp[rn-1] = cy;
      rn -= (cy == 0);
    }
  else
    mpn_copyd (rp + limbs, u->_mp_d, un);

  mpn_zero (rp, limbs);

  r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
}

void
mpz_addmul_ui (mpz_t r, const mpz_t u, unsigned long int v)
{
  mpz_t t;
  mpz_init_set_ui (t, v);
  mpz_mul (t, u, t);
  mpz_add (r, r, t);
  mpz_clear (t);
}

void
mpz_submul_ui (mpz_t r, const mpz_t u, unsigned long int v)
{
  mpz_t t;
  mpz_init_set_ui (t, v);
  mpz_mul (t, u, t);
  mpz_sub (r, r, t);
  mpz_clear (t);
}

void
mpz_addmul (mpz_t r, const mpz_t u, const mpz_t v)
{
  mpz_t t;
  mpz_init (t);
  mpz_mul (t, u, v);
  mpz_add (r, r, t);
  mpz_clear (t);
}

void
mpz_submul (mpz_t r, const mpz_t u, const mpz_t v)
{
  mpz_t t;
  mpz_init (t);
  mpz_mul (t, u, v);
  mpz_sub (r, r, t);
  mpz_clear (t);
}


/* MPZ division */
enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };

/* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */
static int
mpz_div_qr (mpz_t q, mpz_t r,
	    const mpz_t n, const mpz_t d, enum mpz_div_round_mode mode)
{
  mp_size_t ns, ds, nn, dn, qs;
  ns = n->_mp_size;
  ds = d->_mp_size;

  if (ds == 0)
    gmp_die("mpz_div_qr: Divide by zero.");

  if (ns == 0)
    {
      if (q)
	q->_mp_size = 0;
      if (r)
	r->_mp_size = 0;
      return 0;
    }

  nn = GMP_ABS (ns);
  dn = GMP_ABS (ds);

  qs = ds ^ ns;

  if (nn < dn)
    {
      if (mode == GMP_DIV_CEIL && qs >= 0)
	{
	  /* q = 1, r = n - d */
	  if (r)
	    mpz_sub (r, n, d);
	  if (q)
	    mpz_set_ui (q, 1);
	}
      else if (mode == GMP_DIV_FLOOR && qs < 0)
	{
	  /* q = -1, r = n + d */
	  if (r)
	    mpz_add (r, n, d);
	  if (q)
	    mpz_set_si (q, -1);
	}
      else
	{
	  /* q = 0, r = d */
	  if (r)
	    mpz_set (r, n);
	  if (q)
	    q->_mp_size = 0;
	}
      return 1;
    }
  else
    {
      mp_ptr np, qp;
      mp_size_t qn, rn;
      mpz_t tq, tr;

      mpz_init_set (tr, n);
      np = tr->_mp_d;

      qn = nn - dn + 1;

      if (q)
	{
	  mpz_init2 (tq, qn * GMP_LIMB_BITS);
	  qp = tq->_mp_d;
	}
      else
	qp = NULL;

      mpn_div_qr (qp, np, nn, d->_mp_d, dn);

      if (qp)
	{
	  qn -= (qp[qn-1] == 0);

	  tq->_mp_size = qs < 0 ? -qn : qn;
	}
      rn = mpn_normalized_size (np, dn);
      tr->_mp_size = ns < 0 ? - rn : rn;

      if (mode == GMP_DIV_FLOOR && qs < 0 && rn != 0)
	{
	  if (q)
	    mpz_sub_ui (tq, tq, 1);
	  if (r)
	    mpz_add (tr, tr, d);
	}
      else if (mode == GMP_DIV_CEIL && qs >= 0 && rn != 0)
	{
	  if (q)
	    mpz_add_ui (tq, tq, 1);
	  if (r)
	    mpz_sub (tr, tr, d);
	}

      if (q)
	{
	  mpz_swap (tq, q);
	  mpz_clear (tq);
	}
      if (r)
	mpz_swap (tr, r);

      mpz_clear (tr);

      return rn != 0;
    }
}

void
mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
{
  mpz_div_qr (q, r, n, d, GMP_DIV_CEIL);
}

void
mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
{
  mpz_div_qr (q, r, n, d, GMP_DIV_FLOOR);
}

void
mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
{
  mpz_div_qr (q, r, n, d, GMP_DIV_TRUNC);
}

void
mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
{
  mpz_div_qr (q, NULL, n, d, GMP_DIV_CEIL);
}

void
mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
{
  mpz_div_qr (q, NULL, n, d, GMP_DIV_FLOOR);
}

void
mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
{
  mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC);
}

void
mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
{
  mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
}

void
mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
{
  mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
}

void
mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
{
  mpz_div_qr (NULL, r, n, d, GMP_DIV_TRUNC);
}

void
mpz_mod (mpz_t r, const mpz_t n, const mpz_t d)
{
  mpz_div_qr (NULL, r, n, d, d->_mp_size >= 0 ? GMP_DIV_FLOOR : GMP_DIV_CEIL);
}

static void
mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
		enum mpz_div_round_mode mode)
{
  mp_size_t un, qn;
  mp_size_t limb_cnt;
  mp_ptr qp;
  int adjust;

  un = u->_mp_size;
  if (un == 0)
    {
      q->_mp_size = 0;
      return;
    }
  limb_cnt = bit_index / GMP_LIMB_BITS;
  qn = GMP_ABS (un) - limb_cnt;
  bit_index %= GMP_LIMB_BITS;

  if (mode == ((un > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* un != 0 here. */
    /* Note: Below, the final indexing at limb_cnt is valid because at
       that point we have qn > 0. */
    adjust = (qn <= 0
	      || !mpn_zero_p (u->_mp_d, limb_cnt)
	      || (u->_mp_d[limb_cnt]
		  & (((mp_limb_t) 1 << bit_index) - 1)));
  else
    adjust = 0;

  if (qn <= 0)
    qn = 0;
  else
    {
      qp = MPZ_REALLOC (q, qn);

      if (bit_index != 0)
	{
	  mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index);
	  qn -= qp[qn - 1] == 0;
	}
      else
	{
	  mpn_copyi (qp, u->_mp_d + limb_cnt, qn);
	}
    }

  q->_mp_size = qn;

  if (adjust)
    mpz_add_ui (q, q, 1);
  if (un < 0)
    mpz_neg (q, q);
}

static void
mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
		enum mpz_div_round_mode mode)
{
  mp_size_t us, un, rn;
  mp_ptr rp;
  mp_limb_t mask;

  us = u->_mp_size;
  if (us == 0 || bit_index == 0)
    {
      r->_mp_size = 0;
      return;
    }
  rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
  assert (rn > 0);

  rp = MPZ_REALLOC (r, rn);
  un = GMP_ABS (us);

  mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index);

  if (rn > un)
    {
      /* Quotient (with truncation) is zero, and remainder is
	 non-zero */
      if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
	{
	  /* Have to negate and sign extend. */
	  mp_size_t i;

	  gmp_assert_nocarry (! mpn_neg (rp, u->_mp_d, un));
	  for (i = un; i < rn - 1; i++)
	    rp[i] = GMP_LIMB_MAX;

	  rp[rn-1] = mask;
	  us = -us;
	}
      else
	{
	  /* Just copy */
	  if (r != u)
	    mpn_copyi (rp, u->_mp_d, un);

	  rn = un;
	}
    }
  else
    {
      if (r != u)
	mpn_copyi (rp, u->_mp_d, rn - 1);

      rp[rn-1] = u->_mp_d[rn-1] & mask;

      if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
	{
	  /* If r != 0, compute 2^{bit_count} - r. */
	  mpn_neg (rp, rp, rn);

	  rp[rn-1] &= mask;

	  /* us is not used for anything else, so we can modify it
	     here to indicate flipped sign. */
	  us = -us;
	}
    }
  rn = mpn_normalized_size (rp, rn);
  r->_mp_size = us < 0 ? -rn : rn;
}

void
mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
{
  mpz_div_q_2exp (r, u, cnt, GMP_DIV_CEIL);
}

void
mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
{
  mpz_div_q_2exp (r, u, cnt, GMP_DIV_FLOOR);
}

void
mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
{
  mpz_div_q_2exp (r, u, cnt, GMP_DIV_TRUNC);
}

void
mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
{
  mpz_div_r_2exp (r, u, cnt, GMP_DIV_CEIL);
}

void
mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
{
  mpz_div_r_2exp (r, u, cnt, GMP_DIV_FLOOR);
}

void
mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
{
  mpz_div_r_2exp (r, u, cnt, GMP_DIV_TRUNC);
}

void
mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d)
{
  gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC));
}

int
mpz_divisible_p (const mpz_t n, const mpz_t d)
{
  return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
}

int
mpz_congruent_p (const mpz_t a, const mpz_t b, const mpz_t m)
{
  mpz_t t;
  int res;

  /* a == b (mod 0) iff a == b */
  if (mpz_sgn (m) == 0)
    return (mpz_cmp (a, b) == 0);

  mpz_init (t);
  mpz_sub (t, a, b);
  res = mpz_divisible_p (t, m);
  mpz_clear (t);

  return res;
}

static unsigned long
mpz_div_qr_ui (mpz_t q, mpz_t r,
	       const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
{
  unsigned long ret;
  mpz_t rr, dd;

  mpz_init (rr);
  mpz_init_set_ui (dd, d);
  mpz_div_qr (q, rr, n, dd, mode);
  mpz_clear (dd);
  ret = mpz_get_ui (rr);

  if (r)
    mpz_swap (r, rr);
  mpz_clear (rr);

  return ret;
}

unsigned long
mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (q, r, n, d, GMP_DIV_CEIL);
}

unsigned long
mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (q, r, n, d, GMP_DIV_FLOOR);
}

unsigned long
mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (q, r, n, d, GMP_DIV_TRUNC);
}

unsigned long
mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_CEIL);
}

unsigned long
mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_FLOOR);
}

unsigned long
mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC);
}

unsigned long
mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_CEIL);
}
unsigned long
mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
}
unsigned long
mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_TRUNC);
}

unsigned long
mpz_cdiv_ui (const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_CEIL);
}

unsigned long
mpz_fdiv_ui (const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_FLOOR);
}

unsigned long
mpz_tdiv_ui (const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC);
}

unsigned long
mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
}

void
mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d)
{
  gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC));
}

int
mpz_divisible_ui_p (const mpz_t n, unsigned long d)
{
  return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
}


/* GCD */
static mp_limb_t
mpn_gcd_11 (mp_limb_t u, mp_limb_t v)
{
  unsigned shift;

  assert ( (u | v) > 0);

  if (u == 0)
    return v;
  else if (v == 0)
    return u;

  gmp_ctz (shift, u | v);

  u >>= shift;
  v >>= shift;

  if ( (u & 1) == 0)
    MP_LIMB_T_SWAP (u, v);

  while ( (v & 1) == 0)
    v >>= 1;

  while (u != v)
    {
      if (u > v)
	{
	  u -= v;
	  do
	    u >>= 1;
	  while ( (u & 1) == 0);
	}
      else
	{
	  v -= u;
	  do
	    v >>= 1;
	  while ( (v & 1) == 0);
	}
    }
  return u << shift;
}

unsigned long
mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v)
{
  mpz_t t;
  mpz_init_set_ui(t, v);
  mpz_gcd (t, u, t);
  if (v > 0)
    v = mpz_get_ui (t);

  if (g)
    mpz_swap (t, g);

  mpz_clear (t);

  return v;
}

static mp_bitcnt_t
mpz_make_odd (mpz_t r)
{
  mp_bitcnt_t shift;

  assert (r->_mp_size > 0);
  /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */
  shift = mpn_common_scan (r->_mp_d[0], 0, r->_mp_d, 0, 0);
  mpz_tdiv_q_2exp (r, r, shift);

  return shift;
}

void
mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
{
  mpz_t tu, tv;
  mp_bitcnt_t uz, vz, gz;

  if (u->_mp_size == 0)
    {
      mpz_abs (g, v);
      return;
    }
  if (v->_mp_size == 0)
    {
      mpz_abs (g, u);
      return;
    }

  mpz_init (tu);
  mpz_init (tv);

  mpz_abs (tu, u);
  uz = mpz_make_odd (tu);
  mpz_abs (tv, v);
  vz = mpz_make_odd (tv);
  gz = GMP_MIN (uz, vz);

  if (tu->_mp_size < tv->_mp_size)
    mpz_swap (tu, tv);

  mpz_tdiv_r (tu, tu, tv);
  if (tu->_mp_size == 0)
    {
      mpz_swap (g, tv);
    }
  else
    for (;;)
      {
	int c;

	mpz_make_odd (tu);
	c = mpz_cmp (tu, tv);
	if (c == 0)
	  {
	    mpz_swap (g, tu);
	    break;
	  }
	if (c < 0)
	  mpz_swap (tu, tv);

	if (tv->_mp_size == 1)
	  {
	    mp_limb_t vl = tv->_mp_d[0];
	    mp_limb_t ul = mpz_tdiv_ui (tu, vl);
	    mpz_set_ui (g, mpn_gcd_11 (ul, vl));
	    break;
	  }
	mpz_sub (tu, tu, tv);
      }
  mpz_clear (tu);
  mpz_clear (tv);
  mpz_mul_2exp (g, g, gz);
}

void
mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
{
  mpz_t tu, tv, s0, s1, t0, t1;
  mp_bitcnt_t uz, vz, gz;
  mp_bitcnt_t power;

  if (u->_mp_size == 0)
    {
      /* g = 0 u + sgn(v) v */
      signed long sign = mpz_sgn (v);
      mpz_abs (g, v);
      if (s)
	s->_mp_size = 0;
      if (t)
	mpz_set_si (t, sign);
      return;
    }

  if (v->_mp_size == 0)
    {
      /* g = sgn(u) u + 0 v */
      signed long sign = mpz_sgn (u);
      mpz_abs (g, u);
      if (s)
	mpz_set_si (s, sign);
      if (t)
	t->_mp_size = 0;
      return;
    }

  mpz_init (tu);
  mpz_init (tv);
  mpz_init (s0);
  mpz_init (s1);
  mpz_init (t0);
  mpz_init (t1);

  mpz_abs (tu, u);
  uz = mpz_make_odd (tu);
  mpz_abs (tv, v);
  vz = mpz_make_odd (tv);
  gz = GMP_MIN (uz, vz);

  uz -= gz;
  vz -= gz;

  /* Cofactors corresponding to odd gcd. gz handled later. */
  if (tu->_mp_size < tv->_mp_size)
    {
      mpz_swap (tu, tv);
      MPZ_SRCPTR_SWAP (u, v);
      MPZ_PTR_SWAP (s, t);
      MP_BITCNT_T_SWAP (uz, vz);
    }

  /* Maintain
   *
   * u = t0 tu + t1 tv
   * v = s0 tu + s1 tv
   *
   * where u and v denote the inputs with common factors of two
   * eliminated, and det (s0, t0; s1, t1) = 2^p. Then
   *
   * 2^p tu =  s1 u - t1 v
   * 2^p tv = -s0 u + t0 v
   */

  /* After initial division, tu = q tv + tu', we have
   *
   * u = 2^uz (tu' + q tv)
   * v = 2^vz tv
   *
   * or
   *
   * t0 = 2^uz, t1 = 2^uz q
   * s0 = 0,    s1 = 2^vz
   */

  mpz_setbit (t0, uz);
  mpz_tdiv_qr (t1, tu, tu, tv);
  mpz_mul_2exp (t1, t1, uz);

  mpz_setbit (s1, vz);
  power = uz + vz;

  if (tu->_mp_size > 0)
    {
      mp_bitcnt_t shift;
      shift = mpz_make_odd (tu);
      mpz_mul_2exp (t0, t0, shift);
      mpz_mul_2exp (s0, s0, shift);
      power += shift;

      for (;;)
	{
	  int c;
	  c = mpz_cmp (tu, tv);
	  if (c == 0)
	    break;

	  if (c < 0)
	    {
	      /* tv = tv' + tu
	       *
	       * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv'
	       * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */

	      mpz_sub (tv, tv, tu);
	      mpz_add (t0, t0, t1);
	      mpz_add (s0, s0, s1);

	      shift = mpz_make_odd (tv);
	      mpz_mul_2exp (t1, t1, shift);
	      mpz_mul_2exp (s1, s1, shift);
	    }
	  else
	    {
	      mpz_sub (tu, tu, tv);
	      mpz_add (t1, t0, t1);
	      mpz_add (s1, s0, s1);

	      shift = mpz_make_odd (tu);
	      mpz_mul_2exp (t0, t0, shift);
	      mpz_mul_2exp (s0, s0, shift);
	    }
	  power += shift;
	}
    }

  /* Now tv = odd part of gcd, and -s0 and t0 are corresponding
     cofactors. */

  mpz_mul_2exp (tv, tv, gz);
  mpz_neg (s0, s0);

  /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To
     adjust cofactors, we need u / g and v / g */

  mpz_divexact (s1, v, tv);
  mpz_abs (s1, s1);
  mpz_divexact (t1, u, tv);
  mpz_abs (t1, t1);

  while (power-- > 0)
    {
      /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */
      if (mpz_odd_p (s0) || mpz_odd_p (t0))
	{
	  mpz_sub (s0, s0, s1);
	  mpz_add (t0, t0, t1);
	}
      assert (mpz_even_p (t0) && mpz_even_p (s0));
      mpz_tdiv_q_2exp (s0, s0, 1);
      mpz_tdiv_q_2exp (t0, t0, 1);
    }

  /* Arrange so that |s| < |u| / 2g */
  mpz_add (s1, s0, s1);
  if (mpz_cmpabs (s0, s1) > 0)
    {
      mpz_swap (s0, s1);
      mpz_sub (t0, t0, t1);
    }
  if (u->_mp_size < 0)
    mpz_neg (s0, s0);
  if (v->_mp_size < 0)
    mpz_neg (t0, t0);

  mpz_swap (g, tv);
  if (s)
    mpz_swap (s, s0);
  if (t)
    mpz_swap (t, t0);

  mpz_clear (tu);
  mpz_clear (tv);
  mpz_clear (s0);
  mpz_clear (s1);
  mpz_clear (t0);
  mpz_clear (t1);
}

void
mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v)
{
  mpz_t g;

  if (u->_mp_size == 0 || v->_mp_size == 0)
    {
      r->_mp_size = 0;
      return;
    }

  mpz_init (g);

  mpz_gcd (g, u, v);
  mpz_divexact (g, u, g);
  mpz_mul (r, g, v);

  mpz_clear (g);
  mpz_abs (r, r);
}

void
mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v)
{
  if (v == 0 || u->_mp_size == 0)
    {
      r->_mp_size = 0;
      return;
    }

  v /= mpz_gcd_ui (NULL, u, v);
  mpz_mul_ui (r, u, v);

  mpz_abs (r, r);
}

int
mpz_invert (mpz_t r, const mpz_t u, const mpz_t m)
{
  mpz_t g, tr;
  int invertible;

  if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0)
    return 0;

  mpz_init (g);
  mpz_init (tr);

  mpz_gcdext (g, tr, NULL, u, m);
  invertible = (mpz_cmp_ui (g, 1) == 0);

  if (invertible)
    {
      if (tr->_mp_size < 0)
	{
	  if (m->_mp_size >= 0)
	    mpz_add (tr, tr, m);
	  else
	    mpz_sub (tr, tr, m);
	}
      mpz_swap (r, tr);
    }

  mpz_clear (g);
  mpz_clear (tr);
  return invertible;
}


/* Higher level operations (sqrt, pow and root) */

void
mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e)
{
  unsigned long bit;
  mpz_t tr;
  mpz_init_set_ui (tr, 1);

  bit = GMP_ULONG_HIGHBIT;
  do
    {
      mpz_mul (tr, tr, tr);
      if (e & bit)
	mpz_mul (tr, tr, b);
      bit >>= 1;
    }
  while (bit > 0);

  mpz_swap (r, tr);
  mpz_clear (tr);
}

void
mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e)
{
  mpz_t b;

  mpz_init_set_ui (b, blimb);
  mpz_pow_ui (r, b, e);
  mpz_clear (b);
}

void
mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
{
  mpz_t tr;
  mpz_t base;
  mp_size_t en, mn;
  mp_srcptr mp;
  struct gmp_div_inverse minv;
  unsigned shift;
  mp_ptr tp = NULL;

  en = GMP_ABS (e->_mp_size);
  mn = GMP_ABS (m->_mp_size);
  if (mn == 0)
    gmp_die ("mpz_powm: Zero modulo.");

  if (en == 0)
    {
      mpz_set_ui (r, 1);
      return;
    }

  mp = m->_mp_d;
  mpn_div_qr_invert (&minv, mp, mn);
  shift = minv.shift;

  if (shift > 0)
    {
      /* To avoid shifts, we do all our reductions, except the final
	 one, using a *normalized* m. */
      minv.shift = 0;

      tp = gmp_xalloc_limbs (mn);
      gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift));
      mp = tp;
    }

  mpz_init (base);

  if (e->_mp_size < 0)
    {
      if (!mpz_invert (base, b, m))
	gmp_die ("mpz_powm: Negative exponent and non-invertible base.");
    }
  else
    {
      mp_size_t bn;
      mpz_abs (base, b);

      bn = base->_mp_size;
      if (bn >= mn)
	{
	  mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv);
	  bn = mn;
	}

      /* We have reduced the absolute value. Now take care of the
	 sign. Note that we get zero represented non-canonically as
	 m. */
      if (b->_mp_size < 0)
	{
	  mp_ptr bp = MPZ_REALLOC (base, mn);
	  gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn));
	  bn = mn;
	}
      base->_mp_size = mpn_normalized_size (base->_mp_d, bn);
    }
  mpz_init_set_ui (tr, 1);

  while (--en >= 0)
    {
      mp_limb_t w = e->_mp_d[en];
      mp_limb_t bit;

      bit = GMP_LIMB_HIGHBIT;
      do
	{
	  mpz_mul (tr, tr, tr);
	  if (w & bit)
	    mpz_mul (tr, tr, base);
	  if (tr->_mp_size > mn)
	    {
	      mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
	      tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
	    }
	  bit >>= 1;
	}
      while (bit > 0);
    }

  /* Final reduction */
  if (tr->_mp_size >= mn)
    {
      minv.shift = shift;
      mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
      tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
    }
  if (tp)
    gmp_free (tp);

  mpz_swap (r, tr);
  mpz_clear (tr);
  mpz_clear (base);
}

void
mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
{
  mpz_t e;

  mpz_init_set_ui (e, elimb);
  mpz_powm (r, b, e, m);
  mpz_clear (e);
}

/* x=trunc(y^(1/z)), r=y-x^z */
void
mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
{
  int sgn;
  mpz_t t, u;

  sgn = y->_mp_size < 0;
  if ((~z & sgn) != 0)
    gmp_die ("mpz_rootrem: Negative argument, with even root.");
  if (z == 0)
    gmp_die ("mpz_rootrem: Zeroth root.");

  if (mpz_cmpabs_ui (y, 1) <= 0) {
    if (x)
      mpz_set (x, y);
    if (r)
      r->_mp_size = 0;
    return;
  }

  mpz_init (u);
  mpz_init (t);
  mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1);

  if (z == 2) /* simplify sqrt loop: z-1 == 1 */
    do {
      mpz_swap (u, t);			/* u = x */
      mpz_tdiv_q (t, y, u);		/* t = y/x */
      mpz_add (t, t, u);		/* t = y/x + x */
      mpz_tdiv_q_2exp (t, t, 1);	/* x'= (y/x + x)/2 */
    } while (mpz_cmpabs (t, u) < 0);	/* |x'| < |x| */
  else /* z != 2 */ {
    mpz_t v;

    mpz_init (v);
    if (sgn)
      mpz_neg (t, t);

    do {
      mpz_swap (u, t);			/* u = x */
      mpz_pow_ui (t, u, z - 1);		/* t = x^(z-1) */
      mpz_tdiv_q (t, y, t);		/* t = y/x^(z-1) */
      mpz_mul_ui (v, u, z - 1);		/* v = x*(z-1) */
      mpz_add (t, t, v);		/* t = y/x^(z-1) + x*(z-1) */
      mpz_tdiv_q_ui (t, t, z);		/* x'=(y/x^(z-1) + x*(z-1))/z */
    } while (mpz_cmpabs (t, u) < 0);	/* |x'| < |x| */

    mpz_clear (v);
  }

  if (r) {
    mpz_pow_ui (t, u, z);
    mpz_sub (r, y, t);
  }
  if (x)
    mpz_swap (x, u);
  mpz_clear (u);
  mpz_clear (t);
}

int
mpz_root (mpz_t x, const mpz_t y, unsigned long z)
{
  int res;
  mpz_t r;

  mpz_init (r);
  mpz_rootrem (x, r, y, z);
  res = r->_mp_size == 0;
  mpz_clear (r);

  return res;
}

/* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */
void
mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u)
{
  mpz_rootrem (s, r, u, 2);
}

void
mpz_sqrt (mpz_t s, const mpz_t u)
{
  mpz_rootrem (s, NULL, u, 2);
}

int
mpz_perfect_square_p (const mpz_t u)
{
  if (u->_mp_size <= 0)
    return (u->_mp_size == 0);
  else
    return mpz_root (NULL, u, 2);
}

int
mpn_perfect_square_p (mp_srcptr p, mp_size_t n)
{
  mpz_t t;

  assert (n > 0);
  assert (p [n-1] != 0);
  return mpz_root (NULL, mpz_roinit_normal_n (t, p, n), 2);
}

mp_size_t
mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n)
{
  mpz_t s, r, u;
  mp_size_t res;

  assert (n > 0);
  assert (p [n-1] != 0);

  mpz_init (r);
  mpz_init (s);
  mpz_rootrem (s, r, mpz_roinit_normal_n (u, p, n), 2);

  assert (s->_mp_size == (n+1)/2);
  mpn_copyd (sp, s->_mp_d, s->_mp_size);
  mpz_clear (s);
  res = r->_mp_size;
  if (rp)
    mpn_copyd (rp, r->_mp_d, res);
  mpz_clear (r);
  return res;
}

/* Combinatorics */

void
mpz_mfac_uiui (mpz_t x, unsigned long n, unsigned long m)
{
  mpz_set_ui (x, n + (n == 0));
  if (m + 1 < 2) return;
  while (n > m + 1)
    mpz_mul_ui (x, x, n -= m);
}

void
mpz_2fac_ui (mpz_t x, unsigned long n)
{
  mpz_mfac_uiui (x, n, 2);
}

void
mpz_fac_ui (mpz_t x, unsigned long n)
{
  mpz_mfac_uiui (x, n, 1);
}

void
mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
{
  mpz_t t;

  mpz_set_ui (r, k <= n);

  if (k > (n >> 1))
    k = (k <= n) ? n - k : 0;

  mpz_init (t);
  mpz_fac_ui (t, k);

  for (; k > 0; --k)
    mpz_mul_ui (r, r, n--);

  mpz_divexact (r, r, t);
  mpz_clear (t);
}


/* Primality testing */

/* Computes Kronecker (a/b) with odd b, a!=0 and GCD(a,b) = 1 */
/* Adapted from JACOBI_BASE_METHOD==4 in mpn/generic/jacbase.c */
static int
gmp_jacobi_coprime (mp_limb_t a, mp_limb_t b)
{
  int c, bit = 0;

  assert (b & 1);
  assert (a != 0);
  /* assert (mpn_gcd_11 (a, b) == 1); */

  /* Below, we represent a and b shifted right so that the least
     significant one bit is implicit. */
  b >>= 1;

  gmp_ctz(c, a);
  a >>= 1;

  do
    {
      a >>= c;
      /* (2/b) = -1 if b = 3 or 5 mod 8 */
      bit ^= c & (b ^ (b >> 1));
      if (a < b)
	{
	  bit ^= a & b;
	  a = b - a;
	  b -= a;
	}
      else
	{
	  a -= b;
	  assert (a != 0);
	}

      gmp_ctz(c, a);
      ++c;
    }
  while (b > 0);

  return bit & 1 ? -1 : 1;
}

static void
gmp_lucas_step_k_2k (mpz_t V, mpz_t Qk, const mpz_t n)
{
  mpz_mod (Qk, Qk, n);
  /* V_{2k} <- V_k ^ 2 - 2Q^k */
  mpz_mul (V, V, V);
  mpz_submul_ui (V, Qk, 2);
  mpz_tdiv_r (V, V, n);
  /* Q^{2k} = (Q^k)^2 */
  mpz_mul (Qk, Qk, Qk);
}

/* Computes V_k, Q^k (mod n) for the Lucas' sequence */
/* with P=1, Q=Q; k = (n>>b0)|1. */
/* Requires an odd n > 4; b0 > 0; -2*Q must not overflow a long */
/* Returns (U_k == 0) and sets V=V_k and Qk=Q^k. */
static int
gmp_lucas_mod (mpz_t V, mpz_t Qk, long Q,
	       mp_bitcnt_t b0, const mpz_t n)
{
  mp_bitcnt_t bs;
  mpz_t U;
  int res;

  assert (b0 > 0);
  assert (Q <= - (LONG_MIN / 2));
  assert (Q >= - (LONG_MAX / 2));
  assert (mpz_cmp_ui (n, 4) > 0);
  assert (mpz_odd_p (n));

  mpz_init_set_ui (U, 1); /* U1 = 1 */
  mpz_set_ui (V, 1); /* V1 = 1 */
  mpz_set_si (Qk, Q);

  for (bs = mpz_sizeinbase (n, 2) - 1; --bs >= b0;)
    {
      /* U_{2k} <- U_k * V_k */
      mpz_mul (U, U, V);
      /* V_{2k} <- V_k ^ 2 - 2Q^k */
      /* Q^{2k} = (Q^k)^2 */
      gmp_lucas_step_k_2k (V, Qk, n);

      /* A step k->k+1 is performed if the bit in $n$ is 1	*/
      /* mpz_tstbit(n,bs) or the the bit is 0 in $n$ but	*/
      /* should be 1 in $n+1$ (bs == b0)			*/
      if (b0 == bs || mpz_tstbit (n, bs))
	{
	  /* Q^{k+1} <- Q^k * Q */
	  mpz_mul_si (Qk, Qk, Q);
	  /* U_{k+1} <- (U_k + V_k) / 2 */
	  mpz_swap (U, V); /* Keep in V the old value of U_k */
	  mpz_add (U, U, V);
	  /* We have to compute U/2, so we need an even value, */
	  /* equivalent (mod n) */
	  if (mpz_odd_p (U))
	    mpz_add (U, U, n);
	  mpz_tdiv_q_2exp (U, U, 1);
	  /* V_{k+1} <-(D*U_k + V_k) / 2 =
			U_{k+1} + (D-1)/2*U_k = U_{k+1} - 2Q*U_k */
	  mpz_mul_si (V, V, -2*Q);
	  mpz_add (V, U, V);
	  mpz_tdiv_r (V, V, n);
	}
      mpz_tdiv_r (U, U, n);
    }

  res = U->_mp_size == 0;
  mpz_clear (U);
  return res;
}

/* Performs strong Lucas' test on x, with parameters suggested */
/* for the BPSW test. Qk is only passed to recycle a variable. */
/* Requires GCD (x,6) = 1.*/
static int
gmp_stronglucas (const mpz_t x, mpz_t Qk)
{
  mp_bitcnt_t b0;
  mpz_t V, n;
  mp_limb_t maxD, D; /* The absolute value is stored. */
  long Q;
  mp_limb_t tl;

  /* Test on the absolute value. */
  mpz_roinit_normal_n (n, x->_mp_d, GMP_ABS (x->_mp_size));

  assert (mpz_odd_p (n));
  /* assert (mpz_gcd_ui (NULL, n, 6) == 1); */
  if (mpz_root (Qk, n, 2))
    return 0; /* A square is composite. */

  /* Check Ds up to square root (in case, n is prime)
     or avoid overflows */
  maxD = (Qk->_mp_size == 1) ? Qk->_mp_d [0] - 1 : GMP_LIMB_MAX;

  D = 3;
  /* Search a D such that (D/n) = -1 in the sequence 5,-7,9,-11,.. */
  /* For those Ds we have (D/n) = (n/|D|) */
  do
    {
      if (D >= maxD)
	return 1 + (D != GMP_LIMB_MAX); /* (1 + ! ~ D) */
      D += 2;
      tl = mpz_tdiv_ui (n, D);
      if (tl == 0)
	return 0;
    }
  while (gmp_jacobi_coprime (tl, D) == 1);

  mpz_init (V);

  /* n-(D/n) = n+1 = d*2^{b0}, with d = (n>>b0) | 1 */
  b0 = mpz_scan0 (n, 0);

  /* D= P^2 - 4Q; P = 1; Q = (1-D)/4 */
  Q = (D & 2) ? (long) (D >> 2) + 1 : -(long) (D >> 2);

  if (! gmp_lucas_mod (V, Qk, Q, b0, n))	/* If Ud != 0 */
    while (V->_mp_size != 0 && --b0 != 0)	/* while Vk != 0 */
      /* V <- V ^ 2 - 2Q^k */
      /* Q^{2k} = (Q^k)^2 */
      gmp_lucas_step_k_2k (V, Qk, n);

  mpz_clear (V);
  return (b0 != 0);
}

static int
gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y,
		 const mpz_t q, mp_bitcnt_t k)
{
  assert (k > 0);

  /* Caller must initialize y to the base. */
  mpz_powm (y, y, q, n);

  if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
    return 1;

  while (--k > 0)
    {
      mpz_powm_ui (y, y, 2, n);
      if (mpz_cmp (y, nm1) == 0)
	return 1;
      /* y == 1 means that the previous y was a non-trivial square root
	 of 1 (mod n). y == 0 means that n is a power of the base.
	 In either case, n is not prime. */
      if (mpz_cmp_ui (y, 1) <= 0)
	return 0;
    }
  return 0;
}

/* This product is 0xc0cfd797, and fits in 32 bits. */
#define GMP_PRIME_PRODUCT \
  (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL)

/* Bit (p+1)/2 is set, for each odd prime <= 61 */
#define GMP_PRIME_MASK 0xc96996dcUL

int
mpz_probab_prime_p (const mpz_t n, int reps)
{
  mpz_t nm1;
  mpz_t q;
  mpz_t y;
  mp_bitcnt_t k;
  int is_prime;
  int j;

  /* Note that we use the absolute value of n only, for compatibility
     with the real GMP. */
  if (mpz_even_p (n))
    return (mpz_cmpabs_ui (n, 2) == 0) ? 2 : 0;

  /* Above test excludes n == 0 */
  assert (n->_mp_size != 0);

  if (mpz_cmpabs_ui (n, 64) < 0)
    return (GMP_PRIME_MASK >> (n->_mp_d[0] >> 1)) & 2;

  if (mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT) != 1)
    return 0;

  /* All prime factors are >= 31. */
  if (mpz_cmpabs_ui (n, 31*31) < 0)
    return 2;

  mpz_init (nm1);
  mpz_init (q);

  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
  mpz_abs (nm1, n);
  nm1->_mp_d[0] -= 1;
  k = mpz_scan1 (nm1, 0);
  mpz_tdiv_q_2exp (q, nm1, k);

  /* BPSW test */
  mpz_init_set_ui (y, 2);
  is_prime = gmp_millerrabin (n, nm1, y, q, k) && gmp_stronglucas (n, y);
  reps -= 24; /* skip the first 24 repetitions */

  /* Use Miller-Rabin, with a deterministic sequence of bases, a[j] =
     j^2 + j + 41 using Euler's polynomial. We potentially stop early,
     if a[j] >= n - 1. Since n >= 31*31, this can happen only if reps >
     30 (a[30] == 971 > 31*31 == 961). */

  for (j = 0; is_prime & (j < reps); j++)
    {
      mpz_set_ui (y, (unsigned long) j*j+j+41);
      if (mpz_cmp (y, nm1) >= 0)
	{
	  /* Don't try any further bases. This "early" break does not affect
	     the result for any reasonable reps value (<=5000 was tested) */
	  assert (j >= 30);
	  break;
	}
      is_prime = gmp_millerrabin (n, nm1, y, q, k);
    }
  mpz_clear (nm1);
  mpz_clear (q);
  mpz_clear (y);

  return is_prime;
}


/* Logical operations and bit manipulation. */

/* Numbers are treated as if represented in two's complement (and
   infinitely sign extended). For a negative values we get the two's
   complement from -x = ~x + 1, where ~ is bitwise complement.
   Negation transforms

     xxxx10...0

   into

     yyyy10...0

   where yyyy is the bitwise complement of xxxx. So least significant
   bits, up to and including the first one bit, are unchanged, and
   the more significant bits are all complemented.

   To change a bit from zero to one in a negative number, subtract the
   corresponding power of two from the absolute value. This can never
   underflow. To change a bit from one to zero, add the corresponding
   power of two, and this might overflow. E.g., if x = -001111, the
   two's complement is 110001. Clearing the least significant bit, we
   get two's complement 110000, and -010000. */

int
mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index)
{
  mp_size_t limb_index;
  unsigned shift;
  mp_size_t ds;
  mp_size_t dn;
  mp_limb_t w;
  int bit;

  ds = d->_mp_size;
  dn = GMP_ABS (ds);
  limb_index = bit_index / GMP_LIMB_BITS;
  if (limb_index >= dn)
    return ds < 0;

  shift = bit_index % GMP_LIMB_BITS;
  w = d->_mp_d[limb_index];
  bit = (w >> shift) & 1;

  if (ds < 0)
    {
      /* d < 0. Check if any of the bits below is set: If so, our bit
	 must be complemented. */
      if (shift > 0 && (mp_limb_t) (w << (GMP_LIMB_BITS - shift)) > 0)
	return bit ^ 1;
      while (--limb_index >= 0)
	if (d->_mp_d[limb_index] > 0)
	  return bit ^ 1;
    }
  return bit;
}

static void
mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index)
ofs	hex dump	ascii
0000	89 50 4e 47 0d 0a 1a 0a 00 00 00 0d 49 48 44 52 00 00 00 80 00 00 00 80 08 04 00 00 00 69 37 a9	.PNG........IHDR.............i7.
0020	40 00 00 00 02 73 42 49 54 08 08 55 ec 46 04 00 00 00 09 70 48 59 73 00 00 03 b1 00 00 03 b1 01	@....sBIT..U.F.....pHYs.........
0040	f5 83 ed 49 00 00 00 19 74 45 58 74 53 6f 66 74 77 61 72 65 00 77 77 77 2e 69 6e 6b 73 63 61 70	...I....tEXtSoftware.www.inkscap
0060	65 2e 6f 72 67 9b ee 3c 1a 00 00 08 f6 49 44 41 54 78 01 ed c1 0b ac 9d 05 7d 00 f0 ff bd 05 4a	e.org..<.....IDATx.......}.....J
0080	eb a5 0f a0 94 d7 e4 d1 0a ab a9 2e 8e 4c 8c 21 23 68 b6 69 d4 65 63 53 b7 c9 26 9a 39 12 83 c4	.............L.!#h.i.ecS..&.9...
00a0	b9 65 9b d3 0d aa 0b cb cc 66 94 c8 b2 6c 93 54 c4 31 86 89 59 04 25 33 8a 83 89 82 8f e9 80 a1	.e.......f...l.T.1..Y.%3........
00c0	8c 96 22 50 1e 6d ef 39 e7 7b 9f fb 38 bf f5 e3 50 28 ed b9 ed bd a7 f7 b6 bd c9 f7 fb 45 a3 d1	.."P.m.9.{..8...P(...........E..
00e0	68 34 1a 8d 46 a3 d1 68 34 1a 8d 46 a3 d1 68 34 1a 8d 46 a3 d1 68 34 e6 89 13 9d 6b 2c e6 c0 6a	h4..F..h4..F..h4..F..h4....k,..j
0100	e7 39 21 16 3f c7 ba ca 4f 65 c6 e5 1e f6 5e 4b e2 20 1c e3 4a 0f cb 8c cb fd c4 e5 46 63 f1 72	.9!.?...Oe....^K....J.......Fc.r
0120	8a 87 a4 5e 94 b8 c7 ca 38 00 27 fa 8e cc 8b 12 f7 5a 19 8b 93 11 3f 34 e9 a5 2a 77 1b 8d 19 18	...^....8.'......Z....?4..*w....
0140	f5 1d 5d 2f 55 b9 37 16 27 ef d4 b1 bf c4 e5 31 03 1f 92 da 5f db 2f c7 e2 e3 58 db 0d 76 7f cc	..]/U.7.'......1...._./...X..v..
0160	c0 93 06 fb 7e 2c 3e 2e 97 18 ac 8c 81 2c 97 1b 6c dc 2f c4 e2 62 c4 63 66 92 59 11 03 58 ab 63	....~,>......,..l./..b.cf.Y..X.c
0180	26 5f 8f c5 c5 52 b7 ba db 83 76 28 8d 9b b6 b7 b6 d3 63 00 eb 8c db db a4 96 d2 36 df f5 19 bf	&_...R....v(......c........6....
01a0	11 8b 95 63 fc a2 4f 69 49 ec d1 b2 2e 06 f0 6a e3 f6 e8 d8 e9 5a af b5 24 16 2b 6b 9c 63 b5 25	...c..OiI......j.....Z..$.+k.c.%
01c0	f1 1c 4b fd be 5d 4a b5 a7 bc 3c 06 70 a1 5d 6a b9 9d de 63 49 3c c7 b1 56 3b d7 9a 58 2c 6c b0	..K..]J...<.p.]j...cI<..V;..X,l.
01e0	c9 03 0a 2d e3 da 0a f7 fb 1b e7 c5 6e 56 fa 82 04 5b ad 8a 01 5c 64 1c 99 eb 2d 8f dd 6c f4 69	...-........nV...[...\d...-..l.i
0200	3f 96 1b b7 c3 b8 c2 fd 3e 6a 5d 1c cd 6c 74 87 44 a5 a3 ab af ab a3 2b 75 97 0d b1 9b 0f 48 3d	?.......>j]..lt.D......+u.....H=
0220	61 59 0c e0 12 85 d4 7b 62 37 af 71 9f 44 57 6a 4a df 94 44 25 f1 25 eb e3 68 64 d4 87 75 94 26	aY.....{b7.q.DWjJ..D%.%..hd..u.&
0240	ec 6f 5a 29 75 b5 91 08 bf e2 db 46 63 00 bf a4 f2 c6 08 4b fc ad 8e d2 b4 fd 4d 28 b5 bd 3f 8e	.oZ)u......Fc......K......M(..?.
0260	36 96 b8 59 4b 6e 66 85 dc 17 1d 17 61 2c 06 b2 cc 2b 22 2c f3 35 99 c2 cc 32 6d d7 1b 89 a3 89	6..YKnf.....a,...+",.5...2m.....
0280	2f c8 4c 38 b0 49 b9 9b e2 80 8c ba 4d 6e d2 81 95 32 9f 8e a3 87 f7 49 4d 38 b8 09 89 2b e3 79	/.L8.I......Mn...2.....IM8...+.y
02a0	7e c6 1f da ec 06 1f 70 46 3c cf c7 a4 26 1d 5c 29 75 69 1c 1d 9c aa 25 35 3b a9 96 d3 22 8c fa	~......pF<...&.\)ui....%5;..."..
02c0	6b 6d a9 52 29 d1 71 8d d1 08 e7 4b 14 66 27 33 6e 45 1c 0d 5c af 30 7b a5 cd 11 6e 94 c8 ed 91	km.R).q....K.f'3nE..\.0{...n....
02e0	69 fb e7 08 b7 a9 cc 5e 6a 53 1c 79 4e f7 88 c2 ec e5 3a de a5 ad b2 b7 ca b8 cb a4 ba 66 2f f1	i......^jS.yN.....:..........f/.
0300	80 d5 71 a4 f9 07 95 b9 69 d9 26 b3 af dc 56 a5 b9 49 6d 8a 23 cb a9 52 a5 b9 29 4d 9a b2 af ca	..q.....i.&...V..Im.#..R..)M....
0320	b4 cc dc 24 5a 56 c5 91 e4 ef 55 e6 aa 30 6d 90 29 a5 b9 ca 5c 1d 47 8e b5 52 a5 61 65 26 64 52	...$ZV....U..0m.)...\.G..R.ae&dR
0340	5d a9 61 25 5a 56 c5 91 e2 3a 95 e1 74 95 9e f4 c7 2e 76 89 3f f5 b4 ca b4 e1 64 fe 22 8e 0c 27	].a%ZV...:..t.....v.?.....d."..'
0360	4b 94 86 d1 53 d8 ec b8 78 9e a5 fe 4d 65 38 89 71 27 c4 91 e0 53 2a c3 a9 dc 6d 24 f6 62 89 1f	K...S...x...Me8.q'...S*...m$.b..
0380	28 0c 27 f7 e7 71 f8 39 49 aa 30 9c c4 05 b1 0f 17 e9 18 4e aa e5 84 38 dc fc 9d ca b0 76 c5 7e	(.'..q.9I.0........N...8.....v.~
03a0	bc 4c 6e 58 b9 3f 8b c3 cb 49 12 85 61 55 b1 1f 63 32 c3 4a b5 8c c5 e1 e4 13 4a c3 ab e2 39 46	.LnX.?...I..aU..c2.J......J...9F
03c0	ad b5 26 9e 63 85 d4 f0 0a 7f 12 87 8f 55 3a 72 c3 ab 22 ac f2 49 2d 6d 1d 2d 7f 65 a5 55 52 c3	..&.c........U:r.."..I-m.-.e.UR.
03e0	cb ec b2 3c 0e 17 d7 2a 1d 8a ca 06 cf 28 a4 6a 89 d2 76 17 48 1d 8a cc 87 e2 f0 b0 52 5b ee 50	...<...*.....(.j..v.H.......R[.P
0400	54 9e 51 d8 5b e1 29 99 43 91 da 69 79 1c 0e 3e ae b0 bf 9e 4c 2a d3 91 eb 39 b0 29 a5 7d 15 26	T.Q.[.).C..iy..>....L*...9.).}.&
0420	1c 58 4f ae a3 a3 23 35 65 7f 89 0f c6 c2 b3 52 5b 66 6f 85 ae 8e d2 0e 0f b8 d7 43 76 29 24 2a	.XO...#5e......R[fo........Cv)$*
0440	95 99 f4 14 f6 55 98 34 93 ae 52 a6 d0 f2 7f be e5 01 cf 28 75 74 e5 f6 96 da 61 59 2c 34 d7 c8	.....U.4..R........(ut....aY,4..
0460	bc 28 57 68 fb ac 37 3b 3e 5e 60 cc a5 fe 55 a2 50 3a 54 a5 4a ea df fd b6 d5 f1 02 c7 7a 83 eb	.(Wh..7;>^`...U.P:T.J........z..
0480	8d cb 64 5e d4 76 65 2c 2c 2b 8c 4b f4 f5 54 12 1f 33 16 03 39 c9 75 52 a5 43 51 c8 6c 76 46 0c	..d^.ve,,+.K..T..3..9.uR.CQ.lvF.
04a0	64 a9 3f d2 56 99 d2 97 78 d6 b2 58 48 3e 2a d1 37 a9 72 97 53 e2 80 ac f3 88 c2 a4 e1 54 1e b3	d.?.V...x..XH>*.7.r.S........T..
04c0	31 0e c8 4a ff 21 57 ea db e5 fd b1 70 8c d9 29 51 9b 50 f8 47 23 71 50 c6 7c 55 69 da 60 93 26	1..J.!W.....p..)Q.P.G#qP.\|Ui.`.&
04e0	0d d6 53 f9 b6 b1 38 28 a3 3e 21 53 aa a5 9e b2 34 16 8a bf d4 52 eb c9 7d 3e 66 c9 a8 db 54 f6	..S...8(.>!S....4....R..}>f...T.
0500	55 ea ea c8 15 3a 2a 5d fb 2a 3d 68 59 cc 92 6b 65 26 d4 76 b8 2a 16 86 53 4c 49 d5 2a df 73 5c	U....:].=hY..ke&.v...SLI..s\
0520	cc 9a a5 1e 56 d8 5b 2e f7 59 eb 23 8c 38 df 66 b9 cc de 72 3b 9c 16 b3 66 c4 ed 2a b5 cc a4 97	....V.[..Y.#.8.f...r;...f..*....
0540	c5 42 f0 79 a5 5a a6 63 6d cc 89 0d 52 85 3d 2a 99 77 c4 5e 5c 26 57 d9 63 5a ea e2 98 13 63 76	.B.y.Z.cm...R.=*.w.^\&W.cZ....cv
0560	c9 d0 53 fa 64 2c 04 8f 4b d5 72 57 c7 9c f9 9c ae 3d 4a d7 47 18 75 ad 7b dc ed ca d8 cd 0d 0a	..S.d,..K.rW.....=J.G.u.{.......
0580	7b 94 be 11 73 e6 0a 89 5a ee fe 98 7f 2e d4 56 cb b4 8c c5 9c 39 43 6e 42 df 63 ce 8b 70 8d ae	{...s...Z......V.....9CnB.c..p..
05a0	5a db 3b 23 fc bc 67 f5 f5 64 5e 15 73 e6 18 db a5 98 d6 71 5e cc 37 9b 74 d5 2a ff 14 43 71 a7	Z.;#..g..d^.s......q^.7.t.*..Cq.
05c0	4c 5f c7 da 08 df b4 c7 17 23 fc ac b6 be c4 4f 62 28 36 a9 d4 32 57 c5 7c 73 9f 44 ad e5 4d 31	L_.......#.....Ob(6..2W.\|s.D..M1
05e0	14 ef 53 ea 4b ac 89 f0 35 7b dc 1c e1 3c 6d 7d 5d 9b 62 28 5e a3 a3 56 f8 72 cc 37 89 0c 3d a5	..S.K...5{...<m}].b(^..V.r.7..=.
0600	e3 63 28 ce 94 e9 4b 9d 14 e1 37 3d a3 f6 98 8d 11 d6 6b e9 6b 7b 6d 0c c5 88 b6 0a 1d 8f c7 fc	.c(...K...7=......k.k{m.........
0620	b2 5c 6a 0a 95 9d 31 34 b9 69 b5 d4 ea d8 cd c5 6e b2 d9 f9 b1 9b 73 b5 f4 b5 ad 89 21 f9 be 04	.\j...14.i......n.....s.....!...
0640	85 c2 48 cc 27 67 d9 a9 96 7a 28 86 e6 09 99 5a 6a 65 ec c3 d9 c6 d5 7a 4a 23 31 24 5f 96 a2 27	..H.'g...z(....Zje.....zJ#1$_..'
0660	b1 22 e6 93 73 ec 54 4b dc 13 43 b3 55 ae b6 c5 58 ec c3 99 1e 51 9b 92 c4 d0 fc 8b 5c 2d b1 32	."..s.TK..C.U...X....Q......\-.2
0680	e6 93 b3 ec 50 cb fc 28 86 66 bb 02 d3 4a 23 b1 0f 2b 25 fa aa 18 9a db a5 e8 49 8c c5 7c b2 54	....P..(.f...J#..+%.......I..\|.T
06a0	ae 87 ae ed 31 34 85 29 14 b2 d8 8f 51 a5 29 b5 c2 58 0c c9 8f 24 e8 ca 63 be 19 97 a8 95 8e 8f	....14.)....Q.)..X...$..c.......
06c0	a1 78 b9 4c ad e3 7f 63 00 db 25 6a 89 0b 63 28 46 74 54 68 db 16 f3 cd 37 25 6a 85 37 c4 50 fc	.x.L...c..%j..c(FtTh....7%j.7.P.
06e0	9e 52 ad 74 53 0c e0 0e 99 da 84 8f c4 50 bc 52 aa 56 ba 25 e6 9b 8f a8 d4 2a 9f 8b a1 b8 5b ae	.R.tS........P.R.V.%.....*....[.
0700	96 b9 2c 06 f0 41 a5 5a e2 e1 18 8a 8f 9b 50 2b 5c 11 f3 cd 46 6d 3d 94 12 27 c6 9c 59 2f 33 85	..,..A.Z......P+\...Fm=..'..Y/3.
0720	09 99 93 63 00 eb 74 4c ab 65 5e 17 73 66 99 9d 32 4c 4a 9c 19 f3 cf 16 a9 5a e5 33 31 67 6e 57	...c..tL.e^.sf..2LJ......Z.31gnW
0740	a9 15 ee 8c 19 78 40 a6 56 fa 5e cc 99 0f 4b d5 72 f7 c5 42 70 85 54 ad 90 58 17 73 e2 f5 12 13	.....x@.V.^...K.r..Bp.T..X.s....
0760	e8 69 7b 4b cc c0 bb 74 d4 7a 52 97 c6 9c 58 a3 23 53 eb 78 47 2c 04 c7 79 46 aa 56 da 62 2c 66	.i{K...t.zR...X.#S.xG,..yF.V.b,f
0780	cd 5a 4f cb d5 0a 0f 19 8d 19 38 ce 13 32 b5 4c cb fa 98 35 c7 ba 4f a5 96 db 6a 34 16 86 df 95	.ZO.......8..2.L...5..O...j4....
07a0	9a 52 2b dd 61 34 66 c5 52 3f 54 aa 75 b5 5d 1c 07 e0 d7 74 4c aa 95 b6 59 15 b3 64 b3 4c 0f 5d	.R+.a4f.R?T.u.]....tL...Y..d.L.]
07c0	1d 6f 8e 85 e3 4e a5 da a4 d2 37 ac 88 83 b2 d6 7f ab f4 15 6e 8c 83 f0 55 a5 be ca a3 ce 8d 83	.o...N....7.........n...U.......
07e0	b2 d4 cd 0a 5d b5 c2 ad b1 90 9c ea 69 85 5a 4f e9 51 af 8c 03 f2 7a cf 2a f5 d4 4a 8f 5a 1e 07	....].......i.ZO.Q....z.*..J.Z..
0800	61 b5 27 15 fa 4a e3 de 14 07 e4 6c ff 23 d7 55 2b 6c b3 2a 16 96 0b b4 e4 fa 0a 99 1b 9d 1e 03	a.'..J.....l.#.U+l.*............
0820	59 ef 36 a9 5c 5f e1 59 e7 c4 2c d8 68 a7 42 5f 2e f1 9f 5e 1d 03 39 d1 75 52 95 be c2 4e e7 c7	Y.6.\_.Y..,.h.B_...^..9.uR...N..
0840	c2 73 91 96 42 5f a5 2b f7 15 bf e5 b4 78 81 b3 bc db 5d 32 95 49 7d a5 1d 5e 15 b3 e4 75 76 29	.s..B_.+.....x....]2.I}..^...uv)
0860	f4 4d 29 a5 be eb 0f bc c2 48 3c cf c9 7e dd ad 32 95 42 5f 69 97 9f 8b c3 c3 05 9e 52 98 d4 37	.M)......H<..~..2.B_i.......R..7
0880	a9 50 49 a5 9e b0 c5 53 72 a9 4a 6e 5a 5f 25 b7 d5 d9 31 07 36 f8 a9 d2 84 be 9e 4c 25 51 7a d6	.PI....Sr.JnZ_%...1.6......L%Qz.
08a0	36 db e5 32 99 52 57 5f 57 e1 51 1b e2 f0 71 92 af e8 c8 ed ad 92 48 95 a6 bd a8 a7 d0 71 83 13	6..2.RW_W.Q...q.......H......q..
08c0	62 8e 9c e8 4b 3a 32 2f 35 a1 50 79 a9 5c db ad 4e 88 c3 cd af 7a 5c 47 6e ca 60 13 4a 6d 0f ba	b...K:2/5.Py.\..N....z\Gn.`.Jm..
08e0	24 86 e4 ad b6 6a 2b 4c 1a 6c 52 a1 65 8b 37 c6 91 61 89 b7 fb 81 54 ae d0 31 a5 af d0 51 49 a5	$....j+L.lR.e.7..a....T..1...QI.
0900	be ee ad 46 e2 10 18 f5 76 df 92 c8 54 12 b9 da b4 42 47 29 91 fa 2f 6f 33 12 47 96 d3 5d e1 16	...F....v...T....BG)../o3.G..]..
0920	5b a4 2a a5 42 e2 41 37 fa 1d a7 c4 3c 71 ba 77 bb c5 8f a5 0a 85 52 c7 c3 6e f2 5e 6b e2 e8 62	[.*.B.A7....<q.w......R..n.^k..b
0940	89 d1 58 60 46 2d 89 46 a3 d1 68 34 1a 8d 46 a3 d1 68 34 1a 8d 46 a3 d1 68 34 1a 8d 46 a3 d1 98	..X`F-.F..h4..F..h4..F..h4..F...
0960	a3 ff 07 81 51 1c fa 86 78 0e 24 00 00 00 00 49 45 4e 44 ae 42 60 82	....Q...x.$....IEND.B`.