/************************************************************************/
/*									*/
/*		Multi-Precise Arithmetic Package			*/
/*									*/
/*	Author: Robert D. Silverman					*/
/*		MITRE Corp.						*/
/*		Burlington Rd.						*/
/*		Bedford, Mass.						*/
/*									*/
/*	Notes:								*/
/*	All multi-precise numbers are stored in radix "RADIX"		*/
/*	They are stored with their length in word 0. The length		*/
/*	includes word 0.						*/
/*	They are stored with most significant digits to the left and	*/
/*	least to the right.						*/
/*									*/
/*	SIZE(x)  is the size of the multi-precise number x		*/
/*	LASTDIGIT(x)	 is the least significant word			*/
/*	FIRST(x) is the most significant				*/
/*									*/
/************************************************************************/

#include <math.h>
#include <stdio.h> 
#define RADIX		1073741824
#define SQRT_RAD	32768
#define MPDIM		50
#define MPCDIM		200
#define MPBITS		30
#define LOGTWO		50
#define EOS		'\0'
#define BLANK		' '

#define SIZE(x)		x[0] 
#define LASTDIGIT(x)	x[1]
#define FIRST(x)	x[SIZE(x)-1]

#define INIT(x)		for (ini=0; ini<MPDIM; ini++) x[ini] = 0;
 

/************************************************************************/
/*									*/
/*			GLOBALS						*/
/*									*/
/************************************************************************/
 
int ini;
/************************************************************************/
/*									*/
/*	routine to copy A into B					*/
/*									*/
/************************************************************************/

mpcopy(a,b)
int a[],b[];
{   /* start of mpcopy */

register int i,d;
d = SIZE(a);
for (i=0; i<d; i++) b[i] = a[i];

}   /* end of mpcopy */



/************************************************************************/
/*									*/
/*	routine to add A and B	: OUTPUT CAN BE INPUT			*/
/*									*/
/************************************************************************/
 
add(a,b,c)
int a[],b[],c[];

{   /* start of add */

register int i,carry,sa;

carry = 0;
sa = SIZE(a);
if (sa == SIZE(b))
	{
	for (i=1; i<sa; i++)
		{
		if ((c[i] = a[i] + b[i] + carry) < RADIX)
		   {
		   carry = 0;
		   }
		else
		   {
		   c[i] -= RADIX;
		   carry = 1;
		   }
		}
	if (carry == 1)
		{
		SIZE(c) = sa + 1;
		FIRST(c) = 1;
		}
	else SIZE(c) = sa;
	}
 

else if (sa < SIZE(b)) add_single(a,b,c);
else add_single(b,a,c);
 
}   /* end of add */

/************************************************************************/
/*									*/
/*	auxiliary routine for ADD only: assumes SIZE(a) < SIZE(b)	*/
/*	 OUTPUT CAN BE INPUT						*/
/************************************************************************/
 
add_single(a,b,c)
int a[],b[],c[];
 
{   /* start of add_single */
register int i,carry,as,bs;
 
carry = 0;
as = SIZE(a);
bs = SIZE(b);
for (i=1; i<as; i++)
	{
	if ((c[i] = a[i] + b[i] + carry) < RADIX)
		{
		carry = 0;
		}
	else
		{
		c[i] -= RADIX;
		carry = 1;
		}
	}
 
for (; i< bs; i++)
	{
	if ((c[i] = b[i] + carry) < RADIX) carry = 0;
	else
		{
		c[i] -= RADIX;
		carry = 1;
		}
	}
 
if (carry == 1)
	{
	SIZE(c) = bs + 1;
	FIRST(c) = 1;
	}
else SIZE(c) = bs;

}   /* end of add_single */
 
/************************************************************************/
/*									*/
/*	routine to compute C = A - B: assumes A >= B			*/
/*	OUTPUT CAN BE INPUT						*/
/************************************************************************/
 
subt(a,b,c)
int a[],b[],c[];

{   /* start of subt */

register int i,borrow,as,bs;
 
borrow = 0;
as = SIZE(a);
bs = SIZE(b);
for (i=1; i<bs; i++)
	{
	if ((c[i] = a[i] - b[i] - borrow) < 0)
		{
		c[i] += RADIX;
		borrow = 1;
		}
	else borrow = 0;
	}

for (; i<= as; i++)
	{
	if ((c[i] = a[i] - borrow) < 0)
		{
		c[i] += RADIX;
		borrow = 1;
		}
	else borrow = 0;
	}
 
SIZE(c) = as; 
i = as-1;
while (i > 1 && c[i] == 0)
	{
	SIZE(c)--;
	i--;
	}

}   /* end of subt */
 
