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

Annotation of OpenXM/src/kan96xx/Kan/order.c, Revision 1.16

1.16    ! takayama    1: /* $OpenXM: OpenXM/src/kan96xx/Kan/order.c,v 1.15 2005/07/03 11:08:54 ohara Exp $ */
1.1       maekawa     2: #include <stdio.h>
1.15      ohara       3: #include <stdlib.h>
1.1       maekawa     4: #include "datatype.h"
                      5: #include "stackm.h"
                      6: #include "extern.h"
                      7: #include "extern2.h"
                      8:
                      9: /* The format of order.
                     10:    Example:   graded lexicographic order
                     11:    x_{N-1}  x_{N-2}  ...  x_0  D_{N-1}  ....  D_{0}
                     12:     1        1             1    1              1
                     13:     1        0             0    0              0
                     14:     0        1             0    0              0
                     15:     ..............................................
                     16:
                     17:    (ringp->order)[i][j] should be (ringp->order)[i*2*N+j].
                     18:    All order matrix is generated by functions in smacro.sm1
                     19: */
                     20:
                     21: static void warningOrder(char *s);
                     22: static void errorOrder(char *s);
                     23:
                     24: void setOrderByMatrix(order,n,c,l,omsize)
1.4       takayama   25:      int order[];
                     26:      int n,c,l,omsize;
1.1       maekawa    27: {
                     28:   int i,j;
                     29:   int *Order;
                     30:   extern struct ring *CurrentRingp;
                     31:
                     32:   switch_mmLarger("default");
1.4       takayama   33:   /* q-case */
1.1       maekawa    34:   if ( l-c > 0) {
                     35:     switch_mmLarger("qmatrix");
                     36:   }
                     37:
                     38:   Order = (int *)sGC_malloc(sizeof(int)*(2*n)*(omsize));
                     39:   if (Order == (int *)NULL) errorOrder("No memory.");
                     40:   CurrentRingp->order = Order;
                     41:   CurrentRingp->orderMatrixSize = omsize;
                     42:   for (i=0; i<omsize; i++) {
                     43:     for (j=0; j<2*n; j++) {
                     44:       Order[i*2*n+j] = order[i*2*n+j];
                     45:     }
                     46:   }
                     47: }
                     48:
                     49: void showRing(level,ringp)
1.4       takayama   50:      int level;
                     51:      struct ring *ringp;
