home *** CD-ROM | disk | FTP | other *** search
- /*
- * 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);
- }
-