/*
 * FREE GMP - a public domain implementation of a subset of the 
 *           gmp library
 *
 * I hearby place the file in the public domain.
 *
 * Do whatever you want with this code. Change it. Sell it. Claim you
 *  wrote it. 
 * Bugs, complaints, flames, rants: please send email to 
 *    Mark Henderson <markh@wimsey.bc.ca>
 * I'm already aware that fgmp is considerably slower than gmp
 *
 * CREDITS:
 *  Paul Rouse <par@r-cube.demon.co.uk> - generic bug fixes, mpz_sqrt and
 *    mpz_sqrtrem, and modifications to get fgmp to compile on a system 
 *    with int and long of different sizes (specifically MSDOS,286 compiler)
 *  Also see the file "notes" included with the fgmp distribution, for
 *    more credits.
 *
 * VERSION 1.0 - beta 5
 */

#include "gmp.h"
#ifndef NULL
#define NULL ((void *)0)
#endif

#define iabs(x) ((x>0) ? (x) : (-x))
#define imax(x,y) ((x>y)?x:y)
#define LONGBITS (sizeof(mp_limb) * 8)
#define DIGITBITS (LONGBITS - 2)
#define HALFDIGITBITS ((LONGBITS -2)/2)
#ifndef init
#define init(g) { g = (MP_INT *)malloc(sizeof(MP_INT));  mpz_init(g); }
#endif
#ifndef clear
#define clear(g) { mpz_clear(g); free(g); }
#endif
#ifndef even
#define even(a) (!((a)->p[0] & 1))
#endif
#ifndef odd
#define odd(a) (((a)->p[0] & 1))
#endif


#ifndef B64
/* 
 * The values below are for 32 bit machines (i.e. machines with a 
 *  32 bit long type)
 * You'll need to change them, if you're using something else
 * If DIGITBITS is odd, see the comment at the top of mpz_sqrtrem
 */
#define LMAX 0x3fffffffL
#define LC   0xc0000000L
#define OVMASK 0x2
#define CMASK (LMAX+1)
#define FC ((double)CMASK)
#define HLMAX 0x7fffL
#define HCMASK (HLMAX + 1)
#define HIGH(x) (((x) & 0x3fff8000L) >> 15)
#define LOW(x)  ((x) & 0x7fffL)
#else

/* 64 bit long type */
#define LMAX 0x3fffffffffffffffL
#define LC 0xc000000000000000L
#define OVMASK 0x2
#define CMASK (LMAX+1)
#define FC ((double)CMASK)
#define HLMAX 0x7fffffffL
#define HCMASK (HLMAX + 1)
#define HIGH(x) (((x) & 0x3fffffff80000000L) >> 31)
#define LOW(x) ((x) & 0x7fffffffL)

#endif

#define hd(x,i)  (((i)>=2*(x->sz))? 0:(((i)%2) ? HIGH((x)->p[(i)/2]) \
    : LOW((x)->p[(i)/2])))
#define dg(x,i) ((i < x->sz) ? (x->p)[i] : 0)

#ifdef __GNUC__
#define INLINE inline
#else
#define INLINE
#endif

#define lowdigit(x) (((x)->p)[0])

struct is {
    mp_limb v;
    struct is *next;
};


INLINE static void push(i,sp)
mp_limb i;struct is **sp;
{
    struct is *tmp;
    tmp = *sp;
    *sp = (struct is *) malloc(sizeof(struct is));
    (*sp)->v = i;
    (*sp)->next=tmp;
}

INLINE static mp_limb pop(sp)
struct is **sp;
{
    struct is *tmp;
    mp_limb i;
    if (!(*sp))
        return (-1);
    tmp = *sp;
    *sp = (*sp)->next;
    i = tmp->v;
    tmp->v = 0;
    free(tmp);
    return i;
}

void fatal(s)
char *s;
{
    fprintf(stderr,"%s\n",s);
    exit(123);
}

void mpz_init(s)
MP_INT *s;
{
    s->p = (mp_limb *)malloc(sizeof(mp_limb)*2);
#ifdef DEBUG
    printf("malloc: 8 bytes, %08lx\n", (long)(s->p));
#endif
    if (!(s->p))
        fatal("mpz_init: cannot allocate memory");
    (s->p)[0] = 0;
    (s->p)[1] = 0;
    s->sn=0;
    s->sz=2;
}

void mpz_init_set(s,t)
MP_INT *s,*t;
{
    int i;
    s->p = (mp_limb *)malloc(sizeof(mp_limb) * t->sz);
    if (!(s->p))
        fatal("mpz_init: cannot allocate memory");
    for (i=0;i < t->sz ; i++)
        (s->p)[i] = (t->p)[i];

    s->sn = t->sn;
    s->sz = t->sz;
}

void mpz_init_set_ui(s,v)
MP_INT *s;
unsigned long v;
{
    s->p = (mp_limb *)malloc(sizeof(mp_limb)*2);
    if (!(s->p))
        fatal("mpz_init: cannot allocate memory");
    s->p[0] = (v & LMAX);
    s->p[1] = (v & LC) >> DIGITBITS;
    if (v) 
        s->sn = 1;
    else
        s->sn = 0;
    s->sz = 2;
}

void mpz_init_set_si(y,v)
MP_INT *y;
long v;
{
    y->p = (mp_limb *)malloc(sizeof(mp_limb)*2);
    if (!(y->p))
        fatal("mpz_init: cannot allocate memory");
    if (v < 0) {
        y->sn = -1;
        y->p[0] = (-v) & LMAX;
        y->p[1] = ((-v) & LC) >> DIGITBITS;
    }
    else if (v > 0) {
        y->sn = 1;
        y->p[0] = v & LMAX;
        y->p[1] = (v & LC) >> DIGITBITS;
    }
    else {
        y->sn=0;
        y->p[0] = 0;
        y->p[1] = 0;
    }
    y -> sz = 2;
}



void mpz_clear(s)
MP_INT *s;
{
#ifdef DEBUG
    printf("free: %08lx\n", (long)(s->p));
#endif
    if (s->p)
        free(s->p);
    s->p=NULL;
    s->sn=0;
    s->sz=0;
}

/* number of leading zero bits in digit */
static int lzb(a)
mp_limb a;
{
    mp_limb i; int j;
    j = 0;
    for (i = ((mp_limb)1 << (DIGITBITS-1)); i && !(a&i) ; j++,i>>=1) 
        ;
    return j;
}

