Return to int.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Update and addition of some igcd functions, implementation of modern algorithms.
They are igcdbin(binary alg), igcdbmod(bmod algorithm), igcdEuc(Euclidean),
igcdacc(accelerated alg) and igcdcntl for setting default choice.
Also, update for 32-bit version.

/*
 * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED 
 * All rights reserved.
 * 
 * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
 * non-exclusive and royalty-free license to use, copy, modify and
 * redistribute, solely for non-commercial and non-profit purposes, the
 * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
 * conditions of this Agreement. For the avoidance of doubt, you acquire
 * only a limited right to use the SOFTWARE hereunder, and FLL or any
 * third party developer retains all rights, including but not limited to
 * copyrights, in and to the SOFTWARE.
 * 
 * (1) FLL does not grant you a license in any way for commercial
 * purposes. You may use the SOFTWARE only for non-commercial and
 * non-profit purposes only, such as academic, research and internal
 * business use.
 * (2) The SOFTWARE is protected by the Copyright Law of Japan and
 * international copyright treaties. If you make copies of the SOFTWARE,
 * with or without modification, as permitted hereunder, you shall affix
 * to all such copies of the SOFTWARE the above copyright notice.
 * (3) An explicit reference to this SOFTWARE and its copyright owner
 * shall be made on your publication or presentation in any form of the
 * results obtained by use of the SOFTWARE.
 * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
 * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
 * for such modification or the source code of the modified part of the
 * SOFTWARE.
 * 
 * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
 * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
 * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
 * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
 * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
 * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
 * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
 * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
 * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
 * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
 * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
 * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
 * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
 * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
 * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
 * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
 * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
 *
 * $OpenXM: OpenXM_contrib2/asir2000/builtin/int.c,v 1.8 2000/12/21 02:45:16 murao Exp $ 
*/
#include "ca.h"
#include "parse.h"
#include "base.h"

void Pidiv(), Pirem(), Pigcd(), Pilcm(), Pfac(), Prandom(), Pinv();
void Pup2_inv(),Pgf2nton(), Pntogf2n();
void Pup2_init_eg(), Pup2_show_eg();
void Piqr(), Pprime(), Plprime(), Pinttorat();
void Piand(), Pior(), Pixor(), Pishift();
void Pisqrt();
void iand(), ior(), ixor();
void isqrt();
void Plrandom();
void Pset_upkara(), Pset_uptkara(), Pset_up2kara(), Pset_upfft();
void Pmt_save(), Pmt_load();
void Psmall_jacobi();
void Pdp_set_mpi();
void Pntoint32(),Pint32ton();

void Pigcdbin(), Pigcdbmod(), PigcdEuc(), Pigcdacc(), Pigcdcntl();

void Pihex();
void Pimaxrsh(), Pilen();
void Ptype_t_NB();

struct ftab int_tab[] = {
	{"dp_set_mpi",Pdp_set_mpi,-1},
	{"isqrt",Pisqrt,1},
	{"idiv",Pidiv,2},
	{"irem",Pirem,2},
	{"iqr",Piqr,2},
	{"igcd",Pigcd,-2},
	{"ilcm",Pilcm,2},
	{"up2_inv",Pup2_inv,2},
	{"up2_init_eg",Pup2_init_eg,0},
	{"up2_show_eg",Pup2_show_eg,0},
	{"type_t_NB",Ptype_t_NB,2},
	{"gf2nton",Pgf2nton,1},
	{"ntogf2n",Pntogf2n,1},
	{"set_upkara",Pset_upkara,-1},
	{"set_uptkara",Pset_uptkara,-1},
	{"set_up2kara",Pset_up2kara,-1},
	{"set_upfft",Pset_upfft,-1},
	{"inv",Pinv,2},
	{"inttorat",Pinttorat,3},
	{"fac",Pfac,1},
	{"prime",Pprime,1},
	{"lprime",Plprime,1},
	{"random",Prandom,-1},
	{"lrandom",Plrandom,1},
	{"iand",Piand,2},
	{"ior",Pior,2},
	{"ixor",Pixor,2},
	{"ishift",Pishift,2},
	{"small_jacobi",Psmall_jacobi,2},

