Annotation of OpenXM/src/kan96xx/Kan/syz0.c, Revision 1.1
1.1 ! maekawa 1: #include <stdio.h>
! 2: #include "datatype.h"
! 3: #include "extern2.h"
! 4: #include "matrix.h"
! 5: #include "gradedset.h"
! 6:
! 7: extern int KanGBmessage;
! 8: static int Debugsyz0 = 0;
! 9:
! 10: static POLY getFactor0(struct gradedPolySet *grG,int grd,int index);
! 11: static void clearMark(struct gradedPolySet *grG);
! 12: static void printMatrixOfPOLY(struct matrixOfPOLY *mat);
! 13: static void printArrayOfPOLY(struct arrayOfPOLY *mat);
! 14: static int isZeroRow(struct matrixOfPOLY *mat,int i);
! 15:
! 16: static int getSize00OfGradedSet(struct gradedPolySet *g);
! 17: static int positionInPairs(struct pair *top,struct pair *pairs);
! 18: /* if the top is not in the list pairs, it returns -1. */
! 19: static struct pair *oldPairsToNewPairs(struct pair *opairs,
! 20: int *table,int size);
! 21: /* In the function getSyzygy(), reduced Grobner basis *grBases is
! 22: constructed from grG by deleting unnecessary elements.
! 23: The array <<table>> gives the correspondence between the index
! 24: ig (i-grade), ii (i-index) of *grBases and grG.
! 25: When <<opairs>> keeps syzygies for grG, it returns syzygies for
! 26: *grBases.
! 27: Note that the result syzygies are not necessarily COMPLETE.
! 28: Note that the <<prev>> fields are not set properly.
! 29: */
! 30: struct matrixOfPOLY *getSyzygy01(struct gradedPolySet *reducedBasis,
! 31: struct pair *excludedPairs);
! 32: /*
! 33: When <<excludedPairs>> is NULL,
! 34: this function computes all syzygies for the <<reducedBasis>>.
! 35: If not, this function computes all syzygies except those in excludedPairs.
! 36: The <<serial>> field of the <<reducedBasis>> is used to change
! 37: syzygy polynomials into arrays. March 15, 1997
! 38: */
! 39:
! 40:
! 41:
! 42: /* if (mark[j]), then the simplification is done. */
! 43: void getBackwardTransformation(grG)
! 44: struct gradedPolySet *grG;
! 45: /* grG->polys[i]->mark[j],syz[j] are modified. */
! 46: {
! 47: int i,j;
! 48: struct polySet *ps;
! 49:
! 50:
! 51: for (i=0; i<grG->maxGrade; i++) {
! 52: ps = grG->polys[i];
! 53: for (j=0; j<ps->size; j++) {
! 54: if (ps->mark[j] == 0 && ps->del[j] == 0) {
! 55: simplifyBT(i,j,grG);
! 56: }
! 57: }
! 58: }
! 59: }
! 60:
! 61: void simplifyBT(grd,index,grG)
! 62: int grd;
! 63: int index;
! 64: struct gradedPolySet *grG;
! 65: {
! 66: POLY s,r;
! 67: int g0,i0;
! 68:
! 69: s = grG->polys[grd]->syz[index]->syz;
! 70: r = ZERO;
! 71:
! 72: if (KanGBmessage) {printf("."); fflush(stdout); }
! 73:
! 74: while (s != ZERO) {
! 75: g0 = srGrade(s);
! 76: i0 = srIndex(s);
! 77: if (grG->polys[g0]->mark[i0]) {
! 78: r = ppAdd(r,cpMult(s->coeffp,grG->polys[g0]->syz[i0]->syz));
! 79: }else{
! 80: simplifyBT(g0,i0,grG);
! 81: r = ppAdd(r,cpMult(s->coeffp,grG->polys[g0]->syz[i0]->syz));
! 82: }
! 83: s = s->next;
! 84: }
! 85: grG->polys[grd]->syz[index]->syz = r;
! 86: grG->polys[grd]->mark[index] = 1;
! 87: }
! 88:
! 89:
! 90: static POLY getFactor0(grG,grd,index)
! 91: struct gradedPolySet *grG;
! 92: int grd;
! 93: int index;
! 94: {
! 95: POLY ans;
! 96: int i,j;
! 97: struct polySet *ps;
! 98:
! 99: ans = cxx(1,0,0,grG->polys[grd]->g[index]->m->ringp);
! 100: for (i=0; i<grG->maxGrade; i++) {
! 101: ps = grG->polys[i];
! 102: for (j=0; j<ps->size; j++) {
! 103: if (ps->mark[j] && !(i == grd && j == index)) {
! 104: ans = ppMult(ps->syz[j]->cf,ans);
! 105: }
! 106: }
! 107: }
! 108: return(ans);
! 109: }
! 110:
! 111: struct arrayOfPOLY *getSyzygy0(grG,zeroPairs)
! 112: struct gradedPolySet *grG;
! 113: struct pair *zeroPairs;
! 114: {
! 115: struct pair *tmp;
! 116: int size,k,i;
! 117: POLY s;
! 118: struct arrayOfPOLY *r;
! 119: struct arrayOfPOLY *ans;
! 120:
! 121: /* outputGradedPolySet(grG); */
! 122: /* outputNode(zeroPairs); */
! 123: /* outputNode() workds as outputPairs() */
! 124: tmp = zeroPairs; size = 0;
! 125: while (tmp != (struct pair *)NULL) {
! 126: size++;
! 127: tmp = tmp->next;
! 128: }
! 129: ans = newArrayOfPOLY(size+1);
! 130:
! 131: k = 0;
! 132:
! 133: while (zeroPairs != (struct pair *)NULL) {
! 134: s = getSyzPolyFromSp(zeroPairs,grG);
! 135: if (s ISZERO) {
! 136:
! 137: }else {
! 138: getArrayOfPOLY(ans,k) = s;
! 139: k++;
! 140: }
! 141: zeroPairs = zeroPairs->next;
! 142: }
! 143:
! 144: r = newArrayOfPOLY(k);
! 145: for (i=0; i<k; i++) {
! 146: getArrayOfPOLY(r,i) = getArrayOfPOLY(ans,i);
! 147: }
! 148: return(r);
! 149: }
! 150:
! 151: struct matrixOfPOLY *getSyzygy(grG,zeroPairs,grBasesp,backwardMatp)
! 152: struct gradedPolySet *grG;
! 153: struct pair *zeroPairs;
! 154: struct gradedPolySet **grBasesp;
! 155: struct matrixOfPOLY **backwardMatp;
! 156: {
! 157: int serial;
! 158: int i,j,kk;
! 159: struct polySet *ps;
! 160: POLY fi;
! 161: int grd,indx;
! 162: struct gradedPolySet *newG;
! 163: POLY r;
! 164: struct syz0 syz;
! 165: struct arrayOfPOLY *ap,*vec;
! 166: struct matrixOfPOLY *mp;
! 167: struct matrixOfPOLY *b,*nc,*m2,*m1,*ans,*ans2;
! 168: struct arrayOfPOLY *dc;
! 169: int size;
! 170: int *table;
! 171: struct pair *excludePairs;
! 172: struct matrixOfPOLY *mp0;
! 173: struct matrixOfPOLY *m0;
! 174: extern int Sugar;
! 175:
! 176: size = getSize00OfGradedSet(grG);
! 177: table = (int *)sGC_malloc(sizeof(int)*4*(size+1));
! 178: if (table == (int *)NULL) errorSyz0("Out of memory.");
! 179:
! 180: if (Debugsyz0) { printf("grG=");outputGradedPolySet(grG,1);}
! 181:
! 182: /* set grBases */
! 183: *grBasesp = newGradedPolySet(0);
! 184: serial = 0;
! 185: for (i=0; i<grG->maxGrade; i++) {
! 186: ps = grG->polys[i];
! 187: for (j=0; j<ps->size; j++) {
! 188: if (ps->del[j] == 0) {
! 189: fi = ps->g[j];
! 190: grd = -1;whereInG(*grBasesp,fi,&grd,&indx,Sugar);
! 191: *grBasesp =
! 192: putPolyInG(*grBasesp,fi,grd,indx,newSyz0(),0,serial);
! 193: table[4*serial] = i; table[4*serial+1] = j;
! 194: table[4*serial+2] = grd; table[4*serial+3] = indx;
! 195: serial++;
! 196: }
! 197: }
! 198: }
! 199: if (Debugsyz0) {printf("*grBasesp=");outputGradedPolySet(*grBasesp,1);}
! 200:
! 201: if (size != serial) {
! 202: fprintf(stderr," size(%d) != serial(%d) ",size,serial);
! 203: errorSyz0("Internal error size != serial");
! 204: }
! 205:
! 206: /* set newG. Compute the forward transformations.*/
! 207: newG = gradedPolySetCopy(grG);
! 208: for (i=0; i<grG->maxGrade;i++) {
! 209: ps = newG->polys[i];
! 210: for (j=0; j<ps->size; j++) {
! 211: fi = ps->g[j];
! 212: r = (*reduction)(fi,*grBasesp,1,&syz);
! 213: if (KanGBmessage) {printf("."); fflush(stdout);}
! 214: if (! (r ISZERO)) errorSyz0("Internal error; getSyzygy(): groebner basis must be given.");
! 215: ps->syz[j] = newSyz0();
! 216: /* printf("%s : ",POLYToString(syz.cf,'*',1));
! 217: printf("%s\n",POLYToString(syz.syz,'*',1)); */
! 218: ps->syz[j]->cf = syz.cf;
! 219: ps->syz[j]->syz = syz.syz;
! 220: }
! 221: }
! 222:
! 223: ap = getSyzygy0(newG,zeroPairs);
! 224:
! 225: /* */
! 226: excludePairs = oldPairsToNewPairs(zeroPairs,table,size);
! 227: if (Debugsyz0) {
! 228: printf("zeroPairs = "); outputNode(zeroPairs);
! 229: for (i=0; i<size; i++) {
! 230: printf("old gi=%d, old ii=%d, new gi=%d, new ii=%d\n",
! 231: table[4*i],table[4*i+1],table[4*i+2],table[4*i+3]);
! 232: }
! 233: printf("excludePairs = "); outputNode(excludePairs);
! 234: printf("\n");
! 235: }
! 236:
! 237: if (KanGBmessage) {printf("#"); fflush(stdout); }
! 238: mp0 = getSyzygy01(*grBasesp,excludePairs);
! 239: if (KanGBmessage) {printf("#"); fflush(stdout); }
! 240:
! 241: /* We compute E-CB. The number of G (G-basis) is serial. */
! 242: /* F = - C G, G = B F , then F + CB F = 0*/
! 243: b = getBackwardMatrixOfPOLY(grG);
! 244: nc = getNC(newG,serial,*grBasesp);
! 245: dc = getDC(newG);
! 246: /*{
! 247: printf("\nb=\n"); printMatrixOfPOLY(b);
! 248: printf("\n\nnc=\n"); printMatrixOfPOLY(nc);
! 249: printf("\n\ndc=\n"); printArrayOfPOLY(dc);
! 250: }*/
! 251:
! 252: if (KanGBmessage) {printf(":"); fflush(stdout);}
! 253: m2 = getSyzygy1(b,nc,dc);
! 254: *backwardMatp = b;
! 255:
! 256: /* mp is the syzygy of the G-basis. */
! 257: mp = newMatrixOfPOLY(ap->n,serial);
! 258: for (i=0; i < ap->n; i++) {
! 259: vec = syzPolyToArrayOfPOLY(serial,getArrayOfPOLY(ap,i),*grBasesp);
! 260: for (j=0; j < serial; j++) {
! 261: getMatrixOfPOLY(mp,i,j) = getArrayOfPOLY(vec,j);
! 262: }
! 263: }
! 264: /*printf(" ap->n = %d, serial=%d \n",ap->n, serial);
! 265: printf("\nmp = \n"); printMatrixOfPOLY(mp); */
! 266:
! 267: if (KanGBmessage) {printf(";"); fflush(stdout);}
! 268: m1 = aaMult(mp,b);
! 269: m0 = aaMult(mp0,b);
! 270: /* printf("\nm0 = \n"); printMatrixOfPOLY(m0); */
! 271:
! 272: ans = newMatrixOfPOLY((m0->m)+(m1->m)+(m2->m),m2->n);
! 273: kk = 0;
! 274: for (i=0; i<m0->m; i++) {
! 275: if (!isZeroRow(m0,i)) {
! 276: for (j=0; j<m0->n; j++) {
! 277: getMatrixOfPOLY(ans,kk,j) = getMatrixOfPOLY(m0,i,j);
! 278: }
! 279: kk++;
! 280: }
! 281: }
! 282: for (i=0; i<m1->m; i++) {
! 283: if (!isZeroRow(m1,i)) {
! 284: for (j=0; j<m1->n; j++) {
! 285: getMatrixOfPOLY(ans,kk,j) = getMatrixOfPOLY(m1,i,j);
! 286: }
! 287: kk++;
! 288: }
! 289: }
! 290: for (i=0; i<m2->m; i++) {
! 291: if (!isZeroRow(m2,i)) {
! 292: for (j=0; j<m2->n; j++) {
! 293: getMatrixOfPOLY(ans,kk,j) = getMatrixOfPOLY(m2,i,j);
! 294: }
! 295: kk++; if (KanGBmessage) printf("*");
! 296: }
! 297: }
! 298: if (kk != ans->m) {
! 299: ans2 = newMatrixOfPOLY(kk,ans->n);
! 300: for (i=0; i<kk; i++) {
! 301: for (j=0; j<ans->n; j++) {
! 302: getMatrixOfPOLY(ans2,i,j) = getMatrixOfPOLY(ans,i,j);
! 303: }
! 304: }
! 305: return(ans2);
! 306: }else{
! 307: return(ans);
! 308: }
! 309:
! 310: }
! 311:
! 312: POLY getSyzPolyFromSp(spij,grG)
! 313: struct pair *spij;
! 314: struct gradedPolySet *grG;
! 315: {
! 316: int ig,ii,jg,ji;
! 317: POLY dk;
! 318: POLY f;
! 319: POLY t,ans;
! 320: int grd,indx;
! 321:
! 322: ig = spij->ig; ii = spij->ii;
! 323: jg = spij->jg; ji = spij->ji;
! 324: if (grG->polys[ig]->del[ii] != 0 || grG->polys[jg]->del[ji] != 0)
! 325: return(ZERO);
! 326: dk = spij->syz; /* in SyzRing */
! 327: clearMark(grG);
! 328: f = dk;
! 329: while (f != POLYNULL) {
! 330: grd = srGrade(f);
! 331: indx = srIndex(f);
! 332: grG->polys[grd]->mark[indx] = 1;
! 333: f=f->next;
! 334: }
! 335:
! 336: f = dk; ans = ZERO;
! 337: while (f != POLYNULL) {
! 338: grd = srGrade(f);
! 339: indx = srIndex(f);
! 340: t = grG->polys[grd]->syz[indx]->syz;
! 341: t = cpMult(f->coeffp,t);
! 342: t = cpMult(toSyzCoeff(getFactor0(grG,grd,indx)),t);
! 343: ans = ppAdd(ans,t);
! 344: f = f->next;
! 345: }
! 346: return(ans);
! 347: }
! 348:
! 349: static void clearMark(grG)
! 350: struct gradedPolySet *grG;
! 351: {
! 352: int i,j;
! 353: struct polySet *ps;
! 354: for (i=0; i<grG->maxGrade; i++) {
! 355: ps = grG->polys[i];
! 356: for (j=0; j<ps->size; j++) {
! 357: ps->mark[j] = 0;
! 358: }
! 359: }
! 360: }
! 361:
! 362:
! 363: struct arrayOfPOLY *syzPolyToArrayOfPOLY(size,f,grG)
! 364: int size;
! 365: POLY f; /* f is in the SyzRingp */
! 366: struct gradedPolySet *grG;
! 367: {
! 368: struct arrayOfPOLY *ap;
! 369: int i,g0,i0,serial;
! 370:
! 371: ap = newArrayOfPOLY(size);
! 372: for (i=0; i<size; i++) {
! 373: getArrayOfPOLY(ap,i) = ZERO;
! 374: }
! 375:
! 376: while (f != POLYNULL) {
! 377: g0 = srGrade(f);
! 378: i0 = srIndex(f);
! 379: serial = grG->polys[g0]->serial[i0];
! 380: if (serial < 0) {
! 381: errorSyz0("syzPolyToArrayOfPOLY(): invalid serial[i] of grG.");
! 382: }
! 383: if (getArrayOfPOLY(ap,serial) != ZERO) {
! 384: errorSyz0("syzPolyToArrayOfPOLY(): syzygy polynomial is broken.");
! 385: }
! 386: getArrayOfPOLY(ap,serial) = srSyzCoeffToPOLY(f->coeffp);
! 387: f = f->next;
! 388: }
! 389: return(ap);
! 390: }
! 391:
! 392:
! 393:
! 394: struct matrixOfPOLY *getBackwardMatrixOfPOLY(struct gradedPolySet *grG)
! 395: {
! 396: /* use serial, del. cf. getBackwardArray() */
! 397: int inputSize,outputSize;
! 398: int i,j,k,p;
! 399: struct arrayOfPOLY *vec;
! 400: struct matrixOfPOLY *mat;
! 401: struct polySet *ps;
! 402:
! 403: inputSize = 0; outputSize = 0;
! 404: for (i=0; i<grG->maxGrade; i++) {
! 405: ps = grG->polys[i];
! 406: for (j=0; j<ps->size; j++) {
! 407: if (ps->serial[j] >= 0) ++inputSize;
! 408: if (ps->del[j] == 0) ++outputSize;
! 409: }
! 410: }
! 411:
! 412: mat = newMatrixOfPOLY(outputSize,inputSize);
! 413: k = 0;
! 414: for (i=0; i<grG->maxGrade; i++) {
! 415: ps = grG->polys[i];
! 416: for (j=0; j<ps->size; j++) {
! 417: if (ps->del[j] == 0) {
! 418: vec = syzPolyToArrayOfPOLY(inputSize,ps->syz[j]->syz,grG);
! 419: for (p=0; p<inputSize; p++) {
! 420: getMatrixOfPOLY(mat,k,p)=getArrayOfPOLY(vec,p);
! 421: }
! 422: k++;
! 423: }
! 424: }
! 425: }
! 426: return(mat);
! 427: }
! 428:
! 429:
! 430: struct matrixOfPOLY *getNC(newG,n,grBases)
! 431: struct gradedPolySet *newG; /* F is stored and indexed by serial. */
! 432: int n; /* The number of G. */
! 433: struct gradedPolySet *grBases; /* G (G-basis) is stored. */
! 434: {
! 435: int size,i,j,k,ii;
! 436: struct matrixOfPOLY *mat;
! 437: struct polySet *ps;
! 438: struct arrayOfPOLY *vec;
! 439:
! 440: if (newG == (struct gradedPolySet *)NULL) {
! 441: return(newMatrixOfPOLY(0,0));
! 442: }
! 443: size = 0;
! 444: for (i=0; i<newG->maxGrade; i++) {
! 445: ps = newG->polys[i];
! 446: for (j=0; j<ps->size; j++) {
! 447: if (ps->serial[j] >= 0) size++;
! 448: }
! 449: }
! 450: mat = newMatrixOfPOLY(size,n);
! 451: for (i=0; i<newG->maxGrade; i++) {
! 452: ps = newG->polys[i];
! 453: for (j=0; j<ps->size; j++) {
! 454: if (ps->serial[j] >= 0) {
! 455: ii = ps->serial[j];
! 456: vec = syzPolyToArrayOfPOLY(n,ps->syz[j]->syz,grBases);
! 457: for (k=0; k<n; k++) {
! 458: getMatrixOfPOLY(mat,ii,k) = getArrayOfPOLY(vec,k);
! 459: }
! 460: }
! 461: }
! 462: }
! 463: return(mat);
! 464: }
! 465:
! 466: struct arrayOfPOLY *getDC(newG)
! 467: struct gradedPolySet *newG; /* F is stored and indexed by serial. */
! 468: {
! 469: int size,i,j,k,ii;
! 470: struct arrayOfPOLY *mat;
! 471: struct polySet *ps;
! 472: extern struct ring *CurrentRingp;
! 473:
! 474: if (newG == (struct gradedPolySet *)NULL) {
! 475: return(newArrayOfPOLY(0));
! 476: }
! 477: size = 0;
! 478: for (i=0; i<newG->maxGrade; i++) {
! 479: ps = newG->polys[i];
! 480: for (j=0; j<ps->size; j++) {
! 481: if (ps->serial[j] >= 0) size++;
! 482: }
! 483: }
! 484: mat = newArrayOfPOLY(size);
! 485: for (i=0; i<newG->maxGrade; i++) {
! 486: ps = newG->polys[i];
! 487: for (j=0; j<ps->size; j++) {
! 488: if (ps->serial[j] >= 0) {
! 489: ii = ps->serial[j];
! 490: getArrayOfPOLY(mat,ii) = ps->syz[j]->cf;
! 491: }
! 492: }
! 493: }
! 494: return(mat);
! 495: }
! 496:
! 497:
! 498:
! 499: /* Syzygy from E-CB */
! 500: struct matrixOfPOLY *getSyzygy1(b,nc,dc)
! 501: struct matrixOfPOLY *b;
! 502: struct matrixOfPOLY *nc;
! 503: struct arrayOfPOLY *dc;
! 504: {
! 505: int m,n2,n;
! 506: struct matrixOfPOLY *mat;
! 507: int i,j,k;
! 508: POLY r;
! 509: POLY tmp;
! 510:
! 511: m = nc->m;
! 512: n2 = nc->n;
! 513: n = b->n;
! 514: mat = newMatrixOfPOLY(m,n);
! 515: for (i=0; i<m; i++) {
! 516: for (j=0; j<n; j++) {
! 517: r = ZERO;
! 518: if (i == j) {
! 519: r = getArrayOfPOLY(dc,i);
! 520: }
! 521: for (k=0; k<n2; k++) {
! 522: tmp = ppMult(getMatrixOfPOLY(nc,i,k),getMatrixOfPOLY(b,k,j));
! 523: r = ppAdd(r,tmp);
! 524: }
! 525: getMatrixOfPOLY(mat,i,j) = r;
! 526: }
! 527: }
! 528: return(mat);
! 529: }
! 530:
! 531: static int isZeroRow(mat,i)
! 532: struct matrixOfPOLY *mat;
! 533: int i;
! 534: {
! 535: int n,j;
! 536: n = mat->n;
! 537: for (j=0; j<n; j++) {
! 538: if (getMatrixOfPOLY(mat,i,j) != ZERO) return(0);
! 539: }
! 540: return(1);
! 541: }
! 542:
! 543: void errorSyz0(s)
! 544: char *s;
! 545: {
! 546: fprintf(stderr,"Error(syz0.c): %s \n",s);
! 547: exit(10);
! 548: }
! 549:
! 550:
! 551: static void printMatrixOfPOLY(mat)
! 552: struct matrixOfPOLY *mat;
! 553: {
! 554: int n,m,i,j;
! 555: POLY f;
! 556: m = mat->m; n = mat->n;
! 557: for (i=0; i<m; i++) {
! 558: for (j=0; j<n; j++) {
! 559: f = getMatrixOfPOLY(mat,i,j);
! 560: printf("%s, ",POLYToString(f,'*',1));
! 561: }
! 562: printf("\n");
! 563: }
! 564: printf("\n\n");
! 565: }
! 566:
! 567: static void printArrayOfPOLY(mat)
! 568: struct arrayOfPOLY *mat;
! 569: {
! 570: int n,m,i,j;
! 571: POLY f;
! 572: m = mat->n;
! 573: for (i=0; i<m; i++) {
! 574: f = getArrayOfPOLY(mat,i);
! 575: printf("%s, ",POLYToString(f,'*',1));
! 576: }
! 577: printf("\n\n");
! 578: }
! 579:
! 580: struct matrixOfPOLY *getSyzygy01(struct gradedPolySet *reducedBasis,
! 581: struct pair *excludePairs)
! 582: {
! 583: int r;
! 584: struct gradedPolySet *g;
! 585: struct gradedPairs *d;
! 586: int i,j;
! 587: int grade,indx;
! 588: POLY gt;
! 589: struct pair *top;
! 590: int ig,ii,jg,ji;
! 591: POLY gi,gj;
! 592: struct spValue h;
! 593: struct syz0 syz;
! 594: int pgrade = 0;
! 595: POLY rd;
! 596: POLY syzPoly;
! 597: POLY syzCf;
! 598: struct syz0 *syzp;
! 599: int serial;
! 600: struct pair *listP;
! 601: struct pair *listPfirst;
! 602: extern int Verbose;
! 603: extern struct ring *CurrentRingp;
! 604: struct polySet *ps;
! 605: int size;
! 606: int listPsize = 0;
! 607: struct arrayOfPOLY *vec;
! 608: struct matrixOfPOLY *mp;
! 609: extern int Sugar;
! 610:
! 611:
! 612: if (Debugsyz0) printf("--------------------- getSyzygy01 -----------\n");
! 613: if (reducedBasis == (struct gradedPolySet *)NULL)
! 614: return((struct matrixOfPOLY *)NULL);
! 615:
! 616: g = newGradedPolySet(reducedBasis->maxGrade+1);
! 617: for (i=0; i<g->lim; i++) { g->polys[i] = newPolySet(1); }
! 618: d = newGradedPairs((reducedBasis->maxGrade)*(reducedBasis->maxGrade)+1);
! 619: for (i=0; i< reducedBasis->maxGrade; i++) {
! 620: ps = reducedBasis->polys[i];
! 621: for (j=0; j< ps->size; j++) {
! 622: gt = ps->g[j];
! 623: grade=-1;whereInG(g,gt,&grade,&indx,Sugar);
! 624: d = updatePairs(d,gt,grade,indx,g);
! 625: g = putPolyInG(g,gt,grade,indx,newSyz0(),1,ps->serial[j]);
! 626: if (Debugsyz0) {
! 627: outputGradedPairs(d); outputGradedPolySet(g,1);
! 628: }
! 629: }
! 630: }
! 631:
! 632: listP = newPair((struct pair *)NULL); /* A list of syzygies will be stored*/
! 633: listPfirst = listP;
! 634:
! 635: while ((top = getPair(d)) != (struct pair *)NULL) {
! 636: if (positionInPairs(top,excludePairs) < 0) {
! 637: ig = top->ig; ii = top->ii; /* [ig,ii] */
! 638: jg = top->jg; ji = top->ji; /* [jg,ji] */
! 639: gi = g->polys[ig]->g[ii];
! 640: gj = g->polys[jg]->g[ji];
! 641:
! 642: h = (*sp)(gi,gj);
! 643: rd = ppAddv(ppMult(h.a,gi),ppMult(h.b,gj));
! 644: rd = (*reduction)(rd,g,1,&syz);
! 645: syzPoly = syz.syz;
! 646: syzCf = syz.cf;
! 647:
! 648: if (KanGBmessage) {
! 649: if (pgrade != top->grade) {
! 650: pgrade = top->grade;
! 651: printf(" %d",pgrade);
! 652: fflush(stdout);
! 653: }else{
! 654: if (rd ISZERO) {
! 655: printf("o"); fflush(stdout);
! 656: }else{
! 657: printf("."); fflush(stdout);
! 658: }
! 659: }
! 660: }
! 661:
! 662: if (!(rd ISZERO)) {
! 663: fprintf(stderr,"The given argument of getSyzygy01 is not a g-basis.\n");
! 664: return((struct matrixOfPOLY *)NULL);
! 665: }else{
! 666: top->syz = ppAdd(toSyzPoly(h.a,ig,ii),toSyzPoly(h.b,jg,ji));
! 667: top->syz = cpMult(toSyzCoeff(syzCf),top->syz);
! 668: top->syz = ppAdd(top->syz,syzPoly);
! 669: listP->next = top; top->prev = listP; listP = listP->next;
! 670: listPsize++;
! 671: }
! 672: }
! 673: }
! 674:
! 675: if (KanGBmessage) { printf("c"); fflush(stdout); }
! 676:
! 677: size = getSize00OfGradedSet(g);
! 678: listPfirst = listPfirst->next; /* It is the true top of sygyzies. */
! 679: mp = newMatrixOfPOLY(listPsize,size);
! 680: for (i=0; i<listPsize; i++) {
! 681: vec = syzPolyToArrayOfPOLY(size,listPfirst->syz,g);
! 682: for (j=0; j<size; j++) {
! 683: getMatrixOfPOLY(mp,i,j) = getArrayOfPOLY(vec,j);
! 684: }
! 685: listPfirst = listPfirst->next;
! 686: }
! 687: return(mp);
! 688:
! 689:
! 690: }
! 691:
! 692: static int getSize00OfGradedSet(struct gradedPolySet *g) {
! 693: int size;
! 694: int i,j;
! 695: struct polySet *ps;
! 696: size = 0;
! 697: if (g == (struct gradedPolySet *)NULL) return(0);
! 698: for (i=0; i<g->maxGrade; i++) {
! 699: ps = g->polys[i];
! 700: for (j=0; j<ps->size; j++) {
! 701: if (ps->del[j] == 0) {
! 702: size += 1;
! 703: }
! 704: }
! 705: }
! 706: return(size);
! 707: }
! 708:
! 709: static int positionInPairs(struct pair *top, struct pair *pairs) {
! 710: struct pair *tmp;
! 711: int pos;
! 712: tmp = pairs;
! 713: pos = 0;
! 714: if (top == (struct pair *)NULL) return(-1);
! 715: while (tmp != (struct pair *)NULL) {
! 716: if (((top->ig == tmp->ig) && (top->ii == tmp->ii) &&
! 717: (top->jg == tmp->jg) && (top->ji == tmp->ji)) ||
! 718: ((top->ig == tmp->jg) && (top->ii == tmp->ji) &&
! 719: (top->jg == tmp->ig) && (top->ji == tmp->ii))) {
! 720: return(pos);
! 721: }
! 722: pos++;
! 723: tmp = tmp->next;
! 724: }
! 725: return(-1);
! 726: }
! 727:
! 728: static struct pair *oldPairsToNewPairs(struct pair *opairs,
! 729: int *table,int size) {
! 730: /* Never loop up prev field. */
! 731: int ig,ii,jg,ji;
! 732: int p,q;
! 733: struct pair *ans;
! 734:
! 735: if (opairs == (struct pair *)NULL) return((struct pair *)NULL);
! 736: ig = opairs->ig; ii = opairs->ii;
! 737: jg = opairs->jg; ji = opairs->ji;
! 738: for (p=0; p<size; p++) {
! 739: if (table[4*p] == ig && table[4*p+1] == ii ) {
! 740: for (q = 0; q<size; q++) {
! 741: if (table[4*q] == jg && table[4*q+1] == ji) {
! 742: ans = newPair(NULL);
! 743: *ans = *opairs;
! 744: ans->prev = NULL;
! 745: ans->ig = table[4*p+2]; ans->ii = table[4*p+3];
! 746: ans->jg = table[4*q+2]; ans->ji = table[4*q+3];
! 747: ans->next = oldPairsToNewPairs(opairs->next,table,size);
! 748: return(ans);
! 749: }
! 750: }
! 751: }
! 752: }
! 753: return(oldPairsToNewPairs(opairs->next,table,size));
! 754: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>