[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.13

1.13    ! takayama    1: /* $OpenXM: OpenXM/src/kan96xx/Kan/poly4.c,v 1.12 2003/08/24 05:19:42 takayama Exp $ */
1.1       maekawa     2: #include <stdio.h>
                      3: #include "datatype.h"
                      4: #include "stackm.h"
                      5: #include "extern.h"
                      6: #include "extern2.h"
                      7: #include "matrix.h"
                      8: static void shell(int v[],int n);
                      9: static int degreeOfPrincipalPart(POLY f);
                     10: static int degreeOfInitW(POLY f,int w[]);
1.10      takayama   11: static int degreeOfInitWS(POLY f,int w[],int s[]);
1.13    ! takayama   12: static int dDegree(POLY f);
        !            13: static POLY dHomogenize(POLY f);
1.1       maekawa    14:
                     15: static void shell(v,n)
1.3       takayama   16:      int v[];
                     17:      int n;
1.1       maekawa    18: {
                     19:   int gap,i,j,temp;
                     20:
                     21:   for (gap = n/2; gap > 0; gap /= 2) {
                     22:     for (i = gap; i<n; i++) {
                     23:       for (j=i-gap ; j>=0 && v[j]<v[j+gap]; j -= gap) {
1.3       takayama   24:         temp = v[j];
                     25:         v[j] = v[j+gap];
                     26:         v[j+gap] = temp;
1.1       maekawa    27:       }
                     28:     }
                     29:   }
                     30: }
                     31:
                     32:
                     33: struct matrixOfPOLY *parts(f,v)
1.3       takayama   34:      POLY f;
                     35:      POLY v;  /* v must be a single variable, e.g. x */
1.1       maekawa    36: {
                     37:   struct matrixOfPOLY *evPoly;
                     38:   int vi = 0;  /* index of v */
                     39:   int vx = 1;  /* x --> 1, D--> 0 */
                     40:   int n,evSize,i,k,e;
                     41:   int *ev;
                     42:   struct object *evList;
                     43:   struct object *list;
                     44:   struct object ob;
                     45:   POLY ans;
                     46:   POLY h;
                     47:   extern struct ring *CurrentRingp;
                     48:   POLY ft;
                     49:
                     50:
                     51:   if (f ISZERO || v ISZERO) {
                     52:     evPoly = newMatrixOfPOLY(2,1);
                     53:     getMatrixOfPOLY(evPoly,0,0) = ZERO;
                     54:     getMatrixOfPOLY(evPoly,1,0) = ZERO;
                     55:     return(evPoly);
                     56:   }
                     57:   n = v->m->ringp->n;
                     58:   /* get the index of the variable v */
                     59:   for (i=0; i<n; i++) {
                     60:     if (v->m->e[i].x) {
                     61:       vx = 1; vi = i; break;
                     62:     }else if (v->m->e[i].D) {
                     63:       vx = 0; vi = i; break;
                     64:     }
                     65:   }
                     66:   ft = f;
                     67:   /* get the vector of exponents */
                     68:   evList = NULLLIST;
                     69:   while (ft != POLYNULL) {
                     70:     if (vx) {
                     71:       e = ft->m->e[vi].x;
                     72:     }else{
                     73:       e = ft->m->e[vi].D;
                     74:     }
                     75:     ft = ft->next;
                     76:     ob = KpoInteger(e);
                     77:     if (!memberQ(evList,ob)) {
                     78:       list = newList(&ob);
                     79:       evList = vJoin(evList,list);
                     80:     }
                     81:   }
                     82:   /*printf("evList = "); printObjectList(evList);*/
                     83:   evSize = klength(evList);
                     84:   ev = (int *)sGC_malloc(sizeof(int)*(evSize+1));
                     85:   if (ev == (int *)NULL) errorPoly("No more memory.");
                     86:   for (i=0; i<evSize; i++) {
                     87:     ev[i] = KopInteger(car(evList));
                     88:     evList = cdr(evList);
                     89:   }
                     90:   /* sort ev */
                     91:   shell(ev,evSize);
                     92:
                     93:   /* get coefficients */
                     94:   evPoly = newMatrixOfPOLY(2,evSize);
                     95:   for (i=0; i<evSize; i++) {
                     96:     ans = ZERO;
                     97:     getMatrixOfPOLY(evPoly,0,i) = cxx(ev[i],0,0,CurrentRingp);
                     98:     ft = f;
                     99:     while (ft != POLYNULL) {
                    100:       if (vx) {
1.3       takayama  101:         if (ft->m->e[vi].x == ev[i]) {
                    102:           h = newCell(ft->coeffp,monomialCopy(ft->m));
                    103:           xset0(h,vi); /* touch monomial part, so you need to copy it above. */
                    104:           ans = ppAdd(ans,h);
                    105:         }
1.1       maekawa   106:       }else{
1.3       takayama  107:         if (ft->m->e[vi].D == ev[i]) {
                    108:           h = newCell(ft->coeffp,monomialCopy(ft->m));
                    109:           dset0(h,vi);
                    110:           ans = ppAdd(ans,h);
                    111:         }
1.1       maekawa   112:       }
                    113:       ft = ft->next;
                    114:     }
                    115:     getMatrixOfPOLY(evPoly,1,i) = ans;
                    116:   }
                    117:   return(evPoly);
                    118: }
1.3       takayama  119:
1.1       maekawa   120: struct object parts2(f,v)
1.3       takayama  121:      POLY f;
                    122:      POLY v;  /* v must be a single variable, e.g. x */
1.1       maekawa   123: {
                    124:   struct matrixOfPOLY *evPoly;
                    125:   int vi = 0;  /* index of v */
                    126:   int vx = 1;  /* x --> 1, D--> 0 */
                    127:   int n,evSize,i,k,e;
                    128:   int *ev;
                    129:   struct object *evList;
                    130:   struct object *list;
                    131:   struct object ob;
                    132:   POLY ans;
                    133:   POLY h;
                    134:   POLY ft;
                    135:   struct object ob1,ob2,rob;
                    136:
                    137:
                    138:   if (f ISZERO || v ISZERO) {
                    139:     evPoly = newMatrixOfPOLY(2,1);
                    140:     getMatrixOfPOLY(evPoly,0,0) = ZERO;
                    141:     getMatrixOfPOLY(evPoly,1,0) = ZERO;
                    142:     rob = newObjectArray(2);
                    143:     ob1 = newObjectArray(1);
                    144:     ob2 = newObjectArray(1);
                    145:     putoa(ob1,0,KpoInteger(0));
                    146:     putoa(ob2,0,KpoPOLY(POLYNULL));
                    147:     putoa(rob,0,ob1); putoa(rob,1,ob2);
                    148:     return(rob);
                    149:   }
                    150:   n = v->m->ringp->n;
                    151:   /* get the index of the variable v */
                    152:   for (i=0; i<n; i++) {
                    153:     if (v->m->e[i].x) {
                    154:       vx = 1; vi = i; break;
                    155:     }else if (v->m->e[i].D) {
                    156:       vx = 0; vi = i; break;
                    157:     }
                    158:   }
                    159:   ft = f;
                    160:   /* get the vector of exponents */
                    161:   evList = NULLLIST;
                    162:   while (ft != POLYNULL) {
                    163:     if (vx) {
                    164:       e = ft->m->e[vi].x;
                    165:     }else{
                    166:       e = ft->m->e[vi].D;
                    167:     }
                    168:     ft = ft->next;
                    169:     ob = KpoInteger(e);
                    170:     if (!memberQ(evList,ob)) {
                    171:       list = newList(&ob);
                    172:       evList = vJoin(evList,list);
                    173:     }
                    174:   }
                    175:   /*printf("evList = "); printObjectList(evList);*/
                    176:   evSize = klength(evList);
                    177:   ev = (int *)sGC_malloc(sizeof(int)*(evSize+1));
                    178:   if (ev == (int *)NULL) errorPoly("No more memory.");
                    179:   for (i=0; i<evSize; i++) {
                    180:     ev[i] = KopInteger(car(evList));
                    181:     evList = cdr(evList);
                    182:   }
                    183:   /* sort ev */
                    184:   shell(ev,evSize);
                    185:
                    186:   /* get coefficients */
                    187:   evPoly = newMatrixOfPOLY(2,evSize);
                    188:   for (i=0; i<evSize; i++) {
                    189:     ans = ZERO;
                    190:     /* getMatrixOfPOLY(evPoly,0,i) = cxx(ev[i],0,0,CurrentRingp); */
                    191:     getMatrixOfPOLY(evPoly,0,i) = POLYNULL;
                    192:     ft = f;
                    193:     while (ft != POLYNULL) {
                    194:       if (vx) {
1.3       takayama  195:         if (ft->m->e[vi].x == ev[i]) {
                    196:           h = newCell(ft->coeffp,monomialCopy(ft->m));
                    197:           xset0(h,vi); /* touch monomial part, so you need to copy it above. */
                    198:           ans = ppAdd(ans,h);
                    199:         }
1.1       maekawa   200:       }else{
1.3       takayama  201:         if (ft->m->e[vi].D == ev[i]) {
                    202:           h = newCell(ft->coeffp,monomialCopy(ft->m));
                    203:           dset0(h,vi);
                    204:           ans = ppAdd(ans,h);
                    205:         }
1.1       maekawa   206:       }
                    207:       ft = ft->next;
                    208:     }
                    209:     getMatrixOfPOLY(evPoly,1,i) = ans;
                    210:   }
                    211:   rob = newObjectArray(2);
                    212:   ob1 = newObjectArray(evSize);
                    213:   ob2 = newObjectArray(evSize);
                    214:   for (i=0; i<evSize; i++) {
                    215:     putoa(ob2,i,KpoPOLY(getMatrixOfPOLY(evPoly,1,i)));
                    216:     putoa(ob1,i,KpoInteger(ev[i]));
                    217:   }
                    218:   putoa(rob,0,ob1); putoa(rob,1,ob2);
                    219:   return(rob);
                    220: }
1.3       takayama  221:
1.1       maekawa   222: int pDegreeWrtV(f,v)
1.3       takayama  223:      POLY f;
                    224:      POLY v;
1.1       maekawa   225: {
                    226:   int vx = 1;
                    227:   int vi = 0;
                    228:   int i,n;
                    229:   int ans;
                    230:   if (f ISZERO || v ISZERO) return(0);
                    231:   n = f->m->ringp->n;
                    232:   for (i=0; i<n; i++) {
                    233:     if (v->m->e[i].x) {
                    234:       vx = 1; vi = i;
                    235:       break;
                    236:     }else if (v->m->e[i].D) {
                    237:       vx = 0; vi = i;
                    238:       break;
                    239:     }
                    240:   }
                    241:   if (vx) {
                    242:     ans = f->m->e[vi].x;
                    243:   }else{
                    244:     ans = f->m->e[vi].D;
                    245:   }
                    246:   f = f->next;
                    247:   while (f != POLYNULL) {
                    248:     if (vx) {
                    249:       if (f->m->e[vi].x > ans) ans = f->m->e[vi].x;
                    250:     }else{
                    251:       if (f->m->e[vi].D > ans) ans = f->m->e[vi].D;
                    252:     }
                    253:     f = f->next;
                    254:   }
                    255:   return(ans);
                    256: }
                    257:
                    258: int containVectorVariable(POLY f)
                    259: {
                    260:   MONOMIAL tf;
                    261:   static int nn,mm,ll,cc,n,m,l,c;
                    262:   static struct ring *cr = (struct ring *)NULL;
                    263:   int i;
                    264:
                    265:   if (f ISZERO) return(0);
                    266:   tf = f->m;
                    267:   if (tf->ringp != cr) {
                    268:     n = tf->ringp->n;
                    269:     m = tf->ringp->m;
                    270:     l = tf->ringp->l;
                    271:     c = tf->ringp->c;
                    272:     nn = tf->ringp->nn;
                    273:     mm = tf->ringp->mm;
                    274:     ll = tf->ringp->ll;
                    275:     cc = tf->ringp->cc;
                    276:     cr = tf->ringp;
                    277:   }
                    278:
                    279:   while (f != POLYNULL) {
                    280:     tf = f->m;
                    281:     for (i=cc; i<c; i++) {
                    282:       if ( tf->e[i].x ) return(1);
                    283:       if ( tf->e[i].D ) return(1);
                    284:     }
                    285:     for (i=ll; i<l; i++) {
                    286:       if (tf->e[i].x) return(1);
                    287:       if (tf->e[i].D) return(1);
                    288:     }
                    289:     for (i=mm; i<m; i++) {
                    290:       if (tf->e[i].x) return(1);
                    291:       if (tf->e[i].D) return(1);
                    292:     }
                    293:     for (i=nn; i<n; i++) {
                    294:       if (tf->e[i].x) return(1);
                    295:       if (tf->e[i].D) return(1);
                    296:     }
                    297:     f = f->next;
                    298:   }
                    299:   return(0);
                    300:
                    301: }
                    302:
                    303: POLY homogenize(f)
1.3       takayama  304:      POLY f;
                    305:      /* homogenize by using (*grade)(f) */
1.1       maekawa   306: {
                    307:   POLY t;
                    308:   int maxg;
                    309:   int flag,d;
1.13    ! takayama  310:   extern int Homogenize;
1.1       maekawa   311:
                    312:   if (f == ZERO) return(f);
1.13    ! takayama  313:   if (Homogenize == 3) { /* double homogenization Dx x = x Dx + h H */
        !           314:     return dHomogenize(f);
        !           315:   }
1.1       maekawa   316:   t = f; maxg = (*grade)(f); flag = 0;
                    317:   while (t != POLYNULL) {
                    318:     d = (*grade)(t);
                    319:     if (d != maxg) flag = 1;
                    320:     if (d > maxg) {
                    321:       maxg = d;
                    322:     }
                    323:     t = t->next;
                    324:   }
                    325:   if (flag == 0) return(f);
                    326:
                    327:   f = pmCopy(f); /* You can rewrite the monomial parts */
                    328:   t = f;
                    329:   while (t != POLYNULL) {
                    330:     d = (*grade)(t);
                    331:     if (d != maxg) {
                    332:       t->m->e[0].D += maxg-d; /* Multiply h^(maxg-d) */
                    333:     }
                    334:     t = t->next;
                    335:   }
                    336:   return(f);
                    337: }
                    338:
                    339: int isHomogenized(f)
1.3       takayama  340:      POLY f;
1.1       maekawa   341: {
                    342:   POLY t;
                    343:   extern int Homogenize_vec;
                    344:   int maxg;
                    345:   if (!Homogenize_vec) return(isHomogenized_vec(f));
                    346:   if (f == ZERO) return(1);
1.4       takayama  347:   if (f->m->ringp->weightedHomogenization) {
                    348:        return 1; /* BUG: do not chech in case of one-zero homogenization */
                    349:   }
1.1       maekawa   350:   maxg = (*grade)(f);
                    351:   t = f;
                    352:   while (t != POLYNULL) {
                    353:     if (maxg != (*grade)(t)) return(0);
                    354:     t = t->next;
                    355:   }
                    356:   return(1);
                    357: }
                    358:
                    359: int isHomogenized_vec(f)
1.3       takayama  360:      POLY f;
1.1       maekawa   361: {
1.3       takayama  362:   /* This is not efficient version. *grade should be grade_module1v(). */
1.1       maekawa   363:   POLY t;
                    364:   int ggg;
                    365:   if (f == ZERO) return(1);
1.4       takayama  366:   if (f->m->ringp->weightedHomogenization) {
                    367:        return 1; /* BUG: do not chech in case of one-zero homogenization */
                    368:   }
1.1       maekawa   369:   while (f != POLYNULL) {
                    370:     t = f;
                    371:     ggg = (*grade)(f);
                    372:     while (t != POLYNULL) {
                    373:       if ((*isSameComponent)(f,t)) {
1.3       takayama  374:         if (ggg != (*grade)(t)) return(0);
1.1       maekawa   375:       }
                    376:       t = t->next;
                    377:     }
                    378:     f = f->next;
                    379:   }
                    380:   return(1);
                    381: }
                    382:
1.13    ! takayama  383: static POLY dHomogenize(f)
        !           384: POLY f;
        !           385: {
        !           386:   POLY t;
        !           387:   int maxg, maxdg;
        !           388:   int flag,d,dd,neg;
        !           389:
        !           390:   if (f == ZERO) return(f);
        !           391:   t = f; maxg = (*grade)(f); flag = 0;
        !           392:   maxdg = dDegree(f);
        !           393:   while (t != POLYNULL) {
        !           394:     d = (*grade)(t);
        !           395:     if (d != maxg) flag = 1;
        !           396:     if (d > maxg) {
        !           397:       maxg = d;
        !           398:     }
        !           399:     d = dDegree(f);
        !           400:     if (d > maxdg) {
        !           401:       maxdg = d;
        !           402:     }
        !           403:     t = t->next;
        !           404:   }
        !           405:   if (flag == 0) return(f);
        !           406:
        !           407:   t = f; neg = 0;
        !           408:   while (t != POLYNULL) {
        !           409:     d = (*grade)(t);
        !           410:     dd = dDegree(t);
        !           411:     if (maxg-d-(maxdg-dd) < neg) {
        !           412:       neg = maxg-d-(maxdg-dd);
        !           413:     }
        !           414:     t = t->next;
        !           415:   }
        !           416:   neg = -neg;
        !           417:
        !           418:   f = pmCopy(f); /* You can rewrite the monomial parts */
        !           419:   t = f;
        !           420:   while (t != POLYNULL) {
        !           421:     d = (*grade)(t);
        !           422:     dd = dDegree(t);
        !           423:     t->m->e[0].D += maxdg-dd; /* h */
        !           424:     t->m->e[0].x += maxg-d-(maxdg-dd)+neg; /* Multiply H */
        !           425:     /* example Dx^2+Dx+x */
        !           426:     t = t->next;
        !           427:   }
        !           428:   return(f);
        !           429: }
1.1       maekawa   430:
                    431: static int degreeOfPrincipalPart(f)
1.3       takayama  432:      POLY f;
1.1       maekawa   433: {
                    434:   int n,i,dd;
                    435:   if (f ISZERO) return(0);
                    436:   n = f->m->ringp->n; dd = 0;
                    437:   /* D[0] is homogenization var */
                    438:   for (i=1; i<n; i++) {
1.13    ! takayama  439:     dd += f->m->e[i].D;
        !           440:   }
        !           441:   return(dd);
        !           442: }
        !           443:
        !           444: static int dDegree(f)
        !           445:      POLY f;
        !           446: {
        !           447:   int nn,i,dd,m;
        !           448:   if (f ISZERO) return(0);
        !           449:   nn = f->m->ringp->nn; dd = 0;
        !           450:   m = f->m->ringp->m;
        !           451:   for (i=m; i<nn; i++) {
1.1       maekawa   452:     dd += f->m->e[i].D;
                    453:   }
                    454:   return(dd);
                    455: }
                    456:
                    457: POLY POLYToPrincipalPart(f)
1.3       takayama  458:      POLY f;
1.1       maekawa   459: {
                    460:   POLY node;
                    461:   struct listPoly nod;
                    462:   POLY h;
                    463:   POLY g;
                    464:   int maxd = -20000; /* very big negative number */
                    465:   int dd;
                    466:   node = &nod; node->next = POLYNULL; h = node;
                    467:
                    468:   g = pCopy(f); /* shallow copy */
                    469:   while (!(f ISZERO)) {
                    470:     dd = degreeOfPrincipalPart(f);
                    471:     if (dd > maxd) maxd = dd;
                    472:     f = f->next;
                    473:   }
                    474:   while (!(g ISZERO)) {
                    475:     dd = degreeOfPrincipalPart(g);
                    476:     if (dd == maxd) {
                    477:       h->next = g;
                    478:       h = h->next;
                    479:     }
                    480:     g = g->next;
                    481:   }
                    482:   h->next = POLYNULL;
                    483:   return(node->next);
                    484: }
                    485:
                    486: static int degreeOfInitW(f,w)
1.3       takayama  487:      POLY f;
                    488:      int w[];
1.1       maekawa   489: {
                    490:   int n,i,dd;
                    491:   if (f ISZERO) {
                    492:     errorPoly("degreeOfInitW(0,w) ");
                    493:   }
                    494:   n = f->m->ringp->n; dd = 0;
                    495:   for (i=0; i<n; i++) {
                    496:     dd += (f->m->e[i].D)*w[n+i];
                    497:     dd += (f->m->e[i].x)*w[i];
                    498:   }
                    499:   return(dd);
                    500: }
                    501:
                    502: POLY POLYToInitW(f,w)
1.3       takayama  503:      POLY f;
                    504:      int w[]; /* weight vector */
1.1       maekawa   505: {
                    506:   POLY h;
                    507:   POLY g;
                    508:   int maxd;
                    509:   int dd;
1.11      takayama  510:   h = POLYNULL;
1.1       maekawa   511:
1.11      takayama  512:   /*printf("1:%s\n",POLYToString(f,'*',1));*/
1.1       maekawa   513:   if (f ISZERO) return(f);
                    514:   maxd = degreeOfInitW(f,w);
1.11      takayama  515:   g = f;
1.1       maekawa   516:   while (!(f ISZERO)) {
                    517:     dd = degreeOfInitW(f,w);
                    518:     if (dd > maxd) maxd = dd;
                    519:     f = f->next;
                    520:   }
                    521:   while (!(g ISZERO)) {
                    522:     dd = degreeOfInitW(g,w);
                    523:     if (dd == maxd) {
1.11      takayama  524:       h = ppAdd(h,newCell(g->coeffp,g->m)); /* it might be slow. */
1.1       maekawa   525:     }
                    526:     g = g->next;
                    527:   }
1.11      takayama  528:   /*printf("2:%s\n",POLYToString(h,'*',1));*/
                    529:   return(h);
1.10      takayama  530: }
                    531:
                    532: static int degreeOfInitWS(f,w,s)
                    533:      POLY f;
                    534:      int w[];
                    535:         int s[];
                    536: {
                    537:   int n,i,dd;
                    538:   if (f ISZERO) {
                    539:     errorPoly("degreeOfInitWS(0,w) ");
                    540:   }
1.11      takayama  541:   if (s == (int *) NULL) return degreeOfInitW(f,w);
1.10      takayama  542:   n = f->m->ringp->n; dd = 0;
                    543:   for (i=0; i<n-1; i++) {
                    544:     dd += (f->m->e[i].D)*w[n+i];
                    545:     dd += (f->m->e[i].x)*w[i];
                    546:   }
                    547:   dd += s[(f->m->e[n-1].x)];
                    548:   return(dd);
                    549: }
                    550:
                    551: POLY POLYToInitWS(f,w,s)
                    552:      POLY f;
                    553:      int w[]; /* weight vector */
                    554:         int s[]; /* shift vector */
                    555: {
                    556:   POLY h;
                    557:   POLY g;
                    558:   int maxd;
                    559:   int dd;
1.11      takayama  560:   h = POLYNULL;
1.10      takayama  561:
1.11      takayama  562:   /*printf("1s:%s\n",POLYToString(f,'*',1));*/
1.10      takayama  563:   if (f ISZERO) return(f);
                    564:   maxd = degreeOfInitWS(f,w,s);
1.11      takayama  565:   g = f;
1.10      takayama  566:   while (!(f ISZERO)) {
                    567:     dd = degreeOfInitWS(f,w,s);
                    568:     if (dd > maxd) maxd = dd;
                    569:     f = f->next;
                    570:   }
                    571:   while (!(g ISZERO)) {
                    572:     dd = degreeOfInitWS(g,w,s);
                    573:     if (dd == maxd) {
1.11      takayama  574:       h = ppAdd(h,newCell(g->coeffp,g->m)); /* it might be slow. */
1.10      takayama  575:     }
                    576:     g = g->next;
                    577:   }
1.11      takayama  578:   /*printf("2s:%s\n",POLYToString(h,'*',1));*/
                    579:   return(h);
1.10      takayama  580: }
                    581:
                    582: int ordWsAll(f,w,s)
                    583:      POLY f;
                    584:      int w[]; /* weight vector */
                    585:         int s[]; /* shift vector */
                    586: {
                    587:   int maxd;
                    588:   int dd;
                    589:
                    590:   if (f ISZERO)  errorPoly("ordWsAll(0,w,s) ");
                    591:   maxd = degreeOfInitWS(f,w,s);
                    592:   while (!(f ISZERO)) {
                    593:     dd = degreeOfInitWS(f,w,s);
                    594:     if (dd > maxd) maxd = dd;
                    595:     f = f->next;
                    596:   }
                    597:   return maxd;
1.1       maekawa   598: }
                    599:
                    600:
                    601: /*
                    602: 1.The substitution  "ringp->multiplication = ...." is allowed only in
                    603:   KsetUpRing(), so the check in KswitchFunction is not necessary.
                    604: 2.mmLarger != matrix and AvoidTheSameRing==1, then error
                    605: 3.If Schreyer = 1, then the system always generates a new ring.
                    606: 4.The execution of set_order_by_matrix is not allowed when Avoid... == 1.
                    607: 5.When mmLarger == tower  (in tower.sm1, tower-sugar.sm1), we do
                    608:   as follows with our own risk.
                    609: [(AvoidTheSameRing)] pushEnv [ [(AvoidTheSameRing) 0] system_variable (mmLarger) (tower) switch_function ] pop popEnv
                    610: */
                    611: int isTheSameRing(struct ring *rstack[],int rp, struct ring *newRingp)
                    612: {
                    613:   struct ring *rrr;
                    614:   int i,j,k;
                    615:   int a=0;
                    616:   for (k=0; k<rp; k++) {
                    617:     rrr = rstack[k];
                    618:     if (rrr->p != newRingp->p) { a=1; goto bbb ; }
                    619:     if (rrr->n != newRingp->n) { a=2; goto bbb ; }
                    620:     if (rrr->nn != newRingp->nn) { a=3; goto bbb ; }
                    621:     if (rrr->m != newRingp->m) { a=4; goto bbb ; }
                    622:     if (rrr->mm != newRingp->mm) { a=5; goto bbb ; }
                    623:     if (rrr->l != newRingp->l) { a=6; goto bbb ; }
                    624:     if (rrr->ll != newRingp->ll) { a=7; goto bbb ; }
                    625:     if (rrr->c != newRingp->c) { a=8; goto bbb ; }
                    626:     if (rrr->cc != newRingp->cc) { a=9; goto bbb ; }
                    627:     for (i=0; i<rrr->n; i++) {
                    628:       if (strcmp(rrr->x[i],newRingp->x[i])!=0) { a=10; goto bbb ; }
                    629:       if (strcmp(rrr->D[i],newRingp->D[i])!=0) { a=11; goto bbb ; }
                    630:     }
                    631:     if (rrr->orderMatrixSize != newRingp->orderMatrixSize) { a=12; goto bbb ; }
                    632:     for (i=0; i<rrr->orderMatrixSize; i++) {
                    633:       for (j=0; j<2*(rrr->n); j++) {
1.3       takayama  634:         if (rrr->order[i*2*(rrr->n)+j] != newRingp->order[i*2*(rrr->n)+j])
                    635:           { a=13; goto bbb ; }
1.1       maekawa   636:       }
                    637:     }
                    638:     if (rrr->next != newRingp->next) { a=14; goto bbb ; }
                    639:     if (rrr->multiplication != newRingp->multiplication) { a=15; goto bbb ; }
                    640:     /* if (rrr->schreyer != newRingp->schreyer) { a=16; goto bbb ; }*/
                    641:     if (newRingp->schreyer == 1) { a=16; goto bbb; }
1.4       takayama  642:     if (rrr->weightedHomogenization != newRingp->weightedHomogenization) { a=16; goto bbb; }
1.7       takayama  643:     if (rrr->degreeShiftSize != newRingp->degreeShiftSize) {
                    644:       a = 17; goto bbb;
                    645:     }
                    646:     if (rrr->degreeShiftN != newRingp->degreeShiftN) {
                    647:       a = 17; goto bbb;
                    648:     }
                    649:     for (i=0; i < rrr->degreeShiftSize; i++) {
                    650:       for (j=0; j< rrr->degreeShiftN; j++) {
                    651:         if (rrr->degreeShift[i*(rrr->degreeShiftN)+j] !=
                    652:             newRingp->degreeShift[i*(rrr->degreeShiftN)+j]) {
                    653:           a = 17; goto bbb;
                    654:         }
                    655:       }
                    656:     }
                    657:
1.1       maekawa   658:     /* The following fields are ignored.
                    659:        void *gbListTower;
                    660:        int *outputOrder;
                    661:        char *name;
                    662:     */
                    663:     /* All tests are passed. */
                    664:     return(k);
                    665:   bbb: ;
                    666:     /* for debugging. */
                    667:     /* fprintf(stderr," reason=%d, ",a); */
                    668:   }
                    669:   return(-1);
                    670: }
1.6       takayama  671:
                    672: /* s->1 */
                    673: POLY goDeHomogenizeS(POLY f) {
1.9       takayama  674:   POLY lRule[1];
                    675:   POLY rRule[1];
                    676:   struct ring *rp;
                    677:   POLY ans;
                    678:   /* printf("1:[%s]\n",POLYToString(f,'*',1)); */
                    679:   if (f == POLYNULL) return f;
                    680:   rp = f->m->ringp;
                    681:   if (rp->next == NULL) {
                    682:        lRule[0] = cxx(1,0,1,rp);
                    683:        rRule[0] = cxx(1,0,0,rp);
                    684:        ans=replace(f,lRule,rRule,1);
                    685:   }else{
                    686:        struct coeff *cp;
                    687:        POLY t;
                    688:        POLY nc;
                    689:     ans = POLYNULL;
                    690:        while (f != POLYNULL) {
                    691:          cp = f->coeffp;
                    692:          if (cp->tag == POLY_COEFF) {
                    693:                t = goDeHomogenizeS((cp->val).f);
1.11      takayama  694:                nc = newCell(polyToCoeff(t,f->m->ringp),monomialCopy(f->m));
1.9       takayama  695:                ans = ppAddv(ans,nc);
                    696:                f = f->next;
                    697:          }else{
                    698:                ans = f; break;
                    699:          }
                    700:        }
                    701:   }
                    702:   /* printf("2:[%s]\n",POLYToString(ans,'*',1)); */
                    703:   return ans;
                    704: }
                    705:
                    706: POLY goDeHomogenizeS_buggy(POLY f) {
1.6       takayama  707:   POLY node;
                    708:   POLY lastf;
                    709:   struct listPoly nod;
                    710:   POLY h;
                    711:   POLY tf;
                    712:   int gt,first;
                    713:
1.9       takayama  714:   printf("1:[%s]\n",POLYToString(f,'*',1));
1.6       takayama  715:   if (f == POLYNULL) return(POLYNULL);
                    716:   node = &nod;
                    717:   node->next = POLYNULL;
                    718:   lastf = POLYNULL;
                    719:   first = 1;
                    720:   while (f != POLYNULL) {
                    721:     tf = newCell(f->coeffp,monomialCopy(f->m));
                    722:     tf->m->e[0].x = 0;  /* H, s variable in the G-O paper. */
                    723:     if (first) {
                    724:       node->next = tf;
                    725:       lastf = tf;
                    726:       first = 0;
                    727:     }else{
                    728:       gt = (*mmLarger)(lastf,tf);
                    729:       if (gt == 1) {
                    730:         lastf->next = tf;
                    731:         lastf = tf;
                    732:       }else{
                    733:         h = node->next;
                    734:         h = ppAddv(h,tf);
                    735:         node->next = h;
                    736:         lastf = h;
                    737:         while (lastf->next != POLYNULL) {
                    738:           lastf = lastf->next;
                    739:         }
                    740:       }
                    741:     }
                    742:     f = f->next;
                    743:   }
1.9       takayama  744:   printf("2:[%s]\n",POLYToString(node->next,'*',1));
1.6       takayama  745:   return (node->next);
                    746: }
                    747:
1.5       takayama  748: /* Granger-Oaku's homogenization for the ecart tangent cone.
                    749:    Note: 2003.07.10.
                    750:    ds[] is the degree shift.
                    751:    ei ( element index ). If it is < 0, then e[n-1]->x will be used,
                    752:                          else ei is used.
1.6       takayama  753:    if onlyS is set to 1, then input is assumed to be (u,v)-h-homogeneous.
1.5       takayama  754: */
1.6       takayama  755: POLY goHomogenize(POLY f,int u[],int v[],int ds[],int dssize,int ei,int onlyS)
1.5       takayama  756: {
                    757:   POLY node;
                    758:   POLY lastf;
                    759:   struct listPoly nod;
                    760:   POLY h;
                    761:   POLY tf;
                    762:   int gt,first,m,mp,t,tp,dsIdx,message;
1.7       takayama  763:   struct ring *rp;
1.5       takayama  764:
                    765:   message = 1;
                    766:   if (f == POLYNULL) return(POLYNULL);
1.7       takayama  767:   rp = f->m->ringp;
1.12      takayama  768:   /*
1.7       takayama  769:   if ((rp->degreeShiftSize == 0) && (dssize > 0)) {
                    770:        warningPoly("You are trying to homogenize a polynomial with degree shift. However, the polynomial belongs to the ring without degreeShift option. It may cause a trouble in comparison in free module.\n");
                    771:   }
1.12      takayama  772:   */
1.5       takayama  773:   node = &nod;
                    774:   node->next = POLYNULL;
                    775:   lastf = POLYNULL;
                    776:   first = 1;
                    777:   while (f != POLYNULL) {
                    778:     if (first) {
                    779:       t = m = dGrade1(f);
                    780:       tp = mp = uvGrade1(f,u,v,ds,dssize,ei);
                    781:     }else{
                    782:       t =  dGrade1(f);
                    783:       tp = uvGrade1(f,u,v,ds,dssize,ei);
                    784:       if (t > m) m = t;
                    785:       if (tp < mp) mp = tp;
                    786:     }
                    787:
                    788:     tf = newCell(f->coeffp,monomialCopy(f->m));
                    789:        /* Automatic dehomogenize. Not +=  */
                    790:        if (message && ((tf->m->e[0].D != 0) || (tf->m->e[0].x != 0))) {
                    791:       /*go-debug fprintf(stderr,"Automatic dehomogenize and homogenize.\n"); */
                    792:          message = 0;
                    793:        }
1.6       takayama  794:        if (!onlyS) {
                    795:          tf->m->e[0].D = -t;  /* h */
                    796:        }
1.5       takayama  797:     tf->m->e[0].x = tp;  /* H, s variable in the G-O paper. */
                    798:        /*go-debug printf("t(h)=%d, tp(uv+ds)=%d\n",t,tp); */
                    799:     if (first) {
                    800:       node->next = tf;
                    801:       lastf = tf;
                    802:       first = 0;
                    803:     }else{
                    804:       gt = (*mmLarger)(lastf,tf);
                    805:       if (gt == 1) {
                    806:         lastf->next = tf;
                    807:         lastf = tf;
                    808:       }else{
                    809:         /*go-debug printf("?\n"); */
                    810:         h = node->next;
                    811:         h = ppAddv(h,tf);
                    812:         node->next = h;
                    813:         lastf = h;
                    814:         while (lastf->next != POLYNULL) {
                    815:           lastf = lastf->next;
                    816:         }
                    817:       }
                    818:        }
                    819:        f = f->next;
                    820:   }
                    821:   h = node->next;
                    822:   /*go-debug printf("m=%d, mp=%d\n",m,mp); */
                    823:   while (h != POLYNULL) {
1.12      takayama  824:     /*go-debug printf("Old: h=%d, s=%d\n",h->m->e[0].D,h->m->e[0].x);  */
1.6       takayama  825:     if (!onlyS) h->m->e[0].D += m;   /* h */
1.5       takayama  826:     h->m->e[0].x += -mp; /* H, s*/
                    827:     /*go-debug printf("New: h=%d, s=%d\n",h->m->e[0].D,h->m->e[0].x); */
                    828:     h = h->next;
                    829:   }
                    830:   return (node->next);
                    831: }
                    832:
                    833: /* u[] = -1, v[] = 1 */
1.6       takayama  834: POLY goHomogenize11(POLY f,int ds[],int dssize,int ei,int onlyS)
1.5       takayama  835: {
                    836:   int r;
                    837:   int i,t,n,m,nn;
                    838:   MONOMIAL tf;
                    839:   static int *u;
                    840:   static int *v;
                    841:   static struct ring *cr = (struct ring *)NULL;
                    842:
                    843:   if (f == POLYNULL) return POLYNULL;
                    844:
                    845:   tf = f->m;
                    846:   if (tf->ringp != cr) {
                    847:     n = tf->ringp->n;
                    848:     m = tf->ringp->m;
                    849:     nn = tf->ringp->nn;
                    850:     cr = tf->ringp;
                    851:        u = (int *)sGC_malloc(sizeof(int)*n);
                    852:        v = (int *)sGC_malloc(sizeof(int)*n);
                    853:        for (i=0; i<n; i++) u[i]=v[i]=0;
                    854:        for (i=m; i<nn; i++) {
                    855:          u[i] = -1; v[i] = 1;
                    856:        }
                    857:   }
1.6       takayama  858:   return(goHomogenize(f,u,v,ds,dssize,ei,onlyS));
1.5       takayama  859: }
                    860:
1.6       takayama  861: POLY goHomogenize_dsIdx(POLY f,int u[],int v[],int dsIdx,int ei,int onlyS)
1.5       takayama  862: {
                    863:   if (f == POLYNULL) return POLYNULL;
                    864: }
1.6       takayama  865: POLY goHomogenize11_dsIdx(POLY f,int ds[],int dsIdx,int ei,int onlyS)
1.5       takayama  866: {
                    867:   if (f == POLYNULL) return POLYNULL;
1.8       takayama  868: }
                    869:
                    870: /* cf. KsetUpRing() in kanExport0.c */
                    871: struct ring *newRingOverFp(struct ring *rp,int p) {
                    872:   struct ring *newRingp;
                    873:   char *ringName = NULL;
                    874:   char pstr[64];
                    875:   sprintf(pstr,"_%d",p);
                    876:   ringName = (char *)sGC_malloc(128);
                    877:   newRingp = (struct ring *)sGC_malloc(sizeof(struct ring));
                    878:   if (newRingp == NULL) errorPoly("No more memory.\n");
                    879:   strcpy(ringName,rp->name);
                    880:   strcat(ringName,pstr);
                    881:   *newRingp = *rp;
                    882:   newRingp->p = p;
                    883:   newRingp->name = ringName;
                    884:   return newRingp;
                    885: }
                    886:
                    887: /*
                    888:    P = 3001;
                    889:    L = [ ];
                    890:    while (P<10000) {
                    891:      L=cons(P,L);
                    892:      P = pari(nextprime,P+1);
                    893:        }
                    894:    print(L);
                    895: */
                    896: #define N799  799
                    897: static int nextPrime(void) {
                    898:   static int pt = 0;
                    899:   static int tb[N799] =
                    900: {3001,3011,3019,3023,3037,3041,3049,3061,3067,3079,3083,3089,3109,3119,3121,3137,3163,3167,3169,3181,3187,3191,3203,3209,3217,3221,3229,3251,3253,3257,3259,3271,3299,3301,3307,3313,3319,3323,3329,3331,3343,3347,3359,3361,3371,3373,3389,3391,3407,3413,3433,3449,3457,3461,3463,3467,3469,3491,3499,3511,3517,3527,3529,3533,3539,3541,3547,3557,3559,3571,3581,3583,3593,3607,3613,3617,3623,3631,3637,3643,3659,3671,3673,3677,3691,3697,3701,3709,3719,3727,3733,3739,3761,3767,3769,3779,3793,3797,3803,3821,3823,3833,3847,3851,3853,3863,3877,3881,3889,3907,3911,3917,3919,3923,3929,3931,3943,3947,3967,3989,
                    901:   4001,4003,4007,4013,4019,4021,4027,4049,4051,4057,4073,4079,4091,4093,4099,4111,4127,4129,4133,4139,4153,4157,4159,4177,4201,4211,4217,4219,4229,4231,4241,4243,4253,4259,4261,4271,4273,4283,4289,4297,4327,4337,4339,4349,4357,4363,4373,4391,4397,4409,4421,4423,4441,4447,4451,4457,4463,4481,4483,4493,4507,4513,4517,4519,4523,4547,4549,4561,4567,4583,4591,4597,4603,4621,4637,4639,4643,4649,4651,4657,4663,4673,4679,4691,4703,4721,4723,4729,4733,4751,4759,4783,4787,4789,4793,4799,4801,4813,4817,4831,4861,4871,4877,4889,4903,4909,4919,4931,4933,4937,4943,4951,4957,4967,4969,4973,4987,4993,4999,
                    902:   5003,5009,5011,5021,5023,5039,5051,5059,5077,5081,5087,5099,5101,5107,5113,5119,5147,5153,5167,5171,5179,5189,5197,5209,5227,5231,5233,5237,5261,5273,5279,5281,5297,5303,5309,5323,5333,5347,5351,5381,5387,5393,5399,5407,5413,5417,5419,5431,5437,5441,5443,5449,5471,5477,5479,5483,5501,5503,5507,5519,5521,5527,5531,5557,5563,5569,5573,5581,5591,5623,5639,5641,5647,5651,5653,5657,5659,5669,5683,5689,5693,5701,5711,5717,5737,5741,5743,5749,5779,5783,5791,5801,5807,5813,5821,5827,5839,5843,5849,5851,5857,5861,5867,5869,5879,5881,5897,5903,5923,5927,5939,5953,5981,5987,
                    903:   6007,6011,6029,6037,6043,6047,6053,6067,6073,6079,6089,6091,6101,6113,6121,6131,6133,6143,6151,6163,6173,6197,6199,6203,6211,6217,6221,6229,6247,6257,6263,6269,6271,6277,6287,6299,6301,6311,6317,6323,6329,6337,6343,6353,6359,6361,6367,6373,6379,6389,6397,6421,6427,6449,6451,6469,6473,6481,6491,6521,6529,6547,6551,6553,6563,6569,6571,6577,6581,6599,6607,6619,6637,6653,6659,6661,6673,6679,6689,6691,6701,6703,6709,6719,6733,6737,6761,6763,6779,6781,6791,6793,6803,6823,6827,6829,6833,6841,6857,6863,6869,6871,6883,6899,6907,6911,6917,6947,6949,6959,6961,6967,6971,6977,6983,6991,6997,
                    904:   7001,7013,7019,7027,7039,7043,7057,7069,7079,7103,7109,7121,7127,7129,7151,7159,7177,7187,7193,7207,7211,7213,7219,7229,7237,7243,7247,7253,7283,7297,7307,7309,7321,7331,7333,7349,7351,7369,7393,7411,7417,7433,7451,7457,7459,7477,7481,7487,7489,7499,7507,7517,7523,7529,7537,7541,7547,7549,7559,7561,7573,7577,7583,7589,7591,7603,7607,7621,7639,7643,7649,7669,7673,7681,7687,7691,7699,7703,7717,7723,7727,7741,7753,7757,7759,7789,7793,7817,7823,7829,7841,7853,7867,7873,7877,7879,7883,7901,7907,7919,7927,7933,7937,7949,7951,7963,7993,
                    905: 8009,8011,8017,8039,8053,8059,8069,8081,8087,8089,8093,8101,8111,8117,8123,8147,8161,8167,8171,8179,8191,8209,8219,8221,8231,8233,8237,8243,8263,8269,8273,8287,8291,8293,8297,8311,8317,8329,8353,8363,8369,8377,8387,8389,8419,8423,8429,8431,8443,8447,8461,8467,8501,8513,8521,8527,8537,8539,8543,8563,8573,8581,8597,8599,8609,8623,8627,8629,8641,8647,8663,8669,8677,8681,8689,8693,8699,8707,8713,8719,8731,8737,8741,8747,8753,8761,8779,8783,8803,8807,8819,8821,8831,8837,8839,8849,8861,8863,8867,8887,8893,8923,8929,8933,8941,8951,8963,8969,8971,8999,
                    906:   9001,9007,9011,9013,9029,9041,9043,9049,9059,9067,9091,9103,9109,9127,9133,9137,9151,9157,9161,9173,9181,9187,9199,9203,9209,9221,9227,9239,9241,9257,9277,9281,9283,9293,9311,9319,9323,9337,9341,9343,9349,9371,9377,9391,9397,9403,9413,9419,9421,9431,9433,9437,9439,9461,9463,9467,9473,9479,9491,9497,9511,9521,9533,9539,9547,9551,9587,9601,9613,9619,9623,9629,9631,9643,9649,9661,9677,9679,9689,9697,9719,9721,9733,9739,9743,9749,9767,9769,9781,9787,9791,9803,9811,9817,9829,9833,9839,9851,9857,9859,9871,9883,9887,9901,9907,9923,9929,9931,9941,9949,9967,9973};
                    907:
                    908:   if (pt <N799) {
                    909:        return tb[pt++];
                    910:   }else{
                    911:        pt = 0;
                    912:        return tb[pt++];
                    913:   }
                    914: }
                    915:
                    916: int getPrime(int p) {
                    917:   int i;
                    918:   if (p <= 2) return nextPrime();
                    919:   for (i=2; i<p; i++) {
                    920:        if (p % i == 0) {
                    921:          return nextPrime();
                    922:        }
                    923:   }
                    924:   return p;
1.5       takayama  925: }

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