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

1.3     ! noro        1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.2 2000/04/05 08:32:17 noro Exp $ */
1.1       noro        2: #include "ca.h"
                      3:
                      4: #define NV(p) ((p)->nv)
                      5: #define C(p) ((p)->c)
                      6:
                      7: #define ORD_REVGRADLEX 0
                      8: #define ORD_GRADLEX 1
                      9: #define ORD_LEX 2
                     10: #define ORD_BREVGRADLEX 3
                     11: #define ORD_BGRADLEX 4
                     12: #define ORD_BLEX 5
                     13: #define ORD_BREVREV 6
                     14: #define ORD_BGRADREV 7
                     15: #define ORD_BLEXREV 8
                     16: #define ORD_ELIM 9
                     17:
                     18: int (*cmpdl)()=cmpdl_revgradlex;
                     19: int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex};
                     20:
1.2       noro       21: void comm_muld(VL,DP,DP,DP *);
                     22: void weyl_muld(VL,DP,DP,DP *);
                     23: void weyl_muldm(VL,DP,MP,DP *);
                     24: void weyl_mulmm(VL,MP,MP,int,DP *);
                     25: void mkwc(int,int,Q *);
                     26:
                     27: int do_weyl;
                     28:
1.1       noro       29: int dp_nelim,dp_fcoeffs;
                     30: struct order_spec dp_current_spec;
                     31: int *dp_dl_work;
                     32:
                     33: int has_fcoef(DP);
                     34: int has_fcoef_p(P);
                     35:
                     36: int has_fcoef(f)
                     37: DP f;
                     38: {
                     39:        MP t;
                     40:
                     41:        if ( !f )
                     42:                return 0;
                     43:        for ( t = BDY(f); t; t = NEXT(t) )
                     44:                if ( has_fcoef_p(t->c) )
                     45:                        break;
                     46:        return t ? 1 : 0;
                     47: }
                     48:
                     49: int has_fcoef_p(f)
                     50: P f;
                     51: {
                     52:        DCP dc;
                     53:
                     54:        if ( !f )
                     55:                return 0;
                     56:        else if ( NUM(f) )
                     57:                return (NID((Num)f) == N_LM || NID((Num)f) == N_GF2N) ? 1 : 0;
                     58:        else {
                     59:                for ( dc = DC(f); dc; dc = NEXT(dc) )
                     60:                        if ( has_fcoef_p(COEF(dc)) )
                     61:                                return 1;
                     62:                return 0;
                     63:        }
                     64: }
                     65:
                     66: void initd(spec)
                     67: struct order_spec *spec;
                     68: {
                     69:        switch ( spec->id ) {
                     70:                case 2:
                     71:                        cmpdl = cmpdl_matrix;
                     72:                        dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
                     73:                        break;
                     74:                case 1:
                     75:                        cmpdl = cmpdl_order_pair;
                     76:                        break;
                     77:                default:
                     78:                        switch ( spec->ord.simple ) {
                     79:                                case ORD_REVGRADLEX:
                     80:                                        cmpdl = cmpdl_revgradlex; break;
                     81:                                case ORD_GRADLEX:
                     82:                                        cmpdl = cmpdl_gradlex; break;
                     83:                                case ORD_BREVGRADLEX:
                     84:                                        cmpdl = cmpdl_brevgradlex; break;
                     85:                                case ORD_BGRADLEX:
                     86:                                        cmpdl = cmpdl_bgradlex; break;
                     87:                                case ORD_BLEX:
                     88:                                        cmpdl = cmpdl_blex; break;
                     89:                                case ORD_BREVREV:
                     90:                                        cmpdl = cmpdl_brevrev; break;
                     91:                                case ORD_BGRADREV:
                     92:                                        cmpdl = cmpdl_bgradrev; break;
                     93:                                case ORD_BLEXREV:
                     94:                                        cmpdl = cmpdl_blexrev; break;
                     95:                                case ORD_ELIM:
                     96:                                        cmpdl = cmpdl_elim; break;
                     97:                                case ORD_LEX: default:
                     98:                                        cmpdl = cmpdl_lex; break;
                     99:                        }
                    100:                        break;
                    101:        }
                    102:        dp_current_spec = *spec;
                    103: }
                    104:
                    105: void ptod(vl,dvl,p,pr)
                    106: VL vl,dvl;
                    107: P p;
                    108: DP *pr;
                    109: {
                    110:        int isconst = 0;
                    111:        int n,i;
                    112:        VL tvl;
                    113:        V v;
                    114:        DL d;
                    115:        MP m;
                    116:        DCP dc;
                    117:        DP r,s,t,u;
                    118:        P x,c;
                    119:
                    120:        if ( !p )
                    121:                *pr = 0;
                    122:        else {
                    123:                for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
                    124:                if ( NUM(p) ) {
                    125:                        NEWDL(d,n);
                    126:                        NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
                    127:                } else {
                    128:                        for ( i = 0, tvl = dvl, v = VR(p);
                    129:                                tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
                    130:                        if ( !tvl ) {
                    131:                                for ( dc = DC(p), s = 0, MKV(v,x); dc; dc = NEXT(dc) ) {
                    132:                                        ptod(vl,dvl,COEF(dc),&t); pwrp(vl,x,DEG(dc),&c);
                    133:                                        muldc(vl,t,c,&r); addd(vl,r,s,&t); s = t;
                    134:                                }
                    135:                                *pr = s;
                    136:                        } else {
                    137:                                for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
                    138:                                        ptod(vl,dvl,COEF(dc),&t);
                    139:                                        NEWDL(d,n); d->td = QTOS(DEG(dc)); d->d[i] = d->td;
                    140:                                        NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
1.2       noro      141:                                        comm_muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
1.1       noro      142:                                }
                    143:                                *pr = s;
                    144:                        }
                    145:                }
                    146:        }
                    147:        if ( !dp_fcoeffs && has_fcoef(*pr) )
                    148:                dp_fcoeffs = 1;
                    149: }
                    150:
                    151: void dtop(vl,dvl,p,pr)
                    152: VL vl,dvl;
                    153: DP p;
                    154: P *pr;
                    155: {
                    156:        int n,i;
                    157:        DL d;
                    158:        MP m;
                    159:        P r,s,t,u,w;
                    160:        Q q;
                    161:        VL tvl;
                    162:
                    163:        if ( !p )
                    164:                *pr = 0;
                    165:        else {
                    166:                for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
                    167:                        t = C(m);
                    168:                        if ( NUM(t) && NID((Num)t) == N_M ) {
                    169:                                mptop(t,&u); t = u;
                    170:                        }
                    171:                        for ( i = 0, d = m->dl, tvl = dvl;
                    172:                                i < n; tvl = NEXT(tvl), i++ ) {
                    173:                                MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
                    174:                                mulp(vl,t,u,&w); t = w;
                    175:                        }
                    176:                        addp(vl,s,t,&u); s = u;
                    177:                }
                    178:                *pr = s;
                    179:        }
                    180: }
                    181:
                    182: void nodetod(node,dp)
                    183: NODE node;
                    184: DP *dp;
                    185: {
                    186:        NODE t;
                    187:        int len,i,td;
                    188:        Q e;
                    189:        DL d;
                    190:        MP m;
                    191:        DP u;
                    192:
                    193:        for ( t = node, len = 0; t; t = NEXT(t), len++ );
                    194:        NEWDL(d,len);
                    195:        for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
                    196:                e = (Q)BDY(t);
                    197:                if ( !e )
                    198:                        d->d[i] = 0;
                    199:                else if ( !NUM(e) || !RATN(e) || !INT(e) )
                    200:                        error("nodetod : invalid input");
                    201:                else {
                    202:                        d->d[i] = QTOS((Q)e); td += d->d[i];
                    203:                }
                    204:        }
                    205:        d->td = td;
                    206:        NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0;
                    207:        MKDP(len,m,u); u->sugar = td; *dp = u;
                    208: }
                    209:
                    210: int sugard(m)
                    211: MP m;
                    212: {
                    213:        int s;
                    214:
                    215:        for ( s = 0; m; m = NEXT(m) )
                    216:                s = MAX(s,m->dl->td);
                    217:        return s;
                    218: }
                    219:
                    220: void addd(vl,p1,p2,pr)
                    221: VL vl;
                    222: DP p1,p2,*pr;
                    223: {
                    224:        int n;
                    225:        MP m1,m2,mr,mr0;
                    226:        P t;
                    227:
                    228:        if ( !p1 )
                    229:                *pr = p2;
                    230:        else if ( !p2 )
                    231:                *pr = p1;
                    232:        else {
                    233:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
                    234:                        switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
                    235:                                case 0:
                    236:                                        addp(vl,C(m1),C(m2),&t);
                    237:                                        if ( t ) {
                    238:                                                NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
                    239:                                        }
                    240:                                        m1 = NEXT(m1); m2 = NEXT(m2); break;
                    241:                                case 1:
                    242:                                        NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
                    243:                                        m1 = NEXT(m1); break;
                    244:                                case -1:
                    245:                                        NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
                    246:                                        m2 = NEXT(m2); break;
                    247:                        }
                    248:                if ( !mr0 )
                    249:                        if ( m1 )
                    250:                                mr0 = m1;
                    251:                        else if ( m2 )
                    252:                                mr0 = m2;
                    253:                        else {
                    254:                                *pr = 0;
                    255:                                return;
                    256:                        }
                    257:                else if ( m1 )
                    258:                        NEXT(mr) = m1;
                    259:                else if ( m2 )
                    260:                        NEXT(mr) = m2;
                    261:                else
                    262:                        NEXT(mr) = 0;
                    263:                MKDP(NV(p1),mr0,*pr);
                    264:                if ( *pr )
                    265:                        (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                    266:        }
                    267: }
                    268:
                    269: /* for F4 symbolic reduction */
                    270:
                    271: void symb_addd(p1,p2,pr)
                    272: DP p1,p2,*pr;
                    273: {
                    274:        int n;
                    275:        MP m1,m2,mr,mr0;
                    276:        P t;
                    277:
                    278:        if ( !p1 )
                    279:                *pr = p2;
                    280:        else if ( !p2 )
                    281:                *pr = p1;
                    282:        else {
                    283:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
                    284:                        NEXTMP(mr0,mr); C(mr) = (P)ONE;
                    285:                        switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
                    286:                                case 0:
                    287:                                        mr->dl = m1->dl;
                    288:                                        m1 = NEXT(m1); m2 = NEXT(m2); break;
                    289:                                case 1:
                    290:                                        mr->dl = m1->dl;
                    291:                                        m1 = NEXT(m1); break;
                    292:                                case -1:
                    293:                                        mr->dl = m2->dl;
                    294:                                        m2 = NEXT(m2); break;
                    295:                        }
                    296:                }
                    297:                if ( !mr0 )
                    298:                        if ( m1 )
                    299:                                mr0 = m1;
                    300:                        else if ( m2 )
                    301:                                mr0 = m2;
                    302:                        else {
                    303:                                *pr = 0;
                    304:                                return;
                    305:                        }
                    306:                else if ( m1 )
                    307:                        NEXT(mr) = m1;
                    308:                else if ( m2 )
                    309:                        NEXT(mr) = m2;
                    310:                else
                    311:                        NEXT(mr) = 0;
                    312:                MKDP(NV(p1),mr0,*pr);
                    313:                if ( *pr )
                    314:                        (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1.3     ! noro      315:        }
        !           316: }
        !           317:
        !           318: /*
        !           319:  * destructive merge of two list
        !           320:  *
        !           321:  * p1, p2 : list of DL
        !           322:  * return : a merged list
        !           323:  */
        !           324:
        !           325: NODE symb_merge(m1,m2,n)
        !           326: NODE m1,m2;
        !           327: int n;
        !           328: {
        !           329:        NODE top,prev,cur,m,t;
        !           330:
        !           331:        if ( !m1 )
        !           332:                return m2;
        !           333:        else if ( !m2 )
        !           334:                return m1;
        !           335:        else {
        !           336:                switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
        !           337:                        case 0:
        !           338:                                top = m1; m = NEXT(m2);
        !           339:                                break;
        !           340:                        case 1:
        !           341:                                top = m1; m = m2;
        !           342:                                break;
        !           343:                        case -1:
        !           344:                                top = m2; m = m1;
        !           345:                                break;
        !           346:                }
        !           347:                prev = top; cur = NEXT(top);
        !           348:                /* BDY(prev) > BDY(m) always holds */
        !           349:                while ( cur && m ) {
        !           350:                        switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
        !           351:                                case 0:
        !           352:                                        m = NEXT(m);
        !           353:                                        prev = cur; cur = NEXT(cur);
        !           354:                                        break;
        !           355:                                case 1:
        !           356:                                        t = NEXT(cur); NEXT(cur) = m; m = t;
        !           357:                                        prev = cur; cur = NEXT(cur);
        !           358:                                        break;
        !           359:                                case -1:
        !           360:                                        NEXT(prev) = m; m = cur;
        !           361:                                        prev = NEXT(prev); cur = NEXT(prev);
        !           362:                                        break;
        !           363:                        }
        !           364:                }
        !           365:                if ( !cur )
        !           366:                        NEXT(prev) = m;
        !           367:                return top;
