Annotation of OpenXM/src/kan96xx/Kan/red.c, Revision 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>