===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/Mgfs.c,v
retrieving revision 1.11
retrieving revision 1.19
diff -u -p -r1.11 -r1.19
--- OpenXM_contrib2/asir2000/engine/Mgfs.c	2001/10/09 01:36:10	1.11
+++ OpenXM_contrib2/asir2000/engine/Mgfs.c	2015/08/14 13:51:54	1.19
@@ -1,13 +1,15 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.10 2001/09/03 01:04:26 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.18 2015/08/06 10:01:52 fujimoto 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;
 
 #if defined(__GNUC__)
-#define INLINE inline
-#elif defined(VISUAL)
+#define INLINE static inline
+#elif defined(VISUAL) || defined(__MINGW32__)
 #define INLINE __inline
 #else
 #define INLINE
@@ -21,6 +23,18 @@ INLINE int _ADDSF(int a,int 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];
@@ -48,9 +62,15 @@ INLINE int _ADDSF(int a,int 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;
@@ -200,8 +220,8 @@ void kmulsfum(UM n1,UM n2,UM nr)
 	len = (d+1)*d2;
 	r0 = (unsigned int *)ALLOCA(len*sizeof(int));
 	bzero((char *)r0,len*sizeof(int));
-	m = W_UMALLOC(d2+1);
-	carry = W_UMALLOC(d2+1);
+	m = W_UMALLOC(d1+d2+1);
+	carry = W_UMALLOC(d1+d2+1);
 	t = W_UMALLOC(d1+d2+1);
 	s = W_UMALLOC(d1+d2+1);
 	for ( DEG(carry) = -1, i = 0, r = r0; i < d; i++, r += d2 ) {
@@ -225,16 +245,17 @@ void kmulsfum(UM n1,UM n2,UM nr)
 
 void kmulsfummain(UM n1,UM n2,UM nr)
 {
-	int d1,d2,h,len;
+	int d1,d2,h,len,d;
 	UM n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
 
 	d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (d1+1)/2;
-	n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1);
-	n2lo = W_UMALLOC(d2+1); n2hi = W_UMALLOC(d2+1);
-	lo = W_UMALLOC(d1+d2+1); hi = W_UMALLOC(d1+d2+1);
-	mid1 = W_UMALLOC(d1+d2+1); mid2 = W_UMALLOC(d1+d2+1);
-	mid = W_UMALLOC(d1+d2+1);
-	s1 = W_UMALLOC(d1+d2+1); s2 = W_UMALLOC(d1+d2+1);
+	d = d1+d2+1;
+	n1lo = W_UMALLOC(d); n1hi = W_UMALLOC(d);
+	n2lo = W_UMALLOC(d); n2hi = W_UMALLOC(d);
+	lo = W_UMALLOC(d); hi = W_UMALLOC(d);
+	mid1 = W_UMALLOC(d); mid2 = W_UMALLOC(d);
+	mid = W_UMALLOC(d);
+	s1 = W_UMALLOC(d); s2 = W_UMALLOC(d);
 	extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi);
 	extractum(n2,0,h,n2lo); extractum(n2,h,d2-h,n2hi);
 	kmulsfum(n1hi,n2hi,hi); kmulsfum(n1lo,n2lo,lo);
@@ -300,7 +321,7 @@ void diffsfum(UM f,UM fd)
 
 	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);
 }
@@ -537,4 +558,221 @@ void clearbm(int n,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++ ) {
+					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)
+{
+	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;
 }