	{"igcdbin",Pigcdbin,2},		/* HM@CCUT extension */
	{"igcdbmod",Pigcdbmod,2},	/* HM@CCUT extension */
	{"igcdeuc",PigcdEuc,2},		/* HM@CCUT extension */
	{"igcdacc",Pigcdacc,2},		/* HM@CCUT extension */
	{"igcdcntl",Pigcdcntl,-1},	/* HM@CCUT extension */

	{"mt_save",Pmt_save,1},
	{"mt_load",Pmt_load,1},
	{"ntoint32",Pntoint32,1},
	{"int32ton",Pint32ton,1},
	{0,0,0},
};

static int is_prime_small(unsigned int);
static unsigned int gcd_small(unsigned int,unsigned int);
int TypeT_NB_check(unsigned int, unsigned int);
int mpi_mag;

void Pntoint32(arg,rp)
NODE arg;
USINT *rp;
{
	Q q;
	unsigned int t;

	asir_assert(ARG0(arg),O_N,"ntoint32");
	q = (Q)ARG0(arg);
	if ( !q ) {
		MKUSINT(*rp,0);
		return;
	}
	if ( NID(q)!=N_Q || !INT(q) || PL(NM(q))>1 )
		error("ntoint32 : invalid argument");
	t = BD(NM(q))[0];
	if ( SGN(q) < 0 )
		t = -t;
	MKUSINT(*rp,t);
}

void Pint32ton(arg,rp)
NODE arg;
Q *rp;
{
	int t;

	asir_assert(ARG0(arg),O_USINT,"int32ton");
	t = (int)BDY((USINT)ARG0(arg));
	STOQ(t,*rp);
}

void Pdp_set_mpi(arg,rp)
NODE arg;
Q *rp;
{
	if ( arg ) {
		asir_assert(ARG0(arg),O_N,"dp_set_mpi");
		mpi_mag = QTOS((Q)ARG0(arg));
	}
	STOQ(mpi_mag,*rp);
}

void Psmall_jacobi(arg,rp)
NODE arg;
Q *rp;
{
	Q a,m;
	int a0,m0,s;

	a = (Q)ARG0(arg);
	m = (Q)ARG1(arg);
	asir_assert(a,O_N,"small_jacobi");
	asir_assert(m,O_N,"small_jacobi");
	if ( !a )
		 *rp = ONE;
	else if ( !m || !INT(m) || !INT(a) 
		|| PL(NM(m))>1 || PL(NM(a))>1 || SGN(m) < 0 || EVENN(NM(m)) )
		error("small_jacobi : invalid input");
	else {
		a0 = QTOS(a); m0 = QTOS(m);
		s = small_jacobi(a0,m0);
		STOQ(s,*rp);
	}
}

int small_jacobi(a,m)
int a,m;
{
	int m4,m8,a4,j1,i,s;

	a %= m;
	if ( a == 0 || a == 1 )
		return 1;
	else if ( a < 0 ) {
		j1 = small_jacobi(-a,m);
		m4 = m%4;
		return m4==1?j1:-j1;
	} else { 
		for ( i = 0; a && !(a&1); i++, a >>= 1 );
		if ( i&1 ) {
			m8 = m%8;
			s = (m8==1||m8==7) ? 1 : -1;
		} else
			s = 1;
		/* a, m are odd */
		j1 = small_jacobi(m%a,a);
		m4 = m%4; a4 = a%4;
		s *= (m4==1||a4==1) ? 1 : -1; 
		return j1*s;
	}  
}

void Ptype_t_NB(arg,rp)
NODE arg;
Q *rp;
{
	if ( TypeT_NB_check(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg))))
		*rp = ONE;
	else
		*rp = 0;
}

int TypeT_NB_check(unsigned int m, unsigned int t)
{
	unsigned int p,k,u,h,d;

	if ( !(m%8) )
		return 0;
	p = t*m+1;
	if ( !is_prime_small(p) )
		return 0;
	for ( k = 1, u = 2%p; ; k++ )
		if ( u == 1 )
			break;
		else
			u = (2*u)%p;
	h = t*m/k;
	d = gcd_small(h,m);
	return d == 1 ? 1 : 0;
}

