===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/Mgfs.c,v
retrieving revision 1.5
retrieving revision 1.15
diff -u -p -r1.5 -r1.15
--- OpenXM_contrib2/asir2000/engine/Mgfs.c	2001/06/25 05:30:49	1.5
+++ OpenXM_contrib2/asir2000/engine/Mgfs.c	2002/12/18 06:15:40	1.15
@@ -1,15 +1,21 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.4 2001/06/25 04:11:42 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.14 2002/09/27 08:40:49 noro Exp $ */
 
 #include "ca.h"
+#include "inline.h"
 
+extern int current_gfs_p;
 extern int up_kara_mag, current_gfs_q1;
 extern int *current_gfs_plus1;
 
-void mulssfum(UM,int,UM);
-void kmulsfummain(UM,UM,UM);
+#if defined(__GNUC__)
+#define INLINE inline
+#elif defined(VISUAL)
+#define INLINE __inline
+#else
+#define INLINE
+#endif
 
-inline int _ADDSF(a,b)
-int a,b;
+INLINE int _ADDSF(int a,int b)
 {
 	if ( !a )
 		return b;
@@ -17,6 +23,18 @@ int a,b;
 		return a;
 
 	a = IFTOF(a); b = IFTOF(b);
+
+	if ( !current_gfs_ntoi ) {
+		a = a+b-current_gfs_q;
+		if ( a == 0 )
+			return 0;
+		else {
+			if ( a < 0 )
+				a += current_gfs_q;
+			return FTOIF(a);
+		}
+	}
+
 	if ( a > b ) {
 		/* tab[a]+tab[b] = tab[b](tab[a-b]+1) */
 		a = current_gfs_plus1[a-b];
@@ -42,12 +60,17 @@ int a,b;
 	}
 }
 
-inline int _MULSF(a,b)
-int a,b;
+INLINE int _MULSF(int a,int b)
 {
+	int c;
+
 	if ( !a || !b )
 		return 0;
-	else {
+	else if ( !current_gfs_ntoi ) {
+		a = IFTOF(a); b = IFTOF(b);
+		DMAR(a,b,0,current_gfs_q,c);
+		return FTOIF(c);
+	} else {
 		a = IFTOF(a) + IFTOF(b);
 		if ( a >= current_gfs_q1 )
 			a -= current_gfs_q1;
@@ -55,8 +78,7 @@ int a,b;
 	}
 }
 
-void addsfum(p1,p2,pr)
-UM p1,p2,pr;
+void addsfum(UM p1,UM p2,UM pr)
 {
 	int *c1,*c2,*cr,i,dmax,dmin;
 		
@@ -83,8 +105,7 @@ UM p1,p2,pr;
 		DEG(pr) = dmax;
 }
 
-void subsfum(p1,p2,pr)
-UM p1,p2,pr;
+void subsfum(UM p1,UM p2,UM pr)
 {
 	int *c1,*c2,*cr,i;
 	int dmax,dmin;
@@ -119,8 +140,7 @@ UM p1,p2,pr;
 		DEG(pr) = dmax;
 }
 		
-void gcdsfum(p1,p2,pr)
-UM p1,p2,pr;
+void gcdsfum(UM p1,UM p2,UM pr)
 {
 	int inv;
 	UM t1,t2,q,tum;
@@ -144,8 +164,8 @@ UM p1,p2,pr;
 		mulssfum(t2,inv,pr);
 	}
 }
-void mulsfum(p1,p2,pr)
-UM p1,p2,pr;
+
+void mulsfum(UM p1,UM p2,UM pr)
 {
 	int *pc1,*pcr;
 	int *c1,*c2,*cr;
@@ -166,9 +186,7 @@ UM p1,p2,pr;
 	DEG(pr) = d1 + d2;
 }
 
-void mulssfum(p,n,pr)
-int n;
-UM p,pr;
+void mulssfum(UM p,int n,UM pr)
 {
 	int *sp,*dp;
 	int i;
@@ -178,8 +196,7 @@ UM p,pr;
 		*dp = _MULSF(*sp,n);
 }
 	
