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

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

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