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

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

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