/*
 * a simple prime checker
 * return value: 1  ---  prime number
 *               0  ---  composite number
 */

static int is_prime_small(unsigned int a)
{
	unsigned int b,t,i;

	if ( !(a%2) ) return 0;
	for ( t = a, i = 0; t; i++, t >>= 1 );
	/* b >= sqrt(a) */
	b = 1<<((i+1)/2);

	/* divisibility test by all odd numbers <= b */
	for ( i = 3; i <= b; i += 2 )
		if ( !(a%i) )
			return 0;
	return 1;
}

/*
 * gcd for unsigned int as integers
 * return value: GCD(a,b)
 *
 */


static unsigned int gcd_small(unsigned int a,unsigned int b)
{
	unsigned int t;

	if ( b > a ) {
		t = a; a = b; b = t;
	}
	/* Euclid's algorithm */
	while ( 1 )
		if ( !(t = a%b) ) return b;
		else {
			a = b; b = t;
		}
}

void Pmt_save(arg,rp)
NODE arg;
Q *rp;
{
	int ret;

	ret = mt_save(BDY((STRING)ARG0(arg)));
	STOQ(ret,*rp);
}

void Pmt_load(arg,rp)
NODE arg;
Q *rp;
{
	int ret;

	ret = mt_load(BDY((STRING)ARG0(arg)));
	STOQ(ret,*rp);
}

void Pisqrt(arg,rp)
NODE arg;
Q *rp;
{
	Q a;
	N r;

	a = (Q)ARG0(arg);
	asir_assert(a,O_N,"isqrt");
	if ( !a )
		*rp = 0;
	else if ( SGN(a) < 0 )
		error("isqrt : negative argument");
	else {
		isqrt(NM(a),&r);
		NTOQ(r,1,*rp);
	}
}

void Pidiv(arg,rp)
NODE arg;
Obj *rp;
{
	N q,r;
	Q a;
	Q dnd,dvr;
	
	dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg);
	asir_assert(dnd,O_N,"idiv");
	asir_assert(dvr,O_N,"idiv");
	if ( !dvr )
		error("idiv: division by 0");
	else if ( !dnd )
		*rp = 0;
	else {
		divn(NM(dnd),NM(dvr),&q,&r);
		NTOQ(q,SGN(dnd)*SGN(dvr),a); *rp = (Obj)a;
	}
}

void Pirem(arg,rp)
NODE arg;
Obj *rp;
{
	N q,r;
	Q a;
	Q dnd,dvr;
	
	dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg);
	asir_assert(dnd,O_N,"irem");
	asir_assert(dvr,O_N,"irem");
	if ( !dvr )
		error("irem: division by 0");
	else if ( !dnd )
		*rp = 0;
	else {
		divn(NM(dnd),NM(dvr),&q,&r);
		NTOQ(r,SGN(dnd),a); *rp = (Obj)a;
	}
}

void Piqr(arg,rp)
NODE arg;
LIST *rp;
{
	N q,r;
	Q a,b;
	Q dnd,dvr;
	NODE node1,node2;

	dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg);
	if ( !dvr )
		error("iqr: division by 0");
	else if ( !dnd )
		a = b = 0;
	else if ( OID(dnd) == O_VECT ) {
		iqrv((VECT)dnd,dvr,rp); return;
	} else {
		asir_assert(dnd,O_N,"iqr");
		asir_assert(dvr,O_N,"iqr");
		divn(NM(dnd),NM(dvr),&q,&r);
		NTOQ(q,SGN(dnd)*SGN(dvr),a);
		NTOQ(r,SGN(dnd),b);
	}
	MKNODE(node2,b,0); MKNODE(node1,a,node2); MKLIST(*rp,node1);
}