/************************************************************************/
/*									*/
/*	routine to do single precision multiply				*/
/*									*/
/************************************************************************/
 
mul_single(a,b,c)
int a[],b,c[];
 
{   /* start of mul_single */
register int carry,i,d;
int  x[2];

d = SIZE(a);
carry = 0;
for (i=1; i<d; i++)
	{
	muladd(a[i],b,carry,x);		/* x = a[i]*b + carry		*/
	carry = x[1];
	c[i] = x[0];
	}

if (carry == 0) SIZE(c) = d;
else
	{
	SIZE(c) = d+1;
	FIRST(c) = carry;
	}
 
}   /* end of mul_single */


/************************************************************************/
/*									*/
/*	routine to do single precision divide				*/
/*									*/
/************************************************************************/

div_single(a,b,result)
int a[],b,result[];
 
{   /* start of div_single */
register int i,r,as;
int  x[2];
 
r = 0;
as = SIZE(a);
for (i=as-1; i>0; i--)
	{
	divrem(r,a[i],b,x);
	result[i] = x[0];
	r = x[1];
	}
SIZE(result) = as;
i = SIZE(result);
while (i>1 && result[i-1] == 0)
	{
	SIZE(result)--;
	i = SIZE(result);
	}
 
return(r);

}   /* end of div_single */

/************************************************************************/
/*									*/
/*	routine to do single precision remaindering			*/
/*									*/
/************************************************************************/
 
smod(a,b)
int a[],b;
 
{   /* start of smod */
register int i,r,as;
int  x[2];
 
r = 0;
as = SIZE(a);
for (i=as-1; i>0; i--)
	{
	divrem(r,a[i],b,x);
	r = x[1];
	}
 
return(r);

}   /* end of smod */

/************************************************************************/
/*									*/
/*	routine to multiply C = A*B					*/
/*	OUTPUT CAN'T BE INPUT						*/
/************************************************************************/
 
mult(a,b,c)
int a[],b[],c[];
 
{   /* start of mult */
register int carry,i,j,k;
int x[2],da,db;
 
da = SIZE(a);
db = SIZE(b);
for (i=1; i< da; i++) c[i] = 0;		/* initialize output		*/
 
for (j=1; j<db; j++)
	{
	carry = 0;
	for (i=1; i<da; i++)
		{
		k = i + j - 1;
		muladd(a[i],b[j],c[k],x);
		x[0] += carry;
		if (x[0] >= RADIX)
		   {
		   x[0] -= RADIX;
		   x[1]++;
		   }
		carry = x[1];
		c[k] = x[0];
		}
	k = da + j - 1;
	c[k] = carry;
	}
 
SIZE(c) = da + db - 1;
if (FIRST(c) == 0) SIZE(c)--;
 
}   /* end of mult */


/************************************************************************/
/*									*/
/*		routine to perform division				*/
/*									*/
/************************************************************************/
 
div(a,b,q,r)
int a[],b[],q[],r[];
 
