=================================================================== 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 +#include #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; + } +}