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

Annotation of OpenXM_contrib2/asir2000/engine/dist.c, Revision 1.34

1.8       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
1.9       noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.8       noro       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.34    ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.33 2005/11/16 23:42:53 noro Exp $
1.8       noro       49: */
1.1       noro       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
1.12      noro       62: #define ORD_WEYL_ELIM 10
1.13      noro       63: #define ORD_HOMO_WW_DRL 11
1.21      noro       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();
1.1       noro       68:
                     69: int (*cmpdl)()=cmpdl_revgradlex;
                     70: int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex};
                     71:
1.2       noro       72: int do_weyl;
                     73:
1.1       noro       74: int dp_nelim,dp_fcoeffs;
1.27      noro       75: struct order_spec *dp_current_spec;
1.31      noro       76: struct modorder_spec *dp_current_modspec;
1.1       noro       77: int *dp_dl_work;
                     78:
1.24      noro       79: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr);
                     80: void comm_quod(VL vl,DP p1,DP p2,DP *pr);
                     81: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr);
                     82: void muldc_trunc(VL vl,DP p,P c,DL dl,DP *pr);
1.29      noro       83:
                     84: void order_init()
                     85: {
                     86:        struct order_spec *spec;
                     87:
                     88:        create_order_spec(0,0,&spec);
                     89:        initd(spec);
1.31      noro       90:        create_modorder_spec(0,0,&dp_current_modspec);
1.29      noro       91: }
1.24      noro       92:
1.22      noro       93: int has_sfcoef(DP f)
1.1       noro       94: {
                     95:        MP t;
                     96:
                     97:        if ( !f )
                     98:                return 0;
                     99:        for ( t = BDY(f); t; t = NEXT(t) )
1.22      noro      100:                if ( has_sfcoef_p(t->c) )
1.1       noro      101:                        break;
                    102:        return t ? 1 : 0;
                    103: }
                    104:
1.22      noro      105: int has_sfcoef_p(P f)
1.1       noro      106: {
                    107:        DCP dc;
                    108:
                    109:        if ( !f )
                    110:                return 0;
                    111:        else if ( NUM(f) )
1.22      noro      112:                return (NID((Num)f) == N_GFS) ? 1 : 0;
1.1       noro      113:        else {
                    114:                for ( dc = DC(f); dc; dc = NEXT(dc) )
1.22      noro      115:                        if ( has_sfcoef_p(COEF(dc)) )
1.1       noro      116:                                return 1;
                    117:                return 0;
                    118:        }
                    119: }
                    120:
1.19      noro      121: void initd(struct order_spec *spec)
1.1       noro      122: {
                    123:        switch ( spec->id ) {
1.28      noro      124:                case 3:
                    125:                        cmpdl = cmpdl_composite;
                    126:                        dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
                    127:                        break;
1.1       noro      128:                case 2:
                    129:                        cmpdl = cmpdl_matrix;
                    130:                        dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
                    131:                        break;
                    132:                case 1:
                    133:                        cmpdl = cmpdl_order_pair;
                    134:                        break;
                    135:                default:
                    136:                        switch ( spec->ord.simple ) {
                    137:                                case ORD_REVGRADLEX:
                    138:                                        cmpdl = cmpdl_revgradlex; break;
                    139:                                case ORD_GRADLEX:
                    140:                                        cmpdl = cmpdl_gradlex; break;
                    141:                                case ORD_BREVGRADLEX:
                    142:                                        cmpdl = cmpdl_brevgradlex; break;
                    143:                                case ORD_BGRADLEX:
                    144:                                        cmpdl = cmpdl_bgradlex; break;
                    145:                                case ORD_BLEX:
                    146:                                        cmpdl = cmpdl_blex; break;
                    147:                                case ORD_BREVREV:
                    148:                                        cmpdl = cmpdl_brevrev; break;
                    149:                                case ORD_BGRADREV:
                    150:                                        cmpdl = cmpdl_bgradrev; break;
                    151:                                case ORD_BLEXREV:
                    152:                                        cmpdl = cmpdl_blexrev; break;
                    153:                                case ORD_ELIM:
                    154:                                        cmpdl = cmpdl_elim; break;
1.12      noro      155:                                case ORD_WEYL_ELIM:
                    156:                                        cmpdl = cmpdl_weyl_elim; break;
1.13      noro      157:                                case ORD_HOMO_WW_DRL:
                    158:                                        cmpdl = cmpdl_homo_ww_drl; break;
1.21      noro      159:                                case ORD_DRL_ZIGZAG:
                    160:                                        cmpdl = cmpdl_drl_zigzag; break;
                    161:                                case ORD_HOMO_WW_DRL_ZIGZAG:
                    162:                                        cmpdl = cmpdl_homo_ww_drl_zigzag; break;
1.1       noro      163:                                case ORD_LEX: default:
                    164:                                        cmpdl = cmpdl_lex; break;
                    165:                        }
                    166:                        break;
                    167:        }
1.27      noro      168:        dp_current_spec = spec;
1.1       noro      169: }
                    170:
1.19      noro      171: void ptod(VL vl,VL dvl,P p,DP *pr)
1.1       noro      172: {
                    173:        int isconst = 0;
1.16      noro      174:        int n,i,j,k;
1.1       noro      175:        VL tvl;
                    176:        V v;
                    177:        DL d;
                    178:        MP m;
                    179:        DCP dc;
1.16      noro      180:        DCP *w;
1.1       noro      181:        DP r,s,t,u;
                    182:        P x,c;
                    183:
                    184:        if ( !p )
                    185:                *pr = 0;
                    186:        else {
                    187:                for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
                    188:                if ( NUM(p) ) {
                    189:                        NEWDL(d,n);
                    190:                        NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
                    191:                } else {
                    192:                        for ( i = 0, tvl = dvl, v = VR(p);
                    193:                                tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
                    194:                        if ( !tvl ) {
1.16      noro      195:                                for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
                    196:                                w = (DCP *)ALLOCA(k*sizeof(DCP));
                    197:                                for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
                    198:                                        w[j] = dc;
                    199:
                    200:                                for ( j = k-1, s = 0, MKV(v,x); j >= 0; j-- ) {
                    201:                                        ptod(vl,dvl,COEF(w[j]),&t); pwrp(vl,x,DEG(w[j]),&c);
1.1       noro      202:                                        muldc(vl,t,c,&r); addd(vl,r,s,&t); s = t;
                    203:                                }
                    204:                                *pr = s;
                    205:                        } else {
1.16      noro      206:                                for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
                    207:                                w = (DCP *)ALLOCA(k*sizeof(DCP));
                    208:                                for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
                    209:                                        w[j] = dc;
                    210:
                    211:                                for ( j = k-1, s = 0; j >= 0; j-- ) {
                    212:                                        ptod(vl,dvl,COEF(w[j]),&t);
1.20      noro      213:                                        NEWDL(d,n); d->d[i] = QTOS(DEG(w[j]));
                    214:                                        d->td = MUL_WEIGHT(d->d[i],i);
1.1       noro      215:                                        NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
1.2       noro      216:                                        comm_muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
1.1       noro      217:                                }
                    218:                                *pr = s;
                    219:                        }
                    220:                }
                    221:        }
1.17      noro      222: #if 0
1.22      noro      223:        if ( !dp_fcoeffs && has_sfcoef(*pr) )
                    224:                dp_fcoeffs = N_GFS;
1.17      noro      225: #endif
1.1       noro      226: }
                    227:
1.19      noro      228: void dtop(VL vl,VL dvl,DP p,P *pr)
1.1       noro      229: {
1.16      noro      230:        int n,i,j,k;
1.1       noro      231:        DL d;
                    232:        MP m;
1.16      noro      233:        MP *a;
1.1       noro      234:        P r,s,t,u,w;
                    235:        Q q;
                    236:        VL tvl;
                    237:
                    238:        if ( !p )
                    239:                *pr = 0;
                    240:        else {
1.16      noro      241:                for ( k = 0, m = BDY(p); m; m = NEXT(m), k++ );
                    242:                a = (MP *)ALLOCA(k*sizeof(MP));
                    243:                for ( j = 0, m = BDY(p); j < k; m = NEXT(m), j++ )
                    244:                        a[j] = m;
                    245:
                    246:                for ( n = p->nv, j = k-1, s = 0; j >= 0; j-- ) {
                    247:                        m = a[j];
1.1       noro      248:                        t = C(m);
                    249:                        if ( NUM(t) && NID((Num)t) == N_M ) {
                    250:                                mptop(t,&u); t = u;
                    251:                        }
                    252:                        for ( i = 0, d = m->dl, tvl = dvl;
                    253:                                i < n; tvl = NEXT(tvl), i++ ) {
                    254:                                MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
                    255:                                mulp(vl,t,u,&w); t = w;
                    256:                        }
                    257:                        addp(vl,s,t,&u); s = u;
                    258:                }
                    259:                *pr = s;
                    260:        }
                    261: }
                    262:
1.19      noro      263: void nodetod(NODE node,DP *dp)
1.1       noro      264: {
                    265:        NODE t;
                    266:        int len,i,td;
                    267:        Q e;
                    268:        DL d;
                    269:        MP m;
                    270:        DP u;
                    271:
                    272:        for ( t = node, len = 0; t; t = NEXT(t), len++ );
                    273:        NEWDL(d,len);
                    274:        for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
                    275:                e = (Q)BDY(t);
                    276:                if ( !e )
                    277:                        d->d[i] = 0;
                    278:                else if ( !NUM(e) || !RATN(e) || !INT(e) )
                    279:                        error("nodetod : invalid input");
                    280:                else {
1.20      noro      281:                        d->d[i] = QTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
1.1       noro      282:                }
                    283:        }
                    284:        d->td = td;
                    285:        NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0;
                    286:        MKDP(len,m,u); u->sugar = td; *dp = u;
                    287: }
                    288:
1.19      noro      289: int sugard(MP m)
1.1       noro      290: {
                    291:        int s;
                    292:
                    293:        for ( s = 0; m; m = NEXT(m) )
                    294:                s = MAX(s,m->dl->td);
                    295:        return s;
                    296: }
                    297:
1.19      noro      298: void addd(VL vl,DP p1,DP p2,DP *pr)
1.1       noro      299: {
                    300:        int n;
                    301:        MP m1,m2,mr,mr0;
                    302:        P t;
1.30      ohara     303:        DL d;
1.1       noro      304:
                    305:        if ( !p1 )
                    306:                *pr = p2;
                    307:        else if ( !p2 )
                    308:                *pr = p1;
                    309:        else {
1.30      ohara     310:                if ( OID(p1) <= O_R ) {
                    311:                        n = NV(p2);     NEWDL(d,n);
1.31      noro      312:                        NEWMP(m1); m1->dl = d; C(m1) = (P)p1; NEXT(m1) = 0;
1.30      ohara     313:                        MKDP(n,m1,p1); (p1)->sugar = 0;
                    314:                }
                    315:                if ( OID(p2) <= O_R ) {
                    316:                        n = NV(p1);     NEWDL(d,n);
1.31      noro      317:                        NEWMP(m2); m2->dl = d; C(m2) = (P)p2; NEXT(m2) = 0;
1.30      ohara     318:                        MKDP(n,m2,p2); (p2)->sugar = 0;
                    319:                }
1.1       noro      320:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
                    321:                        switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
                    322:                                case 0:
                    323:                                        addp(vl,C(m1),C(m2),&t);
                    324:                                        if ( t ) {
                    325:                                                NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
                    326:                                        }
                    327:                                        m1 = NEXT(m1); m2 = NEXT(m2); break;
                    328:                                case 1:
                    329:                                        NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
                    330:                                        m1 = NEXT(m1); break;
                    331:                                case -1:
                    332:                                        NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
                    333:                                        m2 = NEXT(m2); break;
                    334:                        }
                    335:                if ( !mr0 )
                    336:                        if ( m1 )
                    337:                                mr0 = m1;
                    338:                        else if ( m2 )
                    339:                                mr0 = m2;
                    340:                        else {
                    341:                                *pr = 0;
                    342:                                return;
                    343:                        }
                    344:                else if ( m1 )
                    345:                        NEXT(mr) = m1;
                    346:                else if ( m2 )
                    347:                        NEXT(mr) = m2;
                    348:                else
                    349:                        NEXT(mr) = 0;
                    350:                MKDP(NV(p1),mr0,*pr);
                    351:                if ( *pr )
                    352:                        (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                    353:        }
                    354: }
                    355:
                    356: /* for F4 symbolic reduction */
                    357:
1.19      noro      358: void symb_addd(DP p1,DP p2,DP *pr)
1.1       noro      359: {
                    360:        int n;
                    361:        MP m1,m2,mr,mr0;
                    362:
                    363:        if ( !p1 )
                    364:                *pr = p2;
                    365:        else if ( !p2 )
                    366:                *pr = p1;
                    367:        else {
                    368:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
                    369:                        NEXTMP(mr0,mr); C(mr) = (P)ONE;
                    370:                        switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
                    371:                                case 0:
                    372:                                        mr->dl = m1->dl;
                    373:                                        m1 = NEXT(m1); m2 = NEXT(m2); break;
                    374:                                case 1:
                    375:                                        mr->dl = m1->dl;
                    376:                                        m1 = NEXT(m1); break;
                    377:                                case -1:
                    378:                                        mr->dl = m2->dl;
                    379:                                        m2 = NEXT(m2); break;
                    380:                        }
                    381:                }
                    382:                if ( !mr0 )
                    383:                        if ( m1 )
                    384:                                mr0 = m1;
                    385:                        else if ( m2 )
                    386:                                mr0 = m2;
                    387:                        else {
                    388:                                *pr = 0;
                    389:                                return;
                    390:                        }
                    391:                else if ( m1 )
                    392:                        NEXT(mr) = m1;
                    393:                else if ( m2 )
                    394:                        NEXT(mr) = m2;
                    395:                else
                    396:                        NEXT(mr) = 0;
                    397:                MKDP(NV(p1),mr0,*pr);
                    398:                if ( *pr )
                    399:                        (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1.3       noro      400:        }
                    401: }
                    402:
                    403: /*
                    404:  * destructive merge of two list
                    405:  *
                    406:  * p1, p2 : list of DL
                    407:  * return : a merged list
                    408:  */
                    409:
1.19      noro      410: NODE symb_merge(NODE m1,NODE m2,int n)
1.3       noro      411: {
                    412:        NODE top,prev,cur,m,t;
1.25      noro      413:        int c,i;
                    414:        DL d1,d2;
1.3       noro      415:
                    416:        if ( !m1 )
                    417:                return m2;
                    418:        else if ( !m2 )
                    419:                return m1;
                    420:        else {
                    421:                switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
                    422:                        case 0:
                    423:                                top = m1; m = NEXT(m2);
                    424:                                break;
                    425:                        case 1:
                    426:                                top = m1; m = m2;
                    427:                                break;
                    428:                        case -1:
                    429:                                top = m2; m = m1;
                    430:                                break;
                    431:                }
                    432:                prev = top; cur = NEXT(top);
                    433:                /* BDY(prev) > BDY(m) always holds */
                    434:                while ( cur && m ) {
1.25      noro      435:                        d1 = (DL)BDY(cur);
                    436:                        d2 = (DL)BDY(m);
1.26      noro      437: #if 1
                    438:                        switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
                    439: #else
                    440:                        /* XXX only valid for DRL */
1.25      noro      441:                        if ( d1->td > d2->td )
                    442:                                c = 1;
                    443:                        else if ( d1->td < d2->td )
                    444:                                c = -1;
                    445:                        else {
                    446:                                for ( i = n-1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
                    447:                                if ( i < 0 )
                    448:                                        c = 0;
                    449:                                else if ( d1->d[i] < d2->d[i] )
                    450:                                        c = 1;
                    451:                                else
                    452:                                        c = -1;
                    453:                        }
                    454:                        switch ( c ) {
                    455: #endif
1.3       noro      456:                                case 0:
                    457:                                        m = NEXT(m);
                    458:                                        prev = cur; cur = NEXT(cur);
                    459:                                        break;
                    460:                                case 1:
                    461:                                        t = NEXT(cur); NEXT(cur) = m; m = t;
                    462:                                        prev = cur; cur = NEXT(cur);
                    463:                                        break;
                    464:                                case -1:
                    465:                                        NEXT(prev) = m; m = cur;
                    466:                                        prev = NEXT(prev); cur = NEXT(prev);
                    467:                                        break;
1.18      noro      468:                        }
                    469:                }
                    470:                if ( !cur )
                    471:                        NEXT(prev) = m;
1.23      noro      472:                return top;
                    473:        }
                    474: }
                    475:
                    476: void _adddl(int n,DL d1,DL d2,DL d3)
                    477: {
                    478:        int i;
                    479:
                    480:        d3->td = d1->td+d2->td;
                    481:        for ( i = 0; i < n; i++ )
                    482:                d3->d[i] = d1->d[i]+d2->d[i];
                    483: }
                    484:
                    485: /* m1 <- m1 U dl*f, destructive */
                    486:
                    487: NODE mul_dllist(DL dl,DP f);
                    488:
                    489: NODE symb_mul_merge(NODE m1,DL dl,DP f,int n)
                    490: {
                    491:        NODE top,prev,cur,n1;
                    492:        DP g;
                    493:        DL t,s;
                    494:        MP m;
                    495:
                    496:        if ( !m1 )
                    497:                return mul_dllist(dl,f);
                    498:        else if ( !f )
                    499:                return m1;
                    500:        else {
                    501:                m = BDY(f);
                    502:                NEWDL_NOINIT(t,n);
                    503:                _adddl(n,m->dl,dl,t);
                    504:                top = m1; prev = 0; cur = m1;
                    505:                while ( m ) {
                    506:                        switch ( (*cmpdl)(n,(DL)BDY(cur),t) ) {
                    507:                                case 0:
                    508:                                        prev = cur; cur = NEXT(cur);
                    509:                                        if ( !cur ) {
                    510:                                                MKDP(n,m,g);
                    511:                                                NEXT(prev) = mul_dllist(dl,g);
                    512:                                                return;
                    513:                                        }
                    514:                                        m = NEXT(m);
                    515:                                        if ( m ) _adddl(n,m->dl,dl,t);
                    516:                                        break;
                    517:                                case 1:
                    518:                                        prev = cur; cur = NEXT(cur);
                    519:                                        if ( !cur ) {
                    520:                                                MKDP(n,m,g);
                    521:                                                NEXT(prev) = mul_dllist(dl,g);
                    522:                                                return;
                    523:                                        }
                    524:                                        break;
                    525:                                case -1:
                    526:                                        NEWDL_NOINIT(s,n);
                    527:                                        s->td = t->td;
                    528:                                        bcopy(t->d,s->d,n*sizeof(int));
                    529:                                        NEWNODE(n1);
                    530:                                        n1->body = (pointer)s;
                    531:                                        NEXT(n1) = cur;
                    532:                                        if ( !prev ) {
                    533:                                                top = n1; cur = n1;
                    534:                                        } else {
                    535:                                                NEXT(prev) = n1; prev = n1;
                    536:                                        }
                    537:                                        m = NEXT(m);
                    538:                                        if ( m ) _adddl(n,m->dl,dl,t);
                    539:                                        break;
                    540:                        }
                    541:                }
1.18      noro      542:                return top;
                    543:        }
                    544: }
                    545:
1.19      noro      546: DLBUCKET symb_merge_bucket(DLBUCKET m1,DLBUCKET m2,int n)
1.18      noro      547: {
                    548:        DLBUCKET top,prev,cur,m,t;
                    549:
                    550:        if ( !m1 )
                    551:                return m2;
                    552:        else if ( !m2 )
                    553:                return m1;
                    554:        else {
                    555:                if ( m1->td == m2->td ) {
                    556:                        top = m1;
                    557:                        BDY(top) = symb_merge(BDY(top),BDY(m2),n);
                    558:                        m = NEXT(m2);
                    559:                } else if ( m1->td > m2->td ) {
                    560:                        top = m1; m = m2;
                    561:                } else {
                    562:                        top = m2; m = m1;
                    563:                }
                    564:                prev = top; cur = NEXT(top);
                    565:                /* prev->td > m->td always holds */
                    566:                while ( cur && m ) {
                    567:                        if ( cur->td == m->td ) {
                    568:                                BDY(cur) = symb_merge(BDY(cur),BDY(m),n);
                    569:                                m = NEXT(m);
                    570:                                prev = cur; cur = NEXT(cur);
                    571:                        } else if ( cur->td > m->td ) {
                    572:                                t = NEXT(cur); NEXT(cur) = m; m = t;
                    573:                                prev = cur; cur = NEXT(cur);
                    574:                        } else {
                    575:                                NEXT(prev) = m; m = cur;
                    576:                                prev = NEXT(prev); cur = NEXT(prev);
1.3       noro      577:                        }
                    578:                }
                    579:                if ( !cur )
                    580:                        NEXT(prev) = m;
                    581:                return top;
1.1       noro      582:        }
                    583: }
                    584:
1.19      noro      585: void subd(VL vl,DP p1,DP p2,DP *pr)
1.1       noro      586: {
                    587:        DP t;
                    588:
                    589:        if ( !p2 )
                    590:                *pr = p1;
                    591:        else {
                    592:                chsgnd(p2,&t); addd(vl,p1,t,pr);
                    593:        }
                    594: }
                    595:
1.19      noro      596: void chsgnd(DP p,DP *pr)
1.1       noro      597: {
                    598:        MP m,mr,mr0;
1.33      noro      599:        Obj r;
1.1       noro      600:
                    601:        if ( !p )
                    602:                *pr = 0;
1.33      noro      603:        else if ( OID(p) <= O_R ) {
                    604:                chsgnr((Obj)p,&r); *pr = (DP)r;
                    605:        } else {
1.1       noro      606:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    607:                        NEXTMP(mr0,mr); chsgnp(C(m),&C(mr)); mr->dl = m->dl;
                    608:                }
                    609:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    610:                if ( *pr )
                    611:                        (*pr)->sugar = p->sugar;
                    612:        }
                    613: }
                    614:
1.19      noro      615: void muld(VL vl,DP p1,DP p2,DP *pr)
1.1       noro      616: {
1.2       noro      617:        if ( ! do_weyl )
                    618:                comm_muld(vl,p1,p2,pr);
                    619:        else
                    620:                weyl_muld(vl,p1,p2,pr);
                    621: }
                    622:
1.19      noro      623: void comm_muld(VL vl,DP p1,DP p2,DP *pr)
1.2       noro      624: {
1.1       noro      625:        MP m;
                    626:        DP s,t,u;
1.5       noro      627:        int i,l,l1;
                    628:        static MP *w;
                    629:        static int wlen;
1.1       noro      630:
                    631:        if ( !p1 || !p2 )
                    632:                *pr = 0;
                    633:        else if ( OID(p1) <= O_P )
                    634:                muldc(vl,p2,(P)p1,pr);
                    635:        else if ( OID(p2) <= O_P )
                    636:                muldc(vl,p1,(P)p2,pr);
                    637:        else {
1.5       noro      638:                for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
1.4       noro      639:                for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
1.5       noro      640:                if ( l1 < l ) {
                    641:                        t = p1; p1 = p2; p2 = t;
                    642:                        l = l1;
                    643:                }
                    644:                if ( l > wlen ) {
                    645:                        if ( w ) GC_free(w);
                    646:                        w = (MP *)MALLOC(l*sizeof(MP));
                    647:                        wlen = l;
                    648:                }
1.4       noro      649:                for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
                    650:                        w[i] = m;
                    651:                for ( s = 0, i = l-1; i >= 0; i-- ) {
                    652:                        muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
1.1       noro      653:                }
1.5       noro      654:                bzero(w,l*sizeof(MP));
1.1       noro      655:                *pr = s;
                    656:        }
                    657: }
                    658:
1.24      noro      659: /* discard terms which is not a multiple of dl */
                    660:
                    661: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr)
                    662: {
                    663:        MP m;
                    664:        DP s,t,u;
                    665:        int i,l,l1;
                    666:        static MP *w;
                    667:        static int wlen;
                    668:
                    669:        if ( !p1 || !p2 )
                    670:                *pr = 0;
                    671:        else if ( OID(p1) <= O_P )
                    672:                muldc_trunc(vl,p2,(P)p1,dl,pr);
                    673:        else if ( OID(p2) <= O_P )
                    674:                muldc_trunc(vl,p1,(P)p2,dl,pr);
                    675:        else {
                    676:                for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
                    677:                for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
                    678:                if ( l1 < l ) {
                    679:                        t = p1; p1 = p2; p2 = t;
                    680:                        l = l1;
                    681:                }
                    682:                if ( l > wlen ) {
                    683:                        if ( w ) GC_free(w);
                    684:                        w = (MP *)MALLOC(l*sizeof(MP));
                    685:                        wlen = l;
                    686:                }
                    687:                for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
                    688:                        w[i] = m;
                    689:                for ( s = 0, i = l-1; i >= 0; i-- ) {
                    690:                        muldm_trunc(vl,p1,w[i],dl,&t); addd(vl,s,t,&u); s = u;
                    691:                }
                    692:                bzero(w,l*sizeof(MP));
                    693:                *pr = s;
                    694:        }
                    695: }
                    696:
                    697: void comm_quod(VL vl,DP p1,DP p2,DP *pr)
                    698: {
                    699:        MP m,m0;
                    700:        DP s,t;
                    701:        int i,n,sugar;
                    702:        DL d1,d2,d;
                    703:        Q a,b;
                    704:
                    705:        if ( !p2 )
                    706:                error("comm_quod : invalid input");
                    707:        if ( !p1 )
                    708:                *pr = 0;
                    709:        else {
                    710:                n = NV(p1);
                    711:                d2 = BDY(p2)->dl;
                    712:                m0 = 0;
                    713:                sugar = p1->sugar;
                    714:                while ( p1 ) {
                    715:                        d1 = BDY(p1)->dl;
                    716:                        NEWDL(d,n);
                    717:                        d->td = d1->td - d2->td;
                    718:                        for ( i = 0; i < n; i++ )
                    719:                                d->d[i] = d1->d[i]-d2->d[i];
                    720:                        NEXTMP(m0,m);
                    721:                        m->dl = d;
                    722:                        divq((Q)BDY(p1)->c,(Q)BDY(p2)->c,&a); chsgnq(a,&b);
                    723:                        C(m) = (P)b;
                    724:                        muldm_trunc(vl,p2,m,d2,&t);
                    725:                        addd(vl,p1,t,&s); p1 = s;
                    726:                        C(m) = (P)a;
                    727:                }
                    728:                if ( m0 ) {
                    729:                        NEXT(m) = 0; MKDP(n,m0,*pr);
                    730:                } else
                    731:                        *pr = 0;
                    732:                /* XXX */
                    733:                if ( *pr )
                    734:                        (*pr)->sugar = sugar - d2->td;
                    735:        }
                    736: }
                    737:
1.19      noro      738: void muldm(VL vl,DP p,MP m0,DP *pr)
1.1       noro      739: {
                    740:        MP m,mr,mr0;
                    741:        P c;
                    742:        DL d;
                    743:        int n;
                    744:
                    745:        if ( !p )
                    746:                *pr = 0;
                    747:        else {
                    748:                for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
                    749:                        m; m = NEXT(m) ) {
                    750:                        NEXTMP(mr0,mr);
                    751:                        if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                    752:                                mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                    753:                        else
                    754:                                mulp(vl,C(m),c,&C(mr));
                    755:                        adddl(n,m->dl,d,&mr->dl);
                    756:                }
                    757:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    758:                if ( *pr )
                    759:                        (*pr)->sugar = p->sugar + m0->dl->td;
1.2       noro      760:        }
                    761: }
                    762:
1.24      noro      763: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr)
                    764: {
                    765:        MP m,mr,mr0;
                    766:        P c;
                    767:        DL d,tdl;
                    768:        int n,i;
                    769:
                    770:        if ( !p )
                    771:                *pr = 0;
                    772:        else {
                    773:                n = NV(p);
                    774:                NEWDL(tdl,n);
                    775:                for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl;
                    776:                        m; m = NEXT(m) ) {
                    777:                        _adddl(n,m->dl,d,tdl);
                    778:                        for ( i = 0; i < n; i++ )
                    779:                                if ( tdl->d[i] < dl->d[i] )
                    780:                                        break;
                    781:                        if ( i < n )
                    782:                                continue;
                    783:                        NEXTMP(mr0,mr);
                    784:                        mr->dl = tdl;
                    785:                        NEWDL(tdl,n);
                    786:                        if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                    787:                                mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                    788:                        else
                    789:                                mulp(vl,C(m),c,&C(mr));
                    790:                }
                    791:                if ( mr0 ) {
                    792:                        NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    793:                } else
                    794:                        *pr = 0;
                    795:                if ( *pr )
                    796:                        (*pr)->sugar = p->sugar + m0->dl->td;
                    797:        }
                    798: }
                    799:
1.19      noro      800: void weyl_muld(VL vl,DP p1,DP p2,DP *pr)
1.2       noro      801: {
                    802:        MP m;
                    803:        DP s,t,u;
1.4       noro      804:        int i,l;
1.5       noro      805:        static MP *w;
                    806:        static int wlen;
1.2       noro      807:
                    808:        if ( !p1 || !p2 )
                    809:                *pr = 0;
                    810:        else if ( OID(p1) <= O_P )
                    811:                muldc(vl,p2,(P)p1,pr);
                    812:        else if ( OID(p2) <= O_P )
                    813:                muldc(vl,p1,(P)p2,pr);
                    814:        else {
1.10      noro      815:                for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
1.5       noro      816:                if ( l > wlen ) {
                    817:                        if ( w ) GC_free(w);
                    818:                        w = (MP *)MALLOC(l*sizeof(MP));
                    819:                        wlen = l;
                    820:                }
1.10      noro      821:                for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
1.4       noro      822:                        w[i] = m;
                    823:                for ( s = 0, i = l-1; i >= 0; i-- ) {
1.10      noro      824:                        weyl_muldm(vl,w[i],p2,&t); addd(vl,s,t,&u); s = u;
1.2       noro      825:                }
1.5       noro      826:                bzero(w,l*sizeof(MP));
1.2       noro      827:                *pr = s;
                    828:        }
                    829: }
                    830:
1.10      noro      831: /* monomial * polynomial */
                    832:
1.19      noro      833: void weyl_muldm(VL vl,MP m0,DP p,DP *pr)
1.2       noro      834: {
                    835:        DP r,t,t1;
                    836:        MP m;
1.10      noro      837:        DL d0;
                    838:        int n,n2,l,i,j,tlen;
                    839:        static MP *w,*psum;
                    840:        static struct cdl *tab;
1.5       noro      841:        static int wlen;
1.10      noro      842:        static int rtlen;
1.2       noro      843:
                    844:        if ( !p )
                    845:                *pr = 0;
                    846:        else {
1.4       noro      847:                for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
1.5       noro      848:                if ( l > wlen ) {
                    849:                        if ( w ) GC_free(w);
                    850:                        w = (MP *)MALLOC(l*sizeof(MP));
                    851:                        wlen = l;
                    852:                }
1.4       noro      853:                for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
                    854:                        w[i] = m;
1.10      noro      855:
                    856:                n = NV(p); n2 = n>>1;
                    857:                d0 = m0->dl;
                    858:                for ( i = 0, tlen = 1; i < n2; i++ )
                    859:                        tlen *= d0->d[n2+i]+1;
                    860:                if ( tlen > rtlen ) {
                    861:                        if ( tab ) GC_free(tab);
                    862:                        if ( psum ) GC_free(psum);
                    863:                        rtlen = tlen;
                    864:                        tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
                    865:                        psum = (MP *)MALLOC(rtlen*sizeof(MP));
                    866:                }
                    867:                bzero(psum,tlen*sizeof(MP));
                    868:                for ( i = l-1; i >= 0; i-- ) {
                    869:                        bzero(tab,tlen*sizeof(struct cdl));
                    870:                        weyl_mulmm(vl,m0,w[i],n,tab,tlen);
                    871:                        for ( j = 0; j < tlen; j++ ) {
                    872:                                if ( tab[j].c ) {
                    873:                                        NEWMP(m); m->dl = tab[j].d; C(m) = tab[j].c; NEXT(m) = psum[j];
                    874:                                        psum[j] = m;
                    875:                                }
                    876:                        }
1.2       noro      877:                }
1.10      noro      878:                for ( j = tlen-1, r = 0; j >= 0; j-- )
                    879:                        if ( psum[j] ) {
                    880:                                MKDP(n,psum[j],t); addd(vl,r,t,&t1); r = t1;
                    881:                        }
1.2       noro      882:                if ( r )
                    883:                        r->sugar = p->sugar + m0->dl->td;
                    884:                *pr = r;
                    885:        }
                    886: }
                    887:
1.10      noro      888: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
                    889: /* rtab : array of length (e0+1)*(e1+1)*... */
1.2       noro      890:
1.19      noro      891: void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
1.2       noro      892: {
1.19      noro      893:        P c,c0,c1;
1.10      noro      894:        DL d,d0,d1,dt;
                    895:        int i,j,a,b,k,l,n2,s,min,curlen;
                    896:        struct cdl *p;
                    897:        static Q *ctab;
                    898:        static struct cdl *tab;
1.5       noro      899:        static int tablen;
1.10      noro      900:        static struct cdl *tmptab;
                    901:        static int tmptablen;
1.2       noro      902:
1.10      noro      903:
                    904:        if ( !m0 || !m1 ) {
                    905:                rtab[0].c = 0;
                    906:                rtab[0].d = 0;
                    907:                return;
                    908:        }
                    909:        c0 = C(m0); c1 = C(m1);
                    910:        mulp(vl,c0,c1,&c);
                    911:        d0 = m0->dl; d1 = m1->dl;
                    912:        n2 = n>>1;
                    913:        curlen = 1;
                    914:        NEWDL(d,n);
                    915:        if ( n & 1 )
                    916:                /* offset of h-degree */
                    917:                d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
                    918:        else
                    919:                d->td = 0;
                    920:        rtab[0].c = c;
                    921:        rtab[0].d = d;
                    922:
                    923:        if ( rtablen > tmptablen ) {
                    924:                if ( tmptab ) GC_free(tmptab);
                    925:                tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
                    926:                tmptablen = rtablen;
                    927:        }
                    928:        for ( i = 0; i < n2; i++ ) {
                    929:                a = d0->d[i]; b = d1->d[n2+i];
                    930:                k = d0->d[n2+i]; l = d1->d[i];
1.20      noro      931:
                    932:                /* degree of xi^a*(Di^k*xi^l)*Di^b */
                    933:                a += l;
                    934:                b += k;
                    935:                s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
                    936:
1.10      noro      937:                if ( !k || !l ) {
                    938:                        for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
                    939:                                if ( p->c ) {
                    940:                                        dt = p->d;
                    941:                                        dt->d[i] = a;
                    942:                                        dt->d[n2+i] = b;
                    943:                                        dt->td += s;
1.5       noro      944:                                }
1.10      noro      945:                        }
                    946:                        curlen *= k+1;
                    947:                        continue;
                    948:                }
                    949:                if ( k+1 > tablen ) {
                    950:                        if ( tab ) GC_free(tab);
                    951:                        if ( ctab ) GC_free(ctab);
                    952:                        tablen = k+1;
                    953:                        tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
                    954:                        ctab = (Q *)MALLOC(tablen*sizeof(Q));
                    955:                }
                    956:                /* compute xi^a*(Di^k*xi^l)*Di^b */
                    957:                min = MIN(k,l);
                    958:                mkwc(k,l,ctab);
                    959:                bzero(tab,(k+1)*sizeof(struct cdl));
                    960:                if ( n & 1 )
                    961:                        for ( j = 0; j <= min; j++ ) {
                    962:                                NEWDL(d,n);
1.20      noro      963:                                d->d[i] = a-j; d->d[n2+i] = b-j;
1.10      noro      964:                                d->td = s;
1.20      noro      965:                                d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1.10      noro      966:                                tab[j].d = d;
                    967:                                tab[j].c = (P)ctab[j];
                    968:                        }
                    969:                else
                    970:                        for ( j = 0; j <= min; j++ ) {
                    971:                                NEWDL(d,n);
1.20      noro      972:                                d->d[i] = a-j; d->d[n2+i] = b-j;
                    973:                                d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
1.10      noro      974:                                tab[j].d = d;
                    975:                                tab[j].c = (P)ctab[j];
                    976:                        }
                    977:                bzero(ctab,(min+1)*sizeof(Q));
                    978:                comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
                    979:                curlen *= k+1;
                    980:                bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
                    981:        }
                    982: }
                    983:
                    984: /* direct product of two cdl tables
                    985:   rt[] = [
                    986:     t[0]*t1[0],...,t[n-1]*t1[0],
                    987:     t[0]*t1[1],...,t[n-1]*t1[1],
                    988:     ...
                    989:     t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
                    990:   ]
                    991: */
                    992:
1.19      noro      993: void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
1.10      noro      994: {
                    995:        int i,j;
                    996:        struct cdl *p;
                    997:        P c;
                    998:        DL d;
                    999:
                   1000:        bzero(rt,n*n1*sizeof(struct cdl));
                   1001:        for ( j = 0, p = rt; j < n1; j++ ) {
                   1002:                c = t1[j].c;
                   1003:                d = t1[j].d;
                   1004:                if ( !c )
                   1005:                        break;
                   1006:                for ( i = 0; i < n; i++, p++ ) {
                   1007:                        if ( t[i].c ) {
                   1008:                                mulp(vl,t[i].c,c,&p->c);
                   1009:                                adddl(nv,t[i].d,d,&p->d);
                   1010:                        }
1.6       noro     1011:                }
1.1       noro     1012:        }
                   1013: }
                   1014:
1.19      noro     1015: void muldc(VL vl,DP p,P c,DP *pr)
1.1       noro     1016: {
                   1017:        MP m,mr,mr0;
                   1018:
                   1019:        if ( !p || !c )
                   1020:                *pr = 0;
                   1021:        else if ( NUM(c) && UNIQ((Q)c) )
                   1022:                *pr = p;
                   1023:        else if ( NUM(c) && MUNIQ((Q)c) )
                   1024:                chsgnd(p,pr);
                   1025:        else {
                   1026:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1027:                        NEXTMP(mr0,mr);
                   1028:                        if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                   1029:                                mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                   1030:                        else
                   1031:                                mulp(vl,C(m),c,&C(mr));
                   1032:                        mr->dl = m->dl;
                   1033:                }
                   1034:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                   1035:                if ( *pr )
                   1036:                        (*pr)->sugar = p->sugar;
                   1037:        }
1.24      noro     1038: }
                   1039:
                   1040: void muldc_trunc(VL vl,DP p,P c,DL dl,DP *pr)
                   1041: {
                   1042:        MP m,mr,mr0;
                   1043:        DL mdl;
                   1044:        int i,n;
                   1045:
                   1046:        if ( !p || !c ) {
                   1047:                *pr = 0; return;
                   1048:        }
                   1049:        n = NV(p);
                   1050:        for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1051:                mdl = m->dl;
                   1052:                for ( i = 0; i < n; i++ )
                   1053:                        if ( mdl->d[i] < dl->d[i] )
                   1054:                                break;
                   1055:                if ( i < n )
                   1056:                        break;
                   1057:                NEXTMP(mr0,mr);
                   1058:                if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                   1059:                        mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                   1060:                else
                   1061:                        mulp(vl,C(m),c,&C(mr));
                   1062:                mr->dl = m->dl;
                   1063:        }
                   1064:        NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                   1065:        if ( *pr )
                   1066:                (*pr)->sugar = p->sugar;
1.1       noro     1067: }
                   1068:
1.19      noro     1069: void divsdc(VL vl,DP p,P c,DP *pr)
1.1       noro     1070: {
                   1071:        MP m,mr,mr0;
                   1072:
                   1073:        if ( !c )
                   1074:                error("disvsdc : division by 0");
                   1075:        else if ( !p )
                   1076:                *pr = 0;
                   1077:        else {
                   1078:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1079:                        NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
                   1080:                }
                   1081:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                   1082:                if ( *pr )
                   1083:                        (*pr)->sugar = p->sugar;
                   1084:        }
                   1085: }
                   1086:
1.19      noro     1087: void adddl(int n,DL d1,DL d2,DL *dr)
1.1       noro     1088: {
                   1089:        DL dt;
                   1090:        int i;
                   1091:
                   1092:        if ( !d1->td )
                   1093:                *dr = d2;
                   1094:        else if ( !d2->td )
                   1095:                *dr = d1;
                   1096:        else {
                   1097:                *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
                   1098:                dt->td = d1->td + d2->td;
                   1099:                for ( i = 0; i < n; i++ )
                   1100:                        dt->d[i] = d1->d[i]+d2->d[i];
                   1101:        }
1.11      noro     1102: }
                   1103:
                   1104: /* d1 += d2 */
                   1105:
1.19      noro     1106: void adddl_destructive(int n,DL d1,DL d2)
1.11      noro     1107: {
                   1108:        int i;
                   1109:
                   1110:        d1->td += d2->td;
                   1111:        for ( i = 0; i < n; i++ )
                   1112:                d1->d[i] += d2->d[i];
1.1       noro     1113: }
                   1114:
1.19      noro     1115: int compd(VL vl,DP p1,DP p2)
1.1       noro     1116: {
                   1117:        int n,t;
                   1118:        MP m1,m2;
                   1119:
                   1120:        if ( !p1 )
                   1121:                return p2 ? -1 : 0;
                   1122:        else if ( !p2 )
                   1123:                return 1;
                   1124:        else {
                   1125:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
                   1126:                        m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
                   1127:                        if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
                   1128:                                (t = compp(vl,C(m1),C(m2)) ) )
                   1129:                                return t;
                   1130:                if ( m1 )
                   1131:                        return 1;
                   1132:                else if ( m2 )
                   1133:                        return -1;
                   1134:                else
                   1135:                        return 0;
                   1136:        }
                   1137: }
                   1138:
1.19      noro     1139: int cmpdl_lex(int n,DL d1,DL d2)
1.1       noro     1140: {
                   1141:        int i;
                   1142:
                   1143:        for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
                   1144:        return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
                   1145: }
                   1146:
1.19      noro     1147: int cmpdl_revlex(int n,DL d1,DL d2)
1.1       noro     1148: {
                   1149:        int i;
                   1150:
                   1151:        for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
                   1152:        return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                   1153: }
                   1154:
1.19      noro     1155: int cmpdl_gradlex(int n,DL d1,DL d2)
1.1       noro     1156: {
                   1157:        if ( d1->td > d2->td )
                   1158:                return 1;
                   1159:        else if ( d1->td < d2->td )
                   1160:                return -1;
                   1161:        else
                   1162:                return cmpdl_lex(n,d1,d2);
                   1163: }
                   1164:
1.19      noro     1165: int cmpdl_revgradlex(int n,DL d1,DL d2)
1.1       noro     1166: {
1.25      noro     1167:        register int i,c;
1.7       noro     1168:        register int *p1,*p2;
                   1169:
1.1       noro     1170:        if ( d1->td > d2->td )
                   1171:                return 1;
                   1172:        else if ( d1->td < d2->td )
                   1173:                return -1;
1.7       noro     1174:        else {
1.25      noro     1175:                i = n-1;
                   1176:                p1 = d1->d+n-1;
                   1177:                p2 = d2->d+n-1;
                   1178:                while ( i >= 7 ) {
                   1179:                        c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1180:                        c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1181:                        c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1182:                        c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1183:                        c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1184:                        c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1185:                        c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1186:                        c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1187:                        i -= 8;
                   1188:                }
                   1189:                switch ( i ) {
                   1190:                        case 6:
                   1191:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1192:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1193:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1194:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1195:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1196:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1197:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1198:                                return 0;
                   1199:                        case 5:
                   1200:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1201:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1202:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1203:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1204:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1205:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1206:                                return 0;
                   1207:                        case 4:
                   1208:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1209:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1210:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1211:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1212:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1213:                                return 0;
                   1214:                        case 3:
                   1215:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1216:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1217:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1218:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1219:                                return 0;
                   1220:                        case 2:
                   1221:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1222:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1223:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1224:                                return 0;
                   1225:                        case 1:
                   1226:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1227:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1228:                                return 0;
                   1229:                        case 0:
                   1230:                                c = (*p1--) - (*p2--); if ( c ) goto LAST;
                   1231:                                return 0;
                   1232:                        default:
                   1233:                                return 0;
                   1234:                }
                   1235: LAST:
                   1236:                if ( c > 0 ) return -1;
                   1237:                else return 1;
1.7       noro     1238:        }
1.1       noro     1239: }
                   1240:
1.19      noro     1241: int cmpdl_blex(int n,DL d1,DL d2)
1.1       noro     1242: {
                   1243:        int c;
                   1244:
                   1245:        if ( c = cmpdl_lex(n-1,d1,d2) )
                   1246:                return c;
                   1247:        else {
                   1248:                c = d1->d[n-1] - d2->d[n-1];
                   1249:                return c > 0 ? 1 : c < 0 ? -1 : 0;
                   1250:        }
                   1251: }
                   1252:
1.19      noro     1253: int cmpdl_bgradlex(int n,DL d1,DL d2)
1.1       noro     1254: {
                   1255:        int e1,e2,c;
                   1256:
                   1257:        e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
                   1258:        if ( e1 > e2 )
                   1259:                return 1;
                   1260:        else if ( e1 < e2 )
                   1261:                return -1;
                   1262:        else {
                   1263:                c = cmpdl_lex(n-1,d1,d2);
                   1264:                if ( c )
                   1265:                        return c;
                   1266:                else
                   1267:                        return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
                   1268:        }
                   1269: }
                   1270:
1.19      noro     1271: int cmpdl_brevgradlex(int n,DL d1,DL d2)
1.1       noro     1272: {
                   1273:        int e1,e2,c;
                   1274:
                   1275:        e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
                   1276:        if ( e1 > e2 )
                   1277:                return 1;
                   1278:        else if ( e1 < e2 )
                   1279:                return -1;
                   1280:        else {
                   1281:                c = cmpdl_revlex(n-1,d1,d2);
                   1282:                if ( c )
                   1283:                        return c;
                   1284:                else
                   1285:                        return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
                   1286:        }
                   1287: }
                   1288:
1.19      noro     1289: int cmpdl_brevrev(int n,DL d1,DL d2)
1.1       noro     1290: {
                   1291:        int e1,e2,f1,f2,c,i;
                   1292:
                   1293:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                   1294:                e1 += d1->d[i]; e2 += d2->d[i];
                   1295:        }
                   1296:        f1 = d1->td - e1; f2 = d2->td - e2;
                   1297:        if ( e1 > e2 )
                   1298:                return 1;
                   1299:        else if ( e1 < e2 )
                   1300:                return -1;
                   1301:        else {
                   1302:                c = cmpdl_revlex(dp_nelim,d1,d2);
                   1303:                if ( c )
                   1304:                        return c;
                   1305:                else if ( f1 > f2 )
                   1306:                        return 1;
                   1307:                else if ( f1 < f2 )
                   1308:                        return -1;
                   1309:                else {
                   1310:                        for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                   1311:                        return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                   1312:                }
                   1313:        }
                   1314: }
                   1315:
1.19      noro     1316: int cmpdl_bgradrev(int n,DL d1,DL d2)
1.1       noro     1317: {
                   1318:        int e1,e2,f1,f2,c,i;
                   1319:
                   1320:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                   1321:                e1 += d1->d[i]; e2 += d2->d[i];
                   1322:        }
                   1323:        f1 = d1->td - e1; f2 = d2->td - e2;
                   1324:        if ( e1 > e2 )
                   1325:                return 1;
                   1326:        else if ( e1 < e2 )
                   1327:                return -1;
                   1328:        else {
                   1329:                c = cmpdl_lex(dp_nelim,d1,d2);
                   1330:                if ( c )
                   1331:                        return c;
                   1332:                else if ( f1 > f2 )
                   1333:                        return 1;
                   1334:                else if ( f1 < f2 )
                   1335:                        return -1;
                   1336:                else {
                   1337:                        for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                   1338:                        return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                   1339:                }
                   1340:        }
                   1341: }
                   1342:
1.19      noro     1343: int cmpdl_blexrev(int n,DL d1,DL d2)
1.1       noro     1344: {
                   1345:        int e1,e2,f1,f2,c,i;
                   1346:
                   1347:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                   1348:                e1 += d1->d[i]; e2 += d2->d[i];
                   1349:        }
                   1350:        f1 = d1->td - e1; f2 = d2->td - e2;
                   1351:        c = cmpdl_lex(dp_nelim,d1,d2);
                   1352:        if ( c )
                   1353:                return c;
                   1354:        else if ( f1 > f2 )
                   1355:                return 1;
                   1356:        else if ( f1 < f2 )
                   1357:                return -1;
                   1358:        else {
                   1359:                for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                   1360:                return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                   1361:        }
                   1362: }
                   1363:
1.19      noro     1364: int cmpdl_elim(int n,DL d1,DL d2)
1.1       noro     1365: {
                   1366:        int e1,e2,i;
                   1367:
                   1368:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                   1369:                e1 += d1->d[i]; e2 += d2->d[i];
                   1370:        }
                   1371:        if ( e1 > e2 )
                   1372:                return 1;
                   1373:        else if ( e1 < e2 )
                   1374:                return -1;
                   1375:        else
                   1376:                return cmpdl_revgradlex(n,d1,d2);
1.12      noro     1377: }
                   1378:
1.19      noro     1379: int cmpdl_weyl_elim(int n,DL d1,DL d2)
1.12      noro     1380: {
                   1381:        int e1,e2,i;
                   1382:
                   1383:        for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
                   1384:                e1 += d1->d[n-i]; e2 += d2->d[n-i];
                   1385:        }
                   1386:        if ( e1 > e2 )
                   1387:                return 1;
                   1388:        else if ( e1 < e2 )
                   1389:                return -1;
                   1390:        else if ( d1->td > d2->td )
                   1391:                return 1;
                   1392:        else if ( d1->td < d2->td )
                   1393:                return -1;
                   1394:        else return -cmpdl_revlex(n,d1,d2);
1.13      noro     1395: }
                   1396:
                   1397: /*
                   1398:        a special ordering
                   1399:        1. total order
                   1400:        2. (-w,w) for the first 2*m variables
                   1401:        3. DRL for the first 2*m variables
                   1402: */
                   1403:
1.20      noro     1404: extern int *current_weyl_weight_vector;
1.13      noro     1405:
1.19      noro     1406: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
1.13      noro     1407: {
                   1408:        int e1,e2,m,i;
                   1409:        int *p1,*p2;
                   1410:
                   1411:        if ( d1->td > d2->td )
                   1412:                return 1;
                   1413:        else if ( d1->td < d2->td )
                   1414:                return -1;
                   1415:
                   1416:        m = n>>1;
1.21      noro     1417:        for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
                   1418:                e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
                   1419:                e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1.13      noro     1420:        }
                   1421:        if ( e1 > e2 )
                   1422:                return 1;
                   1423:        else if ( e1 < e2 )
                   1424:                return -1;
                   1425:
                   1426:        e1 = d1->td - d1->d[n-1];
                   1427:        e2 = d2->td - d2->d[n-1];
                   1428:        if ( e1 > e2 )
                   1429:                return 1;
                   1430:        else if ( e1 < e2 )
                   1431:                return -1;
                   1432:
                   1433:        for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
                   1434:                i >= 0 && *p1 == *p2; i--, p1--, p2-- );
                   1435:        return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