{   /* start of div */
register int i,j,k; 
int qhat,d,m,db,t[2],len;
int x[2],borrow,tmp;
int junk1[MPDIM],junk2[MPDIM],junk3[MPDIM];

/*	if divisor is single precision call div_single			*/

if (SIZE(b) < 3)
	{
	r[1] = div_single(a,b[1],q);
	SIZE(r) = 2;
	return;
	}

/*	now we do full blown long division (expensive!)			*/ 
 
m = SIZE(a) - SIZE(b) + 2;
db = SIZE(b);
if (m > 1)
	{   /* outer loop */
	d = RADIX/(1 + FIRST(b));
	mul_single(a,d,junk1);
	mul_single(b,d,junk2);		/* normalize input parameters	*/
	junk2[SIZE(junk2)] = 0;		/* pad with high order zero	*/
	len = SIZE(junk1) - 1;
	if (len == SIZE(a)-1) { len++; junk1[len] = 0; }

	for (j=0; j<m-1; j++)
		{			/* compute qhat first		*/
		if (junk1[len] == junk2[db-1])
		   {
		   qhat = RADIX-1;
		   i = junk1[len-1] + junk2[db-1];
		   }
		else
		   {
		   divrem(junk1[len],junk1[len-1],junk2[db-1],t);
		   qhat = t[0];
		   i = t[1];
		   }
		muladd(qhat,junk2[db-2],0,x);

		while (x[1] > i || (x[1] == i && x[0] > junk1[len-2]))
			{
			qhat--;
			i += junk2[db-1];
			muladd(qhat,junk2[db-2],0,x);
			}
		mul_single(junk2,qhat,junk3);	/* multiply and subtract*/

		if (SIZE(junk2) == SIZE(junk3))
			{
			k = SIZE(junk2);
			SIZE(junk3) = k+1;
			junk3[k] = 0;
			}
		borrow = 0;
		k = len - SIZE(junk2) + 1;

		for (i=1; i<SIZE(junk3); i++)	/* subtract		*/
			{
			tmp = junk1[k] - junk3[i] - borrow;
			if (tmp < 0) { borrow = 1; junk1[k] = tmp+RADIX; }
			else         { borrow = 0; junk1[k] = tmp; }
			k++;
			}

		k = m-j-1;			/* test remainder	*/
		if (borrow == 0)  q[k] = qhat; 
		else
			{			/* add back the borrow	*/
			q[k] = qhat-1;
			borrow = 0;
			k = len - SIZE(junk2) + 1;
			for (i=1; i<SIZE(junk2); i++)
			    {
			    tmp = junk1[k] + junk2[i] + borrow;
			    if (tmp < RADIX) { borrow=0; junk1[k]=tmp; }
			    else { borrow=1; junk1[k] = tmp-RADIX; }
			    k++;
			    }
			}
		len--;
		}   /* end loop on j */

	while (db > 2 && junk1[db-1] == 0) db--;
	SIZE(junk1) = db;
	i = div_single(junk1,d,r);			/* unnormalize	*/
	if (q[m-1] == 0) SIZE(q) = m-1;
	else SIZE(q) = m;
	}
else
	{						/* q = 0, r = a	*/
	mpcopy(a,r);
	q[0] = 2;
	q[1] = 0;
	}

}   /* end of div */

 
/************************************************************************/
/*									*/
/*	routine to convert multi-precision integer to string of		*/
/*	digits of length MPITOC <= size					*/
/*									*/
/************************************************************************/
 
mpitoc(input,str,size)
int input[];
char str[];
int size;
 
{   /* start of mpitoc */

register int i,j,k;
int junk1[MPCDIM],temp,d;
static char digits[10] = { '0','1','2','3','4','5','6','7','8','9' };
char tstr[MPCDIM];

mpcopy(input,junk1);

for (i=0; i<size; i++)
	{
	tstr[i] = EOS;
	str[i] = EOS;
	}
i = 0;
while(i < size)
	{
	d = div_single(junk1,10,junk1);
	tstr[i] = digits[d];
	if ((SIZE(junk1) <= 2 && junk1[1] == 0)) break;
	i++;
	}
 
temp = strlen(tstr);
for (j=0; j<temp; j++) str[j] = tstr[temp-j-1];	/* reverse string	*/
 
return(temp);

}   /* end of mpitoc */

/************************************************************************/
/*									*/
/*	routine to convert string of digits: input to MP	 	*/
/*	intger output.							*/
/*									*/
/************************************************************************/
 
mpctoi(input,output)
int output[];
char input[];

{   /* start of mpctoi */

int d,
    index,
    x[2];

index = 0;
SIZE(x) = 2;
while (input[index] == BLANK) index = index + 1;

SIZE(output) = 2;
for (output[1] = 0; input[index] != EOS; index++)
	{
	d = input[index] - '0';
	mul_single(output,10,output);
	x[1] = d;
	add(output,x,output);
	}
 
}   /* end of mpctoi */

/************************************************************************/
/*									*/
/*	routine to output an MP integer to the terminal			*/
/*	and to an output string						*/
/*									*/
/************************************************************************/
 
write_num(n,numstr)
int n[];
char numstr[];
 
{   /* start of write_num */
register int i;
int nd;
 
nd = mpitoc(n,numstr,MPCDIM);

for (i=0; i<nd; i++)
	{
	printf("%c",numstr[i]);
	}
printf("\n");
 
}   /* end of write_num */

