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