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

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.26    ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2018/engine/dist.c,v 1.25 2021/01/11 08:37:44 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:
1.26    ! noro      213: #if 0
1.1       noro      214: void ptod(VL vl,VL dvl,P p,DP *pr)
                    215: {
                    216:   int n,i,j,k;
                    217:   VL tvl;
                    218:   V v;
                    219:   DL d;
                    220:   MP m;
                    221:   DCP dc;
                    222:   DCP *w;
                    223:   DP r,s,t,u;
                    224:   P x,c;
                    225:
                    226:   if ( !p )
                    227:     *pr = 0;
                    228:   else if ( OID(p) > O_P )
                    229:     error("ptod : only polynomials can be converted.");
                    230:   else {
                    231:     for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
                    232:     if ( NUM(p) ) {
                    233:       NEWDL(d,n);
                    234:       NEWMP(m); m->dl = d; C(m) = (Obj)p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
                    235:     } else {
                    236:       for ( i = 0, tvl = dvl, v = VR(p);
                    237:         tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
                    238:       if ( !tvl ) {
                    239:         for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
                    240:         w = (DCP *)ALLOCA(k*sizeof(DCP));
                    241:         for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
                    242:           w[j] = dc;
                    243:
                    244:         for ( j = k-1, s = 0, MKV(v,x); j >= 0; j-- ) {
                    245:           ptod(vl,dvl,COEF(w[j]),&t); pwrp(vl,x,DEG(w[j]),&c);
                    246:           muldc(vl,t,(Obj)c,&r); addd(vl,r,s,&t); s = t;
                    247:         }
                    248:         *pr = s;
                    249:       } else {
                    250:         for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
                    251:         w = (DCP *)ALLOCA(k*sizeof(DCP));
                    252:         for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
                    253:           w[j] = dc;
                    254:
                    255:         for ( j = k-1, s = 0; j >= 0; j-- ) {
                    256:           ptod(vl,dvl,COEF(w[j]),&t);
1.2       noro      257:           NEWDL(d,n); d->d[i] = ZTOS(DEG(w[j]));
1.1       noro      258:           d->td = MUL_WEIGHT(d->d[i],i);
                    259:           NEWMP(m); m->dl = d; C(m) = (Obj)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
                    260:           comm_muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
                    261:         }
                    262:         *pr = s;
                    263:       }
                    264:     }
                    265:   }
                    266: #if 0
                    267:   if ( !dp_fcoeffs && has_sfcoef(*pr) )
                    268:     dp_fcoeffs = N_GFS;
                    269: #endif
                    270: }
1.26    ! noro      271: #else
        !           272: void ptod(VL vl,VL dvl,P p,DP *pr)
        !           273: {
        !           274:   int n,i,j,k;
        !           275:   VL tvl;
        !           276:   V v;
        !           277:   DL d;
        !           278:   MP m;
        !           279:   DCP dc;
        !           280:   DP *y;
        !           281:   DP r,s,t,u;
        !           282:   P x,c;
        !           283:
        !           284:   if ( !p )
        !           285:     *pr = 0;
        !           286:   else if ( OID(p) > O_P )
        !           287:     error("ptod : only polynomials can be converted.");
        !           288:   else {
        !           289:     for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
        !           290:     if ( NUM(p) ) {
        !           291:       NEWDL(d,n);
        !           292:       NEWMP(m); m->dl = d; C(m) = (Obj)p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
        !           293:     } else {
        !           294:       for ( i = 0, tvl = dvl, v = VR(p); tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
        !           295:       for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
        !           296:       y = (DP *)ALLOCA(k*sizeof(DP));
        !           297:       if ( !tvl ) {
        !           298:         MKV(v,x);
        !           299:         for ( dc = DC(p), j = 0; dc; dc = NEXT(dc), j++ ) {
        !           300:           ptod(vl,dvl,COEF(dc),&t); pwrp(vl,x,DEG(dc),&c);
        !           301:           muldc(vl,t,(Obj)c,&y[j]);
        !           302:         }
        !           303:       } else {
        !           304:         for ( dc = DC(p), j = 0; dc; dc = NEXT(dc), j++ ) {
        !           305:           ptod(vl,dvl,COEF(dc),&t);
        !           306:           NEWDL(d,n);
        !           307:           d->d[i] = ZTOS(DEG(dc));
        !           308:           d->td = MUL_WEIGHT(d->d[i],i);
        !           309:           NEWMP(m); m->dl = d; C(m) = (Obj)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
        !           310:           comm_muld(vl,t,u,&y[j]);
        !           311:         }
        !           312:       }
        !           313:       for ( j = k-1, s = 0; j >= 0; j-- ) {
        !           314:         addd(vl,y[j],s,&t); s = t;
        !           315:       }
        !           316:       *pr = s;
        !           317:     }
        !           318:   }
        !           319: #if 0
        !           320:   if ( !dp_fcoeffs && has_sfcoef(*pr) )
        !           321:     dp_fcoeffs = N_GFS;
        !           322: #endif
        !           323: }
        !           324: #endif
1.1       noro      325:
                    326: void dtop(VL vl,VL dvl,DP p,Obj *pr)
                    327: {
                    328:   int n,i,j,k;
                    329:   DL d;
                    330:   MP m;
                    331:   MP *a;
                    332:   P r;
                    333:   Obj t,w,s,u;
                    334:   Z q;
                    335:   VL tvl;
                    336:
                    337:   if ( !p )
                    338:     *pr = 0;
                    339:   else {
                    340:     for ( k = 0, m = BDY(p); m; m = NEXT(m), k++ );
                    341:     a = (MP *)ALLOCA(k*sizeof(MP));
                    342:     for ( j = 0, m = BDY(p); j < k; m = NEXT(m), j++ )
                    343:       a[j] = m;
                    344:
                    345:     for ( n = p->nv, j = k-1, s = 0; j >= 0; j-- ) {
                    346:       m = a[j];
                    347:       t = C(m);
                    348:       if ( NUM(t) && NID((Num)t) == N_M ) {
                    349:         mptop((P)t,(P *)&u); t = u;
                    350:       }
                    351:       for ( i = 0, d = m->dl, tvl = dvl;
                    352:         i < n; tvl = NEXT(tvl), i++ ) {
1.2       noro      353:         MKV(tvl->v,r); STOZ(d->d[i],q); pwrp(vl,r,q,(P *)&u);
1.1       noro      354:         arf_mul(vl,t,(Obj)u,&w); t = w;
                    355:       }
                    356:       arf_add(vl,s,t,&u); s = u;
                    357:     }
                    358:     *pr = s;
                    359:   }
                    360: }
                    361:
                    362: void nodetod(NODE node,DP *dp)
                    363: {
                    364:   NODE t;
                    365:   int len,i,td;
                    366:   Q e;
                    367:   DL d;
                    368:   MP m;
                    369:   DP u;
                    370:
                    371:   for ( t = node, len = 0; t; t = NEXT(t), len++ );
                    372:   NEWDL(d,len);
                    373:   for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
                    374:     e = (Q)BDY(t);
                    375:     if ( !e )
                    376:       d->d[i] = 0;
                    377:     else if ( !NUM(e) || !RATN(e) || !INT(e) )
                    378:       error("nodetod : invalid input");
                    379:     else {
1.2       noro      380:       d->d[i] = ZTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
1.1       noro      381:     }
                    382:   }
                    383:   d->td = td;
                    384:   NEWMP(m); m->dl = d; C(m) = (Obj)ONE; NEXT(m) = 0;
                    385:   MKDP(len,m,u); u->sugar = td; *dp = u;
                    386: }
                    387:
                    388: void nodetodpm(NODE node,Obj pos,DPM *dp)
                    389: {
                    390:   NODE t;
1.9       noro      391:   int len,i,td,p;
1.1       noro      392:   Q e;
                    393:   DL d;
                    394:   DMM m;
                    395:   DPM u;
                    396:
                    397:   for ( t = node, len = 0; t; t = NEXT(t), len++ );
                    398:   NEWDL(d,len);
                    399:   for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
                    400:     e = (Q)BDY(t);
                    401:     if ( !e )
                    402:       d->d[i] = 0;
                    403:     else if ( !NUM(e) || !RATN(e) || !INT(e) )
                    404:       error("nodetodpm : invalid input");
                    405:     else {
1.2       noro      406:       d->d[i] = ZTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
1.1       noro      407:     }
                    408:   }
                    409:   d->td = td;
1.9       noro      410:   p = ZTOS((Q)pos);
1.19      noro      411:   if ( dp_current_spec->module_top_weight ) {
1.9       noro      412:     if ( p > dp_current_spec->module_rank )
                    413:       error("nodetodpm : inconsistent order spec");
                    414:     d->td += dp_current_spec->module_top_weight[p-1];
                    415:   }
                    416:   NEWDMM(m); m->dl = d; m->pos = p; C(m) = (Obj)ONE; NEXT(m) = 0;
1.1       noro      417:   MKDPM(len,m,u); u->sugar = td; *dp = u;
                    418: }
                    419:
                    420: void dtodpm(DP d,int pos,DPM *dp)
                    421: {
                    422:   DMM mr0,mr;
                    423:   MP m;
1.9       noro      424:   int shift;
1.1       noro      425:
                    426:   if ( !d ) *dp = 0;
                    427:   else {
1.9       noro      428:     shift = 0;
1.19      noro      429:     if ( dp_current_spec->module_top_weight ) {
1.9       noro      430:       if ( pos > dp_current_spec->module_rank )
                    431:         error("nodetodpm : inconsistent order spec");
                    432:       shift = dp_current_spec->module_top_weight[pos-1];
                    433:     }
1.1       noro      434:     for ( m = BDY(d), mr0 = 0; m; m = NEXT(m) ) {
                    435:       NEXTDMM(mr0,mr);
                    436:       mr->dl = m->dl;
1.9       noro      437:       mr->dl->td += shift;
1.1       noro      438:       mr->pos = pos;
                    439:       C(mr) = C(m);
                    440:     }
                    441:     MKDPM(d->nv,mr0,*dp); (*dp)->sugar = d->sugar;
                    442:   }
                    443: }
                    444:
                    445: int sugard(MP m)
                    446: {
                    447:   int s;
                    448:
                    449:   for ( s = 0; m; m = NEXT(m) )
                    450:     s = MAX(s,m->dl->td);
                    451:   return s;
                    452: }
                    453:
                    454: void addd(VL vl,DP p1,DP p2,DP *pr)
                    455: {
                    456:   int n;
                    457:   MP m1,m2,mr=0,mr0;
                    458:   Obj t;
                    459:   DL d;
                    460:
                    461:   if ( !p1 )
                    462:     *pr = p2;
                    463:   else if ( !p2 )
                    464:     *pr = p1;
                    465:   else {
                    466:     if ( OID(p1) <= O_R ) {
                    467:       n = NV(p2);  NEWDL(d,n);
                    468:       NEWMP(m1); m1->dl = d; C(m1) = (Obj)p1; NEXT(m1) = 0;
                    469:       MKDP(n,m1,p1); (p1)->sugar = 0;
                    470:     }
                    471:     if ( OID(p2) <= O_R ) {
                    472:       n = NV(p1);  NEWDL(d,n);
                    473:       NEWMP(m2); m2->dl = d; C(m2) = (Obj)p2; NEXT(m2) = 0;
                    474:       MKDP(n,m2,p2); (p2)->sugar = 0;
                    475:     }
                    476:     for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
                    477:       switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
                    478:         case 0:
                    479:           arf_add(vl,C(m1),C(m2),&t);
                    480:           if ( t ) {
                    481:             NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
                    482:           }
                    483:           m1 = NEXT(m1); m2 = NEXT(m2); break;
                    484:         case 1:
                    485:           NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
                    486:           m1 = NEXT(m1); break;
                    487:         case -1:
                    488:           NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
                    489:           m2 = NEXT(m2); break;
                    490:       }
                    491:     if ( !mr0 )
                    492:       if ( m1 )
                    493:         mr0 = m1;
                    494:       else if ( m2 )
                    495:         mr0 = m2;
                    496:       else {
                    497:         *pr = 0;
                    498:         return;
                    499:       }
                    500:     else if ( m1 )
                    501:       NEXT(mr) = m1;
                    502:     else if ( m2 )
                    503:       NEXT(mr) = m2;
                    504:     else
                    505:       NEXT(mr) = 0;
                    506:     MKDP(NV(p1),mr0,*pr);
                    507:     if ( *pr )
                    508:       (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                    509:   }
                    510: }
                    511:
                    512: /* for F4 symbolic reduction */
                    513:
                    514: void symb_addd(DP p1,DP p2,DP *pr)
                    515: {
                    516:   int n;
                    517:   MP m1,m2,mr=0,mr0;
                    518:
                    519:   if ( !p1 )
                    520:     *pr = p2;
                    521:   else if ( !p2 )
                    522:     *pr = p1;
                    523:   else {
                    524:     for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
                    525:       NEXTMP(mr0,mr); C(mr) = (Obj)ONE;
                    526:       switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
                    527:         case 0:
                    528:           mr->dl = m1->dl;
                    529:           m1 = NEXT(m1); m2 = NEXT(m2); break;
                    530:         case 1:
                    531:           mr->dl = m1->dl;
                    532:           m1 = NEXT(m1); break;
                    533:         case -1:
                    534:           mr->dl = m2->dl;
                    535:           m2 = NEXT(m2); break;
                    536:       }
                    537:     }
                    538:     if ( !mr0 )
                    539:       if ( m1 )
                    540:         mr0 = m1;
                    541:       else if ( m2 )
                    542:         mr0 = m2;
                    543:       else {
                    544:         *pr = 0;
                    545:         return;
                    546:       }
                    547:     else if ( m1 )
                    548:       NEXT(mr) = m1;
                    549:     else if ( m2 )
                    550:       NEXT(mr) = m2;
                    551:     else
                    552:       NEXT(mr) = 0;
                    553:     MKDP(NV(p1),mr0,*pr);
                    554:     if ( *pr )
                    555:       (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                    556:   }
                    557: }
                    558:
                    559: /*
                    560:  * destructive merge of two list
                    561:  *
                    562:  * p1, p2 : list of DL
                    563:  * return : a merged list
                    564:  */
                    565:
                    566: NODE symb_merge(NODE m1,NODE m2,int n)
                    567: {
                    568:   NODE top=0,prev,cur,m=0,t;
                    569:   DL d1,d2;
                    570:
                    571:   if ( !m1 )
                    572:     return m2;
                    573:   else if ( !m2 )
                    574:     return m1;
                    575:   else {
                    576:     switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
                    577:       case 0:
                    578:         top = m1; m = NEXT(m2);
                    579:         break;
                    580:       case 1:
                    581:         top = m1; m = m2;
                    582:         break;
                    583:       case -1:
                    584:         top = m2; m = m1;
                    585:         break;
                    586:     }
                    587:     prev = top; cur = NEXT(top);
                    588:     /* BDY(prev) > BDY(m) always holds */
                    589:     while ( cur && m ) {
                    590:       d1 = (DL)BDY(cur);
                    591:       d2 = (DL)BDY(m);
                    592: #if 1
                    593:       switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
                    594: #else
                    595:       /* XXX only valid for DRL */
                    596:       if ( d1->td > d2->td )
                    597:         c = 1;
                    598:       else if ( d1->td < d2->td )
                    599:         c = -1;
                    600:       else {
                    601:         for ( i = n-1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
                    602:         if ( i < 0 )
                    603:           c = 0;
                    604:         else if ( d1->d[i] < d2->d[i] )
                    605:           c = 1;
                    606:         else
                    607:           c = -1;
                    608:       }
                    609:       switch ( c ) {
                    610: #endif
                    611:         case 0:
                    612:           m = NEXT(m);
                    613:           prev = cur; cur = NEXT(cur);
                    614:           break;
                    615:         case 1:
                    616:           t = NEXT(cur); NEXT(cur) = m; m = t;
                    617:           prev = cur; cur = NEXT(cur);
                    618:           break;
                    619:         case -1:
                    620:           NEXT(prev) = m; m = cur;
                    621:           prev = NEXT(prev); cur = NEXT(prev);
                    622:           break;
                    623:       }
                    624:     }
                    625:     if ( !cur )
                    626:       NEXT(prev) = m;
                    627:     return top;
                    628:   }
                    629: }
                    630:
                    631: void _adddl(int n,DL d1,DL d2,DL d3)
                    632: {
                    633:   int i;
                    634:
                    635:   d3->td = d1->td+d2->td;
                    636:   for ( i = 0; i < n; i++ )
                    637:     d3->d[i] = d1->d[i]+d2->d[i];
                    638: }
                    639:
1.25      noro      640: void _subdl(int n,DL d1,DL d2,DL d3)
                    641: {
                    642:   int i;
                    643:
                    644:   d3->td = d1->td-d2->td;
                    645:   for ( i = 0; i < n; i++ )
                    646:     d3->d[i] = d1->d[i]-d2->d[i];
                    647: }
                    648:
                    649:
1.3       noro      650: void _addtodl(int n,DL d1,DL d2)
                    651: {
                    652:   int i;
                    653:
                    654:   d2->td += d1->td;
                    655:   for ( i = 0; i < n; i++ )
                    656:     d2->d[i] += d1->d[i];
                    657: }
                    658:
1.23      noro      659: void _subfromdl(int n,DL d1,DL d2)
                    660: {
                    661:   int i;
                    662:
                    663:   d2->td -= d1->td;
                    664:   for ( i = 0; i < n; i++ )
                    665:     d2->d[i] -= d1->d[i];
                    666: }
                    667:
1.3       noro      668: void _copydl(int n,DL d1,DL d2)
                    669: {
                    670:   int i;
                    671:
                    672:   d2->td = d1->td;
                    673:   for ( i = 0; i < n; i++ )
                    674:     d2->d[i] = d1->d[i];
                    675: }
                    676:
                    677: int _eqdl(int n,DL d1,DL d2)
                    678: {
                    679:   int i;
                    680:
                    681:   if ( d2->td != d1->td ) return 0;
                    682:   for ( i = 0; i < n; i++ )
                    683:     if ( d2->d[i] != d1->d[i] ) return 0;
                    684:   return 1;
                    685: }
                    686:
1.1       noro      687: /* m1 <- m1 U dl*f, destructive */
                    688:
                    689: NODE mul_dllist(DL dl,DP f);
                    690:
                    691: NODE symb_mul_merge(NODE m1,DL dl,DP f,int n)
                    692: {
                    693:   NODE top,prev,cur,n1;
                    694:   DP g;
                    695:   DL t,s;
                    696:   MP m;
                    697:
                    698:   if ( !m1 )
                    699:     return mul_dllist(dl,f);
                    700:   else if ( !f )
                    701:     return m1;
                    702:   else {
                    703:     m = BDY(f);
                    704:     NEWDL_NOINIT(t,n);
                    705:     _adddl(n,m->dl,dl,t);
                    706:     top = m1; prev = 0; cur = m1;
                    707:     while ( m ) {
                    708:       switch ( (*cmpdl)(n,(DL)BDY(cur),t) ) {
                    709:         case 0:
                    710:           prev = cur; cur = NEXT(cur);
                    711:           if ( !cur ) {
                    712:             MKDP(n,m,g);
                    713:             NEXT(prev) = mul_dllist(dl,g);
                    714:             return top;
                    715:           }
                    716:           m = NEXT(m);
                    717:           if ( m ) _adddl(n,m->dl,dl,t);
                    718:           break;
                    719:         case 1:
                    720:           prev = cur; cur = NEXT(cur);
                    721:           if ( !cur ) {
                    722:             MKDP(n,m,g);
                    723:             NEXT(prev) = mul_dllist(dl,g);
                    724:             return top;
                    725:           }
                    726:           break;
                    727:         case -1:
                    728:           NEWDL_NOINIT(s,n);
                    729:           s->td = t->td;
                    730:           bcopy(t->d,s->d,n*sizeof(int));
                    731:           NEWNODE(n1);
                    732:           n1->body = (pointer)s;
                    733:           NEXT(n1) = cur;
                    734:           if ( !prev ) {
                    735:             top = n1; cur = n1;
                    736:           } else {
                    737:             NEXT(prev) = n1; prev = n1;
                    738:           }
                    739:           m = NEXT(m);
                    740:           if ( m ) _adddl(n,m->dl,dl,t);
                    741:           break;
                    742:       }
                    743:     }
                    744:     return top;
                    745:   }
                    746: }
                    747:
                    748: DLBUCKET symb_merge_bucket(DLBUCKET m1,DLBUCKET m2,int n)
                    749: {
                    750:   DLBUCKET top,prev,cur,m,t;
                    751:
                    752:   if ( !m1 )
                    753:     return m2;
                    754:   else if ( !m2 )
                    755:     return m1;
                    756:   else {
                    757:     if ( m1->td == m2->td ) {
                    758:       top = m1;
                    759:       BDY(top) = symb_merge(BDY(top),BDY(m2),n);
                    760:       m = NEXT(m2);
                    761:     } else if ( m1->td > m2->td ) {
                    762:       top = m1; m = m2;
                    763:     } else {
                    764:       top = m2; m = m1;
                    765:     }
                    766:     prev = top; cur = NEXT(top);
                    767:     /* prev->td > m->td always holds */
                    768:     while ( cur && m ) {
                    769:       if ( cur->td == m->td ) {
                    770:         BDY(cur) = symb_merge(BDY(cur),BDY(m),n);
                    771:         m = NEXT(m);
                    772:         prev = cur; cur = NEXT(cur);
                    773:       } else if ( cur->td > m->td ) {
                    774:         t = NEXT(cur); NEXT(cur) = m; m = t;
                    775:         prev = cur; cur = NEXT(cur);
                    776:       } else {
                    777:         NEXT(prev) = m; m = cur;
                    778:         prev = NEXT(prev); cur = NEXT(prev);
                    779:       }
                    780:     }
                    781:     if ( !cur )
                    782:       NEXT(prev) = m;
                    783:     return top;
                    784:   }
                    785: }
                    786:
                    787: void subd(VL vl,DP p1,DP p2,DP *pr)
                    788: {
                    789:   DP t;
                    790:
                    791:   if ( !p2 )
                    792:     *pr = p1;
                    793:   else {
                    794:     chsgnd(p2,&t); addd(vl,p1,t,pr);
                    795:   }
                    796: }
                    797:
                    798: void chsgnd(DP p,DP *pr)
                    799: {
                    800:   MP m,mr=0,mr0;
                    801:   Obj r;
                    802:
                    803:   if ( !p )
                    804:     *pr = 0;
                    805:   else if ( OID(p) <= O_R ) {
                    806:     arf_chsgn((Obj)p,&r); *pr = (DP)r;
                    807:   } else {
                    808:     for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    809:       NEXTMP(mr0,mr); arf_chsgn(C(m),&C(mr)); mr->dl = m->dl;
                    810:     }
                    811:     NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    812:     if ( *pr )
                    813:       (*pr)->sugar = p->sugar;
                    814:   }
                    815: }
                    816:
                    817: void muld(VL vl,DP p1,DP p2,DP *pr)
                    818: {
                    819:   if ( ! do_weyl )
                    820:     comm_muld(vl,p1,p2,pr);
                    821:   else
                    822:     weyl_muld(vl,p1,p2,pr);
                    823: }
                    824:
                    825: void comm_muld(VL vl,DP p1,DP p2,DP *pr)
                    826: {
                    827:   MP m;
                    828:   DP s,t,u;
                    829:   int i,l,l1;
                    830:   static MP *w;
                    831:   static int wlen;
                    832:
                    833:   if ( !p1 || !p2 )
                    834:     *pr = 0;
                    835:   else if ( OID(p1) != O_DP )
                    836:     muldc(vl,p2,(Obj)p1,pr);
                    837:   else if ( OID(p2) != O_DP )
                    838:     muldc(vl,p1,(Obj)p2,pr);
                    839:   else {
                    840:     for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
                    841:     for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
                    842:     if ( l1 < l ) {
                    843:       t = p1; p1 = p2; p2 = t;
                    844:       l = l1;
                    845:     }
                    846:     if ( l > wlen ) {
                    847:       if ( w ) GCFREE(w);
                    848:       w = (MP *)MALLOC(l*sizeof(MP));
                    849:       wlen = l;
                    850:     }
                    851:     for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
                    852:       w[i] = m;
                    853:     for ( s = 0, i = l-1; i >= 0; i-- ) {
                    854:       muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
                    855:     }
                    856:     bzero(w,l*sizeof(MP));
                    857:     *pr = s;
                    858:   }
                    859: }
                    860:
                    861: /* discard terms which is not a multiple of dl */
                    862:
                    863: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr)
                    864: {
                    865:   MP m;
                    866:   DP s,t,u;
                    867:   int i,l,l1;
                    868:   static MP *w;
                    869:   static int wlen;
                    870:
                    871:   if ( !p1 || !p2 )
                    872:     *pr = 0;
                    873:   else if ( OID(p1) != O_DP )
                    874:     muldc_trunc(vl,p2,(Obj)p1,dl,pr);
                    875:   else if ( OID(p2) != O_DP )
                    876:     muldc_trunc(vl,p1,(Obj)p2,dl,pr);
                    877:   else {
                    878:     for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
                    879:     for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
                    880:     if ( l1 < l ) {
                    881:       t = p1; p1 = p2; p2 = t;
                    882:       l = l1;
                    883:     }
                    884:     if ( l > wlen ) {
                    885:       if ( w ) GCFREE(w);
                    886:       w = (MP *)MALLOC(l*sizeof(MP));
                    887:       wlen = l;
                    888:     }
                    889:     for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
                    890:       w[i] = m;
                    891:     for ( s = 0, i = l-1; i >= 0; i-- ) {
                    892:       muldm_trunc(vl,p1,w[i],dl,&t); addd(vl,s,t,&u); s = u;
                    893:     }
                    894:     bzero(w,l*sizeof(MP));
                    895:     *pr = s;
                    896:   }
                    897: }
                    898:
                    899: void comm_quod(VL vl,DP p1,DP p2,DP *pr)
                    900: {
                    901:   MP m=0,m0;
                    902:   DP s,t;
                    903:   int i,n,sugar;
                    904:   DL d1,d2,d;
                    905:   Q a,b;
                    906:
                    907:   if ( !p2 )
                    908:     error("comm_quod : invalid input");
                    909:   if ( !p1 )
                    910:     *pr = 0;
                    911:   else {
                    912:     n = NV(p1);
                    913:     d2 = BDY(p2)->dl;
                    914:     m0 = 0;
                    915:     sugar = p1->sugar;
                    916:     while ( p1 ) {
                    917:       d1 = BDY(p1)->dl;
                    918:       NEWDL(d,n);
                    919:       d->td = d1->td - d2->td;
                    920:       for ( i = 0; i < n; i++ )
                    921:         d->d[i] = d1->d[i]-d2->d[i];
                    922:       NEXTMP(m0,m);
                    923:       m->dl = d;
                    924:       divq((Q)BDY(p1)->c,(Q)BDY(p2)->c,&a); chsgnq(a,&b);
                    925:       C(m) = (Obj)b;
                    926:       muldm_trunc(vl,p2,m,d2,&t);
                    927:       addd(vl,p1,t,&s); p1 = s;
                    928:       C(m) = (Obj)a;
                    929:     }
                    930:     if ( m0 ) {
                    931:       NEXT(m) = 0; MKDP(n,m0,*pr);
                    932:     } else
                    933:       *pr = 0;
                    934:     /* XXX */
                    935:     if ( *pr )
                    936:       (*pr)->sugar = sugar - d2->td;
                    937:   }
                    938: }
                    939:
                    940: void muldm(VL vl,DP p,MP m0,DP *pr)
                    941: {
                    942:   MP m,mr=0,mr0;
                    943:   Obj c;
                    944:   DL d;
                    945:   int n;
                    946:
                    947:   if ( !p )
                    948:     *pr = 0;
                    949:   else {
                    950:     for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
                    951:       m; m = NEXT(m) ) {
                    952:       NEXTMP(mr0,mr);
                    953:       if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                    954:         mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                    955:       else
                    956:         arf_mul(vl,C(m),c,&C(mr));
                    957:       adddl(n,m->dl,d,&mr->dl);
                    958:     }
                    959:     NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    960:     if ( *pr )
                    961:       (*pr)->sugar = p->sugar + m0->dl->td;
                    962:   }
                    963: }
                    964:
                    965: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr)
                    966: {
                    967:   MP m,mr=0,mr0;
                    968:   Obj c;
                    969:   DL d,tdl;
                    970:   int n,i;
                    971:
                    972:   if ( !p )
                    973:     *pr = 0;
                    974:   else {
                    975:     n = NV(p);
                    976:     NEWDL(tdl,n);
                    977:     for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl;
                    978:       m; m = NEXT(m) ) {
                    979:       _adddl(n,m->dl,d,tdl);
                    980:       for ( i = 0; i < n; i++ )
                    981:         if ( tdl->d[i] < dl->d[i] )
                    982:           break;
                    983:       if ( i < n )
                    984:         continue;
                    985:       NEXTMP(mr0,mr);
                    986:       mr->dl = tdl;
                    987:       NEWDL(tdl,n);
                    988:       if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                    989:         mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                    990:       else
                    991:         arf_mul(vl,C(m),(Obj)c,&C(mr));
                    992:     }
                    993:     if ( mr0 ) {
                    994:       NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    995:     } else
                    996:       *pr = 0;
                    997:     if ( *pr )
                    998:       (*pr)->sugar = p->sugar + m0->dl->td;
                    999:   }
                   1000: }
                   1001:
                   1002: void weyl_muld(VL vl,DP p1,DP p2,DP *pr)
                   1003: {
                   1004:   MP m;
                   1005:   DP s,t,u;
                   1006:   int i,l;
                   1007:   static MP *w;
                   1008:   static int wlen;
                   1009:
                   1010:   if ( !p1 || !p2 )
                   1011:     *pr = 0;
                   1012:   else if ( OID(p1) != O_DP )
                   1013:     muldc(vl,p2,(Obj)p1,pr);
                   1014:   else if ( OID(p2) != O_DP )
                   1015:     muldc(vl,p1,(Obj)p2,pr);
                   1016:   else {
                   1017:     for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
                   1018:     if ( l > wlen ) {
                   1019:       if ( w ) GCFREE(w);
                   1020:       w = (MP *)MALLOC(l*sizeof(MP));
                   1021:       wlen = l;
                   1022:     }
                   1023:     for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
                   1024:       w[i] = m;
                   1025:     for ( s = 0, i = l-1; i >= 0; i-- ) {
                   1026:       weyl_muldm(vl,w[i],p2,&t); addd(vl,s,t,&u); s = u;
                   1027:     }
                   1028:     bzero(w,l*sizeof(MP));
                   1029:     *pr = s;
                   1030:   }
                   1031: }
                   1032:
                   1033: void actm(VL vl,int nv,MP m1,MP m2,DP *pr)
                   1034: {
                   1035:   DL d1,d2,d;
                   1036:   int n2,i,j,k;
                   1037:   Z jq,c,c1;
                   1038:   MP m;
                   1039:   Obj t;
                   1040:
                   1041:   d1 = m1->dl;
                   1042:   d2 = m2->dl;
                   1043:   for ( i = 0; i < nv; i++ )
                   1044:     if ( d1->d[i] > d2->d[i] ) {
                   1045:       *pr = 0; return;
                   1046:     }
                   1047:   NEWDL(d,nv);
                   1048:   c = ONE;
                   1049:   for ( i = 0; i < nv; i++ ) {
                   1050:     for ( j = d2->d[i], k = d1->d[i]; k > 0; k--, j-- ) {
1.2       noro     1051:       STOZ(j,jq); mulz(c,jq,&c1); c = c1;
1.1       noro     1052:     }
                   1053:     d->d[i] = d2->d[i]-d1->d[i];
                   1054:   }
                   1055:   arf_mul(vl,C(m1),C(m2),&t);
                   1056:   NEWMP(m);
                   1057:   arf_mul(vl,(Obj)c,t,&C(m));
                   1058:   m->dl = d;
                   1059:   MKDP(nv,m,*pr);
                   1060: }
                   1061:
                   1062: void weyl_actd(VL vl,DP p1,DP p2,DP *pr)
                   1063: {
                   1064:   int n;
                   1065:   MP m1,m2;
                   1066:   DP d,r,s;
                   1067:
                   1068:   if ( !p1 || !p2 ) *pr = 0;
                   1069:   else {
                   1070:     n = NV(p1);
                   1071:     r = 0;
                   1072:     for ( m1 = BDY(p1); m1; m1 = NEXT(m1) )
                   1073:       for ( m2 = BDY(p2); m2; m2 = NEXT(m2) ) {
                   1074:         actm(vl,n,m1,m2,&d);
                   1075:         addd(vl,r,d,&s); r = s;
                   1076:       }
                   1077:     *pr = r;
                   1078:   }
                   1079: }
                   1080:
                   1081: /* monomial * polynomial */
                   1082:
                   1083: void weyl_muldm(VL vl,MP m0,DP p,DP *pr)
                   1084: {
                   1085:   DP r,t,t1;
                   1086:   MP m;
                   1087:   DL d0;
                   1088:   int n,n2,l,i,j,tlen;
                   1089:   static MP *w,*psum;
                   1090:   static struct cdl *tab;
                   1091:   static int wlen;
                   1092:   static int rtlen;
                   1093:
                   1094:   if ( !p )
                   1095:     *pr = 0;
                   1096:   else {
                   1097:     for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
                   1098:     if ( l > wlen ) {
                   1099:       if ( w ) GCFREE(w);
                   1100:       w = (MP *)MALLOC(l*sizeof(MP));
                   1101:       wlen = l;
                   1102:     }
                   1103:     for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
                   1104:       w[i] = m;
                   1105:
                   1106:     n = NV(p); n2 = n>>1;
                   1107:     d0 = m0->dl;
                   1108:     for ( i = 0, tlen = 1; i < n2; i++ )
                   1109:       tlen *= d0->d[n2+i]+1;
                   1110:     if ( tlen > rtlen ) {
                   1111:       if ( tab ) GCFREE(tab);
                   1112:       if ( psum ) GCFREE(psum);
                   1113:       rtlen = tlen;
                   1114:       tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
                   1115:       psum = (MP *)MALLOC(rtlen*sizeof(MP));
                   1116:     }
                   1117:     bzero(psum,tlen*sizeof(MP));
                   1118:     for ( i = l-1; i >= 0; i-- ) {
                   1119:       bzero(tab,tlen*sizeof(struct cdl));
                   1120:       weyl_mulmm(vl,m0,w[i],n,tab,tlen);
                   1121:       for ( j = 0; j < tlen; j++ ) {
                   1122:         if ( tab[j].c ) {
                   1123:           NEWMP(m); m->dl = tab[j].d; C(m) = (Obj)tab[j].c; NEXT(m) = psum[j];
                   1124:           psum[j] = m;
                   1125:         }
                   1126:       }
                   1127:     }
                   1128:     for ( j = tlen-1, r = 0; j >= 0; j-- )
                   1129:       if ( psum[j] ) {
                   1130:         MKDP(n,psum[j],t); addd(vl,r,t,&t1); r = t1;
                   1131:       }
                   1132:     if ( r )
                   1133:       r->sugar = p->sugar + m0->dl->td;
                   1134:     *pr = r;
                   1135:   }
                   1136: }
                   1137:
                   1138: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
                   1139: /* rtab : array of length (e0+1)*(e1+1)*... */
                   1140:
                   1141: void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
                   1142: {
                   1143:   Obj c,c0,c1;
                   1144:   DL d,d0,d1,dt;
                   1145:   int i,j,a,b,k,l,n2,s,min,curlen;
                   1146:   struct cdl *p;
                   1147:   static Z *ctab;
                   1148:   static struct cdl *tab;
                   1149:   static int tablen;
                   1150:   static struct cdl *tmptab;
                   1151:   static int tmptablen;
                   1152:
                   1153:
                   1154:   if ( !m0 || !m1 ) {
                   1155:     rtab[0].c = 0;
                   1156:     rtab[0].d = 0;
                   1157:     return;
                   1158:   }
                   1159:   c0 = C(m0); c1 = C(m1);
                   1160:   arf_mul(vl,c0,c1,&c);
                   1161:   d0 = m0->dl; d1 = m1->dl;
                   1162:   n2 = n>>1;
                   1163:   curlen = 1;
                   1164:   NEWDL(d,n);
                   1165:   if ( n & 1 )
                   1166:     /* offset of h-degree */
                   1167:      d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
                   1168:   else
                   1169:     d->td = 0;
                   1170:   rtab[0].c = c;
                   1171:   rtab[0].d = d;
                   1172:
                   1173:   if ( rtablen > tmptablen ) {
                   1174:     if ( tmptab ) GCFREE(tmptab);
                   1175:     tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
                   1176:     tmptablen = rtablen;
                   1177:   }
                   1178:   for ( i = 0; i < n2; i++ ) {
                   1179:     a = d0->d[i]; b = d1->d[n2+i];
                   1180:     k = d0->d[n2+i]; l = d1->d[i];
                   1181:
                   1182:     /* degree of xi^a*(Di^k*xi^l)*Di^b */
                   1183:     a += l;
                   1184:     b += k;
                   1185:     s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
                   1186:
                   1187:     if ( !k || !l ) {
                   1188:       for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
                   1189:         if ( p->c ) {
                   1190:           dt = p->d;
                   1191:           dt->d[i] = a;
                   1192:           dt->d[n2+i] = b;
                   1193:           dt->td += s;
                   1194:         }
                   1195:       }
                   1196:       curlen *= k+1;
                   1197:       continue;
                   1198:     }
                   1199:     if ( k+1 > tablen ) {
                   1200:       if ( tab ) GCFREE(tab);
                   1201:       if ( ctab ) GCFREE(ctab);
                   1202:       tablen = k+1;
                   1203:       tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
                   1204:       ctab = (Z *)MALLOC(tablen*sizeof(Q));
                   1205:     }
                   1206:     /* compute xi^a*(Di^k*xi^l)*Di^b */
                   1207:     min = MIN(k,l);
                   1208:     mkwc(k,l,ctab);
                   1209:     bzero(tab,(k+1)*sizeof(struct cdl));
                   1210:     if ( n & 1 )
                   1211:       for ( j = 0; j <= min; j++ ) {
                   1212:         NEWDL(d,n);
                   1213:         d->d[i] = a-j; d->d[n2+i] = b-j;
                   1214:         d->td = s;
                   1215:         d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
                   1216:         tab[j].d = d;
                   1217:         tab[j].c = (Obj)ctab[j];
                   1218:       }
                   1219:     else
                   1220:       for ( j = 0; j <= min; j++ ) {
                   1221:         NEWDL(d,n);
                   1222:         d->d[i] = a-j; d->d[n2+i] = b-j;
                   1223:         d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
                   1224:         tab[j].d = d;
                   1225:         tab[j].c = (Obj)ctab[j];
                   1226:       }
                   1227:     bzero(ctab,(min+1)*sizeof(Q));
                   1228:     comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
                   1229:     curlen *= k+1;
                   1230:     bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
                   1231:   }
                   1232: }
                   1233:
                   1234: /* direct product of two cdl tables
                   1235:   rt[] = [
                   1236:     t[0]*t1[0],...,t[n-1]*t1[0],
                   1237:     t[0]*t1[1],...,t[n-1]*t1[1],
                   1238:     ...
                   1239:     t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
                   1240:   ]
                   1241: */
                   1242:
                   1243: void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
                   1244: {
                   1245:   int i,j;
                   1246:   struct cdl *p;
                   1247:   Obj c;
                   1248:   DL d;
                   1249:
                   1250:   bzero(rt,n*n1*sizeof(struct cdl));
                   1251:   for ( j = 0, p = rt; j < n1; j++ ) {
                   1252:     c = (Obj)t1[j].c;
                   1253:     d = t1[j].d;
                   1254:     if ( !c )
                   1255:       break;
                   1256:     for ( i = 0; i < n; i++, p++ ) {
                   1257:       if ( t[i].c ) {
                   1258:         arf_mul(vl,(Obj)t[i].c,c,(Obj *)&p->c);
                   1259:         adddl(nv,t[i].d,d,&p->d);
                   1260:       }
                   1261:     }
                   1262:   }
                   1263: }
                   1264:
                   1265: void muldc(VL vl,DP p,Obj c,DP *pr)
                   1266: {
                   1267:   MP m,mr=0,mr0;
                   1268:
                   1269:   if ( !p || !c )
                   1270:     *pr = 0;
                   1271:   else if ( NUM(c) && UNIQ((Q)c) )
                   1272:     *pr = p;
                   1273:   else if ( NUM(c) && MUNIQ((Q)c) )
                   1274:     chsgnd(p,pr);
                   1275:   else {
                   1276:     for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1277:       NEXTMP(mr0,mr);
                   1278:       if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                   1279:         mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                   1280:       else
                   1281:         arf_mul(vl,C(m),c,&C(mr));
                   1282:       mr->dl = m->dl;
                   1283:     }
                   1284:     NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                   1285:     if ( *pr )
                   1286:       (*pr)->sugar = p->sugar;
                   1287:   }
                   1288: }
                   1289:
                   1290: void divdc(VL vl,DP p,Obj c,DP *pr)
                   1291: {
                   1292:   Obj inv;
                   1293:
                   1294:   arf_div(vl,(Obj)ONE,c,&inv);
                   1295:   muld(vl,p,(DP)inv,pr);
                   1296: }
                   1297:
                   1298: void muldc_trunc(VL vl,DP p,Obj c,DL dl,DP *pr)
                   1299: {
                   1300:   MP m,mr=0,mr0;
                   1301:   DL mdl;
                   1302:   int i,n;
                   1303:
                   1304:   if ( !p || !c ) {
                   1305:     *pr = 0; return;
                   1306:   }
                   1307:   n = NV(p);
                   1308:   for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1309:     mdl = m->dl;
                   1310:     for ( i = 0; i < n; i++ )
                   1311:       if ( mdl->d[i] < dl->d[i] )
                   1312:         break;
                   1313:     if ( i < n )
                   1314:       break;
                   1315:     NEXTMP(mr0,mr);
                   1316:     if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                   1317:       mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                   1318:     else
                   1319:       arf_mul(vl,C(m),c,&C(mr));
                   1320:     mr->dl = m->dl;
                   1321:   }
                   1322:   NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                   1323:   if ( *pr )
                   1324:     (*pr)->sugar = p->sugar;
                   1325: }
                   1326:
                   1327: void divsdc(VL vl,DP p,P c,DP *pr)
                   1328: {
                   1329:   MP m,mr=0,mr0;
                   1330:
                   1331:   if ( !c )
                   1332:     error("disvsdc : division by 0");
                   1333:   else if ( !p )
                   1334:     *pr = 0;
                   1335:   else if ( OID(c) > O_P )
                   1336:     error("divsdc : invalid argument");
                   1337:   else {
                   1338:     for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1339:       NEXTMP(mr0,mr); divsp(vl,(P)C(m),c,(P *)&C(mr)); mr->dl = m->dl;
                   1340:     }
                   1341:     NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                   1342:     if ( *pr )
                   1343:       (*pr)->sugar = p->sugar;
                   1344:   }
                   1345: }
                   1346:
                   1347: void adddl(int n,DL d1,DL d2,DL *dr)
                   1348: {
                   1349:   DL dt;
                   1350:   int i;
                   1351:
                   1352:   *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
                   1353:   dt->td = d1->td + d2->td;
                   1354:   for ( i = 0; i < n; i++ )
                   1355:     dt->d[i] = d1->d[i]+d2->d[i];
                   1356: }
                   1357:
                   1358: /* d1 += d2 */
                   1359:
                   1360: void adddl_destructive(int n,DL d1,DL d2)
                   1361: {
                   1362:   int i;
                   1363:
                   1364:   d1->td += d2->td;
                   1365:   for ( i = 0; i < n; i++ )
                   1366:     d1->d[i] += d2->d[i];
                   1367: }
                   1368:
                   1369: int compd(VL vl,DP p1,DP p2)
                   1370: {
                   1371:   int n,t;
                   1372:   MP m1,m2;
                   1373:
                   1374:   if ( !p1 )
                   1375:     return p2 ? -1 : 0;
                   1376:   else if ( !p2 )
                   1377:     return 1;
                   1378:   else if ( NV(p1) != NV(p2) ) {
                   1379:     error("compd : size mismatch");
                   1380:     return 0; /* XXX */
                   1381:   } else {
                   1382:     for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
                   1383:       m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
                   1384:       if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
                   1385:         (t = arf_comp(vl,C(m1),C(m2)) ) )
                   1386:         return t;
                   1387:     if ( m1 )
                   1388:       return 1;
                   1389:     else if ( m2 )
                   1390:       return -1;
                   1391:     else
                   1392:       return 0;
                   1393:   }
                   1394: }
                   1395:
                   1396: int cmpdl_lex(int n,DL d1,DL d2)
                   1397: {
                   1398:   int i;
                   1399:
                   1400:   for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
                   1401:   return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
                   1402: }
                   1403:
                   1404: int cmpdl_revlex(int n,DL d1,DL d2)
                   1405: {
                   1406:   int i;
                   1407:
                   1408:   for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
                   1409:   return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                   1410: }
                   1411:
                   1412: int cmpdl_gradlex(int n,DL d1,DL d2)
                   1413: {
                   1414:   if ( d1->td > d2->td )
                   1415:     return 1;
                   1416:   else if ( d1->td < d2->td )
                   1417:     return -1;
                   1418:   else
                   1419:     return cmpdl_lex(n,d1,d2);
                   1420: }
                   1421:
                   1422: int cmpdl_revgradlex(int n,DL d1,DL d2)
                   1423: {
                   1424:   register int i,c;
                   1425:   register int *p1,*p2;
                   1426:
                   1427:   if ( d1->td > d2->td )
                   1428:     return 1;
                   1429:   else if ( d1->td < d2->td )
                   1430:     return -1;
                   1431:   else {
                   1432:     i = n-1;
                   1433:     p1 = d1->d+n-1;
                   1434:     p2 = d2->d+n-1;
                   1435:     while ( i >= 7 ) {
                   1436:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1437:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1438:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1439:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1440:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1441:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1442:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1443:       c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1444:       i -= 8;
                   1445:     }
                   1446:     switch ( i ) {
                   1447:       case 6:
                   1448:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1449:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1450:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1451:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1452:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1453:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1454:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1455:         return 0;
                   1456:       case 5:
                   1457:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1458:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1459:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1460:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1461:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1462:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1463:         return 0;
                   1464:       case 4:
                   1465:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1466:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1467:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1468:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1469:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1470:         return 0;
                   1471:       case 3:
                   1472:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1473:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1474:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1475:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1476:         return 0;
                   1477:       case 2:
                   1478:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1479:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1480:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1481:         return 0;
                   1482:       case 1:
                   1483:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1484:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1485:         return 0;
                   1486:       case 0:
                   1487:         c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1488:         return 0;
                   1489:       default:
                   1490:         return 0;
                   1491:     }
                   1492: LAST:
                   1493:     if ( c > 0 ) return -1;
                   1494:     else return 1;
                   1495:   }
                   1496: }
                   1497:
                   1498: int cmpdl_blex(int n,DL d1,DL d2)
                   1499: {
                   1500:   int c;
                   1501:
                   1502:   if ( (c = cmpdl_lex(n-1,d1,d2)) )
                   1503:     return c;
                   1504:   else {
                   1505:     c = d1->d[n-1] - d2->d[n-1];
                   1506:     return c > 0 ? 1 : c < 0 ? -1 : 0;
                   1507:   }
                   1508: }
                   1509:
                   1510: int cmpdl_bgradlex(int n,DL d1,DL d2)
                   1511: {
                   1512:   int e1,e2,c;
                   1513:
                   1514:   e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
                   1515:   if ( e1 > e2 )
                   1516:     return 1;
                   1517:   else if ( e1 < e2 )
                   1518:     return -1;
                   1519:   else {
                   1520:     c = cmpdl_lex(n-1,d1,d2);
                   1521:     if ( c )
                   1522:       return c;
                   1523:     else
                   1524:       return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
                   1525:   }
                   1526: }
                   1527:
                   1528: int cmpdl_brevgradlex(int n,DL d1,DL d2)
                   1529: {
                   1530:   int e1,e2,c;
                   1531:
                   1532:   e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
                   1533:   if ( e1 > e2 )
                   1534:     return 1;
                   1535:   else if ( e1 < e2 )
                   1536:     return -1;
                   1537:   else {
                   1538:     c = cmpdl_revlex(n-1,d1,d2);
                   1539:     if ( c )
                   1540:       return c;
                   1541:     else
                   1542:       return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
                   1543:   }
                   1544: }
                   1545:
                   1546: int cmpdl_brevrev(int n,DL d1,DL d2)
                   1547: {
                   1548:   int e1,e2,f1,f2,c,i;
                   1549:
                   1550:   for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                   1551:     e1 += d1->d[i]; e2 += d2->d[i];
                   1552:   }
                   1553:   f1 = d1->td - e1; f2 = d2->td - e2;
                   1554:   if ( e1 > e2 )
                   1555:     return 1;
                   1556:   else if ( e1 < e2 )
                   1557:     return -1;
                   1558:   else {
                   1559:     c = cmpdl_revlex(dp_nelim,d1,d2);
                   1560:     if ( c )
                   1561:       return c;
                   1562:     else if ( f1 > f2 )
                   1563:       return 1;
                   1564:     else if ( f1 < f2 )
                   1565:       return -1;
                   1566:     else {
                   1567:       for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                   1568:       return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                   1569:     }
                   1570:   }
                   1571: }
                   1572:
                   1573: int cmpdl_bgradrev(int n,DL d1,DL d2)
                   1574: {
                   1575:   int e1,e2,f1,f2,c,i;
                   1576:
                   1577:   for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                   1578:     e1 += d1->d[i]; e2 += d2->d[i];
                   1579:   }
                   1580:   f1 = d1->td - e1; f2 = d2->td - e2;
                   1581:   if ( e1 > e2 )
                   1582:     return 1;
                   1583:   else if ( e1 < e2 )
                   1584:     return -1;
                   1585:   else {
                   1586:     c = cmpdl_lex(dp_nelim,d1,d2);
                   1587:     if ( c )
                   1588:       return c;
                   1589:     else if ( f1 > f2 )
                   1590:       return 1;
                   1591:     else if ( f1 < f2 )
                   1592:       return -1;
                   1593:     else {
                   1594:       for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                   1595:       return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                   1596:     }
                   1597:   }
                   1598: }
                   1599:
                   1600: int cmpdl_blexrev(int n,DL d1,DL d2)
                   1601: {
                   1602:   int e1,e2,f1,f2,c,i;
                   1603:
                   1604:   for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                   1605:     e1 += d1->d[i]; e2 += d2->d[i];
                   1606:   }
                   1607:   f1 = d1->td - e1; f2 = d2->td - e2;
                   1608:   c = cmpdl_lex(dp_nelim,d1,d2);
                   1609:   if ( c )
                   1610:     return c;
                   1611:   else if ( f1 > f2 )
                   1612:     return 1;
                   1613:   else if ( f1 < f2 )
                   1614:     return -1;
                   1615:   else {
                   1616:     for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                   1617:     return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                   1618:   }
                   1619: }
                   1620:
                   1621: int cmpdl_elim(int n,DL d1,DL d2)
                   1622: {
                   1623:   int e1,e2,i;
                   1624:
                   1625:   for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                   1626:     e1 += d1->d[i]; e2 += d2->d[i];
                   1627:   }
                   1628:   if ( e1 > e2 )
                   1629:     return 1;
                   1630:   else if ( e1 < e2 )
                   1631:     return -1;
                   1632:   else
                   1633:     return cmpdl_revgradlex(n,d1,d2);
                   1634: }
                   1635:
                   1636: int cmpdl_weyl_elim(int n,DL d1,DL d2)
                   1637: {
                   1638:   int e1,e2,i;
                   1639:
                   1640:   for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
                   1641:     e1 += d1->d[n-i]; e2 += d2->d[n-i];
                   1642:   }
                   1643:   if ( e1 > e2 )
                   1644:     return 1;
                   1645:   else if ( e1 < e2 )
                   1646:     return -1;
                   1647:   else if ( d1->td > d2->td )
                   1648:     return 1;
                   1649:   else if ( d1->td < d2->td )
                   1650:     return -1;
                   1651:   else return -cmpdl_revlex(n,d1,d2);
                   1652: }
                   1653:
                   1654: /*
                   1655:   a special ordering
                   1656:   1. total order
                   1657:   2. (-w,w) for the first 2*m variables
                   1658:   3. DRL for the first 2*m variables
                   1659: */
                   1660:
                   1661: extern int *current_weyl_weight_vector;
                   1662:
                   1663: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
                   1664: {
                   1665:   int e1,e2,m,i;
                   1666:   int *p1,*p2;
                   1667:
                   1668:   if ( d1->td > d2->td )
                   1669:     return 1;
                   1670:   else if ( d1->td < d2->td )
                   1671:     return -1;
                   1672:
                   1673:   m = n>>1;
                   1674:   for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
                   1675:     e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
                   1676:     e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
                   1677:   }
                   1678:   if ( e1 > e2 )
                   1679:     return 1;
                   1680:   else if ( e1 < e2 )
                   1681:     return -1;
                   1682:
                   1683:   e1 = d1->td - d1->d[n-1];
                   1684:   e2 = d2->td - d2->d[n-1];
                   1685:   if ( e1 > e2 )
                   1686:     return 1;
                   1687:   else if ( e1 < e2 )
                   1688:     return -1;
                   1689:
                   1690:   for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
                   1691:     i >= 0 && *p1 == *p2; i--, p1--, p2-- );
                   1692:   return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
                   1693: }
                   1694:
                   1695: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
                   1696: {
                   1697:   int i,t,m;
                   1698:   int *p1,*p2;
                   1699:
                   1700:   if ( d1->td > d2->td )
                   1701:     return 1;
                   1702:   else if ( d1->td < d2->td )
                   1703:     return -1;
                   1704:   else {
                   1705:     m = n>>1;
                   1706:     for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
                   1707:       if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1;
                   1708:       if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1;
                   1709:     }
                   1710:     return 0;
                   1711:   }
                   1712: }
                   1713:
                   1714: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
                   1715: {
                   1716:   int e1,e2,m,i,t;
                   1717:   int *p1,*p2;
                   1718:
                   1719:   if ( d1->td > d2->td )
                   1720:     return 1;
                   1721:   else if ( d1->td < d2->td )
                   1722:     return -1;
                   1723:
                   1724:   m = n>>1;
                   1725:   for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
                   1726:     e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
                   1727:     e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
                   1728:   }
                   1729:   if ( e1 > e2 )
                   1730:     return 1;
                   1731:   else if ( e1 < e2 )
                   1732:     return -1;
                   1733:
                   1734:   e1 = d1->td - d1->d[n-1];
                   1735:   e2 = d2->td - d2->d[n-1];
                   1736:   if ( e1 > e2 )
                   1737:     return 1;
                   1738:   else if ( e1 < e2 )
                   1739:     return -1;
                   1740:
                   1741:   for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
                   1742:     if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1;
                   1743:     if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1;
                   1744:   }
                   1745:   return 0;
                   1746: }
                   1747:
                   1748: int cmpdl_order_pair(int n,DL d1,DL d2)
                   1749: {
                   1750:   int e1,e2,i,j,l;
                   1751:   int *t1,*t2;
                   1752:   int len,head;
                   1753:   struct order_pair *pair;
                   1754:
                   1755:   len = dp_current_spec->ord.block.length;
                   1756:   if ( n != dp_current_spec->nv )
                   1757:     error("cmpdl_order_pair : incompatible order specification");
                   1758:   pair = dp_current_spec->ord.block.order_pair;
                   1759:
                   1760:   head = 0;
                   1761:   for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
                   1762:     l = pair[i].length;
                   1763:     switch ( pair[i].order ) {
                   1764:       case 0:
                   1765:         for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
                   1766:           e1 += MUL_WEIGHT(t1[j],head+j);
                   1767:           e2 += MUL_WEIGHT(t2[j],head+j);
                   1768:         }
                   1769:         if ( e1 > e2 )
                   1770:           return 1;
                   1771:         else if ( e1 < e2 )
                   1772:           return -1;
                   1773:         else {
                   1774:           for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
                   1775:           if ( j >= 0 )
                   1776:             return t1[j] < t2[j] ? 1 : -1;
                   1777:         }
                   1778:         break;
                   1779:       case 1:
                   1780:         for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
                   1781:           e1 += MUL_WEIGHT(t1[j],head+j);
                   1782:           e2 += MUL_WEIGHT(t2[j],head+j);
                   1783:         }
                   1784:         if ( e1 > e2 )
                   1785:           return 1;
                   1786:         else if ( e1 < e2 )
                   1787:           return -1;
                   1788:         else {
                   1789:           for ( j = 0; j < l && t1[j] == t2[j]; j++ );
                   1790:           if ( j < l )
                   1791:             return t1[j] > t2[j] ? 1 : -1;
                   1792:         }
                   1793:         break;
                   1794:       case 2:
                   1795:         for ( j = 0; j < l && t1[j] == t2[j]; j++ );
                   1796:         if ( j < l )
                   1797:           return t1[j] > t2[j] ? 1 : -1;
                   1798:         break;
                   1799:       default:
                   1800:         error("cmpdl_order_pair : invalid order"); break;
                   1801:     }
                   1802:     t1 += l; t2 += l; head += l;
                   1803:   }
                   1804:   return 0;
                   1805: }
                   1806:
                   1807: int cmpdl_composite(int nv,DL d1,DL d2)
                   1808: {
                   1809:   int n,i,j,k,start,s,len;
                   1810:   int *dw;
                   1811:   struct sparse_weight *sw;
                   1812:   struct weight_or_block *worb;
                   1813:   int *w,*t1,*t2;
                   1814:
                   1815:   n = dp_current_spec->ord.composite.length;
                   1816:   worb = dp_current_spec->ord.composite.w_or_b;
                   1817:   w = dp_dl_work;
                   1818:   for ( i = 0, t1 = d1->d, t2 = d2->d; i < nv; i++ )
                   1819:     w[i] = t1[i]-t2[i];
                   1820:   for ( i = 0; i < n; i++, worb++ ) {
                   1821:     len = worb->length;
                   1822:     switch ( worb->type ) {
                   1823:       case IS_DENSE_WEIGHT:
                   1824:         dw = worb->body.dense_weight;
                   1825:         for ( j = 0, s = 0; j < len; j++ )
                   1826:           s += dw[j]*w[j];
                   1827:         if ( s > 0 ) return 1;
                   1828:         else if ( s < 0 ) return -1;
                   1829:         break;
                   1830:       case IS_SPARSE_WEIGHT:
                   1831:         sw = worb->body.sparse_weight;
                   1832:         for ( j = 0, s = 0; j < len; j++ )
                   1833:           s += sw[j].value*w[sw[j].pos];
                   1834:         if ( s > 0 ) return 1;
                   1835:         else if ( s < 0 ) return -1;
                   1836:         break;
                   1837:       case IS_BLOCK:
                   1838:         start = worb->body.block.start;
                   1839:         switch ( worb->body.block.order ) {
                   1840:           case 0:
                   1841:             for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
                   1842:               s += MUL_WEIGHT(w[k],k);
                   1843:             }
                   1844:             if ( s > 0 ) return 1;
                   1845:             else if ( s < 0 ) return -1;
                   1846:             else {
                   1847:               for ( j = k-1; j >= start && w[j] == 0; j-- );
                   1848:               if ( j >= start )
                   1849:                 return w[j] < 0 ? 1 : -1;
                   1850:             }
                   1851:             break;
                   1852:           case 1:
                   1853:             for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
                   1854:               s += MUL_WEIGHT(w[k],k);
                   1855:             }
                   1856:             if ( s > 0 ) return 1;
                   1857:             else if ( s < 0 ) return -1;
                   1858:             else {
                   1859:               for ( j = 0, k = start;  j < len && w[j] == 0; j++, k++ );
                   1860:               if ( j < len )
                   1861:                 return w[j] > 0 ? 1 : -1;
                   1862:             }
                   1863:             break;
                   1864:           case 2:
                   1865:             for ( j = 0, k = start;  j < len && w[j] == 0; j++, k++ );
                   1866:             if ( j < len )
                   1867:               return w[j] > 0 ? 1 : -1;
                   1868:             break;
                   1869:         }
                   1870:         break;
                   1871:     }
                   1872:   }
                   1873:   return 0;
                   1874: }
                   1875:
                   1876: int cmpdl_matrix(int n,DL d1,DL d2)
                   1877: {
                   1878:   int *v,*w,*t1,*t2;
                   1879:   int s,i,j,len;
                   1880:   int **matrix;
                   1881:
                   1882:   for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
                   1883:     w[i] = t1[i]-t2[i];
                   1884:   len = dp_current_spec->ord.matrix.row;
                   1885:   matrix = dp_current_spec->ord.matrix.matrix;
                   1886:   for ( j = 0; j < len; j++ ) {
                   1887:     v = matrix[j];
                   1888:     for ( i = 0, s = 0; i < n; i++ )
                   1889:       s += v[i]*w[i];
                   1890:     if ( s > 0 )
                   1891:       return 1;
                   1892:     else if ( s < 0 )
                   1893:       return -1;
                   1894:   }
                   1895:   return 0;
                   1896: }
                   1897:
                   1898: int cmpdl_top_weight(int n,DL d1,DL d2)
                   1899: {
                   1900:   int *w;
                   1901:   Z **mat;
                   1902:   Z *a;
                   1903:   mpz_t sum;
                   1904:   int len,i,sgn,tsgn,row,k;
                   1905:   int *t1,*t2;
                   1906:
                   1907:   w = (int *)ALLOCA(n*sizeof(int));
                   1908:   len = current_top_weight_len+3;
                   1909:   t1 = d1->d; t2 = d2->d;
                   1910:   for ( i = 0; i < n; i++ ) w[i] = t1[i]-t2[i];
                   1911:   mpz_init_set_ui(sum,0);
                   1912:   if ( OID(current_top_weight) == O_VECT ) {
                   1913:       mat = (Z **)&BDY((VECT)current_top_weight);
                   1914:     row = 1;
                   1915:   } else {
                   1916:       mat = (Z **)BDY((MAT)current_top_weight);
                   1917:     row = ((MAT)current_top_weight)->row;
                   1918:   }
                   1919:   for ( k = 0; k < row; k++ ) {
                   1920:     a = mat[k];
                   1921:     for ( i = 0; i < n; i++ ) {
                   1922:       if ( !a[i] || !w[i] ) continue;
                   1923:       if ( w[i] > 0 )
                   1924:         mpz_addmul_ui(sum,BDY(a[i]),(unsigned int)w[i]);
                   1925:       else
                   1926:         mpz_submul_ui(sum,BDY(a[i]),(unsigned int)(-w[i]));
                   1927:     }
                   1928:     sgn = mpz_sgn(sum);
                   1929:     if ( sgn > 0 ) return 1;
                   1930:     else if ( sgn < 0 ) return -1;
                   1931:   }
                   1932:   return (*cmpdl_tie_breaker)(n,d1,d2);
                   1933: }
                   1934:
                   1935: GeoBucket create_bucket()
                   1936: {
                   1937:   GeoBucket g;
                   1938:
                   1939:   g = CALLOC(1,sizeof(struct oGeoBucket));
                   1940:   g->m = 32;
                   1941:   return g;
                   1942: }
                   1943:
                   1944: int length(NODE d);
                   1945:
                   1946: void add_bucket(GeoBucket g,NODE d,int nv)
                   1947: {
                   1948:   int l,k,m;
                   1949:
                   1950:   l = length(d);
                   1951:   for ( k = 0, m = 1; l > m; k++, m <<= 1 );
                   1952:   /* 2^(k-1) < l <= 2^k */
                   1953:   d = symb_merge(g->body[k],d,nv);
                   1954:   for ( ; length(d) > (1<<(k)); k++ ) {
                   1955:     g->body[k] = 0;
                   1956:     d = symb_merge(g->body[k+1],d,nv);
                   1957:   }
                   1958:   g->body[k] = d;
                   1959:   g->m = MAX(g->m,k);
                   1960: }
                   1961:
                   1962: DL remove_head_bucket(GeoBucket g,int nv)
                   1963: {
                   1964:   int j,i,c,m;
                   1965:   DL d;
                   1966:
                   1967:   j = -1;
                   1968:   m = g->m;
                   1969:   for ( i = 0; i <= m; i++ ) {
                   1970:     if ( !g->body[i] )
                   1971:       continue;
                   1972:     if ( j < 0 ) j = i;
                   1973:     else {
                   1974:       c = (*cmpdl)(nv,g->body[i]->body,g->body[j]->body);
                   1975:       if ( c > 0 )
                   1976:         j = i;
                   1977:       else if ( c == 0 )
                   1978:         g->body[i] = NEXT(g->body[i]);
                   1979:     }
                   1980:   }
                   1981:   if ( j < 0 )
                   1982:     return 0;
                   1983:   else {
                   1984:     d = g->body[j]->body;
                   1985:     g->body[j] = NEXT(g->body[j]);
                   1986:     return d;
                   1987:   }
                   1988: }
                   1989:
                   1990: /*  DPV functions */
                   1991:
                   1992: void adddv(VL vl,DPV p1,DPV p2,DPV *pr)
                   1993: {
                   1994:   int i,len;
                   1995:   DP *e;
                   1996:
                   1997:   if ( !p1 || !p2 )
                   1998:     error("adddv : invalid argument");
                   1999:   else if ( p1->len != p2->len )
                   2000:     error("adddv : size mismatch");
                   2001:   else {
                   2002:     len = p1->len;
                   2003:     e = (DP *)MALLOC(p1->len*sizeof(DP));
                   2004:     for ( i = 0; i < len; i++ )
                   2005:       addd(vl,p1->body[i],p2->body[i],&e[i]);
                   2006:     MKDPV(len,e,*pr);
                   2007:     (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                   2008:   }
                   2009: }
                   2010:
                   2011: void subdv(VL vl,DPV p1,DPV p2,DPV *pr)
                   2012: {
                   2013:   int i,len;
                   2014:   DP *e;
                   2015:
                   2016:   if ( !p1 || !p2 )
                   2017:     error("subdv : invalid argument");
                   2018:   else if ( p1->len != p2->len )
                   2019:     error("subdv : size mismatch");
                   2020:   else {
                   2021:     len = p1->len;
                   2022:     e = (DP *)MALLOC(p1->len*sizeof(DP));
                   2023:     for ( i = 0; i < len; i++ )
                   2024:       subd(vl,p1->body[i],p2->body[i],&e[i]);
                   2025:     MKDPV(len,e,*pr);
                   2026:     (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                   2027:   }
                   2028: }
                   2029:
                   2030: void chsgndv(DPV p1,DPV *pr)
                   2031: {
                   2032:   int i,len;
                   2033:   DP *e;
                   2034:
                   2035:   if ( !p1 )
                   2036:     error("subdv : invalid argument");
                   2037:   else {
                   2038:     len = p1->len;
                   2039:     e = (DP *)MALLOC(p1->len*sizeof(DP));
                   2040:     for ( i = 0; i < len; i++ )
                   2041:       chsgnd(p1->body[i],&e[i]);
                   2042:     MKDPV(len,e,*pr);
                   2043:     (*pr)->sugar = p1->sugar;
                   2044:   }
                   2045: }
                   2046:
                   2047: void muldv(VL vl,DP p1,DPV p2,DPV *pr)
                   2048: {
                   2049:   int i,len;
                   2050:   DP *e;
                   2051:
                   2052:   len = p2->len;
                   2053:   e = (DP *)MALLOC(p2->len*sizeof(DP));
                   2054:   if ( !p1 ) {
                   2055:     MKDPV(len,e,*pr);
                   2056:     (*pr)->sugar = 0;
                   2057:   } else {
                   2058:     for ( i = 0; i < len; i++ )
                   2059:       muld(vl,p1,p2->body[i],&e[i]);
                   2060:     MKDPV(len,e,*pr);
                   2061:     (*pr)->sugar = p1->sugar + p2->sugar;
                   2062:   }
                   2063: }
                   2064:
                   2065: int compdv(VL vl,DPV p1,DPV p2)
                   2066: {
                   2067:   int i,t,len;
                   2068:
                   2069:   if ( p1->len != p2->len ) {
                   2070:     error("compdv : size mismatch");
                   2071:     return 0; /* XXX */
                   2072:   } else {
                   2073:     len = p1->len;
                   2074:     for ( i = 0; i < len; i++ )
                   2075:       if ( (t = compd(vl,p1->body[i],p2->body[i])) )
                   2076:         return t;
                   2077:     return 0;
                   2078:   }
                   2079: }
                   2080:
                   2081: int ni_next(int *a,int n)
                   2082: {
                   2083:   int i,j,k,kj;
                   2084:
                   2085:   /* find the first nonzero a[j] */
                   2086:   for ( j = 0; j < n && a[j] == 0; j++ );
                   2087:   /* find the first zero a[k] after a[j] */
                   2088:   for ( k = j; k < n && a[k] == 1; k++ );
                   2089:   if ( k == n ) return 0;
                   2090:   /* a[0] = 0, ... , a[j-1] = 0, a[j] = 1, ..., a[k-1] = 1, a[k] = 0 */
                   2091:   /* a[0] = 1,..., a[k-j-2] = 1, a[k-j-1] = 0, ..., a[k-1] = 0, a[k] = 1 */
                   2092:   kj = k-j-1;
                   2093:   for ( i = 0; i < kj; i++ ) a[i] = 1;
                   2094:   for ( ; i < k; i++ ) a[i] = 0;
                   2095:   a[k] = 1;
                   2096:   return 1;
                   2097: }
                   2098:
                   2099: int comp_nbm(NBM a,NBM b)
                   2100: {
                   2101:   int d,i,ai,bi;
                   2102:   unsigned int *ab,*bb;
                   2103:
                   2104:   if ( a->d > b->d ) return 1;
                   2105:   else if ( a->d < b->d ) return -1;
                   2106:   else {
                   2107:     d = a->d; ab = a->b; bb = b->b;
                   2108: #if 0
                   2109:     w = (d+31)/32;
                   2110:     for ( i = 0; i < w; i++ )
                   2111:       if ( ab[i] > bb[i] ) return 1;
                   2112:       else if ( ab[i] < bb[i] ) return -1;
                   2113: #else
                   2114:     for ( i = 0; i < d; i++ ) {
                   2115:       ai = NBM_GET(ab,i);
                   2116:       bi = NBM_GET(bb,i);
                   2117:       if ( ai > bi ) return 1;
                   2118:       else if ( ai < bi ) return -1;
                   2119:     }
                   2120: #endif
                   2121:     return 0;
                   2122:   }
                   2123: }
                   2124:
                   2125: NBM mul_nbm(NBM a,NBM b)
                   2126: {
                   2127:   int ad,bd,d,i,j;
                   2128:   unsigned int *ab,*bb,*mb;
                   2129:   NBM m;
                   2130:
                   2131:   ad = a->d; bd = b->d; ab = a->b; bb = b->b;
                   2132:   d = ad + bd;
                   2133:   NEWNBM(m); NEWNBMBDY(m,d);
                   2134:   m->d = d; mulp(CO,a->c,b->c,&m->c); mb = m->b;
                   2135:   j = 0;
                   2136:   for ( i = 0; i < ad; i++, j++ )
                   2137:     if ( NBM_GET(ab,i) ) NBM_SET(mb,j);
                   2138:     else NBM_CLR(mb,j);
                   2139:   for ( i = 0; i < bd; i++, j++ )
                   2140:     if ( NBM_GET(bb,i) ) NBM_SET(mb,j);
                   2141:     else NBM_CLR(mb,j);
                   2142:   return m;
                   2143: }
                   2144:
                   2145: NBP nbmtonbp(NBM m)
                   2146: {
                   2147:   NODE n;
                   2148:   NBP u;
                   2149:
                   2150:   MKNODE(n,m,0);
                   2151:   MKNBP(u,n);
                   2152:   return u;
                   2153: }
                   2154:
                   2155: /* a=c*x*rest -> a0= x*rest, ah=x, ar=rest */
                   2156:
                   2157: P separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
                   2158: {
                   2159:   int i,d1;
                   2160:   NBM t;
                   2161:
                   2162:   if ( !a->d ) error("separate_nbm : invalid argument");
                   2163:
                   2164:   if ( a0 ) {
                   2165:     NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
                   2166:     *a0 = nbmtonbp(t);
                   2167:   }
                   2168:
                   2169:   if ( ah ) {
                   2170:     NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
                   2171:     if ( NBM_GET(a->b,0) ) NBM_SET(t->b,0);
                   2172:     else NBM_CLR(t->b,0);
                   2173:     *ah = nbmtonbp(t);
                   2174:   }
                   2175:
                   2176:   if ( ar ) {
                   2177:     d1 = a->d-1;
                   2178:     NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
                   2179:     for ( i = 0; i < d1; i++ ) {
                   2180:       if ( NBM_GET(a->b,i+1) ) NBM_SET(t->b,i);
                   2181:       else NBM_CLR(t->b,i);
                   2182:     }
                   2183:     *ar = nbmtonbp(t);
                   2184:   }
                   2185:
                   2186:   return a->c;
                   2187: }
                   2188:
                   2189: /* a=c*rest*x -> a0= rest*x, ar=rest, at=x */
                   2190:
                   2191: P separate_tail_nbm(NBM a,NBP *a0,NBP *ar,NBP *at)
                   2192: {
                   2193:   int i,d,d1;
                   2194:   NBM t;
                   2195:
                   2196:   if ( !(d=a->d) ) error("separate_tail_nbm : invalid argument");
                   2197:
                   2198:   if ( a0 ) {
                   2199:     NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
                   2200:     *a0 = nbmtonbp(t);
                   2201:   }
                   2202:
                   2203:   d1 = a->d-1;
                   2204:   if ( at ) {
                   2205:     NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
                   2206:     if ( NBM_GET(a->b,d1) ) NBM_SET(t->b,0);
                   2207:     else NBM_CLR(t->b,0);
                   2208:     *at = nbmtonbp(t);
                   2209:   }
                   2210:
                   2211:   if ( ar ) {
                   2212:     NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
                   2213:     for ( i = 0; i < d1; i++ ) {
                   2214:       if ( NBM_GET(a->b,i) ) NBM_SET(t->b,i);
                   2215:       else NBM_CLR(t->b,i);
                   2216:     }
                   2217:     *ar = nbmtonbp(t);
                   2218:   }
                   2219:
                   2220:   return a->c;
                   2221: }
                   2222:
                   2223: NBP make_xky(int k)
                   2224: {
                   2225:   int k1,i;
                   2226:   NBM t;
                   2227:
                   2228:   NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
                   2229:   k1 = k-1;
                   2230:   for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
                   2231:   NBM_CLR(t->b,i);
                   2232:   return nbmtonbp(t);
                   2233: }
                   2234:
                   2235: /* a=c*x^(k-1)*y*rest -> a0= x^(k-1)*y*rest, ah=x^(k-1)*y, ar=rest */
                   2236:
                   2237: P separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
                   2238: {
                   2239:   int i,d1,k,k1;
                   2240:   NBM t;
                   2241:
                   2242:   if ( !a->d )
                   2243:     error("separate_nbm : invalid argument");
                   2244:   for ( i = 0; i < a->d && NBM_GET(a->b,i); i++ );
                   2245:   if ( i == a->d )
                   2246:     error("separate_nbm : invalid argument");
                   2247:   k1 = i;
                   2248:   k = i+1;
                   2249:
                   2250:   if ( a0 ) {
                   2251:     NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
                   2252:     *a0 = nbmtonbp(t);
                   2253:   }
                   2254:
                   2255:   if ( ah ) {
                   2256:     NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
                   2257:     for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
                   2258:     NBM_CLR(t->b,i);
                   2259:     *ah = nbmtonbp(t);
                   2260:   }
                   2261:
                   2262:   if ( ar ) {
                   2263:     d1 = a->d-k;
                   2264:     NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
                   2265:     for ( i = 0; i < d1; i++ ) {
                   2266:       if ( NBM_GET(a->b,i+k) ) NBM_SET(t->b,i);
                   2267:       else NBM_CLR(t->b,i);
                   2268:     }
                   2269:     *ar = nbmtonbp(t);
                   2270:   }
                   2271:
                   2272:   return a->c;
                   2273: }
                   2274:
                   2275: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
                   2276: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
                   2277: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp);
                   2278: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
                   2279:
                   2280: NBP shuffle_mul_nbm(NBM a,NBM b)
                   2281: {
                   2282:   NBP u,a0,ah,ar,b0,bh,br,a1,b1,t;
                   2283:   P ac,bc,c;
                   2284:
                   2285:   if ( !a->d || !b->d )
                   2286:     u = nbmtonbp(mul_nbm(a,b));
                   2287:   else {
                   2288:     ac = separate_nbm(a,&a0,&ah,&ar);
                   2289:     bc = separate_nbm(b,&b0,&bh,&br);
                   2290:     mulp(CO,ac,bc,&c);
                   2291:     shuffle_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
                   2292:     shuffle_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
                   2293:     addnbp(CO,a1,b1,&t); mulnbp(CO,(NBP)c,t,&u);
                   2294:   }
                   2295:   return u;
                   2296: }
                   2297:
                   2298: NBP harmonic_mul_nbm(NBM a,NBM b)
                   2299: {
                   2300:   NBP u,a0,ah,ar,b0,bh,br,a1,b1,t,s,abk,ab1;
                   2301:   P ac,bc,c;
                   2302:
                   2303:   if ( !a->d || !b->d )
                   2304:     u = nbmtonbp(mul_nbm(a,b));
                   2305:   else {
                   2306:     mulp(CO,a->c,b->c,&c);
                   2307:     ac = separate_xky_nbm(a,&a0,&ah,&ar);
                   2308:     bc = separate_xky_nbm(b,&b0,&bh,&br);
                   2309:     mulp(CO,ac,bc,&c);
                   2310:     harmonic_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
                   2311:     harmonic_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
                   2312:     abk = make_xky(((NBM)BDY(BDY(ah)))->d+((NBM)BDY(BDY(bh)))->d);
                   2313:     harmonic_mulnbp(CO,ar,br,&t); mulnbp(CO,abk,t,&ab1);
                   2314:     addnbp(CO,a1,b1,&t); addnbp(CO,t,ab1,&s); mulnbp(CO,(NBP)c,s,&u);
                   2315:   }
                   2316:   return u;
                   2317:
                   2318: }
                   2319:
                   2320: void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2321: {
                   2322:   NODE b1,b2,br=0,br0;
                   2323:   NBM m1,m2,m;
                   2324:   P c;
                   2325:
                   2326:   if ( !p1 )
                   2327:     *rp = p2;
                   2328:   else if ( !p2 )
                   2329:     *rp = p1;
                   2330:   else {
                   2331:     for ( b1 = BDY(p1), b2 = BDY(p2), br0 = 0; b1 && b2; ) {
                   2332:       m1 = (NBM)BDY(b1); m2 = (NBM)BDY(b2);
                   2333:       switch ( comp_nbm(m1,m2) ) {
                   2334:         case 0:
                   2335:           addp(CO,m1->c,m2->c,&c);
                   2336:           if ( c ) {
                   2337:             NEXTNODE(br0,br);
                   2338:             NEWNBM(m); m->d = m1->d; m->c = c; m->b = m1->b;
                   2339:             BDY(br) = (pointer)m;
                   2340:           }
                   2341:           b1 = NEXT(b1); b2 = NEXT(b2); break;
                   2342:         case 1:
                   2343:           NEXTNODE(br0,br); BDY(br) = BDY(b1);
                   2344:           b1 = NEXT(b1); break;
                   2345:         case -1:
                   2346:           NEXTNODE(br0,br); BDY(br) = BDY(b2);
                   2347:           b2 = NEXT(b2); break;
                   2348:       }
                   2349:     }
                   2350:     if ( !br0 )
                   2351:       if ( b1 )
                   2352:         br0 = b1;
                   2353:       else if ( b2 )
                   2354:         br0 = b2;
                   2355:       else {
                   2356:         *rp = 0;
                   2357:         return;
                   2358:       }
                   2359:     else if ( b1 )
                   2360:       NEXT(br) = b1;
                   2361:     else if ( b2 )
                   2362:         NEXT(br) = b2;
                   2363:     else
                   2364:       NEXT(br) = 0;
                   2365:     MKNBP(*rp,br0);
                   2366:   }
                   2367: }
                   2368:
                   2369: void subnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2370: {
                   2371:   NBP t;
                   2372:
                   2373:   chsgnnbp(p2,&t);
                   2374:   addnbp(vl,p1,t,rp);
                   2375: }
                   2376:
                   2377: void chsgnnbp(NBP p,NBP *rp)
                   2378: {
                   2379:   NODE r0,r=0,b;
                   2380:   NBM m,m1;
                   2381:
                   2382:   for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
                   2383:     NEXTNODE(r0,r);
                   2384:     m = (NBM)BDY(b);
                   2385:     NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnp(m->c,&m1->c);
                   2386:     BDY(r) = m1;
                   2387:   }
                   2388:   if ( r0 ) NEXT(r) = 0;
                   2389:   MKNBP(*rp,r0);
                   2390: }
                   2391:
                   2392: void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2393: {
                   2394:   NODE b,n;
                   2395:   NBP r,t,s;
                   2396:   NBM m;
                   2397:
                   2398:   if ( !p1 || !p2 ) {
                   2399:     *rp = 0; return;
                   2400:   }
                   2401:   if ( OID(p1) != O_NBP ) {
                   2402:     if ( !POLY(p1) )
                   2403:       error("mulnbp : invalid argument");
                   2404:     NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
                   2405:     MKNODE(n,m,0); MKNBP(p1,n);
                   2406:   }
                   2407:   if ( OID(p2) != O_NBP ) {
                   2408:     if ( !POLY(p2) )
                   2409:       error("mulnbp : invalid argument");
                   2410:     NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
                   2411:     MKNODE(n,m,0); MKNBP(p2,n);
                   2412:   }
                   2413:   if ( length(BDY(p1)) < length(BDY(p2)) ) {
                   2414:     for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {
                   2415:       mulnbmnbp(vl,(NBM)BDY(b),p2,&t);
                   2416:       addnbp(vl,r,t,&s); r = s;
                   2417:     }
                   2418:     *rp = r;
                   2419:   } else {
                   2420:     for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) {
                   2421:       mulnbpnbm(vl,p1,(NBM)BDY(b),&t);
                   2422:       addnbp(vl,r,t,&s); r = s;
                   2423:     }
                   2424:     *rp = r;
                   2425:   }
                   2426: }
                   2427:
                   2428: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp)
                   2429: {
                   2430:   NODE b,r0,r=0;
                   2431:
                   2432:   if ( !p ) *rp = 0;
                   2433:   else {
                   2434:     for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
                   2435:       NEXTNODE(r0,r);
                   2436:       BDY(r) = mul_nbm(m,(NBM)BDY(b));
                   2437:     }
                   2438:     if ( r0 ) NEXT(r) = 0;
                   2439:     MKNBP(*rp,r0);
                   2440:   }
                   2441: }
                   2442:
                   2443: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp)
                   2444: {
                   2445:   NODE b,r0,r=0;
                   2446:
                   2447:   if ( !p ) *rp = 0;
                   2448:   else {
                   2449:     for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
                   2450:       NEXTNODE(r0,r);
                   2451:       BDY(r) = mul_nbm((NBM)BDY(b),m);
                   2452:     }
                   2453:     if ( r0 ) NEXT(r) = 0;
                   2454:     MKNBP(*rp,r0);
                   2455:   }
                   2456: }
                   2457:
                   2458: void pwrnbp(VL vl,NBP a,Z q,NBP *c)
                   2459: {
                   2460:   NBP a1,a2;
                   2461:   Z q1,r1,two;
                   2462:   NBM m;
                   2463:   NODE r;
                   2464:
                   2465:   if ( !q ) {
                   2466:      NEWNBM(m); m->d = 0; m->c = (P)ONE; m->b = 0;
                   2467:      MKNODE(r,m,0); MKNBP(*c,r);
                   2468:   } else if ( !a )
                   2469:     *c = 0;
                   2470:   else if ( UNIQ(q) )
                   2471:     *c = a;
                   2472:   else {
1.2       noro     2473:     STOZ(2,two);
1.1       noro     2474:     divqrz(q,two,&q1,&r1);
                   2475:     pwrnbp(vl,a,q1,&a1);
                   2476:     mulnbp(vl,a1,a1,&a2);
                   2477:     if ( r1 )
                   2478:       mulnbp(vl,a2,a,c);
                   2479:     else
                   2480:       *c = a2;
                   2481:   }
                   2482: }
                   2483:
                   2484: int compnbp(VL vl,NBP p1,NBP p2)
                   2485: {
                   2486:   NODE n1,n2;
                   2487:   NBM m1,m2;
                   2488:   int t;
                   2489:
                   2490:   if ( !p1 )
                   2491:     return p2 ? -1 : 0;
                   2492:   else if ( !p2 )
                   2493:     return 1;
                   2494:   else {
                   2495:     for ( n1 = BDY(p1), n2 = BDY(p2);
                   2496:       n1 && n2; n1 = NEXT(n1), n2 = NEXT(n2) ) {
                   2497:       m1 = (NBM)BDY(n1); m2 = (NBM)BDY(n2);
                   2498:       if ( (t = comp_nbm(m1,m2)) || (t = compp(CO,m1->c,m2->c) ) )
                   2499:         return t;
                   2500:     }
                   2501:     if ( n1 )
                   2502:       return 1;
                   2503:     else if ( n2 )
                   2504:       return -1;
                   2505:     else
                   2506:       return 0;
                   2507:   }
                   2508: }
                   2509:
                   2510: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2511: {
                   2512:   NODE b1,b2,n;
                   2513:   NBP r,t,s;
                   2514:   NBM m;
                   2515:
                   2516:   if ( !p1 || !p2 ) {
                   2517:     *rp = 0; return;
                   2518:   }
                   2519:   if ( OID(p1) != O_NBP ) {
                   2520:     if ( !POLY(p1) )
                   2521:       error("shuffle_mulnbp : invalid argument");
                   2522:     NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
                   2523:     MKNODE(n,m,0); MKNBP(p1,n);
                   2524:   }
                   2525:   if ( OID(p2) != O_NBP ) {
                   2526:     if ( !POLY(p2) )
                   2527:       error("shuffle_mulnbp : invalid argument");
                   2528:     NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
                   2529:     MKNODE(n,m,0); MKNBP(p2,n);
                   2530:   }
                   2531:   for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
                   2532:     for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
                   2533:       t = shuffle_mul_nbm(m,(NBM)BDY(b2));
                   2534:       addnbp(vl,r,t,&s); r = s;
                   2535:     }
                   2536:   *rp = r;
                   2537: }
                   2538:
                   2539: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2540: {
                   2541:   NODE b1,b2,n;
                   2542:   NBP r,t,s;
                   2543:   NBM m;
                   2544:
                   2545:   if ( !p1 || !p2 ) {
                   2546:     *rp = 0; return;
                   2547:   }
                   2548:   if ( OID(p1) != O_NBP ) {
                   2549:     if ( !POLY(p1) )
                   2550:       error("harmonic_mulnbp : invalid argument");
                   2551:     NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
                   2552:     MKNODE(n,m,0); MKNBP(p1,n);
                   2553:   }
                   2554:   if ( OID(p2) != O_NBP ) {
                   2555:     if ( !POLY(p2) )
                   2556:       error("harmonic_mulnbp : invalid argument");
                   2557:     NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
                   2558:     MKNODE(n,m,0); MKNBP(p2,n);
                   2559:   }
                   2560:   for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
                   2561:     for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
                   2562:       t = harmonic_mul_nbm(m,(NBM)BDY(b2));
                   2563:       addnbp(vl,r,t,&s); r = s;
                   2564:     }
                   2565:   *rp = r;
                   2566: }
                   2567:
                   2568: #if 0
                   2569: NBP shuffle_mul_nbm(NBM a,NBM b)
                   2570: {
                   2571:   int ad,bd,d,i,ai,bi,bit,s;
                   2572:   int *ab,*bb,*wmb,*w;
                   2573:   NBM wm,tm;
                   2574:   P c,c1;
                   2575:   NODE r,t,t1,p;
                   2576:   NBP u;
                   2577:
                   2578:   ad = a->d; bd = b->d; ab = a->b; bb = b->b;
                   2579:   d = ad + bd;
                   2580:   w = (int *)ALLOCA(d*sizeof(int));
                   2581:   NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   2582:   for ( i = 0; i < ad; i++ ) w[i] = 1;
                   2583:   for ( ; i < d; i++ ) w[i] = 0;
                   2584:   mulp(CO,a->c,b->c,&c);
                   2585:   r = 0;
                   2586:   do {
                   2587:     wm->d = d; wm->c = c;
                   2588:     ai = 0; bi = 0;
                   2589:     for ( i = 0; i < d; i++ ) {
                   2590:       if ( w[i] ) { bit = NBM_GET(ab,ai); ai++; }
                   2591:       else { bit = NBM_GET(bb,bi); bi++; }
                   2592:       if ( bit ) NBM_SET(wmb,i);
                   2593:       else NBM_CLR(wmb,i);
                   2594:     }
                   2595:     for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
                   2596:       tm = (NBM)BDY(t);
                   2597:       s = comp_nbm(tm,wm);
                   2598:       if ( s < 0 ) {
                   2599:         /* insert */
                   2600:         MKNODE(t1,wm,t);
                   2601:         if ( !p ) r = t1;
                   2602:         else NEXT(p) = t1;
                   2603:         NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   2604:         break;
                   2605:       } else if ( s == 0 ) {
                   2606:         /* add coefs */
                   2607:         addp(CO,tm->c,c,&c1);
                   2608:         if ( c1 ) tm->c = c1;
                   2609:         else NEXT(p) = NEXT(t);
                   2610:         break;
                   2611:       }
                   2612:     }
                   2613:     if ( !t ) {
                   2614:       /* append */
                   2615:       MKNODE(t1,wm,t);
                   2616:       if ( !p ) r = t1;
                   2617:       else NEXT(p) = t1;
                   2618:       NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   2619:     }
                   2620:   } while ( ni_next(w,d) );
                   2621:   MKNBP(u,r);
                   2622:   return u;
                   2623: }
                   2624:
                   2625: int nbmtoxky(NBM a,int *b)
                   2626: {
                   2627:   int d,i,j,k;
                   2628:   int *p;
                   2629:
                   2630:   d = a->d; p = a->b;
                   2631:   for ( i = j = 0, k = 1; i < d; i++ ) {
                   2632:     if ( !NBM_GET(p,i) ) {
                   2633:       b[j++] = k;
                   2634:       k = 1;
                   2635:     } else k++;
                   2636:   }
                   2637:   return j;
                   2638: }
                   2639:
                   2640: NBP harmonic_mul_nbm(NBM a,NBM b)
                   2641: {
                   2642:   int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1;
                   2643:   int i,j,k,ia,ib,s;
                   2644:   int *wa,*wb,*w,*wab,*wa1,*wmb;
                   2645:   P c,c1;
                   2646:   NBM wm,tm;
                   2647:   NODE r,t1,t,p;
                   2648:   NBP u;
                   2649:
                   2650:   da = a->d; db = b->d; d = da+db;
                   2651:   wa = (int *)ALLOCA(da*sizeof(int));
                   2652:   wb = (int *)ALLOCA(db*sizeof(int));
                   2653:   la = nbmtoxky(a,wa);
                   2654:   lb = nbmtoxky(b,wb);
                   2655:   mulp(CO,a->c,b->c,&c);
                   2656:   /* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */
                   2657:   /* lmax : total length */
                   2658:   lmax = la+lb;
                   2659:   lmin = la>lb?la:lb;
                   2660:   w = (int *)ALLOCA(lmax*sizeof(int));
                   2661:   /* position of a+b */
                   2662:   wab = (int *)ALLOCA(lmax*sizeof(int));
                   2663:   /* position of a */
                   2664:   wa1 = (int *)ALLOCA(lmax*sizeof(int));
                   2665:   NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   2666:   for ( l = lmin, r = 0; l <= lmax; l++ ) {
                   2667:     lab = lmax - l;
                   2668:     la1 = la - lab;
                   2669:     lb1 = lb - lab;
                   2670:     lab1 = l-lab;
                   2671:     /* partion l into three parts: a, b, a+b */
                   2672:     /* initialize wab */
                   2673:     for ( i = 0; i < lab; i++ ) wab[i] = 1;
                   2674:     for ( ; i < l; i++ ) wab[i] = 0;
                   2675:     do {
                   2676:       /* initialize wa1 */
                   2677:       for ( i = 0; i < la1; i++ ) wa1[i] = 1;
                   2678:       for ( ; i < lab1; i++ ) wa1[i] = 0;
                   2679:       do {
                   2680:         ia = 0; ib = 0;
                   2681:         for ( i = j = 0; i < l; i++ )
                   2682:           if ( wab[i] ) w[i] = wa[ia++]+wb[ib++];
                   2683:           else if ( wa1[j++] ) w[i] = wa[ia++];
                   2684:           else w[i] = wb[ib++];
                   2685:         for ( i = j = 0; i < l; i++ ) {
                   2686:           for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j);
                   2687:           NBM_CLR(wmb,j); j++;
                   2688:         }
                   2689:         wm->d = j; wm->c = c;
                   2690:         for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
                   2691:           tm = (NBM)BDY(t);
                   2692:           s = comp_nbm(tm,wm);
                   2693:           if ( s < 0 ) {
                   2694:             /* insert */
                   2695:             MKNODE(t1,wm,t);
                   2696:             if ( !p ) r = t1;
                   2697:             else NEXT(p) = t1;
                   2698:             NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   2699:             break;
                   2700:           } else if ( s == 0 ) {
                   2701:             /* add coefs */
                   2702:             addp(CO,tm->c,c,&c1);
                   2703:             if ( c1 ) tm->c = c1;
                   2704:             else NEXT(p) = NEXT(t);
                   2705:             break;
                   2706:           }
                   2707:         }
                   2708:         if ( !t ) {
                   2709:           /* append */
                   2710:           MKNODE(t1,wm,t);
                   2711:           if ( !p ) r = t1;
                   2712:           else NEXT(p) = t1;
                   2713:           NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   2714:         }
                   2715:       } while ( ni_next(wa1,lab1) );
                   2716:     } while ( ni_next(wab,l) );
                   2717:   }
                   2718:   MKNBP(u,r);
                   2719:   return u;
                   2720: }
                   2721: #endif
                   2722:
                   2723: /* DPM functions */
                   2724:
1.3       noro     2725: DMMstack dmm_stack;
1.9       noro     2726: int dpm_base_ordtype;;
1.3       noro     2727:
1.9       noro     2728: DMMstack push_schreyer_order(LIST data,DMMstack stack)
1.3       noro     2729: {
                   2730:   DMMstack t;
1.14      noro     2731:   DP dp;
                   2732:   MP mp;
                   2733:   DMM m0,m1;
                   2734:   DPM dpm0,dpm1;
                   2735:   int len,i,nv;
1.7       noro     2736:   NODE in,t1;
1.9       noro     2737:   LIST l;
1.4       noro     2738:
                   2739:   /* data = [DPM,...,DPM] */
1.9       noro     2740:   if ( !dmm_stack && ( !dp_current_spec || dp_current_spec->id < 256 ) )
                   2741:       error("push_schreyer_order : base module order is not set");
1.4       noro     2742:   in = BDY(data);
                   2743:   len = length(in);
                   2744:   NEWDMMstack(t);
                   2745:   t->rank = len;
                   2746:   t->in = (DMM *)MALLOC((len+1)*sizeof(DMM));
1.14      noro     2747:   t->sum = (DMM *)MALLOC((len+1)*sizeof(DMM));
1.9       noro     2748:   if ( stack ) {
                   2749:     MKNODE(t1,data,BDY(stack->obj)); MKLIST(l,t1); t->obj = l;
1.14      noro     2750:     for ( i = 1; i <= len; i++, in = NEXT(in) ) {
                   2751:        m1 = t->in[i] = BDY((DPM)BDY(in));
                   2752:        NEWMP(mp); mp->dl = m1->dl; mp->c = m1->c; NEXT(mp) = 0;
                   2753:        nv = ((DPM)BDY(in))->nv;
                   2754:        MKDP(nv,mp,dp); dp->sugar = mp->dl->td;
                   2755:        m0 = stack->sum[m1->pos]; MKDPM(nv,m0,dpm0);
                   2756:        mulobjdpm(CO,(Obj)dp,dpm0,&dpm1);
                   2757:        t->sum[i] = BDY(dpm1);
                   2758:     }
1.9       noro     2759:   } else {
                   2760:     MKNODE(t1,data,0); MKLIST(l,t1); t->obj = l;
1.14      noro     2761:     for ( i = 1; i <= len; i++, in = NEXT(in) ) {
                   2762:       t->sum[i] = t->in[i] = BDY((DPM)BDY(in));
                   2763:     }
1.4       noro     2764:   }
1.9       noro     2765:   t->next = stack;;
                   2766:   dpm_ordtype = 3;
                   2767:   return t;
1.4       noro     2768: }
                   2769:
                   2770: // data=[Ink,...,In0]
                   2771: // Ini = a list of module monomials
                   2772:
                   2773: void set_schreyer_order(LIST data)
                   2774: {
                   2775:   NODE in;
                   2776:   LIST *w;
                   2777:   int i,len;
1.3       noro     2778:
                   2779:   if ( !data ) {
                   2780:     dmm_stack = 0;
                   2781:     if ( dp_current_spec && dp_current_spec->id >= 256 )
1.9       noro     2782:       dpm_ordtype = dp_current_spec->module_ordtype;
1.3       noro     2783:     else
                   2784:       dpm_ordtype = 0;
                   2785:     return;
                   2786:   } else {
1.9       noro     2787:     if ( !dp_current_spec || dp_current_spec->id < 256 )
                   2788:       error("set_schreyer_order : base module order is not set");
1.4       noro     2789:     dmm_stack = 0;
1.9       noro     2790:     dpm_base_ordtype = dp_current_spec->module_ordtype;
1.4       noro     2791:     in = BDY(data);
1.3       noro     2792:     len = length(in);
1.4       noro     2793:     w = (LIST *)MALLOC(len*sizeof(LIST));
                   2794:     for ( i = 0; i < len; i++, in = NEXT(in) ) w[i] = (LIST)BDY(in);
1.9       noro     2795:     for ( i = len-1; i >= 0; i-- ) dmm_stack = push_schreyer_order(w[i],dmm_stack);
1.3       noro     2796:   }
                   2797: }
                   2798:
1.15      noro     2799: void set_schreyer_level(DMMstack_array array,int level)
                   2800: {
                   2801:   if ( !level ) {
                   2802:     dmm_stack = 0;
                   2803:     if ( dp_current_spec && dp_current_spec->id >= 256 )
                   2804:       dpm_ordtype = dp_current_spec->module_ordtype;
                   2805:     else
                   2806:       dpm_ordtype = 0;
                   2807:     return;
                   2808:   } else {
                   2809:     if ( !dp_current_spec || dp_current_spec->id < 256 )
                   2810:       error("set_schreyer_level : base module order is not set");
                   2811:     if ( level > array->len )
                   2812:       error("set_schreyer_level : invalid level");
                   2813:     dmm_stack = array->body[level-1];
                   2814:     dpm_base_ordtype = dp_current_spec->module_ordtype;
                   2815:     dpm_ordtype = 3;
                   2816:   }
                   2817: }
                   2818:
1.4       noro     2819: // construct a base of syz(g)
                   2820: // assuming the schrerer order is properly set
                   2821:
                   2822: DP dpm_sp_hm(DPM p1,DPM p2);
                   2823: void dpm_sp(DPM p1,DPM p2,DPM *sp,DP *t1,DP *t2);
1.12      noro     2824: DPM dpm_nf_and_quotient3(DPM sp,VECT psv,DPM *nf,P *dn);
                   2825: DPM dpm_nf_and_quotient4(DPM sp,DPM *ps,VECT psiv,DPM head,DPM *nf,P *dn);
                   2826: DPM dpm_sp_nf(VECT psv,VECT psiv,int i,int j,DPM *nf);
