Annotation of OpenXM/src/kan96xx/Kan/red.c, Revision 1.2
1.2 ! takayama 1: /* $OpenXM$ */
1.1 maekawa 2: #include <stdio.h>
3: #include "datatype.h"
4: #include "extern2.h"
5: #include "gradedset.h"
6:
7: #define mymax(p,q) (p>q?p:q)
8:
9: int DebugReductionRed = 0;
10: extern int Sugar;
11:
12: struct spValue sp_gen(f,g)
13: POLY f;
14: POLY g;
15: /* the results may be rewritten. */
16: {
17: struct spValue r;
18: POLY a;
19: POLY b;
20: MONOMIAL tf,tg;
21: MONOMIAL ta,tb;
22: int n,i;
23: struct coeff *ac;
24: struct coeff *bc;
25: int amodify,bmodify,c,ell;
26:
27: if ((f ISZERO) || (g ISZERO)) {
28: warningGradedSet("sp_gen(): f and g must not be zero. Returns zero.");
29: r.a = ZERO;
30: r.b = ZERO;
31: return(r);
32: }
33:
34: checkRingSp(f,g,r);
35: tf = f->m; tg = g->m;
36: n = tf->ringp->n;
37: ta = newMonomial(tf->ringp);
38: tb = newMonomial(tf->ringp);
39: for (i=0; i<n; i++) {
40: ta->e[i].x = mymax(tf->e[i].x,tg->e[i].x) - tf->e[i].x;
41: ta->e[i].D = mymax(tf->e[i].D,tg->e[i].D) - tf->e[i].D;
42: tb->e[i].x = mymax(tf->e[i].x,tg->e[i].x) - tg->e[i].x;
43: tb->e[i].D = mymax(tf->e[i].D,tg->e[i].D) - tg->e[i].D;
44: }
45:
46: switch (f->coeffp->tag) {
47: case INTEGER:
48: a = newCell(coeffCopy(g->coeffp),ta);
49: b = newCell(coeffCopy(f->coeffp),tb);
50: b->coeffp->val.i = -(b->coeffp->val.i);
51: break;
52: case MP_INTEGER:
53: {
54: MP_INT *gcd,*ac;
55: gcd = newMP_INT();
56: mpz_gcd(gcd,f->coeffp->val.bigp,g->coeffp->val.bigp);
57: ac = newMP_INT();
58: mpz_mdiv(ac,g->coeffp->val.bigp,gcd);
59: mpz_mdiv(gcd,f->coeffp->val.bigp,gcd);
60: mpz_neg(gcd,gcd);
61:
62: a = newCell(mpintToCoeff(ac,tf->ringp),ta);
63: b = newCell(mpintToCoeff(gcd,tf->ringp),tb);
64: }
65: break;
66: case POLY_COEFF:
67: c = f->m->ringp->c; ell = f->m->ringp->l;
68: if (ell-c > 0) { /* q-case */
69: amodify = 0;
70: for (i=c; i<ell; i++) {
71: amodify += (tb->e[i].D)*(tg->e[i].x);
72: }
73: bmodify = 0;
74: for (i=c; i<ell; i++) {
75: bmodify += (ta->e[i].D)*(tf->e[i].x);
76: }
77: if (amodify > bmodify) {
78: amodify = amodify-bmodify;
79: bmodify = 0;
80: }else{
81: bmodify = bmodify - amodify;
82: amodify = 0;
83: }
84: }
85: if (amodify) {
86: /* a ---> q^amodify a, b ---> - b */
87: ac = polyToCoeff(cxx(1,0,amodify,f->m->ringp->next),f->m->ringp);
88: Cmult(ac,ac,g->coeffp);
89: a = newCell(ac,ta);
90: bc = coeffNeg(f->coeffp,f->m->ringp);
91: b = newCell(bc,tb);
92: }else{
93: /* a ---> a, b ---> -q^bmodify b */
94: a = newCell(coeffCopy(g->coeffp),ta);
95: bc = coeffNeg(f->coeffp,f->m->ringp);
96: Cmult(bc,polyToCoeff(cxx(1,0,bmodify,f->m->ringp->next),f->m->ringp),bc);
97: b = newCell(bc,tb);
98: }
99: break;
100: default:
101: errorGradedSet("sp_gen(): Unsupported tag.");
102: break;
103: }
104: /* r.sp = ppAddv(ppMult(a,f),ppMult(b,g)); */
105: r.a = a;
106: r.b = b;
107: return(r);
108: }
109:
110:
111: POLY reduction1_gen_debug(f,g,needSyz,c,h)
112: POLY f;
113: POLY g;
114: int needSyz;
115: POLY *c; /* set */
116: POLY *h; /* set */
117: /* f must be reducible by g. r = c*f + h*g */
118: {
119: extern struct ring *CurrentRingp;
120: struct ring *rp;
121: struct spValue sv;
122: POLY f2;
123: int showLength = 0;
124:
125: /*DebugReductionRed = 1;*/
126: if (needSyz) {
127: if (f ISZERO) { rp = CurrentRingp; } else {rp = f->m->ringp; }
128: *c = cxx(1,0,0,rp);
129: *h = ZERO;
130: }
131:
132: sv = (*sp)(f,g);
133: f2 = ppAddv(cpMult((sv.a)->coeffp,f),ppMult(sv.b,g));
134: if (showLength) {printf(" [%d] ",pLength(f)); fflush(stdout);}
135: if (!isOrdered(f2) || !isOrdered(f)) {
136: if (!isOrdered(f)) {
137: printf("\nf is not ordered polynomial.");
138: }else if (!isOrdered(f)) {
139: printf("\nf2 is not ordered polynomial.");
140: }
141: printf("f=%s,\nf2=%s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
142: getchar();
143: getchar();
144: }
145:
146: if ((*mmLarger)(f2,f) >= 1) {
147: printf("error in reduction1.");
148: printf("f=%s --> f2=%s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
149: printf("f2 = (%s)*f+(%s)*(%s)\n",POLYToString(sv.a,'*',1),
150: POLYToString(sv.b,'*',1),POLYToString(g,'*',1));
151: getchar();
152: getchar();
153: }
154: if (DebugReductionRed) {
155: printf("%s --> %s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
156: }
157: f = f2;
158: if (needSyz) {
159: *c = ppMult(sv.a,*c);
160: *h = ppAdd(ppMult(sv.a,*h),sv.b);
161: }
162:
163: while ((*isReducible)(f,g)) {
164: sv = (*sp)(f,g);
165: f2 = ppAddv(cpMult((sv.a)->coeffp,f),ppMult(sv.b,g));
166: if (DebugReductionRed) {
167: printf("%s --> %s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
168: }
169: if (showLength) {printf(" [%d] ",pLength(f)); fflush(stdout);}
170: if (!isOrdered(f2) || !isOrdered(f)) {
171: if (!isOrdered(f)) {
172: printf("\nf is not ordered polynomial.");
173: }else if (!isOrdered(f)) {
174: printf("\nf2 is not ordered polynomial.");
175: }
176: printf("f=%s,\nf2=%s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
177: getchar();
178: getchar();
179: }
180:
181: if ((*mmLarger)(f2,f) >= 1) {
182: printf("error in reduction1.");
183: printf("f=%s --> f2=%s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
184: printf("f2 = (%s)*f+(%s)*(%s)\n",POLYToString(sv.a,'*',1),
185: POLYToString(sv.b,'*',1),
186: POLYToString(g,'*',1));
187: getchar();
188: getchar();
189: }
190: f = f2;
191: if (needSyz) {
192: *c = ppMult(sv.a,*c);
193: *h = ppAdd(ppMult(sv.a,*h),sv.b);
194: }
195: }
196: return(f);
197: }
198:
199: POLY reduction1_gen(f,g,needSyz,c,h)
200: POLY f;
201: POLY g;
202: int needSyz;
203: POLY *c; /* set */
204: POLY *h; /* set */
205: /* f must be reducible by g. r = c*f + h*g */
206: {
207: extern struct ring *CurrentRingp;
208: struct ring *rp;
209: struct spValue sv;
210: POLY f2;
211:
212:
213: if (needSyz) {
214: if (f ISZERO) { rp = CurrentRingp; } else {rp = f->m->ringp; }
215: *c = cxx(1,0,0,rp);
216: *h = ZERO;
217: }
218:
219: sv = (*sp)(f,g);
220: f2 = ppAddv(cpMult((sv.a)->coeffp,f),ppMult(sv.b,g));
221: if (DebugReductionRed) {
222: printf("c=%s, d=%s, g=%s: f --> c*f + d*g.\n",
223: POLYToString(sv.a,'*',1),POLYToString(sv.b,'*',1),POLYToString(g,'*',1));
224: printf("%s --> %s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
225: }
226: f = f2;
227: if (needSyz) {
228: *c = ppMult(sv.a,*c);
229: *h = ppAdd(ppMult(sv.a,*h),sv.b);
230: }
231:
232: while ((*isReducible)(f,g)) {
233: sv = (*sp)(f,g);
234: f2 = ppAddv(cpMult((sv.a)->coeffp,f),ppMult(sv.b,g));
235: if (DebugReductionRed) {
236: printf("! c=%s, d=%s, g=%s: f --> c*f + d*g.\n",
237: POLYToString(sv.a,'*',1),POLYToString(sv.b,'*',1),POLYToString(g,'*',1));
238: printf("%s --> %s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
239: }
240: f = f2;
241: if (needSyz) {
242: *c = ppMult(sv.a,*c);
243: *h = ppAdd(ppMult(sv.a,*h),sv.b);
244: }
245: }
246: return(f);
247: }
248:
249: POLY reduction1Cdr_gen(f,fs,g,needSyz,c,h)
250: POLY f;
251: POLY fs;
252: POLY g;
253: int needSyz;
254: POLY *c; /* set */
255: POLY *h; /* set */
256: /* f must be reducible by g. r = c*f + h*g */
257: {
258: extern struct ring *CurrentRingp;
259: struct ring *rp;
260: struct spValue sv;
261: POLY f2;
262:
263:
264: if (needSyz) {
265: if (f ISZERO) { rp = CurrentRingp; } else {rp = f->m->ringp; }
266: *c = cxx(1,0,0,rp);
267: *h = ZERO;
268: }
269:
270: sv = (*sp)(fs,g);
271: f2 = ppAddv(cpMult((sv.a)->coeffp,f),ppMult(sv.b,g));
272: if (DebugReductionRed) {
273: printf("%s --> %s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
274: }
275: f = f2;
276: if (needSyz) {
277: *c = ppMult(sv.a,*c);
278: *h = ppAdd(ppMult(sv.a,*h),sv.b);
279: }
280:
281:
282: while ((fs = (*isCdrReducible)(f,g)) != ZERO) {
283: sv = (*sp)(fs,g);
284: f2 = ppAddv(cpMult((sv.a)->coeffp,f),ppMult(sv.b,g));
285: if (DebugReductionRed) {
286: printf("%s --> %s\n",POLYToString(f,'*',1),POLYToString(f2,'*',1));
287: }
288: f = f2;
289: if (needSyz) {
290: *c = ppMult(sv.a,*c);
291: *h = ppAdd(ppMult(sv.a,*h),sv.b);
292: }
293: }
294: return(f);
295: }
296:
297:
298: /* for debug */
299: int isOrdered(f)
300: POLY f;
301: { POLY g;
302: if (f ISZERO) return(1);
303: g = f->next;
304: while (g != POLYNULL) {
305: if ((*mmLarger)(g,f)>=1) return(0);
306: f = g;
307: g = f->next;
308: }
309: return(1);
310: }
311:
312:
313: POLY reduction_gen(f,gset,needSyz,syzp)
314: POLY f;
315: struct gradedPolySet *gset;
316: int needSyz;
317: struct syz0 *syzp; /* set */
318: {
319: int reduced,reduced1,reduced2;
320: int grd;
321: struct polySet *set;
322: POLY cf,syz;
323: int i;
324: POLY cc,cg;
325:
326: extern struct ring *CurrentRingp;
327: struct ring *rp;
328:
329: if (needSyz) {
330: if (f ISZERO) { rp = CurrentRingp; } else { rp = f->m->ringp; }
331: cf = cxx(1,0,0,rp);
332: syz = ZERO;
333: }
334:
335: reduced = 0; /* no */
336: do {
337: reduced1 = 0; /* no */
338: grd = 0;
339: while (grd < gset->maxGrade) {
340: if (!Sugar) {
341: if (grd > (*grade)(f)) break;
342: }
343: set = gset->polys[grd];
344: do {
345: reduced2 = 0; /* no */
346: for (i=0; i<set->size; i++) {
347: if (f ISZERO) goto ss;
348: if ((*isReducible)(f,set->g[i])) {
349: f = (*reduction1)(f,set->g[i],needSyz,&cc,&cg);
350: if (needSyz) {
351: cf = ppMult(cc,cf);
352: syz = cpMult(toSyzCoeff(cc),syz);
353: syz = ppAddv(syz,toSyzPoly(cg,grd,i));
354: }
355: reduced = reduced1 = reduced2 = 1; /* yes */
356: }
357: }
358: } while (reduced2 != 0);
359: grd++;
360: }
361: }while (reduced1 != 0);
362:
363: ss: ;
364: if (needSyz) {
365: syzp->cf = cf; /* cf is in the CurrentRingp */
366: syzp->syz = syz; /* syz is in the SyzRingp */
367: }
368: return(f);
369: }
370:
371: POLY reduction_gen_rev(f,gset,needSyz,syzp)
372: POLY f;
373: struct gradedPolySet *gset;
374: int needSyz;
375: struct syz0 *syzp; /* set */
376: {
377: int reduced,reduced1,reduced2;
378: int grd;
379: struct polySet *set;
380: POLY cf,syz;
381: int i;
382: POLY cc,cg;
383:
384: extern struct ring *CurrentRingp;
385: struct ring *rp;
386:
387: if (needSyz) {
388: if (f ISZERO) { rp = CurrentRingp; } else { rp = f->m->ringp; }
389: cf = cxx(1,0,0,rp);
390: syz = ZERO;
391: }
392:
393: reduced = 0; /* no */
394: do {
395: reduced1 = 0; /* no */
396: grd = ((*grade)(f) < gset->maxGrade-1?(*grade)(f): gset->maxGrade-1);
397: while (grd >= 0) { /* reverse order for reduction */
398: set = gset->polys[grd];
399: do {
400: reduced2 = 0; /* no */
401: for (i=0; i<set->size; i++) {
402: if (f ISZERO) goto ss;
403: if ((*isReducible)(f,set->g[i])) {
404: f = (*reduction1)(f,set->g[i],needSyz,&cc,&cg);
405: if (needSyz) {
406: cf = ppMult(cc,cf);
407: syz = cpMult(toSyzCoeff(cc),syz);
408: syz = ppAddv(syz,toSyzPoly(cg,grd,i));
409: }
410: reduced = reduced1 = reduced2 = 1; /* yes */
411: }
412: }
413: } while (reduced2 != 0);
414: grd--;
415: }
416: }while (reduced1 != 0);
417:
418: ss: ;
419: if (needSyz) {
420: syzp->cf = cf; /* cf is in the CurrentRingp */
421: syzp->syz = syz; /* syz is in the SyzRingp */
422: }
423: return(f);
424: }
425:
426: POLY reductionCdr_gen(f,gset,needSyz,syzp)
427: POLY f;
428: struct gradedPolySet *gset;
429: int needSyz;
430: struct syz0 *syzp; /* set */
431: {
432: int reduced,reduced1,reduced2;
433: int grd;
434: struct polySet *set;
435: POLY cf,syz;
436: int i;
437: POLY cc,cg;
438: POLY fs;
439:
440: extern struct ring *CurrentRingp;
441: struct ring *rp;
442:
443: if (needSyz) {
444: if (f ISZERO) { rp = CurrentRingp; } else {rp = f->m->ringp; }
445: cf = cxx(1,0,0,rp);
446: syz = ZERO;
447: }
448:
449: reduced = 0; /* no */
450: do {
451: reduced1 = 0; /* no */
452: grd = 0;
453: while (grd < gset->maxGrade) {
454: if (!Sugar) {
455: if (grd > (*grade)(f)) break;
456: }
457: set = gset->polys[grd];
458: do {
459: reduced2 = 0; /* no */
460: for (i=0; i<set->size; i++) {
461: if (f ISZERO) goto ss;
462: if ((fs =(*isCdrReducible)(f,set->g[i])) != ZERO) {
463: f = (*reduction1Cdr)(f,fs,set->g[i],needSyz,&cc,&cg);
464: if (needSyz) {
465: cf = ppMult(cc,cf);
466: syz = cpMult(toSyzCoeff(cc),syz);
467: syz = ppAddv(syz,toSyzPoly(cg,grd,i));
468: }
469: reduced = reduced1 = reduced2 = 1; /* yes */
470: }
471: }
472: } while (reduced2 != 0);
473: grd++;
474: }
475: }while (reduced1 != 0);
476:
477: ss: ;
478: if (needSyz) {
479: syzp->cf = cf;
480: syzp->syz = syz;
481: }
482: return(f);
483: }
484:
485: int isReducible_gen(f,g)
486: POLY f;
487: POLY g;
488: {
489: int n,i;
490: MONOMIAL tf;
491: MONOMIAL tg;
492:
493: if (f ISZERO) return(0);
494: if (g ISZERO) return(0);
495:
496: checkRingIsR(f,g);
497:
498: tf = f->m; tg = g->m;
499: n = tf->ringp->n;
500: for (i=0; i<n; i++) {
501: if (tf->e[i].x < tg->e[i].x) return(0);
502: if (tf->e[i].D < tg->e[i].D) return(0);
503: }
504: return(1);
505: }
506:
507: POLY isCdrReducible_gen(f,g)
508: POLY f;
509: POLY g;
510: {
511: while (f != POLYNULL) {
512: if ((*isReducible)(f,g)) {
513: return(f);
514: }
515: f = f->next;
516: }
517: return(ZERO);
518: }
519:
520: POLY lcm_gen(f,g)
521: POLY f;
522: POLY g;
523: {
524: MONOMIAL tf,tg;
525: MONOMIAL lcm;
526: int n;
527: int i;
528:
529: tf = f->m; tg = g->m;
530: n = tf->ringp->n;
531: lcm = newMonomial(tf->ringp);
532: for (i=0; i<n; i++) {
533: lcm->e[i].x = mymax(tf->e[i].x,tg->e[i].x);
534: lcm->e[i].D = mymax(tf->e[i].D,tg->e[i].D);
535: }
536: return(newCell(intToCoeff(1,tf->ringp),lcm));
537: }
538:
539: int grade_gen(f)
540: POLY f;
541: {
542: int r;
543: int i,n;
544: MONOMIAL tf;
545: if (f ISZERO) return(-1);
546: tf = f->m;
547: n = tf->ringp->n;
548: r = 0;
549: for (i=0; i<n; i++) {
550: r += tf->e[i].x;
551: r += tf->e[i].D;
552: }
553: return(r);
554: }
555:
556: /* constructors */
557: POLY toSyzPoly(cg,grd,index)
558: POLY cg;
559: int grd;
560: int index;
561: /* the result is read only. */
562: {
563: extern struct ring *SyzRingp;
564: POLY r;
565: r = newCell(toSyzCoeff(cg),newMonomial(SyzRingp));
566: r->m->e[0].x = grd;
567: r->m->e[0].D = index;
568: return(r);
569: }
570:
571: struct coeff *toSyzCoeff(f)
572: POLY f;
573: {
574: extern struct ring *SyzRingp;
575: struct coeff *c;
576: c = newCoeff();
577: c->tag = POLY_COEFF;
578: c->val.f = f;
579: c->p = SyzRingp->p;
580: return(c);
581: }
582:
583: void initSyzRingp() {
584: extern struct ring *SyzRingp;
585: extern struct ring *CurrentRingp;
586: static char *x[]={"grade"};
587: static char *d[]={"index"};
588: static int order[]={1,0,
589: 0,1};
590: static int outputOrderForSyzRing[] = {0,1};
591: static int ringSerial = 0;
592: char *ringName = NULL;
593: SyzRingp = (struct ring *)sGC_malloc(sizeof(struct ring));
594: if (SyzRingp == (struct ring *)NULL)
595: errorGradedSet("initSyzRingp(); No more memory");
596: SyzRingp->p = CurrentRingp->p;
597: SyzRingp->n = 1; SyzRingp->m = 1; SyzRingp->l = 1; SyzRingp->c = 1;
598: SyzRingp->nn = 1; SyzRingp->mm = 1; SyzRingp->ll = 1;
599: SyzRingp->cc = 1;
600: SyzRingp->x = x;
601: SyzRingp->D = d;
602: SyzRingp->order = order;
603: SyzRingp->orderMatrixSize = 2;
604: setFromTo(SyzRingp);
605: SyzRingp->next = CurrentRingp;
606: SyzRingp->multiplication = mpMult_poly; /* Multiplication is not used.*/
607: SyzRingp->schreyer = 0;
608: SyzRingp->gbListTower = NULL;
609: SyzRingp->outputOrder = outputOrderForSyzRing;
610: ringName = (char *)sGC_malloc(20);
611: if (ringName == NULL) errorGradedSet("No more memory.");
612: sprintf(ringName,"syzring%05d",ringSerial);
613: SyzRingp->name = ringName;
614: }
615:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to red.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / kan96xx / Kan