conservatory
/
stk
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
stk/Mp/fgmp-1.0b5/gmp.c

1752 lines
38 KiB
C
Raw Normal View History

Commit of 3.1.1 version
1996-09-27 06:29:02 -04:00
/*
* 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);
}