1.19      noro     2827: DPM dpm_sp_nf_zlist(VECT psv,VECT psiv,int i,int j,int top,DPM *nf);
1.13      noro     2828: DPM dpm_sp_nf_asir(VECT psv,int i,int j,DPM *nf);
1.4       noro     2829: void dpm_sort(DPM p,DPM *r);
                   2830:
                   2831: extern int DP_Multiple;
1.22      noro     2832: extern int DP_Print;
1.4       noro     2833:
                   2834: void dpm_nf_z(NODE b,DPM g,VECT psv,int full,int multiple,DPM *rp);
                   2835: NODE dpm_sort_list(NODE l);
                   2836: void dpm_ptozp(DPM p,Z *cont,DPM *r);
                   2837:
                   2838: NODE dpm_reduceall(NODE in)
                   2839: {
                   2840:   int n,i;
                   2841:   VECT psv;
                   2842:   DPM *ps;
                   2843:   NODE t,t1;
                   2844:   DPM g,r;
                   2845:   Z cont;
                   2846:
                   2847:   n = length(in);
                   2848:   MKVECT(psv,n);
                   2849:   ps = (DPM *)BDY(psv);
                   2850:   for ( i = 0, t = in; i < n; i++, t = NEXT(t) ) ps[i] = BDY(t);
                   2851:   for ( i = 0; i < n; i++ ) {
                   2852:     g = ps[i]; ps[i] = 0;
1.8       noro     2853: //    dpm_nf_z(0,g,psv,1,DP_Multiple,&ps[i]);
                   2854:     dpm_nf_z(0,g,psv,1,0,&ps[i]);
1.4       noro     2855:   }
                   2856:   t = 0;
                   2857:   for ( i = n-1; i >= 0; i-- ) {
                   2858:     dpm_ptozp(ps[i],&cont,&r);
                   2859:     MKNODE(t1,r,t); t = t1;
                   2860:   }
                   2861:   return t;
                   2862: }
                   2863:
1.8       noro     2864: struct oEGT egra;
                   2865:
1.9       noro     2866: void dpm_ht(DPM d,DPM *r);
1.16      noro     2867: NODE reverse_node(NODE b);
1.9       noro     2868:
1.12      noro     2869: void dpm_schreyer_base(LIST g,LIST *s)
                   2870: {
                   2871:   NODE nd,t0,t,b0,b;
                   2872:   int n,i,j,k,nv,max,pos;
                   2873:   LIST l;
                   2874:   DP h,t1,t2;
                   2875:   MP d;
                   2876:   DMM r0,r,r1;
                   2877:   DPM sp,nf,dpm;
                   2878:   DPM *ps;
                   2879:   VECT psv,psiv;
                   2880:   DPM quo;
1.16      noro     2881:   DP *mm;
1.12      noro     2882:   NODE *psi;
1.13      noro     2883:   NODE n1,n2,n3;
                   2884:   int p1,p2,p3;
                   2885:   struct oEGT eg0,eg1,egsp,egnf;
1.12      noro     2886:   extern struct oEGT egred;
1.15      noro     2887:   extern int sch_count,schrec_count,schlast_count;
1.12      noro     2888:
1.15      noro     2889:   sch_count = schlast_count= 0;
1.12      noro     2890:   init_eg(&egra);
1.13      noro     2891:   init_eg(&egsp);
                   2892:   init_eg(&egnf);
1.12      noro     2893:   nd = BDY(g);
                   2894:   n = length(nd);
                   2895:   MKVECT(psv,n+1);
                   2896:   ps = (DPM *)BDY(psv);
                   2897:   for ( i = 1, t = nd; i <= n; i++, t = NEXT(t) ) ps[i] = (DPM)BDY(t);
                   2898:   for ( i = 1, max = 0; i <= n; i++ )
                   2899:     if ( (pos=BDY(ps[i])->pos) > max ) max = pos;
                   2900:   MKVECT(psiv,max+1);
                   2901:   psi = (NODE *)BDY(psiv);
                   2902:   for ( i = n; i >= 1; i-- ) {
                   2903:     pos = BDY(ps[i])->pos;
                   2904:     MKNODE(nd,(long)i,psi[pos]); psi[pos] = nd;
                   2905:   }
                   2906:   nv = ps[1]->nv;
1.16      noro     2907:   mm = (DP *)MALLOC((n+1)*sizeof(DP));
1.12      noro     2908:   b0 = 0;
1.13      noro     2909:   get_eg(&eg0);
                   2910:   for ( i = 1; i <= max; i++ ) {
1.16      noro     2911:     memset(mm,0,(n+1)*sizeof(DP));
1.13      noro     2912:     for ( n1 = psi[i]; n1; n1 = NEXT(n1) ) {
                   2913:       p1 = (long)BDY(n1);
                   2914:       for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
                   2915:         p2 = (long)BDY(n2);
1.16      noro     2916:         mm[p2] = dpm_sp_hm(ps[p1],ps[p2]);
1.13      noro     2917:       }
                   2918:       for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
                   2919:         p2 = (long)BDY(n2);
1.16      noro     2920:         if ( !mm[p2] ) continue;
                   2921:         for ( h = mm[p2], n3 = NEXT(n1); n3; n3 = NEXT(n3) ) {
1.13      noro     2922:           p3 = (long)BDY(n3);
1.16      noro     2923:           if ( n3 != n2 && mm[p3] && dp_redble(mm[p3],h) ) mm[p3] = 0;
                   2924:         }
                   2925:       }
                   2926:       for ( j = p1+1; j <= n; j++ ) {
                   2927:         if ( mm[j] ) {
                   2928:           quo = dpm_sp_nf(psv,psiv,p1,j,&nf);
                   2929:           if ( nf )
                   2930:             error("dpm_schreyer_base : cannot happen");
                   2931:           NEXTNODE(b0,b); BDY(b) = (pointer)quo;
1.13      noro     2932:         }
                   2933:       }
                   2934:     }
                   2935:   }
                   2936:   get_eg(&eg1); add_eg(&egsp,&eg0,&eg1); print_eg("SP",&egsp);
                   2937:   get_eg(&eg0);
                   2938:   get_eg(&eg1); add_eg(&egnf,&eg0,&eg1); print_eg("NF",&egnf); printf("\n");
1.12      noro     2939:   if ( b0 ) NEXT(b) = 0;
1.17      noro     2940:   for ( t0 = 0, nd = BDY(g); nd; nd = NEXT(nd) ) {
1.12      noro     2941:     dpm_ht((DPM)BDY(nd),&dpm); NEXTNODE(t0,t); BDY(t) = (pointer)dpm;
                   2942:   }
                   2943:   if ( t0 ) NEXT(t) = 0;
                   2944:   MKLIST(l,t0);
                   2945:   dmm_stack = push_schreyer_order(l,dmm_stack);
1.16      noro     2946: //  b0 = dpm_sort_list(b0);
1.12      noro     2947: //  get_eg(&eg0);
                   2948: //  b0 = dpm_reduceall(b0);
                   2949: //  get_eg(&eg1); add_eg(&egra,&eg0,&eg1); print_eg("RA",&egra);
1.4       noro     2950:   MKLIST(*s,b0);
1.8       noro     2951: //  print_eg("red",&egred); printf("\n");
1.15      noro     2952:   printf("sch_count=%d, schlast_count=%d\n",sch_count,schlast_count);
                   2953: }
                   2954:
1.17      noro     2955: void dpm_list_to_array(LIST g,VECT *psvect,VECT *psindvect)
                   2956: {
                   2957:   NODE nd,t;
                   2958:   int n;
                   2959:   VECT psv,psiv;
                   2960:   DPM *ps;
                   2961:   int i,max,pos;
                   2962:   LIST *psi;
                   2963:   LIST l;
                   2964:   Z iz;
                   2965:
                   2966:   nd = BDY(g);
                   2967:   n = length(nd);
                   2968:   MKVECT(psv,n+1);
                   2969:   ps = (DPM *)BDY(psv);
                   2970:   for ( i = 1, t = nd; i <= n; i++, t = NEXT(t) ) ps[i] = (DPM)BDY(t);
                   2971:   for ( i = 1, max = 0; i <= n; i++ )
                   2972:     if ( (pos=BDY(ps[i])->pos) > max ) max = pos;
                   2973:   MKVECT(psiv,max+1);
                   2974:   psi = (LIST *)BDY(psiv);
                   2975:   for ( i = 0; i <= max; i++ ) {
                   2976:     MKLIST(l,0); psi[i] = l;
                   2977:   }
                   2978:   for ( i = n; i >= 1; i-- ) {
                   2979:     pos = BDY(ps[i])->pos;
                   2980:     STOZ(i,iz); MKNODE(nd,iz,BDY(psi[pos])); BDY(psi[pos]) = nd;
                   2981:   }
                   2982:   *psvect = psv; *psindvect = psiv;
                   2983: }
                   2984:
1.19      noro     2985: #if 0
1.17      noro     2986: void dpm_insert_to_zlist(VECT psiv,int pos,int i)
                   2987: {
                   2988:   LIST l;
                   2989:   NODE prev,cur,nd;
                   2990:   int curi;
                   2991:   Z iz;
                   2992:
                   2993:   l = (LIST)BDY(psiv)[pos];
                   2994:   for ( prev = 0, cur = BDY(l); cur; cur = NEXT(cur) ) {
                   2995:     curi = ZTOS((Q)BDY(cur));
                   2996:     if ( curi == i )
                   2997:       error("dpm_insert_to_list : invalid index");
                   2998:     if ( i < curi ) break;
                   2999:     prev = cur;
                   3000:   }
                   3001:   STOZ(i,iz); MKNODE(nd,iz,cur);
                   3002:   if ( prev == 0 ) BDY(l) = nd;
                   3003:   else NEXT(prev) = nd;
                   3004: }