1.1       maekawa    52: {
                     53:   int i,j;
                     54:   FILE *fp;
                     55:   char tmp[100];
                     56:   int N,M,L,C,NN,MM,LL,CC;
                     57:   char **TransX,**TransD;
                     58:   int *Order;
                     59:   int P;
                     60:   char *mtype;
                     61:   extern char *F_isSameComponent;
1.5       takayama   62:   POLY f;
1.6       takayama   63:   POLY fx;
                     64:   POLY fd;
                     65:   POLY rf;
1.1       maekawa    66:   fp = stdout;
                     67:
                     68:   N=ringp->n; M = ringp->m; L = ringp->l; C = ringp->c;
                     69:   NN=ringp->nn; MM = ringp->mm; LL = ringp->ll; CC = ringp->cc;
                     70:   TransX = ringp->x; TransD = ringp->D;
                     71:   Order = ringp->order;
                     72:   P = ringp->p;
                     73:
                     74:
                     75:   fprintf(fp,"\n----------  the current ring ---- name: %s------\n",ringp->name);
                     76:   fprintf(fp,"Characteristic is %d. ",P);
                     77:   fprintf(fp,"N0=%d N=%d NN=%d M=%d MM=%d L=%d LL=%d C=%d CC=%d omsize=%d\n",N0,N,NN,M,MM,L,LL,C,CC,ringp->orderMatrixSize);
                     78:   fprintf(fp,"\n");
                     79:
                     80:   /* print identifier names */
                     81:   if (N-M >0) {
                     82:     fprintf(fp,"Differential variables: ");
                     83:     for (i=M; i<N; i++) fprintf(fp," %4s ",TransX[i]);
                     84:     for (i=M; i<N; i++) fprintf(fp," %4s ",TransD[i]);
                     85:     fprintf(fp,"\n");
                     86:     fprintf(fp,"where ");
                     87:     for (i=M; i<N; i++) {
1.6       takayama   88:       fx = cxx(1,i,1,ringp); fd = cdd(1,i,1,ringp);
                     89:          rf = ppSub(ppMult(fd,fx),ppMult(fx,fd));
                     90:       fprintf(fp," %s %s - %s %s = %s, ",TransD[i],TransX[i],
                     91:               TransX[i],TransD[i],POLYToString(rf,'*',0));
1.1       maekawa    92:     }
                     93:     fprintf(fp,"\n\n");
                     94:   }
                     95:   if (M-L >0) {
                     96:     fprintf(fp,"Difference  variables: ");
                     97:     for (i=L; i<M; i++) fprintf(fp," %4s ",TransX[i]);
                     98:     for (i=L; i<M; i++) fprintf(fp," %4s ",TransD[i]);
                     99:     fprintf(fp,"\n");
                    100:     fprintf(fp,"where ");
                    101:     for (i=L; i<M; i++) {
1.5       takayama  102:       fprintf(fp," %s %s - %s %s = ",TransD[i],TransX[i],
                    103:               TransX[i],TransD[i]);
                    104:       f=ppSub(ppMult(cdd(1,i,1,ringp),cxx(1,i,1,ringp)),
                    105:               ppMult(cxx(1,i,1,ringp),cdd(1,i,1,ringp)));
                    106:       fprintf(fp," %s, ",POLYToString(f,'*',0));
1.1       maekawa   107:     }
                    108:     fprintf(fp,"\n\n");
                    109:   }
                    110:   if (L-C >0) {
                    111:     fprintf(fp,"q-Difference  variables: ");
                    112:     for (i=C; i<L; i++) fprintf(fp," %4s ",TransX[i]);
                    113:     for (i=C; i<L; i++) fprintf(fp," %4s ",TransD[i]);
                    114:     fprintf(fp,"\n");
                    115:     fprintf(fp,"where ");
                    116:     for (i=C; i<L; i++) {
                    117:       fprintf(fp," %s %s = %s %s %s, ",TransD[i],TransX[i],
1.4       takayama  118:               TransX[0],
                    119:               TransX[i],TransD[i]);
1.1       maekawa   120:     }
                    121:     fprintf(fp,"\n\n");
                    122:   }
                    123:   if (C>0) {
                    124:     fprintf(fp,"Commutative  variables: ");
                    125:     for (i=0; i<C; i++) fprintf(fp," %4s ",TransX[i]);
                    126:     for (i=0; i<C; i++) fprintf(fp," %4s ",TransD[i]);
                    127:     fprintf(fp,"\n\n");
                    128:   }
                    129:
                    130:   if (strcmp(F_isSameComponent,"x") == 0) {
                    131:     fprintf(fp,"Integral or summation or graduation variables are : ");
                    132:     for (i=CC; i<C; i++) fprintf(fp," %4s ",TransX[i]);
                    133:     for (i=LL; i<L; i++) fprintf(fp," %4s ",TransX[i]);
                    134:     for (i=MM; i<M; i++) fprintf(fp," %4s ",TransX[i]);
                    135:     for (i=NN; i<N; i++) fprintf(fp," %4s ",TransX[i]);
                    136:     fprintf(fp,"\n");
                    137:   }else if (strcmp(F_isSameComponent,"xd") == 0) {
                    138:     fprintf(fp,"Graduation variables are : ");
                    139:     for (i=CC; i<C; i++) fprintf(fp," %4s ",TransX[i]);
                    140:     for (i=LL; i<L; i++) fprintf(fp," %4s ",TransX[i]);
                    141:     for (i=MM; i<M; i++) fprintf(fp," %4s ",TransX[i]);
                    142:     for (i=NN; i<N; i++) fprintf(fp," %4s ",TransX[i]);
                    143:     for (i=CC; i<C; i++) fprintf(fp," %4s ",TransD[i]);
                    144:     for (i=LL; i<L; i++) fprintf(fp," %4s ",TransD[i]);
                    145:     for (i=MM; i<M; i++) fprintf(fp," %4s ",TransD[i]);
                    146:     for (i=NN; i<N; i++) fprintf(fp," %4s ",TransD[i]);
                    147:     fprintf(fp,"\n");
                    148:   }else {
                    149:     fprintf(fp,"Unknown graduation variable specification.\n\n");
                    150:   }
                    151:   fprintf(fp,"The homogenization variable is : ");
                    152:   fprintf(fp," %4s ",TransD[0]);
                    153:   fprintf(fp,"\n");
                    154:
                    155:
                    156:
                    157:   fprintf(fp,"-------------------------------------------\n");
                    158:   fprintf(fp,"Output order : ");
                    159:   for (i=0; i<2*N; i++) {
                    160:     if (ringp->outputOrder[i] < N) {
                    161:       fprintf(fp,"%s ",TransX[ringp->outputOrder[i]]);
                    162:     }else{
                    163:       fprintf(fp,"%s ",TransD[(ringp->outputOrder[i])-N]);
                    164:     }
                    165:   }
                    166:   fprintf(fp,"\n");
                    167:
                    168:   if (ringp->multiplication == mpMult_poly) {
                    169:     mtype = "poly";
                    170:   }else if  (ringp->multiplication == mpMult_diff) {
                    171:     mtype = "diff";
                    172:   }else if  (ringp->multiplication == mpMult_difference) {
                    173:     mtype = "difference";
                    174:   }else {
                    175:     mtype = "unknown";
                    176:   }
1.16    ! takayama  177:   fprintf(fp,"Multiplication function --%s(%p).\n",
        !           178:           mtype, ringp->multiplication);
1.1       maekawa   179:   if (ringp->schreyer) {
                    180:     fprintf(fp,"schreyer=1, gbListTower=");
                    181:     printObjectList((struct object *)(ringp->gbListTower));
                    182:     fprintf(fp,"\n");
                    183:   }
1.7       takayama  184:   if (ringp->degreeShiftSize) {
1.8       takayama  185:     fprintf(fp,"degreeShift vector (N=%d,Size=%d)= \n[\n",ringp->degreeShiftN,ringp->degreeShiftSize);
1.7       takayama  186:     {
1.8       takayama  187:       int i,j;
                    188:       for (i=0; i<ringp->degreeShiftN; i++) {
                    189:         fprintf(fp," [");
                    190:         for (j=0; j< ringp->degreeShiftSize; j++) {
                    191:           fprintf(fp," %d ",ringp->degreeShift[i*(ringp->degreeShiftSize)+j]);
                    192:         }
                    193:         fprintf(fp,"]\n");
1.7       takayama  194:       }
                    195:     }
                    196:     fprintf(fp,"]\n");
                    197:   }
                    198:   fprintf(fp,"---  weight vectors ---\n");
1.1       maekawa   199:   if (level) printOrder(ringp);
1.13      takayama  200:
                    201:   if (ringp->partialEcart) {
                    202:     fprintf(fp,"---  partialEcartGlobalVarX ---\n");
                    203:     for (i=0; i<ringp->partialEcart; i++) {
                    204:       fprintf(fp," %4s ",TransX[ringp->partialEcartGlobalVarX[i]]);
                    205:     }
                    206:     fprintf(fp,"\n");
                    207:   }
1.1       maekawa   208:
                    209:   if (ringp->next != (struct ring *)NULL) {
                    210:     fprintf(fp,"\n\n-------- The next ring is .... --------------\n");
                    211:     showRing(level,ringp->next);
                    212:   }
                    213: }
                    214:
                    215: /***************************************************************
                    216:    functions related to order
                    217: ******************************************************************/
                    218: #define xtoi(k) ((N-1)-(k))
                    219: #define dtoi(k) ((2*N-1)-(k))
                    220: #define itox(k) ((N-1)-(k))
                    221: #define itod(k) ((2*N-1)-(k))
                    222: #define isX(i) (i<N? 1: 0)
                    223: #define isD(i) (i<N? 0: 1)
                    224: /****************************************************
                    225: i : 0       1         N-1       N           2N-1
                    226: x :x_{N-1} x_{N-2}   x_0
                    227: d :                          D_{N-1}        D_{0}
                    228: if (isX(i))  x_{itox(i)}
                    229: if (isD(i))  D_{itod(i)}
                    230: ******************************************************/
                    231: /* xtoi(0):N-1   xtoi(1):N-2  ....
                    232:    dtoi(0):2N-1  dtoi(1):2N-2 ...
                    233:    itod(N):N-1   dtoi(N-1):N ...
                    234: */
                    235:
                    236: void printOrder(ringp)
