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

Diff for /OpenXM/src/kan96xx/Kan/hilbert.c between version 1.2 and 1.4

version 1.2, 2000/01/16 07:55:39 version 1.4, 2005/06/16 05:07:23
Line 1 
Line 1 
 /* $OpenXM$ */  /* $OpenXM: OpenXM/src/kan96xx/Kan/hilbert.c,v 1.3 2001/05/04 01:06:23 takayama Exp $ */
 /*   hilbert.c  /*   hilbert.c
     1992/06/16       1992/06/16
     1992/06/18       1992/06/18
     1993/04/28       1993/04/28
         \bibitem{Cocoa-Hilbert}       \bibitem{Cocoa-Hilbert}
          Bigatti, A., Caboara, M., Robianno, L. (1991):       Bigatti, A., Caboara, M., Robianno, L. (1991):
          On the computation of Hilbert Poincar\'e series.       On the computation of Hilbert Poincar\'e series.
          Applicable algebra in Engineering, Communication and Computing       Applicable algebra in Engineering, Communication and Computing
         {\bf 2}, 21--33.       {\bf 2}, 21--33.
     1998/10/31       1998/10/31
 */  */
 #include <stdio.h>  #include <stdio.h>
 #include "datatype.h"  #include "datatype.h"
Line 41  struct object hilberto(struct object obgb,struct objec
Line 41  struct object hilberto(struct object obgb,struct objec
   int n; /* number of variables */    int n; /* number of variables */
   int i,j,k;    int i,j,k;
   int n0; /* number of the variables for the polynomials in base. */    int n0; /* number of the variables for the polynomials in base. */
   struct object obf;    struct object obf = OINIT;
   struct ring *rp;    struct ring *rp;
   struct ring *rr = (struct ring *)NULL;    struct ring *rr = (struct ring *)NULL;
   POLY *base;    POLY *base;
Line 50  struct object hilberto(struct object obgb,struct objec
Line 50  struct object hilberto(struct object obgb,struct objec
   int debug = 0;    int debug = 0;
   POLY f;    POLY f;
   POLY g;    POLY g;
   struct object rob;    struct object rob = OINIT;
   struct object ccc;    struct object ccc = OINIT;
   extern struct ring SmallRing;    extern struct ring SmallRing;
   int worg;    int worg;
   extern int WarningNoVectorVariable;    extern int WarningNoVectorVariable;
Line 79  struct object hilberto(struct object obgb,struct objec
Line 79  struct object hilberto(struct object obgb,struct objec
     f = KopPOLY(obf);      f = KopPOLY(obf);
     if (rr != (struct ring*)NULL) {      if (rr != (struct ring*)NULL) {
       if (rr != f->m->ringp) {        if (rr != f->m->ringp) {
         errorKan1("%s\n","Hhilbert(): the element must belong to the same ring.");          errorKan1("%s\n","Hhilbert(): the element must belong to the same ring.");
       }        }
     }else{      }else{
       rr =  f->m->ringp;        rr =  f->m->ringp;
Line 96  struct object hilberto(struct object obgb,struct objec
Line 96  struct object hilberto(struct object obgb,struct objec
     f = KopPOLY(obf);      f = KopPOLY(obf);
     if (rr != (struct ring*)NULL) {      if (rr != (struct ring*)NULL) {
       if (rr != f->m->ringp) {        if (rr != f->m->ringp) {
         errorKan1("%s\n","Hhilbert(): the element must belong to the same ring.");          errorKan1("%s\n","Hhilbert(): the element must belong to the same ring.");
       }        }
     }else{      }else{
       rr =  f->m->ringp;        rr =  f->m->ringp;
Line 126  struct object hilberto(struct object obgb,struct objec
Line 126  struct object hilberto(struct object obgb,struct objec
     f = KopPOLY(obf);      f = KopPOLY(obf);
     for (j=0; j<n0; j++) {      for (j=0; j<n0; j++) {
       if ((f->m->e)[j].x) {        if ((f->m->e)[j].x) {
         x[j] = 1; break;          x[j] = 1; break;
       }        }
       if ((f->m->e)[j].D) {        if ((f->m->e)[j].D) {
         d[j] = 1; break;          d[j] = 1; break;
       }        }
     }      }
   }    }
Line 158  struct object hilberto(struct object obgb,struct objec
Line 158  struct object hilberto(struct object obgb,struct objec
     g = cxx(1,0,0,rp);      g = cxx(1,0,0,rp);
     for (j=0; j<n0; j++) {      for (j=0; j<n0; j++) {
       if ((f->m->e)[j].x) {        if ((f->m->e)[j].x) {
         if (x[j] != -1) {          if (x[j] != -1) {
           (g->m->e)[x[j]].D = (f->m->e)[j].x;            (g->m->e)[x[j]].D = (f->m->e)[j].x;
         }          }
       }        }
       if ((f->m->e)[j].D) {        if ((f->m->e)[j].D) {
         if (d[j] != -1) {          if (d[j] != -1) {
           (g->m->e)[d[j]].D = (f->m->e)[j].D;            (g->m->e)[d[j]].D = (f->m->e)[j].D;
         }          }
       }        }
     }      }
     base[i] = g;      base[i] = g;
Line 220  static int polyToInt(POLY f) {
Line 220  static int polyToInt(POLY f) {
   
   
 static shell(v,n)  static shell(v,n)
 int v[];       int v[];
 int n;       int n;
 {  {
   int gap,i,j,temp;    int gap,i,j,temp;
   
   for (gap = n/2; gap > 0; gap /= 2) {    for (gap = n/2; gap > 0; gap /= 2) {
     for (i = gap; i<n; i++) {      for (i = gap; i<n; i++) {
       for (j=i-gap ; j>=0 && v[j]>v[j+gap]; j -= gap) {        for (j=i-gap ; j>=0 && v[j]>v[j+gap]; j -= gap) {
         temp = v[j];          temp = v[j];
         v[j] = v[j+gap];          v[j] = v[j+gap];
         v[j+gap] = temp;          v[j+gap] = temp;
       }        }
     }      }
   }    }
 }  }
   
 static POLY ppifac(f,n)  static POLY ppifac(f,n)
 POLY f;       POLY f;
 int n;       int n;
 /*  ppifac returns (f+n) (f+n-1) ... (f+1). */       /*  ppifac returns (f+n) (f+n-1) ... (f+1). */
 {  {
   POLY r;    POLY r;
   int i;    int i;
Line 259  int n;
Line 259  int n;
   
   
 static int isReducibleD1(f,g)  static int isReducibleD1(f,g)
 POLY f;       POLY f;
 POLY g;       POLY g;
 {  {
   int i;    int i;
   for (i = M; i < NN; i++) {    for (i = M; i < NN; i++) {
     if (((f->m->e)[i].x >= (g->m->e)[i].x) &&      if (((f->m->e)[i].x >= (g->m->e)[i].x) &&
         ((f->m->e)[i].D >= (g->m->e)[i].D)) {          ((f->m->e)[i].D >= (g->m->e)[i].D)) {
       /* do nothing */        /* do nothing */
     } else {      } else {
       return(0);        return(0);
Line 277  POLY g;
Line 277  POLY g;
   
 static struct arrayOfPOLYold  static struct arrayOfPOLYold
 reduceMonomials(set)  reduceMonomials(set)
 struct arrayOfPOLYold set;       struct arrayOfPOLYold set;
 {  {
   int *drop;  /* 1--> yes. drop set[i]. 0--> no. */    int *drop;  /* 1--> yes. drop set[i]. 0--> no. */
   int i,j;    int i,j;
Line 295  struct arrayOfPOLYold set;
Line 295  struct arrayOfPOLYold set;
   for (i = 0; i < size; i++) {    for (i = 0; i < size; i++) {
     if (!drop[i]) {      if (!drop[i]) {
       for (j=i+1; j < size; j++) {        for (j=i+1; j < size; j++) {
         if (!drop[j]) {          if (!drop[j]) {
           if (isReducibleD1(getArrayOfPOLYold(set,i),            if (isReducibleD1(getArrayOfPOLYold(set,i),
                             getArrayOfPOLYold(set,j))) {                              getArrayOfPOLYold(set,j))) {
             drop[i] = 1; break;              drop[i] = 1; break;
           }else if (isReducibleD1(getArrayOfPOLYold(set,j),            }else if (isReducibleD1(getArrayOfPOLYold(set,j),
                             getArrayOfPOLYold(set,i))) {                                    getArrayOfPOLYold(set,i))) {
             drop[j] = 1; break;              drop[j] = 1; break;
           }else { }            }else { }
         }          }
       }        }
     }      }
   }    }
Line 326  struct arrayOfPOLYold set;
Line 326  struct arrayOfPOLYold set;
   return(ans);    return(ans);
 }  }
   
   
 static int tryDecompose(set,i,j,vA,vWhich)  static int tryDecompose(set,i,j,vA,vWhich)
 struct arrayOfPOLYold set;    /* input: monomials */       struct arrayOfPOLYold set;    /* input: monomials */
 int     i,j;               /* decompose with respect to the (i,j)-th       int     i,j;               /* decompose with respect to the (i,j)-th
                               variable: i=0-->x_j, i=1--->D_j */                                     variable: i=0-->x_j, i=1--->D_j */
 int vA[];                  /* Return value: vA[0] is a_0, ... */       int vA[];                  /* Return value: vA[0] is a_0, ... */
 int vWhich[];              /* Return value: vWhich[i] is <<a>> of set[i] */       int vWhich[];              /* Return value: vWhich[i] is <<a>> of set[i] */
 /* set ---> x_j^(a_0) I_{a_0} + .... + x_j^{a_{p-1}} I_{a_{p-1}} */       /* set ---> x_j^(a_0) I_{a_0} + .... + x_j^{a_{p-1}} I_{a_{p-1}} */
 /* return value is p */       /* return value is p */
 /* tryDecompose is used to find the best decomposition. */       /* tryDecompose is used to find the best decomposition. */
 {  {
   int size,k,ell,ell2;    int size,k,ell,ell2;
   POLY f;    POLY f;
Line 356  int vWhich[];              /* Return value: vWhich[i] 
Line 356  int vWhich[];              /* Return value: vWhich[i] 
   }    }
 #ifdef DEBUG3  #ifdef DEBUG3
   /*for (i=0; i<size; i++) printf("%3d", vWhich[i]);    /*for (i=0; i<size; i++) printf("%3d", vWhich[i]);
   printf("     vWhich\n");*/      printf("     vWhich\n");*/
 #endif  #endif
   /* sort and prune */    /* sort and prune */
   tmp = (int *)sGC_malloc(sizeof(int)*((size+1)*2));    tmp = (int *)sGC_malloc(sizeof(int)*((size+1)*2));
Line 371  int vWhich[];              /* Return value: vWhich[i] 
Line 371  int vWhich[];              /* Return value: vWhich[i] 
   
 #ifdef DEBUG3  #ifdef DEBUG3
   /*for (i=0; i<p; i++) printf("%3d", vA[i]);    /*for (i=0; i<p; i++) printf("%3d", vA[i]);
   printf("---vA\n");*/      printf("---vA\n");*/
 #endif  #endif
   return(p);    return(p);
 }  }
   
 static struct arrayOfPOLYold getJA(a,set,vWhich,ja,ii,xd,ith)  static struct arrayOfPOLYold getJA(a,set,vWhich,ja,ii,xd,ith)
 /* get J_{a_{i+1}} from J_{a_i}       /* get J_{a_{i+1}} from J_{a_i}
    J_{a_{i+1}} = J_{a_i} \cup I_{a_{i+1}}     J_{a_{i+1}} = J_{a_i} \cup I_{a_{i+1}}
 */  */
 int a;                  /* each a_i */       int a;                  /* each a_i */
 struct arrayOfPOLYold set; /* polynomials */  struct arrayOfPOLYold set; /* polynomials */
 int *vWhich;            /* vWhich[0] is exponent a of set[0], .... */  int *vWhich;            /* vWhich[0] is exponent a of set[0], .... */
 struct arrayOfPOLYold ja;  /* J_i */  struct arrayOfPOLYold ja;  /* J_i */
Line 400  int ith;                /* x_{ith} or D_{ith} is the b
Line 400  int ith;                /* x_{ith} or D_{ith} is the b
   
   size = set.n;    size = set.n;
   start = 0;    start = 0;
 /*  input = newArrayOfPOLYold(size); */    /*  input = newArrayOfPOLYold(size); */
   input = newArrayOfPOLYold(size+ja.n+1);  /* 1993/09/08 */    input = newArrayOfPOLYold(size+ja.n+1);  /* 1993/09/08 */
   
   
Line 414  int ith;                /* x_{ith} or D_{ith} is the b
Line 414  int ith;                /* x_{ith} or D_{ith} is the b
     if (vWhich[i] == a) {      if (vWhich[i] == a) {
       f = pcmCopy(getArrayOfPOLYold(set,i));        f = pcmCopy(getArrayOfPOLYold(set,i));
       if (xd == 0) {        if (xd == 0) {
         (f->m->e)[ith].x = 0;          (f->m->e)[ith].x = 0;
       }else{        }else{
         (f->m->e)[ith].D = 0;          (f->m->e)[ith].D = 0;
       }        }
       getArrayOfPOLYold(input,start) = f;        getArrayOfPOLYold(input,start) = f;
       start++ ;        start++ ;
Line 435  int ith;                /* x_{ith} or D_{ith} is the b
Line 435  int ith;                /* x_{ith} or D_{ith} is the b
 #endif  #endif
   return( input );    return( input );
 }  }
   
   
 static struct arrayOfPOLYold *getDecomposition(set,activeX,activeD)  static struct arrayOfPOLYold *getDecomposition(set,activeX,activeD)
 struct arrayOfPOLYold set;       struct arrayOfPOLYold set;
 int activeX[N0];       int activeX[N0];
 int activeD[N0];       int activeD[N0];
 {  {
   int i;    int i;
   int size;    int size;
Line 462  int activeD[N0];
Line 462  int activeD[N0];
     if (activeX[i]) {      if (activeX[i]) {
       p = tryDecompose(set,0,i,vA,vWhich);        p = tryDecompose(set,0,i,vA,vWhich);
       if (p > xmax) {        if (p > xmax) {
         xmax = p;          xmax = p;
         xi = i;          xi = i;
         if (xmax == size) {          if (xmax == size) {
           break;            break;
         }          }
       }        }
     }      }
   }    }
   if (xmax != size) {    if (xmax != size) {
     for ( i = M; i < NN; i++) {      for ( i = M; i < NN; i++) {
       if (activeD[i]) {        if (activeD[i]) {
         p = tryDecompose(set,1,i,vA,vWhich);          p = tryDecompose(set,1,i,vA,vWhich);
         if (p > dmax) {          if (p > dmax) {
           dmax = p;            dmax = p;
           di = i;            di = i;
           if (dmax == size) {            if (dmax == size) {
             break;              break;
           }            }
         }          }
       }        }
     }      }
   }    }
   /*    /*
     ans[0] = (a_0,...,a_{p-1})      ans[0] = (a_0,...,a_{p-1})
     ans[1] = generators of J_{a_0}      ans[1] = generators of J_{a_0}
      ....      ....
     ans[p] = generators of J_{a_{p-1}},   p = xmax or dmax      ans[p] = generators of J_{a_{p-1}},   p = xmax or dmax
   */    */
   if (xmax > dmax) {    if (xmax > dmax) {
Line 520  int activeD[N0];
Line 520  int activeD[N0];
 }  }
   
 static POLY hilbert1T(set)  static POLY hilbert1T(set)
 struct arrayOfPOLYold set;       struct arrayOfPOLYold set;
 /* <<set>> must be reduced basis and each polynomial must have the length 1 */       /* <<set>> must be reduced basis and each polynomial must have the length 1 */
 /* Unnecessary exponents should be set to zero. For example, f = x_{M-1} x_M       /* Unnecessary exponents should be set to zero. For example, f = x_{M-1} x_M
    is illegal input. It should be x_M ( M <= index <= NN ).     is illegal input. It should be x_M ( M <= index <= NN ).
    cf. hilbert1() and hilbert2()     cf. hilbert1() and hilbert2()
 */  */
Line 566  struct arrayOfPOLYold set;
Line 566  struct arrayOfPOLYold set;
     f = getArrayOfPOLYold(set,i);      f = getArrayOfPOLYold(set,i);
     for (j = M; j < NN; j++) {      for (j = M; j < NN; j++) {
       if ((f->m->e)[j].x) {        if ((f->m->e)[j].x) {
         activeX[j] = 1;          activeX[j] = 1;
       }        }
       if ((f->m->e)[j].D) {        if ((f->m->e)[j].D) {
         activeD[j] = 1;          activeD[j] = 1;
       }        }
     }      }
   }    }
   for (i = M; i < NN; i++) {    for (i = M; i < NN; i++) {
     if (activeX[i]) {      if (activeX[i]) {
       if (active) {        if (active) {
         active = 2;          active = 2;
         break;          break;
       }        }
       active = 1;        active = 1;
       xxdd = 0; xi = i;        xxdd = 0; xi = i;
     }      }
     if (activeD[i]) {      if (activeD[i]) {
       if (active) {        if (active) {
         active = 2;          active = 2;
         break;          break;
       }        }
       active = 1;        active = 1;
       xxdd = 1; di = i;        xxdd = 1; di = i;
Line 642  struct arrayOfPOLYold set;
Line 642  struct arrayOfPOLYold set;
   return(ansp);    return(ansp);
 }  }
   
   
 POLY hilbert2(k,p,pArray)  POLY hilbert2(k,p,pArray)
 POLY k;       POLY k;
 int p;       int p;
 POLY pArray[];       POLY pArray[];
 /* This returns n! H(k,p,a^0, ..., a^{p-1}) */       /* This returns n! H(k,p,a^0, ..., a^{p-1}) */
 /* n = (NN-M); */       /* n = (NN-M); */
 /* Expample: hilbert2(xx(0,1),3,...) */       /* Expample: hilbert2(xx(0,1),3,...) */
 {  {
   struct arrayOfPOLYold set;    struct arrayOfPOLYold set;
   int i,j;    int i,j;
Line 698  POLY pArray[];
Line 698  POLY pArray[];
   
 #ifdef DEBUG2  #ifdef DEBUG2
 checkh(set,i)  checkh(set,i)
 struct arrayOfPOLYold set;       struct arrayOfPOLYold set;
 int i;       int i;
 {  {
   if (pLength(getArrayOfPOLYold(set,i)) != 1) {    if (pLength(getArrayOfPOLYold(set,i)) != 1) {
       printf("Size is %d.",pSize(getArrayOfPOLYold(set,i)));      printf("Size is %d.",pSize(getArrayOfPOLYold(set,i)));
       getchar(); getchar(); getchar(); getchar();      getchar(); getchar(); getchar(); getchar();
       getchar();getchar();      getchar();getchar();
     }    }
 }  }
 #endif  #endif
   
 #ifdef DEBUG3  #ifdef DEBUG3
 outputDecomposition(p,activeX,activeD,ja)  outputDecomposition(p,activeX,activeD,ja)
 int p;       int p;
 int activeX[];       int activeX[];
 int activeD[];       int activeD[];
 struct arrayOfPOLYold ja[];       struct arrayOfPOLYold ja[];
 {  {
   int i;    int i;
   printf("\nActiveX: ");    printf("\nActiveX: ");
Line 740  struct arrayOfPOLYold ja[];
Line 740  struct arrayOfPOLYold ja[];
 }  }
   
 outputarrayOfPOLYold(set)  outputarrayOfPOLYold(set)
 struct arrayOfPOLYold set;       struct arrayOfPOLYold set;
 {  {
   int i;    int i;
   for (i=0; i< set.n ; i++) {    for (i=0; i< set.n ; i++) {
Line 751  struct arrayOfPOLYold set;
Line 751  struct arrayOfPOLYold set;
   
   
 warningHilbert(str)  warningHilbert(str)
 char str[];       char str[];
 {  {
   fprintf(stderr,"Warning (hilbert.c): %s\n",str);    fprintf(stderr,"Warning (hilbert.c): %s\n",str);
 }  }
   
 errorHilbert(str)  errorHilbert(str)
 char str[];       char str[];
 {  {
   errorKan1("%s\n",str);    errorKan1("%s\n",str);
 }  }
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.4
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>