/* gmpmodule.c
 *
 * Python interface to the GMP multiple precision library,
 * written for GMP-2.0.
 *
 * Rewritten by Niels Möller, May 1996
 */

#include <assert.h>
#include <math.h>

#include "Python.h"
#include <gmp.h>
#include "longintrepr.h"

#ifdef __GNUC__
#define USE_ALLOCA 1
#endif

typedef struct
{
  PyObject_HEAD
  mpz_t z;
} PympzObject;

static PyObject *gmp_module = NULL; 

static struct gmp_options
{
  long debug;
} options = { 0 };

staticforward PyTypeObject Pympz_Type;
staticforward PyMethodDef Pympz_methods [];

/* Number of bits that are significant when converting an mpz value to
 * a float. Initialized in initgmp(). */
static unsigned
double_mantissa = 0;

#define Pympz_Check(v) ((v)->ob_type == &Pympz_Type)

static PympzObject *
Pympz_new(void)
{
  PympzObject * self;

  if (!(self = PyObject_NEW(PympzObject, &Pympz_Type)))
    return NULL;
  mpz_init(self->z);
  return self;
}

static void
Pympz_dealloc(PympzObject *self)
{
  if (options.debug)
    fprintf(stderr, "Pympz_dealloc: %p\n", self);
  mpz_clear(self->z);
  PyMem_DEL(self);
} /* Pympz_dealloc */

static PyObject *
Pygmp_set_debug(PyObject *self, PyObject *args)
{
  int old = options.debug;

  if (!PyArg_ParseTuple(args, "l", &options.debug))
    return NULL;
  return Py_BuildValue("i", old);
}


/* CONVERSIONS */
static PympzObject *
int2mpz(PyObject *i)
{
  PympzObject *new;

  assert(PyInt_Check(i));

  if (!(new = Pympz_new()))
    return NULL;
  mpz_set_si(new->z, PyInt_AsLong(i));
  return new;
}

static PympzObject *
long2mpz(PyObject * obj)
{
  PympzObject *new;

  mpz_t digit;
  int len;
  int i;
  int negative;
  PyLongObject *l;
  
  assert(PyLong_Check(obj));

  if (!(new = Pympz_new()))
    return NULL;
  mpz_set_si(new->z, 0);
  mpz_init(digit);
  
  l  = (PyLongObject *) obj;
  if (l->ob_size < 0)
    {
      len = - l->ob_size;
      negative = 1;
    }
  else
    {
      negative = 0;
      len = l->ob_size;
    }
  for (i=0; i<len; i++)
    {
      mpz_set_ui(digit, l->ob_digit[i]);
      mpz_mul_2exp(digit, digit, i * SHIFT);
      mpz_ior(new->z, new->z, digit);
    }
  if (negative)
    mpz_neg(new->z, new->z);
  mpz_clear(digit);
  return new;
}

static PympzObject *
str2mpz(PyObject *s)
{
  /* Least significant octet first */
  
  PympzObject *new;
  unsigned char *cp;
  int len;
  int i;
  mpz_t digit;

  assert(PyString_Check(s));

  if (! (new = Pympz_new()))
    return NULL;

  len = PyString_Size(s);
  cp = PyString_AsString(s);

  mpz_set_si(new->z, 0);
  mpz_init(digit);

  for (i=0; i<len; i++)
    {
      mpz_set_ui(digit, cp[i]);
      mpz_mul_2exp(digit, digit, i * 8);
      mpz_ior(new->z, new->z, digit);
    }
  return new;
}