1.4       takayama  237:      struct ring *ringp;
1.1       maekawa   238: {
                    239:   int i,j;
                    240:   FILE *fp;
                    241:   char tmp[100];
                    242:   int N,M,L,C,NN,MM,LL,CC;
                    243:   char **TransX,**TransD;
                    244:   int *Order;
                    245:   int P;
                    246:   int omsize;
                    247:   extern char *F_isSameComponent;
                    248:
                    249:   N=ringp->n; M = ringp->m; L = ringp->l; C = ringp->c;
                    250:   NN=ringp->nn; MM = ringp->mm; LL = ringp->ll; CC = ringp->cc;
                    251:   TransX = ringp->x; TransD = ringp->D;
                    252:   Order = ringp->order;
                    253:   P = ringp->p;
                    254:   omsize = ringp->orderMatrixSize;
                    255:
                    256:   fp = stdout;
                    257:
                    258:
                    259:   for (i=0; i<2*N; i++) printf("%4d",i);
                    260:   fprintf(fp,"\n");
                    261:
                    262:   /* print variables names */
                    263:   for (i=0; i<N; i++) {
                    264:     sprintf(tmp,"x%d",N-1-i);
                    265:     fprintf(fp,"%4s",tmp);
                    266:   }
                    267:   for (i=0; i<N; i++) {
                    268:     sprintf(tmp,"D%d",N-1-i);
                    269:     fprintf(fp,"%4s",tmp);
                    270:   }
                    271:   fprintf(fp,"\n");
                    272:
                    273:   /* print identifier names */
                    274:   for (i=0; i<N; i++) fprintf(fp,"%4s",TransX[itox(i)]);
                    275:   for (i=N; i<2*N; i++) fprintf(fp,"%4s",TransD[itod(i)]);
                    276:   fprintf(fp,"\n");
                    277:
                    278:   /* print D: differential     DE: differential, should be eliminated
1.4       takayama  279:      E: difference
                    280:      Q: q-difference
                    281:      C: commutative
1.1       maekawa   282:   */
                    283:   if (strcmp(F_isSameComponent,"x")== 0 || strcmp(F_isSameComponent,"xd")==0) {
                    284:     for (i=0; i<N; i++) {
                    285:       if ((NN<=itox(i)) && (itox(i)<N)) fprintf(fp,"%4s","DE");
                    286:       if ((M<=itox(i)) && (itox(i)<NN)) fprintf(fp,"%4s","D");
                    287:       if ((MM<=itox(i)) && (itox(i)<M)) fprintf(fp,"%4s","EE");
                    288:       if ((L<=itox(i)) && (itox(i)<MM)) fprintf(fp,"%4s","E");
                    289:       if ((LL<=itox(i)) && (itox(i)<L)) fprintf(fp,"%4s","QE");
                    290:       if ((C<=itox(i)) && (itox(i)<LL)) fprintf(fp,"%4s","Q");
                    291:       if ((CC<=itox(i)) && (itox(i)<C)) fprintf(fp,"%4s","CE");
                    292:       if ((0<=itox(i)) && (itox(i)<CC)) fprintf(fp,"%4s","C");
                    293:     }
                    294:   }
                    295:   if (strcmp(F_isSameComponent,"x")==0) {
                    296:     for (i=N; i<2*N; i++) {
                    297:       if ((M<=itod(i)) && (itod(i)<N)) fprintf(fp,"%4s","D");
                    298:       if ((L<=itod(i)) && (itod(i)<M)) fprintf(fp,"%4s","E");
                    299:       if ((C<=itod(i)) && (itod(i)<L)) fprintf(fp,"%4s","Q");
                    300:       if ((0<=itod(i)) && (itod(i)<C)) fprintf(fp,"%4s","C");
                    301:     }
                    302:   }else if (strcmp(F_isSameComponent,"xd")==0) {
                    303:     for (i=N; i<2*N; i++) {
                    304:       if ((NN<=itod(i)) && (itod(i)<N)) fprintf(fp,"%4s","DE");
                    305:       if ((M<=itod(i)) && (itod(i)<NN)) fprintf(fp,"%4s","D");
                    306:       if ((MM<=itod(i)) && (itod(i)<M)) fprintf(fp,"%4s","EE");
                    307:       if ((L<=itod(i)) && (itod(i)<MM)) fprintf(fp,"%4s","E");
                    308:       if ((LL<=itod(i)) && (itod(i)<L)) fprintf(fp,"%4s","QE");
                    309:       if ((C<=itod(i)) && (itod(i)<LL)) fprintf(fp,"%4s","Q");
                    310:       if ((CC<=itod(i)) && (itod(i)<C)) fprintf(fp,"%4s","CE");
                    311:       if ((0<=itod(i)) && (itod(i)<CC)) fprintf(fp,"%4s","C");
                    312:     }
                    313:   } else {
                    314:     fprintf(fp,"Unknown graduation variable type.\n");
                    315:   }
                    316:   fprintf(fp,"\n");
                    317:
                    318:   for (i=0; i< omsize; i++) {
                    319:     for (j=0; j<2*N; j++) {
                    320:       fprintf(fp,"%4d", Order[i*2*N+j]);
                    321:     }
                    322:     fprintf(fp,"\n");
                    323:   }
                    324:   fprintf(fp,"\n");
                    325:
                    326: }
                    327:
                    328: struct object oGetOrderMatrix(struct ring *ringp)
                    329: {
1.14      takayama  330:   struct object rob = OINIT;
                    331:   struct object ob2 = OINIT;
1.1       maekawa   332:   int n,i,j,m;
                    333:   int *om;
                    334:   n = ringp->n;
                    335:   m = ringp->orderMatrixSize;
                    336:   om = ringp->order;
                    337:   if (m<=0) m = 1;
                    338:   rob = newObjectArray(m);
                    339:   for (i=0; i<m; i++) {
                    340:     ob2 = newObjectArray(2*n);
                    341:     for (j=0; j<2*n; j++) {
                    342:       putoa(ob2,j,KpoInteger(om[2*n*i+j]));
                    343:     }
                    344:     putoa(rob,i,ob2);
                    345:   }
                    346:   return(rob);
                    347: }
                    348:
                    349:
                    350: int mmLarger_matrix(ff,gg)
