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

Annotation of OpenXM/src/kan96xx/Kan/ecart.c, Revision 1.1

1.1     ! takayama    1: /* $OpenXM$ */
        !             2: #include <stdio.h>
        !             3: #include "datatype.h"
        !             4: #include "extern2.h"
        !             5: #include "gradedset.h"
        !             6:
        !             7: #define outofmemory(a) if (a == NULL) {fprintf(stderr,"ecart.c: out of memory.\n"); exit(10);}
        !             8: /* Local structures */
        !             9: struct ecartReducer {
        !            10:   int ell; /* s^\ell */
        !            11:   int first; /* first =1 ==> gset else gg */
        !            12:   int grade;  int gseti;  /* if first==1,  gradedPolySet */
        !            13:   int ggi;                /* if first==0   ecartPoly [] */
        !            14: };
        !            15: struct ecartPolyArray {
        !            16:   int size;
        !            17:   int limit;
        !            18:   POLY *pa;
        !            19: };
        !            20:
        !            21: static struct ecartReducer ecartFindReducer(POLY r,struct gradedPolySet *gset,struct ecartPolyArray *epa);
        !            22: static int ecartCheckPoly(POLY f);  /* check if it does not contain s-variables */
        !            23: static int ecartCheckEnv();         /* check if the environment is OK for ecart div*/
        !            24: static struct ecartPolyArray *ecartPutPolyInG(POLY g,struct ecartPolyArray *eparray);
        !            25: static int ecartGetEll(POLY r,POLY g);
        !            26: static POLY ecartDivideSv(POLY r);
        !            27:
        !            28: extern int DebugReductionRed;
        !            29:
        !            30: /* This is used for goHomogenization */
        !            31: extern int DegreeShifto_size;
        !            32: extern int *DegreeShifto_vec;
        !            33:
        !            34: #define LARGE 0x7fffffff
        !            35:
        !            36:
        !            37: static POLY ecartDivideSv(POLY r) {
        !            38:   POLY f;
        !            39:   int k;
        !            40:   if (r == POLYNULL) return r;
        !            41:   f = r; k = LARGE;
        !            42:   while (r != POLYNULL) {
        !            43:     if (r->m->e[0].x < k) {
        !            44:       k = r->m->e[0].x;
        !            45:     }
        !            46:     r = r->next;
        !            47:   }
        !            48:   if (k > 0) {
        !            49:     f = r;
        !            50:     while (r != POLYNULL) {
        !            51:       r->m->e[0].x -= k;
        !            52:       r = r->next;
        !            53:     }
        !            54:   }
        !            55:   return f;
        !            56: }
        !            57:
        !            58: static int ecartGetEll(POLY f,POLY g) {
        !            59:   int n,i,p;
        !            60:   MONOMIAL tf;
        !            61:   MONOMIAL tg;
        !            62:
        !            63:   if (f ISZERO) return(-1);
        !            64:   if (g ISZERO) return(-1);
        !            65:
        !            66:   checkRingIsR(f,g);
        !            67:
        !            68:   if (!isSameComponent_x(f,g)) return(-1);
        !            69:   tf = f->m; tg = g->m; n = tf->ringp->n;
        !            70:   for (i=1; i<n; i++) {
        !            71:     if (tf->e[i].x < tg->e[i].x) return(-1);
        !            72:     if (tf->e[i].D < tg->e[i].D) return(-1);
        !            73:   }
        !            74:   if (tf->e[0].D < tg->e[0].D) return(-1);  /*  h  */
        !            75:   p = tf->e[i].x - tg->e[i].x;  /* H,  s */
        !            76:   if (p >=0 ) return 0;
        !            77:   else return(-p);
        !            78: }
        !            79:
        !            80:
        !            81: #define EP_SIZE 10
        !            82: static struct ecartPolyArray *ecartPutPolyInG(POLY g,struct ecartPolyArray *eparray)
        !            83: {
        !            84:   struct ecartPolyArray *a;
        !            85:   POLY *ep;
        !            86:   int i;
        !            87:   if (eparray == (struct ecartPolyArray *)NULL) {
        !            88:     a = (struct ecartPolyArray *) sGC_malloc(sizeof(struct ecartPolyArray));
        !            89:     outofmemory(a);
        !            90:     ep = (POLY *) sGC_malloc(sizeof(POLY)*EP_SIZE);
        !            91:     outofmemory(ep);
        !            92:     a->limit = EP_SIZE;
        !            93:     a->size = 0;
        !            94:     a->pa = ep;
        !            95:     eparray = a;
        !            96:   }
        !            97:   if (eparray->size >= eparray->limit) {
        !            98:     a = (struct ecartPolyArray *) sGC_malloc(sizeof(struct ecartPolyArray));
        !            99:     outofmemory(a);
        !           100:     ep = (POLY *) sGC_malloc(sizeof(POLY)*((eparray->limit)*2));
        !           101:     outofmemory(ep);
        !           102:     a->limit = (eparray->limit)*2;
        !           103:     a->size = eparray->size;
        !           104:     a->pa = ep;
        !           105:     for (i=0; i<eparray->size; i++) {
        !           106:       (a->pa)[i] = (eparray->pa)[i];
        !           107:     }
        !           108:     eparray  = a;
        !           109:   }
        !           110:   i = eparray->size;
        !           111:   (eparray->pa)[i] = g;
        !           112:   (eparray->size)++;
        !           113:   return eparray;
        !           114: }
        !           115:
        !           116: static struct ecartReducer ecartFindReducer(POLY r,struct gradedPolySet *gset,
        !           117:                                             struct ecartPolyArray *epa)
        !           118: {
        !           119:   int grd;
        !           120:   struct polySet *set;
        !           121:   int minGrade = 0;
        !           122:   int minGseti = 0;
        !           123:   int minGgi   = 0;
        !           124:   int ell1 = LARGE;
        !           125:   int ell2 = LARGE;
        !           126:   int ell;
        !           127:   int i;
        !           128:   struct ecartReducer er;
        !           129:   /* Try to find a reducer in gset; */
        !           130:   grd = 0;
        !           131:   while (grd < gset->maxGrade) {
        !           132:     set = gset->polys[grd];
        !           133:     for (i=0; i<set->size; i++) {
        !           134:       ell = ecartGetEll(r,set->g[i]);
        !           135:       if ((ell>=0) && (ell < ell1)) {
        !           136:         ell1 = ell;
        !           137:         minGrade = grd; minGseti=i;
        !           138:       }
        !           139:       if (set->gh[i] == POLYNULL) {
        !           140:         /* goHomogenize set->gh[i] */
        !           141:         set->gh[i] = goHomogenize11(set->g[i],DegreeShifto_vec,DegreeShifto_size,-1,1);
        !           142:       }
        !           143:     }
        !           144:     grd++;
        !           145:   }
        !           146:   if (epa != NULL) {
        !           147:     /* Try to find in the second group. */
        !           148:     for (i=0; i< epa->size; i++) {
        !           149:       ell = ecartGetEll(r,(epa->pa)[i]);
        !           150:       if ((ell>=0) && (ell < ell2)) {
        !           151:         ell2 = ell;
        !           152:         minGgi = i;
        !           153:       }
        !           154:     }
        !           155:   }
        !           156:
        !           157:   if (DebugReductionRed) {
        !           158:     printf("ecartFindReducer(): ell1=%d, ell2=%d, minGrade=%d, minGseti=%d, minGgi=%d\n",ell1,ell2,minGrade,minGseti,minGgi);
        !           159:   }
        !           160:   if (ell1 <= ell2) {
        !           161:     if (ell1 == LARGE) {
        !           162:       er.ell = -1;
        !           163:       return er;
        !           164:     }else{
        !           165:       er.ell = ell1;
        !           166:       er.first = 1;
        !           167:       er.grade = minGrade;
        !           168:       er.gseti = minGseti;
        !           169:       return er;
        !           170:     }
        !           171:   }else{
        !           172:       er.ell = ell2;
        !           173:       er.first = 0;
        !           174:       er.ggi = minGgi;
        !           175:       return er;
        !           176:   }
        !           177: }
        !           178:
        !           179: /*
        !           180:   r and gset are assumed to be (0,1)-homogenized (h-homogenized)
        !           181:  */
        !           182: POLY reduction_ecart(r,gset,needSyz,syzp)
        !           183:      POLY r;
        !           184:      struct gradedPolySet *gset;
        !           185:      int needSyz;
        !           186:      struct syz0 *syzp; /* set */
        !           187: {
        !           188:   int reduced,reduced1,reduced2;
        !           189:   int grd;
        !           190:   struct polySet *set;
        !           191:   POLY cf,syz;
        !           192:   int i;
        !           193:   POLY cc,cg;
        !           194:   struct ecartReducer ells;
        !           195:   struct ecartPolyArray *gg;
        !           196:   POLY pp;
        !           197:   int ell;
        !           198:
        !           199:   extern struct ring *CurrentRingp;
        !           200:   struct ring *rp;
        !           201:
        !           202:   gg = NULL;
        !           203:   if (needSyz) {
        !           204:     if (r ISZERO) { rp = CurrentRingp; } else { rp = r->m->ringp; }
        !           205:     cf = cxx(1,0,0,rp);
        !           206:     syz = ZERO;
        !           207:   }
        !           208:
        !           209:   if (r != POLYNULL) {
        !           210:        if (! r->m->ringp->weightedHomogenization) {
        !           211:          errorKan1("%s\n","ecart.c: the given ring must be declared with [(weightedHomogenization) 1]");
        !           212:        }
        !           213:   }
        !           214:
        !           215:   r = goHomogenize11(r,DegreeShifto_vec,DegreeShifto_size,-1,1);
        !           216:   /* 1 means homogenize only s */
        !           217:
        !           218:   do {
        !           219:     if (DebugReductionRed) printf("r=%s\n",POLYToString(r,'*',1));
        !           220:     ells = ecartFindReducer(r,gset,gg);
        !           221:     ell = ells.ell;
        !           222:     if (ell > 0) {
        !           223:       gg = ecartPutPolyInG(r,gg);
        !           224:     }
        !           225:     if (ell >= 0) {
        !           226:       if (ells.first) {
        !           227:         pp = ((gset->polys[ells.grade])->gh)[ells.gseti];
        !           228:       }else{
        !           229:         pp = (gg->pa)[ells.ggi];
        !           230:       }
        !           231:       if (ell > 0) r = mpMult(cxx(1,0,ell,rp),r); /* r = s^ell r */
        !           232:       r = (*reduction1)(r,pp,needSyz,&cc,&cg);
        !           233:       if (needSyz) {
        !           234:         if (ells.first) {
        !           235:           cf = ppMult(cc,cf);
        !           236:           syz = cpMult(toSyzCoeff(cc),syz);
        !           237:           syz = ppAddv(syz,toSyzPoly(cg,ells.grade,ells.gseti));
        !           238:         }else{
        !           239:           /* BUG: not yet */
        !           240:         }
        !           241:       }
        !           242:       if (r ISZERO) goto ss;
        !           243:       r = ecartDivideSv(r); /* r = r/s^? */
        !           244:     }
        !           245:   }while (ell >= 0);
        !           246:
        !           247:  ss: ;
        !           248:   if (needSyz) {
        !           249:     syzp->cf = cf;   /* cf is in the CurrentRingp */
        !           250:     syzp->syz = syz; /* syz is in the SyzRingp */
        !           251:     /* BUG: dehomogenize the syzygy */
        !           252:   }
        !           253:   r = goDeHomogenizeS(r);
        !           254:   return(r);
        !           255: }

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