void Pinttorat(arg,rp)
NODE arg;
LIST *rp;
{
	Q cq,qq,t,u1,v1,r1,nm;
	N m,b,q,r,c,u2,v2,r2;
	NODE node1,node2;
	P p;

	asir_assert(ARG0(arg),O_N,"inttorat");
	asir_assert(ARG1(arg),O_N,"inttorat");
	asir_assert(ARG2(arg),O_N,"inttorat");
	cq = (Q)ARG0(arg); m = NM((Q)ARG1(arg)); b = NM((Q)ARG2(arg));
	if ( !cq ) {
		MKNODE(node2,(pointer)ONE,0); MKNODE(node1,0,node2); MKLIST(*rp,node1);
		return;
	}
	divn(NM(cq),m,&q,&r);
	if ( !r ) {
		MKNODE(node2,(pointer)ONE,0); MKNODE(node1,0,node2); MKLIST(*rp,node1);
		return;
	} else if ( SGN(cq) < 0 ) {
		subn(m,r,&c);
	} else
		c = r;
	u1 = 0; v1 = ONE; u2 = m; v2 = c;
	while ( cmpn(v2,b) >= 0 ) {
		divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
		NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1; 
	}
	if ( cmpn(NM(v1),b) >= 0 )
		*rp = 0;
	else {
		if ( SGN(v1) < 0 ) {
			chsgnp((P)v1,&p); v1 = (Q)p; NTOQ(v2,-1,nm);
		} else
			NTOQ(v2,1,nm);
		MKNODE(node2,v1,0); MKNODE(node1,nm,node2); MKLIST(*rp,node1);
	}
}

void Pigcd(arg,rp)
NODE arg;
Q *rp;
{
	N g;
	Q n1,n2,a;
	
	if ( argc(arg) == 1 ) {
		igcdv((VECT)ARG0(arg),rp);
		return;
	}
	n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
	asir_assert(n1,O_N,"igcd");
	asir_assert(n2,O_N,"igcd");
	if ( !n1 )
		*rp = n2;
	else if ( !n2 )
		*rp = n1;
	else {
		gcdn(NM(n1),NM(n2),&g);
		NTOQ(g,1,a); *rp = a;
	}
}

int comp_n(a,b)
N *a,*b;
{
	return cmpn(*a,*b);
}

void iqrv(a,dvr,rp)
VECT a;
Q dvr;
LIST *rp;
{
	int i,n;
	VECT q,r;
	Q dnd,qi,ri;
	Q *b;
	N qn,rn;
	NODE n0,n1;

	if ( !dvr )
		error("iqrv: division by 0");
	n = a->len; b = (Q *)BDY(a);
	MKVECT(q,n); MKVECT(r,n);
	for ( i = 0; i < n; i++ ) {
		dnd = b[i];
		if ( !dnd ) {
			qi = ri = 0;
		} else {
			divn(NM(dnd),NM(dvr),&qn,&rn);
			NTOQ(qn,SGN(dnd)*SGN(dvr),qi);
			NTOQ(rn,SGN(dnd),ri);
		}
		BDY(q)[i] = (pointer)qi; BDY(r)[i] = (pointer)ri;
	}
	MKNODE(n1,r,0); MKNODE(n0,q,n1); MKLIST(*rp,n0);
}

/*
 * gcd = GCD(a,b), ca = a/g, cb = b/g
 */

void igcd_cofactor(a,b,gcd,ca,cb)
Q a,b;
Q *gcd,*ca,*cb;
{
	N gn,tn;
	
	if ( !a ) {
		if ( !b )
			error("igcd_cofactor : invalid input");
		else {
			*ca = 0;
			*cb = ONE;
			*gcd = b;
		}
	} else if ( !b ) {
		*ca = ONE;
		*cb = 0;
		*gcd = a;
	} else {
		gcdn(NM(a),NM(b),&gn); NTOQ(gn,1,*gcd);
		divsn(NM(a),gn,&tn); NTOQ(tn,SGN(a),*ca);
		divsn(NM(b),gn,&tn); NTOQ(tn,SGN(b),*cb);
	}
}