1.4       takayama  351:      POLY ff; POLY gg;
1.1       maekawa   352: {
                    353:   int exp[2*N0]; /* exponents */
                    354:   int i,k;
                    355:   int sum,flag;
                    356:   int *Order;
                    357:   int N;
                    358:   MONOMIAL f,g;
                    359:   struct ring *rp;
                    360:   int in2;
                    361:   int *from, *to;
                    362:   int omsize;
1.7       takayama  363:   int dssize;
1.8       takayama  364:   int dsn;
1.7       takayama  365:   int *degreeShiftVector;
1.1       maekawa   366:
                    367:   if (ff == POLYNULL ) {
                    368:     if (gg == POLYNULL) return( 2 );
                    369:     else return( 0 );
                    370:   }
                    371:   if (gg == POLYNULL) {
                    372:     if (ff == POLYNULL) return( 2 );
                    373:     else return( 1 );
                    374:   }
                    375:   f = ff->m; g=gg->m;
                    376:
                    377:   rp = f->ringp;
                    378:   Order = rp->order;
                    379:   N = rp->n;
                    380:   from = rp->from;
                    381:   to = rp->to;
                    382:   omsize = rp->orderMatrixSize;
1.7       takayama  383:   if (dssize = rp->degreeShiftSize) {
                    384:        degreeShiftVector = rp->degreeShift;  /* Note. 2003.06.26 */
1.8       takayama  385:        dsn = rp->degreeShiftN;
1.7       takayama  386:   }
1.1       maekawa   387:
                    388:   flag = 1;
                    389:   for (i=N-1,k=0; i>=0; i--,k++) {
                    390:     exp[k] = (f->e[i].x) - (g->e[i].x);
                    391:     exp[k+N] = (f->e[i].D) - (g->e[i].D);
                    392:     if ((exp[k] != 0) || (exp[k+N] != 0)) flag =0;
                    393:   }
                    394:   if (flag==1) return(2);
                    395:   /* exp > 0   <--->  f>g
                    396:      exp = 0   <--->  f=g
                    397:      exp < 0   <--->  f<g
                    398:   */
                    399:   for (i=0; i< omsize; i++) {
                    400:     sum = 0; in2 = i*2*N;
                    401:     /* for (k=0; k<2*N; k++) sum += exp[k]*Order[in2+k]; */
                    402:     for (k=from[i]; k<to[i]; k++) sum += exp[k]*Order[in2+k];
1.8       takayama  403:     if (dssize && ( i < dsn)) { /* Note, 2003.06.26 */
1.7       takayama  404:       if ((f->e[N-1].x < dssize) && (f->e[N-1].x >= 0) &&
                    405:           (g->e[N-1].x < dssize) && (g->e[N-1].x >= 0)) {
1.8       takayama  406:         sum += degreeShiftVector[i*dssize+ (f->e[N-1].x)]
                    407:               -degreeShiftVector[i*dssize+ (g->e[N-1].x)];
1.7       takayama  408:       }else{
1.9       takayama  409:         /*warningOrder("Size mismatch in the degree shift vector. It is ignored.");*/
1.7       takayama  410:       }
                    411:     }
1.1       maekawa   412:     if (sum > 0) return(1);
                    413:     if (sum < 0) return(0);
                    414:   }
                    415:   return(2);
                    416: }
                    417:
                    418: /* This should be used in case of q */
                    419: int mmLarger_qmatrix(ff,gg)
