Annotation of OpenXM/src/kan96xx/Kan/poly4.c, Revision 1.1
1.1 ! maekawa 1: #include <stdio.h>
! 2: #include "datatype.h"
! 3: #include "stackm.h"
! 4: #include "extern.h"
! 5: #include "extern2.h"
! 6: #include "matrix.h"
! 7: static void shell(int v[],int n);
! 8: static int degreeOfPrincipalPart(POLY f);
! 9: static int degreeOfInitW(POLY f,int w[]);
! 10:
! 11:
! 12: static void shell(v,n)
! 13: int v[];
! 14: int n;
! 15: {
! 16: int gap,i,j,temp;
! 17:
! 18: for (gap = n/2; gap > 0; gap /= 2) {
! 19: for (i = gap; i<n; i++) {
! 20: for (j=i-gap ; j>=0 && v[j]<v[j+gap]; j -= gap) {
! 21: temp = v[j];
! 22: v[j] = v[j+gap];
! 23: v[j+gap] = temp;
! 24: }
! 25: }
! 26: }
! 27: }
! 28:
! 29:
! 30: struct matrixOfPOLY *parts(f,v)
! 31: POLY f;
! 32: POLY v; /* v must be a single variable, e.g. x */
! 33: {
! 34: struct matrixOfPOLY *evPoly;
! 35: int vi = 0; /* index of v */
! 36: int vx = 1; /* x --> 1, D--> 0 */
! 37: int n,evSize,i,k,e;
! 38: int *ev;
! 39: struct object *evList;
! 40: struct object *list;
! 41: struct object ob;
! 42: POLY ans;
! 43: POLY h;
! 44: extern struct ring *CurrentRingp;
! 45: POLY ft;
! 46:
! 47:
! 48: if (f ISZERO || v ISZERO) {
! 49: evPoly = newMatrixOfPOLY(2,1);
! 50: getMatrixOfPOLY(evPoly,0,0) = ZERO;
! 51: getMatrixOfPOLY(evPoly,1,0) = ZERO;
! 52: return(evPoly);
! 53: }
! 54: n = v->m->ringp->n;
! 55: /* get the index of the variable v */
! 56: for (i=0; i<n; i++) {
! 57: if (v->m->e[i].x) {
! 58: vx = 1; vi = i; break;
! 59: }else if (v->m->e[i].D) {
! 60: vx = 0; vi = i; break;
! 61: }
! 62: }
! 63: ft = f;
! 64: /* get the vector of exponents */
! 65: evList = NULLLIST;
! 66: while (ft != POLYNULL) {
! 67: if (vx) {
! 68: e = ft->m->e[vi].x;
! 69: }else{
! 70: e = ft->m->e[vi].D;
! 71: }
! 72: ft = ft->next;
! 73: ob = KpoInteger(e);
! 74: if (!memberQ(evList,ob)) {
! 75: list = newList(&ob);
! 76: evList = vJoin(evList,list);
! 77: }
! 78: }
! 79: /*printf("evList = "); printObjectList(evList);*/
! 80: evSize = klength(evList);
! 81: ev = (int *)sGC_malloc(sizeof(int)*(evSize+1));
! 82: if (ev == (int *)NULL) errorPoly("No more memory.");
! 83: for (i=0; i<evSize; i++) {
! 84: ev[i] = KopInteger(car(evList));
! 85: evList = cdr(evList);
! 86: }
! 87: /* sort ev */
! 88: shell(ev,evSize);
! 89:
! 90: /* get coefficients */
! 91: evPoly = newMatrixOfPOLY(2,evSize);
! 92: for (i=0; i<evSize; i++) {
! 93: ans = ZERO;
! 94: getMatrixOfPOLY(evPoly,0,i) = cxx(ev[i],0,0,CurrentRingp);
! 95: ft = f;
! 96: while (ft != POLYNULL) {
! 97: if (vx) {
! 98: if (ft->m->e[vi].x == ev[i]) {
! 99: h = newCell(ft->coeffp,monomialCopy(ft->m));
! 100: xset0(h,vi); /* touch monomial part, so you need to copy it above. */
! 101: ans = ppAdd(ans,h);
! 102: }
! 103: }else{
! 104: if (ft->m->e[vi].D == ev[i]) {
! 105: h = newCell(ft->coeffp,monomialCopy(ft->m));
! 106: dset0(h,vi);
! 107: ans = ppAdd(ans,h);
! 108: }
! 109: }
! 110: ft = ft->next;
! 111: }
! 112: getMatrixOfPOLY(evPoly,1,i) = ans;
! 113: }
! 114: return(evPoly);
! 115: }
! 116:
! 117: struct object parts2(f,v)
! 118: POLY f;
! 119: POLY v; /* v must be a single variable, e.g. x */
! 120: {
! 121: struct matrixOfPOLY *evPoly;
! 122: int vi = 0; /* index of v */
! 123: int vx = 1; /* x --> 1, D--> 0 */
! 124: int n,evSize,i,k,e;
! 125: int *ev;
! 126: struct object *evList;
! 127: struct object *list;
! 128: struct object ob;
! 129: POLY ans;
! 130: POLY h;
! 131: POLY ft;
! 132: struct object ob1,ob2,rob;
! 133:
! 134:
! 135: if (f ISZERO || v ISZERO) {
! 136: evPoly = newMatrixOfPOLY(2,1);
! 137: getMatrixOfPOLY(evPoly,0,0) = ZERO;
! 138: getMatrixOfPOLY(evPoly,1,0) = ZERO;
! 139: rob = newObjectArray(2);
! 140: ob1 = newObjectArray(1);
! 141: ob2 = newObjectArray(1);
! 142: putoa(ob1,0,KpoInteger(0));
! 143: putoa(ob2,0,KpoPOLY(POLYNULL));
! 144: putoa(rob,0,ob1); putoa(rob,1,ob2);
! 145: return(rob);
! 146: }
! 147: n = v->m->ringp->n;
! 148: /* get the index of the variable v */
! 149: for (i=0; i<n; i++) {
! 150: if (v->m->e[i].x) {
! 151: vx = 1; vi = i; break;
! 152: }else if (v->m->e[i].D) {
! 153: vx = 0; vi = i; break;
! 154: }
! 155: }
! 156: ft = f;
! 157: /* get the vector of exponents */
! 158: evList = NULLLIST;
! 159: while (ft != POLYNULL) {
! 160: if (vx) {
! 161: e = ft->m->e[vi].x;
! 162: }else{
! 163: e = ft->m->e[vi].D;
! 164: }
! 165: ft = ft->next;
! 166: ob = KpoInteger(e);
! 167: if (!memberQ(evList,ob)) {
! 168: list = newList(&ob);
! 169: evList = vJoin(evList,list);
! 170: }
! 171: }
! 172: /*printf("evList = "); printObjectList(evList);*/
! 173: evSize = klength(evList);
! 174: ev = (int *)sGC_malloc(sizeof(int)*(evSize+1));
! 175: if (ev == (int *)NULL) errorPoly("No more memory.");
! 176: for (i=0; i<evSize; i++) {
! 177: ev[i] = KopInteger(car(evList));
! 178: evList = cdr(evList);
! 179: }
! 180: /* sort ev */
! 181: shell(ev,evSize);
! 182:
! 183: /* get coefficients */
! 184: evPoly = newMatrixOfPOLY(2,evSize);
! 185: for (i=0; i<evSize; i++) {
! 186: ans = ZERO;
! 187: /* getMatrixOfPOLY(evPoly,0,i) = cxx(ev[i],0,0,CurrentRingp); */
! 188: getMatrixOfPOLY(evPoly,0,i) = POLYNULL;
! 189: ft = f;
! 190: while (ft != POLYNULL) {
! 191: if (vx) {
! 192: if (ft->m->e[vi].x == ev[i]) {
! 193: h = newCell(ft->coeffp,monomialCopy(ft->m));
! 194: xset0(h,vi); /* touch monomial part, so you need to copy it above. */
! 195: ans = ppAdd(ans,h);
! 196: }
! 197: }else{
! 198: if (ft->m->e[vi].D == ev[i]) {
! 199: h = newCell(ft->coeffp,monomialCopy(ft->m));
! 200: dset0(h,vi);
! 201: ans = ppAdd(ans,h);
! 202: }
! 203: }
! 204: ft = ft->next;
! 205: }
! 206: getMatrixOfPOLY(evPoly,1,i) = ans;
! 207: }
! 208: rob = newObjectArray(2);
! 209: ob1 = newObjectArray(evSize);
! 210: ob2 = newObjectArray(evSize);
! 211: for (i=0; i<evSize; i++) {
! 212: putoa(ob2,i,KpoPOLY(getMatrixOfPOLY(evPoly,1,i)));
! 213: putoa(ob1,i,KpoInteger(ev[i]));
! 214: }
! 215: putoa(rob,0,ob1); putoa(rob,1,ob2);
! 216: return(rob);
! 217: }
! 218:
! 219: int pDegreeWrtV(f,v)
! 220: POLY f;
! 221: POLY v;
! 222: {
! 223: int vx = 1;
! 224: int vi = 0;
! 225: int i,n;
! 226: int ans;
! 227: if (f ISZERO || v ISZERO) return(0);
! 228: n = f->m->ringp->n;
! 229: for (i=0; i<n; i++) {
! 230: if (v->m->e[i].x) {
! 231: vx = 1; vi = i;
! 232: break;
! 233: }else if (v->m->e[i].D) {
! 234: vx = 0; vi = i;
! 235: break;
! 236: }
! 237: }
! 238: if (vx) {
! 239: ans = f->m->e[vi].x;
! 240: }else{
! 241: ans = f->m->e[vi].D;
! 242: }
! 243: f = f->next;
! 244: while (f != POLYNULL) {
! 245: if (vx) {
! 246: if (f->m->e[vi].x > ans) ans = f->m->e[vi].x;
! 247: }else{
! 248: if (f->m->e[vi].D > ans) ans = f->m->e[vi].D;
! 249: }
! 250: f = f->next;
! 251: }
! 252: return(ans);
! 253: }
! 254:
! 255: int containVectorVariable(POLY f)
! 256: {
! 257: MONOMIAL tf;
! 258: static int nn,mm,ll,cc,n,m,l,c;
! 259: static struct ring *cr = (struct ring *)NULL;
! 260: int i;
! 261:
! 262: if (f ISZERO) return(0);
! 263: tf = f->m;
! 264: if (tf->ringp != cr) {
! 265: n = tf->ringp->n;
! 266: m = tf->ringp->m;
! 267: l = tf->ringp->l;
! 268: c = tf->ringp->c;
! 269: nn = tf->ringp->nn;
! 270: mm = tf->ringp->mm;
! 271: ll = tf->ringp->ll;
! 272: cc = tf->ringp->cc;
! 273: cr = tf->ringp;
! 274: }
! 275:
! 276: while (f != POLYNULL) {
! 277: tf = f->m;
! 278: for (i=cc; i<c; i++) {
! 279: if ( tf->e[i].x ) return(1);
! 280: if ( tf->e[i].D ) return(1);
! 281: }
! 282: for (i=ll; i<l; i++) {
! 283: if (tf->e[i].x) return(1);
! 284: if (tf->e[i].D) return(1);
! 285: }
! 286: for (i=mm; i<m; i++) {
! 287: if (tf->e[i].x) return(1);
! 288: if (tf->e[i].D) return(1);
! 289: }
! 290: for (i=nn; i<n; i++) {
! 291: if (tf->e[i].x) return(1);
! 292: if (tf->e[i].D) return(1);
! 293: }
! 294: f = f->next;
! 295: }
! 296: return(0);
! 297:
! 298: }
! 299:
! 300: POLY homogenize(f)
! 301: POLY f;
! 302: /* homogenize by using (*grade)(f) */
! 303: {
! 304: POLY t;
! 305: int maxg;
! 306: int flag,d;
! 307:
! 308: if (f == ZERO) return(f);
! 309: t = f; maxg = (*grade)(f); flag = 0;
! 310: while (t != POLYNULL) {
! 311: d = (*grade)(t);
! 312: if (d != maxg) flag = 1;
! 313: if (d > maxg) {
! 314: maxg = d;
! 315: }
! 316: t = t->next;
! 317: }
! 318: if (flag == 0) return(f);
! 319:
! 320: f = pmCopy(f); /* You can rewrite the monomial parts */
! 321: t = f;
! 322: while (t != POLYNULL) {
! 323: d = (*grade)(t);
! 324: if (d != maxg) {
! 325: t->m->e[0].D += maxg-d; /* Multiply h^(maxg-d) */
! 326: }
! 327: t = t->next;
! 328: }
! 329: return(f);
! 330: }
! 331:
! 332: int isHomogenized(f)
! 333: POLY f;
! 334: {
! 335: POLY t;
! 336: extern int Homogenize_vec;
! 337: int maxg;
! 338: if (!Homogenize_vec) return(isHomogenized_vec(f));
! 339: if (f == ZERO) return(1);
! 340: maxg = (*grade)(f);
! 341: t = f;
! 342: while (t != POLYNULL) {
! 343: if (maxg != (*grade)(t)) return(0);
! 344: t = t->next;
! 345: }
! 346: return(1);
! 347: }
! 348:
! 349: int isHomogenized_vec(f)
! 350: POLY f;
! 351: {
! 352: /* This is not efficient version. *grade should be grade_module1v(). */
! 353: POLY t;
! 354: int ggg;
! 355: if (f == ZERO) return(1);
! 356: while (f != POLYNULL) {
! 357: t = f;
! 358: ggg = (*grade)(f);
! 359: while (t != POLYNULL) {
! 360: if ((*isSameComponent)(f,t)) {
! 361: if (ggg != (*grade)(t)) return(0);
! 362: }
! 363: t = t->next;
! 364: }
! 365: f = f->next;
! 366: }
! 367: return(1);
! 368: }
! 369:
! 370:
! 371: static int degreeOfPrincipalPart(f)
! 372: POLY f;
! 373: {
! 374: int n,i,dd;
! 375: if (f ISZERO) return(0);
! 376: n = f->m->ringp->n; dd = 0;
! 377: /* D[0] is homogenization var */
! 378: for (i=1; i<n; i++) {
! 379: dd += f->m->e[i].D;
! 380: }
! 381: return(dd);
! 382: }
! 383:
! 384: POLY POLYToPrincipalPart(f)
! 385: POLY f;
! 386: {
! 387: POLY node;
! 388: struct listPoly nod;
! 389: POLY h;
! 390: POLY g;
! 391: int maxd = -20000; /* very big negative number */
! 392: int dd;
! 393: node = &nod; node->next = POLYNULL; h = node;
! 394:
! 395: g = pCopy(f); /* shallow copy */
! 396: while (!(f ISZERO)) {
! 397: dd = degreeOfPrincipalPart(f);
! 398: if (dd > maxd) maxd = dd;
! 399: f = f->next;
! 400: }
! 401: while (!(g ISZERO)) {
! 402: dd = degreeOfPrincipalPart(g);
! 403: if (dd == maxd) {
! 404: h->next = g;
! 405: h = h->next;
! 406: }
! 407: g = g->next;
! 408: }
! 409: h->next = POLYNULL;
! 410: return(node->next);
! 411: }
! 412:
! 413: static int degreeOfInitW(f,w)
! 414: POLY f;
! 415: int w[];
! 416: {
! 417: int n,i,dd;
! 418: if (f ISZERO) {
! 419: errorPoly("degreeOfInitW(0,w) ");
! 420: }
! 421: n = f->m->ringp->n; dd = 0;
! 422: for (i=0; i<n; i++) {
! 423: dd += (f->m->e[i].D)*w[n+i];
! 424: dd += (f->m->e[i].x)*w[i];
! 425: }
! 426: return(dd);
! 427: }
! 428:
! 429: POLY POLYToInitW(f,w)
! 430: POLY f;
! 431: int w[]; /* weight vector */
! 432: {
! 433: POLY node;
! 434: struct listPoly nod;
! 435: POLY h;
! 436: POLY g;
! 437: int maxd;
! 438: int dd;
! 439: node = &nod; node->next = POLYNULL; h = node;
! 440:
! 441: if (f ISZERO) return(f);
! 442: maxd = degreeOfInitW(f,w);
! 443: g = pCopy(f); /* shallow copy */
! 444: while (!(f ISZERO)) {
! 445: dd = degreeOfInitW(f,w);
! 446: if (dd > maxd) maxd = dd;
! 447: f = f->next;
! 448: }
! 449: while (!(g ISZERO)) {
! 450: dd = degreeOfInitW(g,w);
! 451: if (dd == maxd) {
! 452: h->next = g;
! 453: h = h->next;
! 454: }
! 455: g = g->next;
! 456: }
! 457: h->next = POLYNULL;
! 458: return(node->next);
! 459: }
! 460:
! 461:
! 462: /*
! 463: 1.The substitution "ringp->multiplication = ...." is allowed only in
! 464: KsetUpRing(), so the check in KswitchFunction is not necessary.
! 465: 2.mmLarger != matrix and AvoidTheSameRing==1, then error
! 466: 3.If Schreyer = 1, then the system always generates a new ring.
! 467: 4.The execution of set_order_by_matrix is not allowed when Avoid... == 1.
! 468: 5.When mmLarger == tower (in tower.sm1, tower-sugar.sm1), we do
! 469: as follows with our own risk.
! 470: [(AvoidTheSameRing)] pushEnv [ [(AvoidTheSameRing) 0] system_variable (mmLarger) (tower) switch_function ] pop popEnv
! 471: */
! 472: int isTheSameRing(struct ring *rstack[],int rp, struct ring *newRingp)
! 473: {
! 474: struct ring *rrr;
! 475: int i,j,k;
! 476: int a=0;
! 477: for (k=0; k<rp; k++) {
! 478: rrr = rstack[k];
! 479: if (rrr->p != newRingp->p) { a=1; goto bbb ; }
! 480: if (rrr->n != newRingp->n) { a=2; goto bbb ; }
! 481: if (rrr->nn != newRingp->nn) { a=3; goto bbb ; }
! 482: if (rrr->m != newRingp->m) { a=4; goto bbb ; }
! 483: if (rrr->mm != newRingp->mm) { a=5; goto bbb ; }
! 484: if (rrr->l != newRingp->l) { a=6; goto bbb ; }
! 485: if (rrr->ll != newRingp->ll) { a=7; goto bbb ; }
! 486: if (rrr->c != newRingp->c) { a=8; goto bbb ; }
! 487: if (rrr->cc != newRingp->cc) { a=9; goto bbb ; }
! 488: for (i=0; i<rrr->n; i++) {
! 489: if (strcmp(rrr->x[i],newRingp->x[i])!=0) { a=10; goto bbb ; }
! 490: if (strcmp(rrr->D[i],newRingp->D[i])!=0) { a=11; goto bbb ; }
! 491: }
! 492: if (rrr->orderMatrixSize != newRingp->orderMatrixSize) { a=12; goto bbb ; }
! 493: for (i=0; i<rrr->orderMatrixSize; i++) {
! 494: for (j=0; j<2*(rrr->n); j++) {
! 495: if (rrr->order[i*2*(rrr->n)+j] != newRingp->order[i*2*(rrr->n)+j])
! 496: { a=13; goto bbb ; }
! 497: }
! 498: }
! 499: if (rrr->next != newRingp->next) { a=14; goto bbb ; }
! 500: if (rrr->multiplication != newRingp->multiplication) { a=15; goto bbb ; }
! 501: /* if (rrr->schreyer != newRingp->schreyer) { a=16; goto bbb ; }*/
! 502: if (newRingp->schreyer == 1) { a=16; goto bbb; }
! 503: /* The following fields are ignored.
! 504: void *gbListTower;
! 505: int *outputOrder;
! 506: char *name;
! 507: */
! 508: /* All tests are passed. */
! 509: return(k);
! 510: bbb: ;
! 511: /* for debugging. */
! 512: /* fprintf(stderr," reason=%d, ",a); */
! 513: }
! 514: return(-1);
! 515: }
! 516:
! 517:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>