void igcdv(a,rp)
VECT a;
Q *rp;
{
	int i,j,n,nz;
	N g,gt,q,r;
	N *c;

	n = a->len;
	c = (N *)ALLOCA(n*sizeof(N));
	for ( i = 0; i < n; i++ )
		c[i] = BDY(a)[i]?NM((Q)(BDY(a)[i])):0;
	qsort(c,n,sizeof(N),(int (*) (const void *,const void *))comp_n);
	for ( ; n && ! *c; n--, c++ );

	if ( !n ) {
		*rp = 0; return;
	} else if ( n == 1 ) {
		NTOQ(*c,1,*rp); return;
	}
	gcdn(c[0],c[1],&g);
	for ( i = 2; i < n; i++ ) {
		divn(c[i],g,&q,&r);
		gcdn(g,r,&gt);
		if ( !cmpn(g,gt) ) {
			for ( j = i+1, nz = 0; j < n; j++ ) {
				divn(c[j],g,&q,&r); c[j] = r;
				if ( r )
					nz = 1;
			}
		} else
			g = gt;
	}
	NTOQ(g,1,*rp);
}

void igcdv_estimate(a,rp)
VECT a;
Q *rp;
{
	int n,i,m;
	N s,t,u,g;
	Q *q;

	n = a->len; q = (Q *)a->body;
	if ( n == 1 ) {
		NTOQ(NM(q[0]),1,*rp); return;		
	}

	m = n/2;
	for ( i = 0 , s = 0; i < m; i++ ) {
		addn(s,NM(q[i]),&u); s = u;
	}
	for ( t = 0; i < n; i++ ) {
		addn(t,NM(q[i]),&u); t = u;
	}
	gcdn(s,t,&g); NTOQ(g,1,*rp);
}

void Pilcm(arg,rp)
NODE arg;
Obj *rp;
{
	N g,q,r,l;
	Q n1,n2,a;
	
	n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
	asir_assert(n1,O_N,"ilcm");
	asir_assert(n2,O_N,"ilcm");
	if ( !n1 || !n2 )
		*rp = 0;
	else {
		gcdn(NM(n1),NM(n2),&g); divn(NM(n1),g,&q,&r);
		muln(q,NM(n2),&l); NTOQ(l,1,a); *rp = (Obj)a;
	}
}

void Piand(arg,rp)
NODE arg;
Q *rp;
{
	N g;
	Q n1,n2,a;
	
	n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
	asir_assert(n1,O_N,"iand");
	asir_assert(n2,O_N,"iand");
	if ( !n1 || !n2 )
		*rp = 0;
	else {
		iand(NM(n1),NM(n2),&g);
		NTOQ(g,1,a); *rp = a;
	}
}

void Pior(arg,rp)
NODE arg;
Q *rp;
{
	N g;
	Q n1,n2,a;
	
	n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
	asir_assert(n1,O_N,"ior");
	asir_assert(n2,O_N,"ior");
	if ( !n1 )
		*rp = n2;
	else if ( !n2 )
		*rp = n1;
	else {
		ior(NM(n1),NM(n2),&g);
		NTOQ(g,1,a); *rp = a;
	}
}

void Pixor(arg,rp)
NODE arg;
Q *rp;
{
	N g;
	Q n1,n2,a;
	
	n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
	asir_assert(n1,O_N,"ixor");
	asir_assert(n2,O_N,"ixor");
	if ( !n1 )
		*rp = n2;
	else if ( !n2 )
		*rp = n1;
	else {
		ixor(NM(n1),NM(n2),&g);
		NTOQ(g,1,a); *rp = a;
	}
}

void Pishift(arg,rp)
NODE arg;
Q *rp;
{
	N g;
	Q n1,s,a;
	
	n1 = (Q)ARG0(arg); s = (Q)ARG1(arg);
	asir_assert(n1,O_N,"ixor");
	asir_assert(s,O_N,"ixor");
	if ( !n1 )
		*rp = 0;
	else if ( !s )
		*rp = n1;
	else {
		bshiftn(NM(n1),QTOS(s),&g);
		NTOQ(g,1,a); *rp = a;
	}
}