1.21      noro     1436: }
                   1437:
                   1438: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
                   1439: {
                   1440:        int i,t,m;
                   1441:        int *p1,*p2;
                   1442:
                   1443:        if ( d1->td > d2->td )
                   1444:                return 1;
                   1445:        else if ( d1->td < d2->td )
                   1446:                return -1;
                   1447:        else {
                   1448:                m = n>>1;
                   1449:                for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
                   1450:                        if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1;
                   1451:                        if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1;
                   1452:                }
                   1453:                return 0;
                   1454:        }
                   1455: }
                   1456:
                   1457: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
                   1458: {
                   1459:        int e1,e2,m,i,t;
                   1460:        int *p1,*p2;
                   1461:
                   1462:        if ( d1->td > d2->td )
                   1463:                return 1;
                   1464:        else if ( d1->td < d2->td )
                   1465:                return -1;
                   1466:
                   1467:        m = n>>1;
                   1468:        for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
                   1469:                e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
                   1470:                e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
                   1471:        }
                   1472:        if ( e1 > e2 )
                   1473:                return 1;
                   1474:        else if ( e1 < e2 )
                   1475:                return -1;
                   1476:
                   1477:        e1 = d1->td - d1->d[n-1];
                   1478:        e2 = d2->td - d2->d[n-1];
                   1479:        if ( e1 > e2 )
                   1480:                return 1;
                   1481:        else if ( e1 < e2 )
                   1482:                return -1;
                   1483:
                   1484:        for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
                   1485:                if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1;
                   1486:                if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1;
                   1487:        }
                   1488:        return 0;
1.1       noro     1489: }
                   1490:
1.19      noro     1491: int cmpdl_order_pair(int n,DL d1,DL d2)
1.1       noro     1492: {
                   1493:        int e1,e2,i,j,l;
                   1494:        int *t1,*t2;
1.20      noro     1495:        int len,head;
1.1       noro     1496:        struct order_pair *pair;
                   1497:
1.27      noro     1498:        len = dp_current_spec->ord.block.length;
                   1499:        pair = dp_current_spec->ord.block.order_pair;
1.1       noro     1500:
1.20      noro     1501:        head = 0;
1.1       noro     1502:        for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
                   1503:                l = pair[i].length;
                   1504:                switch ( pair[i].order ) {
                   1505:                        case 0:
                   1506:                                for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20      noro     1507:                                        e1 += MUL_WEIGHT(t1[j],head+j);
                   1508:                                        e2 += MUL_WEIGHT(t2[j],head+j);
1.1       noro     1509:                                }
                   1510:                                if ( e1 > e2 )
                   1511:                                        return 1;
                   1512:                                else if ( e1 < e2 )
                   1513:                                        return -1;
                   1514:                                else {
                   1515:                                        for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
                   1516:                                        if ( j >= 0 )
                   1517:                                                return t1[j] < t2[j] ? 1 : -1;
                   1518:                                }
                   1519:                                break;
                   1520:                        case 1:
                   1521:                                for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20      noro     1522:                                        e1 += MUL_WEIGHT(t1[j],head+j);
                   1523:                                        e2 += MUL_WEIGHT(t2[j],head+j);
1.1       noro     1524:                                }
                   1525:                                if ( e1 > e2 )
                   1526:                                        return 1;
                   1527:                                else if ( e1 < e2 )
                   1528:                                        return -1;
                   1529:                                else {
                   1530:                                        for ( j = 0; j < l && t1[j] == t2[j]; j++ );
                   1531:                                        if ( j < l )
                   1532:                                                return t1[j] > t2[j] ? 1 : -1;
                   1533:                                }
                   1534:                                break;
                   1535:                        case 2:
                   1536:                                for ( j = 0; j < l && t1[j] == t2[j]; j++ );
                   1537:                                if ( j < l )
                   1538:                                        return t1[j] > t2[j] ? 1 : -1;
                   1539:                                break;
                   1540:                        default:
                   1541:                                error("cmpdl_order_pair : invalid order"); break;
                   1542:                }
1.20      noro     1543:                t1 += l; t2 += l; head += l;
1.28      noro     1544:        }
                   1545:        return 0;
                   1546: }
                   1547:
                   1548: int cmpdl_composite(int nv,DL d1,DL d2)
                   1549: {
                   1550:        int n,i,j,k,start,s,len;
                   1551:        int *dw;
                   1552:        struct sparse_weight *sw;
                   1553:        struct weight_or_block *worb;
                   1554:        int *w,*t1,*t2;
                   1555:
                   1556:        n = dp_current_spec->ord.composite.length;
                   1557:        worb = dp_current_spec->ord.composite.w_or_b;
                   1558:        w = dp_dl_work;
                   1559:        for ( i = 0, t1 = d1->d, t2 = d2->d; i < nv; i++ )
                   1560:                w[i] = t1[i]-t2[i];
                   1561:        for ( i = 0; i < n; i++, worb++ ) {
                   1562:                len = worb->length;
                   1563:                switch ( worb->type ) {
                   1564:                        case IS_DENSE_WEIGHT:
                   1565:                                dw = worb->body.dense_weight;
                   1566:                                for ( j = 0, s = 0; j < len; j++ )
                   1567:                                        s += dw[j]*w[j];
                   1568:                                if ( s > 0 ) return 1;
                   1569:                                else if ( s < 0 ) return -1;
                   1570:                                break;
                   1571:                        case IS_SPARSE_WEIGHT:
                   1572:                                sw = worb->body.sparse_weight;
                   1573:                                for ( j = 0, s = 0; j < len; j++ )
                   1574:                                        s += sw[j].value*w[sw[j].pos];
                   1575:                                if ( s > 0 ) return 1;
                   1576:                                else if ( s < 0 ) return -1;
                   1577:                                break;
                   1578:                        case IS_BLOCK:
                   1579:                                start = worb->body.block.start;
                   1580:                                switch ( worb->body.block.order ) {
                   1581:                                        case 0:
                   1582:                                                for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
                   1583:                                                        s += MUL_WEIGHT(w[k],k);
                   1584:                                                }
                   1585:                                                if ( s > 0 ) return 1;
                   1586:                                                else if ( s < 0 ) return -1;
                   1587:                                                else {
                   1588:                                                        for ( j = k-1; j >= start && w[j] == 0; j-- );
                   1589:                                                        if ( j >= start )
                   1590:                                                                return w[j] < 0 ? 1 : -1;
                   1591:                                                }
                   1592:                                                break;
                   1593:                                        case 1:
                   1594:                                                for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
                   1595:                                                        s += MUL_WEIGHT(w[k],k);
                   1596:                                                }
                   1597:                                                if ( s > 0 ) return 1;
                   1598:                                                else if ( s < 0 ) return -1;
                   1599:                                                else {
                   1600:                                                        for ( j = 0, k = start;  j < len && w[j] == 0; j++, k++ );
                   1601:                                                        if ( j < len )
                   1602:                                                                return w[j] > 0 ? 1 : -1;
                   1603:                                                }
                   1604:                                                break;
                   1605:                                        case 2:
                   1606:                                                for ( j = 0, k = start;  j < len && w[j] == 0; j++, k++ );
                   1607:                                                if ( j < len )
                   1608:                                                        return w[j] > 0 ? 1 : -1;
                   1609:                                                break;
                   1610:                                }
                   1611:                                break;
                   1612:                }
1.1       noro     1613:        }
                   1614:        return 0;
                   1615: }
                   1616:
1.19      noro     1617: int cmpdl_matrix(int n,DL d1,DL d2)
1.1       noro     1618: {
                   1619:        int *v,*w,*t1,*t2;
                   1620:        int s,i,j,len;
                   1621:        int **matrix;
                   1622:
                   1623:        for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
                   1624:                w[i] = t1[i]-t2[i];
1.27      noro     1625:        len = dp_current_spec->ord.matrix.row;
                   1626:        matrix = dp_current_spec->ord.matrix.matrix;
1.1       noro     1627:        for ( j = 0; j < len; j++ ) {
                   1628:                v = matrix[j];
                   1629:                for ( i = 0, s = 0; i < n; i++ )
                   1630:                        s += v[i]*w[i];
                   1631:                if ( s > 0 )
                   1632:                        return 1;
                   1633:                else if ( s < 0 )
                   1634:                        return -1;
                   1635:        }
                   1636:        return 0;
1.25      noro     1637: }
                   1638:
                   1639: GeoBucket create_bucket()
                   1640: {
                   1641:        GeoBucket g;
                   1642:
                   1643:        g = CALLOC(1,sizeof(struct oGeoBucket));
                   1644:        g->m = 32;
                   1645:        return g;
                   1646: }
                   1647:
                   1648: void add_bucket(GeoBucket g,NODE d,int nv)
                   1649: {
                   1650:        int l,k,m;
                   1651:
                   1652:        l = length(d);
                   1653:        for ( k = 0, m = 1; l > m; k++, m <<= 1 );
                   1654:        /* 2^(k-1) < l <= 2^k */
                   1655:        d = symb_merge(g->body[k],d,nv);
                   1656:        for ( ; length(d) > (1<<(k)); k++ ) {
                   1657:                g->body[k] = 0;
                   1658:                d = symb_merge(g->body[k+1],d,nv);
                   1659:        }
                   1660:        g->body[k] = d;
                   1661:        g->m = MAX(g->m,k);
                   1662: }
                   1663:
                   1664: DL remove_head_bucket(GeoBucket g,int nv)
                   1665: {
                   1666:        int j,i,c,m;
                   1667:        DL d;
                   1668:
                   1669:        j = -1;
                   1670:        m = g->m;
                   1671:        for ( i = 0; i <= m; i++ ) {
                   1672:                if ( !g->body[i] )
                   1673:                        continue;
                   1674:                if ( j < 0 ) j = i;
                   1675:                else {
                   1676:                        c = (*cmpdl)(nv,g->body[i]->body,g->body[j]->body);
                   1677:                        if ( c > 0 )
                   1678:                                j = i;
                   1679:                        else if ( c == 0 )
                   1680:                                g->body[i] = NEXT(g->body[i]);
                   1681:                }
                   1682:        }
                   1683:        if ( j < 0 )
                   1684:                return 0;
                   1685:        else {
                   1686:                d = g->body[j]->body;
                   1687:                g->body[j] = NEXT(g->body[j]);
                   1688:                return d;
1.31      noro     1689:        }
                   1690: }
                   1691:
                   1692: /*  DPV functions */
                   1693:
                   1694: void adddv(VL vl,DPV p1,DPV p2,DPV *pr)
                   1695: {
                   1696:        int i,len;
                   1697:        DP *e;
                   1698:
                   1699:        if ( !p1 || !p2 )
                   1700:                error("adddv : invalid argument");
                   1701:        else if ( p1->len != p2->len )
                   1702:                error("adddv : size mismatch");
                   1703:        else {
                   1704:                len = p1->len;
                   1705:                e = (DP *)MALLOC(p1->len*sizeof(DP));
                   1706:                for ( i = 0; i < len; i++ )
                   1707:                        addd(vl,p1->body[i],p2->body[i],&e[i]);
                   1708:                MKDPV(len,e,*pr);
                   1709:                (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                   1710:        }
                   1711: }
                   1712:
                   1713: void subdv(VL vl,DPV p1,DPV p2,DPV *pr)
                   1714: {
                   1715:        int i,len;
                   1716:        DP *e;
                   1717:
                   1718:        if ( !p1 || !p2 )
                   1719:                error("subdv : invalid argument");
                   1720:        else if ( p1->len != p2->len )
                   1721:                error("subdv : size mismatch");
                   1722:        else {
                   1723:                len = p1->len;
                   1724:                e = (DP *)MALLOC(p1->len*sizeof(DP));
                   1725:                for ( i = 0; i < len; i++ )
                   1726:                        subd(vl,p1->body[i],p2->body[i],&e[i]);
                   1727:                MKDPV(len,e,*pr);
                   1728:                (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                   1729:        }
                   1730: }
                   1731:
                   1732: void chsgndv(DPV p1,DPV *pr)
                   1733: {
                   1734:        int i,len;
                   1735:        DP *e;
                   1736:
                   1737:        if ( !p1 )
                   1738:                error("subdv : invalid argument");
                   1739:        else {
                   1740:                len = p1->len;
                   1741:                e = (DP *)MALLOC(p1->len*sizeof(DP));
                   1742:                for ( i = 0; i < len; i++ )
                   1743:                        chsgnd(p1->body[i],&e[i]);
                   1744:                MKDPV(len,e,*pr);
                   1745:                (*pr)->sugar = p1->sugar;
                   1746:        }
                   1747: }
                   1748:
                   1749: void muldv(VL vl,DP p1,DPV p2,DPV *pr)
                   1750: {
                   1751:        int i,len;
                   1752:        DP *e;
                   1753:
                   1754:        len = p2->len;
                   1755:        e = (DP *)MALLOC(p2->len*sizeof(DP));
                   1756:        if ( !p1 ) {
                   1757:                MKDPV(len,e,*pr);
                   1758:                (*pr)->sugar = 0;
                   1759:        } else {
                   1760:                for ( i = 0; i < len; i++ )
                   1761:                        muld(vl,p1,p2->body[i],&e[i]);
                   1762:                MKDPV(len,e,*pr);
                   1763:                (*pr)->sugar = p1->sugar + p2->sugar;
                   1764:        }
                   1765: }
                   1766:
                   1767: int compdv(VL vl,DPV p1,DPV p2)
                   1768: {
                   1769:        int i,t,len;
                   1770:
                   1771:        if ( p1->len != p2->len )
                   1772:                error("compdv : size mismatch");
                   1773:        else {
                   1774:                len = p1->len;
                   1775:                for ( i = 0; i < len; i++ )
                   1776:                        if ( t = compd(vl,p1->body[i],p2->body[i]) )
                   1777:                                return t;
                   1778:                return 0;
1.33      noro     1779:        }
                   1780: }
                   1781:
                   1782: int ni_next(int *a,int n)
                   1783: {
                   1784:        int i,j,k,kj;
                   1785:
                   1786:        /* find the first nonzero a[j] */
                   1787:        for ( j = 0; a[j] == 0; j++ );
                   1788:        /* find the first zero a[k] after a[j] */
                   1789:        for ( k = j; k < n && a[k] == 1; k++ );
                   1790:        if ( k == n ) return 0;
                   1791:        /* a[0] = 0, ... , a[j-1] = 0, a[j] = 1, ..., a[k-1] = 1, a[k] = 0 */
                   1792:        /* a[0] = 1,..., a[k-j-2] = 1, a[k-j-1] = 0, ..., a[k-1] = 0, a[k] = 1 */
                   1793:        kj = k-j-1;
                   1794:        for ( i = 0; i < kj; i++ ) a[i] = 1;
                   1795:        for ( ; i < k; i++ ) a[i] = 0;
                   1796:        a[k] = 1;
                   1797:        return 1;
                   1798: }
                   1799:
                   1800: int comp_nbm(NBM a,NBM b)
                   1801: {
                   1802:        int d,i,w;
                   1803:        int *ab,*bb;
                   1804:
                   1805:        if ( a->d > b->d ) return 1;
                   1806:        else if ( a->d < b->d ) return -1;
                   1807:        else {
                   1808:                d = a->d; ab = a->b; bb = b->b;
                   1809:                w = (d+31)/32;
                   1810:                for ( i = 0; i < w; i++ )
                   1811:                        if ( ab[i] > bb[i] ) return 1;
                   1812:                        else if ( ab[i] < bb[i] ) return -1;
                   1813:                return 0;
                   1814:        }
                   1815: }
                   1816:
                   1817: NBM mul_nbm(NBM a,NBM b)
                   1818: {
                   1819:        int ad,bd,d,i,j;
                   1820:        int *ab,*bb,*mb;
                   1821:        NBM m;
                   1822:        Q c,c1;
                   1823:        NODE r;
                   1824:        NBP u;
                   1825:
                   1826:        ad = a->d; bd = b->d; ab = a->b; bb = b->b;
                   1827:        d = ad + bd;
                   1828:        NEWNBM(m); NEWNBMBDY(m,d);
                   1829:        m->d = d; mulq(a->c,b->c,&m->c); mb = m->b;
                   1830:        j = 0;
                   1831:        for ( i = 0; i < ad; i++, j++ )
                   1832:                if ( NBM_GET(ab,i) ) NBM_SET(mb,j);
                   1833:                else NBM_CLR(mb,j);
                   1834:        for ( i = 0; i < bd; i++, j++ )
                   1835:                if ( NBM_GET(bb,i) ) NBM_SET(mb,j);
                   1836:                else NBM_CLR(mb,j);
                   1837:        return m;
                   1838: }
                   1839:
                   1840: NBP shuffle_mul_nbm(NBM a,NBM b)
                   1841: {
                   1842:        int ad,bd,d,i,ai,bi,bit,s;
                   1843:        int *ab,*bb,*wmb,*w;
                   1844:        NBM wm,tm;
                   1845:        Q c,c1;
                   1846:        NODE r,t,t1,p;
                   1847:        NBP u;
                   1848:
                   1849:        ad = a->d; bd = b->d; ab = a->b; bb = b->b;
                   1850:        d = ad + bd;
                   1851:        w = (int *)ALLOCA(d*sizeof(int));
                   1852:        NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   1853:        for ( i = 0; i < ad; i++ ) w[i] = 1;
                   1854:        for ( ; i < d; i++ ) w[i] = 0;
                   1855:        mulq(a->c,b->c,&c);
                   1856:        r = 0;
                   1857:        do {
                   1858:                wm->d = d; wm->c = c;
                   1859:                ai = 0; bi = 0;
                   1860:                for ( i = 0; i < d; i++ ) {
                   1861:                        if ( w[i] ) { bit = NBM_GET(ab,ai); ai++; }
                   1862:                        else { bit = NBM_GET(bb,bi); bi++; }
                   1863:                        if ( bit ) NBM_SET(wmb,i);
                   1864:                        else NBM_CLR(wmb,i);
                   1865:                }
                   1866:                for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
                   1867:                        tm = (NBM)BDY(t);
                   1868:                        s = comp_nbm(tm,wm);
                   1869:                        if ( s < 0 ) {
                   1870:                                /* insert */
                   1871:                                MKNODE(t1,wm,t);
                   1872:                                if ( !p ) r = t1;
                   1873:                                else NEXT(p) = t1;
                   1874:                                NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   1875:                                break;
                   1876:                        } else if ( s == 0 ) {
                   1877:                                /* add coefs */
                   1878:                                addq(tm->c,c,&c1);
                   1879:                                if ( c1 ) tm->c = c1;
                   1880:                                else NEXT(p) = NEXT(t);
                   1881:                                break;
                   1882:                        }
                   1883:                }
                   1884:                if ( !t ) {
                   1885:                        /* append */
                   1886:                        MKNODE(t1,wm,t);
                   1887:                        if ( !p ) r = t1;
                   1888:                        else NEXT(p) = t1;
                   1889:                        NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                   1890:                }
                   1891:        } while ( ni_next(w,d) );
                   1892:        MKNBP(u,r);
                   1893:        return u;
                   1894: }
                   1895:
1.34    ! noro     1896: int nbmtoxky(NBM a,int *b)
        !          1897: {
        !          1898:        int d,i,j,k;
        !          1899:        int *p;
        !          1900:
        !          1901:        d = a->d; p = a->b;
        !          1902:        for ( i = j = 0, k = 1; i < d; i++ ) {
        !          1903:                if ( !NBM_GET(p,i) ) {
        !          1904:                        b[j++] = k;
        !          1905:                        k = 1;
        !          1906:                        i++;
        !          1907:                } else k++;
        !          1908:        }
        !          1909:        return j;
        !          1910: }
        !          1911:
        !          1912: NBP harmonic_mul_nbm(NBM a,NBM b)
        !          1913: {
        !          1914:        int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1;
        !          1915:        int i,j,k,ia,ib,s;
        !          1916:        int *wa,*wb,*w,*wab,*wa1,*wmb;
        !          1917:        Q c,c1;
        !          1918:        NBM wm,tm;
        !          1919:        NODE r,t1,t,p;
        !          1920:        NBP u;
        !          1921:
        !          1922:        da = a->d; db = b->d; d = da+db;
        !          1923:        wa = (int *)ALLOCA(da*sizeof(int));
        !          1924:        wb = (int *)ALLOCA(db*sizeof(int));
        !          1925:        la = nbmtoxky(a,wa);
        !          1926:        lb = nbmtoxky(b,wb);
        !          1927:        mulq(a->c,b->c,&c);
        !          1928:        /* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */
        !          1929:        /* lmax : total length */
        !          1930:        lmax = la+lb;
        !          1931:        lmin = la>lb?la:lb;
        !          1932:        w = (int *)ALLOCA(lmax*sizeof(int));
        !          1933:        /* position of a+b */
        !          1934:        wab = (int *)ALLOCA(lmax*sizeof(int));
        !          1935:        /* position of a */
        !          1936:        wa1 = (int *)ALLOCA(lmax*sizeof(int));
        !          1937:        NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
        !          1938:        for ( l = lmin, r = 0; l <= lmax; l++ ) {
        !          1939:                lab = lmax - l;
        !          1940:                la1 = la - lab;
        !          1941:                lb1 = lb - lab;
        !          1942:                lab1 = l-lab;
        !          1943:                /* partion l into three parts: a, b, a+b */
        !          1944:                /* initialize wab */
        !          1945:                for ( i = 0; i < lab; i++ ) wab[i] = 1;
        !          1946:                for ( ; i < l; i++ ) wab[i] = 0;
        !          1947:                do {
        !          1948:                        /* initialize wa1 */
        !          1949:                        for ( i = 0; i < la1; i++ ) wa1[i] = 1;
        !          1950:                        for ( ; i < lab1; i++ ) wa1[i] = 0;
        !          1951:                        do {
        !          1952:                                ia = 0; ib = 0;
        !          1953:                                for ( i = j = 0; i < l; i++ )
        !          1954:                                        if ( wab[i] ) w[i] = wa[ia++]+wb[ib++];
        !          1955:                                        else if ( wa1[j++] ) w[i] = wa[ia++];
        !          1956:                                        else w[i] = wb[ib++];
        !          1957:                                for ( i = j = 0; i < l; i++ ) {
        !          1958:                                        for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j);
        !          1959:                                        NBM_CLR(wmb,j); j++;
        !          1960:                                }
        !          1961:                                wm->d = j; wm->c = c;
        !          1962:                                for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
        !          1963:                                        tm = (NBM)BDY(t);
        !          1964:                                        s = comp_nbm(tm,wm);
        !          1965:                                        if ( s < 0 ) {
        !          1966:                                                /* insert */
        !          1967:                                                MKNODE(t1,wm,t);
        !          1968:                                                if ( !p ) r = t1;
        !          1969:                                                else NEXT(p) = t1;
        !          1970:                                                NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
        !          1971:                                                break;
        !          1972:                                        } else if ( s == 0 ) {
        !          1973:                                                /* add coefs */
        !          1974:                                                addq(tm->c,c,&c1);
        !          1975:                                                if ( c1 ) tm->c = c1;
        !          1976:                                                else NEXT(p) = NEXT(t);
        !          1977:                                                break;
        !          1978:                                        }
        !          1979:                                }
        !          1980:                                if ( !t ) {
        !          1981:                                        /* append */
        !          1982:                                        MKNODE(t1,wm,t);
        !          1983:                                        if ( !p ) r = t1;
        !          1984:                                        else NEXT(p) = t1;
        !          1985:                                        NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
        !          1986:                                }
        !          1987:                        } while ( ni_next(wa1,lab1) );
        !          1988:                } while ( ni_next(wab,l) );
        !          1989:        }
        !          1990:        MKNBP(u,r);
        !          1991:        return u;
        !          1992: }
        !          1993:
1.33      noro     1994: void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   1995: {
                   1996:        NODE b1,b2,br,br0;
                   1997:        NBM m1,m2,m;
                   1998:        Q c;
                   1999:
                   2000:        if ( !p1 )
                   2001:                *rp = p2;
                   2002:        else if ( !p2 )
                   2003:                *rp = p1;
                   2004:        else {
                   2005:                for ( b1 = BDY(p1), b2 = BDY(p2), br0 = 0; b1 && b2; ) {
                   2006:                        m1 = (NBM)BDY(b1); m2 = (NBM)BDY(b2);
                   2007:                        switch ( comp_nbm(m1,m2) ) {
                   2008:                                case 0:
                   2009:                                        addq(m1->c,m2->c,&c);
                   2010:                                        if ( c ) {
                   2011:                                                NEXTNODE(br0,br);
                   2012:                                                NEWNBM(m); m->d = m1->d; m->c = c; m->b = m1->b;
                   2013:                                                BDY(br) = (pointer)m;
                   2014:                                        }
                   2015:                                        b1 = NEXT(b1); b2 = NEXT(b2); break;
                   2016:                                case 1:
                   2017:                                        NEXTNODE(br0,br); BDY(br) = BDY(b1);
                   2018:                                        b1 = NEXT(b1); break;
                   2019:                                case -1:
                   2020:                                        NEXTNODE(br0,br); BDY(br) = BDY(b2);
                   2021:                                        b2 = NEXT(b2); break;
                   2022:                        }
1.34    ! noro     2023:                }
        !          2024:                if ( !br0 )
        !          2025:                        if ( b1 )
        !          2026:                                br0 = b1;
1.33      noro     2027:                        else if ( b2 )
1.34    ! noro     2028:                                br0 = b2;
        !          2029:                        else {
        !          2030:                                *rp = 0;
        !          2031:                                return;
        !          2032:                        }
        !          2033:                else if ( b1 )
        !          2034:                        NEXT(br) = b1;
        !          2035:                else if ( b2 )
1.33      noro     2036:                                NEXT(br) = b2;
1.34    ! noro     2037:                else
        !          2038:                        NEXT(br) = 0;
        !          2039:                MKNBP(*rp,br0);
1.33      noro     2040:        }
                   2041: }
                   2042:
                   2043: void subnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2044: {
                   2045:        NBP t;
                   2046:
                   2047:        chsgnnbp(p2,&t);
                   2048:        addnbp(vl,p1,t,rp);
                   2049: }
                   2050:
                   2051: void chsgnnbp(NBP p,NBP *rp)
                   2052: {
                   2053:        NODE r0,r,b;
                   2054:        NBM m,m1;
                   2055:
                   2056:        for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
                   2057:                NEXTNODE(r0,r);
                   2058:                m = (NBM)BDY(b);
                   2059:                NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnq(m->c,&m1->c);
1.34    ! noro     2060:                BDY(r) = m1;
1.33      noro     2061:        }
                   2062:        if ( r0 ) NEXT(r) = 0;
                   2063:        MKNBP(*rp,r0);
                   2064: }
                   2065:
                   2066: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp);
                   2067: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
                   2068:
                   2069: void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2070: {
                   2071:        NODE b;
                   2072:        NBP r,t,s;
                   2073:
                   2074:        if ( !p1 || !p2 ) *rp = 0;
                   2075:        else if ( length(BDY(p1)) < length(BDY(p2)) ) {
                   2076:                for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {
                   2077:                        mulnbmnbp(vl,(NBM)BDY(b),p2,&t);
                   2078:                        addnbp(vl,r,t,&s); r = s;
                   2079:                }
                   2080:                *rp = r;
                   2081:        } else {
                   2082:                for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) {
                   2083:                        mulnbpnbm(vl,p1,(NBM)BDY(b),&t);
                   2084:                        addnbp(vl,r,t,&s); r = s;
                   2085:                }
                   2086:                *rp = r;
                   2087:        }
                   2088: }
                   2089:
                   2090: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp)
                   2091: {
                   2092:        NODE b,r0,r;
                   2093:
                   2094:        if ( !p ) *rp = 0;
                   2095:        else {
                   2096:                for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
                   2097:                        NEXTNODE(r0,r);
                   2098:                        BDY(r) = mul_nbm(m,(NBM)BDY(b));
                   2099:                }
                   2100:                if ( r0 ) NEXT(r) = 0;
                   2101:                MKNBP(*rp,r0);
                   2102:        }
                   2103: }
                   2104:
                   2105: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp)
                   2106: {
                   2107:        NODE b,r0,r;
                   2108:
                   2109:        if ( !p ) *rp = 0;
                   2110:        else {
                   2111:                for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
                   2112:                        NEXTNODE(r0,r);
                   2113:                        BDY(r) = mul_nbm((NBM)BDY(b),m);
                   2114:                }
                   2115:                if ( r0 ) NEXT(r) = 0;
                   2116:                MKNBP(*rp,r0);
                   2117:        }
                   2118: }
                   2119:
                   2120: void pwrnbp(VL vl,NBP a,Q q,NBP *c)
                   2121: {
                   2122:        int t;
                   2123:        NBP a1,a2;
                   2124:        N n1;
                   2125:        Q q1;
                   2126:        NBM m;
                   2127:        NODE r;
                   2128:
                   2129:        if ( !q ) {
                   2130:                 NEWNBM(m); m->d = 0; m->c = ONE; m->b = 0;
                   2131:                 MKNODE(r,m,0); MKNBP(*c,r);
                   2132:        } else if ( !a )
                   2133:                *c = 0;
                   2134:        else if ( UNIQ(q) )
                   2135:                *c = a;
                   2136:        else {
                   2137:                t = divin(NM(q),2,&n1); NTOQ(n1,1,q1);
                   2138:                pwrnbp(vl,a,q1,&a1);
                   2139:                mulnbp(vl,a1,a1,&a2);
                   2140:                if ( t )
                   2141:                        mulnbp(vl,a2,a,c);
                   2142:                else
                   2143:                        *c = a2;
                   2144:        }
                   2145: }
                   2146:
                   2147: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                   2148: {
1.34    ! noro     2149:        NODE b1,b2;
1.33      noro     2150:        NBP r,t,s;
1.34    ! noro     2151:        NBM m;
1.33      noro     2152:
                   2153:        if ( !p1 || !p2 ) *rp = 0;
                   2154:        else {
1.34    ! noro     2155:                for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
        !          2156:                        for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
        !          2157:                                t = shuffle_mul_nbm(m,(NBM)BDY(b2));
        !          2158:                                addnbp(vl,r,t,&s); r = s;
        !          2159:                        }
1.33      noro     2160:                *rp = r;
                   2161:        }
                   2162: }
                   2163:
1.34    ! noro     2164: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
1.33      noro     2165: {
1.34    ! noro     2166:        NODE b1,b2;
        !          2167:        NBP r,t,s;
        !          2168:        NBM m;
1.33      noro     2169:
1.34    ! noro     2170:        if ( !p1 || !p2 ) *rp = 0;
1.33      noro     2171:        else {
1.34    ! noro     2172:                for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
        !          2173:                        for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
        !          2174:                                t = harmonic_mul_nbm(m,(NBM)BDY(b2));
        !          2175:                                addnbp(vl,r,t,&s); r = s;
        !          2176:                        }
1.33      noro     2177:                *rp = r;
1.25      noro     2178:        }
1.1       noro     2179: }

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