void _mpz_realloc(x,size)
MP_INT *x; mp_size size;
{
    if (size > 1 && x->sz < size) {
        int i;
#ifdef DEBUG
    printf("realloc %08lx to size = %ld ", (long)(x->p),(long)(size));
#endif

        if (x->p) 
            x->p=(mp_limb *)realloc(x->p,size * sizeof(mp_limb));
        else
            x->p=(mp_limb *)malloc(size * sizeof(mp_limb));
#ifdef DEBUG
    printf("becomes %08lx\n", (long)(x->p));
#endif

        if (!(x->p))
            fatal("_mpz_realloc: cannot allocate memory");
        for (i=x->sz;i<size;i++)
            (x->p)[i] = 0;
        x->sz = size;
    }
}

void dgset(x,i,n)
MP_INT *x; unsigned int i; mp_limb n;
{
    if (n) {
        if (i >= x->sz) {
            _mpz_realloc(x,(mp_size)(i+1));
            x->sz = i+1;
        }
        (x->p)[i] = n;
    }
}

static unsigned int digits(x)
MP_INT *x;
{
    int i;
    for (i = (x->sz) - 1; i>=0 && (x->p)[i] == 0 ; i--) 
        ;
    return i+1;
}
        
/* y = x */
void mpz_set(y,x)
MP_INT *y; MP_INT *x; 
{
    unsigned int i,k = x->sz;
    if (y->sz < k) {
        k=digits(x);
        _mpz_realloc(y,(mp_size)k);
    }
    if (y->sz > x->sz) {
        mpz_clear(y); mpz_init(y); _mpz_realloc(y,(mp_size)(x->sz));
    }

    for (i=0;i < k; i++)
        (y->p)[i] = (x->p)[i];

    for (;i<y->sz;i++)
        (y->p)[i] = 0;

    y->sn = x->sn;
}

void mpz_set_ui(y,v)
MP_INT *y; unsigned long v; {
    int i;
    for (i=1;i<y->sz;i++)
        y->p[i] = 0;
    y->p[0] = (v & LMAX);
    y->p[1] = (v & LC) >> (DIGITBITS);
    if (v)
        y->sn=1;
    else
        y->sn=0;
}

unsigned long mpz_get_ui(y) 
MP_INT *y; {
    return (y->p[0] | (y->p[1] << DIGITBITS));
}

long mpz_get_si(y)
MP_INT *y; {
    if (y->sn == 0)
        return 0;
    else
        return (y->sn * (y->p[0] | (y->p[1] & 1) << DIGITBITS));
}
        
void mpz_set_si(y,v)
MP_INT *y; long v; {
    int i;
    for (i=1;i<y->sz;i++)
        y->p[i] = 0;
    if (v < 0) {
        y->sn = -1;
        y->p[0] = (-v) & LMAX;
        y->p[1] = ((-v) & LC) >> DIGITBITS;
    }
    else if (v > 0) {
        y->sn = 1;
        y->p[0] = v & LMAX;
        y->p[1] = (v & LC) >> DIGITBITS;
    }
    else {
        y->sn=0;
        y->p[0] = 0;
        y->p[1] = 0;
    }
}

/* z = x + y, without regard for sign */
static void uadd(z,x,y)
MP_INT *z; MP_INT *x; MP_INT *y;
{
    mp_limb c;
    int i;
    MP_INT *t;

    if (y->sz < x->sz) {
        t=x; x=y; y=t;
    }

    /* now y->sz >= x->sz */

    _mpz_realloc(z,(mp_size)((y->sz)+1));

    c=0;
    for (i=0;i<x->sz;i++) {
        if (( z->p[i] = y->p[i] + x->p[i] + c ) & CMASK) {
            c=1; 
            (z->p[i]) &=LMAX;
        }
        else 
            c=0;
    }
    for (;i<y->sz;i++) {
        if ((z->p[i] = (y->p[i] + c)) & CMASK)
            z->p[i]=0;
        else
            c=0;
    }
    (z->p)[y->sz]=c;
}

/* z = y - x, ignoring sign */
/* precondition: abs(y) >= abs(x) */
static void usub(z,y,x)
MP_INT *z; MP_INT *y; MP_INT *x; 
{
    int i;
    mp_limb b,m;
    _mpz_realloc(z,(mp_size)(y->sz));
    b=0;
    for (i=0;i<y->sz;i++) {
        m=((y->p)[i]-b)-dg(x,i);
        if (m < 0) {
            b = 1;
            m = LMAX + 1 + m;
        }
        else
            b = 0;
        z->p[i] = m;
    }
}

/* compare abs(x) and abs(y) */
static int ucmp(y,x)
MP_INT *y;MP_INT *x;
{
    int i;
    for (i=imax(x->sz,y->sz)-1;i>=0;i--) {
        if (dg(y,i) < dg(x,i))
            return (-1);
        else if (dg(y,i) > dg(x,i))
            return 1;
    }
    return 0;
}

int mpz_cmp(x,y)
MP_INT *x; MP_INT *y;
{
    int abscmp;
    if (x->sn < 0 && y->sn > 0)
        return (-1);
    if (x->sn > 0 && y->sn < 0)
        return 1;
    abscmp=ucmp(x,y);
    if (x->sn >=0 && y->sn >=0)
        return abscmp;
    if (x->sn <=0 && y->sn <=0)
        return (-abscmp);
}
    
/* is abs(x) == 0 */
static int uzero(x)
MP_INT *x;
{
    int i;
    for (i=0; i < x->sz; i++)
        if ((x->p)[i] != 0)
            return 0;
    return 1;
}

static void zero(x)
MP_INT *x;
{
    int i;
    x->sn=0;
    for (i=0;i<x->sz;i++)
        (x->p)[i] = 0;
}


/* w = u * v */
void mpz_mul(ww,u,v)
MP_INT *ww;MP_INT *u; MP_INT *v; {
    int i,j;
    mp_limb t0,t1,t2,t3;
    mp_limb cc;
    MP_INT *w;

    w = (MP_INT *)malloc(sizeof(MP_INT));
    mpz_init(w);
    _mpz_realloc(w,(mp_size)(u->sz + v->sz));
    for (j=0; j < 2*u->sz; j++) {
        cc = (mp_limb)0;
        t3 = hd(u,j);
            for (i=0; i < 2*v->sz; i++) {
                t0 = t3 * hd(v,i);
                t1 = HIGH(t0); t0 = LOW(t0);
                if ((i+j)%2) 
                    t2 = HIGH(w->p[(i+j)/2]);
                else
                    t2 = LOW(w->p[(i+j)/2]);
                t2 += cc;
                if (t2 & HCMASK) {
                    cc = 1; t2&=HLMAX;
                }
                else
                    cc = 0;
                t2 += t0;
                if (t2 & HCMASK) {
                    cc++ ; t2&=HLMAX;
                }
                cc+=t1;
                if ((i+j)%2)
                    w->p[(i+j)/2] = LOW(w->p[(i+j)/2]) |
                        (t2 << HALFDIGITBITS);
                else
                    w->p[(i+j)/2] = (HIGH(w->p[(i+j)/2]) << HALFDIGITBITS) |
                        t2;
                    
            }
        if (cc) {
            if ((j+i)%2) 
                w->p[(i+j)/2] += cc << HALFDIGITBITS;
            else
                w->p[(i+j)/2] += cc;
        }
    }
    w->sn = (u->sn) * (v->sn);
    if (w != ww) {
        mpz_set(ww,w);
        mpz_clear(w);
        free(w);
    }
}