static PyObject *
mpz2long(PympzObject *x)
{
  int negative;
  int size;
  int i;
  PyLongObject *new;
  mpz_t temp;

  /* Assume gmp uses larger limbs than the builtin longs */
  assert(mp_bits_per_limb >= SHIFT);
  
  size = (mpz_sizeinbase(x->z, 2) + SHIFT - 1) / SHIFT;
  if (!(new = _PyLong_New(size)))
    return NULL;
  mpz_init(temp);
  if (mpz_cmp_ui(x->z, 0) < 0)
    {
      negative = 1;
      mpz_neg(temp, x->z);
    }
  else
    {
      negative = 0;
      mpz_set(temp, x->z);
    }
  for (i=0; i<size; i++)
    {
      new->ob_digit[i] = mpz_get_ui(temp) & MASK;
      mpz_fdiv_q_2exp(temp, temp, SHIFT);
    }
  assert(mpz_cmp_ui(temp, 0) == 0);
  mpz_clear(temp);

  /* long_normalize() is file-static */
  /* longobjp = long_normalize(longobjp); */
  i = size;
  while ( (i>0) && (new->ob_digit[i-1] == 0))
    i--;
  new->ob_size = i+1;

  if (negative)
    new->ob_size = - new->ob_size;
  return (PyObject *) new;
}

static PyObject *
mpz2int(PympzObject *x)
{
  long res;

  res = mpz_get_si(x->z);
  if (mpz_cmp_si(x->z, res) != 0)
    {
      PyErr_SetString(PyExc_ValueError, "mpz.int() too large integer");
      return NULL;
    }
  return PyInt_FromLong(res);
}

static PyObject *
mpz2float(PympzObject *x)
{
  int negative;
  mpz_t temp;
  int e;
  int bits;
  int i;
  int weight;
  double res;
  
  mpz_init(temp);
  if (mpz_cmp_ui(x->z, 0) < 0)
    {
      negative = 1;
      mpz_neg(temp, x->z);
    }
  else
    {
      negative = 0;
      mpz_set(temp, x->z);
    }
  bits = mpz_sizeinbase(temp, 2);
  if (bits > double_mantissa)
    {
      /* Reduce precision */
      e = bits - double_mantissa;
      mpz_fdiv_q_2exp(temp, temp, e);
    }
  else
    e = 0;
  i = mpz_sizeinbase(temp, 2);
  for (res = 0.0, weight = 0;
       i > 0;
       i -= mp_bits_per_limb)
    {
      res += ldexp( (double) mpz_get_ui(temp), weight);
      weight += mp_bits_per_limb;
      mpz_fdiv_q_2exp(temp, temp, mp_bits_per_limb);
    }
  assert(mpz_cmp_ui(temp, 0) == 0);
  mpz_clear(temp);

  res = ldexp(res, e);
  if (negative)
    res = -res;
  return PyFloat_FromDouble(res);
}

static PyObject *
Pympz_binary(PympzObject *self, PyObject *args)
{
  mpz_t temp;
  int size;
  char *buffer;
  int i;
  PyObject *s;
  
  if (!PyArg_ParseTuple(args, ""))
    return NULL;

  assert(Pympz_Check( (PyObject *) self));

  if (mpz_cmp_ui(self->z, 0) < 0)
    {
      PyErr_SetString(PyExc_ValueError, "mpz.binary() arg must be >= 0");
      return NULL;
    }
  mpz_init_set(temp, self->z);

  size = (mpz_sizeinbase(temp, 2) + 7) / 8;

#ifdef USE_ALLOCA
  buffer = alloca(size);
#else
  if (!(buffer = malloc(size)))
    {
      mpz_clear(temp);
      PyErr_NoMemory();
      return NULL;
    }
#endif
  for (i = 0; i<size; i++)
    {
      buffer[i] = mpz_get_ui(temp) &0xff;
      mpz_fdiv_q_2exp(temp, temp, 8);
    }
  mpz_clear(temp);
  s = PyString_FromStringAndSize(buffer, size);
#ifndef USA_ALLOCA
  free(buffer);
#endif
  return s;
}

