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

1.2     ! noro        1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.1.1.1 1999/12/03 07:39:08 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);
                    315:        }
                    316: }
                    317:
                    318: void subd(vl,p1,p2,pr)
                    319: VL vl;
                    320: DP p1,p2,*pr;
                    321: {
                    322:        DP t;
                    323:
                    324:        if ( !p2 )
                    325:                *pr = p1;
                    326:        else {
                    327:                chsgnd(p2,&t); addd(vl,p1,t,pr);
                    328:        }
                    329: }
                    330:
                    331: void chsgnd(p,pr)
                    332: DP p,*pr;
                    333: {
                    334:        MP m,mr,mr0;
                    335:
                    336:        if ( !p )
                    337:                *pr = 0;
                    338:        else {
                    339:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    340:                        NEXTMP(mr0,mr); chsgnp(C(m),&C(mr)); mr->dl = m->dl;
                    341:                }
                    342:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    343:                if ( *pr )
                    344:                        (*pr)->sugar = p->sugar;
                    345:        }
                    346: }
                    347:
                    348: void muld(vl,p1,p2,pr)
                    349: VL vl;
                    350: DP p1,p2,*pr;
                    351: {
1.2     ! noro      352:        if ( ! do_weyl )
        !           353:                comm_muld(vl,p1,p2,pr);
        !           354:        else
        !           355:                weyl_muld(vl,p1,p2,pr);
        !           356: }
        !           357:
        !           358: void comm_muld(vl,p1,p2,pr)
        !           359: VL vl;
        !           360: DP p1,p2,*pr;
        !           361: {
1.1       noro      362:        MP m;
                    363:        DP s,t,u;
                    364:
                    365:        if ( !p1 || !p2 )
                    366:                *pr = 0;
                    367:        else if ( OID(p1) <= O_P )
                    368:                muldc(vl,p2,(P)p1,pr);
                    369:        else if ( OID(p2) <= O_P )
                    370:                muldc(vl,p1,(P)p2,pr);
                    371:        else {
                    372:                for ( m = BDY(p2), s = 0; m; m = NEXT(m) ) {
                    373:                        muldm(vl,p1,m,&t); addd(vl,s,t,&u); s = u;
                    374:                }
                    375:                *pr = s;
                    376:        }
                    377: }
                    378:
                    379: void muldm(vl,p,m0,pr)
                    380: VL vl;
                    381: DP p;
                    382: MP m0;
                    383: DP *pr;
                    384: {
                    385:        MP m,mr,mr0;
                    386:        P c;
                    387:        DL d;
                    388:        int n;
                    389:
                    390:        if ( !p )
                    391:                *pr = 0;
                    392:        else {
                    393:                for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
                    394:                        m; m = NEXT(m) ) {
                    395:                        NEXTMP(mr0,mr);
                    396:                        if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                    397:                                mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                    398:                        else
                    399:                                mulp(vl,C(m),c,&C(mr));
                    400:                        adddl(n,m->dl,d,&mr->dl);
                    401:                }
                    402:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    403:                if ( *pr )
                    404:                        (*pr)->sugar = p->sugar + m0->dl->td;
1.2     ! noro      405:        }
        !           406: }
        !           407:
        !           408: void weyl_muld(vl,p1,p2,pr)
        !           409: VL vl;
        !           410: DP p1,p2,*pr;
        !           411: {
        !           412:        MP m;
        !           413:        DP s,t,u;
        !           414:
        !           415:        if ( !p1 || !p2 )
        !           416:                *pr = 0;
        !           417:        else if ( OID(p1) <= O_P )
        !           418:                muldc(vl,p2,(P)p1,pr);
        !           419:        else if ( OID(p2) <= O_P )
        !           420:                muldc(vl,p1,(P)p2,pr);
        !           421:        else {
        !           422:                for ( m = BDY(p2), s = 0; m; m = NEXT(m) ) {
        !           423:                        weyl_muldm(vl,p1,m,&t); addd(vl,s,t,&u); s = u;
        !           424:                }
        !           425:                *pr = s;
        !           426:        }
        !           427: }
        !           428:
        !           429: void weyl_muldm(vl,p,m0,pr)
        !           430: VL vl;
        !           431: DP p;
        !           432: MP m0;
        !           433: DP *pr;
        !           434: {
        !           435:        DP r,t,t1;
        !           436:        MP m;
        !           437:        int n;
        !           438:
        !           439:        if ( !p )
        !           440:                *pr = 0;
        !           441:        else {
        !           442:                for ( r = 0, m = BDY(p), n = NV(p); m; m = NEXT(m) ) {
        !           443:                        weyl_mulmm(vl,m,m0,n,&t);
        !           444:                        addd(vl,r,t,&t1); r = t1;
        !           445:                }
        !           446:                if ( r )
        !           447:                        r->sugar = p->sugar + m0->dl->td;
        !           448:                *pr = r;
        !           449:        }
        !           450: }
        !           451:
        !           452: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
        !           453:
        !           454: void weyl_mulmm(vl,m0,m1,n,pr)
        !           455: VL vl;
        !           456: MP m0,m1;
        !           457: int n;
        !           458: DP *pr;
        !           459: {
        !           460:        MP m,mr,mr0;
        !           461:        DP r,t,t1;
        !           462:        P c,c0,c1,cc;
        !           463:        DL d,d0,d1;
        !           464:        int i,j,a,b,k,l,n2,s,min;
        !           465:        Q *tab;
        !           466:
        !           467:        if ( !m0 || !m1 )
        !           468:                *pr = 0;
        !           469:        else {
        !           470:                c0 = C(m0); c1 = C(m1);
        !           471:                mulp(vl,c0,c1,&c);
        !           472:                d0 = m0->dl; d1 = m1->dl;
        !           473:                n2 = n>>1;
        !           474:                for ( i = 0, r = (DP)ONE; i < n2; i++ ) {
        !           475:                        a = d0->d[i]; b = d1->d[n2+i];
        !           476:                        k = d0->d[n2+i]; l = d1->d[i];
        !           477:                        /* compute xi^a*(Di^k*xi^l)*Di^b */
        !           478:                        min = MIN(k,l);
        !           479:                        tab = (Q *)ALLOCA((min+1)*sizeof(Q));
        !           480:                        mkwc(k,l,tab);
        !           481:                        for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
        !           482:                                NEXTMP(mr0,mr);
        !           483:                                NEWDL(d,n);
        !           484:                                d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
        !           485:                                d->td = d->d[i]+d->d[n2+i]; /* XXX */
        !           486:                                s = MAX(s,d->td); /* XXX */
        !           487:                                mr->c = (P)tab[j];
        !           488:                                mr->dl = d;
        !           489:                        }
        !           490:                        if ( mr0 )
        !           491:                                NEXT(mr) = 0;
        !           492:                        MKDP(n,mr0,t);
        !           493:                        if ( t )
        !           494:                                t->sugar = s;
        !           495:                        comm_muld(vl,r,t,&t1); r = t1;
        !           496:                }
        !           497:                muldc(vl,r,c,pr);
1.1       noro      498:        }
                    499: }
                    500:
                    501: void muldc(vl,p,c,pr)
                    502: VL vl;
                    503: DP p;
                    504: P c;
                    505: DP *pr;
                    506: {
                    507:        MP m,mr,mr0;
                    508:
                    509:        if ( !p || !c )
                    510:                *pr = 0;
                    511:        else if ( NUM(c) && UNIQ((Q)c) )
                    512:                *pr = p;
                    513:        else if ( NUM(c) && MUNIQ((Q)c) )
                    514:                chsgnd(p,pr);
                    515:        else {
                    516:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    517:                        NEXTMP(mr0,mr);
                    518:                        if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                    519:                                mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                    520:                        else
                    521:                                mulp(vl,C(m),c,&C(mr));
                    522:                        mr->dl = m->dl;
                    523:                }
                    524:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    525:                if ( *pr )
                    526:                        (*pr)->sugar = p->sugar;
                    527:        }
                    528: }
                    529:
                    530: void divsdc(vl,p,c,pr)
                    531: VL vl;
                    532: DP p;
                    533: P c;
                    534: DP *pr;
                    535: {
                    536:        MP m,mr,mr0;
                    537:
                    538:        if ( !c )
                    539:                error("disvsdc : division by 0");
                    540:        else if ( !p )
                    541:                *pr = 0;
                    542:        else {
                    543:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    544:                        NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
                    545:                }
                    546:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    547:                if ( *pr )
                    548:                        (*pr)->sugar = p->sugar;
                    549:        }
                    550: }
                    551:
                    552: void adddl(n,d1,d2,dr)
                    553: int n;
                    554: DL d1,d2;
                    555: DL *dr;
                    556: {
                    557:        DL dt;
                    558:        int i;
                    559:
                    560:        if ( !d1->td )
                    561:                *dr = d2;
                    562:        else if ( !d2->td )
                    563:                *dr = d1;
                    564:        else {
                    565:                *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
                    566:                dt->td = d1->td + d2->td;
                    567:                for ( i = 0; i < n; i++ )
                    568:                        dt->d[i] = d1->d[i]+d2->d[i];
                    569:        }
                    570: }
                    571:
                    572: int compd(vl,p1,p2)
                    573: VL vl;
                    574: DP p1,p2;
                    575: {
                    576:        int n,t;
                    577:        MP m1,m2;
                    578:
                    579:        if ( !p1 )
                    580:                return p2 ? -1 : 0;
                    581:        else if ( !p2 )
                    582:                return 1;
                    583:        else {
                    584:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
                    585:                        m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
                    586:                        if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
                    587:                                (t = compp(vl,C(m1),C(m2)) ) )
                    588:                                return t;
                    589:                if ( m1 )
                    590:                        return 1;
                    591:                else if ( m2 )
                    592:                        return -1;
                    593:                else
                    594:                        return 0;
                    595:        }
                    596: }
                    597:
                    598: int cmpdl_lex(n,d1,d2)
                    599: int n;
                    600: DL d1,d2;
                    601: {
                    602:        int i;
                    603:
                    604:        for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
                    605:        return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
                    606: }
                    607:
                    608: int cmpdl_revlex(n,d1,d2)
                    609: int n;
                    610: DL d1,d2;
                    611: {
                    612:        int i;
                    613:
                    614:        for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
                    615:        return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                    616: }
                    617:
                    618: int cmpdl_gradlex(n,d1,d2)
                    619: int n;
                    620: DL d1,d2;
                    621: {
                    622:        if ( d1->td > d2->td )
                    623:                return 1;
                    624:        else if ( d1->td < d2->td )
                    625:                return -1;
                    626:        else
                    627:                return cmpdl_lex(n,d1,d2);
                    628: }
                    629:
                    630: int cmpdl_revgradlex(n,d1,d2)
                    631: int n;
                    632: DL d1,d2;
                    633: {
                    634:        if ( d1->td > d2->td )
                    635:                return 1;
                    636:        else if ( d1->td < d2->td )
                    637:                return -1;
                    638:        else
                    639:                return cmpdl_revlex(n,d1,d2);
                    640: }
                    641:
                    642: int cmpdl_blex(n,d1,d2)
                    643: int n;
                    644: DL d1,d2;
                    645: {
                    646:        int c;
                    647:
                    648:        if ( c = cmpdl_lex(n-1,d1,d2) )
                    649:                return c;
                    650:        else {
                    651:                c = d1->d[n-1] - d2->d[n-1];
                    652:                return c > 0 ? 1 : c < 0 ? -1 : 0;
                    653:        }
                    654: }
                    655:
                    656: int cmpdl_bgradlex(n,d1,d2)
                    657: int n;
                    658: DL d1,d2;
                    659: {
                    660:        int e1,e2,c;
                    661:
                    662:        e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
                    663:        if ( e1 > e2 )
                    664:                return 1;
                    665:        else if ( e1 < e2 )
                    666:                return -1;
                    667:        else {
                    668:                c = cmpdl_lex(n-1,d1,d2);
                    669:                if ( c )
                    670:                        return c;
                    671:                else
                    672:                        return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
                    673:        }
                    674: }
                    675:
                    676: int cmpdl_brevgradlex(n,d1,d2)
                    677: int n;
                    678: DL d1,d2;
                    679: {
                    680:        int e1,e2,c;
                    681:
                    682:        e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
                    683:        if ( e1 > e2 )
                    684:                return 1;
                    685:        else if ( e1 < e2 )
                    686:                return -1;
                    687:        else {
                    688:                c = cmpdl_revlex(n-1,d1,d2);
                    689:                if ( c )
                    690:                        return c;
                    691:                else
                    692:                        return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
                    693:        }
                    694: }
                    695:
                    696: int cmpdl_brevrev(n,d1,d2)
                    697: int n;
                    698: DL d1,d2;
                    699: {
                    700:        int e1,e2,f1,f2,c,i;
                    701:
                    702:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                    703:                e1 += d1->d[i]; e2 += d2->d[i];
                    704:        }
                    705:        f1 = d1->td - e1; f2 = d2->td - e2;
                    706:        if ( e1 > e2 )
                    707:                return 1;
                    708:        else if ( e1 < e2 )
                    709:                return -1;
                    710:        else {
                    711:                c = cmpdl_revlex(dp_nelim,d1,d2);
                    712:                if ( c )
                    713:                        return c;
                    714:                else if ( f1 > f2 )
                    715:                        return 1;
                    716:                else if ( f1 < f2 )
                    717:                        return -1;
                    718:                else {
                    719:                        for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                    720:                        return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                    721:                }
                    722:        }
                    723: }
                    724:
                    725: int cmpdl_bgradrev(n,d1,d2)
                    726: int n;
                    727: DL d1,d2;
                    728: {
                    729:        int e1,e2,f1,f2,c,i;
                    730:
                    731:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                    732:                e1 += d1->d[i]; e2 += d2->d[i];
                    733:        }
                    734:        f1 = d1->td - e1; f2 = d2->td - e2;
                    735:        if ( e1 > e2 )
                    736:                return 1;
                    737:        else if ( e1 < e2 )
                    738:                return -1;
                    739:        else {
                    740:                c = cmpdl_lex(dp_nelim,d1,d2);
                    741:                if ( c )
                    742:                        return c;
                    743:                else if ( f1 > f2 )
                    744:                        return 1;
                    745:                else if ( f1 < f2 )
                    746:                        return -1;
                    747:                else {
                    748:                        for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                    749:                        return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                    750:                }
                    751:        }
                    752: }
                    753:
                    754: int cmpdl_blexrev(n,d1,d2)
                    755: int n;
                    756: DL d1,d2;
                    757: {
                    758:        int e1,e2,f1,f2,c,i;
                    759:
                    760:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                    761:                e1 += d1->d[i]; e2 += d2->d[i];
                    762:        }
                    763:        f1 = d1->td - e1; f2 = d2->td - e2;
                    764:        c = cmpdl_lex(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: int cmpdl_elim(n,d1,d2)
                    778: int n;
                    779: DL d1,d2;
                    780: {
                    781:        int e1,e2,i;
                    782:
                    783:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                    784:                e1 += d1->d[i]; e2 += d2->d[i];
                    785:        }
                    786:        if ( e1 > e2 )
                    787:                return 1;
                    788:        else if ( e1 < e2 )
                    789:                return -1;
                    790:        else
                    791:                return cmpdl_revgradlex(n,d1,d2);
                    792: }
                    793:
                    794: int cmpdl_order_pair(n,d1,d2)
                    795: int n;
                    796: DL d1,d2;
                    797: {
                    798:        int e1,e2,i,j,l;
                    799:        int *t1,*t2;
                    800:        int len;
                    801:        struct order_pair *pair;
                    802:
                    803:        len = dp_current_spec.ord.block.length;
                    804:        pair = dp_current_spec.ord.block.order_pair;
                    805:
                    806:        for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
                    807:                l = pair[i].length;
                    808:                switch ( pair[i].order ) {
                    809:                        case 0:
                    810:                                for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
                    811:                                        e1 += t1[j]; e2 += t2[j];
                    812:                                }
                    813:                                if ( e1 > e2 )
                    814:                                        return 1;
                    815:                                else if ( e1 < e2 )
                    816:                                        return -1;
                    817:                                else {
                    818:                                        for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
                    819:                                        if ( j >= 0 )
                    820:                                                return t1[j] < t2[j] ? 1 : -1;
                    821:                                }
                    822:                                break;
                    823:                        case 1:
                    824:                                for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
                    825:                                        e1 += t1[j]; e2 += t2[j];
                    826:                                }
                    827:                                if ( e1 > e2 )
                    828:                                        return 1;
                    829:                                else if ( e1 < e2 )
                    830:                                        return -1;
                    831:                                else {
                    832:                                        for ( j = 0; j < l && t1[j] == t2[j]; j++ );
                    833:                                        if ( j < l )
                    834:                                                return t1[j] > t2[j] ? 1 : -1;
                    835:                                }
                    836:                                break;
                    837:                        case 2:
                    838:                                for ( j = 0; j < l && t1[j] == t2[j]; j++ );
                    839:                                if ( j < l )
                    840:                                        return t1[j] > t2[j] ? 1 : -1;
                    841:                                break;
                    842:                        default:
                    843:                                error("cmpdl_order_pair : invalid order"); break;
                    844:                }
                    845:                t1 += l; t2 += l;
                    846:        }
                    847:        return 0;
                    848: }
                    849:
                    850: int cmpdl_matrix(n,d1,d2)
                    851: int n;
                    852: DL d1,d2;
                    853: {
                    854:        int *v,*w,*t1,*t2;
                    855:        int s,i,j,len;
                    856:        int **matrix;
                    857:
                    858:        for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
                    859:                w[i] = t1[i]-t2[i];
                    860:        len = dp_current_spec.ord.matrix.row;
                    861:        matrix = dp_current_spec.ord.matrix.matrix;
                    862:        for ( j = 0; j < len; j++ ) {
                    863:                v = matrix[j];
                    864:                for ( i = 0, s = 0; i < n; i++ )
                    865:                        s += v[i]*w[i];
                    866:                if ( s > 0 )
                    867:                        return 1;
                    868:                else if ( s < 0 )
                    869:                        return -1;
                    870:        }
                    871:        return 0;
                    872: }

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