void isqrt(a,r)
N a,*r;
{
	int k;
	N x,t,x2,xh,quo,rem;

	if ( !a )
		*r = 0;
	else if ( UNIN(a) )
		*r = ONEN;
	else {
		k = n_bits(a); /* a <= 2^k-1 */
		bshiftn(ONEN,-((k>>1)+(k&1)),&x); /* a <= x^2 */
		while ( 1 ) {
			pwrn(x,2,&t);
			if ( cmpn(t,a) <= 0 ) {
				*r = x; return;
			} else {
				if ( BD(x)[0] & 1 )
					addn(x,a,&t);
				else
					t = a;
				bshiftn(x,-1,&x2); divn(t,x2,&quo,&rem);
				bshiftn(x,1,&xh); addn(quo,xh,&x);
			}
		}
	}
}

void iand(n1,n2,r)
N n1,n2,*r;
{
	int d1,d2,d,i;
	N nr;
	int *p1,*p2,*pr;

	d1 = PL(n1); d2 = PL(n2);
	d = MIN(d1,d2); 
	nr = NALLOC(d);
	for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d; i++ )
		pr[i] = p1[i] & p2[i];
	for ( i = d-1; i >= 0 && !pr[i]; i-- );
	if ( i < 0 )
		*r = 0;
	else {
		PL(nr) = i+1; *r = nr;
	}
}

void ior(n1,n2,r)
N n1,n2,*r;
{
	int d1,d2,i;
	N nr;
	int *p1,*p2,*pr;

	if ( PL(n1) < PL(n2) ) {
		nr = n1; n1 = n2; n2 = nr;
	}
	d1 = PL(n1); d2 = PL(n2);
	*r = nr = NALLOC(d1);
	for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ )
		pr[i] = p1[i] | p2[i];
	for ( ; i < d1; i++ )
		pr[i] = p1[i];
	for ( i = d1-1; i >= 0 && !pr[i]; i-- );
	if ( i < 0 )
		*r = 0;
	else {
		PL(nr) = i+1; *r = nr;
	}
}

void ixor(n1,n2,r)
N n1,n2,*r;
{
	int d1,d2,i;
	N nr;
	int *p1,*p2,*pr;

	if ( PL(n1) < PL(n2) ) {
		nr = n1; n1 = n2; n2 = nr;
	}
	d1 = PL(n1); d2 = PL(n2);
	*r = nr = NALLOC(d1);
	for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ )
		pr[i] = p1[i] ^ p2[i];
	for ( ; i < d1; i++ )
		pr[i] = p1[i];
	for ( i = d1-1; i >= 0 && !pr[i]; i-- );
	if ( i < 0 )
		*r = 0;
	else {
		PL(nr) = i+1; *r = nr;
	}
}

void Pup2_init_eg(rp)
Obj *rp;
{
	up2_init_eg();
	*rp = 0;
}

void Pup2_show_eg(rp)
Obj *rp;
{
	up2_show_eg();
	*rp = 0;
}

void Pgf2nton(arg,rp)
NODE arg;
Q *rp;
{
	if ( !ARG0(arg) )
		*rp = 0;
	else
		up2ton(((GF2N)ARG0(arg))->body,rp);
}

void Pntogf2n(arg,rp)
NODE arg;
GF2N *rp;
{
	UP2 t;

	ntoup2(ARG0(arg),&t);
	MKGF2N(t,*rp);
}

void Pup2_inv(arg,rp)
NODE arg;
P *rp;
{
	UP2 a,b,t;

	ptoup2(ARG0(arg),&a);
	ptoup2(ARG1(arg),&b);
	invup2(a,b,&t);
	up2top(t,rp);
}

void Pinv(arg,rp)
NODE arg;
Num *rp;
{
	Num n;
	Q mod;
	MQ r;	
	int inv;

	n = (Num)ARG0(arg); mod = (Q)ARG1(arg);
	asir_assert(n,O_N,"inv");
	asir_assert(mod,O_N,"inv");
	if ( !n || !mod )
		error("inv: invalid input");
	else
		switch ( NID(n) ) {
			case N_Q:
				invl((Q)n,mod,(Q *)rp);
				break;
			case N_M:
				inv = invm(CONT((MQ)n),QTOS(mod));
				STOMQ(inv,r);
				*rp = (Num)r;
				break;
			default:
				error("inv: invalid input");
		}
}