static PyObject *
Pympz_ascii(PympzObject *self, int base, int with_tag)
{
  PyObject *s;
  char *buffer;
  mpz_t temp;
  int minus;
  int is_long;
  char *p;
  
  assert(Pympz_Check( (PyObject *) self));
  assert( (base >= 2) && (base <= 36));

  mpz_init(temp);
  if (mpz_cmp_ui(self->z, 0) < 0)
    {
      minus = 1;
      mpz_neg(temp, self->z);
    }
  else
    {
      minus = 0;
      mpz_set(temp, self->z);
    }
  
  /* Allocate extra space for:
   *
   * minus sign and trailing NULL byte   (2)
   * 'gmp.mpz()' tag                     (9)
   * '0x' prefix                         (2)
   * 'L' suffix                          (1)
   *                                    -----
   *                                     14
   */
#ifdef USE_ALLOCA
  buffer = alloca(mpz_sizeinbase(self->z, base) + 16);
#else
  if (!(buffer = malloc(mpz_sizeinbase(self->z, base) + 16)))
    {
      mpz_clear(temp);
      PyErr_NoMemory();
      return NULL;
    }
#endif
  is_long = (mpz_cmp_si(self->z, mpz_get_si(self->z)) != 0);
  p = buffer;
  if (with_tag)
    p += sprintf(p, "gmp.mpz(");
  if (minus)
    *(p++) = '-';
  if (base == 8)
    *(p++) = '0';
  else if (base == 16)
    p += sprintf(p, "0x");
  
  mpz_get_str(p, base, temp); /* Doesn't return number of characters */
  p = buffer + strlen(buffer); /* Address of NULL byte */
  if (is_long)
    *(p++) = 'L';
  if (with_tag)
    *(p++) = ')';
  s = PyString_FromStringAndSize(buffer, p - buffer);
#ifndef USE_ALLOCA
  free(buffer);
#endif
  return s;
}

static PyObject *
mpz2str(PympzObject *self)
{
  return Pympz_ascii(self, 10, 1);
}

int
Pympz_convert_arg(PyObject *arg, PyObject **ptr)
{
  if (options.debug)
    fprintf(stderr, "mpz_convert_arg: %p\n", arg);
  
  if (Pympz_Check(arg))
    {
      Py_INCREF(arg);
      *ptr = arg;
      return 1;
    }
  else if (PyInt_Check(arg))
    {
      *ptr = (PyObject *) int2mpz(arg);
      return 1;
    }
  else if (PyLong_Check(arg))
    {
      *ptr = (PyObject *) long2mpz(arg);
      return 1;
    }
  else
    {
      PyErr_SetString(PyExc_TypeError,
		      "argument can not be converted to mpz");
      return 0;
    }
}

/* CONSTRUCTOR */
static PyObject *
Pygmp_mpz(PyObject *self, PyObject *args)
{
  PympzObject *new;
  PyObject *obj;

  if (options.debug)
    fputs("Pygmp_mpz() called...\n", stderr);
  if (!PyArg_ParseTuple(args, "O", &obj))
    return NULL;
  if (PyInt_Check(obj))
    new = int2mpz(obj);
  else if (PyLong_Check(obj))
    new = long2mpz(obj);
  else if PyString_Check(obj)
    new = str2mpz(obj);
  else if (Pympz_Check(obj))
    {
      Py_INCREF(obj);
      new = (PympzObject *) obj;
    }
  else
    {
      PyErr_SetString(PyExc_TypeError,
		      "gmp.mpz() expects integer, long, mpz or string argument");
      return NULL;
    }
  if (options.debug)
    {
      fputs("Pygmp_mpz: created mpz = ", stderr);
      mpz_out_str(stderr, 10, new->z);
      putc('\n', stderr);
    }
  return (PyObject *) new;
} /* Pygmp_mpz() */


/* ARITHMETIC */
#if 0
static PyObject *
Pympz_add(PympzObject *a, PympzObject *b)
{
  PympzObject *s;

  if (!(r = Pympz_new()))
    return NULL;
  Pympz_add(s->z), a->z, b->z);
  return (PyObject *) s;
}
#endif

