===================================================================
RCS file: /home/cvs/OpenXM/src/kan96xx/Kan/poly2.c,v
retrieving revision 1.1.1.1
retrieving revision 1.7
diff -u -p -r1.1.1.1 -r1.7
--- OpenXM/src/kan96xx/Kan/poly2.c	1999/10/08 02:12:01	1.1.1.1
+++ OpenXM/src/kan96xx/Kan/poly2.c	2005/07/03 11:08:54	1.7
@@ -1,4 +1,6 @@
+/* $OpenXM: OpenXM/src/kan96xx/Kan/poly2.c,v 1.6 2005/06/16 05:07:23 takayama Exp $ */
 #include <stdio.h>
+#include <stdlib.h>
 #include "datatype.h"
 #include "stackm.h"
 #include "extern.h"
@@ -7,7 +9,7 @@
 static POLY mapZmonom(POLY f,struct ring *ringp);
 
 POLY ppAdd(f,g)
-POLY f; POLY g;  /* The result is read only. */
+	 POLY f; POLY g;  /* The result is read only. */
 {
   POLY node;
   struct listPoly nod;
@@ -32,8 +34,8 @@ POLY f; POLY g;  /* The result is read only. */
       h = h->next;
       f = f->next;
       if (f == POLYNULL) {
-	h->next = g;
-	return(node->next);
+		h->next = g;
+		return(node->next);
       }
       break;
     case 0: /* f < g */
@@ -41,37 +43,37 @@ POLY f; POLY g;  /* The result is read only. */
       h = h->next;
       g = g->next;
       if (g == POLYNULL) {
-	h->next = f;
-	return(node->next);
+		h->next = f;
+		return(node->next);
       }
       break;
     case 2:/* f == g */
       c = coeffCopy(f->coeffp);
       Cadd(c,c,g->coeffp);
       if (!isZero(c)) {
-	h->next = newCell(c,f->m);
-	h = h->next;
-	f = f->next;
-	g = g->next;
-	if (f == POLYNULL) {
-	  h->next = g;
-	  return(node->next);
-	}
-	if (g == POLYNULL) {
-	  h->next = f;
-	  return(node->next);
-	}
+		h->next = newCell(c,f->m);
+		h = h->next;
+		f = f->next;
+		g = g->next;
+		if (f == POLYNULL) {
+		  h->next = g;
+		  return(node->next);
+		}
+		if (g == POLYNULL) {
+		  h->next = f;
+		  return(node->next);
+		}
       }else{
-	f = f->next;
-	g = g->next;
-	if (f == POLYNULL) {
-	  h->next = g;
-	  return(node->next);
-	}
-	if (g == POLYNULL) {
-	  h->next = f;
-	  return(node->next);
-	}
+		f = f->next;
+		g = g->next;
+		if (f == POLYNULL) {
+		  h->next = g;
+		  return(node->next);
+		}
+		if (g == POLYNULL) {
+		  h->next = f;
+		  return(node->next);
+		}
       }
       break;
     default:
@@ -83,7 +85,7 @@ POLY f; POLY g;  /* The result is read only. */
 }
 
 POLY ppSub(f,g)
