Annotation of OpenXM/src/kan96xx/Kan/hilbert.c, Revision 1.1
1.1 ! maekawa 1: /* hilbert.c
! 2: 1992/06/16
! 3: 1992/06/18
! 4: 1993/04/28
! 5: \bibitem{Cocoa-Hilbert}
! 6: Bigatti, A., Caboara, M., Robianno, L. (1991):
! 7: On the computation of Hilbert Poincar\'e series.
! 8: Applicable algebra in Engineering, Communication and Computing
! 9: {\bf 2}, 21--33.
! 10: 1998/10/31
! 11: */
! 12: #include <stdio.h>
! 13: #include "datatype.h"
! 14: #include "stackm.h"
! 15: #include "extern.h"
! 16: #include "extern2.h"
! 17:
! 18: struct arrayOfPOLYold {
! 19: int n;
! 20: POLY *array;
! 21: };
! 22:
! 23: static int M;
! 24: static int MM;
! 25: static int NN;
! 26: static int N;
! 27:
! 28: static struct ring *hringp;
! 29:
! 30: #define max(a,b) ((a)>(b)? a: b)
! 31:
! 32: POLY hilbert2(POLY k,int p,POLY pArray[]);
! 33: static POLY xx(int i,int k);
! 34: static void pWriteln(POLY f)
! 35: { printf("%s\n",POLYToString(f,' ',0)); fflush(stdout);}
! 36:
! 37: struct object hilberto(struct object obgb,struct object obv)
! 38: {
! 39: int m; /* size of the base. */
! 40: int n; /* number of variables */
! 41: int i,j,k;
! 42: int n0; /* number of the variables for the polynomials in base. */
! 43: struct object obf;
! 44: struct ring *rp;
! 45: struct ring *rr = (struct ring *)NULL;
! 46: POLY *base;
! 47: int *x;
! 48: int *d;
! 49: int debug = 0;
! 50: POLY f;
! 51: POLY g;
! 52: struct object rob;
! 53: struct object ccc;
! 54: extern struct ring SmallRing;
! 55: int worg;
! 56: extern int WarningNoVectorVariable;
! 57:
! 58: rob = NullObject;
! 59:
! 60: if (obgb.tag != Sarray) {
! 61: errorKan1("%s\n","Khilbert(): The first argument must be array of poly.");
! 62: }
! 63: if (obv.tag != Sarray) {
! 64: errorKan1("%s\n","Khilbert(): The second argument must be array of variables.");
! 65: }
! 66: m = getoaSize(obgb);
! 67: n = getoaSize(obv);
! 68: if (n <= 0) errorKan1("%s\n","Khilbert(): The number of variables must be more than 0.");
! 69: if (m <= 0) errorKan1("%s\n","Khilbert(): The number of basis must be more than 0.");
! 70: for (i=0; i<m; i++) {
! 71: obf = getoa(obgb,i);
! 72: if (obf.tag != Spoly) {
! 73: errorKan1("%s\n","Khilbert(): The element of the first array must be a polynomials.\n");
! 74: }
! 75: if (KopPOLY(obf) ISZERO) {
! 76: errorKan1("%s\n","Khilbert(): The element of the first array must not be the zero.\n");
! 77: }
! 78: f = KopPOLY(obf);
! 79: if (rr != (struct ring*)NULL) {
! 80: if (rr != f->m->ringp) {
! 81: errorKan1("%s\n","Hhilbert(): the element must belong to the same ring.");
! 82: }
! 83: }else{
! 84: rr = f->m->ringp;
! 85: }
! 86: }
! 87: for (i=0; i<n; i++) {
! 88: obf = getoa(obv,i);
! 89: if (obf.tag != Spoly) {
! 90: errorKan1("%s\n","Khilbert(): The element of the second array must be a polynomials.\n");
! 91: }
! 92: if (KopPOLY(obf) ISZERO) {
! 93: errorKan1("%s\n","Khilbert(): The element of the second array must not be the zero.\n");
! 94: }
! 95: f = KopPOLY(obf);
! 96: if (rr != (struct ring*)NULL) {
! 97: if (rr != f->m->ringp) {
! 98: errorKan1("%s\n","Hhilbert(): the element must belong to the same ring.");
! 99: }
! 100: }else{
! 101: rr = f->m->ringp;
! 102: }
! 103: }
! 104:
! 105: worg = WarningNoVectorVariable;
! 106: WarningNoVectorVariable = 0;
! 107: rp = KopRingp(KdefaultPolyRing(KpoInteger(n)));
! 108: WarningNoVectorVariable = worg;
! 109:
! 110: hringp = rp;
! 111: MM=0; M=0; NN=n; N = n;
! 112:
! 113: base = (POLY *)sGC_malloc(sizeof(POLY)*(m+1));
! 114: if (base == (POLY *)NULL) errorKan1("%s\n","no more memory.");
! 115: obf = getoa(obgb,0);
! 116: f = KopPOLY(obf);
! 117: n0 = f->m->ringp->n;
! 118: x = (int *)sGC_malloc(sizeof(int)*n0);
! 119: d = (int *)sGC_malloc(sizeof(int)*n0);
! 120: if (d == (int *)NULL) errorKan1("%s\n","no more memory.");
! 121:
! 122: for (i=0; i<n0; i++) {x[i] = d[i] = -1;}
! 123: for (i=0; i<n; i++) {
! 124: obf = getoa(obv,i);
! 125: f = KopPOLY(obf);
! 126: for (j=0; j<n0; j++) {
! 127: if ((f->m->e)[j].x) {
! 128: x[j] = 1; break;
! 129: }
! 130: if ((f->m->e)[j].D) {
! 131: d[j] = 1; break;
! 132: }
! 133: }
! 134: }
! 135:
! 136: j = 0;
! 137: for (i=0; i<n0; i++) {
! 138: if (x[i] != -1) x[i] = j++;
! 139: }
! 140: for (i=0; i<n0; i++) {
! 141: if (d[i] != -1) d[i] = j++;
! 142: }
! 143:
! 144: if (debug) {
! 145: for (i=0; i<n0; i++) {
! 146: printf("x[%d]=%d ",i,x[i]);
! 147: }
! 148: for (i=0; i<n0; i++) {
! 149: printf("d[%d]=%d ",i,d[i]);
! 150: }
! 151: printf("\n");
! 152: }
! 153:
! 154: for (i=0; i<m; i++) {
! 155: obf = getoa(obgb,i);
! 156: f = KopPOLY(obf);
! 157: g = cxx(1,0,0,rp);
! 158: for (j=0; j<n0; j++) {
! 159: if ((f->m->e)[j].x) {
! 160: if (x[j] != -1) {
! 161: (g->m->e)[x[j]].D = (f->m->e)[j].x;
! 162: }
! 163: }
! 164: if ((f->m->e)[j].D) {
! 165: if (d[j] != -1) {
! 166: (g->m->e)[d[j]].D = (f->m->e)[j].D;
! 167: }
! 168: }
! 169: }
! 170: base[i] = g;
! 171: }
! 172:
! 173: if (debug) {
! 174: for (i=0; i<m; i++) {
! 175: printf("%s\n",POLYToString(base[i],'*',1));
! 176: }
! 177: }
! 178:
! 179: rob = newObjectArray(2);
! 180: ccc.tag = SuniversalNumber;
! 181: ccc.lc.universalNumber = intToCoeff(NN-M,&SmallRing);
! 182: putoa(rob,0,ccc);
! 183: putoa(rob,1,KpoPOLY(hilbert2(cdd(1,0,1,hringp),m,base)));
! 184: return(rob);
! 185: }
! 186:
! 187: /*
! 188: static POLY getArrayOfPOLYold(struct arrayOfPOLYold set,int i) {
! 189: if (i < 0 || i >= set.n) {
! 190: errorKan1("%s\n","getArrayOfPOLYold(): index is out of bound.");
! 191: }
! 192: return((set.array)[i]);
! 193: }
! 194: */
! 195: #define getArrayOfPOLYold(set,i) ((set).array)[i]
! 196:
! 197: static struct arrayOfPOLYold newArrayOfPOLYold(int size) {
! 198: static struct arrayOfPOLYold a;
! 199: if (size <= 0) {
! 200: errorKan1("%s\n","newArrayOfPOLYold(): the size must be more than 0.");
! 201: }
! 202: a.n = size;
! 203: a.array = (POLY *)sGC_malloc(sizeof(POLY)*size);
! 204: if (a.array == (POLY *)NULL) errorKan1("%s\n","no more memory.");
! 205: return(a);
! 206: }
! 207:
! 208: static POLY iToPoly(int n) {
! 209: return(cxx(n,0,0,hringp));
! 210: }
! 211: static POLY xx(int i,int k) {
! 212: return(cxx(1,i,k,hringp));
! 213: }
! 214:
! 215: static int polyToInt(POLY f) {
! 216: if (f ISZERO) return(0);
! 217: return(coeffToInt(f->coeffp));
! 218: }
! 219:
! 220:
! 221: static shell(v,n)
! 222: int v[];
! 223: int n;
! 224: {
! 225: int gap,i,j,temp;
! 226:
! 227: for (gap = n/2; gap > 0; gap /= 2) {
! 228: for (i = gap; i<n; i++) {
! 229: for (j=i-gap ; j>=0 && v[j]>v[j+gap]; j -= gap) {
! 230: temp = v[j];
! 231: v[j] = v[j+gap];
! 232: v[j+gap] = temp;
! 233: }
! 234: }
! 235: }
! 236: }
! 237:
! 238: static POLY ppifac(f,n)
! 239: POLY f;
! 240: int n;
! 241: /* ppifac returns (f+n) (f+n-1) ... (f+1). */
! 242: {
! 243: POLY r;
! 244: int i;
! 245: if (n < 0) {
! 246: warningHilbert(" ppifac() is called with negative argument.");
! 247: return(POLYNULL);
! 248: }else if (n==0) {
! 249: return(pcmCopy(f));
! 250: }else {
! 251: r = iToPoly(1);
! 252: for (i=1; i<=n; i++) {
! 253: r = ppMult( ppAdd(f,iToPoly(i)), r);
! 254: }
! 255: return(r);
! 256: }
! 257: }
! 258:
! 259:
! 260: static int isReducibleD1(f,g)
! 261: POLY f;
! 262: POLY g;
! 263: {
! 264: int i;
! 265: for (i = M; i < NN; i++) {
! 266: if (((f->m->e)[i].x >= (g->m->e)[i].x) &&
! 267: ((f->m->e)[i].D >= (g->m->e)[i].D)) {
! 268: /* do nothing */
! 269: } else {
! 270: return(0);
! 271: }
! 272: }
! 273: return(1);
! 274: }
! 275:
! 276:
! 277: static struct arrayOfPOLYold
! 278: reduceMonomials(set)
! 279: struct arrayOfPOLYold set;
! 280: {
! 281: int *drop; /* 1--> yes. drop set[i]. 0--> no. */
! 282: int i,j;
! 283: int size;
! 284: int newsize;
! 285: struct arrayOfPOLYold ans;
! 286: size = set.n; /* get the size of set */
! 287:
! 288: drop = (int *)sGC_malloc(sizeof(int)*(size+1));
! 289: if (drop == (int *)NULL) errorHilbert("No more memory.");
! 290: for (i=0; i < size; i++) {
! 291: drop[i]=0;
! 292: }
! 293: /* O(n^2) algorithm to reduced basis */
! 294: for (i = 0; i < size; i++) {
! 295: if (!drop[i]) {
! 296: for (j=i+1; j < size; j++) {
! 297: if (!drop[j]) {
! 298: if (isReducibleD1(getArrayOfPOLYold(set,i),
! 299: getArrayOfPOLYold(set,j))) {
! 300: drop[i] = 1; break;
! 301: }else if (isReducibleD1(getArrayOfPOLYold(set,j),
! 302: getArrayOfPOLYold(set,i))) {
! 303: drop[j] = 1; break;
! 304: }else { }
! 305: }
! 306: }
! 307: }
! 308: }
! 309: newsize = 0;
! 310: for (i = 0; i < size; i++) {
! 311: if (!drop[i]) ++newsize;
! 312: }
! 313: ans = newArrayOfPOLYold(newsize);
! 314: j = 0;
! 315:
! 316:
! 317: for (i = 0; i < size; i++) {
! 318:
! 319:
! 320: if (!drop[i]) {
! 321: getArrayOfPOLYold(ans,j) = pcmCopy(getArrayOfPOLYold(set,i));
! 322: j++;
! 323: }
! 324: }
! 325: return(ans);
! 326: }
! 327:
! 328:
! 329: static int tryDecompose(set,i,j,vA,vWhich)
! 330: struct arrayOfPOLYold set; /* input: monomials */
! 331: int i,j; /* decompose with respect to the (i,j)-th
! 332: variable: i=0-->x_j, i=1--->D_j */
! 333: int vA[]; /* Return value: vA[0] is a_0, ... */
! 334: int vWhich[]; /* Return value: vWhich[i] is <<a>> of set[i] */
! 335: /* set ---> x_j^(a_0) I_{a_0} + .... + x_j^{a_{p-1}} I_{a_{p-1}} */
! 336: /* return value is p */
! 337: /* tryDecompose is used to find the best decomposition. */
! 338: {
! 339: int size,k,ell,ell2;
! 340: POLY f;
! 341: int p;
! 342: int *tmp;
! 343:
! 344: size = set.n;
! 345: if ( i == 0) { /* focus on x */
! 346: for (k = 0; k < size; k++) {
! 347: f = getArrayOfPOLYold(set,k);
! 348: vWhich[k] = (f->m->e)[j].x;
! 349: }
! 350: }else{ /* focus on D */
! 351: for (k = 0; k < size; k++) {
! 352: f = getArrayOfPOLYold(set,k);
! 353: vWhich[k] = (f->m->e)[j].D;
! 354: }
! 355: }
! 356: #ifdef DEBUG3
! 357: /*for (i=0; i<size; i++) printf("%3d", vWhich[i]);
! 358: printf(" vWhich\n");*/
! 359: #endif
! 360: /* sort and prune */
! 361: tmp = (int *)sGC_malloc(sizeof(int)*((size+1)*2));
! 362: if (tmp == (int *) NULL) errorHilbert("No more memory.");
! 363: for (i=0; i<size; i++) tmp[i] = vWhich[i];
! 364: shell(tmp,size);
! 365: /* prune */
! 366: p = 0; vA[p] = tmp[0]; p++;
! 367: for (i=1; i<size; i++) {
! 368: if (vA[p-1] != tmp[i]) vA[p++] = tmp[i];
! 369: }
! 370:
! 371: #ifdef DEBUG3
! 372: /*for (i=0; i<p; i++) printf("%3d", vA[i]);
! 373: printf("---vA\n");*/
! 374: #endif
! 375: return(p);
! 376: }
! 377:
! 378: static struct arrayOfPOLYold getJA(a,set,vWhich,ja,ii,xd,ith)
! 379: /* get J_{a_{i+1}} from J_{a_i}
! 380: J_{a_{i+1}} = J_{a_i} \cup I_{a_{i+1}}
! 381: */
! 382: int a; /* each a_i */
! 383: struct arrayOfPOLYold set; /* polynomials */
! 384: int *vWhich; /* vWhich[0] is exponent a of set[0], .... */
! 385: struct arrayOfPOLYold ja; /* J_i */
! 386: int ii; /* ii is i */
! 387: int xd; /* xd == 0 --> x, xd == 1 --> D */
! 388: int ith; /* x_{ith} or D_{ith} is the best variable */
! 389: {
! 390: int size;
! 391: struct arrayOfPOLYold input;
! 392: POLY f;
! 393: int start,i;
! 394:
! 395: #ifdef DEBUG3
! 396: printf("var is %d ",a);
! 397: printf("set is "); outputarrayOfPOLYold(set);
! 398: #endif
! 399:
! 400: size = set.n;
! 401: start = 0;
! 402: /* input = newArrayOfPOLYold(size); */
! 403: input = newArrayOfPOLYold(size+ja.n+1); /* 1993/09/08 */
! 404:
! 405:
! 406: if (ii != 0) {
! 407: for (i = 0; i < ja.n ; i++) {
! 408: getArrayOfPOLYold(input,start) = getArrayOfPOLYold(ja,i);
! 409: start++;
! 410: }
! 411: }
! 412: for (i = 0; i < size; i++) {
! 413: if (vWhich[i] == a) {
! 414: f = pcmCopy(getArrayOfPOLYold(set,i));
! 415: if (xd == 0) {
! 416: (f->m->e)[ith].x = 0;
! 417: }else{
! 418: (f->m->e)[ith].D = 0;
! 419: }
! 420: getArrayOfPOLYold(input,start) = f;
! 421: start++ ;
! 422: }
! 423: }
! 424: input.n = start; /* This is not good code. */
! 425: #ifdef DEBUG2
! 426: for (i=0; i<start; i++) {checkh(input,i);} /****/
! 427: #endif
! 428: #ifdef DEBUG3
! 429: printf("input is "); outputarrayOfPOLYold(input);
! 430: #endif
! 431: input= reduceMonomials(input);
! 432: #ifdef DEBUG3
! 433: /*printf("input is "); outputarrayOfPOLYold(input);*/
! 434: #endif
! 435: return( input );
! 436: }
! 437:
! 438:
! 439: static struct arrayOfPOLYold *getDecomposition(set,activeX,activeD)
! 440: struct arrayOfPOLYold set;
! 441: int activeX[N0];
! 442: int activeD[N0];
! 443: {
! 444: int i;
! 445: int size;
! 446: int p;
! 447: int *vA,*vWhich;
! 448: int xmax = 0;
! 449: int dmax = 0;
! 450: int xi,di;
! 451: int resultsize;
! 452: struct arrayOfPOLYold ja;
! 453: struct arrayOfPOLYold *ans;
! 454:
! 455: size = set.n;
! 456: vA = (int *)sGC_malloc(sizeof(int)*(size+1));
! 457: vWhich = (int *)sGC_malloc(sizeof(int)*(size+1));
! 458: if (vA == (int *)NULL || vWhich == (int *)NULL) errorHilbert("No more memory.\n");
! 459: /* find the best decomposition */
! 460: for ( i = M; i < NN; i++) {
! 461: if (activeX[i]) {
! 462: p = tryDecompose(set,0,i,vA,vWhich);
! 463: if (p > xmax) {
! 464: xmax = p;
! 465: xi = i;
! 466: if (xmax == size) {
! 467: break;
! 468: }
! 469: }
! 470: }
! 471: }
! 472: if (xmax != size) {
! 473: for ( i = M; i < NN; i++) {
! 474: if (activeD[i]) {
! 475: p = tryDecompose(set,1,i,vA,vWhich);
! 476: if (p > dmax) {
! 477: dmax = p;
! 478: di = i;
! 479: if (dmax == size) {
! 480: break;
! 481: }
! 482: }
! 483: }
! 484: }
! 485: }
! 486: /*
! 487: ans[0] = (a_0,...,a_{p-1})
! 488: ans[1] = generators of J_{a_0}
! 489: ....
! 490: ans[p] = generators of J_{a_{p-1}}, p = xmax or dmax
! 491: */
! 492: if (xmax > dmax) {
! 493: tryDecompose(set,0,xi,vA,vWhich); /* xi is the best variable */
! 494: ans = (struct arrayOfPOLYold *)sGC_malloc(sizeof(struct arrayOfPOLYold)*(xmax+1));
! 495: if (ans == (struct arrayOfPOLYold *)NULL) errorHilbert("No more memory.");
! 496: ans[0] = newArrayOfPOLYold(xmax);
! 497: for (i = 0; i < xmax; i++) {
! 498: getArrayOfPOLYold(ans[0],i) = iToPoly(vA[i]);
! 499: }
! 500: /* compute J for for each a_i */
! 501: for (i = 0; i < xmax; i++) {
! 502: ans[i+1] = getJA(vA[i],set,vWhich,ans[i],i,0,xi);
! 503: }
! 504: return(ans);
! 505: }else {
! 506: tryDecompose(set,1,di,vA,vWhich); /* di is the best variable */
! 507: ans = (struct arrayOfPOLYold *)sGC_malloc(sizeof(struct arrayOfPOLYold)*(dmax+1));
! 508: if (ans == (struct arrayOfPOLYold *)NULL) errorHilbert("No more memory.");
! 509: ans[0] = newArrayOfPOLYold(dmax);
! 510: for (i = 0; i < dmax; i++) {
! 511: getArrayOfPOLYold((ans[0]),i) = iToPoly(vA[i]);
! 512: }
! 513: /* compute J for for each a_i */
! 514: for (i = 0; i < dmax; i++) {
! 515: ans[i+1] = getJA(vA[i],set,vWhich,ans[i],i,1,di);
! 516: }
! 517: return(ans);
! 518: }
! 519: }
! 520:
! 521: static POLY hilbert1T(set)
! 522: struct arrayOfPOLYold set;
! 523: /* <<set>> must be reduced basis and each polynomial must have the length 1 */
! 524: /* Unnecessary exponents should be set to zero. For example, f = x_{M-1} x_M
! 525: is illegal input. It should be x_M ( M <= index <= NN ).
! 526: cf. hilbert1() and hilbert2()
! 527: */
! 528: /* Example: hilbert1T of x^a is
! 529: x0^0 - x0^(-a) <=== {{\ell+n} \choose n} - {{\ell+n-a} \choose n}
! 530: */
! 531: {
! 532: int activeX[N0];
! 533: int activeD[N0];
! 534: int size;
! 535: int i,j;
! 536: POLY f;
! 537: int active = 0;
! 538: int xxdd,xi,di,a;
! 539: struct arrayOfPOLYold *ja;
! 540: struct arrayOfPOLYold hja;
! 541: int p;
! 542: POLY ansp;
! 543:
! 544: #ifdef DEBUG
! 545: static int cc = 0;
! 546: /* debugging */ /***/
! 547: printf("hilbert1T: size of set is %d\n",set.n);
! 548: for (i=0; i<set.n; i++) {
! 549: printf("%s\n", POLYToString(getArrayOfPOLYold(set,i),' ',0));
! 550: }
! 551: printf("^^%d--------------------------------------------\n",cc++);
! 552: #endif
! 553:
! 554: size = set.n;
! 555: if (size < 1) {
! 556: #ifdef DEBUG
! 557: cc--; printf("<<%d>> value is 1.\n",cc); /***/
! 558: #endif
! 559: return(iToPoly(1));
! 560: }
! 561: for (i = 0; i < N; i++) {
! 562: activeX[i] = activeD[i] = 0;
! 563: }
! 564: for (i = 0; i < size; i++) {
! 565: f = getArrayOfPOLYold(set,i);
! 566: for (j = M; j < NN; j++) {
! 567: if ((f->m->e)[j].x) {
! 568: activeX[j] = 1;
! 569: }
! 570: if ((f->m->e)[j].D) {
! 571: activeD[j] = 1;
! 572: }
! 573: }
! 574: }
! 575: for (i = M; i < NN; i++) {
! 576: if (activeX[i]) {
! 577: if (active) {
! 578: active = 2;
! 579: break;
! 580: }
! 581: active = 1;
! 582: xxdd = 0; xi = i;
! 583: }
! 584: if (activeD[i]) {
! 585: if (active) {
! 586: active = 2;
! 587: break;
! 588: }
! 589: active = 1;
! 590: xxdd = 1; di = i;
! 591: }
! 592: }
! 593: if (active == 0) {
! 594: #ifdef DEBUG
! 595: cc--; printf("<<%d>> value is 0.\n",cc); /***/
! 596: #endif
! 597: return(POLYNULL);
! 598: }else if (active == 1) {
! 599: f = getArrayOfPOLYold(set,0);
! 600: if (xxdd == 0) {
! 601: a = (f->m->e)[xi].x;
! 602: #ifdef BEBUG
! 603: cc--; printf("<<%d>> 1-x0^(%d).\n",cc,a); /***/
! 604: #endif
! 605: return(ppSub(iToPoly(1),xx(0,-a)));
! 606: }else {
! 607: a = (f->m->e)[di].D;
! 608: #ifdef DEBUG
! 609: cc--; printf("<<%d>> 1-x0^(%d).\n",cc,a); /***/
! 610: #endif
! 611: return(ppSub(iToPoly(1),xx(0,-a)));
! 612: }
! 613: }
! 614:
! 615: /* compute <J_{a_0}>, ..., <J_{a_{p-1}}> */
! 616: ja = getDecomposition(set,activeX,activeD);
! 617: p = (ja[0]).n;
! 618: #ifdef DEBUG3
! 619: outputDecomposition(p,activeX,activeD,ja);
! 620: #endif
! 621: hja = newArrayOfPOLYold(p);
! 622: for (i=0; i<p; i++) {
! 623: getArrayOfPOLYold(hja,i) = hilbert1T(ja[i+1]);
! 624: }
! 625:
! 626: a = polyToInt(getArrayOfPOLYold(ja[0],0));
! 627: ansp = ppSub(iToPoly(1),xx(0,-a));
! 628: f = ppMult(getArrayOfPOLYold(hja,0),xx(0,-a));
! 629: ansp = ppAdd(ansp,f);
! 630: for (i=1; i<p; i++) {
! 631: f = ppSub(getArrayOfPOLYold(hja,i),getArrayOfPOLYold(hja,i-1));
! 632: a = polyToInt(getArrayOfPOLYold(ja[0],i));
! 633: f = ppMult(f,xx(0,-a));
! 634: ansp = ppAdd(ansp,f);
! 635: }
! 636:
! 637: #ifdef DEBUG
! 638: printf("<<%d>>",--cc); pWriteln(ansp);/***/
! 639: #endif
! 640:
! 641: return(ansp);
! 642: }
! 643:
! 644:
! 645: POLY hilbert2(k,p,pArray)
! 646: POLY k;
! 647: int p;
! 648: POLY pArray[];
! 649: /* This returns n! H(k,p,a^0, ..., a^{p-1}) */
! 650: /* n = (NN-M); */
! 651: /* Expample: hilbert2(xx(0,1),3,...) */
! 652: {
! 653: struct arrayOfPOLYold set;
! 654: int i,j;
! 655: POLY f;
! 656: POLY ans;
! 657: POLY c;
! 658: POLY g;
! 659: int n;
! 660: int size,a;
! 661:
! 662: set = newArrayOfPOLYold(p);
! 663: for (i=0; i<p; i++) {
! 664: if (pArray[i] ISZERO) {
! 665: warningHilbert("Argument of hilbert1 contains 0.\n");
! 666: return(POLYNULL);
! 667: }
! 668: f = pcmCopy(pArray[i]); /* pArray is already monomial poly */
! 669: /* This is for hilbert2 */
! 670: for (j = 0; j < MM; j++) {
! 671: (f->m->e)[j].x = 0; (f->m->e)[j].D = 0;
! 672: }
! 673: for (j = M; j < NN; j++) {
! 674: (f->m->e)[j].x = 0;
! 675: }
! 676: getArrayOfPOLYold(set,i) = f;
! 677: }
! 678: set = reduceMonomials(set);
! 679: f = hilbert1T(set); /* c x0^a ---> c {{k+n+a} \choose n} */
! 680:
! 681: n = NN-M;
! 682: ans = POLYNULL;
! 683: while (!(f ISZERO)) {
! 684: a = (f->m->e)[0].x;
! 685: c = newCell(f->coeffp,f->m); c = pcmCopy(c);
! 686: (c->m->e)[0].x = 0;
! 687: g = ppifac(ppAdd(k,iToPoly(a)),n);
! 688: ans = ppAdd(ans, ppMult(c,g));
! 689: f = f->next;
! 690: }
! 691: return(ans);
! 692:
! 693: }
! 694:
! 695:
! 696:
! 697:
! 698: #ifdef DEBUG2
! 699: checkh(set,i)
! 700: struct arrayOfPOLYold set;
! 701: int i;
! 702: {
! 703: if (pLength(getArrayOfPOLYold(set,i)) != 1) {
! 704: printf("Size is %d.",pSize(getArrayOfPOLYold(set,i)));
! 705: getchar(); getchar(); getchar(); getchar();
! 706: getchar();getchar();
! 707: }
! 708: }
! 709: #endif
! 710:
! 711: #ifdef DEBUG3
! 712: outputDecomposition(p,activeX,activeD,ja)
! 713: int p;
! 714: int activeX[];
! 715: int activeD[];
! 716: struct arrayOfPOLYold ja[];
! 717: {
! 718: int i;
! 719: printf("\nActiveX: ");
! 720: for (i=0; i<N; i++) {
! 721: printf("%3d",activeX[i]);
! 722: }
! 723: printf("\nActiveD: ");
! 724: for (i=0; i<N; i++) {
! 725: printf("%3d",activeD[i]);
! 726: }
! 727: printf(" Decomposed into J_0 to J_%d\n",p-1);
! 728: if (1) {
! 729: printf("Exponents: ");
! 730: for (i=0; i<p; i++) {
! 731: printf("%s, ",POLYToString((ja[0]).array[i],' ',0));
! 732: }
! 733: printf("\n");
! 734: for (i=0; i<p; i++) {
! 735: outputarrayOfPOLYold(ja[i+1]);
! 736: printf("\n");
! 737: }
! 738: }
! 739: }
! 740:
! 741: outputarrayOfPOLYold(set)
! 742: struct arrayOfPOLYold set;
! 743: {
! 744: int i;
! 745: for (i=0; i< set.n ; i++) {
! 746: printf("%s, ",POLYToString(set.array[i],' ',0));
! 747: }
! 748: }
! 749: #endif
! 750:
! 751:
! 752: warningHilbert(str)
! 753: char str[];
! 754: {
! 755: fprintf(stderr,"Warning (hilbert.c): %s\n",str);
! 756: }
! 757:
! 758: errorHilbert(str)
! 759: char str[];
! 760: {
! 761: errorKan1("%s\n",str);
! 762: }
! 763:
! 764:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to hilbert.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / kan96xx / Kan