[BACK]Return to hilbert.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / Kan

Annotation of OpenXM/src/kan96xx/Kan/hilbert.c, Revision 1.2

1.2     ! takayama    1: /* $OpenXM$ */
1.1       maekawa     2: /*   hilbert.c
                      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
                     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) {
                     82:        errorKan1("%s\n","Hhilbert(): the element must belong to the same ring.");
                     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) {
                     99:        errorKan1("%s\n","Hhilbert(): the element must belong to the same ring.");
                    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) {
                    129:        x[j] = 1; break;
                    130:       }
                    131:       if ((f->m->e)[j].D) {
                    132:        d[j] = 1; break;
                    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) {
                    161:        if (x[j] != -1) {
                    162:          (g->m->e)[x[j]].D = (f->m->e)[j].x;
                    163:        }
                    164:       }
                    165:       if ((f->m->e)[j].D) {
                    166:        if (d[j] != -1) {
                    167:          (g->m->e)[d[j]].D = (f->m->e)[j].D;
                    168:        }
                    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)
                    223: int v[];
                    224: int n;
                    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) {
                    231:        temp = v[j];
                    232:        v[j] = v[j+gap];
                    233:        v[j+gap] = temp;
                    234:       }
                    235:     }
                    236:   }
                    237: }
                    238:
                    239: static POLY ppifac(f,n)
                    240: POLY f;
                    241: int n;
                    242: /*  ppifac returns (f+n) (f+n-1) ... (f+1). */
                    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)
                    262: POLY f;
                    263: POLY g;
                    264: {
                    265:   int i;
                    266:   for (i = M; i < NN; i++) {
                    267:     if (((f->m->e)[i].x >= (g->m->e)[i].x) &&
                    268:        ((f->m->e)[i].D >= (g->m->e)[i].D)) {
                    269:       /* do nothing */
                    270:     } else {
                    271:       return(0);
                    272:     }
                    273:   }
                    274:   return(1);
                    275: }
                    276:
                    277:
                    278: static struct arrayOfPOLYold
                    279: reduceMonomials(set)
                    280: struct arrayOfPOLYold set;
                    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++) {
                    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:        }
                    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:
                    329:
                    330: static int tryDecompose(set,i,j,vA,vWhich)
                    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. */
                    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]);
                    359:   printf("     vWhich\n");*/
                    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]);
                    374:   printf("---vA\n");*/
                    375: #endif
                    376:   return(p);
                    377: }
                    378:
                    379: static struct arrayOfPOLYold getJA(a,set,vWhich,ja,ii,xd,ith)
                    380: /* get J_{a_{i+1}} from J_{a_i}
                    381:    J_{a_{i+1}} = J_{a_i} \cup I_{a_{i+1}}
                    382: */
                    383: int a;                  /* each a_i */
                    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;
                    403: /*  input = newArrayOfPOLYold(size); */
                    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) {
                    417:        (f->m->e)[ith].x = 0;
                    418:       }else{
                    419:        (f->m->e)[ith].D = 0;
                    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: }
                    438:
                    439:
                    440: static struct arrayOfPOLYold *getDecomposition(set,activeX,activeD)
                    441: struct arrayOfPOLYold set;
                    442: int activeX[N0];
                    443: int activeD[N0];
                    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) {
                    465:        xmax = p;
                    466:        xi = i;
                    467:        if (xmax == size) {
                    468:          break;
                    469:        }
                    470:       }
                    471:     }
                    472:   }
                    473:   if (xmax != size) {
                    474:     for ( i = M; i < NN; i++) {
                    475:       if (activeD[i]) {
                    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:        }
                    484:       }
                    485:     }
                    486:   }
                    487:   /*
                    488:     ans[0] = (a_0,...,a_{p-1})
                    489:     ans[1] = generators of J_{a_0}
                    490:      ....
                    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)
                    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
                    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) {
                    569:        activeX[j] = 1;
                    570:       }
                    571:       if ((f->m->e)[j].D) {
                    572:        activeD[j] = 1;
                    573:       }
                    574:     }
                    575:   }
                    576:   for (i = M; i < NN; i++) {
                    577:     if (activeX[i]) {
                    578:       if (active) {
                    579:        active = 2;
                    580:        break;
                    581:       }
                    582:       active = 1;
                    583:       xxdd = 0; xi = i;
                    584:     }
                    585:     if (activeD[i]) {
                    586:       if (active) {
                    587:        active = 2;
                    588:        break;
                    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:
                    645:
                    646: POLY hilbert2(k,p,pArray)
                    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,...) */
                    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)
                    701: struct arrayOfPOLYold set;
                    702: int i;
                    703: {
                    704:   if (pLength(getArrayOfPOLYold(set,i)) != 1) {
                    705:       printf("Size is %d.",pSize(getArrayOfPOLYold(set,i)));
                    706:       getchar(); getchar(); getchar(); getchar();
                    707:       getchar();getchar();
                    708:     }
                    709: }
                    710: #endif
                    711:
                    712: #ifdef DEBUG3
                    713: outputDecomposition(p,activeX,activeD,ja)
                    714: int p;
                    715: int activeX[];
                    716: int activeD[];
                    717: struct arrayOfPOLYold ja[];
                    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)
                    743: struct arrayOfPOLYold set;
                    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)
                    754: char str[];
                    755: {
                    756:   fprintf(stderr,"Warning (hilbert.c): %s\n",str);
                    757: }
                    758:
                    759: errorHilbert(str)
                    760: char str[];
                    761: {
                    762:   errorKan1("%s\n",str);
                    763: }
                    764:
                    765:

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>