-POLY f; POLY g;  /* The result is read only. */
+	 POLY f; POLY g;  /* The result is read only. */
 {
   POLY h;
   struct coeff *c;
@@ -98,8 +100,8 @@ POLY f; POLY g;  /* The result is read only. */
 
 
 POLY cpMult(c,f)
-struct coeff *c;
-POLY f;
+	 struct coeff *c;
+	 POLY f;
 {
   POLY node;
   struct listPoly nod;
@@ -124,7 +126,7 @@ POLY f;
 }
 
 MONOMIAL monomialAdd_poly(m,m2)
-MONOMIAL m,m2;
+	 MONOMIAL m,m2;
 {
   extern int Msize;
   MONOMIAL f;
@@ -144,8 +146,8 @@ MONOMIAL m,m2;
 
 /* Note that mpMult_poly is called from mmLarger_tower! */
 POLY mpMult_poly(f,g)
-POLY f;
-POLY g;
+	 POLY f;
+	 POLY g;
 {
   POLY node;
   struct listPoly nod;
@@ -175,7 +177,7 @@ POLY g;
 }
 
 POLY ppMult_old(f,g)
-POLY f,g;
+	 POLY f,g;
 {
   POLY r;
   POLY tmp;
@@ -189,7 +191,7 @@ POLY f,g;
 }
 
 POLY ppAddv(f,g)
-POLY f; POLY g;  /* It breaks f and g. Use it just after calling mpMult() */
+	 POLY f; POLY g;  /* It breaks f and g. Use it just after calling mpMult() */
 {
   POLY node;
   struct listPoly nod;
@@ -214,44 +216,44 @@ POLY f; POLY g;  /* It breaks f and g. Use it just aft
       h->next = f;
       h = h->next; f = f->next;;
       if (f == POLYNULL) {
-	h->next = g;
-	return(node->next);
+		h->next = g;
+		return(node->next);
       }
       break;
     case 0: /* f < g */
       h->next =  g;
       h = h->next; g = g->next;
       if (g == POLYNULL) {
-	h->next = f;
-	return(node->next);
+		h->next = f;
+		return(node->next);
       }
       break;
     case 2:/* f == g */
       c = f->coeffp;
       Cadd(c,c,g->coeffp);
       if (!isZero(c)) {
-	h->next = f;
-	h = h->next; f = f->next;;
-	g = g->next;
-	if (f == POLYNULL) {
-	  h->next = g;
-	  return(node->next);
-	}
-	if (g == POLYNULL) {
-	  h->next = f;
-	  return(node->next);
-	}
+		h->next = f;
+		h = h->next; f = f->next;;
+		g = g->next;
+		if (f == POLYNULL) {
+		  h->next = g;
+		  return(node->next);
+		}
+		if (g == POLYNULL) {
+		  h->next = f;
+		  return(node->next);
+		}
       }else{
-	f = f->next;
-	g = g->next;
-	if (f == POLYNULL) {
-	  h->next = g;
-	  return(node->next);
-	}
-	if (g == POLYNULL) {
-	  h->next = f;
-	  return(node->next);
-	}
+		f = f->next;
+		g = g->next;
+		if (f == POLYNULL) {
+		  h->next = g;
+		  return(node->next);
+		}
+		if (g == POLYNULL) {
+		  h->next = f;
+		  return(node->next);
+		}
       }
       break;
     default:
@@ -263,8 +265,8 @@ POLY f; POLY g;  /* It breaks f and g. Use it just aft
 }
 
 POLY pPower(f,k)
-POLY f;
-int k;
+	 POLY f;
+	 int k;
 {
   POLY r;
   int i,n;
@@ -295,8 +297,8 @@ int k;
 }
 
 POLY pPower_poly(f,k)
-POLY f;
-int k;
+	 POLY f;
+	 int k;
 {
   POLY r;
   int i,n;
@@ -327,8 +329,8 @@ int k;
 }
 
 POLY modulop_trash(f,ringp)
-POLY f;
-struct ring *ringp;
+	 POLY f;
+	 struct ring *ringp;
 {
   int p;
   POLY h;
@@ -381,12 +383,12 @@ struct ring *ringp;
       mpz_mod_ui(c,f->coeffp->val.bigp,(unsigned long int)p);
       cc = (int) mpz_get_si(c);
       if (cc != 0) {
-	h->next = newCell(newCoeff(),monomialCopy(f->m));
-	h = h->next;
-	h->m->ringp = ringp;
-	h->coeffp->tag = INTEGER;
-	h->coeffp->p = p;
-	h->coeffp->val.i = cc;
+		h->next = newCell(newCoeff(),monomialCopy(f->m));
+		h = h->next;
+		h->m->ringp = ringp;
+		h->coeffp->tag = INTEGER;
+		h->coeffp->p = p;
+		h->coeffp->val.i = cc;
       }
       f = f->next;
     }
@@ -403,9 +405,9 @@ struct ring *ringp;
 }
   
 POLY modulop(f,ringp)
-POLY f;
-struct ring *ringp;
-/* Z[x] ---> R[x] where R=Z, Z/Zp, ringp->next. */
+	 POLY f;
+	 struct ring *ringp;
+	 /* Z[x] ---> R[x] where R=Z, Z/Zp, ringp->next. */
 {
   int p;
   POLY h;
@@ -446,12 +448,12 @@ struct ring *ringp;
       mpz_mod_ui(c,f->coeffp->val.bigp,(unsigned long int)p);
       cc = (int) mpz_get_si(c);
       if (cc != 0) {
-	h->next = newCell(newCoeff(),monomialCopy(f->m));
-	h = h->next;
-	h->m->ringp = ringp;
-	h->coeffp->tag = INTEGER;
-	h->coeffp->p = p;
-	h->coeffp->val.i = cc;
+		h->next = newCell(newCoeff(),monomialCopy(f->m));
+		h = h->next;
+		h->m->ringp = ringp;
+		h->coeffp->tag = INTEGER;
+		h->coeffp->p = p;
+		h->coeffp->val.i = cc;
       }
       f = f->next;
     }
@@ -470,9 +472,9 @@ struct ring *ringp;
 }
   
 POLY modulopZ(f,pcoeff)
-POLY f;
-struct coeff *pcoeff;
-/* Z[x] ---> Z[x] , f ---> f mod pcoeff*/
+	 POLY f;
+	 struct coeff *pcoeff;
+	 /* Z[x] ---> Z[x] , f ---> f mod pcoeff*/
 {
   int p;
   POLY h;
@@ -534,9 +536,9 @@ struct coeff *pcoeff;
 }
   
 struct pairOfPOLY quotientByNumber(f,pcoeff)
-POLY f;
-struct coeff *pcoeff;
-/* Z[x] ---> Z[x],Z[x] ,  f = first*pcoeff + second */
+	 POLY f;
+	 struct coeff *pcoeff;
+	 /* Z[x] ---> Z[x],Z[x] ,  f = first*pcoeff + second */
 {
   int p;
   POLY h;
@@ -630,8 +632,8 @@ struct coeff *pcoeff;
   
 
 POLY modulo0(f,ringp)
-POLY f;
-struct ring *ringp;
+	 POLY f;
+	 struct ring *ringp;
 {
   int p;
   POLY h;
@@ -649,8 +651,8 @@ struct ring *ringp;
       node = pcmCopy(f);
       f = node;
       while (f != POLYNULL) {
-	f->m->ringp = ringp; /* Touch the monomial "ringp" field. */
-	f = f->next;
+		f->m->ringp = ringp; /* Touch the monomial "ringp" field. */
+		f = f->next;
       }
       return(node);
     }
@@ -691,9 +693,9 @@ struct ring *ringp;
   
 
 struct object test(ob)  /* test3 */
-struct object ob;
+	 struct object ob;
 {
-  struct object rob;
+  struct object rob = OINIT;
   int k;
   static POLY f0;
   static POLY f1;
@@ -713,7 +715,7 @@ struct object ob;
     for (i=k; i>=0; i--) {
       f0->next = bxx(BiiPower(-k,i),0,i,CurrentRingp);
       if (f0->next != POLYNULL) {
-	f0 = f0->next;
+		f0 = f0->next;
       }
     }
     f0 = addNode->next;
@@ -725,7 +727,7 @@ struct object ob;
     for (i=k; i>=0; i--) {
       f1->next = bxx(BiiPower(k,i),0,i,CurrentRingp);
       if (f1->next != POLYNULL) {
-	f1 = f1->next;
+		f1 = f1->next;
       }
     }
     f1 = addNode->next;
@@ -745,7 +747,7 @@ struct object ob;
 
 
 int pLength(f)
-POLY f;
+	 POLY f;
 {
   int c=0;
   if (f ISZERO) return(0);
@@ -758,17 +760,17 @@ POLY f;
 
 
 POLY ppAddv2(f,g,top,nexttop)
-POLY f; POLY g;  /* It breaks f and g. Use it just after calling mpMult() */
-POLY top;
-POLY *nexttop;
-/* top is the starting address in the list f.
-   if top == POLYNULL, start from f.
+	 POLY f; POLY g;  /* It breaks f and g. Use it just after calling mpMult() */
+	 POLY top;
+	 POLY *nexttop;
+	 /* top is the starting address in the list f.
+		if top == POLYNULL, start from f.
 
-   *nexttop == 0
-            == g
-            == h or 0
+		*nexttop == 0
+		== g
+		== h or 0
 
-  It must be called as r = ppAddv2(r,g,...); 	    
+		It must be called as r = ppAddv2(r,g,...); 	    
 */
 {
   POLY node;
@@ -792,11 +794,11 @@ POLY *nexttop;
   if (top != POLYNULL) {
     while (f != top) {
       if (f == POLYNULL) {
-	fprintf(stderr,"\nppAddv2(): Internal error.\n");fflush(stderr);
-	fprintf(stderr,"f = %s\n",POLYToString(f0,'*',0));
-	fprintf(stderr,"g = %s\n",POLYToString(g0,'*',0));
-	fprintf(stderr,"top=%s\n",POLYToString(top,'*',0));
-	errorPoly("ppAddv2(). Internal error=1.");
+		fprintf(stderr,"\nppAddv2(): Internal error.\n");fflush(stderr);
+		fprintf(stderr,"f = %s\n",POLYToString(f0,'*',0));
+		fprintf(stderr,"g = %s\n",POLYToString(g0,'*',0));
+		fprintf(stderr,"top=%s\n",POLYToString(top,'*',0));
+		errorPoly("ppAddv2(). Internal error=1.");
       }
       h->next = f;
       h = h->next;
@@ -815,54 +817,54 @@ POLY *nexttop;
       h->next = f;
       h = h->next; f = f->next;;
       if (f == POLYNULL) {
-	h->next = g;
-	return(node->next);
+		h->next = g;
+		return(node->next);
       }
       break;
     case 0: /* f < g */
       h->next =  g;
       h = h->next; g = g->next;
       if (g == POLYNULL) {
-	h->next = f;
-	return(node->next);
+		h->next = f;
+		return(node->next);
       }
       break;
     case 2:/* f == g */
       c = g->coeffp;
       Cadd(c,f->coeffp,c);
       if (!isZero(c)) {
-	h->next = g;
-	h = h->next;
-	f = f->next;;
-	g = g->next;
-	if (f == POLYNULL) {
-	  h->next = g;
-	  return(node->next);
-	}
-	if (g == POLYNULL) {
-	  h->next = f;
-	  return(node->next);
-	}
+		h->next = g;
+		h = h->next;
+		f = f->next;;
+		g = g->next;
+		if (f == POLYNULL) {
+		  h->next = g;
+		  return(node->next);
+		}
+		if (g == POLYNULL) {
+		  h->next = f;
+		  return(node->next);
+		}
       }else{
-	if (g == g0) {
-	  if (h != node) {
-	    *nexttop = h;
-	  }else{
-	    *nexttop = POLYNULL;
-	  }
-	}
+		if (g == g0) {
+		  if (h != node) {
+			*nexttop = h;
+		  }else{
+			*nexttop = POLYNULL;
+		  }
+		}
 	
-	f = f->next;
-	g = g->next;
+		f = f->next;
+		g = g->next;
 
-	if (f == POLYNULL) {
-	  h->next = g;
-	  return(node->next);
-	}
-	if (g == POLYNULL) {
-	  h->next = f;
-	  return(node->next);
-	}
+		if (f == POLYNULL) {
+		  h->next = g;
+		  return(node->next);
+		}
+		if (g == POLYNULL) {
+		  h->next = f;
+		  return(node->next);
+		}
       }
       break;
     default:
@@ -874,7 +876,7 @@ POLY *nexttop;
 }
 
 POLY ppMult(f,g)
-POLY f,g;
+	 POLY f,g;
 {
   POLY r;
   POLY tmp;
@@ -893,7 +895,7 @@ POLY f,g;
 }
 
 POLY ppMult_poly(f,g)
-POLY f,g;
+	 POLY f,g;
 {
   POLY r;
   POLY tmp;
@@ -910,8 +912,8 @@ POLY f,g;
 }
 
 POLY mapZmonom(f,ringp)
-POLY f; /* assumes monomial. f \in Z[x] */
-struct ring *ringp;  /* R[x] */
+	 POLY f; /* assumes monomial. f \in Z[x] */
+	 struct ring *ringp;  /* R[x] */
 {
   struct ring *nextRing;
   struct ring nextRing0;
@@ -946,4 +948,82 @@ struct ring *ringp;  /* R[x] */
   node->coeffp->val.f = ff;
   return(node);
 }
+
+POLY reduceContentOfPoly(POLY f,struct coeff **contp) {
+  struct coeff *cont;
+  struct coeff *cOne = NULL;
+  extern struct ring *SmallRingp;
+  if (cOne == NULL) cOne = intToCoeff(1,SmallRingp);
+  *contp = cOne;
   
+  if (f == POLYNULL) return f;
+  if (f->m->ringp->p != 0) return f;
+  if (f->coeffp->tag != MP_INTEGER) return f;
+  cont = gcdOfCoeff(f);
+  *contp = cont;
+  if (coeffGreater(cont,cOne)) {
+	f = quotientByNumber(f,cont).first;
+  }
+  return f;
+}
+
+int coeffSizeMin(POLY f) {
+  int size;
+  int t;
+  if (f == POLYNULL) return 0;
+  if (f->m->ringp->p != 0) return 0;
+  if (f->coeffp->tag != MP_INTEGER) return 0;
+  size = mpz_size(f->coeffp->val.bigp);
+  while (f != POLYNULL) {
+	t = mpz_size(f->coeffp->val.bigp);
+	if (t < size)  size = t;
+	if (size == 1) return size;
+	f = f->next;
+  }
+}
+
+struct coeff *gcdOfCoeff(POLY f) {
+  extern struct ring *SmallRingp;
+  struct coeff *t;
+  MP_INT *tmp;
+  MP_INT *tmp2;
+  static MP_INT *cOne = NULL;
+  if (cOne == NULL) {
+	cOne = newMP_INT();
+	mpz_set_si(cOne,(long) 1);
+  }
+  if (f == POLYNULL) return intToCoeff(0,SmallRingp);
+  if (f->m->ringp->p != 0) return intToCoeff(0,SmallRingp);
+  if (f->coeffp->tag != MP_INTEGER) return intToCoeff(0,SmallRingp);
+  tmp = f->coeffp->val.bigp;
+  tmp2 = newMP_INT();
+  while (f != POLYNULL) {
+    mpz_gcd(tmp2,tmp,f->coeffp->val.bigp); /* tmp = tmp2 OK? */
+	tmp = tmp2;
+	if (mpz_cmp(tmp,cOne)==0) return intToCoeff(1,SmallRingp);
+	f = f->next;
+  }
+  return mpintToCoeff(tmp,SmallRingp);
+  
+}
+
+int shouldReduceContent(POLY f,int ss) {
+  extern DoCancel;
+  static int prevSize = 1;
+  int size;
+  if (f == POLYNULL) return 0;
+  if (f->m->ringp->p != 0) return 0;
+  if (f->coeffp->tag != MP_INTEGER) return 0;
+  if (DoCancel & 2) return 1;
+  /* Apply the Noro strategy to reduce content. */
+  size = mpz_size(f->coeffp->val.bigp);
+  if (ss > 0) {
+	prevSize = size;
+	return 0;
+  }
+  if (size > 2*prevSize) {
+	return 1;
+  }else{
+	return 0;
+  }
+}