/************************************************************************/
/*									*/
/*	routine to read an MP integer from the terminal			*/
/*									*/
/************************************************************************/
 
get_num(number,numstr)
int number[];
char numstr[];
 
{   /* start of get_num */
 
scanf("%s",numstr);			/* read in string of digits	*/
mpctoi(numstr,number);			/* convert to internal format	*/ 
 
}   /* end of get_num */
 
 

/************************************************************************/
/*									*/
/*	routine to compare two MP integers				*/
/*									*/
/************************************************************************/
 
mpcmp(a,b)
int a[],b[];
 
{    /* start of mpcmp */
register int i;

if (SIZE(a) > SIZE(b)) return(1);
if (SIZE(b) > SIZE(a)) return(-1);
 
for (i=SIZE(a)-1; i > 0; i--)
	{
	if (a[i] > b[i]) return(1);
	if (a[i] < b[i]) return(-1);
	}
 
return(0);

}   /* end of mpcmp */

/************************************************************************/
/*									*/
/*		routine to return [sqrt(n)]				*/
/*									*/
/************************************************************************/
 
mpsqrt(a,b) 
int a[],b[];
 
{   /* start of mpsqrt */

int junkc[MPCDIM],junk3[MPCDIM];
register int i,da,db;
 
da = SIZE(a);
db = 1 + da/2;
SIZE(b) = db;
for (i=1; i<db; i++) b[i] = 0;
if (da % 2 == 0) FIRST(b) = 1+(int)sqrt((double)a[da-1]);
/* da-1, da-2 in these two lines fixed by ssw 4 June 1991 */
else FIRST(b) = 1+(int)(sqrt((double)a[da-1]*(double)RADIX+(double)a[da-2]));
if (FIRST(b) >= RADIX)
	{					/* this is rare		*/
	b[db] = 1;
	FIRST(b) -= RADIX;
	SIZE(b) = db+1;
	}
while(1)					/* loop to infinity	*/
	{
	div(a,b,junkc,junk3);			/* junkc = a/b		*/
	add(b,junkc,junkc);			/* junkc += b		*/
	da = div_single(junkc,2,junkc);		/* divide by 2		*/
	if (mpcmp(junkc,b) == 0) break;		/* converged		*/
	mpcopy(junkc,b);
	}
 
}   /* end of mpsqrt */
 
/************************************************************************/
/*									*/
/*	this routine removes the largest possible power of 2 from N	*/
/*									*/
/************************************************************************/
 
mpodiv(a)
int a[];
 
{   /* start of mpodiv */
register int i,j,index;
int d,accum;
if (SIZE(a) < 2 && a[1] == 0) return(-1);
i=1;
while (a[i] == 0) i++;
i -= 1;

if (i > 0)
	{
	d = SIZE(a) - i;
	for (j=1; j<=d; j++)
		{
		index = j+i;
		a[j] = a[index];
		}
	SIZE(a) = d;
	}

accum = 1; 
j = 0;
d = a[1];
while (d % 2 == 0) { j++; d = d/2; accum = 2*accum; }
if (j > 0) div_single(a,accum,a);

return(MPBITS*i + j);
 
}   /* end of mpodiv */
 
/************************************************************************/
/*									*/
/*		routine to compute GCD(a,b)				*/
/*									*/
/************************************************************************/
 
gcd(a,b,c)
int a[],b[],c[];
 
{   /* start of gcd */
register int i,j;
int junk1[MPCDIM],temp[MPCDIM];
 
if (mpcmp(a,b) == 0)				/* if the two are equal	*/
	{
	mpcopy(a,c);
	return;
	}
 
mpcopy(a,c);
mpcopy(b,junk1);
i = mpodiv(c);
if (SIZE(b) == 2 && b[1] == 0)			/* b is zero		*/
	{
	printf("Attempt to take GCD with 0\n");
	return;
	}

if (SIZE(a) == 2 && a[1] == 0)			/* a is zero		*/
	{
	printf("Attempt to take GCD with 0\n");
	return;
	}

i = mpcmp(c,junk1);
while (i != 0)
	{
	if (i < 0)
		{
		subt(junk1,c,temp);
		mpcopy(temp,junk1);
		j = mpodiv(junk1);
		}
	else
		{
		subt(c,junk1,temp);
		mpcopy(temp,c);
		j = mpodiv(c);
		}
	i = mpcmp(c,junk1);
	}
 
}   /* end of gcd */