1.4       takayama  420:      POLY ff; POLY gg;
1.1       maekawa   421: {
                    422:   int exp[2*N0]; /* exponents */
                    423:   int i,k;
                    424:   int sum,flag;
                    425:   int *Order;
                    426:   int N;
                    427:   MONOMIAL f,g;
                    428:   int omsize;
                    429:
                    430:   if (ff == POLYNULL ) {
                    431:     if (gg == POLYNULL) return( 2 );
                    432:     else return( 0 );
                    433:   }
                    434:   if (gg == POLYNULL) {
                    435:     if (ff == POLYNULL) return( 2 );
                    436:     else return( 1 );
                    437:   }
                    438:   f = ff->m; g = gg->m;
                    439:   Order = f->ringp->order;
                    440:   N = f->ringp->n;
                    441:   omsize = f->ringp->orderMatrixSize;
                    442:
                    443:   flag = 1;
                    444:   for (i=N-1,k=0; i>=0; i--,k++) {
                    445:     exp[k] = (f->e[i].x) - (g->e[i].x);
                    446:     exp[k+N] = (f->e[i].D) - (g->e[i].D);
                    447:     if ((exp[k] != 0) || (exp[k+N] != 0)) flag =0;
                    448:   }
                    449:   if (flag==1) return(2);
                    450:   /* exp > 0   <--->  f>g
                    451:      exp = 0   <--->  f=g
                    452:      exp < 0   <--->  f<g
                    453:   */
                    454:   for (i=0; i< omsize; i++) {
                    455:     sum = 0;
                    456:     /* In case of q, you should do as follows */
                    457:     for (k=0; k<N-1; k++) sum += exp[k]*Order[i*2*N+k]; /* skip k= N-1 -->q */
                    458:     for (k=N; k<2*N-1; k++) sum += exp[k]*Order[i*2*N+k]; /* SKip k= 2*N-1 */
                    459:     if (sum > 0) return(1);
                    460:     else if (sum < 0) return(0);
                    461:   }
                    462:   if (exp[N-1] > 0) return(1);
                    463:   else if (exp[N-1] < 0) return(0);
                    464:   else return(2);
                    465: }
                    466:
                    467: /* x(N-1)>x(N-2)>....>D(N-1)>....>D(0) */
                    468: mmLarger_pureLexicographic(f,g)
1.4       takayama  469:      POLY f;
                    470:      POLY g;
1.1       maekawa   471: {
                    472:   int i,r;
                    473:   int n;
                    474:   MONOMIAL fm,gm;
                    475:   /* Note that this function ignores the order matrix of the given
                    476:      ring. */
                    477:   if (f == POLYNULL ) {
                    478:     if (g == POLYNULL) return( 2 );
                    479:     else return( 0 );
                    480:   }
                    481:   if (g == POLYNULL) {
                    482:     if (f == POLYNULL) return( 2 );
                    483:     else return( 1 );
                    484:   }
                    485:
                    486:
                    487:   fm = f->m; gm = g->m;
                    488:   n = fm->ringp->n;
                    489:   for (i=n-1; i>=0; i--) {
                    490:     r = (fm->e[i].x) - (gm->e[i].x);
                    491:     if (r > 0) return(1);
                    492:     else if (r < 0) return(0);
                    493:     else ;
                    494:   }
                    495:
                    496:   for (i=n-1; i>=0; i--) {
                    497:     r = (fm->e[i].D) - (gm->e[i].D);
                    498:     if (r > 0) return(1);
                    499:     else if (r < 0) return(0);
                    500:     else ;
                    501:   }
                    502:
                    503:   return(2);
                    504:
                    505: }
                    506:
                    507:
                    508: void setFromTo(ringp)
1.4       takayama  509:      struct ring *ringp;
1.1       maekawa   510: {
                    511:   int n;
                    512:   int i,j,oasize;
                    513:   if (ringp->order == (int *)NULL) errorOrder("setFromTo(); no order matrix.");
                    514:   n = (ringp->n)*2;
                    515:   oasize = ringp->orderMatrixSize;
                    516:   ringp->from = (int *)sGC_malloc(sizeof(int)*oasize);
                    517:   ringp->to = (int *)sGC_malloc(sizeof(int)*oasize);
                    518:   if (ringp->from == (int *)NULL  || ringp->to == (int *)NULL) {
                    519:     errorOrder("setFromTo(): No memory.");
                    520:   }
                    521:   for (i=0; i<oasize; i++) {
                    522:     ringp->from[i] = 0; ringp->to[i] = n;
                    523:     for (j=0; j<n; j++) {
                    524:       if (ringp->order[i*n+j] != 0) {
1.4       takayama  525:         ringp->from[i] = j;
                    526:         break;
1.1       maekawa   527:       }
                    528:     }
                    529:     for (j=n-1; j>=0; j--) {
                    530:       if (ringp->order[i*n+j] != 0) {
1.4       takayama  531:         ringp->to[i] = j+1;
                    532:         break;
1.1       maekawa   533:       }
                    534:     }
                    535:   }
                    536: }
                    537:
                    538: /* It ignores h and should be used with mmLarger_tower */
                    539: /* cf. mmLarger_matrix.  h always must be checked at last. */
                    540: static int mmLarger_matrix_schreyer(ff,gg)