1.19      noro     3005: #else
                   3006: void dpm_insert_to_zlist(VECT psiv,int pos,int i)
                   3007: {
                   3008:   LIST l;
                   3009:   NODE prev,cur,nd;
                   3010:   int curi;
                   3011:   Z iz;
                   3012:
                   3013:   l = (LIST)BDY(psiv)[pos];
                   3014:   for ( prev = 0, cur = BDY(l); cur; cur = NEXT(cur) ) {
                   3015:     curi = ZTOS((Q)BDY(cur));
                   3016:     if ( curi == i )
                   3017:       error("dpm_insert_to_list : invalid index");
                   3018:     prev = cur;
                   3019:   }
                   3020:   STOZ(i,iz); MKNODE(nd,iz,cur);
                   3021:   if ( prev == 0 ) BDY(l) = nd;
                   3022:   else NEXT(prev) = nd;
                   3023: }
                   3024: #endif
1.17      noro     3025:
                   3026: void dpm_schreyer_base_zlist(LIST g,LIST *s)
                   3027: {
                   3028:   NODE nd,t0,t,b0,b;
                   3029:   int n,i,j,k,nv,max,pos;
                   3030:   LIST l;
                   3031:   DP h,t1,t2;
                   3032:   MP d;
                   3033:   DMM r0,r,r1;
                   3034:   DPM sp,nf,dpm;
                   3035:   DPM *ps;
                   3036:   VECT psv,psiv;
                   3037:   DPM quo;
                   3038:   DP *mm;
                   3039:   LIST *psi;
                   3040:   NODE n1,n2,n3;
                   3041:   int p1,p2,p3;
                   3042:   Z iz;
                   3043:   struct oEGT eg0,eg1,egsp,egnf;
                   3044:   extern struct oEGT egred;
                   3045:   extern int sch_count,schrec_count,schlast_count;
                   3046:
                   3047:   sch_count = schlast_count= 0;
                   3048:   init_eg(&egra);
                   3049:   init_eg(&egsp);
                   3050:   init_eg(&egnf);
                   3051:   dpm_list_to_array(g,&psv,&psiv);
                   3052:   ps = (DPM *)BDY(psv);
                   3053:   psi = (LIST *)BDY(psiv);
                   3054:   nv = ps[1]->nv;
                   3055:   n = psv->len-1;
                   3056:   max = psiv->len-1;
                   3057:   mm = (DP *)MALLOC((n+1)*sizeof(DP));
                   3058:   b0 = 0;
                   3059:   get_eg(&eg0);
                   3060:   for ( i = 1; i <= max; i++ ) {
                   3061:     memset(mm,0,(n+1)*sizeof(DP));
                   3062:     for ( n1 = BDY((LIST)psi[i]); n1; n1 = NEXT(n1) ) {
                   3063:       p1 = ZTOS((Q)BDY(n1));
                   3064:       for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
                   3065:         p2 = ZTOS((Q)BDY(n2));
                   3066:         mm[p2] = dpm_sp_hm(ps[p1],ps[p2]);
                   3067:       }
                   3068:       for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
                   3069:         p2 = ZTOS((Q)BDY(n2));
                   3070:         if ( !mm[p2] ) continue;
                   3071:         for ( h = mm[p2], n3 = NEXT(n1); n3; n3 = NEXT(n3) ) {
                   3072:           p3 = ZTOS((Q)BDY(n3));
                   3073:           if ( n3 != n2 && mm[p3] && dp_redble(mm[p3],h) ) mm[p3] = 0;
                   3074:         }
                   3075:       }
                   3076:       for ( j = p1+1; j <= n; j++ ) {
                   3077:         if ( mm[j] ) {
1.19      noro     3078:           quo = dpm_sp_nf_zlist(psv,psiv,p1,j,0,&nf);
1.17      noro     3079:           if ( nf )
                   3080:             error("dpm_schreyer_base : cannot happen");
                   3081:           NEXTNODE(b0,b); BDY(b) = (pointer)quo;
                   3082:         }
                   3083:       }
                   3084:     }
                   3085:   }
                   3086:   get_eg(&eg1); add_eg(&egsp,&eg0,&eg1); print_eg("SP",&egsp);
                   3087:   get_eg(&eg0);
                   3088:   get_eg(&eg1); add_eg(&egnf,&eg0,&eg1); print_eg("NF",&egnf); printf("\n");
                   3089:   if ( b0 ) NEXT(b) = 0;
                   3090:   for ( t0 = 0, nd = BDY(g); nd; nd = NEXT(nd) ) {
                   3091:     dpm_ht((DPM)BDY(nd),&dpm); NEXTNODE(t0,t); BDY(t) = (pointer)dpm;
                   3092:   }
                   3093:   if ( t0 ) NEXT(t) = 0;
                   3094:   MKLIST(l,t0);
                   3095:   dmm_stack = push_schreyer_order(l,dmm_stack);
                   3096: //  b0 = dpm_sort_list(b0);
                   3097: //  get_eg(&eg0);
                   3098: //  b0 = dpm_reduceall(b0);
                   3099: //  get_eg(&eg1); add_eg(&egra,&eg0,&eg1); print_eg("RA",&egra);
                   3100:   MKLIST(*s,b0);
                   3101: //  print_eg("red",&egred); printf("\n");
                   3102:   printf("sch_count=%d, schlast_count=%d\n",sch_count,schlast_count);
                   3103: }
1.15      noro     3104:
1.19      noro     3105: static int compdp_nv;
                   3106:
                   3107: int compdp_revgradlex(DP *a,DP *b)
                   3108: {
                   3109:   return -cmpdl_revgradlex(compdp_nv,BDY(*a)->dl,BDY(*b)->dl);
                   3110: }
                   3111:
                   3112: int compdp_lex(DP *a,DP *b)
                   3113: {
                   3114:   return -cmpdl_lex(compdp_nv,BDY(*a)->dl,BDY(*b)->dl);
                   3115: }
                   3116:
1.22      noro     3117: DMMstack_array dpm_schreyer_frame(NODE g,int lex)
1.15      noro     3118: {
                   3119:   LIST l;
                   3120:   NODE nd,in,b0,b,n1,n2,n3,t;
                   3121:   NODE *psi;
                   3122:   long p1,p2,p3;
1.19      noro     3123:   int nv,n,i,j,k,max,pos,level;
1.15      noro     3124:   DMMstack s,s1;
                   3125:   DMM m1,m0,dmm;
                   3126:   MP mp;
                   3127:   DP dp,h;
1.19      noro     3128:   DP *m,*m2;
1.15      noro     3129:   DPM dpm,dpm0,dpm1;
                   3130:   VECT psv,psiv;
                   3131:   DPM *ps;
                   3132:   DMMstack_array dmmstack_array;
                   3133:
                   3134:   nd = g;
1.19      noro     3135:   compdp_nv = nv = ((DPM)BDY(nd))->nv;
1.15      noro     3136:   s = 0;
                   3137:   level = 0;
                   3138:   while ( 1 ) {
                   3139:     /* store the current nd to a DMMstack */
                   3140:     n = length(nd);
                   3141:     NEWDMMstack(s1);
                   3142:     MKLIST(l,nd); s1->obj = l;
                   3143:     s1->rank = n;
                   3144:     s1->in = (DMM *)MALLOC((n+1)*sizeof(DMM));
                   3145:     s1->sum = (DMM *)MALLOC((n+1)*sizeof(DMM));
                   3146:     NEXT(s1) = s;
                   3147:     if ( s ) {
                   3148:       for ( i = 1, in = nd; i <= n; i++, in = NEXT(in) ) {
                   3149:         m1 = s1->in[i] = BDY((DPM)BDY(in));
                   3150:         NEWMP(mp); mp->dl = m1->dl; mp->c = m1->c; NEXT(mp) = 0;
                   3151:         MKDP(nv,mp,dp); dp->sugar = mp->dl->td;
                   3152:         m0 = s->sum[m1->pos]; MKDPM(nv,m0,dpm0);
                   3153:         mulobjdpm(CO,(Obj)dp,dpm0,&dpm1);
                   3154:         s1->sum[i] = BDY(dpm1);
                   3155:       }
                   3156:     } else {
1.18      noro     3157:       for ( i = 1, in = nd; i <= n; i++, in = NEXT(in) )  {
                   3158:         m0 = BDY((DPM)BDY(in));
                   3159:         NEWDMM(m1); *m1 = *m0; m1->c = (Obj)ONE; NEXT(m1) = 0;
                   3160:         s1->sum[i] = s1->in[i] = m1;
                   3161:       }
1.15      noro     3162:     }
                   3163:     s = s1;
                   3164:     level++;
1.22      noro     3165:     if ( DP_Print ) printf("level=%d,len=%d\n",level,n);
1.15      noro     3166:
                   3167:     /* create new list */
                   3168:     MKVECT(psv,n+1);
                   3169:     ps = (DPM *)BDY(psv);
                   3170:     for ( i = 1, t = nd; i <= n; i++, t = NEXT(t) ) ps[i] = (DPM)BDY(t);
                   3171:     for ( i = 1, max = 0; i <= n; i++ )
                   3172:       if ( (pos=BDY(ps[i])->pos) > max ) max = pos;
                   3173:     MKVECT(psiv,max+1);
                   3174:     psi = (NODE *)BDY(psiv);
                   3175:     for ( i = n; i >= 1; i-- ) {
                   3176:       pos = BDY(ps[i])->pos;
                   3177:       MKNODE(nd,(long)i,psi[pos]); psi[pos] = nd;
                   3178:     }
                   3179:     m = (DP *)MALLOC((n+1)*sizeof(DP));
1.19      noro     3180:     m2 = (DP *)MALLOC((n+1)*sizeof(DP));
1.15      noro     3181:     b0 = 0;
                   3182:     for ( i = 1; i <= max; i++ ) {
                   3183:       for ( n1 = psi[i]; n1; n1 = NEXT(n1) ) {
                   3184:         bzero(m,(n+1)*sizeof(DP));
                   3185:         p1 = (long)BDY(n1);
                   3186:         for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
                   3187:           p2 = (long)BDY(n2);
                   3188:           h = dpm_sp_hm(ps[p1],ps[p2]);
                   3189:           for ( n3 = NEXT(n1); n3 != n2; n3 = NEXT(n3) ) {
                   3190:             p3 = (long)BDY(n3);
                   3191:             if ( m[p3] ) {
                   3192:               if ( dp_redble(h,m[p3]) ) {
                   3193:                 h = 0; break;
                   3194:               }
                   3195:               if ( dp_redble(m[p3],h) ) m[p3] = 0;
                   3196:             }
                   3197:           }
                   3198:           if ( h ) m[p2] = h;
                   3199:         }
1.22      noro     3200:         if ( lex ) {
                   3201:           // compress m to m2
                   3202:           for ( j = 0, n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
                   3203:             p2 = (long)BDY(n2);
                   3204:             if ( m[p2] ) m2[j++] = m[p2];
                   3205:           }
                   3206:           qsort(m2,j,sizeof(DP),(int (*)(const void *,const void *))compdp_lex);
                   3207:           for ( k = 0; k < j; k++ ) {
                   3208:             NEWDMM(dmm); dmm->dl = BDY(m2[k])->dl; dmm->pos = p1; dmm->c = (Obj)ONE;
1.15      noro     3209:             MKDPM(nv,dmm,dpm);
                   3210:             NEXTNODE(b0,b); BDY(b) = (pointer)dpm;
                   3211:           }
1.22      noro     3212:         } else {
                   3213:           for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
                   3214:             p2 = (long)BDY(n2);
                   3215:             if ( m[p2] ) {
                   3216:               NEWDMM(dmm); dmm->dl = BDY(m[p2])->dl; dmm->pos = p1; dmm->c = (Obj)ONE;
                   3217:               MKDPM(nv,dmm,dpm);
                   3218:               NEXTNODE(b0,b); BDY(b) = (pointer)dpm;
                   3219:             }
                   3220:           }
1.15      noro     3221:         }
                   3222:       }
                   3223:     }
                   3224:     if ( !b0 ) {
                   3225:       NEWDMMstack_array(dmmstack_array);
                   3226:       dmmstack_array->len = level;
                   3227:       dmmstack_array->body = (DMMstack *)MALLOC(level*sizeof(DMMstack));
                   3228:       for ( i = level-1, s1 = s; i >= 0; i--, s1 = NEXT(s1) )
                   3229:         dmmstack_array->body[i] = s1;
                   3230:       return dmmstack_array;
                   3231:     } else {
                   3232:       NEXT(b) = 0;
                   3233:       nd = b0;
                   3234:     }
                   3235:   }
1.4       noro     3236: }
1.15      noro     3237:
                   3238: int sch_count,schlast_count;
1.4       noro     3239:
1.3       noro     3240: int compdmm_schreyer(int n,DMM m1,DMM m2)
                   3241: {
1.15      noro     3242:   int pos1,pos2,t,npos1,npos2;
                   3243:   DMM *in,*sum;
                   3244:   DMMstack s;
                   3245:   static DL d1=0,d2=0;
                   3246:   static int dlen=0;
                   3247:
                   3248:   sch_count++;
                   3249:   pos1 = m1->pos; pos2 = m2->pos;
                   3250:   if ( pos1 == pos2 ) return (*cmpdl)(n,m1->dl,m2->dl);
                   3251:   if ( n > dlen ) {
                   3252:     NEWDL(d1,n); NEWDL(d2,n); dlen = n;
                   3253:   }
                   3254:   sum = dmm_stack->sum;
                   3255:   _adddl(n,m1->dl,sum[pos1]->dl,d1);
                   3256:   _adddl(n,m2->dl,sum[pos2]->dl,d2);
                   3257:   t = (*cmpdl)(n,d1,d2);
                   3258:   if ( sum[pos1]->pos == sum[pos2]->pos && t == 0 ) {
                   3259:     for ( s = dmm_stack; s; s = NEXT(s) ) {
                   3260:       in = s->in;
                   3261:       npos1 = in[pos1]->pos;
                   3262:       npos2 = in[pos2]->pos;
                   3263:       if ( npos1 == npos2 ) break;
                   3264:       else {
                   3265:         pos1 = npos1;
                   3266:         pos2 = npos2;
                   3267:       }
                   3268:     }
                   3269:     // (pos1,pos2) = the last pair s.t. pos1 != pos2
                   3270:     // simply compare pos1 and pos2
1.3       noro     3271:     if ( pos1 < pos2 ) return 1;
                   3272:     else if ( pos1 > pos2 ) return -1;
                   3273:   } else {
1.15      noro     3274:     schlast_count++;
                   3275:     pos1 = sum[pos1]->pos;
                   3276:     pos2 = sum[pos2]->pos;
                   3277:     if ( dpm_base_ordtype == 1 ) {
                   3278:       if ( pos1 < pos2 ) return 1;
                   3279:       else if ( pos1 > pos2 ) return -1;
                   3280:       else return t;
                   3281:     } else {
                   3282:       if ( t ) return t;
                   3283:       else if ( pos1 < pos2 ) return 1;
                   3284:       else if ( pos1 > pos2 ) return -1;
                   3285:       else return 0;
                   3286:     }
1.3       noro     3287:   }
1.24      noro     3288:   /* XXX */
                   3289:   return 0;
1.3       noro     3290: }
                   3291:
1.20      noro     3292: int compdmm_schreyer_old(int n,DMM m1,DMM m2)
                   3293: {
                   3294:   int pos1,pos2,t,npos1,npos2;
                   3295:   DMM *in,*sum;
                   3296:   DMMstack s;
                   3297:   static DL d1=0,d2=0;
                   3298:   static int dlen=0;
                   3299:
                   3300:   sch_count++;
                   3301:   pos1 = m1->pos; pos2 = m2->pos;
                   3302:   if ( pos1 == pos2 ) return (*cmpdl)(n,m1->dl,m2->dl);
                   3303:   if ( n > dlen ) {
                   3304:     NEWDL(d1,n); NEWDL(d2,n); dlen = n;
                   3305:   }
                   3306:   sum = dmm_stack->sum;
                   3307:   _copydl(n,m1->dl,d1);
                   3308:   _copydl(n,m2->dl,d2);
                   3309:   for ( s = dmm_stack; s; s = NEXT(s) ) {
                   3310:     in = s->in;
                   3311:     _addtodl(n,in[pos1]->dl,d1);
                   3312:     _addtodl(n,in[pos2]->dl,d2);
                   3313:     if ( in[pos1]->pos == in[pos2]->pos && _eqdl(n,d1,d2)) {
                   3314:       if ( pos1 < pos2 ) return 1;
                   3315:       else if ( pos1 > pos2 ) return -1;
                   3316:       else return 0;
                   3317:     }
                   3318:     pos1 = in[pos1]->pos;
                   3319:     pos2 = in[pos2]->pos;
                   3320:     if ( pos1 == pos2 ) return (*cmpdl)(n,d1,d2);
                   3321:   }
                   3322:   if ( dpm_base_ordtype == 1 ) {
                   3323:     if ( pos1 < pos2 ) return 1;
                   3324:     else if ( pos1 > pos2 ) return -1;
                   3325:     else return (*cmpdl)(n,d1,d2);
                   3326:   } else {
                   3327:     t = (*cmpdl)(n,d1,d2);
                   3328:     if ( t ) return t;
                   3329:     else if ( pos1 < pos2 ) return 1;
                   3330:     else if ( pos1 > pos2 ) return -1;
                   3331:     else return 0;
                   3332:   }
                   3333: }
                   3334:
                   3335: extern int NaiveSchreyer;
                   3336:
1.4       noro     3337: int compdmm(int n,DMM m1,DMM m2)
                   3338: {
1.20      noro     3339:   int t,t1;
1.19      noro     3340:   int *base_ord;
1.4       noro     3341:
                   3342:   switch ( dpm_ordtype ) {
1.9       noro     3343:   case 0: /* TOP ord->pos */
1.4       noro     3344:     t = (*cmpdl)(n,m1->dl,m2->dl);
                   3345:     if ( t ) return t;
                   3346:     else if ( m1->pos < m2->pos ) return 1;
                   3347:     else if ( m1->pos > m2->pos ) return -1;
                   3348:     else return 0;
1.9       noro     3349:   case 1: /* POT : pos->ord */
1.4       noro     3350:     if ( m1->pos < m2->pos ) return 1;
                   3351:     else if ( m1->pos > m2->pos ) return -1;
                   3352:     else return (*cmpdl)(n,m1->dl,m2->dl);
1.9       noro     3353:   case 2:  /* wPOT: weight->pos->ord */
                   3354:     if ( m1->dl->td > m2->dl->td ) return 1;
                   3355:     else if ( m1->dl->td < m2->dl->td ) return 1;
                   3356:     else if ( m1->pos < m2->pos ) return 1;
                   3357:     else if ( m1->pos > m2->pos ) return -1;
                   3358:     else return (*cmpdl)(n,m1->dl,m2->dl);
                   3359:   case 3: /* Schreyer */
1.20      noro     3360:     if ( NaiveSchreyer )
                   3361:       t = compdmm_schreyer_old(n,m1,m2);
                   3362:     else
                   3363:       t = compdmm_schreyer(n,m1,m2);
                   3364:     return t;
1.19      noro     3365:   case 4:  /* POT with base_ord */
                   3366:     base_ord = dp_current_spec->module_base_ord;
                   3367:     if ( base_ord[m1->pos] < base_ord[m2->pos] ) return 1;
                   3368:     else if ( base_ord[m1->pos] > base_ord[m2->pos] ) return -1;
                   3369:     else return (*cmpdl)(n,m1->dl,m2->dl);
1.4       noro     3370:   default:
                   3371:     error("compdmm : invalid dpm_ordtype");
1.24      noro     3372:     return 0;
1.4       noro     3373:   }
                   3374: }
1.1       noro     3375:
                   3376: void adddpm(VL vl,DPM p1,DPM p2,DPM *pr)
                   3377: {
1.4       noro     3378:   int n,s;
1.1       noro     3379:   DMM m1,m2,mr=0,mr0;
                   3380:   Obj t;
                   3381:   DL d;
                   3382:
                   3383:   if ( !p1 )
                   3384:     *pr = p2;
                   3385:   else if ( !p2 )
                   3386:     *pr = p1;
                   3387:   else {
1.4       noro     3388:     for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
                   3389:       s = compdmm(n,m1,m2);
                   3390:       switch ( s ) {
1.1       noro     3391:         case 0:
                   3392:           arf_add(vl,C(m1),C(m2),&t);
                   3393:           if ( t ) {
                   3394:             NEXTDMM(mr0,mr); mr->pos = m1->pos; mr->dl = m1->dl; C(mr) = t;
                   3395:           }
                   3396:           m1 = NEXT(m1); m2 = NEXT(m2); break;
                   3397:         case 1:
                   3398:           NEXTDMM(mr0,mr); mr->pos = m1->pos; mr->dl = m1->dl; C(mr) = C(m1);
                   3399:           m1 = NEXT(m1); break;
                   3400:         case -1:
                   3401:           NEXTDMM(mr0,mr); mr->pos = m2->pos; mr->dl = m2->dl; C(mr) = C(m2);
                   3402:           m2 = NEXT(m2); break;
                   3403:       }
1.4       noro     3404:     }
1.1       noro     3405:     if ( !mr0 )
                   3406:       if ( m1 )
                   3407:         mr0 = m1;
                   3408:       else if ( m2 )
                   3409:         mr0 = m2;
                   3410:       else {
                   3411:         *pr = 0;
                   3412:         return;
                   3413:       }
                   3414:     else if ( m1 )
                   3415:       NEXT(mr) = m1;
                   3416:     else if ( m2 )
                   3417:       NEXT(mr) = m2;
                   3418:     else
                   3419:       NEXT(mr) = 0;
                   3420:     MKDPM(NV(p1),mr0,*pr);
                   3421:     if ( *pr )
                   3422:       (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                   3423:   }
                   3424: }
                   3425:
                   3426: void subdpm(VL vl,DPM p1,DPM p2,DPM *pr)
                   3427: {
                   3428:   DPM t;
                   3429:
                   3430:   if ( !p2 )
                   3431:     *pr = p1;
                   3432:   else {
                   3433:     chsgndpm(p2,&t); adddpm(vl,p1,t,pr);
                   3434:   }
                   3435: }
                   3436:
                   3437: void chsgndpm(DPM p,DPM *pr)
                   3438: {
                   3439:   DMM m,mr=0,mr0;
                   3440:   Obj r;
                   3441:
                   3442:   if ( !p )
                   3443:     *pr = 0;
                   3444:   else {
                   3445:     for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   3446:       NEXTDMM(mr0,mr); arf_chsgn(C(m),&C(mr)); mr->pos = m->pos; mr->dl = m->dl;
                   3447:     }
                   3448:     NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
                   3449:     if ( *pr )
                   3450:       (*pr)->sugar = p->sugar;
                   3451:   }
                   3452: }
                   3453:
1.8       noro     3454: void mulcmp(Obj c,MP m)
                   3455: {
                   3456:   MP t;
                   3457:   Obj c1;
                   3458:
                   3459:   for ( t = m; t; t = NEXT(t) ) {
                   3460:     arf_mul(CO,c,C(t),&c1); C(t) = c1;
                   3461:   }
                   3462: }
                   3463:
                   3464: void mulcdmm(Obj c,DMM m)
                   3465: {
                   3466:   DMM t;
                   3467:   Obj c1;
                   3468:
                   3469:   for ( t = m; t; t = NEXT(t) ) {
                   3470:     arf_mul(CO,c,C(t),&c1); C(t) = c1;
                   3471:   }
                   3472: }
                   3473:
1.1       noro     3474: void mulcdpm(VL vl,Obj c,DPM p,DPM *pr)
                   3475: {
                   3476:   DMM m,mr=0,mr0;
                   3477:
                   3478:   if ( !p || !c )
                   3479:     *pr = 0;
                   3480:   else if ( NUM(c) && UNIQ((Q)c) )
                   3481:     *pr = p;
                   3482:   else if ( NUM(c) && MUNIQ((Q)c) )
                   3483:     chsgndpm(p,pr);
                   3484:   else {
                   3485:     for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   3486:       NEXTDMM(mr0,mr);
                   3487:       arf_mul(vl,C(m),c,&C(mr));
                   3488:       mr->pos = m->pos;
                   3489:       mr->dl = m->dl;
                   3490:     }
                   3491:     NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
                   3492:     if ( *pr )
                   3493:       (*pr)->sugar = p->sugar;
                   3494:   }
                   3495: }
                   3496:
                   3497: void comm_mulmpdpm(VL vl,MP m0,DPM p,DPM *pr)
                   3498: {
                   3499:   DMM m,mr=0,mr0;
                   3500:   DL d;
                   3501:   Obj c;
                   3502:   int n;
                   3503:
                   3504:   if ( !p )
                   3505:     *pr = 0;
                   3506:   else {
                   3507:     n = NV(p);
                   3508:     d = m0->dl;
                   3509:     c = C(m0);
                   3510:     for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   3511:       NEXTDMM(mr0,mr);
                   3512:       arf_mul(vl,C(m),c,&C(mr));
                   3513:       mr->pos = m->pos;
                   3514:       adddl(n,m->dl,d,&mr->dl);
                   3515:     }
                   3516:     NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
                   3517:     if ( *pr )
                   3518:       (*pr)->sugar = p->sugar;
                   3519:   }
                   3520: }
                   3521:
                   3522: void weyl_mulmpdpm(VL vl,MP m0,DPM p,DPM *pr)
                   3523: {
                   3524:   DPM r,t,t1;
                   3525:   DMM m;
                   3526:   DL d0;
                   3527:   int n,n2,l,i,j,tlen;
                   3528:   struct oMP mp;
                   3529:   static DMM *w,*psum;
                   3530:   static struct cdl *tab;
                   3531:   static int wlen;
                   3532:   static int rtlen;
                   3533:
                   3534:   if ( !p )
                   3535:     *pr = 0;
                   3536:   else {
                   3537:     for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
                   3538:     if ( l > wlen ) {
                   3539:       if ( w ) GCFREE(w);
                   3540:       w = (DMM *)MALLOC(l*sizeof(DMM));
                   3541:       wlen = l;
                   3542:     }
                   3543:     for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
                   3544:       w[i] = m;
                   3545:
                   3546:     n = NV(p); n2 = n>>1;
                   3547:     d0 = m0->dl;
                   3548:     for ( i = 0, tlen = 1; i < n2; i++ )
                   3549:       tlen *= d0->d[n2+i]+1;
                   3550:     if ( tlen > rtlen ) {
                   3551:       if ( tab ) GCFREE(tab);
                   3552:       if ( psum ) GCFREE(psum);
                   3553:       rtlen = tlen;
                   3554:       tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
                   3555:       psum = (DMM *)MALLOC(rtlen*sizeof(DMM));
                   3556:     }
                   3557:     bzero(psum,tlen*sizeof(DMM));
                   3558:     for ( i = l-1; i >= 0; i-- ) {
                   3559:       bzero(tab,tlen*sizeof(struct cdl));
                   3560:       mp.dl = w[i]->dl; mp.c = C(w[i]); mp.next = 0;
                   3561:       weyl_mulmm(vl,m0,&mp,n,tab,tlen);
                   3562:       for ( j = 0; j < tlen; j++ ) {
                   3563:         if ( tab[j].c ) {
                   3564:           NEWDMM(m); m->dl = tab[j].d; m->pos = w[i]->pos; C(m) = (Obj)tab[j].c; NEXT(m) = psum[j];
                   3565:           psum[j] = m;
                   3566:         }
                   3567:       }
                   3568:     }
                   3569:     for ( j = tlen-1, r = 0; j >= 0; j-- )
                   3570:       if ( psum[j] ) {
                   3571:         MKDPM(n,psum[j],t); adddpm(vl,r,t,&t1); r = t1;
                   3572:       }
                   3573:     if ( r )
                   3574:       r->sugar = p->sugar + m0->dl->td;
                   3575:     *pr = r;
                   3576:   }
                   3577: }
                   3578:
                   3579: void mulobjdpm(VL vl,Obj p1,DPM p2,DPM *pr)
                   3580: {
                   3581:   MP m;
                   3582:   DPM s,t,u;
                   3583:
                   3584:   if ( !p1 || !p2 )
                   3585:     *pr = 0;
                   3586:   else if ( OID(p1) != O_DP )
                   3587:     mulcdpm(vl,p1,p2,pr);
                   3588:   else {
                   3589:     s = 0;
                   3590:     for ( m = BDY((DP)p1); m; m = NEXT(m) ) {
                   3591:       if ( do_weyl )
                   3592:         weyl_mulmpdpm(vl,m,p2,&t);
                   3593:       else
                   3594:         comm_mulmpdpm(vl,m,p2,&t);
                   3595:       adddpm(vl,s,t,&u); s = u;
                   3596:     }
                   3597:     *pr = s;
                   3598:   }
                   3599: }
                   3600:
                   3601: int compdpm(VL vl,DPM p1,DPM p2)
                   3602: {
                   3603:   int n,t;
                   3604:   DMM m1,m2;
                   3605:
                   3606:   if ( !p1 )
                   3607:     return p2 ? -1 : 0;
                   3608:   else if ( !p2 )
                   3609:     return 1;
                   3610:   else if ( NV(p1) != NV(p2) ) {
                   3611:     error("compdpm : size mismatch");
                   3612:     return 0; /* XXX */
                   3613:   } else {
                   3614:     for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
                   3615:       m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
                   3616:       if ( (t = compdmm(n,m1,m2)) ||
                   3617:         (t = arf_comp(vl,C(m1),C(m2)) ) )
                   3618:         return t;
                   3619:     if ( m1 )
                   3620:       return 1;
                   3621:     else if ( m2 )
                   3622:       return -1;
                   3623:     else
                   3624:       return 0;
                   3625:   }
                   3626: }
                   3627:
1.11      noro     3628: void dpm_removecont2(DPM p1,DPM p2,DPM *r1p,DPM *r2p,Z *contp);
                   3629:
1.5       noro     3630: // p = ...+c*<<0,...0:pos>>+...
                   3631: DPM dpm_eliminate_term(DPM a,DPM p,Obj c,int pos)
                   3632: {
                   3633:   MP d0,d;
1.10      noro     3634:   DL dl;
1.5       noro     3635:   DMM m;
                   3636:   DP f;
1.11      noro     3637:   DPM a1,p1,r,r1,dmy;
                   3638:   Z dmyz;
1.5       noro     3639:
                   3640:   if ( !a ) return 0;
                   3641:   d0 = 0;
                   3642:   for ( m = BDY(a); m; m = NEXT(m) )
                   3643:     if ( m->pos == pos ) {
1.10      noro     3644:       NEXTMP(d0,d);
                   3645:       arf_chsgn(m->c,&d->c);
1.19      noro     3646:       if ( !dp_current_spec || !dp_current_spec->module_top_weight )
1.11      noro     3647:         d->dl = m->dl;
1.10      noro     3648:       else {
                   3649:         NEWDL(dl,NV(a));
                   3650:         _copydl(NV(a),m->dl,dl);
                   3651:         dl->td -= dp_current_spec->module_top_weight[pos-1];
                   3652:         d->dl = dl;
                   3653:       }
1.5       noro     3654:     }
                   3655:   if ( d0 ) {
1.6       noro     3656:     NEXT(d) = 0; MKDP(NV(a),d0,f);
                   3657:     mulcdpm(CO,c,a,&a1);
                   3658:     mulobjdpm(CO,(Obj)f,p,&p1);
                   3659:     adddpm(CO,a1,p1,&r);
1.11      noro     3660:     dpm_removecont2(0,r,&dmy,&r1,&dmyz);
                   3661:     return r1;
1.5       noro     3662:   } else
1.6       noro     3663:     return a;
1.5       noro     3664: }
                   3665:
                   3666: // <<...:i>> -> <<...:tab[i]>>
                   3667: DPM dpm_compress(DPM p,int *tab)
                   3668: {
                   3669:   DMM m,mr0,mr;
                   3670:   DPM t;
                   3671:
                   3672:   if ( !p ) return 0;
                   3673:   else {
                   3674:     for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
                   3675:       NEXTDMM(mr0,mr);
1.6       noro     3676:       mr->dl = m->dl; mr->c = m->c; mr->pos = tab[m->pos];
1.5       noro     3677:       if ( mr->pos <= 0 )
                   3678:         error("dpm_compress : invalid position");
                   3679:     }
                   3680:     NEXT(mr) = 0;
                   3681:     MKDPM(p->nv,mr0,t); t->sugar = p->sugar;
                   3682:     return t;
                   3683:   }
                   3684: }
                   3685:
                   3686: // input : s, s = syz(m) output simplified s, m
1.21      noro     3687: // assuming the term order is POT
1.10      noro     3688: void dpm_simplify_syz(LIST s,LIST m,LIST *s1,LIST *m1,LIST *w1)
1.5       noro     3689: {
1.9       noro     3690:   int lm,ls,i,j,k,pos,nv;
1.5       noro     3691:   DPM *am,*as;
                   3692:   DPM p;
                   3693:   DMM d;
                   3694:   Obj c;
1.10      noro     3695:   Z q;
                   3696:   int *tab,*dd,*new_w;
1.5       noro     3697:   NODE t,t1;
                   3698:
                   3699:   lm = length(BDY(m));
                   3700:   am = (DPM *)MALLOC((lm+1)*sizeof(DPM));
                   3701:   ls = length(BDY(s));
                   3702:   as = (DPM *)MALLOC(ls*sizeof(DPM));
                   3703:   for ( i = 1, t = BDY(m); i <= lm; i++, t = NEXT(t) ) am[i] = (DPM)BDY(t);
1.6       noro     3704:   for ( i = 0, t = BDY(s); i < ls; i++, t = NEXT(t) ) as[i] = (DPM)BDY(t);
1.5       noro     3705:
                   3706:   for ( i = 0; i < ls; i++ ) {
                   3707:     p = as[i];
1.6       noro     3708:     if ( p == 0 ) continue;
1.9       noro     3709:     nv = NV(p);
1.21      noro     3710:     for ( d = BDY(p); d; ) {
1.9       noro     3711:       dd = d->dl->d;
                   3712:       for ( k = 0; k < nv; k++ ) if ( dd[k] ) break;
                   3713:       if ( k == nv ) break;
1.21      noro     3714:       pos = d->pos;
                   3715:       while ( d && d->pos == pos ) d = NEXT(d);
1.9       noro     3716:     }
1.5       noro     3717:     if ( d ) {
                   3718:       c = d->c; pos = d->pos;
                   3719:       for ( j = 0; j < ls; j++ )
1.6       noro     3720:         if ( j != i ) {
1.5       noro     3721:           as[j] = dpm_eliminate_term(as[j],p,c,pos);
1.6       noro     3722:         }
1.5       noro     3723:       // remove m[i]
                   3724:       am[pos] = 0;
                   3725:       as[i] = 0;
                   3726:     }
                   3727:   }
                   3728:   // compress s
                   3729:   // create index table from am[]
                   3730:   // (0 0 * 0 * ...) -> (0 0 1 0 2 ... ) which means 2->1, 4->2, ...
                   3731:   tab = (int *)MALLOC((lm+1)*sizeof(int));
                   3732:   for ( j = 0, i = 1; i <= lm; i++ ) {
                   3733:     if ( am[i] ) { j++; tab[i] = j; }
1.10      noro     3734:     else { tab[i] = 0; }
1.5       noro     3735:   }
                   3736:   t = 0;
                   3737:   for ( i = ls-1; i >= 0; i-- )
                   3738:     if ( as[i] ) {
                   3739:       p = dpm_compress(as[i],tab);
                   3740:       MKNODE(t1,(pointer)p,t); t = t1;
                   3741:     }
                   3742:   MKLIST(*s1,t);
1.10      noro     3743:
1.19      noro     3744:   if ( dp_current_spec && dp_current_spec->module_top_weight ) {
1.10      noro     3745:     new_w = (int *)MALLOC(j*sizeof(int));
                   3746:     for ( j = 0, i = 1; i <= lm; i++ )
                   3747:       if ( tab[i] ) { new_w[j++] = dp_current_spec->module_top_weight[i-1]; }
                   3748:     t = 0;
                   3749:     for ( i = j-1; i >= 0; i-- ) {
                   3750:       STOZ(new_w[i],q);
                   3751:       MKNODE(t1,q,t); t = t1;
                   3752:     }
                   3753:   } else
                   3754:     t = 0;
                   3755:   MKLIST(*w1,t);
                   3756:
1.5       noro     3757:   t = 0;
                   3758:   for ( i = lm; i >= 1; i-- )
                   3759:     if ( am[i] ) {
                   3760:       MKNODE(t1,(pointer)am[i],t); t = t1;
                   3761:     }
                   3762:   MKLIST(*m1,t);
                   3763: }

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