/************************************************************************/
/*									*/
/*	routine to compute b^e MOD m for MP integers b,e, and m		*/
/*									*/
/************************************************************************/
 
mpower(b,e,m,a)
int b[],e[],m[],a[];
 
{   /* start of mpower */

int junkc[MPDIM],junk4[MPDIM];
int bits[MPBITS+1],ebit;
register int i,n,j,es;

mpcopy(b,a);
es = SIZE(e)-1;
for (n=es; n>0; n--)
	{
	ebit = e[n];
	i = 0;
	for (j=1; j<=MPBITS; j++)
		{
		i++;
		bits[i] = ebit & 1;
		ebit = ebit/2;
		if (n == SIZE(e)-1 && ebit == 0) { i--; break; }
		}
	for (; i>0; i--)
		{
		mult(a,a,junkc);
		div(junkc,m,junk4,a);
		if (bits[i] != 0)
			{
			mult(a,b,junkc);
			div(junkc,m,junk4,a);
			}
		}
	}
 
}   /* end of mpower */
 
 
/************************************************************************/
/*									*/
/*	ARITHMETIC PRIMITIVES: (A*B + C)/D				*/
/*			       (A*B + C) MOD D				*/
/*									*/
/************************************************************************/
 
mulold(a,b,c,m)
int a,b,c,m[2];
 
{   /* start of muladd */
 
long d;
d = a*b + c;
m[0] = d % RADIX;
m[1] = d/RADIX;
 
}   /* end of muladd */
 
/************************************************************************/
/*									*/
/*	ARITHMETIC PRIMITIVE: q = (A*RADIX + B)/C			*/
/*									*/
/************************************************************************/
 
divold(a,b,q)
int a[2],b,*q;
 
{   /* start of divrem */
long t;
 
t = a[1]*RADIX + a[0];
*q = t/b;
return(t % b);
 
}   /* end of divrem */
  

/************************************************************************/
/*									*/ 
/*	routine to compute A**N for single precision A and N		*/
/*	and a MP result							*/
/*									*/
/************************************************************************/
 
cpower(base,power,result)
int base,
    power,
    result[];
 
{   /* start of cpower */
int temp[MPCDIM],
    temp1[MPCDIM],
    bits[LOGTWO],
    max;
register int i;
 
INIT(result);
INIT(temp);
binary(power,bits,&max);
SIZE(temp) = 2;
temp[1] = base;
SIZE(result) = 2;
result[1] = base;

for (i=max-1; i>=0; i--)
	{
	if (bits[i] == 1)
		{
		mult(result,result,temp1);
		mult(temp1,temp,result);
		}
	else
		{
		mult(result,result,temp1);
		mpcopy(temp1,result);
		}
	}

 
}  /* end of cpower */


/************************************************************************/
/*									*/
/*	routine to return binary representation of an integer		*/
/*									*/
/************************************************************************/
 
binary(n,bits,max)
int n,bits[],*max;
 
{   /* start of binary */
register int i,p,q,r;
 
for (i=0; i<32; i++) bits[i] = 0;
 
p = n;
while (p > 0)
	{
	q = 1;
	r = 0;
	while (q <= p/2)
		{
		q = 2*q;
		r++;
		}
	bits[r] = 1;
	p = p-q;
	}
 
i = 31;
while (bits[i] != 1) i--;
*max = i;
 
}   /* end of binary */
 

/************************************************************************/
/*									*/
/*	routine to display a number in internal format on the sceen	*/
/*									*/
/************************************************************************/
 
internal(num)
int num[];
 
{   /* start of internal */
int i;
 
for (i=SIZE(num)-1; i>0; i--) printf("%d, ",num[i]);
printf("\n");
 
}   /* end of internal */

/************************************************************************/
/*									*/
/*	routine to dump current number into a file			*/
/*									*/
/************************************************************************/
 
fileit(number,ptr)
int number[];
FILE *ptr;
{
register int nd;
char junk[MPCDIM];
 
nd = mpitoc(number,junk,MPCDIM);
fprintf(ptr,"%s\n",junk);
} 

