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

1.4     ! noro        1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.3 2000/04/13 06:01:02 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;
1.4     ! noro      417:        int i,l;
        !           418:        MP *w;
1.1       noro      419:
                    420:        if ( !p1 || !p2 )
                    421:                *pr = 0;
                    422:        else if ( OID(p1) <= O_P )
                    423:                muldc(vl,p2,(P)p1,pr);
                    424:        else if ( OID(p2) <= O_P )
                    425:                muldc(vl,p1,(P)p2,pr);
                    426:        else {
1.4     ! noro      427:                for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
        !           428:                w = ALLOCA(l*sizeof(MP));
        !           429:                for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
        !           430:                        w[i] = m;
        !           431:                for ( s = 0, i = l-1; i >= 0; i-- ) {
        !           432:                        muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
1.1       noro      433:                }
                    434:                *pr = s;
                    435:        }
                    436: }
                    437:
                    438: void muldm(vl,p,m0,pr)
                    439: VL vl;
                    440: DP p;
                    441: MP m0;
                    442: DP *pr;
                    443: {
                    444:        MP m,mr,mr0;
                    445:        P c;
                    446:        DL d;
                    447:        int n;
                    448:
                    449:        if ( !p )
                    450:                *pr = 0;
                    451:        else {
                    452:                for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
                    453:                        m; m = NEXT(m) ) {
                    454:                        NEXTMP(mr0,mr);
                    455:                        if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                    456:                                mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                    457:                        else
                    458:                                mulp(vl,C(m),c,&C(mr));
                    459:                        adddl(n,m->dl,d,&mr->dl);
                    460:                }
                    461:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    462:                if ( *pr )
                    463:                        (*pr)->sugar = p->sugar + m0->dl->td;
1.2       noro      464:        }
                    465: }
                    466:
                    467: void weyl_muld(vl,p1,p2,pr)
                    468: VL vl;
                    469: DP p1,p2,*pr;
                    470: {
                    471:        MP m;
                    472:        DP s,t,u;
1.4     ! noro      473:        int i,l;
        !           474:        MP *w;
1.2       noro      475:
                    476:        if ( !p1 || !p2 )
                    477:                *pr = 0;
                    478:        else if ( OID(p1) <= O_P )
                    479:                muldc(vl,p2,(P)p1,pr);
                    480:        else if ( OID(p2) <= O_P )
                    481:                muldc(vl,p1,(P)p2,pr);
                    482:        else {
1.4     ! noro      483:                for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
        !           484:                w = ALLOCA(l*sizeof(MP));
        !           485:                for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
        !           486:                        w[i] = m;
        !           487:                for ( s = 0, i = l-1; i >= 0; i-- ) {
        !           488:                        weyl_muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
1.2       noro      489:                }
                    490:                *pr = s;
                    491:        }
                    492: }
                    493:
                    494: void weyl_muldm(vl,p,m0,pr)
                    495: VL vl;
                    496: DP p;
                    497: MP m0;
                    498: DP *pr;
                    499: {
                    500:        DP r,t,t1;
                    501:        MP m;
1.4     ! noro      502:        int n,l,i;
        !           503:        MP *w;
1.2       noro      504:
                    505:        if ( !p )
                    506:                *pr = 0;
                    507:        else {
1.4     ! noro      508:                for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
        !           509:                w = ALLOCA(l*sizeof(MP));
        !           510:                for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
        !           511:                        w[i] = m;
        !           512:                for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
        !           513:                        weyl_mulmm(vl,w[i],m0,n,&t);
1.2       noro      514:                        addd(vl,r,t,&t1); r = t1;
                    515:                }
                    516:                if ( r )
                    517:                        r->sugar = p->sugar + m0->dl->td;
                    518:                *pr = r;
                    519:        }
                    520: }
                    521:
                    522: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
                    523:
                    524: void weyl_mulmm(vl,m0,m1,n,pr)
                    525: VL vl;
                    526: MP m0,m1;
                    527: int n;
                    528: DP *pr;
                    529: {
                    530:        MP m,mr,mr0;
                    531:        DP r,t,t1;
                    532:        P c,c0,c1,cc;
                    533:        DL d,d0,d1;
                    534:        int i,j,a,b,k,l,n2,s,min;
                    535:        Q *tab;
                    536:
                    537:        if ( !m0 || !m1 )
                    538:                *pr = 0;
                    539:        else {
                    540:                c0 = C(m0); c1 = C(m1);
                    541:                mulp(vl,c0,c1,&c);
                    542:                d0 = m0->dl; d1 = m1->dl;
                    543:                n2 = n>>1;
                    544:                for ( i = 0, r = (DP)ONE; i < n2; i++ ) {
                    545:                        a = d0->d[i]; b = d1->d[n2+i];
                    546:                        k = d0->d[n2+i]; l = d1->d[i];
                    547:                        /* compute xi^a*(Di^k*xi^l)*Di^b */
                    548:                        min = MIN(k,l);
                    549:                        tab = (Q *)ALLOCA((min+1)*sizeof(Q));
                    550:                        mkwc(k,l,tab);
                    551:                        for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
                    552:                                NEXTMP(mr0,mr);
                    553:                                NEWDL(d,n);
                    554:                                d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
                    555:                                d->td = d->d[i]+d->d[n2+i]; /* XXX */
                    556:                                s = MAX(s,d->td); /* XXX */
                    557:                                mr->c = (P)tab[j];
                    558:                                mr->dl = d;
                    559:                        }
                    560:                        if ( mr0 )
                    561:                                NEXT(mr) = 0;
                    562:                        MKDP(n,mr0,t);
                    563:                        if ( t )
                    564:                                t->sugar = s;
                    565:                        comm_muld(vl,r,t,&t1); r = t1;
                    566:                }
                    567:                muldc(vl,r,c,pr);
1.1       noro      568:        }
                    569: }
                    570:
                    571: void muldc(vl,p,c,pr)
                    572: VL vl;
                    573: DP p;
                    574: P c;
                    575: DP *pr;
                    576: {
                    577:        MP m,mr,mr0;
                    578:
                    579:        if ( !p || !c )
                    580:                *pr = 0;
                    581:        else if ( NUM(c) && UNIQ((Q)c) )
                    582:                *pr = p;
                    583:        else if ( NUM(c) && MUNIQ((Q)c) )
                    584:                chsgnd(p,pr);
                    585:        else {
                    586:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    587:                        NEXTMP(mr0,mr);
                    588:                        if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
                    589:                                mulq((Q)C(m),(Q)c,(Q *)&C(mr));
                    590:                        else
                    591:                                mulp(vl,C(m),c,&C(mr));
                    592:                        mr->dl = m->dl;
                    593:                }
                    594:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    595:                if ( *pr )
                    596:                        (*pr)->sugar = p->sugar;
                    597:        }
                    598: }
                    599:
                    600: void divsdc(vl,p,c,pr)
                    601: VL vl;
                    602: DP p;
                    603: P c;
                    604: DP *pr;
                    605: {
                    606:        MP m,mr,mr0;
                    607:
                    608:        if ( !c )
                    609:                error("disvsdc : division by 0");
                    610:        else if ( !p )
                    611:                *pr = 0;
                    612:        else {
                    613:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    614:                        NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
                    615:                }
                    616:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    617:                if ( *pr )
                    618:                        (*pr)->sugar = p->sugar;
                    619:        }
                    620: }
                    621:
                    622: void adddl(n,d1,d2,dr)
                    623: int n;
                    624: DL d1,d2;
                    625: DL *dr;
                    626: {
                    627:        DL dt;
                    628:        int i;
                    629:
                    630:        if ( !d1->td )
                    631:                *dr = d2;
                    632:        else if ( !d2->td )
                    633:                *dr = d1;
                    634:        else {
                    635:                *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
                    636:                dt->td = d1->td + d2->td;
                    637:                for ( i = 0; i < n; i++ )
                    638:                        dt->d[i] = d1->d[i]+d2->d[i];
                    639:        }
                    640: }
                    641:
                    642: int compd(vl,p1,p2)
                    643: VL vl;
                    644: DP p1,p2;
                    645: {
                    646:        int n,t;
                    647:        MP m1,m2;
                    648:
                    649:        if ( !p1 )
                    650:                return p2 ? -1 : 0;
                    651:        else if ( !p2 )
                    652:                return 1;
                    653:        else {
                    654:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
                    655:                        m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
                    656:                        if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
                    657:                                (t = compp(vl,C(m1),C(m2)) ) )
                    658:                                return t;
                    659:                if ( m1 )
                    660:                        return 1;
                    661:                else if ( m2 )
                    662:                        return -1;
                    663:                else
                    664:                        return 0;
                    665:        }
                    666: }
                    667:
                    668: int cmpdl_lex(n,d1,d2)
                    669: int n;
                    670: DL d1,d2;
                    671: {
                    672:        int i;
                    673:
                    674:        for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
                    675:        return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
                    676: }
                    677:
                    678: int cmpdl_revlex(n,d1,d2)
                    679: int n;
                    680: DL d1,d2;
                    681: {
                    682:        int i;
                    683:
                    684:        for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
                    685:        return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                    686: }
                    687:
                    688: int cmpdl_gradlex(n,d1,d2)
                    689: int n;
                    690: DL d1,d2;
                    691: {
                    692:        if ( d1->td > d2->td )
                    693:                return 1;
                    694:        else if ( d1->td < d2->td )
                    695:                return -1;
                    696:        else
                    697:                return cmpdl_lex(n,d1,d2);
                    698: }
                    699:
                    700: int cmpdl_revgradlex(n,d1,d2)
                    701: int n;
                    702: DL d1,d2;
                    703: {
                    704:        if ( d1->td > d2->td )
                    705:                return 1;
                    706:        else if ( d1->td < d2->td )
                    707:                return -1;
                    708:        else
                    709:                return cmpdl_revlex(n,d1,d2);
                    710: }
                    711:
                    712: int cmpdl_blex(n,d1,d2)
                    713: int n;
                    714: DL d1,d2;
                    715: {
                    716:        int c;
                    717:
                    718:        if ( c = cmpdl_lex(n-1,d1,d2) )
                    719:                return c;
                    720:        else {
                    721:                c = d1->d[n-1] - d2->d[n-1];
                    722:                return c > 0 ? 1 : c < 0 ? -1 : 0;
                    723:        }
                    724: }
                    725:
                    726: int cmpdl_bgradlex(n,d1,d2)
                    727: int n;
                    728: DL d1,d2;
                    729: {
                    730:        int e1,e2,c;
                    731:
                    732:        e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
                    733:        if ( e1 > e2 )
                    734:                return 1;
                    735:        else if ( e1 < e2 )
                    736:                return -1;
                    737:        else {
                    738:                c = cmpdl_lex(n-1,d1,d2);
                    739:                if ( c )
                    740:                        return c;
                    741:                else
                    742:                        return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
                    743:        }
                    744: }
                    745:
                    746: int cmpdl_brevgradlex(n,d1,d2)
                    747: int n;
                    748: DL d1,d2;
                    749: {
                    750:        int e1,e2,c;
                    751:
                    752:        e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
                    753:        if ( e1 > e2 )
                    754:                return 1;
                    755:        else if ( e1 < e2 )
                    756:                return -1;
                    757:        else {
                    758:                c = cmpdl_revlex(n-1,d1,d2);
                    759:                if ( c )
                    760:                        return c;
                    761:                else
                    762:                        return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
                    763:        }
                    764: }
                    765:
                    766: int cmpdl_brevrev(n,d1,d2)
                    767: int n;
                    768: DL d1,d2;
                    769: {
                    770:        int e1,e2,f1,f2,c,i;
                    771:
                    772:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                    773:                e1 += d1->d[i]; e2 += d2->d[i];
                    774:        }
                    775:        f1 = d1->td - e1; f2 = d2->td - e2;
                    776:        if ( e1 > e2 )
                    777:                return 1;
                    778:        else if ( e1 < e2 )
                    779:                return -1;
                    780:        else {
                    781:                c = cmpdl_revlex(dp_nelim,d1,d2);
                    782:                if ( c )
                    783:                        return c;
                    784:                else if ( f1 > f2 )
                    785:                        return 1;
                    786:                else if ( f1 < f2 )
                    787:                        return -1;
                    788:                else {
                    789:                        for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                    790:                        return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                    791:                }
                    792:        }
                    793: }
                    794:
                    795: int cmpdl_bgradrev(n,d1,d2)
                    796: int n;
                    797: DL d1,d2;
                    798: {
                    799:        int e1,e2,f1,f2,c,i;
                    800:
                    801:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                    802:                e1 += d1->d[i]; e2 += d2->d[i];
                    803:        }
                    804:        f1 = d1->td - e1; f2 = d2->td - e2;
                    805:        if ( e1 > e2 )
                    806:                return 1;
                    807:        else if ( e1 < e2 )
                    808:                return -1;
                    809:        else {
                    810:                c = cmpdl_lex(dp_nelim,d1,d2);
                    811:                if ( c )
                    812:                        return c;
                    813:                else if ( f1 > f2 )
                    814:                        return 1;
                    815:                else if ( f1 < f2 )
                    816:                        return -1;
                    817:                else {
                    818:                        for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                    819:                        return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                    820:                }
                    821:        }
                    822: }
                    823:
                    824: int cmpdl_blexrev(n,d1,d2)
                    825: int n;
                    826: DL d1,d2;
                    827: {
                    828:        int e1,e2,f1,f2,c,i;
                    829:
                    830:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                    831:                e1 += d1->d[i]; e2 += d2->d[i];
                    832:        }
                    833:        f1 = d1->td - e1; f2 = d2->td - e2;
                    834:        c = cmpdl_lex(dp_nelim,d1,d2);
                    835:        if ( c )
                    836:                return c;
                    837:        else if ( f1 > f2 )
                    838:                return 1;
                    839:        else if ( f1 < f2 )
                    840:                return -1;
                    841:        else {
                    842:                for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
                    843:                return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
                    844:        }
                    845: }
                    846:
                    847: int cmpdl_elim(n,d1,d2)
                    848: int n;
                    849: DL d1,d2;
                    850: {
                    851:        int e1,e2,i;
                    852:
                    853:        for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
                    854:                e1 += d1->d[i]; e2 += d2->d[i];
                    855:        }
                    856:        if ( e1 > e2 )
                    857:                return 1;
                    858:        else if ( e1 < e2 )
                    859:                return -1;
                    860:        else
                    861:                return cmpdl_revgradlex(n,d1,d2);
                    862: }
                    863:
                    864: int cmpdl_order_pair(n,d1,d2)
                    865: int n;
                    866: DL d1,d2;
                    867: {
                    868:        int e1,e2,i,j,l;
                    869:        int *t1,*t2;
                    870:        int len;
                    871:        struct order_pair *pair;
                    872:
                    873:        len = dp_current_spec.ord.block.length;
                    874:        pair = dp_current_spec.ord.block.order_pair;
                    875:
                    876:        for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
                    877:                l = pair[i].length;
                    878:                switch ( pair[i].order ) {
                    879:                        case 0:
                    880:                                for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
                    881:                                        e1 += t1[j]; e2 += t2[j];
                    882:                                }
                    883:                                if ( e1 > e2 )
                    884:                                        return 1;
                    885:                                else if ( e1 < e2 )
                    886:                                        return -1;
                    887:                                else {
                    888:                                        for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
                    889:                                        if ( j >= 0 )
                    890:                                                return t1[j] < t2[j] ? 1 : -1;
                    891:                                }
                    892:                                break;
                    893:                        case 1:
                    894:                                for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
                    895:                                        e1 += t1[j]; e2 += t2[j];
                    896:                                }
                    897:                                if ( e1 > e2 )
                    898:                                        return 1;
                    899:                                else if ( e1 < e2 )
                    900:                                        return -1;
                    901:                                else {
                    902:                                        for ( j = 0; j < l && t1[j] == t2[j]; j++ );
                    903:                                        if ( j < l )
                    904:                                                return t1[j] > t2[j] ? 1 : -1;
                    905:                                }
                    906:                                break;
                    907:                        case 2:
                    908:                                for ( j = 0; j < l && t1[j] == t2[j]; j++ );
                    909:                                if ( j < l )
                    910:                                        return t1[j] > t2[j] ? 1 : -1;
                    911:                                break;
                    912:                        default:
                    913:                                error("cmpdl_order_pair : invalid order"); break;
                    914:                }
                    915:                t1 += l; t2 += l;
                    916:        }
                    917:        return 0;
                    918: }
                    919:
                    920: int cmpdl_matrix(n,d1,d2)
                    921: int n;
                    922: DL d1,d2;
                    923: {
                    924:        int *v,*w,*t1,*t2;
                    925:        int s,i,j,len;
                    926:        int **matrix;
                    927:
                    928:        for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
                    929:                w[i] = t1[i]-t2[i];
                    930:        len = dp_current_spec.ord.matrix.row;
                    931:        matrix = dp_current_spec.ord.matrix.matrix;
                    932:        for ( j = 0; j < len; j++ ) {
                    933:                v = matrix[j];
                    934:                for ( i = 0, s = 0; i < n; i++ )
                    935:                        s += v[i]*w[i];
                    936:                if ( s > 0 )
                    937:                        return 1;
                    938:                else if ( s < 0 )
                    939:                        return -1;
                    940:        }
                    941:        return 0;
                    942: }

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