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

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

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