/* n must be < DIGITBITS */
static void urshift(c1,a,n)
MP_INT *c1;MP_INT *a;int n;
{
    mp_limb cc = 0;
    if (n >= DIGITBITS)
        fatal("urshift: n >= DIGITBITS");
    if (n == 0)
        mpz_set(c1,a);
    else {
        MP_INT c; int i;
        mp_limb rm = (((mp_limb)1<<n) - 1);
        mpz_init(&c); _mpz_realloc(&c,(mp_size)(a->sz));
        for (i=a->sz-1;i >= 0; i--) {
            c.p[i] = ((a->p[i] >> n) | cc) & LMAX;
            cc = (a->p[i] & rm) << (DIGITBITS - n);
        }
        mpz_set(c1,&c);
        mpz_clear(&c);
    }
}
    
/* n must be < DIGITBITS */
static void ulshift(c1,a,n)
MP_INT *c1;MP_INT *a;int n;
{
    mp_limb cc = 0;
    if (n >= DIGITBITS)
        fatal("ulshift: n >= DIGITBITS");
    if (n == 0)
        mpz_set(c1,a);
    else {
        MP_INT c; int i;
        mp_limb rm = (((mp_limb)1<<n) - 1) << (DIGITBITS -n);
        mpz_init(&c); _mpz_realloc(&c,(mp_size)(a->sz + 1));
        for (i=0;i < a->sz; i++) {
            c.p[i] = ((a->p[i] << n) | cc) & LMAX;
            cc = (a->p[i] & rm) >> (DIGITBITS -n);
        }
        c.p[i] = cc;
        mpz_set(c1,&c);
        mpz_clear(&c);
    }
}

void mpz_div_2exp(z,x,e) 
MP_INT *z; MP_INT *x; unsigned long e;
{
    short sn = x->sn;
    if (e==0)
        mpz_set(z,x);
    else {
        unsigned long i;
        long digs = (e / DIGITBITS);
        unsigned long bs = (e % (DIGITBITS));
        MP_INT y; mpz_init(&y);
        _mpz_realloc(&y,(mp_size)((x->sz) - digs));
        for (i=0; i < (x->sz - digs); i++)
            (y.p)[i] = (x->p)[i+digs];
        if (bs) {
            urshift(z,&y,bs);
        }
        else {
            mpz_set(z,&y);
        }
        if (uzero(z))
            z->sn = 0;
        else
            z->sn = sn;
        mpz_clear(&y);
    }
}

void mpz_mul_2exp(z,x,e) 
MP_INT *z; MP_INT *x; unsigned long e;
{
    short sn = x->sn;
    if (e==0)
        mpz_set(z,x);
    else {
        unsigned long i;
        long digs = (e / DIGITBITS);
        unsigned long bs = (e % (DIGITBITS));
        MP_INT y; mpz_init(&y);
        _mpz_realloc(&y,(mp_size)((x->sz)+digs));
        for (i=digs;i<((x->sz) + digs);i++)
            (y.p)[i] = (x->p)[i - digs];
        if (bs) {
            ulshift(z,&y,bs);
        }
        else {
            mpz_set(z,&y);
        }
        z->sn = sn;
        mpz_clear(&y);
    }
}

void mpz_mod_2exp(z,x,e) 
MP_INT *z; MP_INT *x; unsigned long e;
{
    short sn = x->sn;
    if (e==0)
        mpz_set_ui(z,0L);
    else {
        unsigned long i;
        MP_INT y;
        long digs = (e / DIGITBITS);
        unsigned long bs = (e % (DIGITBITS));
        if (digs > x->sz || (digs == (x->sz) && bs)) {
            mpz_set(z,x);
            return;
        }
            
        mpz_init(&y);
        if (bs)
            _mpz_realloc(&y,(mp_size)(digs+1));
        else
            _mpz_realloc(&y,(mp_size)digs);
        for (i=0; i<digs; i++)
            (y.p)[i] = (x->p)[i];
        if (bs) {
            (y.p)[digs] = ((x->p)[digs]) & (((mp_limb)1 << bs) - 1);
        }
        mpz_set(z,&y);
        if (uzero(z))
            z->sn = 0;
        else
            z->sn = sn;
        mpz_clear(&y);
    }
}
    

