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