-void kmulsfum(n1,n2,nr)
-UM n1,n2,nr;
+void kmulsfum(UM n1,UM n2,UM nr)
 {
 	UM n,t,s,m,carry;
 	int d,d1,d2,len,i,l;
@@ -226,8 +243,7 @@ UM n1,n2,nr;
 	bcopy((char *)r0,(char *)COEF(nr),l*sizeof(int));
 }
 
-void kmulsfummain(n1,n2,nr)
-UM n1,n2,nr;
+void kmulsfummain(UM n1,UM n2,UM nr)
 {
 	int d1,d2,h,len;
 	UM n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
@@ -261,8 +277,7 @@ UM n1,n2,nr;
 		copyum(t1,nr);
 }
 
-int divsfum(p1,p2,pq)
-UM p1,p2,pq;
+int divsfum(UM p1,UM p2,UM pq)
 {
 	int *pc1,*pct;
 	int *c1,*c2,*ct;
@@ -298,21 +313,19 @@ UM p1,p2,pq;
 	return i;
 }
 
-void diffsfum(f,fd)
-UM f,fd;
+void diffsfum(UM f,UM fd)
 {
 	int *dp,*sp;
 	int i;
 
 	for ( i = DEG(f), dp = COEF(fd)+i-1, sp = COEF(f)+i; 
 		i >= 1; i--, dp--, sp-- ) {
-		*dp = _MULSF(*sp,_itosf(i));
+		*dp = _MULSF(*sp,_itosf(i%current_gfs_p));
 	}
 	degum(fd,DEG(f) - 1);
 }
 
-void monicsfum(f)
-UM f;
+void monicsfum(UM f)
 {
 	int *sp;
 	int i,inv;
@@ -323,126 +336,93 @@ UM f;
 		*sp = _MULSF(*sp,inv);
 }
 
