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

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

1.1     ! maekawa     1: #include <stdio.h>
        !             2: #include "datatype.h"
        !             3: #include "stackm.h"
        !             4: #include "extern.h"
        !             5: #include "extern2.h"
        !             6: #include "matrix.h"
        !             7: static void shell(int v[],int n);
        !             8: static int degreeOfPrincipalPart(POLY f);
        !             9: static int degreeOfInitW(POLY f,int w[]);
        !            10:
        !            11:
        !            12: static void shell(v,n)
        !            13: int v[];
        !            14: int n;
        !            15: {
        !            16:   int gap,i,j,temp;
        !            17:
        !            18:   for (gap = n/2; gap > 0; gap /= 2) {
        !            19:     for (i = gap; i<n; i++) {
        !            20:       for (j=i-gap ; j>=0 && v[j]<v[j+gap]; j -= gap) {
        !            21:        temp = v[j];
        !            22:        v[j] = v[j+gap];
        !            23:        v[j+gap] = temp;
        !            24:       }
        !            25:     }
        !            26:   }
        !            27: }
        !            28:
        !            29:
        !            30: struct matrixOfPOLY *parts(f,v)
        !            31: POLY f;
        !            32: POLY v;  /* v must be a single variable, e.g. x */
        !            33: {
        !            34:   struct matrixOfPOLY *evPoly;
        !            35:   int vi = 0;  /* index of v */
        !            36:   int vx = 1;  /* x --> 1, D--> 0 */
        !            37:   int n,evSize,i,k,e;
        !            38:   int *ev;
        !            39:   struct object *evList;
        !            40:   struct object *list;
        !            41:   struct object ob;
        !            42:   POLY ans;
        !            43:   POLY h;
        !            44:   extern struct ring *CurrentRingp;
        !            45:   POLY ft;
        !            46:
        !            47:
        !            48:   if (f ISZERO || v ISZERO) {
        !            49:     evPoly = newMatrixOfPOLY(2,1);
        !            50:     getMatrixOfPOLY(evPoly,0,0) = ZERO;
        !            51:     getMatrixOfPOLY(evPoly,1,0) = ZERO;
        !            52:     return(evPoly);
        !            53:   }
        !            54:   n = v->m->ringp->n;
        !            55:   /* get the index of the variable v */
        !            56:   for (i=0; i<n; i++) {
        !            57:     if (v->m->e[i].x) {
        !            58:       vx = 1; vi = i; break;
        !            59:     }else if (v->m->e[i].D) {
        !            60:       vx = 0; vi = i; break;
        !            61:     }
        !            62:   }
        !            63:   ft = f;
        !            64:   /* get the vector of exponents */
        !            65:   evList = NULLLIST;
        !            66:   while (ft != POLYNULL) {
        !            67:     if (vx) {
        !            68:       e = ft->m->e[vi].x;
        !            69:     }else{
        !            70:       e = ft->m->e[vi].D;
        !            71:     }
        !            72:     ft = ft->next;
        !            73:     ob = KpoInteger(e);
        !            74:     if (!memberQ(evList,ob)) {
        !            75:       list = newList(&ob);
        !            76:       evList = vJoin(evList,list);
        !            77:     }
        !            78:   }
        !            79:   /*printf("evList = "); printObjectList(evList);*/
        !            80:   evSize = klength(evList);
        !            81:   ev = (int *)sGC_malloc(sizeof(int)*(evSize+1));
        !            82:   if (ev == (int *)NULL) errorPoly("No more memory.");
        !            83:   for (i=0; i<evSize; i++) {
        !            84:     ev[i] = KopInteger(car(evList));
        !            85:     evList = cdr(evList);
        !            86:   }
        !            87:   /* sort ev */
        !            88:   shell(ev,evSize);
        !            89:
        !            90:   /* get coefficients */
        !            91:   evPoly = newMatrixOfPOLY(2,evSize);
        !            92:   for (i=0; i<evSize; i++) {
        !            93:     ans = ZERO;
        !            94:     getMatrixOfPOLY(evPoly,0,i) = cxx(ev[i],0,0,CurrentRingp);
        !            95:     ft = f;
        !            96:     while (ft != POLYNULL) {
        !            97:       if (vx) {
        !            98:        if (ft->m->e[vi].x == ev[i]) {
        !            99:          h = newCell(ft->coeffp,monomialCopy(ft->m));
        !           100:          xset0(h,vi); /* touch monomial part, so you need to copy it above. */
        !           101:          ans = ppAdd(ans,h);
        !           102:        }
        !           103:       }else{
        !           104:        if (ft->m->e[vi].D == ev[i]) {
        !           105:          h = newCell(ft->coeffp,monomialCopy(ft->m));
        !           106:          dset0(h,vi);
        !           107:          ans = ppAdd(ans,h);
        !           108:        }
        !           109:       }
        !           110:       ft = ft->next;
        !           111:     }
        !           112:     getMatrixOfPOLY(evPoly,1,i) = ans;
        !           113:   }
        !           114:   return(evPoly);
        !           115: }
        !           116:
        !           117: struct object parts2(f,v)
        !           118: POLY f;
        !           119: POLY v;  /* v must be a single variable, e.g. x */
        !           120: {
        !           121:   struct matrixOfPOLY *evPoly;
        !           122:   int vi = 0;  /* index of v */
        !           123:   int vx = 1;  /* x --> 1, D--> 0 */
        !           124:   int n,evSize,i,k,e;
        !           125:   int *ev;
        !           126:   struct object *evList;
        !           127:   struct object *list;
        !           128:   struct object ob;
        !           129:   POLY ans;
        !           130:   POLY h;
        !           131:   POLY ft;
        !           132:   struct object ob1,ob2,rob;
        !           133:
        !           134:
        !           135:   if (f ISZERO || v ISZERO) {
        !           136:     evPoly = newMatrixOfPOLY(2,1);
        !           137:     getMatrixOfPOLY(evPoly,0,0) = ZERO;
        !           138:     getMatrixOfPOLY(evPoly,1,0) = ZERO;
        !           139:     rob = newObjectArray(2);
        !           140:     ob1 = newObjectArray(1);
        !           141:     ob2 = newObjectArray(1);
        !           142:     putoa(ob1,0,KpoInteger(0));
        !           143:     putoa(ob2,0,KpoPOLY(POLYNULL));
        !           144:     putoa(rob,0,ob1); putoa(rob,1,ob2);
        !           145:     return(rob);
        !           146:   }
        !           147:   n = v->m->ringp->n;
        !           148:   /* get the index of the variable v */
        !           149:   for (i=0; i<n; i++) {
        !           150:     if (v->m->e[i].x) {
        !           151:       vx = 1; vi = i; break;
        !           152:     }else if (v->m->e[i].D) {
        !           153:       vx = 0; vi = i; break;
        !           154:     }
        !           155:   }
        !           156:   ft = f;
        !           157:   /* get the vector of exponents */
        !           158:   evList = NULLLIST;
        !           159:   while (ft != POLYNULL) {
        !           160:     if (vx) {
        !           161:       e = ft->m->e[vi].x;
        !           162:     }else{
        !           163:       e = ft->m->e[vi].D;
        !           164:     }
        !           165:     ft = ft->next;
        !           166:     ob = KpoInteger(e);
        !           167:     if (!memberQ(evList,ob)) {
        !           168:       list = newList(&ob);
        !           169:       evList = vJoin(evList,list);
        !           170:     }
        !           171:   }
        !           172:   /*printf("evList = "); printObjectList(evList);*/
        !           173:   evSize = klength(evList);
        !           174:   ev = (int *)sGC_malloc(sizeof(int)*(evSize+1));
        !           175:   if (ev == (int *)NULL) errorPoly("No more memory.");
        !           176:   for (i=0; i<evSize; i++) {
        !           177:     ev[i] = KopInteger(car(evList));
        !           178:     evList = cdr(evList);
        !           179:   }
        !           180:   /* sort ev */
        !           181:   shell(ev,evSize);
        !           182:
        !           183:   /* get coefficients */
        !           184:   evPoly = newMatrixOfPOLY(2,evSize);
        !           185:   for (i=0; i<evSize; i++) {
        !           186:     ans = ZERO;
        !           187:     /* getMatrixOfPOLY(evPoly,0,i) = cxx(ev[i],0,0,CurrentRingp); */
        !           188:     getMatrixOfPOLY(evPoly,0,i) = POLYNULL;
        !           189:     ft = f;
        !           190:     while (ft != POLYNULL) {
        !           191:       if (vx) {
        !           192:        if (ft->m->e[vi].x == ev[i]) {
        !           193:          h = newCell(ft->coeffp,monomialCopy(ft->m));
        !           194:          xset0(h,vi); /* touch monomial part, so you need to copy it above. */
        !           195:          ans = ppAdd(ans,h);
        !           196:        }
        !           197:       }else{
        !           198:        if (ft->m->e[vi].D == ev[i]) {
        !           199:          h = newCell(ft->coeffp,monomialCopy(ft->m));
        !           200:          dset0(h,vi);
        !           201:          ans = ppAdd(ans,h);
        !           202:        }
        !           203:       }
        !           204:       ft = ft->next;
        !           205:     }
        !           206:     getMatrixOfPOLY(evPoly,1,i) = ans;
        !           207:   }
        !           208:   rob = newObjectArray(2);
        !           209:   ob1 = newObjectArray(evSize);
        !           210:   ob2 = newObjectArray(evSize);
        !           211:   for (i=0; i<evSize; i++) {
        !           212:     putoa(ob2,i,KpoPOLY(getMatrixOfPOLY(evPoly,1,i)));
        !           213:     putoa(ob1,i,KpoInteger(ev[i]));
        !           214:   }
        !           215:   putoa(rob,0,ob1); putoa(rob,1,ob2);
        !           216:   return(rob);
        !           217: }
        !           218:
        !           219: int pDegreeWrtV(f,v)
        !           220: POLY f;
        !           221: POLY v;
        !           222: {
        !           223:   int vx = 1;
        !           224:   int vi = 0;
        !           225:   int i,n;
        !           226:   int ans;
        !           227:   if (f ISZERO || v ISZERO) return(0);
        !           228:   n = f->m->ringp->n;
        !           229:   for (i=0; i<n; i++) {
        !           230:     if (v->m->e[i].x) {
        !           231:       vx = 1; vi = i;
        !           232:       break;
        !           233:     }else if (v->m->e[i].D) {
        !           234:       vx = 0; vi = i;
        !           235:       break;
        !           236:     }
        !           237:   }
        !           238:   if (vx) {
        !           239:     ans = f->m->e[vi].x;
        !           240:   }else{
        !           241:     ans = f->m->e[vi].D;
        !           242:   }
        !           243:   f = f->next;
        !           244:   while (f != POLYNULL) {
        !           245:     if (vx) {
        !           246:       if (f->m->e[vi].x > ans) ans = f->m->e[vi].x;
        !           247:     }else{
        !           248:       if (f->m->e[vi].D > ans) ans = f->m->e[vi].D;
        !           249:     }
        !           250:     f = f->next;
        !           251:   }
        !           252:   return(ans);
        !           253: }
        !           254:
        !           255: int containVectorVariable(POLY f)
        !           256: {
        !           257:   MONOMIAL tf;
        !           258:   static int nn,mm,ll,cc,n,m,l,c;
        !           259:   static struct ring *cr = (struct ring *)NULL;
        !           260:   int i;
        !           261:
        !           262:   if (f ISZERO) return(0);
        !           263:   tf = f->m;
        !           264:   if (tf->ringp != cr) {
        !           265:     n = tf->ringp->n;
        !           266:     m = tf->ringp->m;
        !           267:     l = tf->ringp->l;
        !           268:     c = tf->ringp->c;
        !           269:     nn = tf->ringp->nn;
        !           270:     mm = tf->ringp->mm;
        !           271:     ll = tf->ringp->ll;
        !           272:     cc = tf->ringp->cc;
        !           273:     cr = tf->ringp;
        !           274:   }
        !           275:
        !           276:   while (f != POLYNULL) {
        !           277:     tf = f->m;
        !           278:     for (i=cc; i<c; i++) {
        !           279:       if ( tf->e[i].x ) return(1);
        !           280:       if ( tf->e[i].D ) return(1);
        !           281:     }
        !           282:     for (i=ll; i<l; i++) {
        !           283:       if (tf->e[i].x) return(1);
        !           284:       if (tf->e[i].D) return(1);
        !           285:     }
        !           286:     for (i=mm; i<m; i++) {
        !           287:       if (tf->e[i].x) return(1);
        !           288:       if (tf->e[i].D) return(1);
        !           289:     }
        !           290:     for (i=nn; i<n; i++) {
        !           291:       if (tf->e[i].x) return(1);
        !           292:       if (tf->e[i].D) return(1);
        !           293:     }
        !           294:     f = f->next;
        !           295:   }
        !           296:   return(0);
        !           297:
        !           298: }
        !           299:
        !           300: POLY homogenize(f)
        !           301: POLY f;
        !           302: /* homogenize by using (*grade)(f) */
        !           303: {
        !           304:   POLY t;
        !           305:   int maxg;
        !           306:   int flag,d;
        !           307:
        !           308:   if (f == ZERO) return(f);
        !           309:   t = f; maxg = (*grade)(f); flag = 0;
        !           310:   while (t != POLYNULL) {
        !           311:     d = (*grade)(t);
        !           312:     if (d != maxg) flag = 1;
        !           313:     if (d > maxg) {
        !           314:       maxg = d;
        !           315:     }
        !           316:     t = t->next;
        !           317:   }
        !           318:   if (flag == 0) return(f);
        !           319:
        !           320:   f = pmCopy(f); /* You can rewrite the monomial parts */
        !           321:   t = f;
        !           322:   while (t != POLYNULL) {
        !           323:     d = (*grade)(t);
        !           324:     if (d != maxg) {
        !           325:       t->m->e[0].D += maxg-d; /* Multiply h^(maxg-d) */
        !           326:     }
        !           327:     t = t->next;
        !           328:   }
        !           329:   return(f);
        !           330: }
        !           331:
        !           332: int isHomogenized(f)
        !           333: POLY f;
        !           334: {
        !           335:   POLY t;
        !           336:   extern int Homogenize_vec;
        !           337:   int maxg;
        !           338:   if (!Homogenize_vec) return(isHomogenized_vec(f));
        !           339:   if (f == ZERO) return(1);
        !           340:   maxg = (*grade)(f);
        !           341:   t = f;
        !           342:   while (t != POLYNULL) {
        !           343:     if (maxg != (*grade)(t)) return(0);
        !           344:     t = t->next;
        !           345:   }
        !           346:   return(1);
        !           347: }
        !           348:
        !           349: int isHomogenized_vec(f)
        !           350: POLY f;
        !           351: {
        !           352: /* This is not efficient version. *grade should be grade_module1v(). */
        !           353:   POLY t;
        !           354:   int ggg;
        !           355:   if (f == ZERO) return(1);
        !           356:   while (f != POLYNULL) {
        !           357:     t = f;
        !           358:     ggg = (*grade)(f);
        !           359:     while (t != POLYNULL) {
        !           360:       if ((*isSameComponent)(f,t)) {
        !           361:        if (ggg != (*grade)(t)) return(0);
        !           362:       }
        !           363:       t = t->next;
        !           364:     }
        !           365:     f = f->next;
        !           366:   }
        !           367:   return(1);
        !           368: }
        !           369:
        !           370:
        !           371: static int degreeOfPrincipalPart(f)
        !           372: POLY f;
        !           373: {
        !           374:   int n,i,dd;
        !           375:   if (f ISZERO) return(0);
        !           376:   n = f->m->ringp->n; dd = 0;
        !           377:   /* D[0] is homogenization var */
        !           378:   for (i=1; i<n; i++) {
        !           379:     dd += f->m->e[i].D;
        !           380:   }
        !           381:   return(dd);
        !           382: }
        !           383:
        !           384: POLY POLYToPrincipalPart(f)
        !           385: POLY f;
        !           386: {
        !           387:   POLY node;
        !           388:   struct listPoly nod;
        !           389:   POLY h;
        !           390:   POLY g;
        !           391:   int maxd = -20000; /* very big negative number */
        !           392:   int dd;
        !           393:   node = &nod; node->next = POLYNULL; h = node;
        !           394:
        !           395:   g = pCopy(f); /* shallow copy */
        !           396:   while (!(f ISZERO)) {
        !           397:     dd = degreeOfPrincipalPart(f);
        !           398:     if (dd > maxd) maxd = dd;
        !           399:     f = f->next;
        !           400:   }
        !           401:   while (!(g ISZERO)) {
        !           402:     dd = degreeOfPrincipalPart(g);
        !           403:     if (dd == maxd) {
        !           404:       h->next = g;
        !           405:       h = h->next;
        !           406:     }
        !           407:     g = g->next;
        !           408:   }
        !           409:   h->next = POLYNULL;
        !           410:   return(node->next);
        !           411: }
        !           412:
        !           413: static int degreeOfInitW(f,w)
        !           414: POLY f;
        !           415: int w[];
        !           416: {
        !           417:   int n,i,dd;
        !           418:   if (f ISZERO) {
        !           419:     errorPoly("degreeOfInitW(0,w) ");
        !           420:   }
        !           421:   n = f->m->ringp->n; dd = 0;
        !           422:   for (i=0; i<n; i++) {
        !           423:     dd += (f->m->e[i].D)*w[n+i];
        !           424:     dd += (f->m->e[i].x)*w[i];
        !           425:   }
        !           426:   return(dd);
        !           427: }
        !           428:
        !           429: POLY POLYToInitW(f,w)
        !           430: POLY f;
        !           431: int w[]; /* weight vector */
        !           432: {
        !           433:   POLY node;
        !           434:   struct listPoly nod;
        !           435:   POLY h;
        !           436:   POLY g;
        !           437:   int maxd;
        !           438:   int dd;
        !           439:   node = &nod; node->next = POLYNULL; h = node;
        !           440:
        !           441:   if (f ISZERO) return(f);
        !           442:   maxd = degreeOfInitW(f,w);
        !           443:   g = pCopy(f); /* shallow copy */
        !           444:   while (!(f ISZERO)) {
        !           445:     dd = degreeOfInitW(f,w);
        !           446:     if (dd > maxd) maxd = dd;
        !           447:     f = f->next;
        !           448:   }
        !           449:   while (!(g ISZERO)) {
        !           450:     dd = degreeOfInitW(g,w);
        !           451:     if (dd == maxd) {
        !           452:       h->next = g;
        !           453:       h = h->next;
        !           454:     }
        !           455:     g = g->next;
        !           456:   }
        !           457:   h->next = POLYNULL;
        !           458:   return(node->next);
        !           459: }
        !           460:
        !           461:
        !           462: /*
        !           463: 1.The substitution  "ringp->multiplication = ...." is allowed only in
        !           464:   KsetUpRing(), so the check in KswitchFunction is not necessary.
        !           465: 2.mmLarger != matrix and AvoidTheSameRing==1, then error
        !           466: 3.If Schreyer = 1, then the system always generates a new ring.
        !           467: 4.The execution of set_order_by_matrix is not allowed when Avoid... == 1.
        !           468: 5.When mmLarger == tower  (in tower.sm1, tower-sugar.sm1), we do
        !           469:   as follows with our own risk.
        !           470: [(AvoidTheSameRing)] pushEnv [ [(AvoidTheSameRing) 0] system_variable (mmLarger) (tower) switch_function ] pop popEnv
        !           471: */
        !           472: int isTheSameRing(struct ring *rstack[],int rp, struct ring *newRingp)
        !           473: {
        !           474:   struct ring *rrr;
        !           475:   int i,j,k;
        !           476:   int a=0;
        !           477:   for (k=0; k<rp; k++) {
        !           478:     rrr = rstack[k];
        !           479:     if (rrr->p != newRingp->p) { a=1; goto bbb ; }
        !           480:     if (rrr->n != newRingp->n) { a=2; goto bbb ; }
        !           481:     if (rrr->nn != newRingp->nn) { a=3; goto bbb ; }
        !           482:     if (rrr->m != newRingp->m) { a=4; goto bbb ; }
        !           483:     if (rrr->mm != newRingp->mm) { a=5; goto bbb ; }
        !           484:     if (rrr->l != newRingp->l) { a=6; goto bbb ; }
        !           485:     if (rrr->ll != newRingp->ll) { a=7; goto bbb ; }
        !           486:     if (rrr->c != newRingp->c) { a=8; goto bbb ; }
        !           487:     if (rrr->cc != newRingp->cc) { a=9; goto bbb ; }
        !           488:     for (i=0; i<rrr->n; i++) {
        !           489:       if (strcmp(rrr->x[i],newRingp->x[i])!=0) { a=10; goto bbb ; }
        !           490:       if (strcmp(rrr->D[i],newRingp->D[i])!=0) { a=11; goto bbb ; }
        !           491:     }
        !           492:     if (rrr->orderMatrixSize != newRingp->orderMatrixSize) { a=12; goto bbb ; }
        !           493:     for (i=0; i<rrr->orderMatrixSize; i++) {
        !           494:       for (j=0; j<2*(rrr->n); j++) {
        !           495:        if (rrr->order[i*2*(rrr->n)+j] != newRingp->order[i*2*(rrr->n)+j])
        !           496:          { a=13; goto bbb ; }
        !           497:       }
        !           498:     }
        !           499:     if (rrr->next != newRingp->next) { a=14; goto bbb ; }
        !           500:     if (rrr->multiplication != newRingp->multiplication) { a=15; goto bbb ; }
        !           501:     /* if (rrr->schreyer != newRingp->schreyer) { a=16; goto bbb ; }*/
        !           502:     if (newRingp->schreyer == 1) { a=16; goto bbb; }
        !           503:     /* The following fields are ignored.
        !           504:        void *gbListTower;
        !           505:        int *outputOrder;
        !           506:        char *name;
        !           507:     */
        !           508:     /* All tests are passed. */
        !           509:     return(k);
        !           510:   bbb: ;
        !           511:     /* for debugging. */
        !           512:     /* fprintf(stderr," reason=%d, ",a); */
        !           513:   }
        !           514:   return(-1);
        !           515: }
        !           516:
        !           517:

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