/************************************************************************/
/*									*/
/*			prime testing routine				*/
/*									*/
/************************************************************************/
 
prime_test(number)
int number[];

{   /* start of prime_test */
 
int temp[MPDIM],temp1[MPDIM],base[MPDIM],q[MPDIM],nminus1[MPDIM],
    quo[MPDIM],k,x,j,i,limit;
double mplog( int number[]);

/*	if We assume the generalized Riemann Hypothesis then we can do 	*/
/*	a very fast tes							*/

x = 101;				/* base for exponentiation	*/
mpcopy(number,temp);			/* make temporary copy		*/
 
SIZE(base) = 2;
FIRST(base) = x;			/* use x as base for Fermat test*/

if (LASTDIGIT(number) % 2 == 0) return(1);	/* make sure it's odd	*/
/* if (LASTDIGIT(number)<2 && SIZE(number)==2) return(1); /* check for 0 or 1 */
if (SIZE(number)==2) { /* special code for small numbers */
	k = LASTDIGIT(number)<2;
	if (k<2) return(1); /* check for 0 or 1 */
	if (k<8) return(0); /* check for 3, 5, 7 */
	for (i=3; i<50000; i++) {
		if (k%i == 0) return(1);
		if (i*i > k) return(0);
		}
	}

temp[1]--;				/* use N-1 as the exponent	*/
mpcopy(temp,nminus1);			/* make copy of N-1		*/

mpower(base,temp,number,temp1);		/* base^(N-1) MOD N		*/

/*	check the residue: if it isn't 1 the number is definitely	*/
/*	composite. Otherwise proceed with the test			*/

if (SIZE(temp1) > 2 || LASTDIGIT(temp1) != 1) return(1);

/*	start by removing largest possible power of 2 from N-1		*/
/*	compute k and q for N = 1 + q*2^k				*/

k = 0;
while (LASTDIGIT(temp) % 2 == 0)
	{
	div_single(temp,2,temp1);
	mpcopy(temp1,temp);
	k++;
	}
mpcopy(temp,q);
 
/*	now we try to find a primitive root of N. Once we do compute  	*/
/*	x^((N-1)/p) until = -1 MOD N.					*/
/*	We are guaranteed to find a primitive root quickly because	*/
/*	of the Riemann Hypothesis.					*/
 
limit = (int)mplog(number);		/* need only limit	*/
					/* iterations, maximum		*/
for (i=0; i<limit; i++)
   {
   x = (x*x + 1) % 15511;		/* choose x			*/
   INIT(temp1);
   FIRST(base) = x;
   j = 0;
   mpower(base,q,number,temp1);		/* compute x^q MOD N		*/
   while (j < k)			/* Rabin/Adelman test		*/
	{
	if ( (j == 0 && SIZE(temp1) == 2 && LASTDIGIT(temp1) == 1) ||
	    mpcmp(temp1,nminus1) == 0) goto outer_loop;
	j++;
	mult(temp1,temp1,temp);
	div(temp,number,quo,temp1);
	}
   return(1);				/* N is not prime		*/
 
outer_loop:;
   }
 
return(0);				/* 	N is  prime		*/
 
}   /* end of prime_test */
 

/************************************************************************/
/*									*/
/*	routine to solve x^2 = N mod p, N and p are given		*/
/*	Note: first pass at this routine: use naiive approach		*/
/*									*/
/************************************************************************/

modsqrt(N,p)
int N,p;
 
{   /* start of modsqrt */
int base[MPDIM];
 
int i;
 
if (N < 0) N += p;
 
if (p % 4 == 3)					/* special case 	*/
	{
	SIZE(base) = 2;
	LASTDIGIT(base) = N;
	return(spower(base,(p+1)/4,p));
	}
 
for (i=0; i<p; i++)
	{
	if (i <= SQRT_RAD && (i*i % p) == N) return(i);
	if (i > SQRT_RAD && mod_mult(i,i,p) == N) return(i);
	}
 
return(0);					/* no solution		*/
 
}   /* end of modsqrt */

/************************************************************************/
/*									*/
/*	routine to compute base^exp MOD mod for MP integer base		*/
/*	and single precision exp and mod.	 			*/
/*									*/
/************************************************************************/
 
spower(base,exp,mod)
int base[],exp,mod;
 
