#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of shell archive."
# Contents:  notes gmp.c gmp.h
# Wrapped by mch@squid on Fri Feb  4 16:02:06 1994
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'notes' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'notes'\"
else
echo shar: Extracting \"'notes'\" \(3670 characters\)
sed "s/^X//" >'notes' <<'END_OF_FILE'
XWELCOME TO FGMP.
X
XFGMP is a public domain implementation of a subset of the GNU gmp library
Xwith the same API. 
X
XFor instance, you can link the following trivial program with either 
Xthis code, or libgmp.a and get the same results.
X------------
X#include <stdio.h>
X#include "gmp.h"
Xmain()
X{
X    MP_INT a; MP_INT b; MP_INT c;
X
X    mpz_init_set_ui(&a,1L); mpz_init_set_ui(&b,2L); mpz_init(&c);
X    mpz_add(&c,&a,&b);
X    printf("\n%s\n", mpz_get_str(NULL,10,&c));
X}
X
X------------
X
XFGMP is really in the public domain. You can do whatever you want with
Xit.
X
XI wrote FGMP so that we would all have access to a (truly free)
Ximplementation of this subset of the API of GNU libgmp. I encourage
Xeveryone to distribute this as widely as possible.
X
XIf you need more documentation, I suggest you look at the file
Xgmp.texi which is included with the GNU gmp library. 
X
XYou can send me bug reports, implementations of missing functions, flames
Xand rants by Email. 
X
XAny submissions of new code to be integrated into fgmp must also be 
Xplaced in the public domain (For the particularly dense, you can 
Xrelease a new fgmp yourself under different licensing terms. This 
Xis a condition for including a submission in a release of FGMP that 
XI personally prepare).
X
XMark Henderson <markh@wimsey.bc.ca>
X
X---
XThis is FGMP 1.0.1
X
XI hearby place this file and all of FGMP in the public domain.
X
XThanks to Paul Rouse <par@r-cube.demon.co.uk> for changes to get fgmp 
Xto work on a 286 MSDOS compiler, the functions mpz_sqrt and 
Xmpz_sqrtrem, plus other general bug fixes. 
X
XThanks also to Erick Gallesio <eg@kaolin.unice.fr> for a fix
Xto mpz_init_set_str
X
XDefine B64 if your "long" type is 64 bits. Otherwise we assume 32
Xbit longs. (The 64 bit version hasn't been tested enough)
X
X
X
XPlatforms:
X
XFGMP has been ported to the following platforms.:
X
XLinux 0.99 (gcc 2.3, i486)
XIBM RS6000/AIX 3.2 (IBM cc, gcc 2.4)
XHP/UX 8.0 (HP cc, gcc 2.3)
XSun OS 4.1, Sun 3/4 (sun cc, gcc 2.2)
XSGI Irix 4.0.5 (SGI cc, gcc 2.3)
XDEC Alpha OSF/1 (DEC cc) 
X    (only lightly tested, 64 bit longs do make a difference,
X     thanks to DEC for providing access on axposf.pa.dec.com). 
X     Define B64 for this platform.
XMSDOS 286 C compiler (see credits above)
X
X---
XSome differences between gmp and fgmp
X
X1. fgmp is considerably slower than gmp
X2. fgmp does not implement the following:
X    all mpq_*
X    internal mpn_* functions
X    mpz_perfect_square_p
X    mpz_inp_raw, mpz_out_raw
X    mp_set_memory_functions, mpz_out_str, mpz_inp_str
X3. fgmp implements the following in addition to the routines in GNU gmp.
X    int mpz_jacobi(MP_INT *a, MP_INT *b)
X    - finds the jacobi symbol (a/b)
X4. mpz_sizeinbase often overestimates the exact value
X
X5. To convert your gmp based program to fgmp (subject to the
Xabove)
X
X- recompile your source. Make sure to include the gmp.h file included
X  with fgmp rather than that included with gmp. (The point is to recompile
X  all files which include gmp.h)
X- link with gmp.o instead of libgmp.a
X
XHere's a complete sorted list of function implemented in fgmp:
X
X_mpz_realloc
Xmpz_abs
Xmpz_add
Xmpz_add_ui
Xmpz_and
Xmpz_clear
Xmpz_cmp
Xmpz_cmp_si
Xmpz_cmp_ui
Xmpz_div
Xmpz_div_2exp
Xmpz_div_ui
Xmpz_divmod
Xmpz_divmod_ui
Xmpz_fac_ui
Xmpz_gcd
Xmpz_gcdext
Xmpz_get_si
Xmpz_get_str
Xmpz_get_ui
Xmpz_init
Xmpz_init_set
Xmpz_init_set_si
Xmpz_init_set_str
Xmpz_init_set_ui
Xmpz_jacobi
Xmpz_mdiv
Xmpz_mdiv_ui
Xmpz_mdivmod
Xmpz_mdivmod_ui
Xmpz_mmod
Xmpz_mmod_ui
Xmpz_mod
Xmpz_mod_2exp
Xmpz_mod_ui
Xmpz_mul
Xmpz_mul_2exp
Xmpz_mul_ui
Xmpz_neg
Xmpz_or
Xmpz_pow_ui
Xmpz_powm
Xmpz_powm_ui
Xmpz_probab_prime_p
Xmpz_random
Xmpz_random2
Xmpz_set
Xmpz_set_si
Xmpz_set_str
Xmpz_set_ui
Xmpz_size
Xmpz_sizeinbase
Xmpz_sqrt
Xmpz_sqrtrem
Xmpz_sub
Xmpz_sub_ui
Xmpz_xor
END_OF_FILE
if test 3670 -ne `wc -c <'notes'`; then
    echo shar: \"'notes'\" unpacked with wrong size!
fi
# end of 'notes'
fi
if test -f 'gmp.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'gmp.c'\"
else
echo shar: Extracting \"'gmp.c'\" \(39198 characters\)
sed "s/^X//" >'gmp.c' <<'END_OF_FILE'
X/*
X * FREE GMP - a public domain implementation of a subset of the 
X *           gmp library
X *
X * I hearby place the file in the public domain.
X *
X * Do whatever you want with this code. Change it. Sell it. Claim you
X *  wrote it. 
X * Bugs, complaints, flames, rants: please send email to 
X *    Mark Henderson <markh@wimsey.bc.ca>
X * I'm already aware that fgmp is considerably slower than gmp
X *
X * CREDITS:
X *  Paul Rouse <par@r-cube.demon.co.uk> - generic bug fixes, mpz_sqrt and
X *    mpz_sqrtrem, and modifications to get fgmp to compile on a system 
X *    with int and long of different sizes (specifically MSDOS,286 compiler).
X *  Also see the file "notes" included with the fgmp distribution, for
X *    more credits.
X *
X * VERSION 1.0.1
X */
X
X#include "gmp.h"
X#ifndef NULL
X#define NULL ((void *)0)
X#endif
X
X#define iabs(x) ((x>0) ? (x) : (-x))
X#define imax(x,y) ((x>y)?x:y)
X#define LONGBITS (sizeof(mp_limb) * 8)
X#define DIGITBITS (LONGBITS - 2)
X#define HALFDIGITBITS ((LONGBITS -2)/2)
X#ifndef init
X#define init(g) { g = (MP_INT *)malloc(sizeof(MP_INT));  mpz_init(g); }
X#endif
X#ifndef clear
X#define clear(g) { mpz_clear(g); free(g); }
X#endif
X#ifndef even
X#define even(a) (!((a)->p[0] & 1))
X#endif
X#ifndef odd
X#define odd(a) (((a)->p[0] & 1))
X#endif
X
X
X#ifndef B64
X/* 
X * The values below are for 32 bit machines (i.e. machines with a 
X *  32 bit long type)
X * You'll need to change them, if you're using something else
X * If DIGITBITS is odd, see the comment at the top of mpz_sqrtrem
X */
X#define LMAX 0x3fffffffL
X#define LC   0xc0000000L
X#define OVMASK 0x2
X#define CMASK (LMAX+1)
X#define FC ((double)CMASK)
X#define HLMAX 0x7fffL
X#define HCMASK (HLMAX + 1)
X#define HIGH(x) (((x) & 0x3fff8000L) >> 15)
X#define LOW(x)  ((x) & 0x7fffL)
X#else
X
X/* 64 bit long type */
X#define LMAX 0x3fffffffffffffffL
X#define LC 0xc000000000000000L
X#define OVMASK 0x2
X#define CMASK (LMAX+1)
X#define FC ((double)CMASK)
X#define HLMAX 0x7fffffffL
X#define HCMASK (HLMAX + 1)
X#define HIGH(x) (((x) & 0x3fffffff80000000L) >> 31)
X#define LOW(x) ((x) & 0x7fffffffL)
X
X#endif
X
X#define hd(x,i)  (((i)>=2*(x->sz))? 0:(((i)%2) ? HIGH((x)->p[(i)/2]) \
X    : LOW((x)->p[(i)/2])))
X#define dg(x,i) ((i < x->sz) ? (x->p)[i] : 0)
X
X#ifdef __GNUC__
X#define INLINE inline
X#else
X#define INLINE
X#endif
X
X#define lowdigit(x) (((x)->p)[0])
X
Xstruct is {
X    mp_limb v;
X    struct is *next;
X};
X
X
XINLINE static void push(i,sp)
Xmp_limb i;struct is **sp;
X{
X    struct is *tmp;
X    tmp = *sp;
X    *sp = (struct is *) malloc(sizeof(struct is));
X    (*sp)->v = i;
X    (*sp)->next=tmp;
X}
X
XINLINE static mp_limb pop(sp)
Xstruct is **sp;
X{
X    struct is *tmp;
X    mp_limb i;
X    if (!(*sp))
X        return (-1);
X    tmp = *sp;
X    *sp = (*sp)->next;
X    i = tmp->v;
X    tmp->v = 0;
X    free(tmp);
X    return i;
X}
X
Xvoid fatal(s)
Xchar *s;
X{
X    fprintf(stderr,"%s\n",s);
X    exit(123);
X}
X
Xvoid mpz_init(s)
XMP_INT *s;
X{
X    s->p = (mp_limb *)malloc(sizeof(mp_limb)*2);
X#ifdef DEBUG
X    printf("malloc: 8 bytes, %08lx\n", (long)(s->p));
X#endif
X    if (!(s->p))
X        fatal("mpz_init: cannot allocate memory");
X    (s->p)[0] = 0;
X    (s->p)[1] = 0;
X    s->sn=0;
X    s->sz=2;
X}
X
Xvoid mpz_init_set(s,t)
XMP_INT *s,*t;
X{
X    int i;
X    s->p = (mp_limb *)malloc(sizeof(mp_limb) * t->sz);
X    if (!(s->p))
X        fatal("mpz_init: cannot allocate memory");
X    for (i=0;i < t->sz ; i++)
X        (s->p)[i] = (t->p)[i];
X
X    s->sn = t->sn;
X    s->sz = t->sz;
X}
X
Xvoid mpz_init_set_ui(s,v)
XMP_INT *s;
Xunsigned long v;
X{
X    s->p = (mp_limb *)malloc(sizeof(mp_limb)*2);
X    if (!(s->p))
X        fatal("mpz_init: cannot allocate memory");
X    s->p[0] = (v & LMAX);
X    s->p[1] = (v & LC) >> DIGITBITS;
X    if (v) 
X        s->sn = 1;
X    else
X        s->sn = 0;
X    s->sz = 2;
X}
X
Xvoid mpz_init_set_si(y,v)
XMP_INT *y;
Xlong v;
X{
X    y->p = (mp_limb *)malloc(sizeof(mp_limb)*2);
X    if (!(y->p))
X        fatal("mpz_init: cannot allocate memory");
X    if (v < 0) {
X        y->sn = -1;
X        y->p[0] = (-v) & LMAX;
X        y->p[1] = ((-v) & LC) >> DIGITBITS;
X    }
X    else if (v > 0) {
X        y->sn = 1;
X        y->p[0] = v & LMAX;
X        y->p[1] = (v & LC) >> DIGITBITS;
X    }
X    else {
X        y->sn=0;
X        y->p[0] = 0;
X        y->p[1] = 0;
X    }
X    y -> sz = 2;
X}
X
X
X
Xvoid mpz_clear(s)
XMP_INT *s;
X{
X#ifdef DEBUG
X    printf("free: %08lx\n", (long)(s->p));
X#endif
X    if (s->p)
X        free(s->p);
X    s->p=NULL;
X    s->sn=0;
X    s->sz=0;
X}
X
X/* number of leading zero bits in digit */
Xstatic int lzb(a)
Xmp_limb a;
X{
X    mp_limb i; int j;
X    j = 0;
X    for (i = ((mp_limb)1 << (DIGITBITS-1)); i && !(a&i) ; j++,i>>=1) 
X        ;
X    return j;
X}
X
Xvoid _mpz_realloc(x,size)
XMP_INT *x; mp_size size;
X{
X    int i;
X    if (size > 1 && x->sz < size) {
X#ifdef DEBUG
X    printf("realloc %08lx to size = %ld ", (long)(x->p),(long)(size));
X#endif
X
X        if (x->p) 
X            x->p=(mp_limb *)realloc(x->p,size * sizeof(mp_limb));
X        else
X            x->p=(mp_limb *)malloc(size * sizeof(mp_limb));
X#ifdef DEBUG
X    printf("becomes %08lx\n", (long)(x->p));
X#endif
X
X        if (!(x->p))
X            fatal("_mpz_realloc: cannot allocate memory");
X        for (i=x->sz;i<size;i++)
X            (x->p)[i] = 0;
X        x->sz = size;
X    } else if (size < x->sz)
X        for (i=size; i<x->sz; i++)
X            (x->p)[i] = 0;
X}
X
Xvoid dgset(x,i,n)
XMP_INT *x; unsigned int i; mp_limb n;
X{
X    if (n) {
X        if (i >= x->sz) {
X            _mpz_realloc(x,(mp_size)(i+1));
X            x->sz = i+1;
X        }
X        (x->p)[i] = n;
X    }
X}
X
Xstatic unsigned int digits(x)
XMP_INT *x;
X{
X    int i;
X    for (i = (x->sz) - 1; i>=0 && (x->p)[i] == 0 ; i--) 
X        ;
X    return i+1;
X}
X        
X/* y = x */
Xvoid mpz_set(y,x)
XMP_INT *y; MP_INT *x; 
X{
X    unsigned int i,k = x->sz;
X    if (y->sz < k) {
X        k=digits(x);
X        _mpz_realloc(y,(mp_size)k);
X    }
X    if (y->sz > x->sz) {
X        mpz_clear(y); mpz_init(y); _mpz_realloc(y,(mp_size)(x->sz));
X    }
X
X    for (i=0;i < k; i++)
X        (y->p)[i] = (x->p)[i];
X
X    for (;i<y->sz;i++)
X        (y->p)[i] = 0;
X
X    y->sn = x->sn;
X}
X
Xvoid mpz_set_ui(y,v)
XMP_INT *y; unsigned long v; {
X    int i;
X    for (i=1;i<y->sz;i++)
X        y->p[i] = 0;
X    y->p[0] = (v & LMAX);
X    y->p[1] = (v & LC) >> (DIGITBITS);
X    if (v)
X        y->sn=1;
X    else
X        y->sn=0;
X}
X
Xunsigned long mpz_get_ui(y) 
XMP_INT *y; {
X    return (y->p[0] | (y->p[1] << DIGITBITS));
X}
X
Xlong mpz_get_si(y)
XMP_INT *y; {
X    if (y->sn == 0)
X        return 0;
X    else
X        return (y->sn * (y->p[0] | (y->p[1] & 1) << DIGITBITS));
X}
X        
Xvoid mpz_set_si(y,v)
XMP_INT *y; long v; {
X    int i;
X    for (i=1;i<y->sz;i++)
X        y->p[i] = 0;
X    if (v < 0) {
X        y->sn = -1;
X        y->p[0] = (-v) & LMAX;
X        y->p[1] = ((-v) & LC) >> DIGITBITS;
X    }
X    else if (v > 0) {
X        y->sn = 1;
X        y->p[0] = v & LMAX;
X        y->p[1] = (v & LC) >> DIGITBITS;
X    }
X    else {
X        y->sn=0;
X        y->p[0] = 0;
X        y->p[1] = 0;
X    }
X}
X
X/* z = x + y, without regard for sign */
Xstatic void uadd(z,x,y)
XMP_INT *z; MP_INT *x; MP_INT *y;
X{
X    mp_limb c;
X    int i;
X    MP_INT *t;
X
X    if (y->sz < x->sz) {
X        t=x; x=y; y=t;
X    }
X
X    /* now y->sz >= x->sz */
X
X    _mpz_realloc(z,(mp_size)((y->sz)+1));
X
X    c=0;
X    for (i=0;i<x->sz;i++) {
X        if (( z->p[i] = y->p[i] + x->p[i] + c ) & CMASK) {
X            c=1; 
X            (z->p[i]) &=LMAX;
X        }
X        else 
X            c=0;
X    }
X    for (;i<y->sz;i++) {
X        if ((z->p[i] = (y->p[i] + c)) & CMASK)
X            z->p[i]=0;
X        else
X            c=0;
X    }
X    (z->p)[y->sz]=c;
X}
X
X/* z = y - x, ignoring sign */
X/* precondition: abs(y) >= abs(x) */
Xstatic void usub(z,y,x)
XMP_INT *z; MP_INT *y; MP_INT *x; 
X{
X    int i;
X    mp_limb b,m;
X    _mpz_realloc(z,(mp_size)(y->sz));
X    b=0;
X    for (i=0;i<y->sz;i++) {
X        m=((y->p)[i]-b)-dg(x,i);
X        if (m < 0) {
X            b = 1;
X            m = LMAX + 1 + m;
X        }
X        else
X            b = 0;
X        z->p[i] = m;
X    }
X}
X
X/* compare abs(x) and abs(y) */
Xstatic int ucmp(y,x)
XMP_INT *y;MP_INT *x;
X{
X    int i;
X    for (i=imax(x->sz,y->sz)-1;i>=0;i--) {
X        if (dg(y,i) < dg(x,i))
X            return (-1);
X        else if (dg(y,i) > dg(x,i))
X            return 1;
X    }
X    return 0;
X}
X
Xint mpz_cmp(x,y)
XMP_INT *x; MP_INT *y;
X{
X    int abscmp;
X    if (x->sn < 0 && y->sn > 0)
X        return (-1);
X    if (x->sn > 0 && y->sn < 0)
X        return 1;
X    abscmp=ucmp(x,y);
X    if (x->sn >=0 && y->sn >=0)
X        return abscmp;
X    if (x->sn <=0 && y->sn <=0)
X        return (-abscmp);
X}
X    
X/* is abs(x) == 0 */
Xstatic int uzero(x)
XMP_INT *x;
X{
X    int i;
X    for (i=0; i < x->sz; i++)
X        if ((x->p)[i] != 0)
X            return 0;
X    return 1;
X}
X
Xstatic void zero(x)
XMP_INT *x;
X{
X    int i;
X    x->sn=0;
X    for (i=0;i<x->sz;i++)
X        (x->p)[i] = 0;
X}
X
X
X/* w = u * v */
Xvoid mpz_mul(ww,u,v)
XMP_INT *ww;MP_INT *u; MP_INT *v; {
X    int i,j;
X    mp_limb t0,t1,t2,t3;
X    mp_limb cc;
X    MP_INT *w;
X
X    w = (MP_INT *)malloc(sizeof(MP_INT));
X    mpz_init(w);
X    _mpz_realloc(w,(mp_size)(u->sz + v->sz));
X    for (j=0; j < 2*u->sz; j++) {
X        cc = (mp_limb)0;
X        t3 = hd(u,j);
X            for (i=0; i < 2*v->sz; i++) {
X                t0 = t3 * hd(v,i);
X                t1 = HIGH(t0); t0 = LOW(t0);
X                if ((i+j)%2) 
X                    t2 = HIGH(w->p[(i+j)/2]);
X                else
X                    t2 = LOW(w->p[(i+j)/2]);
X                t2 += cc;
X                if (t2 & HCMASK) {
X                    cc = 1; t2&=HLMAX;
X                }
X                else
X                    cc = 0;
X                t2 += t0;
X                if (t2 & HCMASK) {
X                    cc++ ; t2&=HLMAX;
X                }
X                cc+=t1;
X                if ((i+j)%2)
X                    w->p[(i+j)/2] = LOW(w->p[(i+j)/2]) |
X                        (t2 << HALFDIGITBITS);
X                else
X                    w->p[(i+j)/2] = (HIGH(w->p[(i+j)/2]) << HALFDIGITBITS) |
X                        t2;
X                    
X            }
X        if (cc) {
X            if ((j+i)%2) 
X                w->p[(i+j)/2] += cc << HALFDIGITBITS;
X            else
X                w->p[(i+j)/2] += cc;
X        }
X    }
X    w->sn = (u->sn) * (v->sn);
X    if (w != ww) {
X        mpz_set(ww,w);
X        mpz_clear(w);
X        free(w);
X    }
X}
X
X/* n must be < DIGITBITS */
Xstatic void urshift(c1,a,n)
XMP_INT *c1;MP_INT *a;int n;
X{
X    mp_limb cc = 0;
X    if (n >= DIGITBITS)
X        fatal("urshift: n >= DIGITBITS");
X    if (n == 0)
X        mpz_set(c1,a);
X    else {
X        MP_INT c; int i;
X        mp_limb rm = (((mp_limb)1<<n) - 1);
X        mpz_init(&c); _mpz_realloc(&c,(mp_size)(a->sz));
X        for (i=a->sz-1;i >= 0; i--) {
X            c.p[i] = ((a->p[i] >> n) | cc) & LMAX;
X            cc = (a->p[i] & rm) << (DIGITBITS - n);
X        }
X        mpz_set(c1,&c);
X        mpz_clear(&c);
X    }
X}
X    
X/* n must be < DIGITBITS */
Xstatic void ulshift(c1,a,n)
XMP_INT *c1;MP_INT *a;int n;
X{
X    mp_limb cc = 0;
X    if (n >= DIGITBITS)
X        fatal("ulshift: n >= DIGITBITS");
X    if (n == 0)
X        mpz_set(c1,a);
X    else {
X        MP_INT c; int i;
X        mp_limb rm = (((mp_limb)1<<n) - 1) << (DIGITBITS -n);
X        mpz_init(&c); _mpz_realloc(&c,(mp_size)(a->sz + 1));
X        for (i=0;i < a->sz; i++) {
X            c.p[i] = ((a->p[i] << n) | cc) & LMAX;
X            cc = (a->p[i] & rm) >> (DIGITBITS -n);
X        }
X        c.p[i] = cc;
X        mpz_set(c1,&c);
X        mpz_clear(&c);
X    }
X}
X
Xvoid mpz_div_2exp(z,x,e) 
XMP_INT *z; MP_INT *x; unsigned long e;
X{
X    short sn = x->sn;
X    if (e==0)
X        mpz_set(z,x);
X    else {
X        unsigned long i;
X        long digs = (e / DIGITBITS);
X        unsigned long bs = (e % (DIGITBITS));
X        MP_INT y; mpz_init(&y);
X        _mpz_realloc(&y,(mp_size)((x->sz) - digs));
X        for (i=0; i < (x->sz - digs); i++)
X            (y.p)[i] = (x->p)[i+digs];
X        if (bs) {
X            urshift(z,&y,bs);
X        }
X        else {
X            mpz_set(z,&y);
X        }
X        if (uzero(z))
X            z->sn = 0;
X        else
X            z->sn = sn;
X        mpz_clear(&y);
X    }
X}
X
Xvoid mpz_mul_2exp(z,x,e) 
XMP_INT *z; MP_INT *x; unsigned long e;
X{
X    short sn = x->sn;
X    if (e==0)
X        mpz_set(z,x);
X    else {
X        unsigned long i;
X        long digs = (e / DIGITBITS);
X        unsigned long bs = (e % (DIGITBITS));
X        MP_INT y; mpz_init(&y);
X        _mpz_realloc(&y,(mp_size)((x->sz)+digs));
X        for (i=digs;i<((x->sz) + digs);i++)
X            (y.p)[i] = (x->p)[i - digs];
X        if (bs) {
X            ulshift(z,&y,bs);
X        }
X        else {
X            mpz_set(z,&y);
X        }
X        z->sn = sn;
X        mpz_clear(&y);
X    }
X}
X
Xvoid mpz_mod_2exp(z,x,e) 
XMP_INT *z; MP_INT *x; unsigned long e;
X{
X    short sn = x->sn;
X    if (e==0)
X        mpz_set_ui(z,0L);
X    else {
X        unsigned long i;
X        MP_INT y;
X        long digs = (e / DIGITBITS);
X        unsigned long bs = (e % (DIGITBITS));
X        if (digs > x->sz || (digs == (x->sz) && bs)) {
X            mpz_set(z,x);
X            return;
X        }
X            
X        mpz_init(&y);
X        if (bs)
X            _mpz_realloc(&y,(mp_size)(digs+1));
X        else
X            _mpz_realloc(&y,(mp_size)digs);
X        for (i=0; i<digs; i++)
X            (y.p)[i] = (x->p)[i];
X        if (bs) {
X            (y.p)[digs] = ((x->p)[digs]) & (((mp_limb)1 << bs) - 1);
X        }
X        mpz_set(z,&y);
X        if (uzero(z))
X            z->sn = 0;
X        else
X            z->sn = sn;
X        mpz_clear(&y);
X    }
X}
X    
X
X/* internal routine to compute x/y and x%y ignoring signs */
Xstatic void udiv(qq,rr,xx,yy)
XMP_INT *qq; MP_INT *rr; MP_INT *xx; MP_INT *yy;
X{
X    MP_INT *q, *x, *y, *r;
X    int ns,xd,yd,j,f,i,ccc;
X    mp_limb zz,z,qhat,b,u,m;
X    ccc=0;
X
X    if (uzero(yy))
X        return;
X    q = (MP_INT *)malloc(sizeof(MP_INT));
X    r = (MP_INT *)malloc(sizeof(MP_INT));
X    x = (MP_INT *)malloc(sizeof(MP_INT)); 
X    y = (MP_INT *)malloc(sizeof(MP_INT));
X    if (!x || !y || !q || !r)
X        fatal("udiv: cannot allocate memory");
X    mpz_init(q); mpz_init(x);mpz_init(y);mpz_init(r);
X    _mpz_realloc(x,(mp_size)((xx->sz)+1));
X    yd = digits(yy);
X    ns = lzb(yy->p[yd-1]);
X    ulshift(x,xx,ns);
X    ulshift(y,yy,ns);
X    xd = digits(x); 
X    _mpz_realloc(q,(mp_size)xd);
X    xd*=2; yd*=2;
X    z = hd(y,yd-1);;
X    for (j=(xd-yd);j>=0;j--) {
X#ifdef DEBUG
X    printf("y=");
X    for (i=yd-1;i>=0;i--)
X        printf(" %04lx", (long)hd(y,i));
X    printf("\n");
X    printf("x=");
X    for (i=xd-1;i>=0;i--)
X        printf(" %04lx", (long)hd(x,i));
X    printf("\n");
X    printf("z=%08lx\n",(long)z);
X#endif
X        
X        if (z == LMAX)
X            qhat = hd(x,j+yd);
X        else {
X            qhat = ((hd(x,j+yd)<< HALFDIGITBITS) + hd(x,j+yd-1)) / (z+1);
X        }
X#ifdef DEBUG
X    printf("qhat=%08lx\n",(long)qhat);
X    printf("hd=%04lx/%04lx\n",(long)hd(x,j+yd),(long)hd(x,j+yd-1));
X#endif
X        b = 0; zz=0;
X        if (qhat) {
X            for (i=0; i < yd; i++) {
X                zz = qhat * hd(y,i);
X                u = hd(x,i+j);
X                u-=b;
X                if (u<0) {
X                    b=1; u+=HLMAX+1;
X                }
X                else
X                    b=0;
X                u-=LOW(zz);
X                if (u < 0) {
X                    b++;
X                    u+=HLMAX+1;
X                }
X                b+=HIGH(zz);
X                if ((i+j)%2) 
X                    x->p[(i+j)/2] = LOW(x->p[(i+j)/2]) | (u << HALFDIGITBITS);
X                else
X                    x->p[(i+j)/2] = (HIGH(x->p[(i+j)/2]) 
X                        << HALFDIGITBITS) | u;
X            }
X            if (b) {
X                if ((j+i)%2) 
X                    x->p[(i+j)/2] -= b << HALFDIGITBITS;
X                else
X                    x->p[(i+j)/2] -= b;
X            }
X        }
X#ifdef DEBUG
X        printf("x after sub=");
X        for (i=xd-1;i>=0;i--)
X            printf(" %04lx", (long)hd(x,i));
X        printf("\n");
X#endif
X        for(;;zz++) {
X            f=1;
X            if (!hd(x,j+yd)) {
X                for(i=yd-1; i>=0; i--) {
X                    if (hd(x,j+i) > hd(y,i)) {
X                        f=1;
X                        break;
X                    }
X                    if (hd(x,j+i) < hd(y,i)) {
X                        f=0;
X                        break;
X                    }
X                }
X            }
X            if (!f)
X                break;
X            qhat++;
X            ccc++;
X#ifdef DEBUG
X            printf("corrected qhat = %08lx\n", (long)qhat);
X#endif
X            b=0;
X            for (i=0;i<yd;i++) {
X                m = hd(x,i+j)-hd(y,i)-b;
X                if (m < 0) {
X                    b = 1;
X                    m = HLMAX + 1 + m;
X                }
X                else
X                    b = 0;
X                if ((i+j)%2) 
X                    x->p[(i+j)/2] = LOW(x->p[(i+j)/2]) | (m << HALFDIGITBITS);
X                else
X                    x->p[(i+j)/2] 
X                        = (HIGH(x->p[(i+j)/2]) << HALFDIGITBITS) | m;
X            }
X            if (b) {
X                if ((j+i)%2) 
X                    x->p[(i+j)/2] -= b << HALFDIGITBITS;
X                else
X                    x->p[(i+j)/2] -= b;
X            }
X#ifdef DEBUG
X            printf("x after cor=");
X            for (i=2*(x->sz)-1;i>=0;i--)
X                printf(" %04lx", (long)hd(x,i));
X            printf("\n");
X#endif
X        }
X        if (j%2) 
X            q->p[j/2] |= qhat << HALFDIGITBITS;
X        else
X            q->p[j/2] |= qhat;
X#ifdef DEBUG
X            printf("x after cor=");
X            for (i=xd - 1;i>=0;i--)
X                printf(" %04lx", (long)hd(q,i));
X            printf("\n");
X#endif
X    }
X    _mpz_realloc(r,(mp_size)(yy->sz));
X    zero(r);
X    urshift(r,x,ns);
X    mpz_set(rr,r);
X    mpz_set(qq,q);
X    mpz_clear(x); mpz_clear(y);
X    mpz_clear(q);  mpz_clear(r);
X    free(q); free(x); free(y); free(r);
X}
X
Xvoid mpz_add(zz,x,y)
XMP_INT *zz;MP_INT *x;MP_INT *y;
X{
X    int mg;
X    MP_INT *z;
X    if (x->sn == 0) {
X        mpz_set(zz,y);
X        return;
X    }
X    if (y->sn == 0) {
X        mpz_set(zz,x);
X        return;
X    }
X    z = (MP_INT *)malloc(sizeof(MP_INT));
X    mpz_init(z);
X
X    if (x->sn > 0 && y->sn > 0) {
X        uadd(z,x,y); z->sn = 1;
X    }
X    else if (x->sn < 0 && y->sn < 0) {
X        uadd(z,x,y); z->sn = -1;
X    }
X    else {
X        /* signs differ */
X        if ((mg = ucmp(x,y)) == 0) {
X            zero(z);
X        }
X        else if (mg > 0) {  /* abs(y) < abs(x) */
X            usub(z,x,y);    
X            z->sn = (x->sn > 0 && y->sn < 0) ? 1 : (-1);
X        }
X        else { /* abs(y) > abs(x) */
X            usub(z,y,x);    
X            z->sn = (x->sn < 0 && y->sn > 0) ? 1 : (-1);
X        }
X    }
X    mpz_set(zz,z);
X    mpz_clear(z);
X    free(z);
X}
X
Xvoid mpz_add_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_add(x,y,&z);
X    mpz_clear(&z);
X}
X
Xint mpz_cmp_ui(x,n)
XMP_INT *x; unsigned long n;
X{
X    MP_INT z; int ret;
X    mpz_init_set_ui(&z,n);
X    ret=mpz_cmp(x,&z);
X    mpz_clear(&z);
X    return ret;
X}
X
Xint mpz_cmp_si(x,n)
XMP_INT *x; long n;
X{
X    MP_INT z; int ret;
X    mpz_init(&z);
X    mpz_set_si(&z,n);
X    ret=mpz_cmp(x,&z);
X    mpz_clear(&z);
X    return ret;
X}
X
X
Xvoid mpz_mul_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_mul(x,y,&z);
X    mpz_clear(&z);
X}
X    
X    
X/* z = x - y  -- just use mpz_add - I'm lazy */
Xvoid mpz_sub(z,x,y)
XMP_INT *z;MP_INT *x; MP_INT *y;
X{
X    MP_INT u;
X    mpz_init(&u);
X    mpz_set(&u,y);
X    u.sn = -(u.sn);
X    mpz_add(z,x,&u);
X    mpz_clear(&u);
X}
X
Xvoid mpz_sub_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_sub(x,y,&z);
X    mpz_clear(&z);
X}
X
Xvoid mpz_div(q,x,y)
XMP_INT *q; MP_INT *x; MP_INT *y;
X{
X    MP_INT r; 
X    short sn1 = x->sn, sn2 = y->sn;
X    mpz_init(&r);
X    udiv(q,&r,x,y);
X    q->sn = sn1*sn2;
X    if (uzero(q))
X        q->sn = 0;
X    mpz_clear(&r);
X}
X
Xvoid mpz_mdiv(q,x,y)
XMP_INT *q; MP_INT *x; MP_INT *y;
X{
X    MP_INT r; 
X    short sn1 = x->sn, sn2 = y->sn, qsign;
X    mpz_init(&r);
X    udiv(q,&r,x,y);
X    qsign = q->sn = sn1*sn2;
X    if (uzero(q))
X        q->sn = 0;
X    /* now if r != 0 and q < 0 we need to round q towards -inf */
X    if (!uzero(&r) && qsign < 0) 
X        mpz_sub_ui(q,q,1L);
X    mpz_clear(&r);
X}
X
Xvoid mpz_mod(r,x,y)
XMP_INT *r; MP_INT *x; MP_INT *y;
X{
X    MP_INT q;
X    short sn = x->sn;
X    mpz_init(&q);
X    if (x->sn == 0) {
X        zero(r);
X        return;
X    }
X    udiv(&q,r,x,y);
X    r->sn = sn;
X    if (uzero(r))
X        r->sn = 0;
X    mpz_clear(&q);
X}
X
Xvoid mpz_divmod(q,r,x,y)
XMP_INT *q; MP_INT *r; MP_INT *x; MP_INT *y;
X{
X    short sn1 = x->sn, sn2 = y->sn;
X    if (x->sn == 0) {
X        zero(q);
X        zero(r);
X        return;
X    }
X    udiv(q,r,x,y);
X    q->sn = sn1*sn2;
X    if (uzero(q))
X        q->sn = 0;
X    r->sn = sn1;
X    if (uzero(r))
X        r->sn = 0;
X}
X
Xvoid mpz_mmod(r,x,y)
XMP_INT *r; MP_INT *x; MP_INT *y;
X{
X    MP_INT q;
X    short sn1 = x->sn, sn2 = y->sn;
X    mpz_init(&q);
X    if (sn1 == 0) {
X        zero(r);
X        return;
X    }
X    udiv(&q,r,x,y);
X    if (uzero(r)) {
X        r->sn = 0;
X        return;
X    }
X    q.sn = sn1*sn2;
X    if (q.sn > 0) 
X        r->sn = sn1;
X    else if (sn1 < 0 && sn2 > 0) {
X        r->sn = 1;
X        mpz_sub(r,y,r);
X    }
X    else {
X        r->sn = 1;
X        mpz_add(r,y,r);
X    }
X}
X
Xvoid mpz_mdivmod(q,r,x,y)
XMP_INT *q;MP_INT *r; MP_INT *x; MP_INT *y;
X{
X    short sn1 = x->sn, sn2 = y->sn, qsign;
X    if (sn1 == 0) {
X        zero(q);
X        zero(r);
X        return;
X    }
X    udiv(q,r,x,y);
X    qsign = q->sn = sn1*sn2;
X    if (uzero(r)) {
X        /* q != 0, since q=r=0 would mean x=0, which was tested above */
X        r->sn = 0;
X        return;
X    }
X    if (q->sn > 0) 
X        r->sn = sn1;
X    else if (sn1 < 0 && sn2 > 0) {
X        r->sn = 1;
X        mpz_sub(r,y,r);
X    }
X    else {
X        r->sn = 1;
X        mpz_add(r,y,r);
X    }
X    if (uzero(q))
X        q->sn = 0;
X    /* now if r != 0 and q < 0 we need to round q towards -inf */
X    if (!uzero(r) && qsign < 0) 
X        mpz_sub_ui(q,q,1L);
X}
X
Xvoid mpz_mod_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init(&z);
X    mpz_set_ui(&z,n);
X    mpz_mod(x,y,&z);
X    mpz_clear(&z);
X}
X
Xvoid mpz_mmod_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init(&z);
X    mpz_set_ui(&z,n);
X    mpz_mmod(x,y,&z);
X    mpz_clear(&z);
X}
X
Xvoid mpz_div_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_div(x,y,&z);
X    mpz_clear(&z);
X}
X
Xvoid mpz_mdiv_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_mdiv(x,y,&z);
X    mpz_clear(&z);
X}
Xvoid mpz_divmod_ui(q,x,y,n)
XMP_INT *q;MP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_divmod(q,x,y,&z);
X    mpz_clear(&z);
X}
X
Xvoid mpz_mdivmod_ui(q,x,y,n)
XMP_INT *q;MP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_mdivmod(q,x,y,&z);
X    mpz_clear(&z);
X}
X
X
X/* 2<=base <=36 - this overestimates the optimal value, which is OK */
Xunsigned int mpz_sizeinbase(x,base)
XMP_INT *x; int base;
X{
X    int i,j;
X    int bits = digits(x) * DIGITBITS;
X    if (base < 2 || base > 36)
X        fatal("mpz_sizeinbase: invalid base");
X    for (j=0,i=1; i<=base;i*=2,j++)
X        ;
X    return ((bits)/(j-1)+1);
X}
X
Xchar *mpz_get_str(s,base,x)
Xchar *s;  int base; MP_INT *x; {
X    MP_INT xx,q,r,bb;
X    char *p,*t,*ps;
X    int sz = mpz_sizeinbase(x,base);
X    mp_limb d;
X    if (base < 2 || base > 36)
X        return s;
X    t = (char *)malloc(sz+2);
X    if (!t) 
X        fatal("cannot allocate memory in mpz_get_str");
X    if (!s) {
X        s=(char *)malloc(sz+2);
X        if (!s) 
X            fatal("cannot allocate memory in mpz_get_str");
X    }
X    if (uzero(x)) {
X        *s='0';
X        *(s+1)='\0';
X        return s;
X    }
X    mpz_init(&xx); mpz_init(&q); mpz_init(&r); mpz_init(&bb);
X    mpz_set(&xx,x);
X    mpz_set_ui(&bb,(unsigned long)base);
X    ps = s;
X    if (x->sn < 0) {
X        *ps++= '-';
X        xx.sn = 1;
X    }
X    p = t;
X    while (!uzero(&xx)) {
X        udiv(&xx,&r,&xx,&bb);
X        d = r.p[0];
X        if (d < 10) 
X            *p++  = (char) (r.p[0] + '0');
X        else 
X            *p++  = (char) (r.p[0] + -10 + 'a');
X    }
X
X    p--;
X    for (;p>=t;p--,ps++)
X        *ps = *p;
X    *ps='\0';
X    
X    mpz_clear(&q); mpz_clear(&r); free(t);  
X    return s;
X}
X
X
Xint mpz_set_str(x,s,base)
XMP_INT *x; char *s; int base;
X{
X    int l,i;
X    int retval = 0;
X    MP_INT t,m,bb;
X    short sn;
X    unsigned int k;
X    mpz_init(&m);
X    mpz_init(&t);
X    mpz_init(&bb);
X    mpz_set_ui(&m,1L);
X    zero(x);
X    while (*s==' ' || *s=='\t' || *s=='\n')
X        s++;
X    if (*s == '-') {
X        sn = -1; s++;
X    }
X    else
X        sn = 1;
X    if (base == 0) {
X        if (*s == '0') {
X            if (*(s+1) == 'x' || *(s+1) == 'X') {
X                base = 16;
X                s+=2;   /* toss 0x */
X            }
X            else {
X                base = 8;
X                s++;    /* toss 0 */
X            }
X        }
X        else
X            base=10;
X    }
X    if (base < 2 || base > 36)
X        fatal("mpz_set_str: invalid base");
X    mpz_set_ui(&bb,(unsigned long)base);
X    l=strlen(s);
X    for (i = l-1; i>=0; i--) {
X        if (s[i]==' ' || s[i]=='\t' || s[i]=='\n') 
X            continue;
X        if (s[i] >= '0' && s[i] <= '9')
X            k = (unsigned int)s[i] - (unsigned int)'0';
X        else if (s[i] >= 'A' && s[i] <= 'Z')
X            k = (unsigned int)s[i] - (unsigned int)'A'+10;
X        else if (s[i] >= 'a' && s[i] <= 'z')
X            k = (unsigned int)s[i] - (unsigned int)'a'+10;
X        else {
X            retval = (-1);
X            break;
X        }
X        if (k >= base) {
X            retval = (-1);
X            break;
X        }
X        mpz_mul_ui(&t,&m,(unsigned long)k);
X        mpz_add(x,x,&t);
X        mpz_mul(&m,&m,&bb);
X#ifdef DEBUG
X        printf("k=%d\n",k);
X        printf("t=%s\n",mpz_get_str(NULL,10,&t));
X        printf("x=%s\n",mpz_get_str(NULL,10,x));
X        printf("m=%s\n",mpz_get_str(NULL,10,&m));
X#endif
X    }
X    if (x->sn)
X        x->sn = sn;
X    mpz_clear(&m);
X    mpz_clear(&bb);
X    mpz_clear(&t);
X    return retval;
X}
X
Xint mpz_init_set_str(x,s,base)
XMP_INT *x; char *s; int base;
X{
X    mpz_init(x); return mpz_set_str(x,s,base);
X}
X
X#define mpz_randombyte (rand()& 0xff)
X
Xvoid mpz_random(x,size)
XMP_INT *x; unsigned int size;
X{
X    unsigned int bits = size * LONGBITS;
X    unsigned int digits = bits/DIGITBITS;
X    unsigned int oflow = bits % DIGITBITS;
X    unsigned int i,j;
X    long t;
X    if (oflow)
X        digits++;
X    _mpz_realloc(x,(mp_size)digits);
X
X    for (i=0; i<digits; i++) {
X        t = 0;
X        for (j=0; j < DIGITBITS; j+=8)
X            t = (t << 8) | mpz_randombyte; 
X        (x->p)[i] = t & LMAX;
X    }
X    if (oflow) 
X        (x->p)[digits-1] &= (((mp_limb)1 << oflow) - 1);
X    x->sn = 0;
X    for (i=0; i<digits; i++)
X        if ((x->p)[i]) {
X            x->sn = 1;
X            break;
X        }
X}
Xvoid mpz_random2(x,size)
XMP_INT *x; unsigned int size;
X{
X    unsigned int bits = size * LONGBITS;
X    unsigned int digits = bits/DIGITBITS;
X    unsigned int oflow = bits % DIGITBITS;
X    unsigned int i,j;
X    long t;
X    if (oflow)
X        digits++;
X    _mpz_realloc(x,(mp_size)digits);
X
X    for (i=0; i<digits; i++) {
X        t = 0;
X        for (j=0; j < DIGITBITS; j+=8)
X            t = (t << 8) | mpz_randombyte; 
X        (x->p)[i] = (t & LMAX) % 2;
X    }
X    if (oflow) 
X        (x->p)[digits-1] &= (((mp_limb)1 << oflow) - 1);
X    x->sn = 0;
X    for (i=0; i<digits; i++)
X        if ((x->p)[i]) {
X            x->sn = 1;
X            break;
X        }
X}
X
Xsize_t mpz_size(x)
XMP_INT *x;
X{
X    int bits = (x->sz)*DIGITBITS;
X    size_t r;
X    
X    r = bits/LONGBITS;
X    if (bits % LONGBITS)
X        r++;
X    return r;
X}
X
Xvoid mpz_abs(x,y)
XMP_INT *x; MP_INT *y;
X{
X    if (x!=y)
X        mpz_set(x,y);
X    if (uzero(x))
X        x->sn = 0;
X    else
X        x->sn = 1;
X}
X
Xvoid mpz_neg(x,y)
XMP_INT *x; MP_INT *y;
X{
X    if (x!=y)
X        mpz_set(x,y);
X    x->sn = -(y->sn);
X}
X
Xvoid mpz_fac_ui(x,n)
XMP_INT *x; unsigned long n;
X{
X    mpz_set_ui(x,1L);
X    if (n==0 || n == 1)
X        return;
X    for (;n>1;n--)
X        mpz_mul_ui(x,x,n);
X}
X
Xvoid mpz_gcd(gg,aa,bb)
XMP_INT *gg;
XMP_INT *aa;
XMP_INT *bb;
X{
X    MP_INT *g,*t,*a,*b;
X    int freeg;
X    
X    t = (MP_INT *)malloc(sizeof(MP_INT));
X    g = (MP_INT *)malloc(sizeof(MP_INT));
X    a = (MP_INT *)malloc(sizeof(MP_INT));
X    b = (MP_INT *)malloc(sizeof(MP_INT));
X    if (!a || !b || !g || !t)
X        fatal("mpz_gcd: cannot allocate memory");
X    mpz_init(g); mpz_init(t); mpz_init(a); mpz_init(b);
X    mpz_abs(a,aa); mpz_abs(b,bb); 
X
X    while (!uzero(b)) {
X        mpz_mod(t,a,b);
X        mpz_set(a,b);
X        mpz_set(b,t);
X    }
X    
X    mpz_set(gg,a);
X    mpz_clear(g); mpz_clear(t); mpz_clear(a); mpz_clear(b);
X    free(g); free(t); free(a); free(b);
X}
X
Xvoid mpz_gcdext(gg,xx,yy,aa,bb)
XMP_INT *gg,*xx,*yy,*aa,*bb;
X{
X    MP_INT *x, *y, *g, *u, *q;
X
X    if (uzero(bb)) {
X        mpz_set(gg,aa);
X        mpz_set_ui(xx,1L);
X        if (yy)
X            mpz_set_ui(yy,0L);
X        return;
X    }
X    
X    g = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(g);
X    q = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(q);
X    y = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(y);
X    x = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(x);
X    u = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(u);
X
X    if (!g || !q || !y || !x || !u)
X        fatal("mpz_gcdext: cannot allocate memory");
X
X    mpz_divmod(q,u,aa,bb);
X    mpz_gcdext(g,x,y,bb,u);
X    if (yy) {
X        mpz_mul(u,q,y);
X        mpz_sub(yy,x,u);
X    }
X    mpz_set(xx,y);
X    mpz_set(gg,g);
X    mpz_clear(g); mpz_clear(q); mpz_clear(y); mpz_clear(x); mpz_clear(u);
X    free(g); free(q); free(y); free(x); free(u);
X}
X
X
X/*
X *    a is a quadratic residue mod b if
X *    x^2 = a mod b      has an integer solution x.
X *
X *    J(a,b) = if a==1 then 1 else
X *             if a is even then J(a/2,b) * ((-1)^(b^2-1)/8))
X *             else J(b mod a, a) * ((-1)^((a-1)*(b-1)/4)))
X *
X *  b>0  b odd
X *
X *
X */
X
Xint mpz_jacobi(ac,bc)
XMP_INT *ac, *bc;
X{
X    int sgn = 1;
X    unsigned long c;
X    MP_INT *t,*a,*b; 
X    if (bc->sn <=0)
X        fatal("mpz_jacobi call with b <= 0");
X    if (even(bc))
X        fatal("mpz_jacobi call with b even");
X
X    init(t); init(a); init(b);
X
X    /* if ac < 0, then we use the fact that 
X     *  (-1/b)= -1 if b mod 4 == 3
X     *          +1 if b mod 4 == 1
X     */
X
X    if (mpz_cmp_ui(ac,0L) < 0 && (lowdigit(bc) % 4) == 3)
X        sgn=-sgn;
X
X    mpz_abs(a,ac); mpz_set(b,bc);
X
X    /* while (a > 1) { */
X    while(mpz_cmp_ui(a,1L) > 0) {
X
X        if (even(a)) {
X
X            /* c = b % 8 */
X            c = lowdigit(b) & 0x07;
X
X            /* b odd, then (b^2-1)/8 is even iff (b%8 == 3,5) */
X
X            /* if b % 8 == 3 or 5 */
X            if (c == 3 || c == 5)
X                sgn = -sgn;
X
X            /* a = a / 2 */
X            mpz_div_2exp(a,a,1L); 
X
X        } 
X        else {
X            /* note: (a-1)(b-1)/4 odd iff a mod 4==b mod 4==3 */
X
X            /* if a mod 4 == 3 and b mod 4 == 3 */
X            if (((lowdigit(a) & 3) == 3) && ((lowdigit(b) & 3) == 3))
X                sgn = -sgn;
X
X            /* set (a,b) = (b mod a,a) */
X            mpz_set(t,a); mpz_mmod(a,b,t); mpz_set(b,t);
X        } 
X    }
X
X    /* if a == 0 then sgn = 0 */
X    if (uzero(a))
X        sgn=0;
X
X    clear(t); clear(a); clear(b);
X    return (sgn);
X}
X
Xvoid mpz_and(z,x,y) /* not the most efficient way to do this */
XMP_INT *z,*x,*y;
X{
X    int i,sz;
X    sz = imax(x->sz, y->sz);
X    _mpz_realloc(z,(mp_size)sz);
X    for (i=0; i < sz; i++)
X        (z->p)[i] = dg(x,i) & dg(y,i);
X    if (x->sn < 0 && y->sn < 0)
X        z->sn = (-1);
X    else
X        z->sn = 1;
X    if (uzero(z))
X        z->sn = 0;
X}
X
Xvoid mpz_or(z,x,y)  /* not the most efficient way to do this */
XMP_INT *z,*x,*y;
X{
X    int i,sz;
X    sz = imax(x->sz, y->sz);
X    _mpz_realloc(z,(mp_size)sz);
X    for (i=0; i < sz; i++)
X        (z->p)[i] = dg(x,i) | dg(y,i);
X    if (x->sn < 0 || y->sn < 0)
X        z->sn = (-1);
X    else
X        z->sn = 1;
X    if (uzero(z))
X        z->sn = 0;
X}
X
Xvoid mpz_xor(z,x,y)  /* not the most efficient way to do this */
XMP_INT *z,*x,*y;
X{
X    int i,sz;
X    sz = imax(x->sz, y->sz);
X    _mpz_realloc(z,(mp_size)sz);
X    for (i=0; i < sz; i++)
X        (z->p)[i] = dg(x,i) ^ dg(y,i);
X    if ((x->sn <= 0 && y->sn > 0) || (x->sn > 0 && y->sn <=0))
X        z->sn = (-1);
X    else
X        z->sn = 1;
X    if (uzero(z))
X        z->sn = 0;
X}
Xvoid mpz_pow_ui(zz,x,e)
XMP_INT *zz, *x;
Xunsigned long e;
X{
X    MP_INT *t;
X    unsigned long mask = (1L<< (LONGBITS-1));
X    
X    if (e==0) {
X        mpz_set_ui(zz,1L);
X        return;
X    }
X
X    init(t);
X    mpz_set(t,x);
X    for (;!(mask &e); mask>>=1) 
X        ;
X    mask>>=1;
X    for (;mask!=0; mask>>=1) {
X        mpz_mul(t,t,t);
X        if (e & mask)
X            mpz_mul(t,t,x);
X    }
X    mpz_set(zz,t);
X    clear(t);
X}
X
Xvoid mpz_powm(zz,x,ee,n)
XMP_INT *zz,*x,*ee,*n;
X{
X    MP_INT *t,*e;
X    struct is *stack = NULL;
X    int k,i;
X    
X    if (uzero(ee)) {
X        mpz_set_ui(zz,1L);
X        return;
X    }
X
X    if (ee->sn < 0) {
X        return;
X    }
X    init(e); init(t); mpz_set(e,ee);
X
X    for (k=0;!uzero(e);k++,mpz_div_2exp(e,e,1L))
X        push(lowdigit(e) & 1,&stack);
X    k--;
X    i=pop(&stack);
X
X    mpz_mod(t,x,n);  /* t=x%n */
X
X    for (i=k-1;i>=0;i--) {
X        mpz_mul(t,t,t); 
X        mpz_mod(t,t,n);  
X        if (pop(&stack)) {
X            mpz_mul(t,t,x); 
X            mpz_mod(t,t,n);
X        }
X    }
X    mpz_set(zz,t);
X    clear(t);
X}
X
Xvoid mpz_powm_ui(z,x,e,n)
XMP_INT *z,*x,*n; unsigned long e;
X{
X    MP_INT f;
X    mpz_init(&f);
X    mpz_set_ui(&f,e);
X    mpz_powm(z,x,&f,n);
X    mpz_clear(&f);
X}
X    
X    
X
X/* Miller-Rabin */
Xstatic int witness(x,n)
XMP_INT *x, *n;
X{
X    MP_INT *t,*e, *e1;
X    struct is *stack = NULL;
X    int trivflag = 0;
X    int k,i;
X    
X    init(e1); init(e); init(t); mpz_sub_ui(e,n,1L); mpz_set(e1,e);
X
X    for (k=0;!uzero(e);k++,mpz_div_2exp(e,e,1L))
X        push(lowdigit(e) & 1,&stack);
X    k--;
X    i=pop(&stack);
X
X    mpz_mod(t,x,n);  /* t=x%n */
X
X    for (i=k-1;i>=0;i--) {
X        trivflag = !mpz_cmp_ui(t,1L) || !mpz_cmp(t,e1);
X        mpz_mul(t,t,t); mpz_mod(t,t,n);  
X        if (!trivflag && !mpz_cmp_ui(t,1L)) {
X            clear(t); clear(e); clear(e1);
X            return 1;
X        }
X            
X        if (pop(&stack)) {
X            mpz_mul(t,t,x); 
X            mpz_mod(t,t,n);
X        }
X    }
X    if (mpz_cmp_ui(t,1L))  { /* t!=1 */
X        clear(t); clear(e); clear(e1);
X        return 1;
X    }
X    else {
X        clear(t); clear(e); clear(e1);
X        return 0;
X    }
X}
X
Xunsigned long smallp[] = {2,3,5,7,11,13,17};
Xint mpz_probab_prime_p(nn,s)
XMP_INT *nn; int s;
X{   
X    MP_INT *a,*n;
X    int j,k,i;
X    long t;
X
X    if (uzero(nn))
X        return 0;
X    init(n); mpz_abs(n,nn);
X    if (!mpz_cmp_ui(n,1L)) {
X        clear(n);
X        return 0;
X    }
X    init(a);
X    for (i=0;i<sizeof(smallp)/sizeof(unsigned long); i++) {
X        mpz_mod_ui(a,n,smallp[i]);
X        if (uzero(a)) {
X            clear(a); clear(n); return 0;
X        }
X    }
X    _mpz_realloc(a,(mp_size)(n->sz));
X    for (k=0; k<s; k++) {
X
X        /* generate a "random" a */
X            for (i=0; i<n->sz; i++) {
X                t = 0;
X                for (j=0; j < DIGITBITS; j+=8)
X                    t = (t << 8) | mpz_randombyte; 
X                (a->p)[i] = t & LMAX;
X            }
X            a->sn = 1;
X            mpz_mod(a,a,n);
X
X        if (witness(a,n)) {
X            clear(a); clear(n);
X            return 0;
X        }
X    }
X    clear(a);clear(n);
X    return 1;
X}   
X
X
X/* Square roots are found by Newton's method, but since we are working
X * with integers we have to consider carefully the termination conditions.
X * If we are seeking x = floor(sqrt(a)), the iteration is
X *     x_{n+1} = floor ((floor (a/x_n) + x_n)/2) == floor ((a + x_n^2)/(2*x))
X * If eps_n represents the error (exactly, eps_n and sqrt(a) real) in the 
X *  form:
X *     x_n = (1 + eps_n) sqrt(a)
X * then it is easy to show that for a >= 4
X *     if 0 <= eps_n, then either 0 <= eps_{n+1} <= (eps_n^2)/2
X *                         or x_{n+1} = floor (sqrt(a))
X *     if x_n = floor (sqrt(a)), then either x_{n+1} = floor (sqrt(a))
X *                               or x_{n+1} = floor (sqrt(a)) + 1
X * Therefore, if the initial approximation x_0 is chosen so that
X *     0 <= eps_0 <= 1/2
X * then the error remains postive and monotonically decreasing with
X *     eps_n <= 2 ^ (-(2^(n+1) - 1))
X * unless we reach the correct floor(sqrt(a)), after which we may oscillate
X * between it and the value one higher.
X * We choose the number of iterations, k, to guarantee
X *     eps_k sqrt(a) < 1,  so x_k <= floor (sqrt (a)) + 1
X * so finally x_k is either the required answer or one too big (which we
X * check by a simple multiplication and comparison).
X *
X * The calculation of the initial approximation *assumes* DIGITBITS is even.
X * It probably always will be, so for now let's leave the code simple,
X * clear and all reachable in testing!
X */
Xvoid mpz_sqrtrem (root, rem, a)
X    MP_INT *root;
X    MP_INT *rem;
X    MP_INT *a;
X{
X    MP_INT tmp;
X    MP_INT r;
X    mp_limb z;
X    unsigned long j, h;
X    int k;
X
X    if (a->sn < 0)
X        /* Negative operand, nothing correct we can do */
X        return;
X
X    /* Now, a good enough approximation to the root is obtained by writing
X     *     a = z * 2^(2j) + y,  4 <= z <= 15 and 0 <= y < 2^(2j)
X     * then taking
X     *     root = ciel (sqrt(z+1)) * 2^j
X     */
X
X    for (j = a->sz - 1; j != 0 && (a->p)[j] == 0; j--);
X    z = (a->p)[j];
X    if (z < 4) {
X        if (j == 0) {
X            /* Special case for small operand (main interation not valid) */
X            mpz_set_ui (root, (z>0) ? 1L : 0L);
X            if (rem)
X                mpz_set_ui (rem, (z>1) ? z-1L : 0L);
X            return;
X        } else {
X            z = (z << 2) + ((a->p)[j-1] >> (DIGITBITS-2));
X            j = j * DIGITBITS - 2;
X        }
X    } else {
X        j = j * DIGITBITS;
X        while (z & ~0xf) {
X            z >>= 2;
X            j += 2;
X        }
X    }
X    j >>= 1;                            /* z and j now as described above */
X    for (k=1, h=4; h < j+3; k++, h*=2); 
X        /* 2^(k+1) >= j+3, since a < 2^(2j+4) */
X    mpz_init_set_ui (&r, (z>8) ? 4L : 3L);
X    mpz_mul_2exp (&r, &r, j);
X
X#ifdef DEBUG
X    printf ("mpz_sqrtrem: z = %lu, j = %lu, k = %d\n", (long)z, j, k);
X#endif
X
X    /* Main iteration */
X    mpz_init (&tmp);
X    for ( ; k > 0; k--) {
X        mpz_div (&tmp, a, &r);
X        mpz_add (&r, &tmp, &r);
X        mpz_div_2exp (&r, &r, 1L);
X    }
X
X    /* Adjust result (which may be one too big) and set remainder */
X    mpz_mul (&tmp, &r, &r);
X    if (mpz_cmp (&tmp, a) > 0) {
X        mpz_sub_ui (root, &r, 1L);
X        if (rem) {
X            mpz_sub (rem, a, &tmp);
X            mpz_add (rem, rem, &r);
X            mpz_add (rem, rem, &r);
X            mpz_sub_ui (rem, rem, 1L);
X        }
X    } else {
X        mpz_set (root, &r);
X        if (rem)
X            mpz_sub (rem, a, &tmp);
X    }
X
X    mpz_clear (&tmp);
X    mpz_clear (&r);
X}
X
X
Xvoid mpz_sqrt (root, a)
X    MP_INT *root;
X    MP_INT *a;
X{
X    mpz_sqrtrem (root, NULL, a);
X}
END_OF_FILE
if test 39198 -ne `wc -c <'gmp.c'`; then
    echo shar: \"'gmp.c'\" unpacked with wrong size!