void Pfac(arg,rp)
NODE arg;
Q *rp;
{ 
	asir_assert(ARG0(arg),O_N,"fac");
	factorial(QTOS((Q)ARG0(arg)),rp); 
}

void Plrandom(arg,rp)
NODE arg;
Q *rp;
{
	N r;

	asir_assert(ARG0(arg),O_N,"lrandom");
	randomn(QTOS((Q)ARG0(arg)),&r);
	NTOQ(r,1,*rp);
}

void Prandom(arg,rp)
NODE arg;
Q *rp;
{
	unsigned int r;

#if 0
#if defined(_PA_RISC1_1)
	r = mrand48()&BMASK;
#else
	if ( arg )
		srandom(QTOS((Q)ARG0(arg)));
	r = random()&BMASK;
#endif
#endif
	if ( arg )
		mt_sgenrand(QTOS((Q)ARG0(arg)));
	r = mt_genrand();
	UTOQ(r,*rp);
}

#if defined(VISUAL)
void srandom(unsigned int);

static unsigned int R_Next;

unsigned int random() {
        if ( !R_Next )
                R_Next = 1;
        return R_Next = (R_Next * 1103515245 + 12345);
}

void srandom(s)
unsigned int s;
{
		if ( s )
			R_Next = s;
        else if ( !R_Next )
            R_Next = 1;
}
#endif

void Pprime(arg,rp)
NODE arg;
Q *rp;
{
	int i;

	asir_assert(ARG0(arg),O_N,"prime");
	i = QTOS((Q)ARG0(arg));
	if ( i < 0 || i >= 1900 )
		*rp = 0;
	else
		STOQ(sprime[i],*rp);
}

void Plprime(arg,rp)
NODE arg;
Q *rp;
{
	int i;

	asir_assert(ARG0(arg),O_N,"lprime");
	i = QTOS((Q)ARG0(arg));
	if ( i < 0 || i >= 1000 )
		*rp = 0;
	else
		STOQ(lprime[i],*rp);
}

extern int up_kara_mag, up_tkara_mag, up_fft_mag;

void Pset_upfft(arg,rp)
NODE arg;
Q *rp;
{
	if ( arg ) {
		asir_assert(ARG0(arg),O_N,"set_upfft");
		up_fft_mag = QTOS((Q)ARG0(arg));
	}
	STOQ(up_fft_mag,*rp);
}

void Pset_upkara(arg,rp)
NODE arg;
Q *rp;
{
	if ( arg ) {
		asir_assert(ARG0(arg),O_N,"set_upkara");
		up_kara_mag = QTOS((Q)ARG0(arg));
	}
	STOQ(up_kara_mag,*rp);
}

void Pset_uptkara(arg,rp)
NODE arg;
Q *rp;
{
	if ( arg ) {
		asir_assert(ARG0(arg),O_N,"set_uptkara");
		up_tkara_mag = QTOS((Q)ARG0(arg));
	}
	STOQ(up_tkara_mag,*rp);
}

extern int up2_kara_mag;

void Pset_up2kara(arg,rp)
NODE arg;
Q *rp;
{
	if ( arg ) {
		asir_assert(ARG0(arg),O_N,"set_up2kara");
		up2_kara_mag = QTOS((Q)ARG0(arg));
	}
	STOQ(up2_kara_mag,*rp);
}

void Pigcdbin(arg,rp)
NODE arg;
Obj *rp;
{
	N g;
	Q n1,n2,a;
	
	n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
	asir_assert(n1,O_N,"igcd");
	asir_assert(n2,O_N,"igcd");
	if ( !n1 )
		*rp = (Obj)n2;
	else if ( !n2 )
		*rp = (Obj)n1;
	else {
		gcdbinn(NM(n1),NM(n2),&g);
		NTOQ(g,1,a); *rp = (Obj)a;
	}
}