/* internal routine to compute x/y and x%y ignoring signs */
static void udiv(qq,rr,xx,yy)
MP_INT *qq; MP_INT *rr; MP_INT *xx; MP_INT *yy;
{
    MP_INT *q, *x, *y, *r;
    int ns,xd,yd,j,f,i,ccc;
    mp_limb zz,z,qhat,b,u,m;
    ccc=0;

    if (uzero(yy))
        return;
    q = (MP_INT *)malloc(sizeof(MP_INT));
    r = (MP_INT *)malloc(sizeof(MP_INT));
    x = (MP_INT *)malloc(sizeof(MP_INT)); y = (MP_INT *)malloc(sizeof(MP_INT));
    if (!x || !y || !q || !r)
        fatal("udiv: cannot allocate memory");
    mpz_init(q); mpz_init(x);mpz_init(y);mpz_init(r);
    _mpz_realloc(x,(mp_size)((xx->sz)+1));
    yd = digits(yy);
    ns = lzb(yy->p[yd-1]);
    ulshift(x,xx,ns);
    ulshift(y,yy,ns);
    xd = digits(x); 
    _mpz_realloc(q,(mp_size)xd);
    xd*=2; yd*=2;
    z = hd(y,yd-1);;
    for (j=(xd-yd);j>=0;j--) {
#ifdef DEBUG
    printf("y=");
    for (i=yd-1;i>=0;i--)
        printf(" %04lx", (long)hd(y,i));
    printf("\n");
    printf("x=");
    for (i=xd-1;i>=0;i--)
        printf(" %04lx", (long)hd(x,i));
    printf("\n");
    printf("z=%08lx\n",(long)z);
#endif
        
        if (z == LMAX)
            qhat = hd(x,j+yd);
        else {
            qhat = ((hd(x,j+yd)<< HALFDIGITBITS) + hd(x,j+yd-1)) / (z+1);
        }
#ifdef DEBUG
    printf("qhat=%08lx\n",(long)qhat);
    printf("hd=%04lx/%04lx\n",(long)hd(x,j+yd),(long)hd(x,j+yd-1));
#endif
        b = 0; zz=0;
        if (qhat) {
            for (i=0; i < yd; i++) {
                zz = qhat * hd(y,i);
                u = hd(x,i+j);
                u-=b;
                if (u<0) {
                    b=1; u+=HLMAX+1;
                }
                else
                    b=0;
                u-=LOW(zz);
                if (u < 0) {
                    b++;
                    u+=HLMAX+1;
                }
                b+=HIGH(zz);
                if ((i+j)%2) 
                    x->p[(i+j)/2] = LOW(x->p[(i+j)/2]) | (u << HALFDIGITBITS);
                else
                    x->p[(i+j)/2] = (HIGH(x->p[(i+j)/2]) << HALFDIGITBITS) | u;
            }
            if (b) {
                if ((j+i)%2) 
                    x->p[(i+j)/2] -= b << HALFDIGITBITS;
                else
                    x->p[(i+j)/2] -= b;
            }
        }
#ifdef DEBUG
        printf("x after sub=");
        for (i=xd-1;i>=0;i--)
            printf(" %04lx", (long)hd(x,i));
        printf("\n");
#endif
        for(;;zz++) {
            f=1;
            if (!hd(x,j+yd)) {
                for(i=yd-1; i>=0; i--) {
                    if (hd(x,j+i) > hd(y,i)) {
                        f=1;
                        break;
                    }
                    if (hd(x,j+i) < hd(y,i)) {
                        f=0;
                        break;
                    }
                }
            }
            if (!f)
                break;
            qhat++;
            ccc++;
#ifdef DEBUG
            printf("corrected qhat = %08lx\n", (long)qhat);
#endif
            b=0;
            for (i=0;i<yd;i++) {
                m = hd(x,i+j)-hd(y,i)-b;
                if (m < 0) {
                    b = 1;
                    m = HLMAX + 1 + m;
                }
                else
                    b = 0;
                if ((i+j)%2) 
                    x->p[(i+j)/2] = LOW(x->p[(i+j)/2]) | (m << HALFDIGITBITS);
                else
                    x->p[(i+j)/2] = (HIGH(x->p[(i+j)/2]) << HALFDIGITBITS) | m;
            }
            if (b) {
                if ((j+i)%2) 
                    x->p[(i+j)/2] -= b << HALFDIGITBITS;
                else
                    x->p[(i+j)/2] -= b;
            }
#ifdef DEBUG
            printf("x after cor=");
            for (i=2*(x->sz)-1;i>=0;i--)
                printf(" %04lx", (long)hd(x,i));
            printf("\n");
#endif
        }
        if (j%2) 
            q->p[j/2] |= qhat << HALFDIGITBITS;
        else
            q->p[j/2] |= qhat;
#ifdef DEBUG
            printf("x after cor=");
            for (i=xd - 1;i>=0;i--)
                printf(" %04lx", (long)hd(q,i));
            printf("\n");
#endif
    }
    _mpz_realloc(r,(mp_size)(yy->sz));
    zero(r);
    urshift(r,x,ns);
    mpz_set(rr,r);
    mpz_set(qq,q);
    mpz_clear(x); mpz_clear(y);
    mpz_clear(q);  mpz_clear(r);
    free(q); free(x); free(y); free(r);
}

void mpz_add(zz,x,y)
MP_INT *zz;MP_INT *x;MP_INT *y;
{
    int mg;
    MP_INT *z;
    if (x->sn == 0) {
        mpz_set(zz,y);
        return;
    }
    if (y->sn == 0) {
        mpz_set(zz,x);
        return;
    }
    z = (MP_INT *)malloc(sizeof(MP_INT));
    mpz_init(z);

    if (x->sn > 0 && y->sn > 0) {
        uadd(z,x,y); z->sn = 1;
    }
    else if (x->sn < 0 && y->sn < 0) {
        uadd(z,x,y); z->sn = -1;
    }
    else {
        /* signs differ */
        if ((mg = ucmp(x,y)) == 0) {
            zero(z);
        }
        else if (mg > 0) {  /* abs(y) < abs(x) */
            usub(z,x,y);    
            z->sn = (x->sn > 0 && y->sn < 0) ? 1 : (-1);
        }
        else { /* abs(y) > abs(x) */
            usub(z,y,x);    
            z->sn = (x->sn < 0 && y->sn > 0) ? 1 : (-1);
        }
    }
    mpz_set(zz,z);
    mpz_clear(z);
    free(z);
}

void mpz_add_ui(x,y,n)
MP_INT *x;MP_INT *y; unsigned long n;
{
    MP_INT z;
    mpz_init_set_ui(&z,n);
    mpz_add(x,y,&z);
    mpz_clear(&z);
}

int mpz_cmp_ui(x,n)
MP_INT *x; unsigned long n;
{
    MP_INT z; int ret;
    mpz_init_set_ui(&z,n);
    ret=mpz_cmp(x,&z);
    mpz_clear(&z);
    return ret;
}

int mpz_cmp_si(x,n)
MP_INT *x; long n;
{
    MP_INT z; int ret;
    mpz_init(&z);
    mpz_set_si(&z,n);
    ret=mpz_cmp(x,&z);
    mpz_clear(&z);
    return ret;
}


void mpz_mul_ui(x,y,n)
MP_INT *x;MP_INT *y; unsigned long n;
{
    MP_INT z;
    mpz_init_set_ui(&z,n);
    mpz_mul(x,y,&z);
    mpz_clear(&z);
}
    
    
/* z = x - y  -- just use mpz_add - I'm lazy */
void mpz_sub(z,x,y)
MP_INT *z;MP_INT *x; MP_INT *y;
{
    MP_INT u;
    mpz_init(&u);
    mpz_set(&u,y);
    u.sn = -(u.sn);
    mpz_add(z,x,&u);
    mpz_clear(&u);
}

void mpz_sub_ui(x,y,n)
MP_INT *x;MP_INT *y; unsigned long n;
{
    MP_INT z;
    mpz_init_set_ui(&z,n);
    mpz_sub(x,y,&z);
    mpz_clear(&z);
}

