[BACK]Return to dist.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2018 / engine

Annotation of OpenXM_contrib2/asir2018/engine/dist.c, Revision 1.8

1.1       noro        1: /*
                      2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
                      3:  * All rights reserved.
                      4:  *
                      5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
                      6:  * non-exclusive and royalty-free license to use, copy, modify and
                      7:  * redistribute, solely for non-commercial and non-profit purposes, the
                      8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
                      9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
                     10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
                     11:  * third party developer retains all rights, including but not limited to
                     12:  * copyrights, in and to the SOFTWARE.
                     13:  *
                     14:  * (1) FLL does not grant you a license in any way for commercial
                     15:  * purposes. You may use the SOFTWARE only for non-commercial and
                     16:  * non-profit purposes only, such as academic, research and internal
                     17:  * business use.
                     18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
                     19:  * international copyright treaties. If you make copies of the SOFTWARE,
                     20:  * with or without modification, as permitted hereunder, you shall affix
                     21:  * to all such copies of the SOFTWARE the above copyright notice.
                     22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
                     23:  * shall be made on your publication or presentation in any form of the
                     24:  * results obtained by use of the SOFTWARE.
                     25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
                     26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
                     27:  * for such modification or the source code of the modified part of the
                     28:  * SOFTWARE.
                     29:  *
                     30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
                     31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
                     32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
                     33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
                     34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
                     35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
                     36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
                     37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
                     38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
                     39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
                     40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
                     41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
                     42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
                     43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
                     44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
                     45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
                     46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
                     47:  *
1.8     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2018/engine/dist.c,v 1.7 2019/09/06 00:11:59 noro Exp $
1.1       noro       49: */
                     50: #include "ca.h"
                     51:
                     52: #define ORD_REVGRADLEX 0
                     53: #define ORD_GRADLEX 1
                     54: #define ORD_LEX 2
                     55: #define ORD_BREVGRADLEX 3
                     56: #define ORD_BGRADLEX 4
                     57: #define ORD_BLEX 5
                     58: #define ORD_BREVREV 6
                     59: #define ORD_BGRADREV 7
                     60: #define ORD_BLEXREV 8
                     61: #define ORD_ELIM 9
                     62: #define ORD_WEYL_ELIM 10
                     63: #define ORD_HOMO_WW_DRL 11
                     64: #define ORD_DRL_ZIGZAG 12
                     65: #define ORD_HOMO_WW_DRL_ZIGZAG 13
                     66:
                     67: int cmpdl_drl_zigzag(), cmpdl_homo_ww_drl_zigzag();
                     68: int cmpdl_top_weight();
                     69:
                     70: int (*cmpdl)()=cmpdl_revgradlex;
                     71: int (*cmpdl_tie_breaker)();
                     72: int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex};
                     73:
                     74: Obj current_top_weight;
                     75: int current_top_weight_len;
                     76:
                     77: int do_weyl;
                     78:
                     79: int dp_nelim,dp_fcoeffs;
                     80: struct order_spec *dp_current_spec;
                     81: struct modorder_spec *dp_current_modspec;
                     82: int *dp_dl_work;
                     83:
                     84: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr);
                     85: void comm_quod(VL vl,DP p1,DP p2,DP *pr);
                     86: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr);
                     87: void muldc_trunc(VL vl,DP p,Obj c,DL dl,DP *pr);
                     88: int create_order_spec(VL vl,Obj obj,struct order_spec **specp);
                     89: void create_modorder_spec(int id,LIST shift,struct modorder_spec **s);
                     90:
                     91: void order_init()
                     92: {
                     93:   struct order_spec *spec;
                     94:
                     95:   create_order_spec(0,0,&spec);
                     96:   initd(spec);
                     97:   create_modorder_spec(0,0,&dp_current_modspec);
                     98: }
                     99:
                    100: int has_sfcoef_p(Obj f);
                    101:
                    102: int has_sfcoef(DP f)
                    103: {
                    104:   MP t;
                    105:
                    106:   if ( !f )
                    107:     return 0;
                    108:   for ( t = BDY(f); t; t = NEXT(t) )
                    109:     if ( has_sfcoef_p(t->c) )
                    110:       break;
                    111:   return t ? 1 : 0;
                    112: }
                    113:
                    114: int has_sfcoef_p(Obj f)
                    115: {
                    116:   DCP dc;
                    117:
                    118:   if ( !f )
                    119:     return 0;
                    120:   else if ( NUM(f) )
                    121:     return (NID((Num)f) == N_GFS) ? 1 : 0;
                    122:   else if ( POLY(f) ) {
                    123:     for ( dc = DC((P)f); dc; dc = NEXT(dc) )
                    124:       if ( has_sfcoef_p((Obj)COEF(dc)) )
                    125:         return 1;
                    126:     return 0;
                    127:   } else
                    128:     return 0;
                    129: }
                    130:
                    131: extern Obj nd_top_weight;
                    132:
                    133: void reset_top_weight()
                    134: {
                    135:   cmpdl = cmpdl_tie_breaker;
                    136:   cmpdl_tie_breaker = 0;
                    137:   nd_top_weight = 0;
                    138:   current_top_weight = 0;
                    139:   current_top_weight_len = 0;
                    140: }
                    141:
                    142: void initd(struct order_spec *spec)
                    143: {
                    144:   int len,i,k,row;
                    145:   Q **mat;
                    146:
                    147:   switch ( spec->id ) {
                    148:     case 3:
                    149:       cmpdl = cmpdl_composite;
                    150:       dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
                    151:       break;
                    152:     case 2:
                    153:       cmpdl = cmpdl_matrix;
                    154:       dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
                    155:       break;
                    156:     case 1:
                    157:       cmpdl = cmpdl_order_pair;
                    158:       break;
                    159:     default:
                    160:       switch ( spec->ord.simple ) {
                    161:         case ORD_REVGRADLEX:
                    162:           cmpdl = cmpdl_revgradlex; break;
                    163:         case ORD_GRADLEX:
                    164:           cmpdl = cmpdl_gradlex; break;
                    165:         case ORD_BREVGRADLEX:
                    166:           cmpdl = cmpdl_brevgradlex; break;
                    167:         case ORD_BGRADLEX:
                    168:           cmpdl = cmpdl_bgradlex; break;
                    169:         case ORD_BLEX:
                    170:           cmpdl = cmpdl_blex; break;
                    171:         case ORD_BREVREV:
                    172:           cmpdl = cmpdl_brevrev; break;
                    173:         case ORD_BGRADREV:
                    174:           cmpdl = cmpdl_bgradrev; break;
                    175:         case ORD_BLEXREV:
                    176:           cmpdl = cmpdl_blexrev; break;
                    177:         case ORD_ELIM:
                    178:           cmpdl = cmpdl_elim; break;
                    179:         case ORD_WEYL_ELIM:
                    180:           cmpdl = cmpdl_weyl_elim; break;
                    181:         case ORD_HOMO_WW_DRL:
                    182:           cmpdl = cmpdl_homo_ww_drl; break;
                    183:         case ORD_DRL_ZIGZAG:
                    184:           cmpdl = cmpdl_drl_zigzag; break;
                    185:         case ORD_HOMO_WW_DRL_ZIGZAG:
                    186:           cmpdl = cmpdl_homo_ww_drl_zigzag; break;
                    187:         case ORD_LEX: default:
                    188:           cmpdl = cmpdl_lex; break;
                    189:       }
                    190:       break;
                    191:   }
                    192:   if ( current_top_weight ) {
                    193:     cmpdl_tie_breaker = cmpdl;
                    194:     cmpdl = cmpdl_top_weight;
                    195:     if ( OID(current_top_weight) == O_VECT ) {
                    196:        mat = (Q **)&BDY((VECT)current_top_weight);
                    197:        row = 1;
                    198:     } else {
                    199:        mat = (Q **)BDY((MAT)current_top_weight);
                    200:        row = ((MAT)current_top_weight)->row;
                    201:     }
                    202:     for ( k = 0, len = 0; k < row; k++ )
                    203:       for ( i = 0; i < spec->nv; i++ )
                    204:         if ( mat[k][i] )
                    205:           len = MAX(mpz_size(BDY((Z)mat[k][i])),len);
                    206:     current_top_weight_len = len;
                    207:   }
                    208:   dp_current_spec = spec;
                    209: }
                    210:
1.3       noro      211: int dpm_ordtype;
1.1       noro      212:
                    213: void ptod(VL vl,VL dvl,P p,DP *pr)
                    214: {
                    215:   int n,i,j,k;
                    216:   VL tvl;
                    217:   V v;
                    218:   DL d;
                    219:   MP m;
                    220:   DCP dc;
                    221:   DCP *w;
                    222:   DP r,s,t,u;
                    223:   P x,c;
                    224:
                    225:   if ( !p )
                    226:     *pr = 0;
                    227:   else if ( OID(p) > O_P )
                    228:     error("ptod : only polynomials can be converted.");
                    229:   else {
                    230:     for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
                    231:     if ( NUM(p) ) {
                    232:       NEWDL(d,n);
                    233:       NEWMP(m); m->dl = d; C(m) = (Obj)p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
                    234:     } else {
                    235:       for ( i = 0, tvl = dvl, v = VR(p);
                    236:         tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
                    237:       if ( !tvl ) {
                    238:         for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
                    239:         w = (DCP *)ALLOCA(k*sizeof(DCP));
                    240:         for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
                    241:           w[j] = dc;
                    242:
                    243:         for ( j = k-1, s = 0, MKV(v,x); j >= 0; j-- ) {
                    244:           ptod(vl,dvl,COEF(w[j]),&t); pwrp(vl,x,DEG(w[j]),&c);
                    245:           muldc(vl,t,(Obj)c,&r); addd(vl,r,s,&t); s = t;
                    246:         }
                    247:         *pr = s;
                    248:       } else {
                    249:         for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
                    250:         w = (DCP *)ALLOCA(k*sizeof(DCP));
                    251:         for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
                    252:           w[j] = dc;
                    253:
                    254:         for ( j = k-1, s = 0; j >= 0; j-- ) {
                    255:           ptod(vl,dvl,COEF(w[j]),&t);
1.2       noro      256:           NEWDL(d,n); d->d[i] = ZTOS(DEG(w[j]));
1.1       noro      257:           d->td = MUL_WEIGHT(d->d[i],i);
                    258:           NEWMP(m); m->dl = d; C(m) = (Obj)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
                    259:           comm_muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
                    260:         }
                    261:         *pr = s;
                    262:       }
                    263:     }
                    264:   }
                    265: #if 0
                    266:   if ( !dp_fcoeffs && has_sfcoef(*pr) )
                    267:     dp_fcoeffs = N_GFS;
                    268: #endif
                    269: }
                    270:
                    271: void dtop(VL vl,VL dvl,DP p,Obj *pr)
                    272: {
                    273:   int n,i,j,k;
                    274:   DL d;
                    275:   MP m;
                    276:   MP *a;
                    277:   P r;
                    278:   Obj t,w,s,u;
                    279:   Z q;
                    280:   VL tvl;
                    281:
                    282:   if ( !p )
                    283:     *pr = 0;
                    284:   else {
                    285:     for ( k = 0, m = BDY(p); m; m = NEXT(m), k++ );
                    286:     a = (MP *)ALLOCA(k*sizeof(MP));
                    287:     for ( j = 0, m = BDY(p); j < k; m = NEXT(m), j++ )
                    288:       a[j] = m;
                    289:
                    290:     for ( n = p->nv, j = k-1, s = 0; j >= 0; j-- ) {
                    291:       m = a[j];
                    292:       t = C(m);
                    293:       if ( NUM(t) && NID((Num)t) == N_M ) {
                    294:         mptop((P)t,(P *)&u); t = u;
                    295:       }
                    296:       for ( i = 0, d = m->dl, tvl = dvl;
                    297:         i < n; tvl = NEXT(tvl), i++ ) {
1.2       noro      298:         MKV(tvl->v,r); STOZ(d->d[i],q); pwrp(vl,r,q,(P *)&u);
1.1       noro      299:         arf_mul(vl,t,(Obj)u,&w); t = w;
                    300:       }
                    301:       arf_add(vl,s,t,&u); s = u;
                    302:     }
                    303:     *pr = s;
                    304:   }
                    305: }
                    306:
                    307: void nodetod(NODE node,DP *dp)
                    308: {
                    309:   NODE t;
                    310:   int len,i,td;
                    311:   Q e;
                    312:   DL d;
                    313:   MP m;
                    314:   DP u;
                    315:
                    316:   for ( t = node, len = 0; t; t = NEXT(t), len++ );
                    317:   NEWDL(d,len);
                    318:   for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
                    319:     e = (Q)BDY(t);
                    320:     if ( !e )
                    321:       d->d[i] = 0;
                    322:     else if ( !NUM(e) || !RATN(e) || !INT(e) )
                    323:       error("nodetod : invalid input");
                    324:     else {
1.2       noro      325:       d->d[i] = ZTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
1.1       noro      326:     }
                    327:   }
                    328:   d->td = td;
                    329:   NEWMP(m); m->dl = d; C(m) = (Obj)ONE; NEXT(m) = 0;
                    330:   MKDP(len,m,u); u->sugar = td; *dp = u;
                    331: }
                    332:
                    333: void nodetodpm(NODE node,Obj pos,DPM *dp)
                    334: {
                    335:   NODE t;
                    336:   int len,i,td;
                    337:   Q e;
                    338:   DL d;
                    339:   DMM m;
                    340:   DPM u;
                    341:
                    342:   for ( t = node, len = 0; t; t = NEXT(t), len++ );
                    343:   NEWDL(d,len);
                    344:   for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
                    345:     e = (Q)BDY(t);
                    346:     if ( !e )
                    347:       d->d[i] = 0;
                    348:     else if ( !NUM(e) || !RATN(e) || !INT(e) )
                    349:       error("nodetodpm : invalid input");
                    350:     else {
1.2       noro      351:       d->d[i] = ZTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
1.1       noro      352:     }
                    353:   }
                    354:   d->td = td;