1.4       takayama  541:      POLY ff; POLY gg;
1.1       maekawa   542: {
                    543:   int exp[2*N0]; /* exponents */
                    544:   int i,k;
                    545:   int sum,flag;
                    546:   int *Order;
                    547:   int N;
                    548:   MONOMIAL f,g;
                    549:   struct ring *rp;
                    550:   int in2;
                    551:   int *from, *to;
                    552:   int omsize;
                    553:
                    554:   if (ff == POLYNULL ) {
                    555:     if (gg == POLYNULL) return( 2 );
                    556:     else return( 0 );
                    557:   }
                    558:   if (gg == POLYNULL) {
                    559:     if (ff == POLYNULL) return( 2 );
                    560:     else return( 1 );
                    561:   }
                    562:   f = ff->m; g=gg->m;
                    563:
                    564:   rp = f->ringp;
                    565:   Order = rp->order;
                    566:   N = rp->n;
                    567:   from = rp->from;
                    568:   to = rp->to;
                    569:   omsize = rp->orderMatrixSize;
                    570:
                    571:   flag = 1;
                    572:   for (i=N-1,k=0; i>0; i--,k++) {
                    573:     exp[k] = (f->e[i].x) - (g->e[i].x);
                    574:     exp[k+N] = (f->e[i].D) - (g->e[i].D);
                    575:     if ((exp[k] != 0) || (exp[k+N] != 0)) flag =0;
                    576:   }
                    577:   exp[N-1] = (f->e[0].x) - (g->e[0].x);
                    578:   exp[2*N-1] = 0;  /* f->e[0].D - g->e[0].D.  Ignore h! */
                    579:   if ((exp[N-1] != 0) || (exp[2*N-1] != 0)) flag =0;
                    580:
                    581:   if (flag==1) return(2);
                    582:   /* exp > 0   <--->  f>g
                    583:      exp = 0   <--->  f=g
                    584:      exp < 0   <--->  f<g
                    585:   */
                    586:   for (i=0; i< omsize; i++) {
                    587:     sum = 0; in2 = i*2*N;
                    588:     /* for (k=0; k<2*N; k++) sum += exp[k]*Order[in2+k]; */
                    589:     for (k=from[i]; k<to[i]; k++) sum += exp[k]*Order[in2+k];
                    590:     if (sum > 0) return(1);
                    591:     if (sum < 0) return(0);
                    592:   }
                    593:   return(2);
                    594: }
                    595:
                    596: int mmLarger_tower(POLY f,POLY g) {
                    597:   struct object *gbList;
                    598:   int r;
                    599:   if (f == POLYNULL) {
                    600:     if (g == POLYNULL)  return(2);
                    601:     else return(0);
                    602:   }
                    603:   if (g == POLYNULL) {
                    604:     if (f == POLYNULL) return(2);
                    605:     else return(1);
                    606:   }
                    607:   if (!(f->m->ringp->schreyer) || !(g->m->ringp->schreyer))
                    608:     return(mmLarger_matrix(f,g));
1.4       takayama  609:   /* modifiable: mmLarger_qmatrix */
1.1       maekawa   610:   gbList = (struct object *)(g->m->ringp->gbListTower);
                    611:   if (gbList == NULL) return(mmLarger_matrix(f,g));
1.4       takayama  612:   /* modifiable: mmLarger_qmatrix */
1.1       maekawa   613:   if (gbList->tag != Slist) {
                    614:     warningOrder("mmLarger_tower(): gbList must be in Slist.\n");
                    615:     return(1);
                    616:   }
                    617:   if (klength(gbList) ==0) return(mmLarger_matrix(f,g));
1.4       takayama  618:   /* modifiable: mmLarger_qmatrix */
1.1       maekawa   619:
                    620:   r = mmLarger_tower3(f,g,gbList);
                    621:   /* printf("mmLarger_tower3(%s,%s) -->  %d\n",POLYToString(head(f),'*',1),POLYToString(head(g),'*',1),r); */
                    622:   if (r == 2) { /* Now, compare by h */
                    623:     if (f->m->e[0].D > g->m->e[0].D) return(1);
                    624:     else if (f->m->e[0].D < g->m->e[0].D) return(0);
                    625:     else return(2);
                    626:   }else{
                    627:     return(r);
                    628:   }
                    629: }
                    630:
                    631: int mmLarger_tower3(POLY f,POLY g,struct object *gbList)
                    632: { /* gbList is assumed to be Slist */
                    633:   int n,fv,gv,t,r,nn;
                    634:   POLY fm;
                    635:   POLY gm;
1.14      takayama  636:   struct object gb = OINIT;
1.1       maekawa   637:
                    638:   if (f == POLYNULL) {
                    639:     if (g == POLYNULL)  return(2);
                    640:     else return(0);
                    641:   }
                    642:   if (g == POLYNULL) {
                    643:     if (f == POLYNULL) return(2);
                    644:     else return(1);   /* It assumes the zero is the minimum element!! */
                    645:   }
                    646:   n = f->m->ringp->n;
                    647:   nn = f->m->ringp->nn;
                    648:   /* critical and modifiable */  /* m e_u > m e_v <==> m g_u > m g_v */
1.4       takayama  649:   /*                  or equal and u < v */
1.1       maekawa   650:   fv = f->m->e[nn].x ; /* extract component (vector) number of f! */
                    651:   gv = g->m->e[nn].x ;
                    652:   if (fv == gv) { /* They have the same component number. */
                    653:     return(mmLarger_matrix_schreyer(f,g));
                    654:   }
                    655:
                    656:   if (gbList == NULL) return(mmLarger_matrix_schreyer(f,g));
1.4       takayama  657:   /* modifiable: mmLarger_qmatrix */
1.1       maekawa   658:   if (gbList->tag != Slist) {
                    659:     warningOrder("mmLarger_tower(): gbList must be in Slist.\n");
                    660:     return(1);
                    661:   }
                    662:   if (klength(gbList) ==0) return(mmLarger_matrix(f,g));
1.4       takayama  663:   /* modifiable: mmLarger_qmatrix */
1.1       maekawa   664:   gb = car(gbList);  /* each entry must be monomials */
                    665:   if (gb.tag != Sarray) {
                    666:     warningOrder("mmLarger_tower3(): car(gbList) must be an array.\n");
                    667:     return(1);
                    668:   }
                    669:   t = getoaSize(gb);
                    670:   if (t == 0) return(mmLarger_tower3(f,g,cdr(gbList)));
                    671:
                    672:   fm = pmCopy(head(f)); fm->m->e[nn].x = 0; /* f is not modified. */
                    673:   gm = pmCopy(head(g)); gm->m->e[nn].x = 0;
                    674:   if (fv >= t || gv >= t) {
                    675:     warningOrder("mmLarger_tower3(): incompatible input and gbList.\n");
                    676:     printf("Length of gb is %d, f is %s, g is %s\n",t,KPOLYToString(f),
1.4       takayama  677:            KPOLYToString(g));
1.3       takayama  678:     KSexecuteString(" show_ring ");
1.1       maekawa   679:     return(1);
                    680:   }
                    681:   /* mpMult_poly is too expensive to call. @@@*/
                    682:   r = mmLarger_tower3(mpMult_poly(fm,KopPOLY(getoa(gb,fv))),
                    683:                       mpMult_poly(gm,KopPOLY(getoa(gb,gv))),
                    684:                       cdr(gbList));
                    685:   if (r != 2) return(r);
                    686:   else if (fv == gv) return(2);
                    687:   else if (fv > gv) return(0); /* modifiable */
                    688:   else if (fv < gv) return(1); /* modifiable */
                    689: }