void mpz_div(q,x,y)
MP_INT *q; MP_INT *x; MP_INT *y;
{
    MP_INT r; 
    short sn1 = x->sn, sn2 = y->sn;
    mpz_init(&r);
    udiv(q,&r,x,y);
    q->sn = sn1*sn2;
    if (uzero(q))
        q->sn = 0;
    mpz_clear(&r);
}

void mpz_mdiv(q,x,y)
MP_INT *q; MP_INT *x; MP_INT *y;
{
    MP_INT r; 
    short sn1 = x->sn, sn2 = y->sn, qsign;
    mpz_init(&r);
    udiv(q,&r,x,y);
    qsign = q->sn = sn1*sn2;
    if (uzero(q))
        q->sn = 0;
    /* now if r != 0 and q < 0 we need to round q towards -inf */
    if (!uzero(&r) && qsign < 0) 
        mpz_sub_ui(q,q,1L);
    mpz_clear(&r);
}

void mpz_mod(r,x,y)
MP_INT *r; MP_INT *x; MP_INT *y;
{
    MP_INT q;
    short sn = x->sn;
    mpz_init(&q);
    if (x->sn == 0) {
        zero(r);
        return;
    }
    udiv(&q,r,x,y);
    r->sn = sn;
    if (uzero(r))
        r->sn = 0;
    mpz_clear(&q);
}

void mpz_divmod(q,r,x,y)
MP_INT *q; MP_INT *r; MP_INT *x; MP_INT *y;
{
    short sn1 = x->sn, sn2 = y->sn;
    if (x->sn == 0) {
        zero(q);
        zero(r);
        return;
    }
    udiv(q,r,x,y);
    q->sn = sn1*sn2;
    if (uzero(q))
        q->sn = 0;
    r->sn = sn1;
    if (uzero(r))
        r->sn = 0;
}

void mpz_mmod(r,x,y)
MP_INT *r; MP_INT *x; MP_INT *y;
{
    MP_INT q;
    short sn1 = x->sn, sn2 = y->sn;
    mpz_init(&q);
    if (sn1 == 0) {
        zero(r);
        return;
    }
    udiv(&q,r,x,y);
    if (uzero(r)) {
        r->sn = 0;
        return;
    }
    q.sn = sn1*sn2;
    if (q.sn > 0) 
        r->sn = sn1;
    else if (sn1 < 0 && sn2 > 0) {
        r->sn = 1;
        mpz_sub(r,y,r);
    }
    else {
        r->sn = 1;
        mpz_add(r,y,r);
    }
}

void mpz_mdivmod(q,r,x,y)
MP_INT *q;MP_INT *r; MP_INT *x; MP_INT *y;
{
    short sn1 = x->sn, sn2 = y->sn, qsign;
    if (sn1 == 0) {
        zero(q);
        zero(r);
        return;
    }
    udiv(q,r,x,y);
    qsign = q->sn = sn1*sn2;
    if (uzero(r)) {
        /* q != 0, since q=r=0 would mean x=0, which was tested above */
        r->sn = 0;
        return;
    }
    if (q->sn > 0) 
        r->sn = sn1;
    else if (sn1 < 0 && sn2 > 0) {
        r->sn = 1;
        mpz_sub(r,y,r);
    }
    else {
        r->sn = 1;
        mpz_add(r,y,r);
    }
    if (uzero(q))
        q->sn = 0;
    /* now if r != 0 and q < 0 we need to round q towards -inf */
    if (!uzero(r) && qsign < 0) 
        mpz_sub_ui(q,q,1L);
}

void mpz_mod_ui(x,y,n)
MP_INT *x;MP_INT *y; unsigned long n;
{
    MP_INT z;
    mpz_init(&z);
    mpz_set_ui(&z,n);
    mpz_mod(x,y,&z);
    mpz_clear(&z);
}

void mpz_mmod_ui(x,y,n)
MP_INT *x;MP_INT *y; unsigned long n;
{
    MP_INT z;
    mpz_init(&z);
    mpz_set_ui(&z,n);
    mpz_mmod(x,y,&z);
    mpz_clear(&z);
}

void mpz_div_ui(x,y,n)
MP_INT *x;MP_INT *y; unsigned long n;
{
    MP_INT z;
    mpz_init_set_ui(&z,n);
    mpz_div(x,y,&z);
    mpz_clear(&z);
}

void mpz_mdiv_ui(x,y,n)
MP_INT *x;MP_INT *y; unsigned long n;
{
    MP_INT z;
    mpz_init_set_ui(&z,n);
    mpz_mdiv(x,y,&z);
    mpz_clear(&z);
}
void mpz_divmod_ui(q,x,y,n)
MP_INT *q;MP_INT *x;MP_INT *y; unsigned long n;
{
    MP_INT z;
    mpz_init_set_ui(&z,n);
    mpz_divmod(q,x,y,&z);
    mpz_clear(&z);
}

void mpz_mdivmod_ui(q,x,y,n)
MP_INT *q;MP_INT *x;MP_INT *y; unsigned long n;
{
    MP_INT z;
    mpz_init_set_ui(&z,n);
    mpz_mdivmod(q,x,y,&z);
    mpz_clear(&z);
}


/* 2<=base <=36 - this overestimates the optimal value, which is OK */
unsigned int mpz_sizeinbase(x,base)
MP_INT *x; int base;
{
    int i,j;
    int bits = digits(x) * DIGITBITS;
    if (base < 2 || base > 36)
        fatal("mpz_sizeinbase: invalid base");
    for (j=0,i=1; i<=base;i*=2,j++)
        ;
    return ((bits)/(j-1)+1);
}

char *mpz_get_str(s,base,x)
char *s;  int base; MP_INT *x; {
    MP_INT xx,q,r,bb;
    char *p,*t,*ps;
    int sz = mpz_sizeinbase(x,base);
    mp_limb d;
    if (base < 2 || base > 36)
        return s;
    t = (char *)malloc(sz+2);
    if (!t) 
        fatal("cannot allocate memory in mpz_get_str");
    if (!s) {
        s=(char *)malloc(sz+2);
        if (!s) 
            fatal("cannot allocate memory in mpz_get_str");
    }
    if (uzero(x)) {
        *s='0';
        *(s+1)='\0';
        return s;
    }
    mpz_init(&xx); mpz_init(&q); mpz_init(&r); mpz_init(&bb);
    mpz_set(&xx,x);
    mpz_set_ui(&bb,(unsigned long)base);
    ps = s;
    if (x->sn < 0) {
        *ps++= '-';
        xx.sn = 1;
    }
    p = t;
    while (!uzero(&xx)) {
        udiv(&xx,&r,&xx,&bb);
        d = r.p[0];
        if (d < 10) 
            *p++  = (char) (r.p[0] + '0');
        else 
            *p++  = (char) (r.p[0] + -10 + 'a');
    }

    p--;
    for (;p>=t;p--,ps++)
        *ps = *p;
    *ps='\0';
    
    mpz_clear(&q); mpz_clear(&r); free(t);  
    return s;
}