fi
# end of 'gmp.c'
fi
if test -f 'gmp.h' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'gmp.h'\"
else
echo shar: Extracting \"'gmp.h'\" \(4203 characters\)
sed "s/^X//" >'gmp.h' <<'END_OF_FILE'
X/*
X * FREE GMP - a public domain implementation of a subset of the 
X *           gmp library
X *
X * I hearby place the file in the public domain.
X *
X * Do whatever you want with this code. Change it. Sell it. Claim you
X *  wrote it. 
X * Bugs, complaints, flames, rants: please send email to 
X *    Mark Henderson <markh@wimsey.bc.ca>
X * I'm already aware that fgmp is considerably slower than gmp
X *
X * CREDITS:
X *  Paul Rouse <par@r-cube.demon.co.uk> - generic bug fixes, mpz_sqrt and
X *    mpz_sqrtrem, and modifications to get fgmp to compile on a system 
X *    with int and long of different sizes (specifically MSDOS,286 compiler)
X *  Also see the file "notes" included with the fgmp distribution, for
X *    more credits.
X *
X * VERSION 1.0.1
X */
X
X#include <stdio.h>
X#include <sys/types.h>
X
X/* for malloc and free */
X#include <stdlib.h>
X
X#ifndef NULL
X#define NULL ((void *)0)
X#endif
X
Xtypedef long mp_limb;
Xtypedef unsigned mp_size;
X
Xtypedef struct mp_int {
X    mp_limb *p;
X    short sn;
X    mp_size sz;
X} MP_INT;
X
X#ifdef __STDC__
X#define PROTO(x) x
X#else
X#define PROTO(x) ()
X#endif
X
Xvoid mpz_init PROTO((MP_INT *s));
Xvoid mpz_init_set PROTO((MP_INT *s, MP_INT *t));
Xvoid mpz_init_set_ui PROTO((MP_INT *s, unsigned long v));
Xvoid mpz_init_set_si PROTO((MP_INT *y, long v));
Xvoid mpz_clear PROTO((MP_INT *s));
Xvoid _mpz_realloc PROTO((MP_INT *x, mp_size size));
Xvoid mpz_set PROTO((MP_INT *y, MP_INT *x));
Xvoid mpz_set_ui PROTO((MP_INT *y, unsigned long v));
Xunsigned long mpz_get_ui PROTO((MP_INT *y));
Xlong mpz_get_si PROTO((MP_INT *y));
Xvoid mpz_set_si PROTO((MP_INT *y, long v));
Xint mpz_cmp PROTO((MP_INT *x, MP_INT *y));
Xvoid mpz_mul PROTO((MP_INT *ww,MP_INT *u, MP_INT *v));
Xvoid mpz_mul_2exp PROTO((MP_INT *z, MP_INT *x, unsigned long e));
Xvoid mpz_div_2exp PROTO((MP_INT *z, MP_INT *x, unsigned long e));
Xvoid mpz_mod_2exp PROTO((MP_INT *z, MP_INT *x, unsigned long e));
Xvoid mpz_add PROTO((MP_INT *zz,MP_INT *x,MP_INT *y));
Xvoid mpz_add_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_mul_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_sub PROTO((MP_INT *z,MP_INT *x, MP_INT *y));
Xvoid mpz_sub_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_div PROTO((MP_INT *q, MP_INT *x, MP_INT *y));
Xvoid mpz_mdiv PROTO((MP_INT *q, MP_INT *x, MP_INT *y));
Xvoid mpz_mod PROTO((MP_INT *r, MP_INT *x, MP_INT *y));
Xvoid mpz_divmod PROTO((MP_INT *q, MP_INT *r, MP_INT *x, MP_INT *y));
Xvoid mpz_mmod PROTO((MP_INT *r, MP_INT *x, MP_INT *y));
Xvoid mpz_mdivmod PROTO((MP_INT *q,MP_INT *r, MP_INT *x, MP_INT *y));
Xvoid mpz_mod_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_mmod_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_div_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_mdiv_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_divmod_ui PROTO((MP_INT *q,MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_mdivmod_ui PROTO((MP_INT *q,MP_INT *x,MP_INT *y, unsigned long n));
Xunsigned int mpz_sizeinbase PROTO((MP_INT *x, int base));
Xchar *mpz_get_str PROTO((char *s,  int base, MP_INT *x));
Xint mpz_set_str PROTO((MP_INT *x, char *s, int base));
Xint mpz_init_set_str PROTO((MP_INT *x, char *s, int base));
Xvoid mpz_random PROTO((MP_INT *x, mp_size size));
Xvoid mpz_random2 PROTO((MP_INT *x, mp_size size));
Xsize_t mpz_size PROTO((MP_INT *x));
Xvoid mpz_abs PROTO((MP_INT *, MP_INT *));
Xvoid mpz_neg PROTO((MP_INT *, MP_INT *));
Xvoid mpz_fac_ui PROTO((MP_INT *, unsigned long));
Xvoid mpz_gcd PROTO((MP_INT *, MP_INT *, MP_INT *));
Xvoid mpz_gcdext PROTO((MP_INT *, MP_INT *, MP_INT *, MP_INT *, MP_INT *));
Xint mpz_jacobi PROTO((MP_INT *, MP_INT *));
Xint mpz_cmp_ui PROTO((MP_INT *, unsigned long));
Xint mpz_cmp_si PROTO((MP_INT *, long));
Xvoid mpz_and PROTO((MP_INT *, MP_INT *, MP_INT *));
Xvoid mpz_or PROTO((MP_INT *, MP_INT *, MP_INT *));
Xvoid mpz_xor PROTO((MP_INT *, MP_INT *, MP_INT *));
Xvoid mpz_pow_ui PROTO((MP_INT *, MP_INT *, unsigned long));
Xvoid mpz_powm PROTO((MP_INT *, MP_INT *, MP_INT *, MP_INT *));
Xvoid mpz_powm_ui PROTO((MP_INT *, MP_INT *, unsigned long, MP_INT *));
Xint mpz_probab_prime_p  PROTO((MP_INT *, int));
Xvoid mpz_sqrtrem PROTO((MP_INT *, MP_INT *, MP_INT *));
Xvoid mpz_sqrt PROTO((MP_INT *, MP_INT *));
END_OF_FILE
if test 4203 -ne `wc -c <'gmp.h'`; then
    echo shar: \"'gmp.h'\" unpacked with wrong size!
fi
# end of 'gmp.h'
fi
echo shar: End of shell archive.
exit 0