1.2       noro      355:   NEWDMM(m); m->dl = d; m->pos = ZTOS((Q)pos); C(m) = (Obj)ONE; NEXT(m) = 0;
1.1       noro      356:   MKDPM(len,m,u); u->sugar = td; *dp = u;
                    357: }
                    358:
                    359: void dtodpm(DP d,int pos,DPM *dp)
                    360: {
                    361:   DMM mr0,mr;
                    362:   MP m;
                    363:
                    364:   if ( !d ) *dp = 0;
                    365:   else {
                    366:     for ( m = BDY(d), mr0 = 0; m; m = NEXT(m) ) {
                    367:       NEXTDMM(mr0,mr);
                    368:       mr->dl = m->dl;
                    369:       mr->pos = pos;
                    370:       C(mr) = C(m);
                    371:     }
                    372:     MKDPM(d->nv,mr0,*dp); (*dp)->sugar = d->sugar;
                    373:   }
                    374: }
                    375:
                    376: int sugard(MP m)
                    377: {
                    378:   int s;
                    379:
                    380:   for ( s = 0; m; m = NEXT(m) )
                    381:     s = MAX(s,m->dl->td);
                    382:   return s;
                    383: }
                    384:
                    385: void addd(VL vl,DP p1,DP p2,DP *pr)
                    386: {
                    387:   int n;
                    388:   MP m1,m2,mr=0,mr0;
                    389:   Obj t;
                    390:   DL d;
                    391:
                    392:   if ( !p1 )
                    393:     *pr = p2;
                    394:   else if ( !p2 )
                    395:     *pr = p1;
                    396:   else {
                    397:     if ( OID(p1) <= O_R ) {
                    398:       n = NV(p2);  NEWDL(d,n);
                    399:       NEWMP(m1); m1->dl = d; C(m1) = (Obj)p1; NEXT(m1) = 0;
                    400:       MKDP(n,m1,p1); (p1)->sugar = 0;
                    401:     }
                    402:     if ( OID(p2) <= O_R ) {
                    403:       n = NV(p1);  NEWDL(d,n);
                    404:       NEWMP(m2); m2->dl = d; C(m2) = (Obj)p2; NEXT(m2) = 0;
                    405:       MKDP(n,m2,p2); (p2)->sugar = 0;
                    406:     }
                    407:     for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
                    408:       switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
                    409:         case 0:
                    410:           arf_add(vl,C(m1),C(m2),&t);
                    411:           if ( t ) {
                    412:             NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
                    413:           }
                    414:           m1 = NEXT(m1); m2 = NEXT(m2); break;
                    415:         case 1:
                    416:           NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
                    417:           m1 = NEXT(m1); break;
                    418:         case -1:
                    419:           NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
                    420:           m2 = NEXT(m2); break;
                    421:       }
                    422:     if ( !mr0 )
                    423:       if ( m1 )
                    424:         mr0 = m1;
                    425:       else if ( m2 )
                    426:         mr0 = m2;
                    427:       else {
                    428:         *pr = 0;
                    429:         return;
                    430:       }
                    431:     else if ( m1 )
                    432:       NEXT(mr) = m1;
                    433:     else if ( m2 )
                    434:       NEXT(mr) = m2;
                    435:     else
                    436:       NEXT(mr) = 0;
                    437:     MKDP(NV(p1),mr0,*pr);
                    438:     if ( *pr )
                    439:       (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                    440:   }
                    441: }
                    442:
                    443: /* for F4 symbolic reduction */
                    444:
                    445: void symb_addd(DP p1,DP p2,DP *pr)
                    446: {
                    447:   int n;
                    448:   MP m1,m2,mr=0,mr0;
                    449:
                    450:   if ( !p1 )
                    451:     *pr = p2;
                    452:   else if ( !p2 )
                    453:     *pr = p1;
                    454:   else {
                    455:     for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
                    456:       NEXTMP(mr0,mr); C(mr) = (Obj)ONE;
                    457:       switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
                    458:         case 0:
                    459:           mr->dl = m1->dl;
                    460:           m1 = NEXT(m1); m2 = NEXT(m2); break;
                    461:         case 1:
                    462:           mr->dl = m1->dl;
                    463:           m1 = NEXT(m1); break;
                    464:         case -1:
                    465:           mr->dl = m2->dl;
                    466:           m2 = NEXT(m2); break;
                    467:       }
                    468:     }
                    469:     if ( !mr0 )
                    470:       if ( m1 )
                    471:         mr0 = m1;
                    472:       else if ( m2 )
                    473:         mr0 = m2;
                    474:       else {
                    475:         *pr = 0;
                    476:         return;
                    477:       }
                    478:     else if ( m1 )
                    479:       NEXT(mr) = m1;
                    480:     else if ( m2 )
                    481:       NEXT(mr) = m2;
                    482:     else
                    483:       NEXT(mr) = 0;
                    484:     MKDP(NV(p1),mr0,*pr);
                    485:     if ( *pr )
                    486:       (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                    487:   }
                    488: }
                    489:
                    490: /*
                    491:  * destructive merge of two list
                    492:  *
                    493:  * p1, p2 : list of DL
                    494:  * return : a merged list
                    495:  */
                    496:
                    497: NODE symb_merge(NODE m1,NODE m2,int n)
                    498: {
                    499:   NODE top=0,prev,cur,m=0,t;
                    500:   DL d1,d2;
                    501:
                    502:   if ( !m1 )
                    503:     return m2;
                    504:   else if ( !m2 )
                    505:     return m1;
                    506:   else {
                    507:     switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
                    508:       case 0:
                    509:         top = m1; m = NEXT(m2);
                    510:         break;
                    511:       case 1:
                    512:         top = m1; m = m2;
                    513:         break;
                    514:       case -1:
                    515:         top = m2; m = m1;
                    516:         break;
                    517:     }
                    518:     prev = top; cur = NEXT(top);
                    519:     /* BDY(prev) > BDY(m) always holds */
                    520:     while ( cur && m ) {
                    521:       d1 = (DL)BDY(cur);
                    522:       d2 = (DL)BDY(m);
                    523: #if 1
                    524:       switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
                    525: #else
                    526:       /* XXX only valid for DRL */
                    527:       if ( d1->td > d2->td )
                    528:         c = 1;
                    529:       else if ( d1->td < d2->td )
                    530:         c = -1;
                    531:       else {
                    532:         for ( i = n-1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
                    533:         if ( i < 0 )
                    534:           c = 0;
                    535:         else if ( d1->d[i] < d2->d[i] )
                    536:           c = 1;
                    537:         else
                    538:           c = -1;
                    539:       }
                    540:       switch ( c ) {
                    541: #endif
                    542:         case 0:
                    543:           m = NEXT(m);
                    544:           prev = cur; cur = NEXT(cur);
                    545:           break;
                    546:         case 1:
                    547:           t = NEXT(cur); NEXT(cur) = m; m = t;
                    548:           prev = cur; cur = NEXT(cur);
                    549:           break;
                    550:         case -1:
                    551:           NEXT(prev) = m; m = cur;
                    552:           prev = NEXT(prev); cur = NEXT(prev);
                    553:           break;
                    554:       }
                    555:     }
                    556:     if ( !cur )
                    557:       NEXT(prev) = m;
                    558:     return top;
                    559:   }
                    560: }
                    561:
                    562: void _adddl(int n,DL d1,DL d2,DL d3)
                    563: {
                    564:   int i;
                    565:
                    566:   d3->td = d1->td+d2->td;
                    567:   for ( i = 0; i < n; i++ )
                    568:     d3->d[i] = d1->d[i]+d2->d[i];
                    569: }
                    570:
1.3       noro      571: void _addtodl(int n,DL d1,DL d2)
                    572: {
                    573:   int i;
                    574:
                    575:   d2->td += d1->td;
                    576:   for ( i = 0; i < n; i++ )
                    577:     d2->d[i] += d1->d[i];
                    578: }
                    579:
                    580: void _copydl(int n,DL d1,DL d2)
                    581: {
                    582:   int i;
                    583:
                    584:   d2->td = d1->td;
                    585:   for ( i = 0; i < n; i++ )
                    586:     d2->d[i] = d1->d[i];
                    587: }
                    588:
                    589: int _eqdl(int n,DL d1,DL d2)
                    590: {
                    591:   int i;
                    592:
                    593:   if ( d2->td != d1->td ) return 0;
                    594:   for ( i = 0; i < n; i++ )
                    595:     if ( d2->d[i] != d1->d[i] ) return 0;
                    596:   return 1;
                    597: }
                    598:
1.1       noro      599: /* m1 <- m1 U dl*f, destructive */
                    600:
                    601: NODE mul_dllist(DL dl,DP f);
                    602:
                    603: NODE symb_mul_merge(NODE m1,DL dl,DP f,int n)
                    604: {
                    605:   NODE top,prev,cur,n1;
                    606:   DP g;
                    607:   DL t,s;
                    608:   MP m;
                    609:
                    610:   if ( !m1 )
                    611:     return mul_dllist(dl,f);
                    612:   else if ( !f )
                    613:     return m1;
                    614:   else {
                    615:     m = BDY(f);
                    616:     NEWDL_NOINIT(t,n);
                    617:     _adddl(n,m->dl,dl,t);
                    618:     top = m1; prev = 0; cur = m1;
                    619:     while ( m ) {
                    620:       switch ( (*cmpdl)(n,(DL)BDY(cur),t) ) {
                    621:         case 0:
                    622:           prev = cur; cur = NEXT(cur);
                    623:           if ( !cur ) {
                    624:             MKDP(n,m,g);
                    625:             NEXT(prev) = mul_dllist(dl,g);
                    626:             return top;
                    627:           }
                    628:           m = NEXT(m);
                    629:           if ( m ) _adddl(n,m->dl,dl,t);
                    630:           break;
                    631:         case 1:
                    632:           prev = cur; cur = NEXT(cur);
                    633:           if ( !cur ) {
                    634:             MKDP(n,m,g);
                    635:             NEXT(prev) = mul_dllist(dl,g);
                    636:             return top;
                    637:           }
                    638:           break;
                    639:         case -1:
                    640:           NEWDL_NOINIT(s,n);
                    641:           s->td = t->td;
                    642:           bcopy(t->d,s->d,n*sizeof(int));
                    643:           NEWNODE(n1);
                    644:           n1->body = (pointer)s;
                    645:           NEXT(n1) = cur;
                    646:           if ( !prev ) {
                    647:             top = n1; cur = n1;
                    648:           } else {
                    649:             NEXT(prev) = n1; prev = n1;
                    650:           }
                    651:           m = NEXT(m);
                    652:           if ( m ) _adddl(n,m->dl,dl,t);
                    653:           break;
                    654:       }
                    655:     }
                    656:     return top;
                    657:   }
                    658: }
                    659:
                    660: DLBUCKET symb_merge_bucket(DLBUCKET m1,DLBUCKET m2,int n)
                    661: {
                    662:   DLBUCKET top,prev,cur,m,t;
                    663:
                    664:   if ( !m1 )
                    665:     return m2;
                    666:   else if ( !m2 )
                    667:     return m1;
                    668:   else {
                    669:     if ( m1->td == m2->td ) {
                    670:       top = m1;
                    671:       BDY(top) = symb_merge(BDY(top),BDY(m2),n);
                    672:       m = NEXT(m2);
                    673:     } else if ( m1->td > m2->td ) {
                    674:       top = m1; m = m2;
                    675:     } else {
                    676:       top = m2; m = m1;
                    677:     }
                    678:     prev = top; cur = NEXT(top);
                    679:     /* prev->td > m->td always holds */
                    680:     while ( cur && m ) {
                    681:       if ( cur->td == m->td ) {
                    682:         BDY(cur) = symb_merge(BDY(cur),BDY(m),n);
                    683:         m = NEXT(m);
                    684:         prev = cur; cur = NEXT(cur);
                    685:       } else if ( cur->td > m->td ) {
                    686:         t = NEXT(cur); NEXT(cur) = m; m = t;
                    687:         prev = cur; cur = NEXT(cur);
                    688:       } else {
                    689:         NEXT(prev) = m; m = cur;
                    690:         prev = NEXT(prev); cur = NEXT(prev);
                    691:       }
                    692:     }
                    693:     if ( !cur )
                    694:       NEXT(prev) = m;
                    695:     return top;
                    696:   }
                    697: }
                    698:
                    699: void subd(VL vl,DP p1,DP p2,DP *pr)
                    700: {
                    701:   DP t;
                    702:
                    703:   if ( !p2 )
                    704:     *pr = p1;
                    705:   else {
                    706:     chsgnd(p2,&t); addd(vl,p1,t,pr);
                    707:   }
                    708: }
                    709:
                    710: void chsgnd(DP p,DP *pr)
                    711: {
                    712:   MP m,mr=0,mr0;
                    713:   Obj r;
                    714:
                    715:   if ( !p )
                    716:     *pr = 0;
                    717:   else if ( OID(p) <= O_R ) {
                    718:     arf_chsgn((Obj)p,&r); *pr = (DP)r;
                    719:   } else {
                    720:     for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    721:       NEXTMP(mr0,mr); arf_chsgn(C(m),&C(mr)); mr->dl = m->dl;
                    722:     }
                    723:     NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    724:     if ( *pr )
                    725:       (*pr)->sugar = p->sugar;
                    726:   }
                    727: }
                    728:
                    729: void muld(VL vl,DP p1,DP p2,DP *pr)
                    730: {
                    731:   if ( ! do_weyl )
                    732:     comm_muld(vl,p1,p2,pr);
                    733:   else
                    734:     weyl_muld(vl,p1,p2,pr);
                    735: }
                    736:
                    737: void comm_muld(VL vl,DP p1,DP p2,DP *pr)
                    738: {
                    739:   MP m;
                    740:   DP s,t,u;
                    741:   int i,l,l1;
                    742:   static MP *w;
                    743:   static int wlen;
                    744:
                    745:   if ( !p1 || !p2 )
                    746:     *pr = 0;
                    747:   else if ( OID(p1) != O_DP )
                    748:     muldc(vl,p2,(Obj)p1,pr);
                    749:   else if ( OID(p2) != O_DP )
                    750:     muldc(vl,p1,(Obj)p2,pr);
                    751:   else {
                    752:     for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
                    753:     for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
                    754:     if ( l1 < l ) {
                    755:       t = p1; p1 = p2; p2 = t;
                    756:       l = l1;
                    757:     }
                    758:     if ( l > wlen ) {
                    759:       if ( w ) GCFREE(w);
                    760:       w = (MP *)MALLOC(l*sizeof(MP));
                    761:       wlen = l;
                    762:     }
                    763:     for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
                    764:       w[i] = m;
                    765:     for ( s = 0, i = l-1; i >= 0; i-- ) {
                    766:       muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
                    767:     }
                    768:     bzero(w,l*sizeof(MP));
                    769:     *pr = s;
                    770:   }
                    771: }
                    772:
                    773: /* discard terms which is not a multiple of dl */
                    774:
                    775: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr)
                    776: {
                    777:   MP m;
                    778:   DP s,t,u;
                    779:   int i,l,l1;
                    780:   static MP *w;
                    781:   static int wlen;
                    782:
                    783:   if ( !p1 || !p2 )
                    784:     *pr = 0;
                    785:   else if ( OID(p1) != O_DP )
                    786:     muldc_trunc(vl,p2,(Obj)p1,dl,pr);
                    787:   else if ( OID(p2) != O_DP )
                    788:     muldc_trunc(vl,p1,(Obj)p2,dl,pr);
                    789:   else {
                    790:     for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
                    791:     for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
                    792:     if ( l1 < l ) {
                    793:       t = p1; p1 = p2; p2 = t;
                    794:       l = l1;
                    795:     }
                    796:     if ( l > wlen ) {
                    797:       if ( w ) GCFREE(w);
                    798:       w = (MP *)MALLOC(l*sizeof(MP));
                    799:       wlen = l;
                    800:     }
                    801:     for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
                    802:       w[i] = m;
                    803:     for ( s = 0, i = l-1; i >= 0; i-- ) {
                    804:       muldm_trunc(vl,p1,w[i],dl,&t); addd(vl,s,t,&u); s = u;
                    805:     }
                    806:     bzero(w,l*sizeof(MP));
                    807:     *pr = s;
                    808:   }
                    809: }
                    810:
                    811: void comm_quod(VL vl,DP p1,DP p2,DP *pr)
                    812: {
                    813:   MP m=0,m0;
                    814:   DP s,t;
                    815:   int i,n,sugar;
                    816:   DL d1,d2,d;
                    817:   Q a,b;
                    818:
                    819:   if ( !p2 )
                    820:     error("comm_quod : invalid input");
                    821:   if ( !p1 )
                    822:     *pr = 0;
                    823:   else {
                    824:     n = NV(p1);
                    825:     d2 = BDY(p2)->dl;
                    826:     m0 = 0;
                    827:     sugar = p1->sugar;
                    828:     while ( p1 ) {
                    829:       d1 = BDY(p1)->dl;
                    830:       NEWDL(d,n);
                    831:       d->td = d1->td - d2->td;
                    832:       for ( i = 0; i < n; i++ )
                    833:         d->d[i] = d1->d[i]-d2->d[i];
                    834:       NEXTMP(m0,m);
                    835:       m->dl = d;
                    836:       divq((Q)BDY(p1)->c,(Q)BDY(p2)->c,&a); chsgnq(a,&b);
                    837:       C(m) = (Obj)b;
                    838:       muldm_trunc(vl,p2,m,d2,&t);
                    839:       addd(vl,p1,t,&s); p1 = s;
                    840:       C(m) = (Obj)a;
                    841:     }
                    842:     if ( m0 ) {
                    843:       NEXT(m) = 0; MKDP(n,m0,*pr);
                    844:     } else
                    845:       *pr = 0;
                    846:     /* XXX */
                    847:     if ( *pr )
                    848:       (*pr)->sugar = sugar - d2->td;
                    849:   }
                    850: }
                    851:
                    852: void muldm(VL vl,DP p,MP m0,DP *pr)
                    853: {
                    854:   MP m,mr=0,mr0;
                    855:   Obj c;
                    856:   DL d;
                    857:   int n;
                    858:
                    859:   if ( !p )
                    860:     *pr = 0;
                    861:   else {
                    862:     for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
                    863:       m; m = NEXT(m) ) {
                    864:       NEXTMP(mr0,mr);
                    865:       if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                    866:         mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                    867:       else
                    868:         arf_mul(vl,C(m),c,&C(mr));
                    869:       adddl(n,m->dl,d,&mr->dl);
                    870:     }
                    871:     NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    872:     if ( *pr )
                    873:       (*pr)->sugar = p->sugar + m0->dl->td;
                    874:   }
                    875: }
                    876:
                    877: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr)
                    878: {
                    879:   MP m,mr=0,mr0;
                    880:   Obj c;
                    881:   DL d,tdl;
                    882:   int n,i;
                    883:
                    884:   if ( !p )
                    885:     *pr = 0;
                    886:   else {
                    887:     n = NV(p);
                    888:     NEWDL(tdl,n);
                    889:     for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl;
                    890:       m; m = NEXT(m) ) {
                    891:       _adddl(n,m->dl,d,tdl);
                    892:       for ( i = 0; i < n; i++ )
                    893:         if ( tdl->d[i] < dl->d[i] )
                    894:           break;
                    895:       if ( i < n )
                    896:         continue;
                    897:       NEXTMP(mr0,mr);
                    898:       mr->dl = tdl;
                    899:       NEWDL(tdl,n);
                    900:       if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                    901:         mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                    902:       else
                    903:         arf_mul(vl,C(m),(Obj)c,&C(mr));
                    904:     }
                    905:     if ( mr0 ) {
                    906:       NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    907:     } else
                    908:       *pr = 0;
                    909:     if ( *pr )
                    910:       (*pr)->sugar = p->sugar + m0->dl->td;
                    911:   }
                    912: }
                    913:
                    914: void weyl_muld(VL vl,DP p1,DP p2,DP *pr)
                    915: {
                    916:   MP m;
                    917:   DP s,t,u;
                    918:   int i,l;
                    919:   static MP *w;
                    920:   static int wlen;
                    921:
                    922:   if ( !p1 || !p2 )
                    923:     *pr = 0;
                    924:   else if ( OID(p1) != O_DP )
                    925:     muldc(vl,p2,(Obj)p1,pr);
                    926:   else if ( OID(p2) != O_DP )
                    927:     muldc(vl,p1,(Obj)p2,pr);
                    928:   else {
                    929:     for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
                    930:     if ( l > wlen ) {
                    931:       if ( w ) GCFREE(w);
                    932:       w = (MP *)MALLOC(l*sizeof(MP));
                    933:       wlen = l;
                    934:     }
                    935:     for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
                    936:       w[i] = m;
                    937:     for ( s = 0, i = l-1; i >= 0; i-- ) {
                    938:       weyl_muldm(vl,w[i],p2,&t); addd(vl,s,t,&u); s = u;
                    939:     }
                    940:     bzero(w,l*sizeof(MP));
                    941:     *pr = s;
                    942:   }
                    943: }
                    944:
                    945: void actm(VL vl,int nv,MP m1,MP m2,DP *pr)
                    946: {
                    947:   DL d1,d2,d;
                    948:   int n2,i,j,k;
                    949:   Z jq,c,c1;
                    950:   MP m;
                    951:   Obj t;
                    952:
                    953:   d1 = m1->dl;
                    954:   d2 = m2->dl;
                    955:   for ( i = 0; i < nv; i++ )
                    956:     if ( d1->d[i] > d2->d[i] ) {
                    957:       *pr = 0; return;
                    958:     }
                    959:   NEWDL(d,nv);
                    960:   c = ONE;
                    961:   for ( i = 0; i < nv; i++ ) {
                    962:     for ( j = d2->d[i], k = d1->d[i]; k > 0; k--, j-- ) {
1.2       noro      963:       STOZ(j,jq); mulz(c,jq,&c1); c = c1;
1.1       noro      964:     }
                    965:     d->d[i] = d2->d[i]-d1->d[i];
                    966:   }
                    967:   arf_mul(vl,C(m1),C(m2),&t);
                    968:   NEWMP(m);
                    969:   arf_mul(vl,(Obj)c,t,&C(m));
                    970:   m->dl = d;
                    971:   MKDP(nv,m,*pr);
                    972: }
                    973:
                    974: void weyl_actd(VL vl,DP p1,DP p2,DP *pr)
                    975: {
                    976:   int n;
                    977:   MP m1,m2;
                    978:   DP d,r,s;
                    979:
                    980:   if ( !p1 || !p2 ) *pr = 0;
                    981:   else {
                    982:     n = NV(p1);
                    983:     r = 0;
                    984:     for ( m1 = BDY(p1); m1; m1 = NEXT(m1) )
                    985:       for ( m2 = BDY(p2); m2; m2 = NEXT(m2) ) {
                    986:         actm(vl,n,m1,m2,&d);
                    987:         addd(vl,r,d,&s); r = s;
                    988:       }
                    989:     *pr = r;
                    990:   }
                    991: }
                    992:
                    993: /* monomial * polynomial */
                    994:
                    995: void weyl_muldm(VL vl,MP m0,DP p,DP *pr)
                    996: {
                    997:   DP r,t,t1;
                    998:   MP m;
                    999:   DL d0;
                   1000:   int n,n2,l,i,j,tlen;
                   1001:   static MP *w,*psum;
                   1002:   static struct cdl *tab;
                   1003:   static int wlen;
                   1004:   static int rtlen;
                   1005:
                   1006:   if ( !p )
                   1007:     *pr = 0;
                   1008:   else {
                   1009:     for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
                   1010:     if ( l > wlen ) {
                   1011:       if ( w ) GCFREE(w);
                   1012:       w = (MP *)MALLOC(l*sizeof(MP));
                   1013:       wlen = l;
                   1014:     }
                   1015:     for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
                   1016:       w[i] = m;
                   1017:
                   1018:     n = NV(p); n2 = n>>1;
                   1019:     d0 = m0->dl;
                   1020:     for ( i = 0, tlen = 1; i < n2; i++ )
                   1021:       tlen *= d0->d[n2+i]+1;
                   1022:     if ( tlen > rtlen ) {
                   1023:       if ( tab ) GCFREE(tab);
                   1024:       if ( psum ) GCFREE(psum);
                   1025:       rtlen = tlen;
                   1026:       tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
                   1027:       psum = (MP *)MALLOC(rtlen*sizeof(MP));
                   1028:     }
                   1029:     bzero(psum,tlen*sizeof(MP));
                   1030:     for ( i = l-1; i >= 0; i-- ) {
                   1031:       bzero(tab,tlen*sizeof(struct cdl));
                   1032:       weyl_mulmm(vl,m0,w[i],n,tab,tlen);
                   1033:       for ( j = 0; j < tlen; j++ ) {
                   1034:         if ( tab[j].c ) {
                   1035:           NEWMP(m); m->dl = tab[j].d; C(m) = (Obj)tab[j].c; NEXT(m) = psum[j];
                   1036:           psum[j] = m;
                   1037:         }
                   1038:       }
                   1039:     }
                   1040:     for ( j = tlen-1, r = 0; j >= 0; j-- )
                   1041:       if ( psum[j] ) {
                   1042:         MKDP(n,psum[j],t); addd(vl,r,t,&t1); r = t1;
                   1043:       }
                   1044:     if ( r )
                   1045:       r->sugar = p->sugar + m0->dl->td;
                   1046:     *pr = r;
                   1047:   }
                   1048: }
                   1049:
                   1050: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
                   1051: /* rtab : array of length (e0+1)*(e1+1)*... */
                   1052:
                   1053: void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
                   1054: {
                   1055:   Obj c,c0,c1;
                   1056:   DL d,d0,d1,dt;
                   1057:   int i,j,a,b,k,l,n2,s,min,curlen;
                   1058:   struct cdl *p;
                   1059:   static Z *ctab;
                   1060:   static struct cdl *tab;
                   1061:   static int tablen;
                   1062:   static struct cdl *tmptab;
                   1063:   static int tmptablen;
                   1064:
                   1065:
                   1066:   if ( !m0 || !m1 ) {
                   1067:     rtab[0].c = 0;
                   1068:     rtab[0].d = 0;
                   1069:     return;
                   1070:   }
                   1071:   c0 = C(m0); c1 = C(m1);
                   1072:   arf_mul(vl,c0,c1,&c);
                   1073:   d0 = m0->dl; d1 = m1->dl;
                   1074:   n2 = n>>1;
                   1075:   curlen = 1;
                   1076:   NEWDL(d,n);
                   1077:   if ( n & 1 )
                   1078:     /* offset of h-degree */
                   1079:      d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
                   1080:   else
                   1081:     d->td = 0;
                   1082:   rtab[0].c = c;
                   1083:   rtab[0].d = d;
                   1084:
                   1085:   if ( rtablen > tmptablen ) {
                   1086:     if ( tmptab ) GCFREE(tmptab);
                   1087:     tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
                   1088:     tmptablen = rtablen;
                   1089:   }
                   1090:   for ( i = 0; i < n2; i++ ) {
                   1091:     a = d0->d[i]; b = d1->d[n2+i];
                   1092:     k = d0->d[n2+i]; l = d1->d[i];
                   1093:
                   1094:     /* degree of xi^a*(Di^k*xi^l)*Di^b */
                   1095:     a += l;
                   1096:     b += k;
                   1097:     s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
                   1098:
                   1099:     if ( !k || !l ) {
                   1100:       for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
                   1101:         if ( p->c ) {
                   1102:           dt = p->d;
                   1103:           dt->d[i] = a;
                   1104:           dt->d[n2+i] = b;
                   1105:           dt->td += s;
                   1106:         }
                   1107:       }
                   1108:       curlen *= k+1;
                   1109:       continue;
                   1110:     }
                   1111:     if ( k+1 > tablen ) {
                   1112:       if ( tab ) GCFREE(tab);
                   1113:       if ( ctab ) GCFREE(ctab);
                   1114:       tablen = k+1;
                   1115:       tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
                   1116:       ctab = (Z *)MALLOC(tablen*sizeof(Q));
                   1117:     }
                   1118:     /* compute xi^a*(Di^k*xi^l)*Di^b */
                   1119:     min = MIN(k,l);
                   1120:     mkwc(k,l,ctab);
                   1121:     bzero(tab,(k+1)*sizeof(struct cdl));
                   1122:     if ( n & 1 )
                   1123:       for ( j = 0; j <= min; j++ ) {
                   1124:         NEWDL(d,n);
                   1125:         d->d[i] = a-j; d->d[n2+i] = b-j;
                   1126:         d->td = s;
                   1127:         d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
                   1128:         tab[j].d = d;
                   1129:         tab[j].c = (Obj)ctab[j];
                   1130:       }
                   1131:     else
                   1132:       for ( j = 0; j <= min; j++ ) {
                   1133:         NEWDL(d,n);
                   1134:         d->d[i] = a-j; d->d[n2+i] = b-j;
                   1135:         d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
                   1136:         tab[j].d = d;
                   1137:         tab[j].c = (Obj)ctab[j];
                   1138:       }
                   1139:     bzero(ctab,(min+1)*sizeof(Q));
                   1140:     comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
                   1141:     curlen *= k+1;
                   1142:     bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
                   1143:   }
                   1144: }
                   1145:
                   1146: /* direct product of two cdl tables
                   1147:   rt[] = [
                   1148:     t[0]*t1[0],...,t[n-1]*t1[0],
                   1149:     t[0]*t1[1],...,t[n-1]*t1[1],
                   1150:     ...
                   1151:     t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
                   1152:   ]
                   1153: */
                   1154:
                   1155: void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
                   1156: {
                   1157:   int i,j;
                   1158:   struct cdl *p;
                   1159:   Obj c;
                   1160:   DL d;
                   1161:
                   1162:   bzero(rt,n*n1*sizeof(struct cdl));
                   1163:   for ( j = 0, p = rt; j < n1; j++ ) {
                   1164:     c = (Obj)t1[j].c;
                   1165:     d = t1[j].d;
                   1166:     if ( !c )
                   1167:       break;
                   1168:     for ( i = 0; i < n; i++, p++ ) {
                   1169:       if ( t[i].c ) {
                   1170:         arf_mul(vl,(Obj)t[i].c,c,(Obj *)&p->c);
                   1171:         adddl(nv,t[i].d,d,&p->d);
                   1172:       }
                   1173:     }
                   1174:   }
                   1175: }
                   1176:
                   1177: void muldc(VL vl,DP p,Obj c,DP *pr)
                   1178: {
                   1179:   MP m,mr=0,mr0;
                   1180:
                   1181:   if ( !p || !c )
                   1182:     *pr = 0;
                   1183:   else if ( NUM(c) && UNIQ((Q)c) )
                   1184:     *pr = p;
                   1185:   else if ( NUM(c) && MUNIQ((Q)c) )
                   1186:     chsgnd(p,pr);
                   1187:   else {
                   1188:     for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1189:       NEXTMP(mr0,mr);
                   1190:       if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                   1191:         mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                   1192:       else
                   1193:         arf_mul(vl,C(m),c,&C(mr));
                   1194:       mr->dl = m->dl;
                   1195:     }
                   1196:     NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                   1197:     if ( *pr )
                   1198:       (*pr)->sugar = p->sugar;
                   1199:   }
                   1200: }
                   1201:
                   1202: void divdc(VL vl,DP p,Obj c,DP *pr)
                   1203: {
                   1204:   Obj inv;
                   1205:
                   1206:   arf_div(vl,(Obj)ONE,c,&inv);
                   1207:   muld(vl,p,(DP)inv,pr);
                   1208: }
                   1209:
                   1210: void muldc_trunc(VL vl,DP p,Obj c,DL dl,DP *pr)
                   1211: {
                   1212:   MP m,mr=0,mr0;
                   1213:   DL mdl;
                   1214:   int i,n;
                   1215:
                   1216:   if ( !p || !c ) {
                   1217:     *pr = 0; return;
                   1218:   }
                   1219:   n = NV(p);
                   1220:   for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1221:     mdl = m->dl;
                   1222:     for ( i = 0; i < n; i++ )
                   1223:       if ( mdl->d[i] < dl->d[i] )
                   1224:         break;
                   1225:     if ( i < n )
                   1226:       break;
                   1227:     NEXTMP(mr0,mr);
                   1228:     if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                   1229:       mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                   1230:     else
                   1231:       arf_mul(vl,C(m),c,&C(mr));
                   1232:     mr->dl = m->dl;
                   1233:   }
                   1234:   NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                   1235:   if ( *pr )
                   1236:     (*pr)->sugar = p->sugar;
                   1237: }
                   1238:
                   1239: void divsdc(VL vl,DP p,P c,DP *pr)
                   1240: {
                   1241:   MP m,mr=0,mr0;
                   1242:
                   1243:   if ( !c )
                   1244:     error("disvsdc : division by 0");
                   1245:   else if ( !p )
                   1246:     *pr = 0;
                   1247:   else if ( OID(c) > O_P )
                   1248:     error("divsdc : invalid argument");
                   1249:   else {
                   1250:     for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1251:       NEXTMP(mr0,mr); divsp(vl,(P)C(m),c,(P *)&C(mr)); mr->dl = m->dl;
                   1252:     }
                   1253:     NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                   1254:     if ( *pr )
                   1255:       (*pr)->sugar = p->sugar;
                   1256:   }
                   1257: }
                   1258:
                   1259: void adddl(int n,DL d1,DL d2,DL *dr)
                   1260: {
                   1261:   DL dt;
                   1262:   int i;
                   1263:
                   1264:   *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
                   1265:   dt->td = d1->td + d2->td;
                   1266:   for ( i = 0; i < n; i++ )
                   1267:     dt->d[i] = d1->d[i]+d2->d[i];
                   1268: }
                   1269:
                   1270: /* d1 += d2 */
                   1271:
                   1272: void adddl_destructive(int n,DL d1,DL d2)
                   1273: {
                   1274:   int i;
                   1275:
                   1276:   d1->td += d2->td;
                   1277:   for ( i = 0; i < n; i++ )
                   1278:     d1->d[i] += d2->d[i];
                   1279: }
                   1280:
                   1281: int compd(VL vl,DP p1,DP p2)
                   1282: {
                   1283:   int n,t;
                   1284:   MP m1,m2;
                   1285:
                   1286:   if ( !p1 )
                   1287:     return p2 ? -1 : 0;
                   1288:   else if ( !p2 )
                   1289:     return 1;
                   1290:   else if ( NV(p1) != NV(p2) ) {
                   1291:     error("compd : size mismatch");
                   1292:     return 0; /* XXX */
                   1293:   } else {
                   1294:     for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
                   1295:       m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
                   1296:       if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
                   1297:         (t = arf_comp(vl,C(m1),C(m2)) ) )
                   1298:         return t;
                   1299:     if ( m1 )
                   1300:       return 1;
                   1301:     else if ( m2 )
                   1302:       return -1;
                   1303:     else
                   1304:       return 0;
                   1305:   }
                   1306: }
                   1307:
                   1308: int cmpdl_lex(int n,DL d1,DL d2)
                   1309: {
                   1310:   int i;
                   1311:
                   1312:   for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
                   1313:   return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
                   1314: }
                   1315:
                   1316: int cmpdl_revlex(int n,DL d1,DL d2)
                   1317: {
                   1318:   int i;
                   1319:
                   1320:   for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
                   1321:   return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                   1322: }
                   1323:
                   1324: int cmpdl_gradlex(int n,DL d1,DL d2)
                   1325: {
                   1326:   if ( d1->td > d2->td )
                   1327:     return 1;
                   1328:   else if ( d1->td < d2->td )
                   1329:     return -1;
                   1330:   else
                   1331:     return cmpdl_lex(n,d1,d2);
                   1332: }
                   1333:
                   1334: int cmpdl_revgradlex(int n,DL d1,DL d2)
                   1335: {
                   1336:   register int i,c;
                   1337:   register int *p1,*p2;
                   1338:
                   1339:   if ( d1->td > d2->td )
                   1340:     return 1;
                   1341:   else if ( d1->td < d2->td )
                   1342:     return -1;
                   1343:   else {
                   1344:     i = n-1;
                   1345:     p1 = d1->d+n-1;
                   1346:     p2 = d2->d+n-1;
                   1347:     while ( i >= 7 ) {
                   1348:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1349:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1350:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1351:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1352:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1353:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1354:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1355:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1356:       i -= 8;
                   1357:     }
                   1358:     switch ( i ) {
                   1359:       case 6:
                   1360:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1361:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1362:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1363:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1364:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1365:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1366:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1367:         return 0;
                   1368:       case 5:
                   1369:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1370:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1371:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1372:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1373:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1374:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1375:         return 0;
                   1376:       case 4:
                   1377:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1378:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1379:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1380:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1381:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1382:         return 0;
                   1383:       case 3:
                   1384:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1385:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1386:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1387:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1388:         return 0;
                   1389:       case 2:
                   1390:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1391:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1392:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1393:         return 0;
                   1394:       case 1:
                   1395:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1396:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1397:         return 0;
                   1398:       case 0:
                   1399:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1400:         return 0;
                   1401:       default:
                   1402:         return 0;
                   1403:     }
                   1404: LAST:
                   1405:     if ( c > 0 ) return -1;
                   1406:     else return 1;
                   1407:   }
                   1408: }
                   1409:
                   1410: int cmpdl_blex(int n,DL d1,DL d2)
                   1411: {
                   1412:   int c;
                   1413:
                   1414:   if ( (c = cmpdl_lex(n-1,d1,d2)) )
                   1415:     return c;
                   1416:   else {
                   1417:     c = d1->d[n-1] - d2->d[n-1];
                   1418:     return c > 0 ? 1 : c < 0 ? -1 : 0;
                   1419:   }
                   1420: }
                   1421:
                   1422: int cmpdl_bgradlex(int n,DL d1,DL d2)
                   1423: {
                   1424:   int e1,e2,c;
                   1425:
                   1426:   e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
                   1427:   if ( e1 > e2 )
                   1428:     return 1;
                   1429:   else if ( e1 < e2 )
                   1430:     return -1;
                   1431:   else {
                   1432:     c = cmpdl_lex(n-1,d1,d2);
                   1433:     if ( c )
                   1434:       return c;
                   1435:     else
                   1436:       return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
                   1437:   }
                   1438: }
                   1439:
                   1440: int cmpdl_brevgradlex(int n,DL d1,DL d2)
                   1441: {
                   1442:   int e1,e2,c;
                   1443:
                   1444:   e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
                   1445:   if ( e1 > e2 )
                   1446:     return 1;
                   1447:   else if ( e1 < e2 )
                   1448:     return -1;
                   1449:   else {
                   1450:     c = cmpdl_revlex(n-1,d1,d2);
                   1451:     if ( c )
                   1452:       return c;
                   1453:     else
                   1454:       return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
                   1455:   }
                   1456: }
                   1457:
                   1458: int cmpdl_brevrev(int n,DL d1,DL d2)
                   1459: {
                   1460:   int e1,e2,f1,f2,c,i;
                   1461:
                   1462:   for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                   1463:     e1 += d1->d[i]; e2 += d2->d[i];
                   1464:   }
                   1465:   f1 = d1->td - e1; f2 = d2->td - e2;
                   1466:   if ( e1 > e2 )
                   1467:     return 1;
                   1468:   else if ( e1 < e2 )
                   1469:     return -1;
                   1470:   else {
                   1471:     c = cmpdl_revlex(dp_nelim,d1,d2);
                   1472:     if ( c )
                   1473:       return c;
                   1474:     else if ( f1 > f2 )
                   1475:       return 1;
                   1476:     else if ( f1 < f2 )
                   1477:       return -1;
                   1478:     else {
                   1479:       for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                   1480:       return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                   1481:     }
                   1482:   }
                   1483: }
                   1484:
                   1485: int cmpdl_bgradrev(int n,DL d1,DL d2)
                   1486: {
                   1487:   int e1,e2,f1,f2,c,i;
                   1488:
                   1489:   for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                   1490:     e1 += d1->d[i]; e2 += d2->d[i];
                   1491:   }
                   1492:   f1 = d1->td - e1; f2 = d2->td - e2;
                   1493:   if ( e1 > e2 )
                   1494:     return 1;
                   1495:   else if ( e1 < e2 )
                   1496:     return -1;
                   1497:   else {
                   1498:     c = cmpdl_lex(dp_nelim,d1,d2);
                   1499:     if ( c )
                   1500:       return c;
                   1501:     else if ( f1 > f2 )
                   1502:       return 1;
                   1503:     else if ( f1 < f2 )
                   1504:       return -1;
                   1505:     else {
                   1506:       for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                   1507:       return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                   1508:     }
                   1509:   }
                   1510: }
                   1511:
                   1512: int cmpdl_blexrev(int n,DL d1,DL d2)
                   1513: {
                   1514:   int e1,e2,f1,f2,c,i;
                   1515:
                   1516:   for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                   1517:     e1 += d1->d[i]; e2 += d2->d[i];
                   1518:   }
                   1519:   f1 = d1->td - e1; f2 = d2->td - e2;
                   1520:   c = cmpdl_lex(dp_nelim,d1,d2);
                   1521:   if ( c )
                   1522:     return c;
                   1523:   else if ( f1 > f2 )
                   1524:     return 1;
                   1525:   else if ( f1 < f2 )
                   1526:     return -1;
                   1527:   else {
                   1528:     for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                   1529:     return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                   1530:   }
                   1531: }
                   1532:
                   1533: int cmpdl_elim(int n,DL d1,DL d2)
                   1534: {
                   1535:   int e1,e2,i;
                   1536:
                   1537:   for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                   1538:     e1 += d1->d[i]; e2 += d2->d[i];
                   1539:   }
                   1540:   if ( e1 > e2 )
                   1541:     return 1;
                   1542:   else if ( e1 < e2 )
                   1543:     return -1;
                   1544:   else
                   1545:     return cmpdl_revgradlex(n,d1,d2);
                   1546: }
                   1547:
                   1548: int cmpdl_weyl_elim(int n,DL d1,DL d2)
                   1549: {
                   1550:   int e1,e2,i;
                   1551:
                   1552:   for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
                   1553:     e1 += d1->d[n-i]; e2 += d2->d[n-i];
                   1554:   }
                   1555:   if ( e1 > e2 )
                   1556:     return 1;
                   1557:   else if ( e1 < e2 )
                   1558:     return -1;
                   1559:   else if ( d1->td > d2->td )
                   1560:     return 1;
                   1561:   else if ( d1->td < d2->td )
                   1562:     return -1;
                   1563:   else return -cmpdl_revlex(n,d1,d2);
                   1564: }
                   1565:
                   1566: /*
                   1567:   a special ordering
                   1568:   1. total order
                   1569:   2. (-w,w) for the first 2*m variables
                   1570:   3. DRL for the first 2*m variables
                   1571: */
                   1572:
                   1573: extern int *current_weyl_weight_vector;
                   1574:
                   1575: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
                   1576: {
                   1577:   int e1,e2,m,i;
                   1578:   int *p1,*p2;
                   1579:
                   1580:   if ( d1->td > d2->td )
                   1581:     return 1;
                   1582:   else if ( d1->td < d2->td )
                   1583:     return -1;
                   1584:
                   1585:   m = n>>1;
                   1586:   for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
                   1587:     e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
                   1588:     e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
                   1589:   }
                   1590:   if ( e1 > e2 )
                   1591:     return 1;
                   1592:   else if ( e1 < e2 )
                   1593:     return -1;
                   1594:
                   1595:   e1 = d1->td - d1->d[n-1];
                   1596:   e2 = d2->td - d2->d[n-1];
                   1597:   if ( e1 > e2 )
                   1598:     return 1;
                   1599:   else if ( e1 < e2 )
                   1600:     return -1;
                   1601:
                   1602:   for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
                   1603:     i >= 0 && *p1 == *p2; i--, p1--, p2-- );
                   1604:   return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
                   1605: }
                   1606:
                   1607: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
                   1608: {
                   1609:   int i,t,m;
                   1610:   int *p1,*p2;
                   1611:
                   1612:   if ( d1->td > d2->td )
                   1613:     return 1;
                   1614:   else if ( d1->td < d2->td )
                   1615:     return -1;
                   1616:   else {
                   1617:     m = n>>1;
                   1618:     for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
                   1619:       if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1;
                   1620:       if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1;
                   1621:     }
                   1622:     return 0;
                   1623:   }
                   1624: }
                   1625:
                   1626: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
                   1627: {
                   1628:   int e1,e2,m,i,t;
                   1629:   int *p1,*p2;
                   1630:
                   1631:   if ( d1->td > d2->td )
                   1632:     return 1;
                   1633:   else if ( d1->td < d2->td )
                   1634:     return -1;
                   1635:
                   1636:   m = n>>1;
                   1637:   for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
                   1638:     e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
                   1639:     e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
                   1640:   }
                   1641:   if ( e1 > e2 )
                   1642:     return 1;
                   1643:   else if ( e1 < e2 )
                   1644:     return -1;
                   1645:
                   1646:   e1 = d1->td - d1->d[n-1];
                   1647:   e2 = d2->td - d2->d[n-1];
                   1648:   if ( e1 > e2 )
                   1649:     return 1;
                   1650:   else if ( e1 < e2 )
                   1651:     return -1;
                   1652:
                   1653:   for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
                   1654:     if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1;
                   1655:     if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1;
                   1656:   }
                   1657:   return 0;
                   1658: }
                   1659:
                   1660: int cmpdl_order_pair(int n,DL d1,DL d2)
                   1661: {
                   1662:   int e1,e2,i,j,l;
                   1663:   int *t1,*t2;
                   1664:   int len,head;
                   1665:   struct order_pair *pair;
                   1666:
                   1667:   len = dp_current_spec->ord.block.length;
                   1668:   if ( n != dp_current_spec->nv )
                   1669:     error("cmpdl_order_pair : incompatible order specification");
                   1670:   pair = dp_current_spec->ord.block.order_pair;
                   1671:
                   1672:   head = 0;
                   1673:   for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
                   1674:     l = pair[i].length;
                   1675:     switch ( pair[i].order ) {
                   1676:       case 0:
                   1677:         for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
                   1678:           e1 += MUL_WEIGHT(t1[j],head+j);
                   1679:           e2 += MUL_WEIGHT(t2[j],head+j);
                   1680:         }
                   1681:         if ( e1 > e2 )
                   1682:           return 1;
                   1683:         else if ( e1 < e2 )
                   1684:           return -1;
                   1685:         else {
                   1686:           for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
                   1687:           if ( j >= 0 )
                   1688:             return t1[j] < t2[j] ? 1 : -1;
                   1689:         }
                   1690:         break;
                   1691:       case 1:
                   1692:         for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
                   1693:           e1 += MUL_WEIGHT(t1[j],head+j);
                   1694:           e2 += MUL_WEIGHT(t2[j],head+j);
                   1695:         }
                   1696:         if ( e1 > e2 )
                   1697:           return 1;
                   1698:         else if ( e1 < e2 )
                   1699:           return -1;
                   1700:         else {
                   1701:           for ( j = 0; j < l && t1[j] == t2[j]; j++ );
                   1702:           if ( j < l )
                   1703:             return t1[j] > t2[j] ? 1 : -1;
                   1704:         }
                   1705:         break;
                   1706:       case 2:
                   1707:         for ( j = 0; j < l && t1[j] == t2[j]; j++ );
                   1708:         if ( j < l )
                   1709:           return t1[j] > t2[j] ? 1 : -1;
                   1710:         break;
                   1711:       default:
                   1712:         error("cmpdl_order_pair : invalid order"); break;
                   1713:     }
                   1714:     t1 += l; t2 += l; head += l;
                   1715:   }
                   1716:   return 0;
                   1717: }
                   1718:
                   1719: int cmpdl_composite(int nv,DL d1,DL d2)
                   1720: {
                   1721:   int n,i,j,k,start,s,len;
                   1722:   int *dw;
                   1723:   struct sparse_weight *sw;
                   1724:   struct weight_or_block *worb;
                   1725:   int *w,*t1,*t2;
                   1726:
                   1727:   n = dp_current_spec->ord.composite.length;
                   1728:   worb = dp_current_spec->ord.composite.w_or_b;
                   1729:   w = dp_dl_work;
                   1730:   for ( i = 0, t1 = d1->d, t2 = d2->d; i < nv; i++ )
                   1731:     w[i] = t1[i]-t2[i];
                   1732:   for ( i = 0; i < n; i++, worb++ ) {
                   1733:     len = worb->length;
                   1734:     switch ( worb->type ) {
                   1735:       case IS_DENSE_WEIGHT:
                   1736:         dw = worb->body.dense_weight;
                   1737:         for ( j = 0, s = 0; j < len; j++ )
                   1738:           s += dw[j]*w[j];
                   1739:         if ( s > 0 ) return 1;
                   1740:         else if ( s < 0 ) return -1;
                   1741:         break;
                   1742:       case IS_SPARSE_WEIGHT:
                   1743:         sw = worb->body.sparse_weight;
                   1744:         for ( j = 0, s = 0; j < len; j++ )
                   1745:           s += sw[j].value*w[sw[j].pos];
                   1746:         if ( s > 0 ) return 1;
                   1747:         else if ( s < 0 ) return -1;
                   1748:         break;
                   1749:       case IS_BLOCK:
                   1750:         start = worb->body.block.start;
                   1751:         switch ( worb->body.block.order ) {
                   1752:           case 0:
                   1753:             for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
                   1754:               s += MUL_WEIGHT(w[k],k);
                   1755:             }
                   1756:             if ( s > 0 ) return 1;
                   1757:             else if ( s < 0 ) return -1;
                   1758:             else {
                   1759:               for ( j = k-1; j >= start && w[j] == 0; j-- );
                   1760:               if ( j >= start )
                   1761:                 return w[j] < 0 ? 1 : -1;
                   1762:             }
                   1763:             break;
                   1764:           case 1:
                   1765:             for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
                   1766:               s += MUL_WEIGHT(w[k],k);
                   1767:             }
                   1768:             if ( s > 0 ) return 1;
                   1769:             else if ( s < 0 ) return -1;
                   1770:             else {
                   1771:               for ( j = 0, k = start;  j < len && w[j] == 0; j++, k++ );
                   1772:               if ( j < len )
                   1773:                 return w[j] > 0 ? 1 : -1;
                   1774:             }
                   1775:             break;
                   1776:           case 2:
                   1777:             for ( j = 0, k = start;  j < len && w[j] == 0; j++, k++ );
                   1778:             if ( j < len )
                   1779:               return w[j] > 0 ? 1 : -1;
                   1780:             break;
                   1781:         }
                   1782:         break;
                   1783:     }
                   1784:   }
                   1785:   return 0;
                   1786: }
                   1787:
                   1788: int cmpdl_matrix(int n,DL d1,DL d2)
                   1789: {
                   1790:   int *v,*w,*t1,*t2;
                   1791:   int s,i,j,len;
                   1792:   int **matrix;
                   1793:
                   1794:   for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
                   1795:     w[i] = t1[i]-t2[i];
                   1796:   len = dp_current_spec->ord.matrix.row;
                   1797:   matrix = dp_current_spec->ord.matrix.matrix;
                   1798:   for ( j = 0; j < len; j++ ) {
                   1799:     v = matrix[j];
                   1800:     for ( i = 0, s = 0; i < n; i++ )
                   1801:       s += v[i]*w[i];
                   1802:     if ( s > 0 )
                   1803:       return 1;
                   1804:     else if ( s < 0 )
                   1805:       return -1;
                   1806:   }
                   1807:   return 0;
                   1808: }
                   1809:
                   1810: int cmpdl_top_weight(int n,DL d1,DL d2)
                   1811: {
                   1812:   int *w;
                   1813:   Z **mat;
                   1814:   Z *a;
                   1815:   mpz_t sum;
                   1816:   int len,i,sgn,tsgn,row,k;
                   1817:   int *t1,*t2;
                   1818:
                   1819:   w = (int *)ALLOCA(n*sizeof(int));
                   1820:   len = current_top_weight_len+3;
                   1821:   t1 = d1->d; t2 = d2->d;
                   1822:   for ( i = 0; i < n; i++ ) w[i] = t1[i]-t2[i];
                   1823:   mpz_init_set_ui(sum,0);
                   1824:   if ( OID(current_top_weight) == O_VECT ) {
                   1825:       mat = (Z **)&BDY((VECT)current_top_weight);
                   1826:     row = 1;
                   1827:   } else {
                   1828:       mat = (Z **)BDY((MAT)current_top_weight);
                   1829:     row = ((MAT)current_top_weight)->row;
                   1830:   }
                   1831:   for ( k = 0; k < row; k++ ) {
                   1832:     a = mat[k];
                   1833:     for ( i = 0; i < n; i++ ) {
                   1834:       if ( !a[i] || !w[i] ) continue;
                   1835:       if ( w[i] > 0 )
                   1836:         mpz_addmul_ui(sum,BDY(a[i]),(unsigned int)w[i]);
                   1837:       else
                   1838:         mpz_submul_ui(sum,BDY(a[i]),(unsigned int)(-w[i]));
                   1839:     }
                   1840:     sgn = mpz_sgn(sum);
                   1841:     if ( sgn > 0 ) return 1;
                   1842:     else if ( sgn < 0 ) return -1;
                   1843:   }
                   1844:   return (*cmpdl_tie_breaker)(n,d1,d2);
                   1845: }
                   1846:
                   1847: GeoBucket create_bucket()
                   1848: {
                   1849:   GeoBucket g;
                   1850:
                   1851:   g = CALLOC(1,sizeof(struct oGeoBucket));
                   1852:   g->m = 32;
                   1853:   return g;
                   1854: }
                   1855:
                   1856: int length(NODE d);
                   1857:
                   1858: void add_bucket(GeoBucket g,NODE d,int nv)
                   1859: {
                   1860:   int l,k,m;
                   1861:
                   1862:   l = length(d);
                   1863:   for ( k = 0, m = 1; l > m; k++, m <<= 1 );
                   1864:   /* 2^(k-1) < l <= 2^k */
                   1865:   d = symb_merge(g->body[k],d,nv);
                   1866:   for ( ; length(d) > (1<<(k)); k++ ) {
                   1867:     g->body[k] = 0;
                   1868:     d = symb_merge(g->body[k+1],d,nv);
                   1869:   }
                   1870:   g->body[k] = d;
                   1871:   g->m = MAX(g->m,k);
                   1872: }
                   1873:
                   1874: DL remove_head_bucket(GeoBucket g,int nv)
                   1875: {
                   1876:   int j,i,c,m;
                   1877:   DL d;
                   1878:
                   1879:   j = -1;
                   1880:   m = g->m;
                   1881:   for ( i = 0; i <= m; i++ ) {
                   1882:     if ( !g->body[i] )
                   1883:       continue;
                   1884:     if ( j < 0 ) j = i;
                   1885:     else {
                   1886:       c = (*cmpdl)(nv,g->body[i]->body,g->body[j]->body);
                   1887:       if ( c > 0 )
                   1888:         j = i;
                   1889:       else if ( c == 0 )
                   1890:         g->body[i] = NEXT(g->body[i]);
                   1891:     }
                   1892:   }
                   1893:   if ( j < 0 )
                   1894:     return 0;
                   1895:   else {
                   1896:     d = g->body[j]->body;
                   1897:     g->body[j] = NEXT(g->body[j]);
                   1898:     return d;
                   1899:   }
                   1900: }
                   1901:
                   1902: /*  DPV functions */
                   1903:
                   1904: void adddv(VL vl,DPV p1,DPV p2,DPV *pr)
                   1905: {
                   1906:   int i,len;
                   1907:   DP *e;
                   1908:
                   1909:   if ( !p1 || !p2 )
                   1910:     error("adddv : invalid argument");
                   1911:   else if ( p1->len != p2->len )
                   1912:     error("adddv : size mismatch");
                   1913:   else {
                   1914:     len = p1->len;
                   1915:     e = (DP *)MALLOC(p1->len*sizeof(DP));
                   1916:     for ( i = 0; i < len; i++ )
                   1917:       addd(vl,p1->body[i],p2->body[i],&e[i]);
                   1918:     MKDPV(len,e,*pr);
                   1919:     (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                   1920:   }
                   1921: }
                   1922:
                   1923: void subdv(VL vl,DPV p1,DPV p2,DPV *pr)
                   1924: {
                   1925:   int i,len;
                   1926:   DP *e;
                   1927:
                   1928:   if ( !p1 || !p2 )
                   1929:     error("subdv : invalid argument");
                   1930:   else if ( p1->len != p2->len )
                   1931:     error("subdv : size mismatch");
                   1932:   else {
                   1933:     len = p1->len;
                   1934:     e = (DP *)MALLOC(p1->len*sizeof(DP));
                   1935:     for ( i = 0; i < len; i++ )
                   1936:       subd(vl,p1->body[i],p2->body[i],&e[i]);
                   1937:     MKDPV(len,e,*pr);
                   1938:     (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                   1939:   }
                   1940: }
                   1941:
                   1942: void chsgndv(DPV p1,DPV *pr)
                   1943: {
                   1944:   int i,len;
                   1945:   DP *e;
                   1946:
                   1947:   if ( !p1 )
                   1948:     error("subdv : invalid argument");
                   1949:   else {
                   1950:     len = p1->len;
                   1951:     e = (DP *)MALLOC(p1->len*sizeof(DP));
                   1952:     for ( i = 0; i < len; i++ )
                   1953:       chsgnd(p1->body[i],&e[i]);
                   1954:     MKDPV(len,e,*pr);
                   1955:     (*pr)->sugar = p1->sugar;
                   1956:   }
                   1957: }
                   1958:
                   1959: void muldv(VL vl,DP p1,DPV p2,DPV *pr)
                   1960: {
                   1961:   int i,len;
                   1962:   DP *e;
                   1963:
                   1964:   len = p2->len;
                   1965:   e = (DP *)MALLOC(p2->len*sizeof(DP));
                   1966:   if ( !p1 ) {
                   1967:     MKDPV(len,e,*pr);
                   1968:     (*pr)->sugar = 0;
                   1969:   } else {
                   1970:     for ( i = 0; i < len; i++ )
                   1971:       muld(vl,p1,p2->body[i],&e[i]);
                   1972:     MKDPV(len,e,*pr);
                   1973:     (*pr)->sugar = p1->sugar + p2->sugar;
                   1974:   }
                   1975: }
                   1976:
                   1977: int compdv(VL vl,DPV p1,DPV p2)
                   1978: {
                   1979:   int i,t,len;
                   1980:
                   1981:   if ( p1->len != p2->len ) {
                   1982:     error("compdv : size mismatch");
                   1983:     return 0; /* XXX */
                   1984:   } else {
                   1985:     len = p1->len;
                   1986:     for ( i = 0; i < len; i++ )
                   1987:       if ( (t = compd(vl,p1->body[i],p2->body[i])) )
                   1988:         return t;
                   1989:     return 0;
                   1990:   }
                   1991: }
                   1992:
                   1993: int ni_next(int *a,int n)
                   1994: {
                   1995:   int i,j,k,kj;
                   1996:
                   1997:   /* find the first nonzero a[j] */
                   1998:   for ( j = 0; j < n && a[j] == 0; j++ );
                   1999:   /* find the first zero a[k] after a[j] */
                   2000:   for ( k = j; k < n && a[k] == 1; k++ );
                   2001:   if ( k == n ) return 0;
                   2002:   /* a[0] = 0, ... , a[j-1] = 0, a[j] = 1, ..., a[k-1] = 1, a[k] = 0 */
                   2003:   /* a[0] = 1,..., a[k-j-2] = 1, a[k-j-1] = 0, ..., a[k-1] = 0, a[k] = 1 */
                   2004:   kj = k-j-1;
                   2005:   for ( i = 0; i < kj; i++ ) a[i] = 1;
                   2006:   for ( ; i < k; i++ ) a[i] = 0;
                   2007:   a[k] = 1;
                   2008:   return 1;
                   2009: }
                   2010:
                   2011: int comp_nbm(NBM a,NBM b)
                   2012: {
                   2013:   int d,i,ai,bi;
                   2014:   unsigned int *ab,*bb;
                   2015:
                   2016:   if ( a->d > b->d ) return 1;
                   2017:   else if ( a->d < b->d ) return -1;
                   2018:   else {
                   2019:     d = a->d; ab = a->b; bb = b->b;
                   2020: #if 0
                   2021:     w = (d+31)/32;
                   2022:     for ( i = 0; i < w; i++ )
                   2023:       if ( ab[i] > bb[i] ) return 1;
                   2024:       else if ( ab[i] < bb[i] ) return -1;
                   2025: #else
                   2026:     for ( i = 0; i < d; i++ ) {
                   2027:       ai = NBM_GET(ab,i);
                   2028:       bi = NBM_GET(bb,i);
                   2029:       if ( ai > bi ) return 1;
                   2030:       else if ( ai < bi ) return -1;
                   2031:     }
                   2032: #endif
                   2033:     return 0;
                   2034:   }
                   2035: }
                   2036:
                   2037: NBM mul_nbm(NBM a,NBM b)
                   2038: {
                   2039:   int ad,bd,d,i,j;
                   2040:   unsigned int *ab,*bb,*mb;
                   2041:   NBM m;
                   2042:
                   2043:   ad = a->d; bd = b->d; ab = a->b; bb = b->b;
                   2044:   d = ad + bd;
                   2045:   NEWNBM(m); NEWNBMBDY(m,d);
                   2046:   m->d = d; mulp(CO,a->c,b->c,&m->c); mb = m->b;
                   2047:   j = 0;
                   2048:   for ( i = 0; i < ad; i++, j++ )
                   2049:     if ( NBM_GET(ab,i) ) NBM_SET(mb,j);
                   2050:     else NBM_CLR(mb,j);
                   2051:   for ( i = 0; i < bd; i++, j++ )
                   2052:     if ( NBM_GET(bb,i) ) NBM_SET(mb,j);
                   2053:     else NBM_CLR(mb,j);
                   2054:   return m;
                   2055: }
                   2056:
                   2057: NBP nbmtonbp(NBM m)
                   2058: {
                   2059:   NODE n;
                   2060:   NBP u;
                   2061:
                   2062:   MKNODE(n,m,0);
                   2063:   MKNBP(u,n);
                   2064:   return u;
                   2065: }
                   2066:
                   2067: /* a=c*x*rest -> a0= x*rest, ah=x, ar=rest */
                   2068:
                   2069: P separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
                   2070: {
                   2071:   int i,d1;
                   2072:   NBM t;
                   2073:
                   2074:   if ( !a->d ) error("separate_nbm : invalid argument");
                   2075:
                   2076:   if ( a0 ) {
                   2077:     NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
                   2078:     *a0 = nbmtonbp(t);
                   2079:   }
                   2080:
                   2081:   if ( ah ) {
                   2082:     NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
                   2083:     if ( NBM_GET(a->b,0) ) NBM_SET(t->b,0);
                   2084:     else NBM_CLR(t->b,0);
                   2085:     *ah = nbmtonbp(t);
                   2086:   }
                   2087:
                   2088:   if ( ar ) {
                   2089:     d1 = a->d-1;
                   2090:     NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
                   2091:     for ( i = 0; i < d1; i++ ) {
                   2092:       if ( NBM_GET(a->b,i+1) ) NBM_SET(t->b,i);
                   2093:       else NBM_CLR(t->b,i);
                   2094:     }
                   2095:     *ar = nbmtonbp(t);
                   2096:   }
                   2097:
                   2098:   return a->c;
                   2099: }
                   2100:
                   2101: /* a=c*rest*x -> a0= rest*x, ar=rest, at=x */
                   2102:
                   2103: P separate_tail_nbm(NBM a,NBP *a0,NBP *ar,NBP *at)
                   2104: {
                   2105:   int i,d,d1;
                   2106:   NBM t;
                   2107:
                   2108:   if ( !(d=a->d) ) error("separate_tail_nbm : invalid argument");
                   2109:
                   2110:   if ( a0 ) {
                   2111:     NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
                   2112:     *a0 = nbmtonbp(t);
                   2113:   }
                   2114:
                   2115:   d1 = a->d-1;
                   2116:   if ( at ) {
                   2117:     NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
                   2118:     if ( NBM_GET(a->b,d1) ) NBM_SET(t->b,0);
                   2119:     else NBM_CLR(t->b,0);
                   2120:     *at = nbmtonbp(t);
                   2121:   }
                   2122:
                   2123:   if ( ar ) {
                   2124:     NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
                   2125:     for ( i = 0; i < d1; i++ ) {
                   2126:       if ( NBM_GET(a->b,i) ) NBM_SET(t->b,i);
                   2127:       else NBM_CLR(t->b,i);
                   2128:     }
                   2129:     *ar = nbmtonbp(t);
                   2130:   }
                   2131:
                   2132:   return a->c;
                   2133: }
                   2134:
                   2135: NBP make_xky(int k)
                   2136: {
                   2137:   int k1,i;
                   2138:   NBM t;
                   2139:
                   2140:   NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
                   2141:   k1 = k-1;
                   2142:   for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
                   2143:   NBM_CLR(t->b,i);
                   2144:   return nbmtonbp(t);
                   2145: }
                   2146:
                   2147: /* a=c*x^(k-1)*y*rest -> a0= x^(k-1)*y*rest, ah=x^(k-1)*y, ar=rest */
                   2148:
                   2149: P separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
                   2150: {
                   2151:   int i,d1,k,k1;
                   2152:   NBM t;
                   2153:
                   2154:   if ( !a->d )
                   2155:     error("separate_nbm : invalid argument");
                   2156:   for ( i = 0; i < a->d && NBM_GET(a->b,i); i++ );
                   2157:   if ( i == a->d )
                   2158:     error("separate_nbm : invalid argument");
                   2159:   k1 = i;
                   2160:   k = i+1;
                   2161:
                   2162:   if ( a0 ) {
                   2163:     NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
                   2164:     *a0 = nbmtonbp(t);
                   2165:   }
                   2166:
                   2167:   if ( ah ) {
                   2168:     NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
                   2169:     for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
                   2170:     NBM_CLR(t->b,i);
                   2171:     *ah = nbmtonbp(t);
                   2172:   }
                   2173:
                   2174:   if ( ar ) {
                   2175:     d1 = a->d-k;
                   2176:     NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
                   2177:     for ( i = 0; i < d1; i++ ) {
                   2178:       if ( NBM_GET(a->b,i+k) ) NBM_SET(t->b,i);
                   2179:       else NBM_CLR(t->b,i);
                   2180:     }
                   2181:     *ar = nbmtonbp(t);
                   2182:   }
                   2183:
                   2184:   return a->c;
                   2185: }
                   2186:
                   2187: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
                   2188: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
                   2189: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp);
                   2190: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
                   2191:
                   2192: NBP shuffle_mul_nbm(NBM a,NBM b)
                   2193: {
                   2194:   NBP u,a0,ah,ar,b0,bh,br,a1,b1,t;
                   2195:   P ac,bc,c;
                   2196:
                   2197:   if ( !a->d || !b->d )
                   2198:     u = nbmtonbp(mul_nbm(a,b));
                   2199:   else {
                   2200:     ac = separate_nbm(a,&a0,&ah,&ar);
                   2201:     bc = separate_nbm(b,&b0,&bh,&br);
                   2202:     mulp(CO,ac,bc,&c);
                   2203:     shuffle_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
                   2204:     shuffle_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
                   2205:     addnbp(CO,a1,b1,&t); mulnbp(CO,(NBP)c,t,&u);
                   2206:   }
                   2207:   return u;
                   2208: }
                   2209:
                   2210: NBP harmonic_mul_nbm(NBM a,NBM b)
                   2211: {
                   2212:   NBP u,a0,ah,ar,b0,bh,br,a1,b1,t,s,abk,ab1;
                   2213:   P ac,bc,c;
                   2214:
                   2215:   if ( !a->d || !b->d )
                   2216:     u = nbmtonbp(mul_nbm(a,b));
                   2217:   else {
                   2218:     mulp(CO,a->c,b->c,&c);
                   2219:     ac = separate_xky_nbm(a,&a0,&ah,&ar);
                   2220:     bc = separate_xky_nbm(b,&b0,&bh,&br);
                   2221:     mulp(CO,ac,bc,&c);
                   2222:     harmonic_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
                   2223:     harmonic_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
                   2224:     abk = make_xky(((NBM)BDY(BDY(ah)))->d+((NBM)BDY(BDY(bh)))->d);
                   2225:     harmonic_mulnbp(CO,ar,br,&t); mulnbp(CO,abk,t,&ab1);
                   2226:     addnbp(CO,a1,b1,&t); addnbp(CO,t,ab1,&s); mulnbp(CO,(NBP)c,s,&u);
                   2227:   }
                   2228:   return u;
                   2229:
                   2230: }
                   2231:
                   2232: void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2233: {
                   2234:   NODE b1,b2,br=0,br0;
                   2235:   NBM m1,m2,m;
                   2236:   P c;
                   2237:
                   2238:   if ( !p1 )
                   2239:     *rp = p2;
                   2240:   else if ( !p2 )
                   2241:     *rp = p1;
                   2242:   else {
                   2243:     for ( b1 = BDY(p1), b2 = BDY(p2), br0 = 0; b1 && b2; ) {
                   2244:       m1 = (NBM)BDY(b1); m2 = (NBM)BDY(b2);
                   2245:       switch ( comp_nbm(m1,m2) ) {
                   2246:         case 0:
                   2247:           addp(CO,m1->c,m2->c,&c);
                   2248:           if ( c ) {
                   2249:             NEXTNODE(br0,br);
                   2250:             NEWNBM(m); m->d = m1->d; m->c = c; m->b = m1->b;
                   2251:             BDY(br) = (pointer)m;
                   2252:           }
                   2253:           b1 = NEXT(b1); b2 = NEXT(b2); break;
                   2254:         case 1:
                   2255:           NEXTNODE(br0,br); BDY(br) = BDY(b1);
                   2256:           b1 = NEXT(b1); break;
                   2257:         case -1:
                   2258:           NEXTNODE(br0,br); BDY(br) = BDY(b2);
                   2259:           b2 = NEXT(b2); break;
                   2260:       }
                   2261:     }
                   2262:     if ( !br0 )
                   2263:       if ( b1 )
                   2264:         br0 = b1;
                   2265:       else if ( b2 )
                   2266:         br0 = b2;
                   2267:       else {
                   2268:         *rp = 0;
                   2269:         return;
                   2270:       }
                   2271:     else if ( b1 )
                   2272:       NEXT(br) = b1;
                   2273:     else if ( b2 )
                   2274:         NEXT(br) = b2;
                   2275:     else
                   2276:       NEXT(br) = 0;
                   2277:     MKNBP(*rp,br0);
                   2278:   }
                   2279: }
                   2280:
                   2281: void subnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2282: {
                   2283:   NBP t;
                   2284:
                   2285:   chsgnnbp(p2,&t);
                   2286:   addnbp(vl,p1,t,rp);
                   2287: }
                   2288:
                   2289: void chsgnnbp(NBP p,NBP *rp)
                   2290: {
                   2291:   NODE r0,r=0,b;
                   2292:   NBM m,m1;
                   2293:
                   2294:   for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
                   2295:     NEXTNODE(r0,r);
                   2296:     m = (NBM)BDY(b);
                   2297:     NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnp(m->c,&m1->c);
                   2298:     BDY(r) = m1;
                   2299:   }
                   2300:   if ( r0 ) NEXT(r) = 0;
                   2301:   MKNBP(*rp,r0);
                   2302: }
                   2303:
                   2304: void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2305: {
                   2306:   NODE b,n;
                   2307:   NBP r,t,s;
                   2308:   NBM m;
                   2309:
                   2310:   if ( !p1 || !p2 ) {
                   2311:     *rp = 0; return;
                   2312:   }
                   2313:   if ( OID(p1) != O_NBP ) {
                   2314:     if ( !POLY(p1) )
                   2315:       error("mulnbp : invalid argument");
                   2316:     NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
                   2317:     MKNODE(n,m,0); MKNBP(p1,n);
                   2318:   }
                   2319:   if ( OID(p2) != O_NBP ) {
                   2320:     if ( !POLY(p2) )
                   2321:       error("mulnbp : invalid argument");
                   2322:     NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
                   2323:     MKNODE(n,m,0); MKNBP(p2,n);
                   2324:   }
                   2325:   if ( length(BDY(p1)) < length(BDY(p2)) ) {
                   2326:     for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {
                   2327:       mulnbmnbp(vl,(NBM)BDY(b),p2,&t);
                   2328:       addnbp(vl,r,t,&s); r = s;
                   2329:     }
                   2330:     *rp = r;
                   2331:   } else {
                   2332:     for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) {
                   2333:       mulnbpnbm(vl,p1,(NBM)BDY(b),&t);
                   2334:       addnbp(vl,r,t,&s); r = s;
                   2335:     }
                   2336:     *rp = r;
                   2337:   }
                   2338: }
                   2339:
                   2340: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp)
                   2341: {
                   2342:   NODE b,r0,r=0;
                   2343:
                   2344:   if ( !p ) *rp = 0;
                   2345:   else {
                   2346:     for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
                   2347:       NEXTNODE(r0,r);
                   2348:       BDY(r) = mul_nbm(m,(NBM)BDY(b));
                   2349:     }
                   2350:     if ( r0 ) NEXT(r) = 0;
                   2351:     MKNBP(*rp,r0);
                   2352:   }
                   2353: }
                   2354:
                   2355: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp)
                   2356: {
                   2357:   NODE b,r0,r=0;
                   2358:
                   2359:   if ( !p ) *rp = 0;
                   2360:   else {
                   2361:     for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
                   2362:       NEXTNODE(r0,r);
                   2363:       BDY(r) = mul_nbm((NBM)BDY(b),m);
                   2364:     }
                   2365:     if ( r0 ) NEXT(r) = 0;
                   2366:     MKNBP(*rp,r0);
                   2367:   }
                   2368: }
                   2369:
                   2370: void pwrnbp(VL vl,NBP a,Z q,NBP *c)
                   2371: {
                   2372:   NBP a1,a2;
                   2373:   Z q1,r1,two;
                   2374:   NBM m;
                   2375:   NODE r;
                   2376:
                   2377:   if ( !q ) {
                   2378:      NEWNBM(m); m->d = 0; m->c = (P)ONE; m->b = 0;
                   2379:      MKNODE(r,m,0); MKNBP(*c,r);
                   2380:   } else if ( !a )
                   2381:     *c = 0;
                   2382:   else if ( UNIQ(q) )
                   2383:     *c = a;
                   2384:   else {
1.2       noro     2385:     STOZ(2,two);
1.1       noro     2386:     divqrz(q,two,&q1,&r1);
                   2387:     pwrnbp(vl,a,q1,&a1);
                   2388:     mulnbp(vl,a1,a1,&a2);
                   2389:     if ( r1 )
                   2390:       mulnbp(vl,a2,a,c);
                   2391:     else
                   2392:       *c = a2;
                   2393:   }
                   2394: }
                   2395:
                   2396: int compnbp(VL vl,NBP p1,NBP p2)
                   2397: {
                   2398:   NODE n1,n2;
                   2399:   NBM m1,m2;
                   2400:   int t;
                   2401:
                   2402:   if ( !p1 )
                   2403:     return p2 ? -1 : 0;
                   2404:   else if ( !p2 )
                   2405:     return 1;
                   2406:   else {
                   2407:     for ( n1 = BDY(p1), n2 = BDY(p2);
                   2408:       n1 && n2; n1 = NEXT(n1), n2 = NEXT(n2) ) {
                   2409:       m1 = (NBM)BDY(n1); m2 = (NBM)BDY(n2);
                   2410:       if ( (t = comp_nbm(m1,m2)) || (t = compp(CO,m1->c,m2->c) ) )
                   2411:         return t;
                   2412:     }
                   2413:     if ( n1 )
                   2414:       return 1;
                   2415:     else if ( n2 )
                   2416:       return -1;
                   2417:     else
                   2418:       return 0;
                   2419:   }
                   2420: }
                   2421:
                   2422: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2423: {
                   2424:   NODE b1,b2,n;
                   2425:   NBP r,t,s;
                   2426:   NBM m;
                   2427:
                   2428:   if ( !p1 || !p2 ) {
                   2429:     *rp = 0; return;
                   2430:   }
                   2431:   if ( OID(p1) != O_NBP ) {
                   2432:     if ( !POLY(p1) )
                   2433:       error("shuffle_mulnbp : invalid argument");
                   2434:     NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
                   2435:     MKNODE(n,m,0); MKNBP(p1,n);
                   2436:   }
                   2437:   if ( OID(p2) != O_NBP ) {
                   2438:     if ( !POLY(p2) )
                   2439:       error("shuffle_mulnbp : invalid argument");
                   2440:     NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
                   2441:     MKNODE(n,m,0); MKNBP(p2,n);
                   2442:   }
                   2443:   for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
                   2444:     for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
                   2445:       t = shuffle_mul_nbm(m,(NBM)BDY(b2));
                   2446:       addnbp(vl,r,t,&s); r = s;
                   2447:     }
                   2448:   *rp = r;
                   2449: }
                   2450:
                   2451: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2452: {
                   2453:   NODE b1,b2,n;
                   2454:   NBP r,t,s;
                   2455:   NBM m;
                   2456:
                   2457:   if ( !p1 || !p2 ) {
                   2458:     *rp = 0; return;
                   2459:   }
                   2460:   if ( OID(p1) != O_NBP ) {
                   2461:     if ( !POLY(p1) )
                   2462:       error("harmonic_mulnbp : invalid argument");
                   2463:     NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
                   2464:     MKNODE(n,m,0); MKNBP(p1,n);
                   2465:   }
                   2466:   if ( OID(p2) != O_NBP ) {
                   2467:     if ( !POLY(p2) )
                   2468:       error("harmonic_mulnbp : invalid argument");
                   2469:     NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
                   2470:     MKNODE(n,m,0); MKNBP(p2,n);
                   2471:   }
                   2472:   for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
                   2473:     for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
                   2474:       t = harmonic_mul_nbm(m,(NBM)BDY(b2));
                   2475:       addnbp(vl,r,t,&s); r = s;
                   2476:     }
                   2477:   *rp = r;
                   2478: }
                   2479:
                   2480: #if 0
                   2481: NBP shuffle_mul_nbm(NBM a,NBM b)
                   2482: {
                   2483:   int ad,bd,d,i,ai,bi,bit,s;
                   2484:   int *ab,*bb,*wmb,*w;
                   2485:   NBM wm,tm;
                   2486:   P c,c1;
                   2487:   NODE r,t,t1,p;
                   2488:   NBP u;
                   2489:
                   2490:   ad = a->d; bd = b->d; ab = a->b; bb = b->b;
                   2491:   d = ad + bd;
                   2492:   w = (int *)ALLOCA(d*sizeof(int));
                   2493:   NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   2494:   for ( i = 0; i < ad; i++ ) w[i] = 1;
                   2495:   for ( ; i < d; i++ ) w[i] = 0;
                   2496:   mulp(CO,a->c,b->c,&c);
                   2497:   r = 0;
                   2498:   do {
                   2499:     wm->d = d; wm->c = c;
                   2500:     ai = 0; bi = 0;
                   2501:     for ( i = 0; i < d; i++ ) {
                   2502:       if ( w[i] ) { bit = NBM_GET(ab,ai); ai++; }
                   2503:       else { bit = NBM_GET(bb,bi); bi++; }
                   2504:       if ( bit ) NBM_SET(wmb,i);
                   2505:       else NBM_CLR(wmb,i);
                   2506:     }
                   2507:     for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
                   2508:       tm = (NBM)BDY(t);
                   2509:       s = comp_nbm(tm,wm);
                   2510:       if ( s < 0 ) {
                   2511:         /* insert */
                   2512:         MKNODE(t1,wm,t);
                   2513:         if ( !p ) r = t1;
                   2514:         else NEXT(p) = t1;
                   2515:         NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   2516:         break;
                   2517:       } else if ( s == 0 ) {
                   2518:         /* add coefs */
                   2519:         addp(CO,tm->c,c,&c1);
                   2520:         if ( c1 ) tm->c = c1;
                   2521:         else NEXT(p) = NEXT(t);
                   2522:         break;
                   2523:       }
                   2524:     }
                   2525:     if ( !t ) {
                   2526:       /* append */
                   2527:       MKNODE(t1,wm,t);
                   2528:       if ( !p ) r = t1;
                   2529:       else NEXT(p) = t1;
                   2530:       NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   2531:     }
                   2532:   } while ( ni_next(w,d) );
                   2533:   MKNBP(u,r);
                   2534:   return u;
                   2535: }
                   2536:
                   2537: int nbmtoxky(NBM a,int *b)
                   2538: {
                   2539:   int d,i,j,k;
                   2540:   int *p;
                   2541:
                   2542:   d = a->d; p = a->b;
                   2543:   for ( i = j = 0, k = 1; i < d; i++ ) {
                   2544:     if ( !NBM_GET(p,i) ) {
                   2545:       b[j++] = k;
                   2546:       k = 1;
                   2547:     } else k++;
                   2548:   }
                   2549:   return j;
                   2550: }
                   2551:
                   2552: NBP harmonic_mul_nbm(NBM a,NBM b)
                   2553: {
                   2554:   int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1;
                   2555:   int i,j,k,ia,ib,s;
                   2556:   int *wa,*wb,*w,*wab,*wa1,*wmb;
                   2557:   P c,c1;
                   2558:   NBM wm,tm;
                   2559:   NODE r,t1,t,p;
                   2560:   NBP u;
                   2561:
                   2562:   da = a->d; db = b->d; d = da+db;
                   2563:   wa = (int *)ALLOCA(da*sizeof(int));
                   2564:   wb = (int *)ALLOCA(db*sizeof(int));
                   2565:   la = nbmtoxky(a,wa);
                   2566:   lb = nbmtoxky(b,wb);
                   2567:   mulp(CO,a->c,b->c,&c);
                   2568:   /* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */
                   2569:   /* lmax : total length */
                   2570:   lmax = la+lb;
                   2571:   lmin = la>lb?la:lb;
                   2572:   w = (int *)ALLOCA(lmax*sizeof(int));
                   2573:   /* position of a+b */
                   2574:   wab = (int *)ALLOCA(lmax*sizeof(int));
                   2575:   /* position of a */
                   2576:   wa1 = (int *)ALLOCA(lmax*sizeof(int));
                   2577:   NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   2578:   for ( l = lmin, r = 0; l <= lmax; l++ ) {
                   2579:     lab = lmax - l;
                   2580:     la1 = la - lab;
                   2581:     lb1 = lb - lab;
                   2582:     lab1 = l-lab;
                   2583:     /* partion l into three parts: a, b, a+b */
                   2584:     /* initialize wab */
                   2585:     for ( i = 0; i < lab; i++ ) wab[i] = 1;
                   2586:     for ( ; i < l; i++ ) wab[i] = 0;
                   2587:     do {
                   2588:       /* initialize wa1 */
                   2589:       for ( i = 0; i < la1; i++ ) wa1[i] = 1;
                   2590:       for ( ; i < lab1; i++ ) wa1[i] = 0;
                   2591:       do {
                   2592:         ia = 0; ib = 0;
                   2593:         for ( i = j = 0; i < l; i++ )
                   2594:           if ( wab[i] ) w[i] = wa[ia++]+wb[ib++];
                   2595:           else if ( wa1[j++] ) w[i] = wa[ia++];
                   2596:           else w[i] = wb[ib++];
                   2597:         for ( i = j = 0; i < l; i++ ) {
                   2598:           for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j);
                   2599:           NBM_CLR(wmb,j); j++;
                   2600:         }
                   2601:         wm->d = j; wm->c = c;
                   2602:         for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
                   2603:           tm = (NBM)BDY(t);
                   2604:           s = comp_nbm(tm,wm);
                   2605:           if ( s < 0 ) {
                   2606:             /* insert */
                   2607:             MKNODE(t1,wm,t);
                   2608:             if ( !p ) r = t1;
                   2609:             else NEXT(p) = t1;
                   2610:             NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   2611:             break;
                   2612:           } else if ( s == 0 ) {
                   2613:             /* add coefs */
                   2614:             addp(CO,tm->c,c,&c1);
                   2615:             if ( c1 ) tm->c = c1;
                   2616:             else NEXT(p) = NEXT(t);
                   2617:             break;
                   2618:           }
                   2619:         }
                   2620:         if ( !t ) {
                   2621:           /* append */
                   2622:           MKNODE(t1,wm,t);
                   2623:           if ( !p ) r = t1;
                   2624:           else NEXT(p) = t1;
                   2625:           NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   2626:         }
                   2627:       } while ( ni_next(wa1,lab1) );
                   2628:     } while ( ni_next(wab,l) );
                   2629:   }
                   2630:   MKNBP(u,r);
                   2631:   return u;
                   2632: }
                   2633: #endif
                   2634:
                   2635: /* DPM functions */
                   2636:
1.3       noro     2637: DMMstack dmm_stack;
1.7       noro     2638: extern LIST schreyer_obj;
1.3       noro     2639:
1.4       noro     2640: void push_schreyer_order(LIST data)
1.3       noro     2641: {
                   2642:   DMMstack t;
                   2643:   int len,i;
1.7       noro     2644:   NODE in,t1;
1.4       noro     2645:
                   2646:   /* data = [DPM,...,DPM] */
                   2647:   in = BDY(data);
                   2648:   len = length(in);
                   2649:   NEWDMMstack(t);
                   2650:   t->rank = len;
                   2651:   t->in = (DMM *)MALLOC((len+1)*sizeof(DMM));
                   2652:   t->ordtype = 0;
                   2653:   for ( i = 1; i <= len; i++, in = NEXT(in) ) {
                   2654:     t->in[i] = BDY((DPM)BDY(in));
                   2655:   }
                   2656:   t->next = dmm_stack;
                   2657:   dmm_stack = t;
                   2658:   dpm_ordtype = 2;
1.7       noro     2659:   MKNODE(t1,data,schreyer_obj?BDY(schreyer_obj):0);
                   2660:   MKLIST(schreyer_obj,t1);
1.4       noro     2661: }
                   2662:
                   2663: // data=[Ink,...,In0]
                   2664: // Ini = a list of module monomials
                   2665:
                   2666: void set_schreyer_order(LIST data)
                   2667: {
                   2668:   NODE in;
                   2669:   LIST *w;
                   2670:   int i,len;
1.3       noro     2671:
                   2672:   if ( !data ) {
                   2673:     dmm_stack = 0;
                   2674:     if ( dp_current_spec && dp_current_spec->id >= 256 )
                   2675:       dpm_ordtype = dp_current_spec->ispot;
                   2676:     else
                   2677:       dpm_ordtype = 0;
                   2678:     return;
                   2679:   } else {
1.4       noro     2680:     dmm_stack = 0;
                   2681:     in = BDY(data);
1.3       noro     2682:     len = length(in);
1.4       noro     2683:     w = (LIST *)MALLOC(len*sizeof(LIST));
                   2684:     for ( i = 0; i < len; i++, in = NEXT(in) ) w[i] = (LIST)BDY(in);
                   2685:     for ( i = len-1; i >= 0; i-- ) push_schreyer_order(w[i]);
                   2686:     dpm_ordtype = 2;
1.3       noro     2687:   }
                   2688: }
                   2689:
1.4       noro     2690: // construct a base of syz(g)
                   2691: // assuming the schrerer order is properly set
                   2692:
                   2693: DP dpm_sp_hm(DPM p1,DPM p2);
                   2694: void dpm_sp(DPM p1,DPM p2,DPM *sp,DP *t1,DP *t2);
                   2695: DP *dpm_nf_and_quotient(NODE b,DPM sp,VECT psv,DPM *nf,P *dn);
                   2696: void dpm_sort(DPM p,DPM *r);
                   2697:
                   2698: extern int DP_Multiple;
                   2699:
                   2700: void dpm_nf_z(NODE b,DPM g,VECT psv,int full,int multiple,DPM *rp);
                   2701: NODE dpm_sort_list(NODE l);
                   2702: void dpm_ptozp(DPM p,Z *cont,DPM *r);
                   2703:
                   2704: NODE dpm_reduceall(NODE in)
                   2705: {
                   2706:   int n,i;
                   2707:   VECT psv;
                   2708:   DPM *ps;
                   2709:   NODE t,t1;
                   2710:   DPM g,r;
                   2711:   Z cont;
                   2712:
                   2713:   n = length(in);
                   2714:   MKVECT(psv,n);
                   2715:   ps = (DPM *)BDY(psv);
                   2716:   for ( i = 0, t = in; i < n; i++, t = NEXT(t) ) ps[i] = BDY(t);
                   2717:   for ( i = 0; i < n; i++ ) {
                   2718:     g = ps[i]; ps[i] = 0;
1.8     ! noro     2719: //    dpm_nf_z(0,g,psv,1,DP_Multiple,&ps[i]);
        !          2720:     dpm_nf_z(0,g,psv,1,0,&ps[i]);
1.4       noro     2721:   }
                   2722:   t = 0;
                   2723:   for ( i = n-1; i >= 0; i-- ) {
                   2724:     dpm_ptozp(ps[i],&cont,&r);
                   2725:     MKNODE(t1,r,t); t = t1;
                   2726:   }
                   2727:   return t;
                   2728: }
                   2729:
1.8     ! noro     2730: struct oEGT egra;
        !          2731:
1.4       noro     2732: void dpm_schreyer_base(LIST g,LIST *s)
                   2733: {
                   2734:   NODE nd,t,b0,b;
                   2735:   int n,i,j,k,nv;
                   2736:   Z cont;
                   2737:   P dn,c;
                   2738:   DP h,t1,t2;
                   2739:   MP d;
                   2740:   DMM r0,r;
                   2741:   DPM sp,nf,dpm;
                   2742:   DPM *ps;
                   2743:   VECT psv;
                   2744:   DP **m,*quo;
1.8     ! noro     2745:   struct oEGT eg0,eg1;
        !          2746:   extern struct oEGT egred;
1.4       noro     2747:
1.8     ! noro     2748:   init_eg(&egra);
1.4       noro     2749:   nd = BDY(g);
                   2750:   n = length(nd);
                   2751:   MKVECT(psv,n);
                   2752:   ps = (DPM *)BDY(psv);
                   2753:   for ( i = 0, t = nd; i < n; i++, t = NEXT(t) ) ps[i] = (DPM)BDY(t);
                   2754:   nv = ps[0]->nv;
                   2755:   m = (DP **)almat_pointer(n,n);
                   2756:   b0 = 0;
                   2757:   for ( i = 0; i < n; i++ ) {
                   2758:     // sp(ps[i],ps[j]) = ti*ps[i]-tj*ps[j] => m[i][j] = ti
                   2759:     for ( j = i+1; j < n; j++ ) m[i][j] = dpm_sp_hm(ps[i],ps[j]);
                   2760:     for ( j = i+1; j < n; j++ ) {
                   2761:       if ( !m[i][j] ) continue;
                   2762:       for ( h = m[i][j], k = i+1; k < n; k++ )
                   2763:         if ( k != j && m[i][k] && dp_redble(m[i][k],h) ) m[i][k] = 0;
                   2764:     }
                   2765:     for ( j = i+1; j < n; j++ ) {
                   2766:       if ( m[i][j] ) {
                   2767:         dpm_sp(ps[i],ps[j],&sp,&t1,&t2);
                   2768:         quo = dpm_nf_and_quotient(0,sp,psv,&nf,&dn);
                   2769:         if ( nf )
                   2770:           error("dpm_schreyer_base : cannot happen");
                   2771:         NEWDMM(r0); r = r0;
                   2772:         mulp(CO,(P)BDY(t1)->c,dn,(P *)&r->c); r->pos = i+1; r->dl = BDY(t1)->dl;
                   2773:         NEWDMM(NEXT(r)); r=NEXT(r);
                   2774:         mulp(CO,(P)BDY(t2)->c,dn,&c); chsgnp(c,(P *)&r->c); r->pos = j+1; r->dl = BDY(t2)->dl;
                   2775:         if ( quo ) {
                   2776:           for ( k = 0; k < n; k++ ) {
                   2777:             if ( !quo[k] ) continue;
                   2778:             for ( d = BDY(quo[k]); d; d = NEXT(d) ) {
                   2779:               NEXTDMM(r0,r); chsgnp((P)d->c,(P *)&r->c); r->pos = k+1; r->dl = d->dl;
                   2780:             }
                   2781:           }
                   2782:         }
                   2783:         NEXT(r) = 0;
                   2784:         MKDPM(nv,r0,dpm); // XXX : sugar is not set
                   2785:         NEXTNODE(b0,b);
                   2786:         BDY(b) = (pointer)dpm;
                   2787:       }
                   2788:     }
                   2789:     if ( b0 ) NEXT(b) = 0;
                   2790:   }
                   2791:   push_schreyer_order(g);
                   2792:   for ( t = b0; t; t = NEXT(t) ) {
                   2793:     dpm_sort((DPM)BDY(t),&dpm);
                   2794:     BDY(t) = (pointer)dpm;
                   2795:   }
                   2796:   b0 = dpm_sort_list(b0);
1.8     ! noro     2797: // get_eg(&eg0);
        !          2798:   b0 = dpm_reduceall(b0);
        !          2799: // get_eg(&eg1); add_eg(&egra,&eg0,&eg1); print_eg("RA",&egra);
1.4       noro     2800:   MKLIST(*s,b0);
1.8     ! noro     2801: //  print_eg("red",&egred); printf("\n");
1.4       noro     2802: }
                   2803:
1.3       noro     2804: int compdmm_schreyer(int n,DMM m1,DMM m2)
                   2805: {
                   2806:    int pos1,pos2,t;
                   2807:    DMM *in;
                   2808:    DMMstack s;
1.4       noro     2809:    static DL d1=0,d2=0;
                   2810:    static int dlen=0;
1.3       noro     2811:
1.4       noro     2812:    pos1 = m1->pos; pos2 = m2->pos;
                   2813:    if ( pos1 == pos2 ) return (*cmpdl)(n,m1->dl,m2->dl);
                   2814:    if ( n > dlen ) {
                   2815:      NEWDL(d1,n); NEWDL(d2,n); dlen = n;
                   2816:    }
                   2817:    _copydl(n,m1->dl,d1);
                   2818:    _copydl(n,m2->dl,d2);
                   2819:    for ( s = dmm_stack; s; s = NEXT(s) ) {
1.3       noro     2820:      in = s->in;
                   2821:      _addtodl(n,in[pos1]->dl,d1);
                   2822:      _addtodl(n,in[pos2]->dl,d2);
1.4       noro     2823:      if ( in[pos1]->pos == in[pos2]->pos && _eqdl(n,d1,d2)) {
1.3       noro     2824:        if ( pos1 < pos2 ) return 1;
                   2825:        else if ( pos1 > pos2 ) return -1;
                   2826:        else return 0;
                   2827:      }
                   2828:      pos1 = in[pos1]->pos;
                   2829:      pos2 = in[pos2]->pos;
1.4       noro     2830:      if ( pos1 == pos2 ) return (*cmpdl)(n,d1,d2);
1.3       noro     2831:    }
                   2832:    // comparison by the bottom order
1.4       noro     2833: LAST:
                   2834:   if ( dpm_ordtype == 1 ) {
1.3       noro     2835:     if ( pos1 < pos2 ) return 1;
                   2836:     else if ( pos1 > pos2 ) return -1;
                   2837:     else return (*cmpdl)(n,d1,d2);
                   2838:   } else {
                   2839:     t = (*cmpdl)(n,d1,d2);
                   2840:     if ( t ) return t;
                   2841:     else if ( pos1 < pos2 ) return 1;
                   2842:     else if ( pos1 > pos2 ) return -1;
                   2843:     else return 0;
                   2844:   }
                   2845: }
                   2846:
1.4       noro     2847: #if 1
                   2848: int compdmm(int n,DMM m1,DMM m2)
                   2849: {
                   2850:   int t;
                   2851:
                   2852:   switch ( dpm_ordtype ) {
                   2853:   case 0:
                   2854:     t = (*cmpdl)(n,m1->dl,m2->dl);
                   2855:     if ( t ) return t;
                   2856:     else if ( m1->pos < m2->pos ) return 1;
                   2857:     else if ( m1->pos > m2->pos ) return -1;
                   2858:     else return 0;
                   2859:   case 1:
                   2860:     if ( m1->pos < m2->pos ) return 1;
                   2861:     else if ( m1->pos > m2->pos ) return -1;
                   2862:     else return (*cmpdl)(n,m1->dl,m2->dl);
                   2863:   case 2:
                   2864:     return compdmm_schreyer(n,m1,m2);
                   2865:   default:
                   2866:     error("compdmm : invalid dpm_ordtype");
                   2867:   }
                   2868: }
                   2869: #else