#define MPZ_BINOP(NAME) \
static PyObject * \
Py##NAME(PympzObject *a, PympzObject *b) \
{ \
  PympzObject *r; \
  if (options.debug) fprintf(stderr, "Py" #NAME ": %p, %p", a, b); \
  if (!(r = Pympz_new())) return NULL; \
  NAME(r->z, a->z, b->z); \
  if (options.debug) fprintf(stderr, "Py" #NAME "-> %p", r); \
  return (PyObject *) r; \
}

MPZ_BINOP(mpz_add)
MPZ_BINOP(mpz_sub)
MPZ_BINOP(mpz_mul)

#define MPZ_DIVOP(NAME, OP) \
static PyObject * \
NAME(PympzObject *a, PympzObject *b) \
{ \
  PympzObject *r; \
  if (options.debug) fprintf(stderr, #NAME ": %p, %p", a, b); \
  if (mpz_cmp_ui(b->z, 0) == 0) { \
    PyErr_SetString(PyExc_ZeroDivisionError, #NAME "by zero"); \
    return NULL; } \
  if (!(r = Pympz_new())) return NULL; \
  OP(r->z, a->z, b->z); \
  if (options.debug) fprintf(stderr, #NAME "-> %p", r); \
  return (PyObject *) r; \
}

MPZ_DIVOP(Pympz_div, mpz_fdiv_q)
MPZ_DIVOP(Pympz_rem, mpz_fdiv_r)

static PyObject *
Pympz_divmod(PympzObject *a, PympzObject *b) 
{ 
  PympzObject *q, *r; 
  if (options.debug)
    fprintf(stderr, "Pympz_divmod: %p, %p", a, b); 
  if (mpz_cmp_ui(b->z, 0) == 0)
    { 
      PyErr_SetString(PyExc_ZeroDivisionError, "mpz.divmod by zero"); 
      return NULL;
    }
  if (!(q = Pympz_new()))
    return NULL;
  if (!(r = Pympz_new()))
    {
      Pympz_dealloc(q);
      return NULL;
    }
  mpz_fdiv_qr(q->z, r->z, a->z, b->z);
  if (options.debug)
    fprintf(stderr, "Pympz_divmod -> %p, %p", q, r); 
  return Py_BuildValue("(OO)", q, r);
}

#define MPZ_MONOP(NAME) \
static PyObject * \
Py##NAME(PympzObject *x) \
{ \
  PympzObject *r; \
  if (options.debug) fprintf(stderr, "Py" #NAME ": %p", x); \
  if (!(r = Pympz_new())) return NULL; \
  NAME(r->z, x->z); \
  if (options.debug) fprintf(stderr, "Py" #NAME "-> %p", r); \
  return (PyObject *) r; \
}

MPZ_MONOP(mpz_abs)
MPZ_MONOP(mpz_neg)

static PyObject *
Pympz_pos(PympzObject *x) { return (PyObject *) x; }

static PyObject *
Pympz_pow(PympzObject *b,PympzObject *e, PympzObject *m)
{
  PympzObject *r;

  if (options.debug)
    fprintf(stderr, "Pympz_pow: %p, %p, %p\n", b, e, m);
  
  if (mpz_cmp_ui(e->z, 0) < 0)
    {
      PyErr_SetString(PyExc_ValueError, "mpz.pow with negative power");
      return NULL;
    }
  if ( (PyObject *) m == Py_None)
    {
      /* When no modulo is present, the exponent must fit in C long */
      unsigned long el;
      el = mpz_get_ui(e->z);
      if (mpz_cmp_ui(e->z, el) != 0)
	{
	  PyErr_SetString(PyExc_ValueError, "mpz.pow outrageous exponent");
	  return NULL;
	}
      if (!(r = Pympz_new()))
	return NULL;
      mpz_pow_ui(r->z, b->z, el);
      if (options.debug)
	fprintf(stderr, "Pympz_pow (ui) -> %p\n", r);
      return (PyObject *) r;
    }
  else
    { /* Modulo exponentiation */
      if (mpz_cmp_ui(m->z, 0) == 0)
	{
	  PyErr_SetString(PyExc_ValueError, "mpz.pow divide by zero");
	    return NULL;
	}
      /* What happens with negative modulo? */
      if (!(r = Pympz_new()))
	return NULL;
      mpz_powm(r->z, b->z, e->z, m->z);
      if (options.debug)
	fprintf(stderr, "Pympz_pow -> %p\n", r);
      return (PyObject *) r;
    }
}


/* COMPARING */
static int
Pympz_cmp(PympzObject *a, PympzObject *b)
{
  int res = mpz_cmp(a->z, b->z);
  return res ? ( res < 0 ? -1: 1) : 0;
}

static int
Pympz_nonzero(PympzObject *x)
{
  return (mpz_cmp_ui(x->z, 0) != 0);
}


/* BIT OPERATIONS */
MPZ_MONOP(mpz_com)
MPZ_BINOP(mpz_and)
MPZ_BINOP(mpz_ior)

static PyObject *
Pympz_xor(PympzObject *a, PympzObject *b)
{ /* Non-optimal */
  mpz_t temp;
  PympzObject *r;

  if (!(r = Pympz_new()))
    return NULL;
  mpz_init(temp);
  mpz_and(temp, a->z, b->z);
  mpz_com(temp, temp);
  mpz_ior(r->z, a->z, b->z);
  mpz_and(r->z, r->z, temp);
  mpz_clear(temp);
  return (PyObject *) r;
}
  
#define MPZ_SHIFTOP(NAME, OP) \
static PyObject * \
NAME(PympzObject *a, PympzObject *b) \
{ \
    PympzObject *r; \
    unsigned long count; \
    if (mpz_cmp_ui(b->z, 0) < 0) { \
       PyErr_SetString(PyExc_ValueError, #NAME " negative shift count"); \
       return NULL; } \
    count = mpz_get_ui(b->z); \
    if (mpz_cmp_ui(b->z, count) != 0) { \
       PyErr_SetString(PyExc_ValueError, "mpz.pow outrageous shift count"); \
       return NULL; } \
    if (! (r = Pympz_new())) return NULL; \
    OP(r->z, a->z, count); \
    return (PyObject *) r; \
}

MPZ_SHIFTOP(Pympz_rshift, mpz_fdiv_q_2exp)
MPZ_SHIFTOP(Pympz_lshift, mpz_mul_2exp)

static PyObject *
Pympz_oct(PympzObject *self)
{
  return Pympz_ascii(self, 8, 1);
}

static PyObject *
Pympz_hex(PympzObject *self)
{
  return Pympz_ascii(self, 16, 1);
}

static int
Pympz_coerce(PyObject **pv, PyObject **pw)
{
  PyObject *z;

  if (options.debug)
    fprintf(stderr, "Pympz.coerce(%p, %p) called...\n", *pv, *pw);
  assert(Pympz_Check(*pv));

  /* Convert all types except floats to mpz */
  if (PyFloat_Check(*pw))
    {
      if (options.debug)
	fprintf(stderr, "Pympz.coerce(): float \n");
      if (!(z = mpz2float( (PympzObject *) *pv)))
	return -1;
      *pv = z;
      Py_INCREF(*pw);
      return 0;
    }
  if (PyInt_Check(*pw))
    {
      if (!(z = (PyObject *) int2mpz(*pw)))
	return -1;
      Py_INCREF(*pv);
      *pw = z;
      return 0;
    }
  if (PyLong_Check(*pw))
    {
      if (!(z = (PyObject *) long2mpz(*pw)))
	return -1;
      Py_INCREF(*pv);
      *pw = z;
      return 0;
    }
  else
    {
      PyErr_SetString(PyExc_TypeError, "coercion to gmp.mpz type failed");
      return -1;
    }
}

static long
Pympz_hash(PympzObject *self)
{
  int i;
  int size;
  long hash;
  mpz_t temp;
  
  size = (mpz_sizeinbase(self->z, 2) + mp_bits_per_limb - 1) / mp_bits_per_limb;
  mpz_init_set(temp, self->z);
  mpz_abs(temp, temp);

  for (i=0, hash = 0; i<size; i++, mpz_fdiv_q_2exp(temp, temp, mp_bits_per_limb))
    hash += mpz_get_ui(temp);
  mpz_clear(temp);
  return hash;
}

static PyObject *
Pympz_getattr(PympzObject *self, char *name)
{
  return Py_FindMethod(Pympz_methods, (PyObject *) self, name);
}


/* Miscellaneous GMP functions */

static PyObject *
Pygmp_gcd(PyObject *self, PyObject *args)
{
  PympzObject *a, *b, *c;

  if (!PyArg_ParseTuple(args, "O&O&", Pympz_convert_arg, &a, Pympz_convert_arg, &b))
    return NULL;

  assert(Pympz_Check( (PyObject *) a));
  assert(Pympz_Check( (PyObject *) b));

  if (!(c = Pympz_new()))
    {
      Py_DECREF(a); Py_DECREF(b);
      return NULL;
    }
  mpz_gcd(c->z, a->z, b->z);
  Py_DECREF(a); Py_DECREF(b);
  return (PyObject *) c;
}

static PyObject *
Pygmp_gcdext(PyObject *self, PyObject *args)
{
  PympzObject *a, *b, *g, *s, *t;
  PyObject *res;
  if (!PyArg_ParseTuple(args, "O&O&",
			Pympz_convert_arg, &a,
			Pympz_convert_arg, &b))
    return NULL;

  assert(Pympz_Check( (PyObject *) a));
  assert(Pympz_Check( (PyObject *) b));

  if (!(g = Pympz_new()))
    {
      Py_DECREF(a); Py_DECREF(b);
      return NULL;
    }
  if (!(s = Pympz_new()))
    {
      Py_DECREF(a); Py_DECREF(b);
      Py_DECREF(g);
      return NULL;
    }
  if (!(t = Pympz_new()))
    {
      Py_DECREF(a); Py_DECREF(b);
      Py_DECREF(g);
      Py_DECREF(s);
      return NULL;
    }
  mpz_gcdext(g->z, s->z, t->z, a->z, b->z);
  Py_DECREF(a); Py_DECREF(b);
  res = Py_BuildValue("OOO", g, s, t); /* Increments refcounts */
  Py_DECREF(g); Py_DECREF(s); Py_DECREF(t);
  
  return (PyObject *) res;
}

static PyObject *
Pygmp_divm(PyObject *self, PyObject *args)
{
  PympzObject *num, *den, *mod, *res;
  
  
  if (! (res = Pympz_new()))
    return NULL;
  if (!PyArg_ParseTuple(args, "O&O&O&",
			Pympz_convert_arg, &num,
			Pympz_convert_arg, &den,
			Pympz_convert_arg, &mod))
    {
      Py_DECREF(res);
      return NULL; /* I hope no more DECREF:s is needed here (?) */
    }
  if (mpz_invert(res->z, den->z, mod->z))
    { /* inverse exists */
      mpz_mul(res->z, res->z, num->z);
      return (PyObject *) res;
    }
  else
    {
      Py_DECREF(res);
      PyErr_SetString(PyExc_ZeroDivisionError, "not invertible");
      return NULL;
    }
}

static PyObject *
Pygmp_fac(PyObject *self, PyObject *args)
{
  PympzObject *fac;
  long n;
  
  if (!PyArg_ParseTuple(args, "l", &n))
    return NULL;

  if (n <= 0)
    {
      PyErr_SetString(PyExc_ValueError, "factorial of negative number");
      return NULL;
    }
  if (! (fac = Pympz_new()))
    return NULL;
  mpz_fac_ui(fac->z, n);

  return (PyObject *) fac;
}

static PyObject *
Pympz_sqrt(PympzObject *self, PyObject *args)
{
  PympzObject *root;

  if (self)   /* Called as method */
    {
      if (!PyArg_ParseTuple(args, ""))
	return NULL;
      Py_INCREF(self);
    }
  else        /* Called as function */
    {
      if (!PyArg_ParseTuple(args, "O&", Pympz_convert_arg, &self))
	return NULL;
    }

  if (! (root = Pympz_new()))
    {
      Py_DECREF(self);
      return NULL;
    }
  mpz_sqrt(root->z, self->z);
  return (PyObject *) root;
}

static PyObject *
Pympz_sqrtrem(PympzObject *self, PyObject *args)
{
  PympzObject *root, *rem;
  PyObject *res;

  if (self)   /* Called as method */
    {
      if (!PyArg_ParseTuple(args, ""))
	return NULL;
      Py_INCREF(self);
    }
  else        /* Called as function */
    {
      if (!PyArg_ParseTuple(args, "O&", Pympz_convert_arg, &self))
	return NULL;
    }

  if (! (root = Pympz_new()))
    {
      Py_DECREF(self);
      return NULL;
    }
  if (! (rem = Pympz_new()))
    {
      Py_DECREF(root); Py_DECREF(self);
      return NULL;
    }
  mpz_sqrtrem(root->z, rem->z, self->z);

  res = Py_BuildValue("OO", root, rem);
  Py_DECREF(self); Py_DECREF(root); Py_DECREF(rem); 
  return res;
}

static PyObject *
Pympz_is_square(PympzObject *self, PyObject *args)
{
  PyObject *res;
  
  if (self)   /* Called as method */
    {
      if (!PyArg_ParseTuple(args, ""))
	return NULL;
      Py_INCREF(self);
    }
  else        /* Called as function */
    {
      if (!PyArg_ParseTuple(args, "O&", Pympz_convert_arg, &self))
	return NULL;
    }
  res = Py_BuildValue("i", mpz_perfect_square_p(self->z));
  Py_DECREF(self);
  return res;
}

static PyObject *
Pympz_is_prime(PympzObject *self, PyObject *args)
{
  PyObject *res;
  int reps;

  reps = 25;  /* Default value */
  if (self)       /* method */
    {
      if (!PyArg_ParseTuple(args, "|i", &reps))
	return NULL;
      Py_INCREF(self);
    } 
  else            /* function */
    {
      if (!PyArg_ParseTuple(args, "O&|i", Pympz_convert_arg, &self, &reps))
	return NULL;
    }
  if (reps <= 0)
    {
      PyErr_SetString(PyExc_ValueError,
		      "repetition count for is_prime must be positive");
      return NULL;
    }
  res = Py_BuildValue("i", mpz_probab_prime_p(self->z, reps));
  Py_DECREF(self);
  return res;
}

static PyNumberMethods mpz_number_methods =
{
  (binaryfunc) Pympz_add,
  (binaryfunc) Pympz_sub,
  (binaryfunc) Pympz_mul,
  (binaryfunc) Pympz_div,
  (binaryfunc) Pympz_rem,
  (binaryfunc) Pympz_divmod,
  (ternaryfunc) Pympz_pow,
  (unaryfunc) Pympz_neg,
  (unaryfunc) Pympz_pos,  /* nb_positive */
  (unaryfunc) Pympz_abs,
  (inquiry) Pympz_nonzero,
  (unaryfunc) Pympz_com,
  (binaryfunc) Pympz_lshift,
  (binaryfunc) Pympz_rshift,
  (binaryfunc) Pympz_and,
  (binaryfunc) Pympz_xor,
  (binaryfunc) Pympz_ior,
  (coercion) Pympz_coerce,
  (unaryfunc) mpz2int,
  (unaryfunc) mpz2long,
  (unaryfunc) mpz2float,
  (unaryfunc) Pympz_oct,
  (unaryfunc) Pympz_hex
};

static PyMethodDef Pygmp_methods [] =
{
  { "set_debug", Pygmp_set_debug, 1 },
  { "mpz", Pygmp_mpz, 1 },
  { "gcd", Pygmp_gcd, 1 },
  { "gcdext", Pygmp_gcdext, 1 },
  { "divm", Pygmp_divm, 1 },
  { "fac", Pygmp_fac, 1},

  { "sqrt", Pympz_sqrt, 1 },
  { "sqrtrem", Pympz_sqrtrem, 1 },
  { "is_square", Pympz_is_square, 1},
  { "is_prime", Pympz_is_prime, 1},
  { "binary", Pympz_binary, 1 },
  { NULL, NULL, 1}
};

statichere PyMethodDef Pympz_methods [] =
{
  { "sqrt", Pympz_sqrt, 1 },
  { "sqrtrem", Pympz_sqrtrem, 1 },
  { "is_square", Pympz_is_square, 1},
  { "is_prime", Pympz_is_prime, 1},
  { "binary", Pympz_binary, 1 },
  { NULL, NULL, 1 }
};

statichere PyTypeObject Pympz_Type =
{
  PyObject_HEAD_INIT(&PyType_Type)
  0,                            /* ob_size */
  "mpz",                        /* tp_name */
  sizeof(PympzObject),		/* tp_basicsize */
  0,				/* tp_itemsize */
  /* methods */
  (destructor) Pympz_dealloc,	/* tp_dealloc */
  0,				/* tp_print */
  (getattrfunc) Pympz_getattr,	/* tp_getattr */
  0,			        /* tp_setattr */
  (cmpfunc) Pympz_cmp,		/* tp_compare */
  (reprfunc) mpz2str,		/* tp_repr */
  &mpz_number_methods,          /* tp_as_number */
  0,                            /* tp_as_sequense */
  0,                            /* tp_as_mapping */
  (hashfunc) Pympz_hash,        /* tp_hash */
  0,                            /* tp_call */
  (reprfunc) mpz2str,           /* tp_str */
  0, 0, 0, 0,                   /* reserved */
  "GNU Multi Precision signed integer",
};

static void *
gmp_allocate(size_t size)
{
  void *res;

  if (options.debug)
    fprintf(stderr, "mp_allocate( %ld )\n", size);
  if ( ! (res = malloc(size)) )
    Py_FatalError("mp_allocate failure");

  if (options.debug)
    fprintf(stderr, "mp_allocate ->%8p\n", res);

  /* MP_SET_TEST(res,alloc_size); */
	
  return res;
} /* mp_allocate() */


static void *
gmp_reallocate(void *ptr, size_t old_size, size_t new_size)
{
  void *res;

  if (options.debug)
    fprintf(stderr, "mp_reallocate: old address %8p, old size %ld\n",
	    ptr, old_size);

  /* MP_DO_TEST(ptr, old_size); */
	
  if ( !(res = realloc(ptr, new_size)) )
    Py_FatalError("mp_reallocate failure");
  
  if (options.debug)
    fprintf(stderr, "mp_reallocate: new address %8p, new size %ld\n",
	    res, new_size);

  /* MP_SET_TEST(res, new_size); */

  return res;
} /* mp_reallocate() */

static void
gmp_free( void *ptr, size_t size)
{
  if (options.debug)
    fprintf(stderr, "mp_free      : old address %8p, old size %ld\n",
	    ptr, size);
  
  /* MP_DO_TEST(ptr, size); */
  free(ptr);
} /* mp_free() */


/* Find out how many limbs are significant in a double */

static unsigned get_precision(void)
{
  int bit;
  double eps;
  
  for (bit = 0, eps = 1.0; 1.0 != (1.0 + eps); bit++)
    eps /= 2;
  return bit;
}
  
void
initgmp(void)
{
  if (options.debug)
    fputs( "initgmp() called...\n", stderr );

  mp_set_memory_functions( gmp_allocate, gmp_reallocate, gmp_free );

  double_mantissa = get_precision();
  gmp_module = Py_InitModule("gmp", Pygmp_methods);
}