void Pigcdbmod(arg,rp)
NODE arg;
Obj *rp;
{
	N g;
	Q n1,n2,a;
	
	n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
	asir_assert(n1,O_N,"igcdbmod");
	asir_assert(n2,O_N,"igcdbmod");
	if ( !n1 )
		*rp = (Obj)n2;
	else if ( !n2 )
		*rp = (Obj)n1;
	else {
		gcdbmodn(NM(n1),NM(n2),&g);
		NTOQ(g,1,a); *rp = (Obj)a;
	}
}

void Pigcdacc(arg,rp)
NODE arg;
Obj *rp;
{
	N g;
	Q n1,n2,a;
	
	n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
	asir_assert(n1,O_N,"igcdacc");
	asir_assert(n2,O_N,"igcdacc");
	if ( !n1 )
		*rp = (Obj)n2;
	else if ( !n2 )
		*rp = (Obj)n1;
	else {
		gcdaccn(NM(n1),NM(n2),&g);
		NTOQ(g,1,a); *rp = (Obj)a;
	}
}

void PigcdEuc(arg,rp)
NODE arg;
Obj *rp;
{
	N g;
	Q n1,n2,a;
	
	n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
	asir_assert(n1,O_N,"igcdbmod");
	asir_assert(n2,O_N,"igcdbmod");
	if ( !n1 )
		*rp = (Obj)n2;
	else if ( !n2 )
		*rp = (Obj)n1;
	else {
		gcdEuclidn(NM(n1),NM(n2),&g);
		NTOQ(g,1,a); *rp = (Obj)a;
	}
}

extern int igcd_algorithm;
    /*	== 0 : Euclid,
     *	== 1 : binary,
     *	== 2 : bmod,
     *	>= 3 : (Weber's accelerated)/(Jebelean's generalized binary) algorithm,
     */
extern int igcd_thre_inidiv;
    /*
     *  In the non-Euclidean algorithms, if the ratio of the lengths (number
     *  of words) of two integers is >= igcd_thre_inidiv, we first perform
     *  remainder calculation.
     *  If == 0, this remainder calculation is not performed.
     */
extern int igcdacc_thre;
    /*
     *  In the accelerated algorithm, if the bit-lengths of two integers is
     *  > igcdacc_thre, "bmod" reduction is done.
     */

void Pigcdcntl(arg,rp)
NODE arg;
Obj *rp;
{
	Obj p;
	Q a;
	int k, m;

	if ( arg ) {
		p = (Obj)ARG0(arg);
		if ( !p ) {
			igcd_algorithm = 0;
			*rp = p;
			return;
		} else if ( OID(p) == O_N ) {
			k = QTOS((Q)p);
			a = (Q)p;
			if ( k >= 0 ) igcd_algorithm = k;
			else if ( k == -1 ) {
			ret_thre:
				k = - igcd_thre_inidiv - igcdacc_thre*10000;
				STOQ(k,a);
				*rp = (Obj)a;
				return;
			} else {
				if ( (m = (-k)%10000) != 0 ) igcd_thre_inidiv = m;
				if ( (m = -k/10000) != 0 ) igcdacc_thre = m;
				goto ret_thre;
			}
		} else if ( OID(p) == O_STR ) {
			char *n = BDY((STRING) p);

			if ( !strcmp( n, "binary" ) || !strcmp( n, "Binary" )
			     || !strcmp( n, "bin" ) || !strcmp( n, "Bin" ) )
				k = igcd_algorithm = 1;
			else if ( !strcmp( n, "bmod" ) || !strcmp( n, "Bmod" ) )
				igcd_algorithm = 2;
			else if ( !strcmp( n, "euc" ) || !strcmp( n, "Euc" )
				  || !strcmp( n, "euclid" ) || !strcmp( n, "Euclid" ) )
				igcd_algorithm = 0;
			else if ( !strcmp( n, "acc" ) || !strcmp( n, "Acc" )
			     || !strcmp( n, "gen" ) || !strcmp( n, "Gen" )
			     || !strcmp( n, "genbin" ) || !strcmp( n, "GenBin" ) )
				igcd_algorithm = 3;
			else error( "igcdhow : invalid algorithm specification" );
		}
	}
	STOQ(igcd_algorithm,a);
	*rp = (Obj)a;
	return;
}