1.1       noro      368:        }
                    369: }
                    370:
                    371: void subd(vl,p1,p2,pr)
                    372: VL vl;
                    373: DP p1,p2,*pr;
                    374: {
                    375:        DP t;
                    376:
                    377:        if ( !p2 )
                    378:                *pr = p1;
                    379:        else {
                    380:                chsgnd(p2,&t); addd(vl,p1,t,pr);
                    381:        }
                    382: }
                    383:
                    384: void chsgnd(p,pr)
                    385: DP p,*pr;
                    386: {
                    387:        MP m,mr,mr0;
                    388:
                    389:        if ( !p )
                    390:                *pr = 0;
                    391:        else {
                    392:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    393:                        NEXTMP(mr0,mr); chsgnp(C(m),&C(mr)); mr->dl = m->dl;
                    394:                }
                    395:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    396:                if ( *pr )
                    397:                        (*pr)->sugar = p->sugar;
                    398:        }
                    399: }
                    400:
                    401: void muld(vl,p1,p2,pr)
                    402: VL vl;
                    403: DP p1,p2,*pr;
                    404: {
1.2       noro      405:        if ( ! do_weyl )
                    406:                comm_muld(vl,p1,p2,pr);
                    407:        else
                    408:                weyl_muld(vl,p1,p2,pr);
                    409: }
                    410:
                    411: void comm_muld(vl,p1,p2,pr)
                    412: VL vl;
                    413: DP p1,p2,*pr;
                    414: {
1.1       noro      415:        MP m;
                    416:        DP s,t,u;
                    417:
                    418:        if ( !p1 || !p2 )
                    419:                *pr = 0;
                    420:        else if ( OID(p1) <= O_P )
                    421:                muldc(vl,p2,(P)p1,pr);
                    422:        else if ( OID(p2) <= O_P )
                    423:                muldc(vl,p1,(P)p2,pr);
                    424:        else {
                    425:                for ( m = BDY(p2), s = 0; m; m = NEXT(m) ) {
                    426:                        muldm(vl,p1,m,&t); addd(vl,s,t,&u); s = u;
                    427:                }
                    428:                *pr = s;
                    429:        }
                    430: }
                    431:
                    432: void muldm(vl,p,m0,pr)
                    433: VL vl;
                    434: DP p;
                    435: MP m0;
                    436: DP *pr;
                    437: {
                    438:        MP m,mr,mr0;
                    439:        P c;
                    440:        DL d;
                    441:        int n;
                    442:
                    443:        if ( !p )
                    444:                *pr = 0;
                    445:        else {
                    446:                for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
                    447:                        m; m = NEXT(m) ) {
                    448:                        NEXTMP(mr0,mr);
                    449:                        if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                    450:                                mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                    451:                        else
                    452:                                mulp(vl,C(m),c,&C(mr));
                    453:                        adddl(n,m->dl,d,&mr->dl);
                    454:                }
                    455:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    456:                if ( *pr )
                    457:                        (*pr)->sugar = p->sugar + m0->dl->td;
1.2       noro      458:        }
                    459: }
                    460:
                    461: void weyl_muld(vl,p1,p2,pr)
                    462: VL vl;
                    463: DP p1,p2,*pr;
                    464: {
                    465:        MP m;
                    466:        DP s,t,u;
                    467:
                    468:        if ( !p1 || !p2 )
                    469:                *pr = 0;
                    470:        else if ( OID(p1) <= O_P )
                    471:                muldc(vl,p2,(P)p1,pr);
                    472:        else if ( OID(p2) <= O_P )
                    473:                muldc(vl,p1,(P)p2,pr);
                    474:        else {
                    475:                for ( m = BDY(p2), s = 0; m; m = NEXT(m) ) {
                    476:                        weyl_muldm(vl,p1,m,&t); addd(vl,s,t,&u); s = u;
                    477:                }
                    478:                *pr = s;
                    479:        }
                    480: }
                    481:
                    482: void weyl_muldm(vl,p,m0,pr)
                    483: VL vl;
                    484: DP p;
                    485: MP m0;
                    486: DP *pr;
                    487: {
                    488:        DP r,t,t1;
                    489:        MP m;
                    490:        int n;
                    491:
                    492:        if ( !p )
                    493:                *pr = 0;
                    494:        else {
                    495:                for ( r = 0, m = BDY(p), n = NV(p); m; m = NEXT(m) ) {
                    496:                        weyl_mulmm(vl,m,m0,n,&t);
                    497:                        addd(vl,r,t,&t1); r = t1;
                    498:                }
                    499:                if ( r )
                    500:                        r->sugar = p->sugar + m0->dl->td;
                    501:                *pr = r;
                    502:        }
                    503: }
                    504:
                    505: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
                    506:
                    507: void weyl_mulmm(vl,m0,m1,n,pr)
                    508: VL vl;
                    509: MP m0,m1;
                    510: int n;
                    511: DP *pr;
                    512: {
                    513:        MP m,mr,mr0;
                    514:        DP r,t,t1;
                    515:        P c,c0,c1,cc;
                    516:        DL d,d0,d1;
                    517:        int i,j,a,b,k,l,n2,s,min;
                    518:        Q *tab;
                    519:
                    520:        if ( !m0 || !m1 )
                    521:                *pr = 0;
                    522:        else {
                    523:                c0 = C(m0); c1 = C(m1);
                    524:                mulp(vl,c0,c1,&c);
                    525:                d0 = m0->dl; d1 = m1->dl;
                    526:                n2 = n>>1;
                    527:                for ( i = 0, r = (DP)ONE; i < n2; i++ ) {
                    528:                        a = d0->d[i]; b = d1->d[n2+i];
                    529:                        k = d0->d[n2+i]; l = d1->d[i];
                    530:                        /* compute xi^a*(Di^k*xi^l)*Di^b */
                    531:                        min = MIN(k,l);
                    532:                        tab = (Q *)ALLOCA((min+1)*sizeof(Q));
                    533:                        mkwc(k,l,tab);
                    534:                        for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
                    535:                                NEXTMP(mr0,mr);
                    536:                                NEWDL(d,n);
                    537:                                d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
                    538:                                d->td = d->d[i]+d->d[n2+i]; /* XXX */
                    539:                                s = MAX(s,d->td); /* XXX */
                    540:                                mr->c = (P)tab[j];
                    541:                                mr->dl = d;
                    542:                        }
                    543:                        if ( mr0 )
                    544:                                NEXT(mr) = 0;
                    545:                        MKDP(n,mr0,t);
                    546:                        if ( t )
                    547:                                t->sugar = s;
                    548:                        comm_muld(vl,r,t,&t1); r = t1;
                    549:                }
                    550:                muldc(vl,r,c,pr);
1.1       noro      551:        }
                    552: }
                    553:
                    554: void muldc(vl,p,c,pr)
                    555: VL vl;
                    556: DP p;
                    557: P c;
                    558: DP *pr;
                    559: {
                    560:        MP m,mr,mr0;
                    561:
                    562:        if ( !p || !c )
                    563:                *pr = 0;
                    564:        else if ( NUM(c) && UNIQ((Q)c) )
                    565:                *pr = p;
                    566:        else if ( NUM(c) && MUNIQ((Q)c) )
                    567:                chsgnd(p,pr);
                    568:        else {
                    569:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    570:                        NEXTMP(mr0,mr);
                    571:                        if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                    572:                                mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                    573:                        else
                    574:                                mulp(vl,C(m),c,&C(mr));
                    575:                        mr->dl = m->dl;
                    576:                }
                    577:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    578:                if ( *pr )
                    579:                        (*pr)->sugar = p->sugar;
                    580:        }
                    581: }
                    582:
                    583: void divsdc(vl,p,c,pr)
                    584: VL vl;
                    585: DP p;
                    586: P c;
                    587: DP *pr;
                    588: {
                    589:        MP m,mr,mr0;
                    590:
                    591:        if ( !c )
                    592:                error("disvsdc : division by 0");
                    593:        else if ( !p )
                    594:                *pr = 0;
                    595:        else {
                    596:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    597:                        NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
                    598:                }
                    599:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    600:                if ( *pr )
                    601:                        (*pr)->sugar = p->sugar;
                    602:        }
                    603: }
                    604:
                    605: void adddl(n,d1,d2,dr)
                    606: int n;
                    607: DL d1,d2;
                    608: DL *dr;
                    609: {
                    610:        DL dt;
                    611:        int i;
                    612:
                    613:        if ( !d1->td )
                    614:                *dr = d2;
                    615:        else if ( !d2->td )
                    616:                *dr = d1;
                    617:        else {
                    618:                *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
                    619:                dt->td = d1->td + d2->td;
                    620:                for ( i = 0; i < n; i++ )
                    621:                        dt->d[i] = d1->d[i]+d2->d[i];
                    622:        }
                    623: }
                    624:
                    625: int compd(vl,p1,p2)
                    626: VL vl;
                    627: DP p1,p2;
                    628: {
                    629:        int n,t;
                    630:        MP m1,m2;
                    631:
                    632:        if ( !p1 )
                    633:                return p2 ? -1 : 0;
                    634:        else if ( !p2 )
                    635:                return 1;
                    636:        else {
                    637:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
                    638:                        m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
                    639:                        if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
                    640:                                (t = compp(vl,C(m1),C(m2)) ) )
                    641:                                return t;
                    642:                if ( m1 )
                    643:                        return 1;
                    644:                else if ( m2 )
                    645:                        return -1;
                    646:                else
                    647:                        return 0;
                    648:        }
                    649: }
                    650:
                    651: int cmpdl_lex(n,d1,d2)
                    652: int n;
                    653: DL d1,d2;
                    654: {
                    655:        int i;
                    656:
                    657:        for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
                    658:        return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
                    659: }
                    660:
                    661: int cmpdl_revlex(n,d1,d2)
                    662: int n;
                    663: DL d1,d2;
                    664: {
                    665:        int i;
                    666:
                    667:        for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
                    668:        return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                    669: }
                    670:
                    671: int cmpdl_gradlex(n,d1,d2)
                    672: int n;
                    673: DL d1,d2;
                    674: {
                    675:        if ( d1->td > d2->td )
                    676:                return 1;
                    677:        else if ( d1->td < d2->td )
                    678:                return -1;
                    679:        else
                    680:                return cmpdl_lex(n,d1,d2);
                    681: }
                    682:
                    683: int cmpdl_revgradlex(n,d1,d2)
                    684: int n;
                    685: DL d1,d2;
                    686: {
                    687:        if ( d1->td > d2->td )
                    688:                return 1;
                    689:        else if ( d1->td < d2->td )
                    690:                return -1;
                    691:        else
                    692:                return cmpdl_revlex(n,d1,d2);
                    693: }
                    694:
                    695: int cmpdl_blex(n,d1,d2)
                    696: int n;
                    697: DL d1,d2;
                    698: {
                    699:        int c;
                    700:
                    701:        if ( c = cmpdl_lex(n-1,d1,d2) )
                    702:                return c;
                    703:        else {
                    704:                c = d1->d[n-1] - d2->d[n-1];
                    705:                return c > 0 ? 1 : c < 0 ? -1 : 0;
                    706:        }
                    707: }
                    708:
                    709: int cmpdl_bgradlex(n,d1,d2)
                    710: int n;
                    711: DL d1,d2;
                    712: {
                    713:        int e1,e2,c;
                    714:
                    715:        e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
                    716:        if ( e1 > e2 )
                    717:                return 1;
                    718:        else if ( e1 < e2 )
                    719:                return -1;
                    720:        else {
                    721:                c = cmpdl_lex(n-1,d1,d2);
                    722:                if ( c )
                    723:                        return c;
                    724:                else
                    725:                        return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
                    726:        }
                    727: }
                    728:
                    729: int cmpdl_brevgradlex(n,d1,d2)
                    730: int n;
                    731: DL d1,d2;
                    732: {
                    733:        int e1,e2,c;
                    734:
                    735:        e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
                    736:        if ( e1 > e2 )
                    737:                return 1;
                    738:        else if ( e1 < e2 )
                    739:                return -1;
                    740:        else {
                    741:                c = cmpdl_revlex(n-1,d1,d2);
                    742:                if ( c )
                    743:                        return c;
                    744:                else
                    745:                        return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
                    746:        }
                    747: }
                    748:
                    749: int cmpdl_brevrev(n,d1,d2)
                    750: int n;
                    751: DL d1,d2;
                    752: {
                    753:        int e1,e2,f1,f2,c,i;
                    754:
                    755:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                    756:                e1 += d1->d[i]; e2 += d2->d[i];
                    757:        }
                    758:        f1 = d1->td - e1; f2 = d2->td - e2;
                    759:        if ( e1 > e2 )
                    760:                return 1;
                    761:        else if ( e1 < e2 )
                    762:                return -1;
                    763:        else {
                    764:                c = cmpdl_revlex(dp_nelim,d1,d2);
                    765:                if ( c )
                    766:                        return c;
                    767:                else if ( f1 > f2 )
                    768:                        return 1;
                    769:                else if ( f1 < f2 )
                    770:                        return -1;
                    771:                else {
                    772:                        for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                    773:                        return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                    774:                }
                    775:        }
                    776: }
                    777:
                    778: int cmpdl_bgradrev(n,d1,d2)
                    779: int n;
                    780: DL d1,d2;
                    781: {
                    782:        int e1,e2,f1,f2,c,i;
                    783:
                    784:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                    785:                e1 += d1->d[i]; e2 += d2->d[i];
                    786:        }
                    787:        f1 = d1->td - e1; f2 = d2->td - e2;
                    788:        if ( e1 > e2 )
                    789:                return 1;
                    790:        else if ( e1 < e2 )
                    791:                return -1;
                    792:        else {
                    793:                c = cmpdl_lex(dp_nelim,d1,d2);
                    794:                if ( c )
                    795:                        return c;
                    796:                else if ( f1 > f2 )
                    797:                        return 1;
                    798:                else if ( f1 < f2 )
                    799:                        return -1;
                    800:                else {
                    801:                        for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                    802:                        return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                    803:                }
                    804:        }
                    805: }
                    806:
                    807: int cmpdl_blexrev(n,d1,d2)
                    808: int n;
                    809: DL d1,d2;
                    810: {
                    811:        int e1,e2,f1,f2,c,i;
                    812:
                    813:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                    814:                e1 += d1->d[i]; e2 += d2->d[i];
                    815:        }
                    816:        f1 = d1->td - e1; f2 = d2->td - e2;
                    817:        c = cmpdl_lex(dp_nelim,d1,d2);
                    818:        if ( c )
                    819:                return c;
                    820:        else if ( f1 > f2 )
                    821:                return 1;
                    822:        else if ( f1 < f2 )
                    823:                return -1;
                    824:        else {
                    825:                for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                    826:                return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                    827:        }
                    828: }
                    829:
                    830: int cmpdl_elim(n,d1,d2)
                    831: int n;
                    832: DL d1,d2;
                    833: {
                    834:        int e1,e2,i;
                    835:
                    836:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                    837:                e1 += d1->d[i]; e2 += d2->d[i];
                    838:        }
                    839:        if ( e1 > e2 )
                    840:                return 1;
                    841:        else if ( e1 < e2 )
                    842:                return -1;
                    843:        else
                    844:                return cmpdl_revgradlex(n,d1,d2);
                    845: }
                    846:
                    847: int cmpdl_order_pair(n,d1,d2)
                    848: int n;
                    849: DL d1,d2;
                    850: {
                    851:        int e1,e2,i,j,l;
                    852:        int *t1,*t2;
                    853:        int len;
                    854:        struct order_pair *pair;
                    855:
                    856:        len = dp_current_spec.ord.block.length;
                    857:        pair = dp_current_spec.ord.block.order_pair;
                    858:
                    859:        for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
                    860:                l = pair[i].length;
                    861:                switch ( pair[i].order ) {
                    862:                        case 0:
                    863:                                for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
                    864:                                        e1 += t1[j]; e2 += t2[j];
                    865:                                }
                    866:                                if ( e1 > e2 )
                    867:                                        return 1;
                    868:                                else if ( e1 < e2 )
                    869:                                        return -1;
                    870:                                else {
                    871:                                        for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
                    872:                                        if ( j >= 0 )
                    873:                                                return t1[j] < t2[j] ? 1 : -1;
                    874:                                }
                    875:                                break;
                    876:                        case 1:
                    877:                                for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
                    878:                                        e1 += t1[j]; e2 += t2[j];
                    879:                                }
                    880:                                if ( e1 > e2 )
                    881:                                        return 1;
                    882:                                else if ( e1 < e2 )
                    883:                                        return -1;
                    884:                                else {
                    885:                                        for ( j = 0; j < l && t1[j] == t2[j]; j++ );
                    886:                                        if ( j < l )
                    887:                                                return t1[j] > t2[j] ? 1 : -1;
                    888:                                }
                    889:                                break;
                    890:                        case 2:
                    891:                                for ( j = 0; j < l && t1[j] == t2[j]; j++ );
                    892:                                if ( j < l )
                    893:                                        return t1[j] > t2[j] ? 1 : -1;
                    894:                                break;
                    895:                        default:
                    896:                                error("cmpdl_order_pair : invalid order"); break;
                    897:                }
                    898:                t1 += l; t2 += l;
                    899:        }
                    900:        return 0;
                    901: }
                    902:
                    903: int cmpdl_matrix(n,d1,d2)
                    904: int n;
                    905: DL d1,d2;
                    906: {
                    907:        int *v,*w,*t1,*t2;
                    908:        int s,i,j,len;
                    909:        int **matrix;
                    910:
                    911:        for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
                    912:                w[i] = t1[i]-t2[i];
                    913:        len = dp_current_spec.ord.matrix.row;
                    914:        matrix = dp_current_spec.ord.matrix.matrix;
                    915:        for ( j = 0; j < len; j++ ) {
                    916:                v = matrix[j];
                    917:                for ( i = 0, s = 0; i < n; i++ )
                    918:                        s += v[i]*w[i];
                    919:                if ( s > 0 )
                    920:                        return 1;
                    921:                else if ( s < 0 )
                    922:                        return -1;
                    923:        }
                    924:        return 0;
                    925: }

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