{   /* start of spower */

int junk[MPDIM];
int bits[MPBITS+1];
register int i,j,base0,base1;

base1 = div_single(base,mod,junk);
base0 = base1;
 
i = 0;
for (j=1; j<=MPBITS; j++)	/* find binary representation of exp	*/
	{
	i++;
	bits[i] = exp & 1;			/* exp MOD 2		*/
	exp = exp/2;
	if (exp == 0) { i--; break; }
	}
for (; i>0; i--)
	{
	if (base1 <= SQRT_RAD)
		base1 = base1*base1 % mod;
	else
		base1 = mod_mult(base1,base1,mod);
	if (bits[i] != 0)
		{
		if (base1 <= SQRT_RAD && base0 <= SQRT_RAD)
			base1 = base0*base1 % mod;
		else base1 = mod_mult(base1,base0,mod);
		}
	}
 
return(base1);
 
}   /* end of spower */

/************************************************************************/
/*									*/
/*	routine to compute base^exp MOD mod for SP integer base		*/
/*	and single precision exp and mod.	 			*/
/*									*/
/************************************************************************/
 
fastpow(base,exp,mod)
int base,exp,mod;
 
{   /* start of fastpow */

int junk[MPDIM];
int bits[MPBITS+1];
register int i,j,base0,base1;

base1 = base % mod;
base0 = base1;
 
i = 0;
for (j=1; j<=MPBITS; j++)	/* find binary representation of exp	*/
	{
	i++;
	bits[i] = exp & 1;			/* exp MOD 2		*/
	exp = exp/2;
	if (exp == 0) { i--; break; }
	}
for (; i>0; i--)
	{
	if (base1 <= SQRT_RAD)
		base1 = base1*base1 % mod;
	else
		base1 = mod_mult(base1,base1,mod);
	if (bits[i] != 0)
		{
		if (base1 <= SQRT_RAD && base0 <= SQRT_RAD)
			base1 = base0*base1 % mod;
		else base1 = mod_mult(base1,base0,mod);
		}
	}
 
return(base1);
 
}   /* end of fastpow */
 
/************************************************************************/
/*									*/
/*	returns Multiplicative inverse of a mod modulus			*/
/*									*/
/************************************************************************/
 
single_modinv(a,modulus)
int a,modulus;
 
{   /* start of single_modinv */
 
register int ps,ps1,ps2=1,parity,dividend,divisor,rem,q;
 
q = modulus/a;
rem = modulus - a*q;
dividend = a;
divisor = rem;
ps1 = q;
parity = 0;
 
while (rem != 0)
	{
	q = 1;
	rem = dividend - divisor;
	if (rem >= divisor)
	   {
	   q = 2;
	   rem -= divisor;
	   if (rem >= divisor)
	      {
	      q = 3;
	      rem -= divisor;
	      if (rem >= divisor)
		{
		q = 4;
		rem -= divisor;
		if (rem >= divisor)
		   {
		   q = dividend/divisor;
		   rem = dividend - q*divisor;
		   }
		}
	      }
	   }
	ps = q*ps1 + ps2;
	ps2 = ps1;
	ps1 = ps;
	dividend = divisor;
	divisor = rem;
	parity = 1 - parity;
	}
 
if (parity == 0) return(ps2);
else return(modulus-ps2);

}   /* end of single_modinv */

 
/************************************************************************/
/*									*/
/*	routine to return the next prime after N			*/
/*									*/
/************************************************************************/
 
next_prime(N)
int N;
 
{   /* start of next_prime */
register int p,q,r;
 
if (N == 2) return(3);
if (N == 3) return(5);
 
p = N;
while (1)
	{
	p += 2;
	if ((p % 3) == 0) continue;
	if ((p % 5) == 0) continue;
	q = 5;
	r = (int)sqrt((double)p) + 1;
	while (q <= r)
		{
		q += 2;
		if (p % q == 0) goto nextp;
		q += 4;
		if (p % q == 0) goto nextp;
		}
	return(p);
nextp:;
	}
 
}   /* end of next_prime */

double mplog(int num[MPDIM])
{
if (SIZE(num) <= 2)
	{
return log((double)FIRST (num));
	} else {
return (double)(SIZE(num)-3)*log((double)RADIX) +
log( ((double)FIRST (num))*(double)RADIX + (double)(num[SIZE(num)-2]) );
	}
}