Annotation of OpenXM/src/kan96xx/Kan/ecart.c, Revision 1.5
1.5 ! takayama 1: /* $OpenXM: OpenXM/src/kan96xx/Kan/ecart.c,v 1.4 2003/07/17 23:47:44 takayama Exp $ */
1.1 takayama 2: #include <stdio.h>
3: #include "datatype.h"
4: #include "extern2.h"
5: #include "gradedset.h"
6:
7: #define outofmemory(a) if (a == NULL) {fprintf(stderr,"ecart.c: out of memory.\n"); exit(10);}
8: /* Local structures */
9: struct ecartReducer {
10: int ell; /* s^\ell */
11: int first; /* first =1 ==> gset else gg */
12: int grade; int gseti; /* if first==1, gradedPolySet */
13: int ggi; /* if first==0 ecartPoly [] */
14: };
15: struct ecartPolyArray {
16: int size;
17: int limit;
18: POLY *pa;
1.5 ! takayama 19: POLY *cf;
! 20: POLY *syz;
1.1 takayama 21: };
22:
23: static struct ecartReducer ecartFindReducer(POLY r,struct gradedPolySet *gset,struct ecartPolyArray *epa);
24: static int ecartCheckPoly(POLY f); /* check if it does not contain s-variables */
25: static int ecartCheckEnv(); /* check if the environment is OK for ecart div*/
1.5 ! takayama 26: static struct ecartPolyArray *ecartPutPolyInG(POLY g,struct ecartPolyArray *eparray,POLY cf, POLY syz);
1.1 takayama 27: static int ecartGetEll(POLY r,POLY g);
1.5 ! takayama 28: static POLY ecartDivideSv(POLY r,int *d);
! 29: /* No automatic homogenization and s is used as a standart var. */
! 30: static POLY reduction_ecart0(POLY r,struct gradedPolySet *gset,
! 31: int needSyz, struct syz0 *syzp);
! 32: /* Automatic homogenization and s->1 */
! 33: static POLY reduction_ecart1(POLY r,struct gradedPolySet *gset,
! 34: int needSyz, struct syz0 *syzp);
! 35: static POLY ecartCheckSyz0(POLY cf,POLY r_0,POLY syz,
! 36: struct gradedPolySet *gg,POLY r);
1.1 takayama 37:
38: extern int DebugReductionRed;
39:
40: /* This is used for goHomogenization */
41: extern int DegreeShifto_size;
42: extern int *DegreeShifto_vec;
43:
1.3 takayama 44: /* It is used reduction_ecart() and ecartFindReducer()
45: to determine if we homogenize in this function */
46: extern int EcartAutomaticHomogenization;
47:
1.1 takayama 48: #define LARGE 0x7fffffff
49:
50:
1.5 ! takayama 51: static POLY ecartDivideSv(POLY r,int *d) {
1.1 takayama 52: POLY f;
53: int k;
1.5 ! takayama 54: *d = 0;
1.1 takayama 55: if (r == POLYNULL) return r;
56: f = r; k = LARGE;
57: while (r != POLYNULL) {
58: if (r->m->e[0].x < k) {
59: k = r->m->e[0].x;
60: }
61: r = r->next;
62: }
63: if (k > 0) {
1.5 ! takayama 64: *d = k;
1.1 takayama 65: f = r;
66: while (r != POLYNULL) {
67: r->m->e[0].x -= k;
68: r = r->next;
69: }
70: }
71: return f;
72: }
73:
74: static int ecartGetEll(POLY f,POLY g) {
75: int n,i,p;
76: MONOMIAL tf;
77: MONOMIAL tg;
78:
79: if (f ISZERO) return(-1);
80: if (g ISZERO) return(-1);
81:
82: checkRingIsR(f,g);
83:
84: if (!isSameComponent_x(f,g)) return(-1);
85: tf = f->m; tg = g->m; n = tf->ringp->n;
86: for (i=1; i<n; i++) {
87: if (tf->e[i].x < tg->e[i].x) return(-1);
88: if (tf->e[i].D < tg->e[i].D) return(-1);
89: }
90: if (tf->e[0].D < tg->e[0].D) return(-1); /* h */
1.2 takayama 91: p = tf->e[0].x - tg->e[0].x; /* H, s */
1.1 takayama 92: if (p >=0 ) return 0;
93: else return(-p);
94: }
95:
96:
97: #define EP_SIZE 10
1.5 ! takayama 98: static struct ecartPolyArray *ecartPutPolyInG(POLY g,struct ecartPolyArray *eparray,POLY cf,POLY syz)
1.1 takayama 99: {
100: struct ecartPolyArray *a;
101: POLY *ep;
102: int i;
103: if (eparray == (struct ecartPolyArray *)NULL) {
104: a = (struct ecartPolyArray *) sGC_malloc(sizeof(struct ecartPolyArray));
105: outofmemory(a);
106: ep = (POLY *) sGC_malloc(sizeof(POLY)*EP_SIZE);
107: outofmemory(ep);
1.5 ! takayama 108: a->cf = (POLY *) sGC_malloc(sizeof(POLY)*EP_SIZE);
! 109: outofmemory(a->cf);
! 110: a->syz = (POLY *) sGC_malloc(sizeof(POLY)*EP_SIZE);
! 111: outofmemory(a->syz);
1.1 takayama 112: a->limit = EP_SIZE;
113: a->size = 0;
114: a->pa = ep;
115: eparray = a;
116: }
117: if (eparray->size >= eparray->limit) {
118: a = (struct ecartPolyArray *) sGC_malloc(sizeof(struct ecartPolyArray));
119: outofmemory(a);
120: ep = (POLY *) sGC_malloc(sizeof(POLY)*((eparray->limit)*2));
121: outofmemory(ep);
1.5 ! takayama 122: a->cf = (POLY *) sGC_malloc(sizeof(POLY)*((eparray->limit)*2));
! 123: outofmemory(a->cf);
! 124: a->syz = (POLY *) sGC_malloc(sizeof(POLY)*((eparray->limit)*2));
! 125: outofmemory(a->syz);
1.1 takayama 126: a->limit = (eparray->limit)*2;
127: a->size = eparray->size;
128: a->pa = ep;
129: for (i=0; i<eparray->size; i++) {
130: (a->pa)[i] = (eparray->pa)[i];
1.5 ! takayama 131: (a->cf)[i] = (eparray->cf)[i];
! 132: (a->syz)[i] = (eparray->syz)[i];
1.1 takayama 133: }
134: eparray = a;
135: }
136: i = eparray->size;
137: (eparray->pa)[i] = g;
1.5 ! takayama 138: (eparray->cf)[i] = cf;
! 139: (eparray->syz)[i] = syz;
1.1 takayama 140: (eparray->size)++;
141: return eparray;
142: }
143:
144: static struct ecartReducer ecartFindReducer(POLY r,struct gradedPolySet *gset,
145: struct ecartPolyArray *epa)
146: {
147: int grd;
148: struct polySet *set;
149: int minGrade = 0;
150: int minGseti = 0;
151: int minGgi = 0;
152: int ell1 = LARGE;
153: int ell2 = LARGE;
154: int ell;
155: int i;
156: struct ecartReducer er;
157: /* Try to find a reducer in gset; */
158: grd = 0;
159: while (grd < gset->maxGrade) {
160: set = gset->polys[grd];
161: for (i=0; i<set->size; i++) {
1.2 takayama 162: if (set->gh[i] == POLYNULL) {
163: /* goHomogenize set->gh[i] */
1.4 takayama 164: if (EcartAutomaticHomogenization) {
165: set->gh[i] = goHomogenize11(set->g[i],DegreeShifto_vec,DegreeShifto_size,-1,1);
166: }else{
167: set->gh[i] = set->g[i];
168: }
1.2 takayama 169: }
170: ell = ecartGetEll(r,set->gh[i]);
1.1 takayama 171: if ((ell>=0) && (ell < ell1)) {
172: ell1 = ell;
173: minGrade = grd; minGseti=i;
174: }
175: }
176: grd++;
177: }
178: if (epa != NULL) {
179: /* Try to find in the second group. */
180: for (i=0; i< epa->size; i++) {
181: ell = ecartGetEll(r,(epa->pa)[i]);
182: if ((ell>=0) && (ell < ell2)) {
183: ell2 = ell;
184: minGgi = i;
185: }
186: }
187: }
188:
189: if (DebugReductionRed) {
190: printf("ecartFindReducer(): ell1=%d, ell2=%d, minGrade=%d, minGseti=%d, minGgi=%d\n",ell1,ell2,minGrade,minGseti,minGgi);
191: }
192: if (ell1 <= ell2) {
193: if (ell1 == LARGE) {
194: er.ell = -1;
195: return er;
196: }else{
197: er.ell = ell1;
198: er.first = 1;
199: er.grade = minGrade;
200: er.gseti = minGseti;
201: return er;
202: }
203: }else{
204: er.ell = ell2;
205: er.first = 0;
206: er.ggi = minGgi;
207: return er;
208: }
209: }
210:
1.5 ! takayama 211: static POLY ecartCheckSyz0(POLY cf,POLY r_0,POLY syz,
! 212: struct gradedPolySet *gg,POLY r)
! 213: {
! 214: POLY f;
! 215: int grd,i;
! 216: POLY q;
! 217: struct coeff *c;
! 218: f = ppMult(cf,r_0);
! 219: while (syz != POLYNULL) {
! 220: grd = syz->m->e[0].x;
! 221: i = syz->m->e[0].D;
! 222: c = syz->coeffp;
! 223: if (c->tag == POLY_COEFF) {
! 224: q = c->val.f;
! 225: }else{
! 226: q = POLYNULL;
! 227: }
! 228: f = ppAdd(f,ppMult(q,((gg->polys[grd])->g)[i]));
! 229: /* not gh, works only for _ecart0 */
! 230: syz = syz->next;
! 231: }
! 232: f = ppSub(f,r);
! 233: return f;
! 234: }
! 235:
! 236:
! 237: POLY reduction_ecart(r,gset,needSyz,syzp)
! 238: POLY r;
! 239: struct gradedPolySet *gset;
! 240: int needSyz;
! 241: struct syz0 *syzp; /* set */
! 242: {
! 243: if (EcartAutomaticHomogenization) {
! 244: return reduction_ecart1(r,gset,needSyz,syzp);
! 245: }else{
! 246: return reduction_ecart0(r,gset,needSyz,syzp);
! 247: }
! 248: }
! 249:
1.1 takayama 250: /*
251: r and gset are assumed to be (0,1)-homogenized (h-homogenized)
1.5 ! takayama 252: Polynomials r and gset are assumed
1.3 takayama 253: to be double homogenized (h-homogenized and s(H)-homogenized)
1.1 takayama 254: */
1.5 ! takayama 255: static POLY reduction_ecart0(r,gset,needSyz,syzp)
1.1 takayama 256: POLY r;
257: struct gradedPolySet *gset;
258: int needSyz;
259: struct syz0 *syzp; /* set */
260: {
261: int reduced,reduced1,reduced2;
262: int grd;
263: struct polySet *set;
264: POLY cf,syz;
265: int i;
266: POLY cc,cg;
267: struct ecartReducer ells;
268: struct ecartPolyArray *gg;
269: POLY pp;
270: int ell;
1.5 ! takayama 271: POLY cf_o;
! 272: POLY syz_o;
! 273: POLY r_0;
1.1 takayama 274:
275: extern struct ring *CurrentRingp;
276: struct ring *rp;
277:
1.5 ! takayama 278: r_0 = r;
1.1 takayama 279: gg = NULL;
280: if (needSyz) {
281: if (r ISZERO) { rp = CurrentRingp; } else { rp = r->m->ringp; }
282: cf = cxx(1,0,0,rp);
283: syz = ZERO;
284: }
285:
286: if (r != POLYNULL) {
1.4 takayama 287: rp = r->m->ringp;
288: if (! rp->weightedHomogenization) {
289: errorKan1("%s\n","ecart.c: the given ring must be declared with [(weightedHomogenization) 1]");
290: }
1.1 takayama 291: }
292:
1.5 ! takayama 293: do {
! 294: if (DebugReductionRed) printf("r=%s\n",POLYToString(r,'*',1));
! 295: ells = ecartFindReducer(r,gset,gg);
! 296: ell = ells.ell;
! 297: if (ell > 0) {
! 298: if (needSyz) {
! 299: gg = ecartPutPolyInG(r,gg,cf,syz);
! 300: }else{
! 301: gg = ecartPutPolyInG(r,gg,POLYNULL,POLYNULL);
! 302: }
! 303: }
! 304: if (ell >= 0) {
! 305: if (ells.first) {
! 306: pp = ((gset->polys[ells.grade])->gh)[ells.gseti];
! 307: }else{
! 308: pp = (gg->pa)[ells.ggi];
! 309: }
! 310: if (ell > 0) r = mpMult(cxx(1,0,ell,rp),r); /* r = s^ell r */
! 311: r = (*reduction1)(r,pp,needSyz,&cc,&cg);
! 312: if (needSyz) {
! 313: if (ells.first) {
! 314: if (ell >0) cc = ppMult(cc,cxx(1,0,ell,rp));
! 315: cf = ppMult(cc,cf);
! 316: syz = cpMult(toSyzCoeff(cc),syz);
! 317: syz = ppAdd(syz,toSyzPoly(cg,ells.grade,ells.gseti));
! 318: }else{
! 319: if (ell >0) cc = ppMult(cc,cxx(1,0,ell,rp));
! 320: cf_o = (gg->cf)[ells.ggi];
! 321: syz_o = (gg->syz)[ells.ggi];
! 322: cf = ppMult(cc,cf);
! 323: cf = ppAdd(cf,ppMult(cg,cf_o));
! 324: syz = cpMult(toSyzCoeff(cc),syz);
! 325: syz = ppAdd(syz,cpMult(toSyzCoeff(cg),syz_o));
! 326: /* Note. 2003.07.19 */
! 327: }
! 328: if (DebugReductionRed) {
! 329: POLY tp;
! 330: tp = ecartCheckSyz0(cf,r_0,syz,gset,r);
! 331: if (tp != POLYNULL) {
! 332: fprintf(stderr,"reduction_ecart0(): sygyzy is broken. Return the Current values.\n");
! 333: fprintf(stderr,"%s\n",POLYToString(tp,'*',1));
! 334: syzp->cf = cf;
! 335: syzp->syz = syz;
! 336: return r;
! 337: }
! 338: }
! 339: }
! 340: if (r ISZERO) goto ss;
! 341: /*r = ecartDivideSv(r,&se); r = r/s^? Don't do this. */
! 342: }
! 343: }while (ell >= 0);
! 344:
! 345: ss: ;
! 346: if (needSyz) {
! 347: syzp->cf = cf; /* cf is in the CurrentRingp */
! 348: syzp->syz = syz; /* syz is in the SyzRingp */
! 349: }
! 350:
! 351: return(r);
! 352: }
! 353:
! 354: /*
! 355: r and gset are assumed to be (0,1)-homogenized (h-homogenized)
! 356: r and gset are not assumed
! 357: to be double homogenized (h-homogenized and s(H)-homogenized)
! 358: They are automatically s(H)-homogenized, and s is set to 1
! 359: when this function returns.
! 360: */
! 361: static POLY reduction_ecart1(r,gset,needSyz,syzp)
! 362: POLY r;
! 363: struct gradedPolySet *gset;
! 364: int needSyz;
! 365: struct syz0 *syzp; /* set */
! 366: {
! 367: int reduced,reduced1,reduced2;
! 368: int grd;
! 369: struct polySet *set;
! 370: POLY cf,syz;
! 371: int i;
! 372: POLY cc,cg;
! 373: struct ecartReducer ells;
! 374: struct ecartPolyArray *gg;
! 375: POLY pp;
! 376: int ell;
! 377: int se;
! 378:
! 379: extern struct ring *CurrentRingp;
! 380: struct ring *rp;
! 381:
! 382: gg = NULL;
! 383: if (needSyz) {
! 384: if (r ISZERO) { rp = CurrentRingp; } else { rp = r->m->ringp; }
! 385: cf = cxx(1,0,0,rp);
! 386: syz = ZERO;
! 387: }
! 388:
! 389: if (r != POLYNULL) {
! 390: rp = r->m->ringp;
! 391: if (! rp->weightedHomogenization) {
! 392: errorKan1("%s\n","ecart.c: the given ring must be declared with [(weightedHomogenization) 1]");
! 393: }
1.3 takayama 394: }
1.5 ! takayama 395:
! 396: r = goHomogenize11(r,DegreeShifto_vec,DegreeShifto_size,-1,1);
1.1 takayama 397: /* 1 means homogenize only s */
398:
399: do {
400: if (DebugReductionRed) printf("r=%s\n",POLYToString(r,'*',1));
401: ells = ecartFindReducer(r,gset,gg);
402: ell = ells.ell;
403: if (ell > 0) {
1.5 ! takayama 404: gg = ecartPutPolyInG(r,gg,POLYNULL,POLYNULL);
1.1 takayama 405: }
406: if (ell >= 0) {
407: if (ells.first) {
408: pp = ((gset->polys[ells.grade])->gh)[ells.gseti];
409: }else{
410: pp = (gg->pa)[ells.ggi];
411: }
412: if (ell > 0) r = mpMult(cxx(1,0,ell,rp),r); /* r = s^ell r */
413: r = (*reduction1)(r,pp,needSyz,&cc,&cg);
414: if (needSyz) {
415: if (ells.first) {
416: cf = ppMult(cc,cf);
417: syz = cpMult(toSyzCoeff(cc),syz);
418: syz = ppAddv(syz,toSyzPoly(cg,ells.grade,ells.gseti));
419: }else{
1.5 ! takayama 420: fprintf(stderr,"It has not yet implemented.\n");
! 421: exit(10);
1.1 takayama 422: /* BUG: not yet */
423: }
424: }
1.5 ! takayama 425: if (r ISZERO) goto ss1;
! 426: r = ecartDivideSv(r,&se); /* r = r/s^? */
1.1 takayama 427: }
428: }while (ell >= 0);
429:
1.5 ! takayama 430: ss1: ;
1.1 takayama 431: if (needSyz) {
432: syzp->cf = cf; /* cf is in the CurrentRingp */
433: syzp->syz = syz; /* syz is in the SyzRingp */
434: /* BUG: dehomogenize the syzygy */
1.5 ! takayama 435: fprintf(stderr,"It has not yet implemented.\n");
! 436: exit(10);
1.1 takayama 437: }
1.5 ! takayama 438:
1.1 takayama 439: r = goDeHomogenizeS(r);
1.5 ! takayama 440:
1.1 takayama 441: return(r);
442: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>