===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/Mgfs.c,v
retrieving revision 1.9
retrieving revision 1.13
diff -u -p -r1.9 -r1.13
--- OpenXM_contrib2/asir2000/engine/Mgfs.c	2001/07/23 05:05:41	1.9
+++ OpenXM_contrib2/asir2000/engine/Mgfs.c	2002/09/27 04:24:04	1.13
@@ -1,14 +1,11 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.8 2001/07/03 01:41:26 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.12 2002/01/13 07:11:46 noro Exp $ */
 
 #include "ca.h"
+#include "inline.h"
 
 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)
@@ -17,8 +14,7 @@ void kmulsfummain(UM,UM,UM);
 #define INLINE
 #endif
 
-INLINE int _ADDSF(a,b)
-int a,b;
+INLINE int _ADDSF(int a,int b)
 {
 	if ( !a )
 		return b;
@@ -26,6 +22,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];
@@ -51,12 +59,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;
@@ -64,8 +77,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;
 		
@@ -92,8 +104,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;
@@ -128,8 +139,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;
@@ -153,8 +163,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;
@@ -175,9 +185,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;
@@ -187,8 +195,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;
@@ -235,8 +242,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;
@@ -270,8 +276,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;
@@ -307,8 +312,7 @@ UM p1,p2,pq;
 	return i;
 }
 
-void diffsfum(f,fd)
-UM f,fd;
+void diffsfum(UM f,UM fd)
 {
 	int *dp,*sp;
 	int i;
@@ -320,8 +324,7 @@ UM f,fd;
 	degum(fd,DEG(f) - 1);
 }
 
-void monicsfum(f)
-UM f;
+void monicsfum(UM f)
 {
 	int *sp;
 	int i,inv;
@@ -332,16 +335,38 @@ UM f;
 		*sp = _MULSF(*sp,inv);
 }
 
+int compsfum(UM a,UM b)
+{
+	int i,da,db;
+
+	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;
+	}
+}
+
 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(f1,f2,fr)
-BM f1,f2,fr;
+void mulsfbm(BM f1,BM f2,BM fr)
 {
 	UM mul,t,s;
-	int i,j,h,d1,d2,dy;
+	int i,j,d1,d2,dy;
 
 	dy = MIN(DEG(f1),DEG(f2));
 
@@ -364,8 +389,7 @@ BM f1,f2,fr;
 	DEG(fr) = dy;
 }
 
-int degbm(f)
-BM f;
+int degbm(BM f)
 {
 	int d,i,dy;
 
@@ -378,8 +402,7 @@ BM f;
 
 /* g += f */
 
-void addtosfbm(f,g)
-BM f,g;
+void addtosfbm(BM f,BM g)
 {
 	int i,d1,d2,dy;
 	UM t;
@@ -397,8 +420,7 @@ BM f,g;
 	}
 }
 
-void eucsfum(f1,f2,a,b)
-UM f1,f2,a,b;
+void eucsfum(UM f1,UM f2,UM a,UM b)
 {
 	UM g1,g2,a1,a2,a3,wm,q,tum;
 	int d,dr;
@@ -432,10 +454,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;
@@ -459,10 +478,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;
@@ -499,9 +515,7 @@ UM g;
 
 /* f(y) -> f(y+a) */
 
-void shiftsfbm(f,a)
-BM f;
-int a;
+void shiftsfbm(BM f,int a)
 {
 	int i,j,d,dy;
 	UM pwr,ya,w,t,s;
@@ -534,9 +548,7 @@ int a;
 	}
 }
 
-void clearbm(n,f)
-int n;
-BM f;
+void clearbm(int n,BM f)
 {
 	int i,dy;
 	UM *c;
@@ -544,4 +556,226 @@ BM 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)
+{
+	NODE prev,cur,n;
+
+	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;
+}
+
+static void lnf(int *r,int *h,int n,int len,NODE g)
+{
+	struct lb *t;
+	int pos,i,j,len1,c;
+	int *r1,*h1;
+
+	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);	
+	}
+}
+
+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;
+
+	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++ ) {
+					if ( prev[k] ) {
+						u = IFTOF(prev[k]);
+						MKGFS(u,s);
+					} else
+						s = 0;
+					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] ) {
+							u = IFTOF(hist[j]);
+							MKGFS(u,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)
+{
+	P t,s;
+
+	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;
 }