int mpz_set_str(x,s,base)
MP_INT *x; char *s; int base;
{
    int l,i;
    int retval = 0;
    MP_INT t,m,bb;
    short sn;
    unsigned int k;
    mpz_init(&m);
    mpz_init(&t);
    mpz_init(&bb);
    mpz_set_ui(&m,1L);
    zero(x);
    while (*s==' ' || *s=='\t' || *s=='\n')
        s++;
    if (*s == '-') {
        sn = -1; s++;
    }
    else
        sn = 1;
    if (base == 0) {
        if (*s == '0') {
            if (*(s+1) == 'x' || *(s+1) == 'X') {
                base = 16;
                s+=2;   /* toss 0x */
            }
            else {
                base = 8;
                s++;    /* toss 0 */
            }
        }
        else
            base=10;
    }
    if (base < 2 || base > 36)
        fatal("mpz_set_str: invalid base");
    mpz_set_ui(&bb,(unsigned long)base);
    l=strlen(s);
    for (i = l-1; i>=0; i--) {
        if (s[i]==' ' || s[i]=='\t' || s[i]=='\n') 
            continue;
        if (s[i] >= '0' && s[i] <= '9')
            k = (unsigned int)s[i] - (unsigned int)'0';
        else if (s[i] >= 'A' && s[i] <= 'Z')
            k = (unsigned int)s[i] - (unsigned int)'A'+10;
        else if (s[i] >= 'a' && s[i] <= 'z')
            k = (unsigned int)s[i] - (unsigned int)'a'+10;
        else {
            retval = (-1);
            break;
        }
        if (k >= base) {
            retval = (-1);
            break;
        }
        mpz_mul_ui(&t,&m,(unsigned long)k);
        mpz_add(x,x,&t);
        mpz_mul(&m,&m,&bb);
#ifdef DEBUG
        printf("k=%d\n",k);
        printf("t=%s\n",mpz_get_str(NULL,10,&t));
        printf("x=%s\n",mpz_get_str(NULL,10,x));
        printf("m=%s\n",mpz_get_str(NULL,10,&m));
#endif
    }
    if (x->sn)
        x->sn = sn;
    mpz_clear(&m);
    mpz_clear(&bb);
    mpz_clear(&t);
    return retval;
}

int mpz_init_set_str(x,s,base)
MP_INT *x; char *s; int base;
{
    mpz_init(x); return mpz_set_str(x,s,base);
}

#define mpz_randombyte (rand()& 0xff)

void mpz_random(x,size)
MP_INT *x; unsigned int size;
{
    unsigned int bits = size * LONGBITS;
    unsigned int digits = bits/DIGITBITS;
    unsigned int oflow = bits % DIGITBITS;
    unsigned int i,j;
    long t;
    if (oflow)
        digits++;
    _mpz_realloc(x,(mp_size)digits);

    for (i=0; i<digits; i++) {
        t = 0;
        for (j=0; j < DIGITBITS; j+=8)
            t = (t << 8) | mpz_randombyte; 
        (x->p)[i] = t & LMAX;
    }
    if (oflow) 
        (x->p)[digits-1] &= (((mp_limb)1 << oflow) - 1);
}
void mpz_random2(x,size)
MP_INT *x; unsigned int size;
{
    unsigned int bits = size * LONGBITS;
    unsigned int digits = bits/DIGITBITS;
    unsigned int oflow = bits % DIGITBITS;
    unsigned int i,j;
    long t;
    if (oflow)
        digits++;
    _mpz_realloc(x,(mp_size)digits);

    for (i=0; i<digits; i++) {
        t = 0;
        for (j=0; j < DIGITBITS; j+=8)
            t = (t << 8) | mpz_randombyte; 
        (x->p)[i] = (t & LMAX) % 2;
    }
    if (oflow) 
        (x->p)[digits-1] &= (((mp_limb)1 << oflow) - 1);
}

size_t mpz_size(x)
MP_INT *x;
{
    int bits = (x->sz)*DIGITBITS;
    size_t r;
    
    r = bits/LONGBITS;
    if (bits % LONGBITS)
        r++;
    return r;
}

void mpz_abs(x,y)
MP_INT *x; MP_INT *y;
{
    if (x!=y)
        mpz_set(x,y);
    if (uzero(x))
        x->sn = 0;
    else
        x->sn = 1;
}

void mpz_neg(x,y)
MP_INT *x; MP_INT *y;
{
    if (x!=y)
        mpz_set(x,y);
    x->sn = -(y->sn);
}

void mpz_fac_ui(x,n)
MP_INT *x; unsigned long n;
{
    mpz_set_ui(x,1L);
    if (n==0 || n == 1)
        return;
    for (;n>1;n--)
        mpz_mul_ui(x,x,n);
}

void mpz_gcd(gg,aa,bb)
MP_INT *gg;
MP_INT *aa;
MP_INT *bb;
{
    MP_INT *g,*t,*a,*b;
    int freeg;
    
    t = (MP_INT *)malloc(sizeof(MP_INT));
    g = (MP_INT *)malloc(sizeof(MP_INT));
    a = (MP_INT *)malloc(sizeof(MP_INT));
    b = (MP_INT *)malloc(sizeof(MP_INT));
    if (!a || !b || !g || !t)
        fatal("mpz_gcd: cannot allocate memory");
    mpz_init(g); mpz_init(t); mpz_init(a); mpz_init(b);
    mpz_abs(a,aa); mpz_abs(b,bb); 

    while (!uzero(b)) {
        mpz_mod(t,a,b);
        mpz_set(a,b);
        mpz_set(b,t);
    }
    
    mpz_set(gg,a);
    mpz_clear(g); mpz_clear(t); mpz_clear(a); mpz_clear(b);
    free(g); free(t); free(a); free(b);
}