1.1       noro     2870: int compdmm(int n,DMM m1,DMM m2)
                   2871: {
                   2872:   int t;
                   2873:
1.3       noro     2874:   if ( dpm_ordtype == 1 ) {
1.1       noro     2875:     if ( m1->pos < m2->pos ) return 1;
                   2876:     else if ( m1->pos > m2->pos ) return -1;
                   2877:     else return (*cmpdl)(n,m1->dl,m2->dl);
1.3       noro     2878:   } else if ( dpm_ordtype == 0 ) {
1.1       noro     2879:     t = (*cmpdl)(n,m1->dl,m2->dl);
                   2880:     if ( t ) return t;
                   2881:     else if ( m1->pos < m2->pos ) return 1;
                   2882:     else if ( m1->pos > m2->pos ) return -1;
                   2883:     else return 0;
1.3       noro     2884:   } else if ( dpm_ordtype == 2 ) {
                   2885:     return compdmm_schreyer(n,m1,m2);
1.1       noro     2886:   }
                   2887: }
1.4       noro     2888: #endif
1.1       noro     2889:
                   2890: void adddpm(VL vl,DPM p1,DPM p2,DPM *pr)
                   2891: {
1.4       noro     2892:   int n,s;
1.1       noro     2893:   DMM m1,m2,mr=0,mr0;
                   2894:   Obj t;
                   2895:   DL d;
                   2896:
                   2897:   if ( !p1 )
                   2898:     *pr = p2;
                   2899:   else if ( !p2 )
                   2900:     *pr = p1;
                   2901:   else {
1.4       noro     2902:     for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
                   2903:       s = compdmm(n,m1,m2);
                   2904:       switch ( s ) {
1.1       noro     2905:         case 0:
                   2906:           arf_add(vl,C(m1),C(m2),&t);
                   2907:           if ( t ) {
                   2908:             NEXTDMM(mr0,mr); mr->pos = m1->pos; mr->dl = m1->dl; C(mr) = t;
                   2909:           }
                   2910:           m1 = NEXT(m1); m2 = NEXT(m2); break;
                   2911:         case 1:
                   2912:           NEXTDMM(mr0,mr); mr->pos = m1->pos; mr->dl = m1->dl; C(mr) = C(m1);
                   2913:           m1 = NEXT(m1); break;
                   2914:         case -1:
                   2915:           NEXTDMM(mr0,mr); mr->pos = m2->pos; mr->dl = m2->dl; C(mr) = C(m2);
                   2916:           m2 = NEXT(m2); break;
                   2917:       }
1.4       noro     2918:     }
1.1       noro     2919:     if ( !mr0 )
                   2920:       if ( m1 )
                   2921:         mr0 = m1;
                   2922:       else if ( m2 )
                   2923:         mr0 = m2;
                   2924:       else {
                   2925:         *pr = 0;
                   2926:         return;
                   2927:       }
                   2928:     else if ( m1 )
                   2929:       NEXT(mr) = m1;
                   2930:     else if ( m2 )
                   2931:       NEXT(mr) = m2;
                   2932:     else
                   2933:       NEXT(mr) = 0;
                   2934:     MKDPM(NV(p1),mr0,*pr);
                   2935:     if ( *pr )
                   2936:       (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                   2937:   }
                   2938: }
                   2939:
                   2940: void subdpm(VL vl,DPM p1,DPM p2,DPM *pr)
                   2941: {
                   2942:   DPM t;
                   2943:
                   2944:   if ( !p2 )
                   2945:     *pr = p1;
                   2946:   else {
                   2947:     chsgndpm(p2,&t); adddpm(vl,p1,t,pr);
                   2948:   }
                   2949: }
                   2950:
                   2951: void chsgndpm(DPM p,DPM *pr)
                   2952: {
                   2953:   DMM m,mr=0,mr0;
                   2954:   Obj r;
                   2955:
                   2956:   if ( !p )
                   2957:     *pr = 0;
                   2958:   else {
                   2959:     for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   2960:       NEXTDMM(mr0,mr); arf_chsgn(C(m),&C(mr)); mr->pos = m->pos; mr->dl = m->dl;
                   2961:     }
                   2962:     NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
                   2963:     if ( *pr )
                   2964:       (*pr)->sugar = p->sugar;
                   2965:   }
                   2966: }
                   2967:
1.8     ! noro     2968: void mulcmp(Obj c,MP m)
        !          2969: {
        !          2970:   MP t;
        !          2971:   Obj c1;
        !          2972:
        !          2973:   for ( t = m; t; t = NEXT(t) ) {
        !          2974:     arf_mul(CO,c,C(t),&c1); C(t) = c1;
        !          2975:   }
        !          2976: }
        !          2977:
        !          2978: void mulcdmm(Obj c,DMM m)
        !          2979: {
        !          2980:   DMM t;
        !          2981:   Obj c1;
        !          2982:
        !          2983:   for ( t = m; t; t = NEXT(t) ) {
        !          2984:     arf_mul(CO,c,C(t),&c1); C(t) = c1;
        !          2985:   }
        !          2986: }
        !          2987:
1.1       noro     2988: void mulcdpm(VL vl,Obj c,DPM p,DPM *pr)
                   2989: {
                   2990:   DMM m,mr=0,mr0;
                   2991:
                   2992:   if ( !p || !c )
                   2993:     *pr = 0;
                   2994:   else if ( NUM(c) && UNIQ((Q)c) )
                   2995:     *pr = p;
                   2996:   else if ( NUM(c) && MUNIQ((Q)c) )
                   2997:     chsgndpm(p,pr);
                   2998:   else {
                   2999:     for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   3000:       NEXTDMM(mr0,mr);
                   3001:       arf_mul(vl,C(m),c,&C(mr));
                   3002:       mr->pos = m->pos;
                   3003:       mr->dl = m->dl;
                   3004:     }
                   3005:     NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
                   3006:     if ( *pr )
                   3007:       (*pr)->sugar = p->sugar;
                   3008:   }
                   3009: }
                   3010:
                   3011: void comm_mulmpdpm(VL vl,MP m0,DPM p,DPM *pr)
                   3012: {
                   3013:   DMM m,mr=0,mr0;
                   3014:   DL d;
                   3015:   Obj c;
                   3016:   int n;
                   3017:
                   3018:   if ( !p )
                   3019:     *pr = 0;
                   3020:   else {
                   3021:     n = NV(p);
                   3022:     d = m0->dl;
                   3023:     c = C(m0);
                   3024:     for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   3025:       NEXTDMM(mr0,mr);
                   3026:       arf_mul(vl,C(m),c,&C(mr));
                   3027:       mr->pos = m->pos;
                   3028:       adddl(n,m->dl,d,&mr->dl);
                   3029:     }
                   3030:     NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
                   3031:     if ( *pr )
                   3032:       (*pr)->sugar = p->sugar;
                   3033:   }
                   3034: }
                   3035:
                   3036: void weyl_mulmpdpm(VL vl,MP m0,DPM p,DPM *pr)
                   3037: {
                   3038:   DPM r,t,t1;
                   3039:   DMM m;
                   3040:   DL d0;
                   3041:   int n,n2,l,i,j,tlen;
                   3042:   struct oMP mp;
                   3043:   static DMM *w,*psum;
                   3044:   static struct cdl *tab;
                   3045:   static int wlen;
                   3046:   static int rtlen;
                   3047:
                   3048:   if ( !p )
                   3049:     *pr = 0;
                   3050:   else {
                   3051:     for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
                   3052:     if ( l > wlen ) {
                   3053:       if ( w ) GCFREE(w);
                   3054:       w = (DMM *)MALLOC(l*sizeof(DMM));
                   3055:       wlen = l;
                   3056:     }
                   3057:     for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
                   3058:       w[i] = m;
                   3059:
                   3060:     n = NV(p); n2 = n>>1;
                   3061:     d0 = m0->dl;
                   3062:     for ( i = 0, tlen = 1; i < n2; i++ )
                   3063:       tlen *= d0->d[n2+i]+1;
                   3064:     if ( tlen > rtlen ) {
                   3065:       if ( tab ) GCFREE(tab);
                   3066:       if ( psum ) GCFREE(psum);
                   3067:       rtlen = tlen;
                   3068:       tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
                   3069:       psum = (DMM *)MALLOC(rtlen*sizeof(DMM));
                   3070:     }
                   3071:     bzero(psum,tlen*sizeof(DMM));
                   3072:     for ( i = l-1; i >= 0; i-- ) {
                   3073:       bzero(tab,tlen*sizeof(struct cdl));
                   3074:       mp.dl = w[i]->dl; mp.c = C(w[i]); mp.next = 0;
                   3075:       weyl_mulmm(vl,m0,&mp,n,tab,tlen);
                   3076:       for ( j = 0; j < tlen; j++ ) {
                   3077:         if ( tab[j].c ) {
                   3078:           NEWDMM(m); m->dl = tab[j].d; m->pos = w[i]->pos; C(m) = (Obj)tab[j].c; NEXT(m) = psum[j];
                   3079:           psum[j] = m;
                   3080:         }
                   3081:       }
                   3082:     }
                   3083:     for ( j = tlen-1, r = 0; j >= 0; j-- )
                   3084:       if ( psum[j] ) {
                   3085:         MKDPM(n,psum[j],t); adddpm(vl,r,t,&t1); r = t1;
                   3086:       }
                   3087:     if ( r )
                   3088:       r->sugar = p->sugar + m0->dl->td;
                   3089:     *pr = r;
                   3090:   }
                   3091: }
                   3092:
                   3093: void mulobjdpm(VL vl,Obj p1,DPM p2,DPM *pr)
                   3094: {
                   3095:   MP m;
                   3096:   DPM s,t,u;
                   3097:
                   3098:   if ( !p1 || !p2 )
                   3099:     *pr = 0;
                   3100:   else if ( OID(p1) != O_DP )
                   3101:     mulcdpm(vl,p1,p2,pr);
                   3102:   else {
                   3103:     s = 0;
                   3104:     for ( m = BDY((DP)p1); m; m = NEXT(m) ) {
                   3105:       if ( do_weyl )
                   3106:         weyl_mulmpdpm(vl,m,p2,&t);
                   3107:       else
                   3108:         comm_mulmpdpm(vl,m,p2,&t);
                   3109:       adddpm(vl,s,t,&u); s = u;
                   3110:     }
                   3111:     *pr = s;
                   3112:   }
                   3113: }
                   3114:
                   3115: int compdpm(VL vl,DPM p1,DPM p2)
                   3116: {
                   3117:   int n,t;
                   3118:   DMM m1,m2;
                   3119:
                   3120:   if ( !p1 )
                   3121:     return p2 ? -1 : 0;
                   3122:   else if ( !p2 )
                   3123:     return 1;
                   3124:   else if ( NV(p1) != NV(p2) ) {
                   3125:     error("compdpm : size mismatch");
                   3126:     return 0; /* XXX */
                   3127:   } else {
                   3128:     for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
                   3129:       m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
                   3130:       if ( (t = compdmm(n,m1,m2)) ||
                   3131:         (t = arf_comp(vl,C(m1),C(m2)) ) )
                   3132:         return t;
                   3133:     if ( m1 )
                   3134:       return 1;
                   3135:     else if ( m2 )
                   3136:       return -1;
                   3137:     else
                   3138:       return 0;
                   3139:   }
                   3140: }
                   3141:
1.5       noro     3142: // p = ...+c*<<0,...0:pos>>+...
                   3143: DPM dpm_eliminate_term(DPM a,DPM p,Obj c,int pos)
                   3144: {
                   3145:   MP d0,d;
                   3146:   DMM m;
                   3147:   DP f;
                   3148:   DPM a1,p1,r;
                   3149:
                   3150:   if ( !a ) return 0;
                   3151:   d0 = 0;
                   3152:   for ( m = BDY(a); m; m = NEXT(m) )
                   3153:     if ( m->pos == pos ) {
                   3154:       NEXTMP(d0,d); d->dl = m->dl; arf_chsgn(m->c,&d->c);
                   3155:     }
                   3156:   if ( d0 ) {
1.6       noro     3157:     NEXT(d) = 0; MKDP(NV(a),d0,f);
                   3158:     mulcdpm(CO,c,a,&a1);
                   3159:     mulobjdpm(CO,(Obj)f,p,&p1);
                   3160:     adddpm(CO,a1,p1,&r);
                   3161:     return r;
1.5       noro     3162:   } else
1.6       noro     3163:     return a;
1.5       noro     3164: }
                   3165:
                   3166: // <<...:i>> -> <<...:tab[i]>>
                   3167: DPM dpm_compress(DPM p,int *tab)
                   3168: {
                   3169:   DMM m,mr0,mr;
                   3170:   DPM t;
                   3171:
                   3172:   if ( !p ) return 0;
                   3173:   else {
                   3174:     for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
                   3175:       NEXTDMM(mr0,mr);
1.6       noro     3176:       mr->dl = m->dl; mr->c = m->c; mr->pos = tab[m->pos];
1.5       noro     3177:       if ( mr->pos <= 0 )
                   3178:         error("dpm_compress : invalid position");
                   3179:     }
                   3180:     NEXT(mr) = 0;
                   3181:     MKDPM(p->nv,mr0,t); t->sugar = p->sugar;
                   3182:     return t;
                   3183:   }
                   3184: }
                   3185:
                   3186: // input : s, s = syz(m) output simplified s, m
1.6       noro     3187: void dpm_simplify_syz(LIST s,LIST m,LIST *s1,LIST *m1)
1.5       noro     3188: {
                   3189:   int lm,ls,i,j,pos;
                   3190:   DPM *am,*as;
                   3191:   DPM p;
                   3192:   DMM d;
                   3193:   Obj c;
                   3194:   int *tab;
                   3195:   NODE t,t1;
                   3196:
                   3197:   lm = length(BDY(m));
                   3198:   am = (DPM *)MALLOC((lm+1)*sizeof(DPM));
                   3199:   ls = length(BDY(s));
                   3200:   as = (DPM *)MALLOC(ls*sizeof(DPM));
                   3201:   for ( i = 1, t = BDY(m); i <= lm; i++, t = NEXT(t) ) am[i] = (DPM)BDY(t);
1.6       noro     3202:   for ( i = 0, t = BDY(s); i < ls; i++, t = NEXT(t) ) as[i] = (DPM)BDY(t);
1.5       noro     3203:
                   3204:   for ( i = 0; i < ls; i++ ) {
                   3205:     p = as[i];
1.6       noro     3206:     if ( p == 0 ) continue;
1.5       noro     3207:     for ( d = BDY(p); d; d = NEXT(d) ) if ( d->dl->td == 0 ) break;
                   3208:     if ( d ) {
                   3209:       c = d->c; pos = d->pos;
                   3210:       for ( j = 0; j < ls; j++ )
1.6       noro     3211:         if ( j != i ) {
1.5       noro     3212:           as[j] = dpm_eliminate_term(as[j],p,c,pos);
1.6       noro     3213:         }
1.5       noro     3214:       // remove m[i]
                   3215:       am[pos] = 0;
                   3216:       as[i] = 0;
                   3217:     }
                   3218:   }
                   3219:   // compress s
                   3220:   // create index table from am[]
                   3221:   // (0 0 * 0 * ...) -> (0 0 1 0 2 ... ) which means 2->1, 4->2, ...
                   3222:   tab = (int *)MALLOC((lm+1)*sizeof(int));
                   3223:   for ( j = 0, i = 1; i <= lm; i++ ) {
                   3224:     if ( am[i] ) { j++; tab[i] = j; }
                   3225:     else tab[i] = 0;
                   3226:   }
                   3227:   t = 0;
                   3228:   for ( i = ls-1; i >= 0; i-- )
                   3229:     if ( as[i] ) {
                   3230:       p = dpm_compress(as[i],tab);
                   3231:       MKNODE(t1,(pointer)p,t); t = t1;
                   3232:     }
                   3233:   MKLIST(*s1,t);
                   3234:   t = 0;
                   3235:   for ( i = lm; i >= 1; i-- )
                   3236:     if ( am[i] ) {
                   3237:       MKNODE(t1,(pointer)am[i],t); t = t1;
                   3238:     }
                   3239:   MKLIST(*m1,t);
                   3240: }

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