-void addsfarray(int,int *,int *);
-void mulsfarray_trunc(int,int *,int *,int *);
-
-void mulsflum(n,f1,f2,fr)
-int n;
-LUM f1,f2,fr;
+int compsfum(UM a,UM b)
 {
-	int max;
-	int i,j,**p1,**p2,*px;
-	int *w,*w1,*w2;
+	int i,da,db;
 
-	p1 = (int **)COEF(f1); p2 = (int **)COEF(f2);
-	w = W_ALLOC(2*(n+1)); w1 = W_ALLOC(DEG(f1)); w2 = W_ALLOC(DEG(f2));
-	for ( i = DEG(f1); i >= 0; i-- ) {
-		for ( j = n - 1, px = p1[i]; ( j >= 0 ) && ( px[j] == 0 ); j-- );
-		w1[i] = ( j == -1 ? 0 : 1 );
+	if ( !a )
+		return !b?0:1;
+	else if ( !b )
+		return 1;
+	else if ( (da = DEG(a)) > (db = DEG(b)) )
+		return 1;
+	else if ( da < db )
+		return -1;
+	else {
+		for ( i = da; i >= 0 && COEF(a)[i] == COEF(b)[i]; i-- );
+		if ( i < 0 )
+			return 0;
+		else if ( (unsigned int)COEF(a)[i] > (unsigned int)COEF(b)[i] )
+			return 1;
+		else
+			return -1;
 	}
-	for ( i = DEG(f2); i >= 0; i-- ) {
-		for ( j = n - 1, px = p2[i]; ( j >= 0 ) && ( px[j] == 0 ); j-- );
-		w2[i] = ( j == -1 ? 0 : 1 );
-	}
-	for ( j = DEG(fr) = DEG(f1) + DEG(f2); j >= 0; j-- ) {
-		for ( i = n - 1, px = COEF(fr)[j]; i >= 0; i-- ) 
-			px[i] = 0;
-		for ( max = MIN(DEG(f1),j), i = MAX(0,j-DEG(f2)); i <= max; i++ ) 
-			if ( w1[i] != 0 && w2[j - i] != 0 ) {
-				mulsfarray_trunc(n,p1[i],p2[j - i],w); addsfarray(n,w,px);
-			}
-	}
 }
 
+void addsfarray(int,int *,int *);
+void mulsfarray_trunc(int,int *,int *,int *);
+
 /* f1 = f1->c[0]+f1->c[1]*y+...,  f2 = f2->c[0]+f2->c[1]*y+... mod y^n */
 
-void mulsfbm(bound,f1,f2,fr)
-int bound;
-BM f1,f2,fr;
+void mulsfbm(BM f1,BM f2,BM fr)
 {
 	UM mul,t,s;
-	int i,j,h,d1,d2;
+	int i,j,d1,d2,dy;
 
-	if ( DEG(f1) < bound || DEG(f2) < bound )
-		error("mulsfbm : invalid input");
+	dy = MIN(DEG(f1),DEG(f2));
 
-	d1 = DEG(COEF(f1)[0]);
-	for ( i = 1; i < bound; i++ )
-		d1 = MAX(DEG(COEF(f1)[i]),d1);
-	d2 = DEG(COEF(f2)[0]);
-	for ( i = 1; i < bound; i++ )
-		d2 = MAX(DEG(COEF(f2)[i]),d2);
+	d1 = degbm(f1);
+	d2 = degbm(f2);
 	t = W_UMALLOC(d1+d2);
 	s = W_UMALLOC(d1+d2);
 
-	for ( i = 0; i < bound; i++ ) {
+	for ( i = 0; i <= dy; i++ ) {
 		mul = COEF(f2)[i];
 		if ( DEG(mul) >= 0 )
-		for ( j = 0; i+j < bound; j++ ) {
+		for ( j = 0; i+j <= dy; j++ ) {
 			if ( COEF(f1)[j] ) {
-				mulsfum(COEF(f1)[j],mul,t);
+				kmulsfum(COEF(f1)[j],mul,t);
 				addsfum(t,COEF(fr)[i+j],s);
 				cpyum(s,COEF(fr)[i+j]);
 			}
 		}
 	}
-	DEG(fr) = bound;
+	DEG(fr) = dy;
 }
 
+int degbm(BM f)
+{
+	int d,i,dy;
+
+	dy = DEG(f);
+	d = DEG(COEF(f)[0]);
+	for ( i = 1; i <= dy; i++ )
+		d = MAX(DEG(COEF(f)[i]),d);
+	return d;
+}
+
 /* g += f */
 
-void addtosfbm(bound,f,g)
-int bound;
-BM f,g;
+void addtosfbm(BM f,BM g)
 {
-	int i,d1,d2;
+	int i,d1,d2,dy;
 	UM t;
 
-	d1 = DEG(COEF(f)[0]);
-	for ( i = 1; i < bound; i++ )
-		d1 = MAX(DEG(COEF(f)[i]),d1);
-	d2 = DEG(COEF(g)[0]);
-	for ( i = 1; i < bound; i++ )
-		d2 = MAX(DEG(COEF(g)[i]),d2);
+	dy = DEG(g);
+	if ( DEG(f) > dy )
+		error("addtosfbm : invalid input");
+	dy = MIN(DEG(f),dy);
+	d1 = degbm(f);
+	d2 = degbm(g);
 	t = W_UMALLOC(MAX(d1,d2));
-	for ( i = 0; i < bound; i++ ) {
+	for ( i = 0; i <= dy; i++ ) {
 		addsfum(COEF(f)[i],COEF(g)[i],t);
 		cpyum(t,COEF(g)[i]);
 	}
 }
 
-void addsfarray(n,a1,a2)
-int n;
-int *a1,*a2;
+void eucsfum(UM f1,UM f2,UM a,UM b)
 {
-	int i;
-
-	for ( i = 0; i < n; i++, a1++, a2++ )
-		if ( *a1 )
-			*a2 = _ADDSF(*a1,*a2);
-}
-
-void mulsfarray_trunc(n,a1,a2,r)
-int n;
-int *a1,*a2,*r;
-{
-	int mul,i,j;
-	int *c1,*c2,*cr;
-	int *pc1,*pc2,*pcr;
-
-	bzero(r,n*sizeof(int));
-	c1 = a1; c2 = a2; cr = r;
-	for ( i = 0; i < n; i++, cr++ ) {
-		if ( mul = *c2++ )
-		for ( j = 0, pc1 = c1, pcr = cr; j+i < n; j++, pc1++, pcr++ )
-			if ( *pc1 )
-				*pcr = _ADDSF(_MULSF(*pc1,mul),*pcr);
-	}
-}
-
-void eucsfum(f1,f2,a,b)
-UM f1,f2,a,b;
-{
 	UM g1,g2,a1,a2,a3,wm,q,tum;
 	int d,dr;
 
@@ -475,10 +455,7 @@ UM f1,f2,a,b;
 
 void shiftsfum(UM,int,UM);
 
-void shiftsflum(n,f,ev)
-int n;
-LUM f;
-int ev;
+void shiftsflum(int n,LUM f,int ev)
 {
 	int d,i,j;
 	UM w,w1;
@@ -502,10 +479,7 @@ int ev;
 
 /* f(x) -> g(x) = f(x+a) */
 
-void shiftsfum(f,a,g)
-UM f;
-int a;
-UM g;
+void shiftsfum(UM f,int a,UM g)
 {
 	int n,i;
 	UM pwr,xa,w,t;
@@ -542,86 +516,262 @@ UM g;
 
 /* f(y) -> f(y+a) */
 
-void shiftsfbm(bound,f,a)
-int bound;
-BM f;
-int a;
+void shiftsfbm(BM f,int a)
 {
-	int i,j,d;
+	int i,j,d,dy;
 	UM pwr,ya,w,t,s;
 	UM *c;
 
-	if ( bound <= 0 )
-		return;
-	else {
-		c = COEF(f);
-		d = DEG(c[0]);
-		for ( i = 1; i < bound; i++ )
-			d = MAX(DEG(c[i]),d);
+	dy = DEG(f);
+	c = COEF(f);
+	d = degbm(f);
 
-		w = W_UMALLOC(d);
-		t = W_UMALLOC(d);
-		s = W_UMALLOC(bound);
+	w = W_UMALLOC(d);
+	t = W_UMALLOC(d);
+	s = W_UMALLOC(dy);
 
-		/* pwr = 1 */
-		pwr = W_UMALLOC(bound); DEG(pwr) = 0; COEF(pwr)[0] = _onesf();
+	/* pwr = 1 */
+	pwr = W_UMALLOC(dy); DEG(pwr) = 0; COEF(pwr)[0] = _onesf();
 
-		/* ya = y+a */
-		ya = W_UMALLOC(1); DEG(ya) = 1; 
-		COEF(ya)[0] = a; COEF(ya)[1] = _onesf();
+	/* ya = y+a */
+	ya = W_UMALLOC(1); DEG(ya) = 1; 
+	COEF(ya)[0] = a; COEF(ya)[1] = _onesf();
 
-		for ( i = 0; i < bound; i++ ) {
-			/* c[i] does not change */
-			for ( j = 0; j < i; j++ ) {
-				mulssfum(c[i],COEF(pwr)[j],w);
-				addsfum(w,c[j],t); cpyum(t,c[j]);
-			}
-			if ( i < bound-1 ) {
-				mulsfum(pwr,ya,s); cpyum(s,pwr);
-			}
+	for ( i = 0; i <= dy; i++ ) {
+		/* c[i] does not change */
+		for ( j = 0; j < i; j++ ) {
+			mulssfum(c[i],COEF(pwr)[j],w);
+			addsfum(w,c[j],t); cpyum(t,c[j]);
 		}
+		if ( i < dy ) {
+			mulsfum(pwr,ya,s); cpyum(s,pwr);
+		}
 	}
 }
 
-/* clear the body of f */
+void clearbm(int n,BM f)
+{
+	int i,dy;
+	UM *c;
 
-void clearsflum(bound,n,f)
-int bound,n;
-LUM f;
+	dy = DEG(f);
+	for ( c = COEF(f), i = 0; i <= dy; i++ )
+		clearum(c[i],n);
+}
+
+static void create_bmono(P c,V x,int i,V y,int j,P *mono);
+
+struct lb {
+	int pos,len;
+	int *r;
+	int *hist;
+};
+
+static NODE insert_lb(NODE g,struct lb *a)
 {
-	int **c;
-	int i;
+	NODE prev,cur,n;
 
-	DEG(f) = 0;
-	for ( c = COEF(f), i = 0; i <= n; i++ )
-		bzero(c[i],(bound+1)*sizeof(int));
+	prev = 0; cur = g;
+	while ( cur ) {
+		if ( a->pos < ((struct lb *)BDY(cur))->pos ) {
+			MKNODE(n,a,cur);
+			if ( !prev )
+				return n;
+			else {
+				NEXT(prev) = n;
+				return g;
+			}
+		} else {
+			prev = cur;
+			cur = NEXT(cur);
+		}
+	}
+	MKNODE(n,a,0);
+	NEXT(prev) = n;
+	return g;
 }
 
-void clearsfbm(bound,n,f)
-int bound,n;
-BM f;
+static void lnf(int *r,int *h,int n,int len,NODE g)
 {
-	int i;
-	UM *c;
+	struct lb *t;
+	int pos,i,j,len1,c;
+	int *r1,*h1;
 
-	DEG(f) = bound;
-	for ( c = COEF(f), i = 0; i < bound; i++ )
-		clearum(c[i],n);
+	for ( ; g; g = NEXT(g) ) {
+		t = (struct lb *)BDY(g);
+		pos = t->pos;
+		if ( c = r[pos] ) {
+			c = _chsgnsf(c);
+			r1 = t->r;
+			h1 = t->hist;
+			len1 = t->len;
+			for ( i = pos; i < n; i++ )
+				if ( r1[i] )
+					r[i] = _ADDSF(r[i],_MULSF(r1[i],c));
+			for ( i = 0; i < len1; i++ )
+				if ( h1[i] )
+					h[i] = _ADDSF(h[i],_MULSF(h1[i],c));
+		}
+	}
+	for ( i = 0; i < n && !r[i]; i++ );
+	if ( i < n ) {
+		c = _invsf(r[i]);
+		for ( j = i; j < n; j++ )
+			if ( r[j] )
+				r[j] = _MULSF(r[j],c);	
+		for ( j = 0; j < len; j++ )
+			if ( h[j] ) 
+				h[j] = _MULSF(h[j],c);	
+	}
 }
 
-/* f += g */
+void sfmintdeg(VL vl,P fx,int dy,int c,P *fr)
+{
+	V x,y;
+	int dx,dxdy,i,j,k,l,d,len,len0,u,dyk;
+	UP *rx;
+	DCP dc;
+	P t,f,mono,f1;
+	UP ut,h;
+	int ***nf;
+	int *r,*hist,*prev,*r1;
+	struct lb *lb;
+	GFS s;
+	NODE g;
 
-void addtosflum(bound,g,f)
-int bound;
-LUM g,f;
+	x = vl->v;
+	y = NEXT(vl)->v;
+	dx = getdeg(x,fx);
+	dxdy = dx*dy;
+	/* rx = -(fx-x^dx) */
+	rx = (UP *)CALLOC(dx,sizeof(UP));
+	for ( dc = DC(fx); dc; dc = NEXT(dc)) {
+		chsgnp(COEF(dc),&t);
+		ptoup(t,&ut);
+		rx[QTOS(DEG(dc))] = ut;
+	}
+	/* nf[d] = normal form table of monomials with total degree d */
+	nf = (int ***)CALLOC(dx+dy+1,sizeof(int **)); /* xxx */
+	nf[0] = (int **)CALLOC(1,sizeof(int *));
+
+	/* nf[0][0] = 1 */
+	r = (int *)CALLOC(dxdy,sizeof(int));
+	r[0] = _onesf();
+	nf[0][0] = r;
+
+	hist = (int *)CALLOC(1,sizeof(int));
+	hist[0] = _onesf();
+
+	lb = (struct lb *)CALLOC(1,sizeof(struct lb));
+	lb->pos = 0;
+	lb->r = r;
+	lb->hist = hist;
+	lb->len = 1;
+
+	/* g : table of normal form as linear form */
+	MKNODE(g,lb,0);
+
+	len = 1;
+	h = UPALLOC(dy);
+	for ( d = 1; ; d++ ) {
+		if ( d > c ){		
+			return;
+		}
+		nf[d] = (int **)CALLOC(d+1,sizeof(int *));
+		len0 = len;
+		len += d+1;
+
+		for ( i = d; i >= 0; i-- ) {
+			/* nf(x^(d-i)*y^i) = nf(y*nf(x^(d-i)*y^(i-1))) */
+			/* nf(x^d) = nf(nf(x^(d-1))*x) */
+			r = (int *)CALLOC(dxdy,sizeof(int));
+			if ( i == 0 ) {
+				prev = nf[d-1][0];
+				bcopy(prev,r+dy,(dxdy-dy)*sizeof(int));
+
+				/* create the head coeff */
+				for ( l = 0, k = dxdy-dy; l < dy; l++, k++ ) {
+					iftogfs(prev[k],&s);
+					COEF(h)[l] = (Num)s;
+				}
+				for ( l = dy-1; l >= 0 && !COEF(h)[l]; l--);
+				DEG(h) = l;
+
+				for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {
+					tmulup(rx[k],h,dy,&ut);
+					if ( ut ) 
+						for ( l = 0; l < dy; l++ ) {
+							s = (GFS)COEF(ut)[l];
+							if ( s ) {
+								u = CONT(s);
+								r[dyk+l] = _ADDSF(r[dyk+l],FTOIF(u));
+							}
+						}
+				}
+			} else {
+				prev = nf[d-1][i-1];
+				for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {
+					for ( l = 1; l < dy; l++ )
+						r[dyk+l] = prev[dyk+l-1];
+				}
+			}
+			nf[d][i] = r;
+			hist = (int *)CALLOC(len,sizeof(int));
+			hist[len0+i] = _onesf();
+			r1 = (int *)CALLOC(dxdy,sizeof(int));
+			bcopy(r,r1,dxdy*sizeof(int));
+			lnf(r1,hist,dxdy,len,g);
+			for ( k = 0; k < dxdy && !r1[k]; k++ );
+			if ( k == dxdy ) {
+				f = 0;
+				for ( k = j = 0; k <= d; k++ )
+					for ( i = 0; i <= k; i++, j++ )
+						if ( hist[j] ) {
+							iftogfs(hist[j],&s);
+							/* mono = s*x^(k-i)*y^i */
+							create_bmono((P)s,x,k-i,y,i,&mono);
+							addp(vl,f,mono,&f1);
+							f = f1;
+						}
+				*fr = f;
+				return;
+			}	else {
+				lb = (struct lb *)CALLOC(1,sizeof(struct lb));
+				lb->pos = k;
+				lb->r = r1;
+				lb->hist = hist;
+				lb->len = len;
+				g = insert_lb(g,lb);
+			}
+		}
+	}
+}
+
+static void create_bmono(P c,V x,int i,V y,int j,P *mono)
 {
-	int dg,i,j;
-	int **cg,**cf;
+	P t,s;
 
-	dg = DEG(g);
-	cg = COEF(g);
-	cf = COEF(f);
-	for ( i = 0; i <= dg; i++ )
-		for ( j = 0; j <= bound; j++ )
-			cf[i][j] = _ADDSF(cf[i][j],cg[i][j]);
+	if ( !i )
+		if ( !j )
+			t = c;
+		else {
+			/* c*y^j */
+			MKV(y,t);
+			COEF(DC(t)) = c;
+			STOQ(j,DEG(DC(t)));
+		}
+	else if ( !j ) {
+		/* c*x^i */
+		MKV(x,t);
+		COEF(DC(t)) = c;
+		STOQ(i,DEG(DC(t)));
+	} else {
+		MKV(y,s);
+		COEF(DC(s)) = c;
+		STOQ(j,DEG(DC(s)));
+		MKV(x,t);
+		COEF(DC(t)) = s;
+		STOQ(i,DEG(DC(t)));
+	}
+	*mono = t;
 }