1.11      takayama  690:
                    691: static struct object auxPruneZeroRow(struct object ob) {
                    692:   int i,m,size;
1.14      takayama  693:   struct object obt = OINIT;
                    694:   struct object rob = OINIT;
1.11      takayama  695:   m = getoaSize(ob);
                    696:   size=0;
                    697:   for (i=0; i<m; i++) {
                    698:        obt = getoa(ob,i);
                    699:        if (getoaSize(obt) != 0) size++;
                    700:   }
                    701:   if (size == m) return ob;
                    702:   rob = newObjectArray(size);
                    703:   for (i=0, size=0; i<m; i++) {
                    704:        obt = getoa(ob,i);
                    705:        if (getoaSize(obt) != 0) {
                    706:          putoa(rob,size,obt); size++;
                    707:        }
                    708:   }
                    709:   return rob;
                    710: }
1.12      takayama  711: static struct object oRingToOXringStructure_long(struct ring *ringp)
1.10      takayama  712: {
1.14      takayama  713:   struct object rob = OINIT;
                    714:   struct object ob2 = OINIT;
                    715:   struct object obMat = OINIT;
                    716:   struct object obV = OINIT;
                    717:   struct object obShift = OINIT;
                    718:   struct object obt = OINIT;
1.10      takayama  719:   char **TransX; char **TransD;
                    720:   int n,i,j,m,p,nonzero;
                    721:   int *om;
                    722:   n = ringp->n;
                    723:   m = ringp->orderMatrixSize;
                    724:   om = ringp->order;
                    725:   TransX = ringp->x; TransD = ringp->D;
                    726:   if (m<=0) m = 1;
                    727:   /*test: (1). getRing /rr set rr (oxRingStructure) dc  */
                    728:   obMat = newObjectArray(m);
                    729:   for (i=0; i<m; i++) {
                    730:     nonzero = 0;
                    731:     for (j=0; j<2*n; j++) {
                    732:       if (om[2*n*i+j] != 0) nonzero++;
                    733:     }
                    734:     ob2 = newObjectArray(nonzero*2);
                    735:     nonzero=0;
                    736:     for (j=0; j<2*n; j++) {
                    737:       /* fprintf(stderr,"%d, ",nonzero); */
                    738:       if (om[2*n*i+j] != 0) {
                    739:         if (j < n) {
                    740:           putoa(ob2,nonzero,KpoString(TransX[n-1-j])); nonzero++;
                    741:         }else{
                    742:           putoa(ob2,nonzero,KpoString(TransD[n-1-(j-n)])); nonzero++;
                    743:         }
                    744:         putoa(ob2,nonzero,KpoUniversalNumber(newUniversalNumber(om[2*n*i+j]))); nonzero++;
                    745:       }
                    746:     }
                    747:     /* printObject(ob2,0,stderr); fprintf(stderr,".\n"); */
                    748:     putoa(obMat,i,ob2);
                    749:   }
1.11      takayama  750:   obMat = auxPruneZeroRow(obMat);
1.10      takayama  751:   /* printObject(obMat,0,stderr); */
                    752:
                    753:   obV = newObjectArray(2*n);
                    754:   for (i=0; i<n; i++) putoa(obV,i,KpoString(TransX[n-1-i]));
                    755:   for (i=0; i<n; i++) putoa(obV,i+n,KpoString(TransD[n-1-i]));
                    756:   /* printObject(obV,0,stderr); */
                    757:
                    758:   if (ringp->degreeShiftSize) {
                    759:     /*test:
                    760:     [(x) ring_of_differential_operators [[(x)]] weight_vector 0
                    761:       [(weightedHomogenization) 1 (degreeShift) [[1 2 1]]] ] define_ring ;
                    762:      (1). getRing /rr set rr (oxRingStructure) dc message
                    763:     */
                    764:     obShift = newObjectArray(ringp->degreeShiftN);
                    765:     for (i=0; i<ringp->degreeShiftN; i++) {
                    766:       obt = newObjectArray(ringp->degreeShiftSize);
                    767:       for (j=0; j< ringp->degreeShiftSize; j++) {
                    768:         putoa(obt,j,KpoUniversalNumber(newUniversalNumber(ringp->degreeShift[i*(ringp->degreeShiftSize)+j])));
                    769:       }
                    770:       putoa(obShift,i,obt);
                    771:     }
                    772:     /* printObject(obShift,0,stderr); */
                    773:   }
                    774:
                    775:   p = 0;
                    776:   if (ringp->degreeShiftSize) {
                    777:     rob = newObjectArray(3);
                    778:     obt = newObjectArray(2);
                    779:     putoa(obt,0,KpoString("degreeShift"));
                    780:     putoa(obt,1,obShift);
                    781:     putoa(rob,p, obt); p++;
                    782:   }else {
                    783:     rob = newObjectArray(2);
                    784:   }
                    785:
                    786:   obt = newObjectArray(2);
                    787:   putoa(obt,0,KpoString("v"));
                    788:   putoa(obt,1,obV);
                    789:   putoa(rob,p, obt); p++;
                    790:
                    791:   obt = newObjectArray(2);
                    792:   putoa(obt,0,KpoString("order"));
                    793:   putoa(obt,1,obMat);
                    794:   putoa(rob,p, obt); p++;
                    795:
1.12      takayama  796:   return(rob);
                    797: }
                    798: static int auxEffectiveVar(int idx,int n) {
                    799:   int x;
                    800:   if (idx < n) x=1; else x=0;
                    801:   if (x) {
                    802:        if ((idx >= 1) && (idx < n-1)) return 1;
                    803:        else return 0;
                    804:   }else{
                    805:        if ( 1 <= idx-n )  return 1;
                    806:        else return 0;
                    807:   }
                    808: }
                    809: /*test:
                    810:    [(x,y) ring_of_differential_operators [[(Dx) 1 (Dy)  1]]
                    811:     weight_vector 0] define_ring
                    812:     (x). getRing (oxRingStructure) dc ::
                    813:  */
                    814: static struct object oRingToOXringStructure_short(struct ring *ringp)
                    815: {
1.14      takayama  816:   struct object rob = OINIT;
                    817:   struct object ob2 = OINIT;
                    818:   struct object obMat = OINIT;
                    819:   struct object obV = OINIT;
                    820:   struct object obShift = OINIT;
                    821:   struct object obt = OINIT;
1.12      takayama  822:   char **TransX; char **TransD;
                    823:   int n,i,j,m,p,nonzero;
                    824:   int *om;
                    825:   n = ringp->n;
                    826:   m = ringp->orderMatrixSize;
                    827:   om = ringp->order;
                    828:   TransX = ringp->x; TransD = ringp->D;
                    829:   if (m<=0) m = 1;
                    830:   /*test: (1). getRing /rr set rr (oxRingStructure) dc  */
                    831:   obMat = newObjectArray(m);
                    832:   for (i=0; i<m; i++) {
                    833:     nonzero = 0;
                    834:     for (j=0; j<2*n; j++) {
                    835:       if ((om[2*n*i+j] != 0) && auxEffectiveVar(j,n)) nonzero++;
                    836:     }
                    837:     ob2 = newObjectArray(nonzero*2);
                    838:     nonzero=0;
                    839:     for (j=0; j<2*n; j++) {
                    840:       /* fprintf(stderr,"%d, ",nonzero); */
                    841:       if ((om[2*n*i+j] != 0) && auxEffectiveVar(j,n)) {
                    842:         if (j < n) {
                    843:           putoa(ob2,nonzero,KpoString(TransX[n-1-j])); nonzero++;
                    844:         }else{
                    845:           putoa(ob2,nonzero,KpoString(TransD[n-1-(j-n)])); nonzero++;
                    846:         }
                    847:         putoa(ob2,nonzero,KpoUniversalNumber(newUniversalNumber(om[2*n*i+j]))); nonzero++;
                    848:       }
                    849:     }
                    850:     /* printObject(ob2,0,stderr); fprintf(stderr,".\n"); */
                    851:     putoa(obMat,i,ob2);
                    852:   }
                    853:   obMat = auxPruneZeroRow(obMat);
                    854:   /* printObject(obMat,0,stderr); */
                    855:
                    856:   obV = newObjectArray(2*n-3);
                    857:   for (i=0; i<n-2; i++) putoa(obV,i,KpoString(TransX[n-1-i-1]));
                    858:   for (i=0; i<n-1; i++) putoa(obV,i+n-2,KpoString(TransD[n-1-i-1]));
                    859:   /* printObject(obV,0,stderr); */
                    860:
                    861:   if (ringp->degreeShiftSize) {
                    862:     /*test:
                    863:     [(x) ring_of_differential_operators [[(x)]] weight_vector 0
                    864:       [(weightedHomogenization) 1 (degreeShift) [[1 2 1]]] ] define_ring ;
                    865:      (1). getRing /rr set rr (oxRingStructure) dc message
                    866:     */
                    867:     obShift = newObjectArray(ringp->degreeShiftN);
                    868:     for (i=0; i<ringp->degreeShiftN; i++) {
                    869:       obt = newObjectArray(ringp->degreeShiftSize);
                    870:       for (j=0; j< ringp->degreeShiftSize; j++) {
                    871:         putoa(obt,j,KpoUniversalNumber(newUniversalNumber(ringp->degreeShift[i*(ringp->degreeShiftSize)+j])));
                    872:       }
                    873:       putoa(obShift,i,obt);
                    874:     }
                    875:     /* printObject(obShift,0,stderr); */
                    876:   }
                    877:
                    878:   p = 0;
                    879:   if (ringp->degreeShiftSize) {
                    880:     rob = newObjectArray(3);
                    881:     obt = newObjectArray(2);
                    882:     putoa(obt,0,KpoString("degreeShift"));
                    883:     putoa(obt,1,obShift);
                    884:     putoa(rob,p, obt); p++;
                    885:   }else {
                    886:     rob = newObjectArray(2);
                    887:   }
                    888:
                    889:   obt = newObjectArray(2);
                    890:   putoa(obt,0,KpoString("v"));
                    891:   putoa(obt,1,obV);
                    892:   putoa(rob,p, obt); p++;
                    893:
                    894:   obt = newObjectArray(2);
                    895:   putoa(obt,0,KpoString("order"));
                    896:   putoa(obt,1,obMat);
                    897:   putoa(rob,p, obt); p++;
                    898:
                    899:   return(rob);
                    900: }
                    901: struct object oRingToOXringStructure(struct ring *ringp)
                    902: {
1.14      takayama  903:   struct object rob = OINIT;
                    904:   struct object tob = OINIT;
1.12      takayama  905:   rob = newObjectArray(2);
                    906:   tob = oRingToOXringStructure_short(ringp);
                    907:   putoa(rob,0,tob);
                    908:   tob = oRingToOXringStructure_long(ringp);
                    909:   putoa(rob,1,tob);
1.10      takayama  910:   return(rob);
                    911: }
                    912:
1.1       maekawa   913: static void warningOrder(s)
1.4       takayama  914:      char *s;
1.1       maekawa   915: {
                    916:   fprintf(stderr,"Warning in order.c: %s\n",s);
                    917: }
                    918:
                    919: static void errorOrder(s)
1.4       takayama  920:      char *s;
1.1       maekawa   921: {
                    922:   fprintf(stderr,"order.c: %s\n",s);
                    923:   exit(14);
                    924: }
                    925:
                    926:

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