void mpz_gcdext(gg,xx,yy,aa,bb)
MP_INT *gg,*xx,*yy,*aa,*bb;
{
    MP_INT *x, *y, *g, *u, *q;

    if (uzero(bb)) {
        mpz_set(gg,aa);
        mpz_set_ui(xx,1L);
        if (yy)
            mpz_set_ui(yy,0L);
        return;
    }
    
    g = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(g);
    q = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(q);
    y = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(y);
    x = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(x);
    u = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(u);

    if (!g || !q || !y || !x || !u)
        fatal("mpz_gcdext: cannot allocate memory");

    mpz_divmod(q,u,aa,bb);
    mpz_gcdext(g,x,y,bb,u);
    if (yy) {
        mpz_mul(u,q,y);
        mpz_sub(yy,x,u);
    }
    mpz_set(xx,y);
    mpz_set(gg,g);
    mpz_clear(g); mpz_clear(q); mpz_clear(y); mpz_clear(x); mpz_clear(u);
    free(g); free(q); free(y); free(x); free(u);
}


/*
 *    a is a quadratic residue mod b if
 *    x^2 = a mod b      has an integer solution x.
 *
 *    J(a,b) = if a==1 then 1 else
 *             if a is even then J(a/2,b) * ((-1)^(b^2-1)/8))
 *             else J(b mod a, a) * ((-1)^((a-1)*(b-1)/4)))
 *
 *  b>0  b odd
 *
 *
 */

int mpz_jacobi(ac,bc)
MP_INT *ac, *bc;
{
    int sgn = 1;
    unsigned long c;
    MP_INT *t,*a,*b; 
    if (bc->sn <=0)
        fatal("mpz_jacobi call with b <= 0");
    if (even(bc))
        fatal("mpz_jacobi call with b even");

    init(t); init(a); init(b);

    /* if ac < 0, then we use the fact that 
     *  (-1/b)= -1 if b mod 4 == 3
     *          +1 if b mod 4 == 1
     */

    if (mpz_cmp_ui(ac,0L) < 0 && (lowdigit(bc) % 4) == 3)
        sgn=-sgn;

    mpz_abs(a,ac); mpz_set(b,bc);

    /* while (a > 1) { */
    while(mpz_cmp_ui(a,1L) > 0) {

        if (even(a)) {

            /* c = b % 8 */
            c = lowdigit(b) & 0x07;

            /* b odd, then (b^2-1)/8 is even iff (b%8 == 3,5) */

            /* if b % 8 == 3 or 5 */
            if (c == 3 || c == 5)
                sgn = -sgn;

            /* a = a / 2 */
            mpz_div_2exp(a,a,1L); 

        } 
        else {
            /* note: (a-1)(b-1)/4 odd iff a mod 4==b mod 4==3 */

            /* if a mod 4 == 3 and b mod 4 == 3 */
            if (((lowdigit(a) & 3) == 3) && ((lowdigit(b) & 3) == 3))
                sgn = -sgn;

            /* set (a,b) = (b mod a,a) */
            mpz_set(t,a); mpz_mmod(a,b,t); mpz_set(b,t);
        } 
    }

    /* if a == 0 then sgn = 0 */
    if (uzero(a))
        sgn=0;

    clear(t); clear(a); clear(b);
    return (sgn);
}

void mpz_and(z,x,y) /* not the most efficient way to do this */
MP_INT *z,*x,*y;
{
    int i,sz;
    sz = imax(x->sz, y->sz);
    _mpz_realloc(z,(mp_size)sz);
    for (i=0; i < sz; i++)
        (z->p)[i] = dg(x,i) & dg(y,i);
    if (x->sn < 0 && y->sn < 0)
        z->sn = (-1);
    else
        z->sn = 1;
    if (uzero(z))
        z->sn = 0;
}

void mpz_or(z,x,y)  /* not the most efficient way to do this */
MP_INT *z,*x,*y;
{
    int i,sz;
    sz = imax(x->sz, y->sz);
    _mpz_realloc(z,(mp_size)sz);
    for (i=0; i < sz; i++)
        (z->p)[i] = dg(x,i) | dg(y,i);
    if (x->sn < 0 || y->sn < 0)
        z->sn = (-1);
    else
        z->sn = 1;
    if (uzero(z))
        z->sn = 0;
}

void mpz_xor(z,x,y)  /* not the most efficient way to do this */
MP_INT *z,*x,*y;
{
    int i,sz;
    sz = imax(x->sz, y->sz);
    _mpz_realloc(z,(mp_size)sz);
    for (i=0; i < sz; i++)
        (z->p)[i] = dg(x,i) ^ dg(y,i);
    if ((x->sn <= 0 && y->sn > 0) || (x->sn > 0 && y->sn <=0))
        z->sn = (-1);
    else
        z->sn = 1;
    if (uzero(z))
        z->sn = 0;
}
void mpz_pow_ui(zz,x,e)
MP_INT *zz, *x;
unsigned long e;
{
    MP_INT *t;
    unsigned long mask = (1L<< (LONGBITS-1));
    
    if (e==0) {
        mpz_set_ui(zz,1L);
        return;
    }

    init(t);
    mpz_set(t,x);
    for (;!(mask &e); mask>>=1) 
        ;
    mask>>=1;
    for (;mask!=0; mask>>=1) {
        mpz_mul(t,t,t);
        if (e & mask)
            mpz_mul(t,t,x);
    }
    mpz_set(zz,t);
    clear(t);
}

void mpz_powm(zz,x,ee,n)
MP_INT *zz,*x,*ee,*n;
{
    MP_INT *t,*e;
    struct is *stack = NULL;
    int k,i;
    
    if (uzero(ee)) {
        mpz_set_ui(zz,1L);
        return;
    }

    if (ee->sn < 0) {
        return;
    }
    init(e); init(t); mpz_set(e,ee);

    for (k=0;!uzero(e);k++,mpz_div_2exp(e,e,1L))
        push(lowdigit(e) & 1,&stack);
    k--;
    i=pop(&stack);

    mpz_mod(t,x,n);  /* t=x%n */

    for (i=k-1;i>=0;i--) {
        mpz_mul(t,t,t); 
        mpz_mod(t,t,n);  
        if (pop(&stack)) {
            mpz_mul(t,t,x); 
            mpz_mod(t,t,n);
        }
    }
    mpz_set(zz,t);
    clear(t);
}

void mpz_powm_ui(z,x,e,n)
MP_INT *z,*x,*n; unsigned long e;
{
    MP_INT f;
    mpz_init(&f);
    mpz_set_ui(&f,e);
    mpz_powm(z,x,&f,n);
    mpz_clear(&f);
}
    
    

