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