/* Miller-Rabin */
static int witness(x,n)
MP_INT *x, *n;
{
    MP_INT *t,*e, *e1;
    struct is *stack = NULL;
    int trivflag = 0;
    int k,i;
    
    init(e1); init(e); init(t); mpz_sub_ui(e,n,1L); mpz_set(e1,e);

    for (k=0;!uzero(e);k++,mpz_div_2exp(e,e,1L))
        push(lowdigit(e) & 1,&stack);
    k--;
    i=pop(&stack);

    mpz_mod(t,x,n);  /* t=x%n */

    for (i=k-1;i>=0;i--) {
        trivflag = !mpz_cmp_ui(t,1L) || !mpz_cmp(t,e1);
        mpz_mul(t,t,t); mpz_mod(t,t,n);  
        if (!trivflag && !mpz_cmp_ui(t,1L)) {
            clear(t); clear(e); clear(e1);
            return 1;
        }
            
        if (pop(&stack)) {
            mpz_mul(t,t,x); 
            mpz_mod(t,t,n);
        }
    }
    if (mpz_cmp_ui(t,1L))  { /* t!=1 */
        clear(t); clear(e); clear(e1);
        return 1;
    }
    else {
        clear(t); clear(e); clear(e1);
        return 0;
    }
}

unsigned long smallp[] = {2,3,5,7,11,13,17};
int mpz_probab_prime_p(nn,s)
MP_INT *nn; int s;
{   
    MP_INT *a,*n;
    int j,k,i;
    long t;

    if (uzero(nn))
        return 0;
    init(n); mpz_abs(n,nn);
    if (!mpz_cmp_ui(n,1L)) {
        clear(n);
        return 0;
    }
    init(a);
    for (i=0;i<sizeof(smallp)/sizeof(unsigned long); i++) {
        mpz_mod_ui(a,n,smallp[i]);
        if (uzero(a)) {
            clear(a); clear(n); return 0;
        }
    }
    _mpz_realloc(a,(mp_size)(n->sz));
    for (k=0; k<s; k++) {

        /* generate a "random" a */
            for (i=0; i<n->sz; i++) {
                t = 0;
                for (j=0; j < DIGITBITS; j+=8)
                    t = (t << 8) | mpz_randombyte; 
                (a->p)[i] = t & LMAX;
            }
            a->sn = 1;
            mpz_mod(a,a,n);

        if (witness(a,n)) {
            clear(a); clear(n);
            return 0;
        }
    }
    clear(a);clear(n);
    return 1;
}   


/* Square roots are found by Newton's method, but since we are working
 * with integers we have to consider carefully the termination conditions.
 * If we are seeking x = floor(sqrt(a)), the iteration is
 *     x_{n+1} = floor ((floor (a/x_n) + x_n)/2) == floor ((a + x_n^2)/(2*x))
 * If eps_n represents the error (exactly, eps_n and sqrt(a) real) in the form:
 *     x_n = (1 + eps_n) sqrt(a)
 * then it is easy to show that for a >= 4
 *     if 0 <= eps_n, then either 0 <= eps_{n+1} <= (eps_n^2)/2
 *                         or x_{n+1} = floor (sqrt(a))
 *     if x_n = floor (sqrt(a)), then either x_{n+1} = floor (sqrt(a))
 *                               or x_{n+1} = floor (sqrt(a)) + 1
 * Therefore, if the initial approximation x_0 is chosen so that
 *     0 <= eps_0 <= 1/2
 * then the error remains postive and monotonically decreasing with
 *     eps_n <= 2 ^ (-(2^(n+1) - 1))
 * unless we reach the correct floor(sqrt(a)), after which we may oscillate
 * between it and the value one higher.
 * We choose the number of iterations, k, to guarantee
 *     eps_k sqrt(a) < 1,  so x_k <= floor (sqrt (a)) + 1
 * so finally x_k is either the required answer or one too big (which we
 * check by a simple multiplication and comparison).
 *
 * The calculation of the initial approximation *assumes* DIGITBITS is even.
 * It probably always will be, so for now let's leave the code simple,
 * clear and all reachable in testing!
 */
void mpz_sqrtrem (root, rem, a)
    MP_INT *root;
    MP_INT *rem;
    MP_INT *a;
{
    MP_INT tmp;
    MP_INT r;
    mp_limb z;
    unsigned long j, h;
    int k;

    if (a->sn < 0)
        /* Negative operand, nothing correct we can do */
        return;

    /* Now, a good enough approximation to the root is obtained by writing
     *     a = z * 2^(2j) + y,  4 <= z <= 15 and 0 <= y < 2^(2j)
     * then taking
     *     root = ciel (sqrt(z+1)) * 2^j
     */

    for (j = a->sz - 1; j != 0 && (a->p)[j] == 0; j--);
    z = (a->p)[j];
    if (z < 4) {
        if (j == 0) {
            /* Special case for small operand (main interation not valid) */
            mpz_set_ui (root, (z>0) ? 1L : 0L);
            if (rem)
                mpz_set_ui (rem, (z>1) ? z-1L : 0L);
            return;
        } else {
            z = (z << 2) + ((a->p)[j-1] >> (DIGITBITS-2));
            j = j * DIGITBITS - 2;
        }
    } else {
        j = j * DIGITBITS;
        while (z & ~0xf) {
            z >>= 2;
            j += 2;
        }
    }
    j >>= 1;                            /* z and j now as described above */
    for (k=1, h=4; h < j+3; k++, h*=2); /* 2^(k+1) >= j+3, since a < 2^(2j+4) */
    mpz_init_set_ui (&r, (z>8) ? 4L : 3L);
    mpz_mul_2exp (&r, &r, j);

#ifdef DEBUG
    printf ("mpz_sqrtrem: z = %lu, j = %lu, k = %d\n", (long)z, j, k);
#endif

    /* Main iteration */
    mpz_init (&tmp);
    for ( ; k > 0; k--) {
        mpz_div (&tmp, a, &r);
        mpz_add (&r, &tmp, &r);
        mpz_div_2exp (&r, &r, 1L);
    }

    /* Adjust result (which may be one too big) and set remainder */
    mpz_mul (&tmp, &r, &r);
    if (mpz_cmp (&tmp, a) > 0) {
        mpz_sub_ui (root, &r, 1L);
        if (rem) {
            mpz_sub (rem, a, &tmp);
            mpz_add (rem, rem, &r);
            mpz_add (rem, rem, &r);
            mpz_sub_ui (rem, rem, 1L);
        }
    } else {
        mpz_set (root, &r);
        if (rem)
            mpz_sub (rem, a, &tmp);
    }

    mpz_clear (&tmp);
    mpz_clear (&r);
}


void mpz_sqrt (root, a)
    MP_INT *root;
    MP_INT *a;
{
    mpz_sqrtrem (root, NULL, a);
}