[BACK]Return to dp-supp.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Annotation of OpenXM_contrib2/asir2000/builtin/dp-supp.c, Revision 1.47

1.2       noro        1: /*
                      2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
                      3:  * All rights reserved.
                      4:  *
                      5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
                      6:  * non-exclusive and royalty-free license to use, copy, modify and
                      7:  * redistribute, solely for non-commercial and non-profit purposes, the
                      8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
                      9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
                     10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
                     11:  * third party developer retains all rights, including but not limited to
                     12:  * copyrights, in and to the SOFTWARE.
                     13:  *
                     14:  * (1) FLL does not grant you a license in any way for commercial
                     15:  * purposes. You may use the SOFTWARE only for non-commercial and
                     16:  * non-profit purposes only, such as academic, research and internal
                     17:  * business use.
                     18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
                     19:  * international copyright treaties. If you make copies of the SOFTWARE,
                     20:  * with or without modification, as permitted hereunder, you shall affix
                     21:  * to all such copies of the SOFTWARE the above copyright notice.
                     22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
                     23:  * shall be made on your publication or presentation in any form of the
                     24:  * results obtained by use of the SOFTWARE.
                     25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3       noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2       noro       27:  * for such modification or the source code of the modified part of the
                     28:  * SOFTWARE.
                     29:  *
                     30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
                     31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
                     32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
                     33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
                     34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
                     35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
                     36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
                     37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
                     38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
                     39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
                     40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
                     41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
                     42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
                     43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
                     44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
                     45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
                     46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
                     47:  *
1.47    ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.46 2007/09/19 05:42:59 noro Exp $
1.2       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include "base.h"
1.16      noro       52: #include "inline.h"
1.1       noro       53: #include "parse.h"
                     54: #include "ox.h"
                     55:
1.5       noro       56: #define HMAG(p) (p_mag(BDY(p)->c))
                     57:
1.1       noro       58: extern int (*cmpdl)();
1.5       noro       59: extern double pz_t_e,pz_t_d,pz_t_d1,pz_t_c;
                     60: extern int dp_nelim,dp_fcoeffs;
1.7       noro       61: extern int NoGCD;
                     62: extern int GenTrace;
                     63: extern NODE TraceList;
                     64:
1.37      noro       65: int show_orderspec;
                     66:
                     67: void print_composite_order_spec(struct order_spec *spec);
                     68:
1.7       noro       69: /*
                     70:  * content reduction
                     71:  *
                     72:  */
                     73:
1.46      noro       74: static NODE RatDenomList;
                     75:
                     76: void init_denomlist()
                     77: {
                     78:        RatDenomList = 0;
                     79: }
                     80:
                     81: void add_denomlist(P f)
                     82: {
                     83:        NODE n;
                     84:
                     85:        if ( OID(f)==O_P ) {
                     86:                MKNODE(n,f,RatDenomList); RatDenomList = n;
                     87:        }
                     88: }
                     89:
                     90: LIST get_denomlist()
                     91: {
                     92:        LIST l;
                     93:
                     94:        MKLIST(l,RatDenomList); RatDenomList = 0;
                     95:        return l;
                     96: }
                     97:
1.20      noro       98: void dp_ptozp(DP p,DP *rp)
1.7       noro       99: {
                    100:        MP m,mr,mr0;
                    101:        int i,n;
                    102:        Q *w;
                    103:        Q dvr;
                    104:        P t;
                    105:
                    106:        if ( !p )
                    107:                *rp = 0;
                    108:        else {
                    109:                for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ );
                    110:                w = (Q *)ALLOCA(n*sizeof(Q));
                    111:                for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
                    112:                        if ( NUM(m->c) )
                    113:                                w[i] = (Q)m->c;
                    114:                        else
                    115:                                ptozp(m->c,1,&w[i],&t);
                    116:                sortbynm(w,n);
                    117:                qltozl(w,n,&dvr);
                    118:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    119:                        NEXTMP(mr0,mr); divsp(CO,m->c,(P)dvr,&mr->c); mr->dl = m->dl;
                    120:                }
                    121:                NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                    122:        }
                    123: }
                    124:
1.20      noro      125: void dp_ptozp2(DP p0,DP p1,DP *hp,DP *rp)
1.7       noro      126: {
                    127:        DP t,s,h,r;
                    128:        MP m,mr,mr0,m0;
                    129:
                    130:        addd(CO,p0,p1,&t); dp_ptozp(t,&s);
                    131:        if ( !p0 ) {
                    132:                h = 0; r = s;
                    133:        } else if ( !p1 ) {
                    134:                h = s; r = 0;
                    135:        } else {
                    136:                for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
                    137:                        m = NEXT(m), m0 = NEXT(m0) ) {
                    138:                        NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
                    139:                }
                    140:                NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
                    141:        }
                    142:        if ( h )
                    143:                h->sugar = p0->sugar;
                    144:        if ( r )
                    145:                r->sugar = p1->sugar;
                    146:        *hp = h; *rp = r;
1.39      ohara     147: }
                    148:
                    149: void dp_ptozp3(DP p,Q *dvr,DP *rp)
                    150: {
                    151:        MP m,mr,mr0;
                    152:        int i,n;
                    153:        Q *w;
                    154:        P t;
                    155:
                    156:        if ( !p ) {
                    157:                *rp = 0; *dvr = 0;
                    158:        }else {
                    159:                for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ );
                    160:                w = (Q *)ALLOCA(n*sizeof(Q));
                    161:                for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
                    162:                        if ( NUM(m->c) )
                    163:                                w[i] = (Q)m->c;
                    164:                        else
                    165:                                ptozp(m->c,1,&w[i],&t);
                    166:                sortbynm(w,n);
                    167:                qltozl(w,n,dvr);
                    168:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    169:                        NEXTMP(mr0,mr); divsp(CO,m->c,(P)(*dvr),&mr->c); mr->dl = m->dl;
                    170:                }
                    171:                NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                    172:        }
1.7       noro      173: }
1.1       noro      174:
1.20      noro      175: void dp_idiv(DP p,Q c,DP *rp)
1.1       noro      176: {
                    177:        Q t;
                    178:        N nm,q;
                    179:        int sgn,s;
                    180:        MP mr0,m,mr;
                    181:
                    182:        if ( !p )
                    183:                *rp = 0;
                    184:        else if ( MUNIQ((Q)c) )
                    185:                *rp = p;
                    186:        else if ( MUNIQ((Q)c) )
                    187:                chsgnd(p,rp);
                    188:        else {
                    189:                nm = NM(c); sgn = SGN(c);
                    190:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    191:                        NEXTMP(mr0,mr);
                    192:
                    193:                        divsn(NM((Q)(m->c)),nm,&q);
                    194:                        s = sgn*SGN((Q)(m->c));
                    195:                        NTOQ(q,s,t);
                    196:                        mr->c = (P)t;
                    197:                        mr->dl = m->dl;
                    198:                }
                    199:                NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
                    200:                if ( *rp )
                    201:                        (*rp)->sugar = p->sugar;
                    202:        }
                    203: }
                    204:
1.20      noro      205: void dp_mbase(NODE hlist,NODE *mbase)
1.1       noro      206: {
                    207:        DL *dl;
                    208:        DL d;
                    209:        int i,j,n,nvar,td;
                    210:
                    211:        n = length(hlist); nvar = ((DP)BDY(hlist))->nv;
                    212:        dl = (DL *)MALLOC(n*sizeof(DL));
                    213:        for ( i = 0; i < n; i++, hlist = NEXT(hlist) )
                    214:                dl[i] = BDY((DP)BDY(hlist))->dl;
                    215:        NEWDL(d,nvar); *mbase = 0;
                    216:        while ( 1 ) {
                    217:                insert_to_node(d,mbase,nvar);
                    218:                for ( i = nvar-1; i >= 0; ) {
1.21      noro      219:                        d->d[i]++;
                    220:                        d->td += MUL_WEIGHT(1,i);
1.1       noro      221:                        for ( j = 0; j < n; j++ ) {
                    222:                                if ( _dl_redble(dl[j],d,nvar) )
                    223:                                        break;
                    224:                        }
                    225:                        if ( j < n ) {
                    226:                                for ( j = nvar-1; j >= i; j-- )
                    227:                                        d->d[j] = 0;
                    228:                                for ( j = 0, td = 0; j < i; j++ )
1.21      noro      229:                                        td += MUL_WEIGHT(d->d[j],j);
1.1       noro      230:                                d->td = td;
                    231:                                i--;
                    232:                        } else
                    233:                                break;
                    234:                }
                    235:                if ( i < 0 )
                    236:                        break;
                    237:        }
                    238: }
                    239:
1.20      noro      240: int _dl_redble(DL d1,DL d2,int nvar)
1.1       noro      241: {
                    242:        int i;
                    243:
                    244:        if ( d1->td > d2->td )
                    245:                return 0;
                    246:        for ( i = 0; i < nvar; i++ )
                    247:                if ( d1->d[i] > d2->d[i] )
                    248:                        break;
                    249:        if ( i < nvar )
                    250:                return 0;
                    251:        else
                    252:                return 1;
                    253: }
                    254:
1.20      noro      255: void insert_to_node(DL d,NODE *n,int nvar)
1.1       noro      256: {
                    257:        DL d1;
                    258:        MP m;
                    259:        DP dp;
                    260:        NODE n0,n1,n2;
                    261:
                    262:        NEWDL(d1,nvar); d1->td = d->td;
                    263:        bcopy((char *)d->d,(char *)d1->d,nvar*sizeof(int));
                    264:        NEWMP(m); m->dl = d1; m->c = (P)ONE; NEXT(m) = 0;
                    265:        MKDP(nvar,m,dp); dp->sugar = d->td;
                    266:        if ( !(*n) ) {
                    267:                MKNODE(n1,dp,0); *n = n1;
                    268:        } else {
                    269:                for ( n1 = *n, n0 = 0; n1; n0 = n1, n1 = NEXT(n1) )
                    270:                        if ( (*cmpdl)(nvar,d,BDY((DP)BDY(n1))->dl) > 0 ) {
                    271:                                MKNODE(n2,dp,n1);
                    272:                                if ( !n0 )
                    273:                                        *n = n2;
                    274:                                else
                    275:                                        NEXT(n0) = n2;
                    276:                                break;
                    277:                        }
                    278:                if ( !n1 ) {
                    279:                        MKNODE(n2,dp,0); NEXT(n0) = n2;
                    280:                }
                    281:        }
                    282: }
                    283:
1.20      noro      284: void dp_vtod(Q *c,DP p,DP *rp)
1.1       noro      285: {
                    286:        MP mr0,m,mr;
                    287:        int i;
                    288:
                    289:        if ( !p )
                    290:                *rp = 0;
                    291:        else {
                    292:                for ( mr0 = 0, m = BDY(p), i = 0; m; m = NEXT(m), i++ ) {
                    293:                        NEXTMP(mr0,mr); mr->c = (P)c[i]; mr->dl = m->dl;
                    294:                }
                    295:                NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
                    296:                (*rp)->sugar = p->sugar;
                    297:        }
                    298: }
                    299:
1.8       noro      300: extern int mpi_mag;
                    301: extern int PCoeffs;
                    302:
1.20      noro      303: void dp_ptozp_d(DP p,DP *rp)
1.1       noro      304: {
                    305:        int i,j,k,l,n,nsep;
                    306:        MP m;
                    307:        NODE tn,n0,n1,n2,n3;
                    308:        struct oVECT v;
                    309:        VECT c,cs;
                    310:        VECT qi,ri;
                    311:        LIST *qr;
                    312:        Obj dmy;
                    313:        Q d0,d1,gcd,a,u,u1;
                    314:        Q *q,*r;
                    315:        STRING iqr_v;
                    316:        pointer *b;
                    317:        N qn,gn;
                    318:        double get_rtime();
                    319:        int blen;
1.8       noro      320:        NODE dist;
                    321:        int ndist;
1.1       noro      322:        double t0;
                    323:        double t_e,t_d,t_d1,t_c;
1.8       noro      324:        extern int DP_NFStat;
                    325:        extern LIST Dist;
1.20      noro      326:        void Pox_rpc();
                    327:        void Pox_pop_local();
1.1       noro      328:
                    329:        if ( !p )
                    330:                *rp = 0;
                    331:        else {
1.8       noro      332:                if ( PCoeffs ) {
                    333:                        dp_ptozp(p,rp); return;
                    334:                }
1.9       noro      335:                if ( !Dist || p_mag(BDY(p)->c) <= mpi_mag ) {
1.8       noro      336:                        dist = 0; ndist = 0;
                    337:                        if ( DP_NFStat ) fprintf(asir_out,"L");
                    338:                } else {
                    339:                        dist = BDY(Dist); ndist = length(dist);
                    340:                        if ( DP_NFStat ) fprintf(asir_out,"D");
                    341:                }
1.1       noro      342:                for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
                    343:                nsep = ndist + 1;
                    344:                if ( n <= nsep ) {
                    345:                        dp_ptozp(p,rp); return;
                    346:                }
                    347:                t0 = get_rtime();
                    348:                dp_dtov(p,&c);
                    349:                igcdv_estimate(c,&d0);
                    350:                t_e = get_rtime()-t0;
                    351:                t0 = get_rtime();
                    352:                dp_dtov(p,&c);
                    353:                sepvect(c,nsep,&cs);
                    354:                MKSTR(iqr_v,"iqr");
                    355:                qr = (LIST *)CALLOC(nsep,sizeof(LIST));
                    356:                q = (Q *)CALLOC(n,sizeof(Q));
                    357:                r = (Q *)CALLOC(n,sizeof(Q));
                    358:                for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
                    359:                        MKNODE(n3,d0,0); MKNODE(n2,b[i],n3);
                    360:                        MKNODE(n1,iqr_v,n2); MKNODE(n0,BDY(tn),n1);
                    361:                        Pox_rpc(n0,&dmy);
                    362:                }
                    363:                iqrv(b[i],d0,&qr[i]);
                    364:                dp_dtov(p,&c);
                    365:                for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
                    366:                        Pox_pop_local(tn,&qr[i]);
                    367:                        if ( OID(qr[i]) == O_ERR ) {
                    368:                                printexpr(CO,(Obj)qr[i]);
                    369:                                error("dp_ptozp_d : aborted");
                    370:                        }
                    371:                }
                    372:                t_d = get_rtime()-t0;
                    373:                t_d1 = t_d/n;
                    374:                t0 = get_rtime();
                    375:                for ( i = j = 0; i < nsep; i++ ) {
                    376:                        tn = BDY(qr[i]); qi = (VECT)BDY(tn); ri = (VECT)BDY(NEXT(tn));
                    377:                        for ( k = 0, l = qi->len; k < l; k++, j++ ) {
                    378:                                q[j] = (Q)BDY(qi)[k]; r[j] = (Q)BDY(ri)[k];
                    379:                        }
                    380:                }
                    381:                v.id = O_VECT; v.len = n; v.body = (pointer *)r; igcdv(&v,&d1);
                    382:                if ( d1 ) {
                    383:                        gcdn(NM(d0),NM(d1),&gn); NTOQ(gn,1,gcd);
                    384:                        divsn(NM(d0),gn,&qn); NTOQ(qn,1,a);
                    385:                        for ( i = 0; i < n; i++ ) {
                    386:                                mulq(a,q[i],&u);
                    387:                                if ( r[i] ) {
                    388:                                        divsn(NM(r[i]),gn,&qn); NTOQ(qn,SGN(r[i]),u1);
                    389:                                        addq(u,u1,&q[i]);
                    390:                                } else
                    391:                                        q[i] = u;
                    392:                        }
                    393:                } else
                    394:                        gcd = d0;
                    395:                dp_vtod(q,p,rp);
                    396:                t_c = get_rtime()-t0;
                    397:                blen=p_mag((P)gcd);
                    398:                pz_t_e += t_e; pz_t_d += t_d; pz_t_d1 += t_d1; pz_t_c += t_c;
                    399:                if ( 0 )
                    400:                        fprintf(stderr,"(%d,%d)",p_mag((P)d0)-blen,blen);
                    401:        }
                    402: }
                    403:
1.20      noro      404: void dp_ptozp2_d(DP p0,DP p1,DP *hp,DP *rp)
1.1       noro      405: {
                    406:        DP t,s,h,r;
                    407:        MP m,mr,mr0,m0;
                    408:
1.8       noro      409:        addd(CO,p0,p1,&t); dp_ptozp_d(t,&s);
1.1       noro      410:        if ( !p0 ) {
                    411:                h = 0; r = s;
                    412:        } else if ( !p1 ) {
                    413:                h = s; r = 0;
                    414:        } else {
                    415:                for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
                    416:                        m = NEXT(m), m0 = NEXT(m0) ) {
                    417:                        NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
                    418:                }
                    419:                NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
                    420:        }
                    421:        if ( h )
                    422:                h->sugar = p0->sugar;
                    423:        if ( r )
                    424:                r->sugar = p1->sugar;
                    425:        *hp = h; *rp = r;
1.5       noro      426: }
                    427:
1.22      noro      428: int have_sf_coef(P p)
                    429: {
                    430:        DCP dc;
                    431:
                    432:        if ( !p )
                    433:                return 0;
                    434:        else if ( NUM(p) )
                    435:                return NID((Num)p) == N_GFS ? 1 : 0;
                    436:        else {
                    437:                for ( dc = DC(p); dc; dc = NEXT(dc) )
                    438:                        if ( have_sf_coef(COEF(dc)) )
                    439:                                return 1;
                    440:                return 0;
                    441:        }
                    442: }
                    443:
1.25      noro      444: void head_coef(P p,Num *c)
                    445: {
                    446:        if ( !p )
                    447:                *c = 0;
                    448:        else if ( NUM(p) )
                    449:                *c = (Num)p;
                    450:        else
                    451:                head_coef(COEF(DC(p)),c);
                    452: }
                    453:
                    454: void dp_monic_sf(DP p,DP *rp)
                    455: {
                    456:        Num c;
                    457:
                    458:        if ( !p )
                    459:                *rp = 0;
                    460:        else {
                    461:                head_coef(BDY(p)->c,&c);
                    462:                divsdc(CO,p,(P)c,rp);
                    463:        }
                    464: }
                    465:
1.20      noro      466: void dp_prim(DP p,DP *rp)
1.5       noro      467: {
1.7       noro      468:        P t,g;
                    469:        DP p1;
                    470:        MP m,mr,mr0;
                    471:        int i,n;
                    472:        P *w;
                    473:        Q *c;
                    474:        Q dvr;
1.46      noro      475:        NODE tn;
1.5       noro      476:
1.7       noro      477:        if ( !p )
                    478:                *rp = 0;
1.23      noro      479:        else if ( dp_fcoeffs == N_GFS ) {
                    480:                for ( m = BDY(p); m; m = NEXT(m) )
1.22      noro      481:                        if ( OID(m->c) == O_N ) {
                    482:                                /* GCD of coeffs = 1 */
1.25      noro      483:                                dp_monic_sf(p,rp);
1.22      noro      484:                                return;
1.23      noro      485:                        } else break;
                    486:                /* compute GCD over the finite fieid */
                    487:                for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
                    488:                w = (P *)ALLOCA(n*sizeof(P));
                    489:                for ( m = BDY(p), i = 0; i < n; m = NEXT(m), i++ )
                    490:                        w[i] = m->c;
                    491:                gcdsf(CO,w,n,&g);
                    492:                if ( NUM(g) )
1.25      noro      493:                        dp_monic_sf(p,rp);
1.23      noro      494:                else {
                    495:                        for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    496:                                NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl;
1.22      noro      497:                        }
1.25      noro      498:                        NEXT(mr) = 0; MKDP(p->nv,mr0,p1); p1->sugar = p->sugar;
                    499:                        dp_monic_sf(p1,rp);
1.22      noro      500:                }
1.23      noro      501:                return;
                    502:        } else if ( dp_fcoeffs )
1.7       noro      503:                *rp = p;
1.23      noro      504:        else if ( NoGCD )
1.7       noro      505:                dp_ptozp(p,rp);
                    506:        else {
                    507:                dp_ptozp(p,&p1); p = p1;
                    508:                for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
                    509:                if ( n == 1 ) {
                    510:                        m = BDY(p);
                    511:                        NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
                    512:                        MKDP(p->nv,mr,*rp); (*rp)->sugar = p->sugar;
                    513:                        return;
                    514:                }
                    515:                w = (P *)ALLOCA(n*sizeof(P));
                    516:                c = (Q *)ALLOCA(n*sizeof(Q));
                    517:                for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
                    518:                        if ( NUM(m->c) ) {
                    519:                                c[i] = (Q)m->c; w[i] = (P)ONE;
                    520:                        } else
                    521:                                ptozp(m->c,1,&c[i],&w[i]);
                    522:                qltozl(c,n,&dvr); heu_nezgcdnpz(CO,w,n,&t); mulp(CO,t,(P)dvr,&g);
                    523:                if ( NUM(g) )
                    524:                        *rp = p;
                    525:                else {
                    526:                        for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    527:                                NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl;
                    528:                        }
                    529:                        NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.46      noro      530:                        add_denomlist(g);
1.5       noro      531:                }
1.7       noro      532:        }
1.5       noro      533: }
                    534:
1.20      noro      535: void heu_nezgcdnpz(VL vl,P *pl,int m,P *pr)
1.5       noro      536: {
                    537:        int i,r;
                    538:        P gcd,t,s1,s2,u;
                    539:        Q rq;
1.40      noro      540:        DCP dc;
                    541:        extern int DP_Print;
                    542:
1.5       noro      543:        while ( 1 ) {
                    544:                for ( i = 0, s1 = 0; i < m; i++ ) {
                    545:                        r = random(); UTOQ(r,rq);
                    546:                        mulp(vl,pl[i],(P)rq,&t); addp(vl,s1,t,&u); s1 = u;
                    547:                }
                    548:                for ( i = 0, s2 = 0; i < m; i++ ) {
                    549:                        r = random(); UTOQ(r,rq);
                    550:                        mulp(vl,pl[i],(P)rq,&t); addp(vl,s2,t,&u); s2 = u;
                    551:                }
                    552:                ezgcdp(vl,s1,s2,&gcd);
1.40      noro      553:                if ( DP_Print > 2 )
                    554:                        { fprintf(asir_out,"(%d)",nmonop(gcd)); fflush(asir_out); }
1.5       noro      555:                for ( i = 0; i < m; i++ ) {
                    556:                        if ( !divtpz(vl,pl[i],gcd,&t) )
                    557:                                break;
                    558:                }
                    559:                if ( i == m )
                    560:                        break;
                    561:        }
                    562:        *pr = gcd;
                    563: }
                    564:
1.20      noro      565: void dp_prim_mod(DP p,int mod,DP *rp)
1.5       noro      566: {
                    567:        P t,g;
                    568:        MP m,mr,mr0;
                    569:
                    570:        if ( !p )
                    571:                *rp = 0;
                    572:        else if ( NoGCD )
                    573:                *rp = p;
                    574:        else {
                    575:                for ( m = BDY(p), g = m->c, m = NEXT(m); m; m = NEXT(m) ) {
                    576:                        gcdprsmp(CO,mod,g,m->c,&t); g = t;
                    577:                }
                    578:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    579:                        NEXTMP(mr0,mr); divsmp(CO,mod,m->c,g,&mr->c); mr->dl = m->dl;
                    580:                }
                    581:                NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                    582:        }
                    583: }
                    584:
1.20      noro      585: void dp_cont(DP p,Q *rp)
1.5       noro      586: {
1.7       noro      587:        VECT v;
1.5       noro      588:
1.7       noro      589:        dp_dtov(p,&v); igcdv(v,rp);
1.5       noro      590: }
                    591:
1.20      noro      592: void dp_dtov(DP dp,VECT *rp)
1.5       noro      593: {
1.7       noro      594:        MP m,t;
                    595:        int i,n;
                    596:        VECT v;
                    597:        pointer *p;
1.5       noro      598:
1.7       noro      599:        m = BDY(dp);
                    600:        for ( t = m, n = 0; t; t = NEXT(t), n++ );
                    601:        MKVECT(v,n);
                    602:        for ( i = 0, p = BDY(v), t = m; i < n; t = NEXT(t), i++ )
                    603:                p[i] = (pointer)(t->c);
                    604:        *rp = v;
1.5       noro      605: }
                    606:
1.7       noro      607: /*
                    608:  * s-poly computation
                    609:  *
                    610:  */
1.5       noro      611:
1.20      noro      612: void dp_sp(DP p1,DP p2,DP *rp)
1.5       noro      613: {
1.7       noro      614:        int i,n,td;
                    615:        int *w;
                    616:        DL d1,d2,d;
                    617:        MP m;
                    618:        DP t,s1,s2,u;
                    619:        Q c,c1,c2;
                    620:        N gn,tn;
1.5       noro      621:
1.7       noro      622:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    623:        w = (int *)ALLOCA(n*sizeof(int));
                    624:        for ( i = 0, td = 0; i < n; i++ ) {
1.21      noro      625:                w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.5       noro      626:        }
1.7       noro      627:
                    628:        NEWDL(d,n); d->td = td - d1->td;
                    629:        for ( i = 0; i < n; i++ )
                    630:                d->d[i] = w[i] - d1->d[i];
                    631:        c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
                    632:        if ( INT(c1) && INT(c2) ) {
                    633:                gcdn(NM(c1),NM(c2),&gn);
                    634:                if ( !UNIN(gn) ) {
                    635:                        divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
                    636:                        divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
1.5       noro      637:                }
                    638:        }
1.7       noro      639:
                    640:        NEWMP(m); m->dl = d; m->c = (P)c2; NEXT(m) = 0;
                    641:        MKDP(n,m,s1); s1->sugar = d->td; muld(CO,s1,p1,&t);
                    642:
                    643:        NEWDL(d,n); d->td = td - d2->td;
                    644:        for ( i = 0; i < n; i++ )
                    645:                d->d[i] = w[i] - d2->d[i];
                    646:        NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0;
                    647:        MKDP(n,m,s2); s2->sugar = d->td; muld(CO,s2,p2,&u);
                    648:
                    649:        subd(CO,t,u,rp);
1.14      noro      650:        if ( GenTrace ) {
                    651:                LIST hist;
                    652:                NODE node;
                    653:
                    654:                node = mknode(4,ONE,0,s1,ONE);
                    655:                MKLIST(hist,node);
                    656:                MKNODE(TraceList,hist,0);
                    657:
                    658:                node = mknode(4,ONE,0,0,ONE);
                    659:                chsgnd(s2,(DP *)&ARG2(node));
                    660:                MKLIST(hist,node);
                    661:                MKNODE(node,hist,TraceList); TraceList = node;
                    662:        }
                    663: }
                    664:
1.20      noro      665: void _dp_sp_dup(DP p1,DP p2,DP *rp)
1.14      noro      666: {
                    667:        int i,n,td;
                    668:        int *w;
                    669:        DL d1,d2,d;
                    670:        MP m;
                    671:        DP t,s1,s2,u;
                    672:        Q c,c1,c2;
                    673:        N gn,tn;
                    674:
                    675:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    676:        w = (int *)ALLOCA(n*sizeof(int));
                    677:        for ( i = 0, td = 0; i < n; i++ ) {
1.21      noro      678:                w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.14      noro      679:        }
                    680:
                    681:        _NEWDL(d,n); d->td = td - d1->td;
                    682:        for ( i = 0; i < n; i++ )
                    683:                d->d[i] = w[i] - d1->d[i];
                    684:        c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
                    685:        if ( INT(c1) && INT(c2) ) {
                    686:                gcdn(NM(c1),NM(c2),&gn);
                    687:                if ( !UNIN(gn) ) {
                    688:                        divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
                    689:                        divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
                    690:                }
                    691:        }
                    692:
                    693:        _NEWMP(m); m->dl = d; m->c = (P)c2; NEXT(m) = 0;
                    694:        _MKDP(n,m,s1); s1->sugar = d->td; _muld_dup(CO,s1,p1,&t); _free_dp(s1);
                    695:
                    696:        _NEWDL(d,n); d->td = td - d2->td;
                    697:        for ( i = 0; i < n; i++ )
                    698:                d->d[i] = w[i] - d2->d[i];
                    699:        _NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0;
                    700:        _MKDP(n,m,s2); s2->sugar = d->td; _muld_dup(CO,s2,p2,&u); _free_dp(s2);
                    701:
                    702:        _addd_destructive(CO,t,u,rp);
1.7       noro      703:        if ( GenTrace ) {
                    704:                LIST hist;
                    705:                NODE node;
                    706:
                    707:                node = mknode(4,ONE,0,s1,ONE);
                    708:                MKLIST(hist,node);
                    709:                MKNODE(TraceList,hist,0);
                    710:
                    711:                node = mknode(4,ONE,0,0,ONE);
                    712:                chsgnd(s2,(DP *)&ARG2(node));
                    713:                MKLIST(hist,node);
                    714:                MKNODE(node,hist,TraceList); TraceList = node;
                    715:        }
                    716: }
                    717:
1.20      noro      718: void dp_sp_mod(DP p1,DP p2,int mod,DP *rp)
1.7       noro      719: {
                    720:        int i,n,td;
                    721:        int *w;
                    722:        DL d1,d2,d;
                    723:        MP m;
                    724:        DP t,s,u;
                    725:
                    726:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    727:        w = (int *)ALLOCA(n*sizeof(int));
                    728:        for ( i = 0, td = 0; i < n; i++ ) {
1.21      noro      729:                w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.7       noro      730:        }
1.18      noro      731:        NEWDL_NOINIT(d,n); d->td = td - d1->td;
1.7       noro      732:        for ( i = 0; i < n; i++ )
                    733:                d->d[i] = w[i] - d1->d[i];
                    734:        NEWMP(m); m->dl = d; m->c = (P)BDY(p2)->c; NEXT(m) = 0;
                    735:        MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p1,s,&t);
1.18      noro      736:        NEWDL_NOINIT(d,n); d->td = td - d2->td;
1.7       noro      737:        for ( i = 0; i < n; i++ )
                    738:                d->d[i] = w[i] - d2->d[i];
                    739:        NEWMP(m); m->dl = d; m->c = (P)BDY(p1)->c; NEXT(m) = 0;
                    740:        MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p2,s,&u);
                    741:        submd(CO,mod,t,u,rp);
                    742: }
                    743:
1.20      noro      744: void _dp_sp_mod_dup(DP p1,DP p2,int mod,DP *rp)
1.7       noro      745: {
                    746:        int i,n,td;
                    747:        int *w;
                    748:        DL d1,d2,d;
                    749:        MP m;
                    750:        DP t,s,u;
                    751:
                    752:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    753:        w = (int *)ALLOCA(n*sizeof(int));
                    754:        for ( i = 0, td = 0; i < n; i++ ) {
1.21      noro      755:                w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.7       noro      756:        }
                    757:        _NEWDL(d,n); d->td = td - d1->td;
                    758:        for ( i = 0; i < n; i++ )
                    759:                d->d[i] = w[i] - d1->d[i];
                    760:        _NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
                    761:        _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p1,&t); _free_dp(s);
                    762:        _NEWDL(d,n); d->td = td - d2->td;
                    763:        for ( i = 0; i < n; i++ )
                    764:                d->d[i] = w[i] - d2->d[i];
                    765:        _NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
                    766:        _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p2,&u); _free_dp(s);
                    767:        _addmd_destructive(mod,t,u,rp);
                    768: }
                    769:
1.20      noro      770: void _dp_sp_mod(DP p1,DP p2,int mod,DP *rp)
1.7       noro      771: {
                    772:        int i,n,td;
                    773:        int *w;
                    774:        DL d1,d2,d;
                    775:        MP m;
                    776:        DP t,s,u;
                    777:
                    778:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    779:        w = (int *)ALLOCA(n*sizeof(int));
                    780:        for ( i = 0, td = 0; i < n; i++ ) {
1.21      noro      781:                w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.7       noro      782:        }
                    783:        NEWDL(d,n); d->td = td - d1->td;
                    784:        for ( i = 0; i < n; i++ )
                    785:                d->d[i] = w[i] - d1->d[i];
                    786:        NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
                    787:        MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p1,&t);
                    788:        NEWDL(d,n); d->td = td - d2->td;
                    789:        for ( i = 0; i < n; i++ )
                    790:                d->d[i] = w[i] - d2->d[i];
                    791:        NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
                    792:        MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p2,&u);
                    793:        addmd_destructive(mod,t,u,rp);
                    794: }
                    795:
                    796: /*
                    797:  * m-reduction
1.13      noro      798:  * do content reduction over Z or Q(x,...)
                    799:  * do nothing over finite fields
1.7       noro      800:  *
                    801:  */
                    802:
1.20      noro      803: void dp_red(DP p0,DP p1,DP p2,DP *head,DP *rest,P *dnp,DP *multp)
1.7       noro      804: {
                    805:        int i,n;
                    806:        DL d1,d2,d;
                    807:        MP m;
                    808:        DP t,s,r,h;
                    809:        Q c,c1,c2;
                    810:        N gn,tn;
                    811:        P g,a;
1.23      noro      812:        P p[2];
1.7       noro      813:
                    814:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    815:        NEWDL(d,n); d->td = d1->td - d2->td;
                    816:        for ( i = 0; i < n; i++ )
                    817:                d->d[i] = d1->d[i]-d2->d[i];
                    818:        c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
1.23      noro      819:        if ( dp_fcoeffs == N_GFS ) {
                    820:                p[0] = (P)c1; p[1] = (P)c2;
                    821:                gcdsf(CO,p,2,&g);
                    822:                divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
                    823:        } else if ( dp_fcoeffs ) {
1.7       noro      824:                /* do nothing */
                    825:        } else if ( INT(c1) && INT(c2) ) {
                    826:                gcdn(NM(c1),NM(c2),&gn);
                    827:                if ( !UNIN(gn) ) {
                    828:                        divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
                    829:                        divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
                    830:                }
                    831:        } else {
                    832:                ezgcdpz(CO,(P)c1,(P)c2,&g);
                    833:                divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
1.46      noro      834:                add_denomlist(g);
1.7       noro      835:        }
                    836:        NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
                    837:        *multp = s;
                    838:        muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); addd(CO,s,t,&r);
                    839:        muldc(CO,p0,(P)c2,&h);
                    840:        *head = h; *rest = r; *dnp = (P)c2;
                    841: }
                    842:
1.41      noro      843: /*
                    844:  * m-reduction by a marked poly
                    845:  * do content reduction over Z or Q(x,...)
                    846:  * do nothing over finite fields
                    847:  *
                    848:  */
                    849:
                    850:
                    851: void dp_red_marked(DP p0,DP p1,DP p2,DP hp2,DP *head,DP *rest,P *dnp,DP *multp)
                    852: {
                    853:        int i,n;
                    854:        DL d1,d2,d;
                    855:        MP m;
                    856:        DP t,s,r,h;
                    857:        Q c,c1,c2;
                    858:        N gn,tn;
                    859:        P g,a;
                    860:        P p[2];
                    861:
                    862:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(hp2)->dl;
                    863:        NEWDL(d,n); d->td = d1->td - d2->td;
                    864:        for ( i = 0; i < n; i++ )
                    865:                d->d[i] = d1->d[i]-d2->d[i];
                    866:        c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(hp2)->c;
                    867:        if ( dp_fcoeffs == N_GFS ) {
                    868:                p[0] = (P)c1; p[1] = (P)c2;
                    869:                gcdsf(CO,p,2,&g);
                    870:                divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
                    871:        } else if ( dp_fcoeffs ) {
                    872:                /* do nothing */
                    873:        } else if ( INT(c1) && INT(c2) ) {
                    874:                gcdn(NM(c1),NM(c2),&gn);
                    875:                if ( !UNIN(gn) ) {
                    876:                        divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
                    877:                        divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
                    878:                }
                    879:        } else {
                    880:                ezgcdpz(CO,(P)c1,(P)c2,&g);
                    881:                divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
                    882:        }
1.47    ! noro      883:        NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
1.41      noro      884:        *multp = s;
1.47    ! noro      885:        muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); subd(CO,s,t,&r);
1.41      noro      886:        muldc(CO,p0,(P)c2,&h);
                    887:        *head = h; *rest = r; *dnp = (P)c2;
                    888: }
                    889:
1.44      noro      890: void dp_red_marked_mod(DP p0,DP p1,DP p2,DP hp2,int mod,DP *head,DP *rest,P *dnp)
                    891: {
                    892:        int i,n;
                    893:        DL d1,d2,d;
                    894:        MP m;
                    895:        DP t,s,r,h;
                    896:        P c1,c2,g,u;
                    897:
                    898:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(hp2)->dl;
                    899:        NEWDL(d,n); d->td = d1->td - d2->td;
                    900:        for ( i = 0; i < n; i++ )
                    901:                d->d[i] = d1->d[i]-d2->d[i];
                    902:        c1 = (P)BDY(p1)->c; c2 = (P)BDY(hp2)->c;
                    903:        gcdprsmp(CO,mod,c1,c2,&g);
                    904:        divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
                    905:        if ( NUM(c2) ) {
                    906:                divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
                    907:        }
                    908:        NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0;
                    909:        MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,s,p2,&t);
                    910:        if ( NUM(c2) ) {
                    911:                addmd(CO,mod,p1,t,&r); h = p0;
                    912:        } else {
                    913:                mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
                    914:        }
                    915:        *head = h; *rest = r; *dnp = c2;
                    916: }
                    917:
1.13      noro      918: /* m-reduction over a field */
                    919:
1.20      noro      920: void dp_red_f(DP p1,DP p2,DP *rest)
1.13      noro      921: {
                    922:        int i,n;
                    923:        DL d1,d2,d;
                    924:        MP m;
1.20      noro      925:        DP t,s;
1.13      noro      926:        Obj a,b;
                    927:
                    928:        n = p1->nv;
                    929:        d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    930:
                    931:        NEWDL(d,n); d->td = d1->td - d2->td;
                    932:        for ( i = 0; i < n; i++ )
                    933:                d->d[i] = d1->d[i]-d2->d[i];
                    934:
                    935:        NEWMP(m); m->dl = d;
                    936:        divr(CO,(Obj)BDY(p1)->c,(Obj)BDY(p2)->c,&a); chsgnr(a,&b);
                    937:        C(m) = (P)b;
                    938:        NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
                    939:
                    940:        muld(CO,s,p2,&t); addd(CO,p1,t,rest);
                    941: }
                    942:
1.20      noro      943: void dp_red_mod(DP p0,DP p1,DP p2,int mod,DP *head,DP *rest,P *dnp)
1.7       noro      944: {
                    945:        int i,n;
                    946:        DL d1,d2,d;
                    947:        MP m;
                    948:        DP t,s,r,h;
                    949:        P c1,c2,g,u;
                    950:
                    951:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    952:        NEWDL(d,n); d->td = d1->td - d2->td;
                    953:        for ( i = 0; i < n; i++ )
                    954:                d->d[i] = d1->d[i]-d2->d[i];
                    955:        c1 = (P)BDY(p1)->c; c2 = (P)BDY(p2)->c;
                    956:        gcdprsmp(CO,mod,c1,c2,&g);
                    957:        divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
                    958:        if ( NUM(c2) ) {
                    959:                divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
                    960:        }
                    961:        NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0;
1.11      noro      962:        MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,s,p2,&t);
1.7       noro      963:        if ( NUM(c2) ) {
                    964:                addmd(CO,mod,p1,t,&r); h = p0;
                    965:        } else {
                    966:                mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
                    967:        }
                    968:        *head = h; *rest = r; *dnp = c2;
                    969: }
                    970:
1.10      noro      971: struct oEGT eg_red_mod;
                    972:
1.20      noro      973: void _dp_red_mod_destructive(DP p1,DP p2,int mod,DP *rp)
1.7       noro      974: {
                    975:        int i,n;
                    976:        DL d1,d2,d;
                    977:        MP m;
                    978:        DP t,s;
1.16      noro      979:        int c,c1,c2;
                    980:        extern int do_weyl;
1.7       noro      981:
                    982:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    983:        _NEWDL(d,n); d->td = d1->td - d2->td;
                    984:        for ( i = 0; i < n; i++ )
                    985:                d->d[i] = d1->d[i]-d2->d[i];
1.16      noro      986:        c = invm(ITOS(BDY(p2)->c),mod);
                    987:        c2 = ITOS(BDY(p1)->c);
                    988:        DMAR(c,c2,0,mod,c1);
1.7       noro      989:        _NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
1.16      noro      990: #if 0
1.7       noro      991:        _MKDP(n,m,s); s->sugar = d->td;
                    992:        _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16      noro      993: #else
                    994:        if ( do_weyl ) {
1.19      noro      995:                _MKDP(n,m,s); s->sugar = d->td;
                    996:                _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16      noro      997:        } else {
                    998:                _mulmdm_dup(mod,p2,m,&t); _FREEMP(m);
                    999:        }
                   1000: #endif
1.10      noro     1001: /* get_eg(&t0); */
1.7       noro     1002:        _addmd_destructive(mod,p1,t,rp);
1.10      noro     1003: /* get_eg(&t1); add_eg(&eg_red_mod,&t0,&t1); */
1.7       noro     1004: }
                   1005:
                   1006: /*
                   1007:  * normal form computation
                   1008:  *
                   1009:  */
1.5       noro     1010:
1.20      noro     1011: void dp_true_nf(NODE b,DP g,DP *ps,int full,DP *rp,P *dnp)
1.5       noro     1012: {
                   1013:        DP u,p,d,s,t,dmy;
                   1014:        NODE l;
                   1015:        MP m,mr;
                   1016:        int i,n;
                   1017:        int *wb;
                   1018:        int sugar,psugar;
                   1019:        P dn,tdn,tdn1;
                   1020:
                   1021:        dn = (P)ONE;
                   1022:        if ( !g ) {
                   1023:                *rp = 0; *dnp = dn; return;
                   1024:        }
                   1025:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1026:        wb = (int *)ALLOCA(n*sizeof(int));
                   1027:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1028:                wb[i] = QTOS((Q)BDY(l));
                   1029:        sugar = g->sugar;
                   1030:        for ( d = 0; g; ) {
                   1031:                for ( u = 0, i = 0; i < n; i++ ) {
                   1032:                        if ( dp_redble(g,p = ps[wb[i]]) ) {
                   1033:                                dp_red(d,g,p,&t,&u,&tdn,&dmy);
                   1034:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1035:                                sugar = MAX(sugar,psugar);
                   1036:                                if ( !u ) {
                   1037:                                        if ( d )
                   1038:                                                d->sugar = sugar;
                   1039:                                        *rp = d; *dnp = dn; return;
                   1040:                                } else {
                   1041:                                        d = t;
                   1042:                                        mulp(CO,dn,tdn,&tdn1); dn = tdn1;
                   1043:                                }
                   1044:                                break;
                   1045:                        }
                   1046:                }
                   1047:                if ( u )
                   1048:                        g = u;
                   1049:                else if ( !full ) {
                   1050:                        if ( g ) {
                   1051:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                   1052:                        }
                   1053:                        *rp = g; *dnp = dn; return;
                   1054:                } else {
                   1055:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1056:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1057:                        addd(CO,d,t,&s); d = s;
                   1058:                        dp_rest(g,&t); g = t;
                   1059:                }
                   1060:        }
                   1061:        if ( d )
                   1062:                d->sugar = sugar;
                   1063:        *rp = d; *dnp = dn;
                   1064: }
                   1065:
1.43      noro     1066: void dp_removecont2(DP p1,DP p2,DP *r1p,DP *r2p,Q *contp)
                   1067: {
                   1068:        struct oVECT v;
                   1069:        int i,n1,n2,n;
                   1070:        MP m,m0,t;
                   1071:        Q *w;
                   1072:        Q h;
                   1073:
                   1074:        if ( p1 ) {
                   1075:                for ( i = 0, m = BDY(p1); m; m = NEXT(m), i++ );
                   1076:                n1 = i;
                   1077:        } else
                   1078:                n1 = 0;
                   1079:        if ( p2 ) {
                   1080:                for ( i = 0, m = BDY(p2); m; m = NEXT(m), i++ );
                   1081:                n2 = i;
                   1082:        } else
                   1083:                n2 = 0;
                   1084:        n = n1+n2;
                   1085:        if ( !n ) {
                   1086:                *r1p = 0; *r2p = 0; *contp = ONE; return;
                   1087:        }
                   1088:        w = (Q *)ALLOCA(n*sizeof(Q));
                   1089:        v.len = n;
                   1090:        v.body = (pointer *)w;
                   1091:        i = 0;
                   1092:        if ( p1 )
                   1093:                for ( m = BDY(p1); i < n1; m = NEXT(m), i++ ) w[i] = (Q)m->c;
                   1094:        if ( p2 )
                   1095:                for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) w[i] = (Q)m->c;
                   1096:        h = w[0]; removecont_array((P *)w,n,1); divq(h,w[0],contp);
                   1097:        i = 0;
                   1098:        if ( p1 ) {
                   1099:                for ( m0 = 0, t = BDY(p1); i < n1; i++, t = NEXT(t) ) {
                   1100:                        NEXTMP(m0,m); m->c = (P)w[i]; m->dl = t->dl;
                   1101:                }
                   1102:                NEXT(m) = 0;
                   1103:                MKDP(p1->nv,m0,*r1p); (*r1p)->sugar = p1->sugar;
                   1104:        } else
                   1105:                *r1p = 0;
                   1106:        if ( p2 ) {
                   1107:                for ( m0 = 0, t = BDY(p2); i < n; i++, t = NEXT(t) ) {
                   1108:                        NEXTMP(m0,m); m->c = (P)w[i]; m->dl = t->dl;
                   1109:                }
                   1110:                NEXT(m) = 0;
                   1111:                MKDP(p2->nv,m0,*r2p); (*r2p)->sugar = p2->sugar;
                   1112:        } else
                   1113:                *r2p = 0;
                   1114: }
                   1115:
1.41      noro     1116: /* true nf by a marked GB */
                   1117:
1.43      noro     1118: void dp_true_nf_marked(NODE b,DP g,DP *ps,DP *hps,DP *rp,P *nmp,P *dnp)
1.41      noro     1119: {
                   1120:        DP u,p,d,s,t,dmy,hp;
                   1121:        NODE l;
                   1122:        MP m,mr;
1.43      noro     1123:        int i,n,hmag;
1.41      noro     1124:        int *wb;
1.43      noro     1125:        int sugar,psugar,multiple;
                   1126:        P nm,tnm1,dn,tdn,tdn1;
                   1127:        Q cont;
1.41      noro     1128:
1.43      noro     1129:        multiple = 0;
                   1130:        hmag = multiple*HMAG(g);
                   1131:        nm = (P)ONE;
1.41      noro     1132:        dn = (P)ONE;
                   1133:        if ( !g ) {
                   1134:                *rp = 0; *dnp = dn; return;
                   1135:        }
                   1136:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1137:        wb = (int *)ALLOCA(n*sizeof(int));
                   1138:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1139:                wb[i] = QTOS((Q)BDY(l));
                   1140:        sugar = g->sugar;
                   1141:        for ( d = 0; g; ) {
                   1142:                for ( u = 0, i = 0; i < n; i++ ) {
                   1143:                        if ( dp_redble(g,hp = hps[wb[i]]) ) {
                   1144:                                p = ps[wb[i]];
                   1145:                                dp_red_marked(d,g,p,hp,&t,&u,&tdn,&dmy);
                   1146:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1147:                                sugar = MAX(sugar,psugar);
                   1148:                                if ( !u ) {
1.43      noro     1149:                                        goto last;
1.41      noro     1150:                                } else {
                   1151:                                        d = t;
                   1152:                                        mulp(CO,dn,tdn,&tdn1); dn = tdn1;
                   1153:                                }
                   1154:                                break;
                   1155:                        }
                   1156:                }
1.43      noro     1157:                if ( u ) {
1.41      noro     1158:                        g = u;
1.43      noro     1159:                        if ( multiple && ((d && HMAG(d)>hmag) || (HMAG(g)>hmag)) ) {
                   1160:                                dp_removecont2(d,g,&t,&u,&cont); d = t; g = u;
                   1161:                                mulp(CO,nm,(P)cont,&tnm1); nm = tnm1;
                   1162:                                if ( d )
                   1163:                                        hmag = multiple*HMAG(d);
                   1164:                                else
                   1165:                                        hmag = multiple*HMAG(g);
                   1166:                        }
                   1167:                } else {
1.41      noro     1168:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1169:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1170:                        addd(CO,d,t,&s); d = s;
                   1171:                        dp_rest(g,&t); g = t;
                   1172:                }
                   1173:        }
1.43      noro     1174: last:
                   1175:        if ( d ) {
                   1176:                dp_removecont2(d,0,&t,&u,&cont); d = t;
                   1177:                mulp(CO,nm,(P)cont,&tnm1); nm = tnm1;
1.41      noro     1178:                d->sugar = sugar;
1.43      noro     1179:        }
                   1180:        *rp = d; *nmp = nm; *dnp = dn;
1.41      noro     1181: }
                   1182:
1.44      noro     1183: void dp_true_nf_marked_mod(NODE b,DP g,DP *ps,DP *hps,int mod,DP *rp,P *dnp)
                   1184: {
                   1185:        DP hp,u,p,d,s,t;
                   1186:        NODE l;
                   1187:        MP m,mr;
                   1188:        int i,n;
                   1189:        int *wb;
                   1190:        int sugar,psugar;
                   1191:        P dn,tdn,tdn1;
                   1192:
                   1193:        dn = (P)ONEM;
                   1194:        if ( !g ) {
                   1195:                *rp = 0; *dnp = dn; return;
                   1196:        }
                   1197:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1198:                wb = (int *)ALLOCA(n*sizeof(int));
                   1199:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1200:                wb[i] = QTOS((Q)BDY(l));
                   1201:        sugar = g->sugar;
                   1202:        for ( d = 0; g; ) {
                   1203:                for ( u = 0, i = 0; i < n; i++ ) {
                   1204:                        if ( dp_redble(g,hp = hps[wb[i]]) ) {
                   1205:                                p = ps[wb[i]];
                   1206:                                dp_red_marked_mod(d,g,p,hp,mod,&t,&u,&tdn);
                   1207:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1208:                                sugar = MAX(sugar,psugar);
                   1209:                                if ( !u ) {
                   1210:                                        if ( d )
                   1211:                                                d->sugar = sugar;
                   1212:                                        *rp = d; *dnp = dn; return;
                   1213:                                } else {
                   1214:                                        d = t;
                   1215:                                        mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
                   1216:                                }
                   1217:                                break;
                   1218:                        }
                   1219:                }
                   1220:                if ( u )
                   1221:                        g = u;
                   1222:                else {
                   1223:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1224:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1225:                        addmd(CO,mod,d,t,&s); d = s;
                   1226:                        dp_rest(g,&t); g = t;
                   1227:                }
                   1228:        }
                   1229:        if ( d )
                   1230:                d->sugar = sugar;
                   1231:        *rp = d; *dnp = dn;
                   1232: }
                   1233:
1.47    ! noro     1234: /* true nf by a marked GB and collect quotients */
        !          1235:
        !          1236: DP *dp_true_nf_and_quotient_marked (NODE b,DP g,DP *ps,DP *hps,DP *rp,P *dnp)
        !          1237: {
        !          1238:        DP u,p,d,s,t,dmy,hp,mult;
        !          1239:        DP *q;
        !          1240:        NODE l;
        !          1241:        MP m,mr;
        !          1242:        int i,n,j;
        !          1243:        int *wb;
        !          1244:        int sugar,psugar,multiple;
        !          1245:        P nm,tnm1,dn,tdn,tdn1;
        !          1246:        Q cont;
        !          1247:
        !          1248:        dn = (P)ONE;
        !          1249:        if ( !g ) {
        !          1250:                *rp = 0; *dnp = dn; return;
        !          1251:        }
        !          1252:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
        !          1253:        wb = (int *)ALLOCA(n*sizeof(int));
        !          1254:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
        !          1255:                wb[i] = QTOS((Q)BDY(l));
        !          1256:        q = (DP *)MALLOC(n*sizeof(DP));
        !          1257:        for ( i = 0; i < n; i++ ) q[i] = 0;
        !          1258:        sugar = g->sugar;
        !          1259:        for ( d = 0; g; ) {
        !          1260:                for ( u = 0, i = 0; i < n; i++ ) {
        !          1261:                        if ( dp_redble(g,hp = hps[wb[i]]) ) {
        !          1262:                                p = ps[wb[i]];
        !          1263:                                dp_red_marked(d,g,p,hp,&t,&u,&tdn,&mult);
        !          1264:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
        !          1265:                                sugar = MAX(sugar,psugar);
        !          1266:                                for ( j = 0; j < n; j++ ) {
        !          1267:                                        muldc(CO,q[j],(P)tdn,&dmy); q[j] = dmy;
        !          1268:                                }
        !          1269:                                addd(CO,q[wb[i]],mult,&dmy); q[wb[i]] = dmy;
        !          1270:                                mulp(CO,dn,tdn,&tdn1); dn = tdn1;
        !          1271:                                d = t;
        !          1272:                                if ( !u ) goto last;
        !          1273:                                break;
        !          1274:                        }
        !          1275:                }
        !          1276:                if ( u ) {
        !          1277:                        g = u;
        !          1278:                } else {
        !          1279:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
        !          1280:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
        !          1281:                        addd(CO,d,t,&s); d = s;
        !          1282:                        dp_rest(g,&t); g = t;
        !          1283:                }
        !          1284:        }
        !          1285: last:
        !          1286:        if ( d ) d->sugar = sugar;
        !          1287:        *rp = d; *dnp = dn;
        !          1288:        return q;
        !          1289: }
        !          1290:
1.13      noro     1291: /* nf computation over Z */
                   1292:
1.20      noro     1293: void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *rp)
1.5       noro     1294: {
                   1295:        DP u,p,d,s,t,dmy1;
                   1296:        P dmy;
                   1297:        NODE l;
                   1298:        MP m,mr;
                   1299:        int i,n;
                   1300:        int *wb;
                   1301:        int hmag;
                   1302:        int sugar,psugar;
                   1303:
                   1304:        if ( !g ) {
                   1305:                *rp = 0; return;
                   1306:        }
                   1307:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1308:        wb = (int *)ALLOCA(n*sizeof(int));
                   1309:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1310:                wb[i] = QTOS((Q)BDY(l));
1.12      noro     1311:
1.13      noro     1312:        hmag = multiple*HMAG(g);
1.5       noro     1313:        sugar = g->sugar;
1.12      noro     1314:
1.5       noro     1315:        for ( d = 0; g; ) {
                   1316:                for ( u = 0, i = 0; i < n; i++ ) {
                   1317:                        if ( dp_redble(g,p = ps[wb[i]]) ) {
                   1318:                                dp_red(d,g,p,&t,&u,&dmy,&dmy1);
                   1319:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1320:                                sugar = MAX(sugar,psugar);
                   1321:                                if ( !u ) {
                   1322:                                        if ( d )
                   1323:                                                d->sugar = sugar;
                   1324:                                        *rp = d; return;
                   1325:                                }
                   1326:                                d = t;
                   1327:                                break;
                   1328:                        }
                   1329:                }
                   1330:                if ( u ) {
                   1331:                        g = u;
                   1332:                        if ( d ) {
1.13      noro     1333:                                if ( multiple && HMAG(d) > hmag ) {
1.5       noro     1334:                                        dp_ptozp2(d,g,&t,&u); d = t; g = u;
                   1335:                                        hmag = multiple*HMAG(d);
                   1336:                                }
                   1337:                        } else {
1.13      noro     1338:                                if ( multiple && HMAG(g) > hmag ) {
1.5       noro     1339:                                        dp_ptozp(g,&t); g = t;
                   1340:                                        hmag = multiple*HMAG(g);
                   1341:                                }
                   1342:                        }
                   1343:                }
                   1344:                else if ( !full ) {
                   1345:                        if ( g ) {
                   1346:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                   1347:                        }
                   1348:                        *rp = g; return;
                   1349:                } else {
                   1350:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1351:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1352:                        addd(CO,d,t,&s); d = s;
                   1353:                        dp_rest(g,&t); g = t;
                   1354:
                   1355:                }
                   1356:        }
                   1357:        if ( d )
                   1358:                d->sugar = sugar;
                   1359:        *rp = d;
                   1360: }
                   1361:
1.13      noro     1362: /* nf computation over a field */
                   1363:
1.20      noro     1364: void dp_nf_f(NODE b,DP g,DP *ps,int full,DP *rp)
1.13      noro     1365: {
                   1366:        DP u,p,d,s,t;
                   1367:        NODE l;
                   1368:        MP m,mr;
                   1369:        int i,n;
                   1370:        int *wb;
                   1371:        int sugar,psugar;
                   1372:
                   1373:        if ( !g ) {
                   1374:                *rp = 0; return;
                   1375:        }
                   1376:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1377:        wb = (int *)ALLOCA(n*sizeof(int));
                   1378:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1379:                wb[i] = QTOS((Q)BDY(l));
                   1380:
                   1381:        sugar = g->sugar;
                   1382:        for ( d = 0; g; ) {
                   1383:                for ( u = 0, i = 0; i < n; i++ ) {
                   1384:                        if ( dp_redble(g,p = ps[wb[i]]) ) {
                   1385:                                dp_red_f(g,p,&u);
                   1386:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1387:                                sugar = MAX(sugar,psugar);
                   1388:                                if ( !u ) {
                   1389:                                        if ( d )
                   1390:                                                d->sugar = sugar;
                   1391:                                        *rp = d; return;
                   1392:                                }
                   1393:                                break;
                   1394:                        }
                   1395:                }
                   1396:                if ( u )
                   1397:                        g = u;
                   1398:                else if ( !full ) {
                   1399:                        if ( g ) {
                   1400:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                   1401:                        }
                   1402:                        *rp = g; return;
                   1403:                } else {
                   1404:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1405:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1406:                        addd(CO,d,t,&s); d = s;
                   1407:                        dp_rest(g,&t); g = t;
                   1408:                }
                   1409:        }
                   1410:        if ( d )
                   1411:                d->sugar = sugar;
                   1412:        *rp = d;
                   1413: }
                   1414:
                   1415: /* nf computation over GF(mod) (only for internal use) */
                   1416:
1.20      noro     1417: void dp_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5       noro     1418: {
                   1419:        DP u,p,d,s,t;
                   1420:        P dmy;
                   1421:        NODE l;
                   1422:        MP m,mr;
                   1423:        int sugar,psugar;
                   1424:
                   1425:        if ( !g ) {
                   1426:                *rp = 0; return;
                   1427:        }
                   1428:        sugar = g->sugar;
                   1429:        for ( d = 0; g; ) {
                   1430:                for ( u = 0, l = b; l; l = NEXT(l) ) {
                   1431:                        if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
                   1432:                                dp_red_mod(d,g,p,mod,&t,&u,&dmy);
                   1433:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1434:                                sugar = MAX(sugar,psugar);
                   1435:                                if ( !u ) {
                   1436:                                        if ( d )
                   1437:                                                d->sugar = sugar;
                   1438:                                        *rp = d; return;
                   1439:                                }
                   1440:                                d = t;
                   1441:                                break;
                   1442:                        }
                   1443:                }
                   1444:                if ( u )
                   1445:                        g = u;
                   1446:                else if ( !full ) {
                   1447:                        if ( g ) {
                   1448:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                   1449:                        }
                   1450:                        *rp = g; return;
                   1451:                } else {
                   1452:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1453:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1454:                        addmd(CO,mod,d,t,&s); d = s;
                   1455:                        dp_rest(g,&t); g = t;
                   1456:                }
                   1457:        }
                   1458:        if ( d )
                   1459:                d->sugar = sugar;
                   1460:        *rp = d;
                   1461: }
                   1462:
1.20      noro     1463: void dp_true_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp,P *dnp)
1.5       noro     1464: {
                   1465:        DP u,p,d,s,t;
                   1466:        NODE l;
                   1467:        MP m,mr;
                   1468:        int i,n;
                   1469:        int *wb;
                   1470:        int sugar,psugar;
                   1471:        P dn,tdn,tdn1;
                   1472:
                   1473:        dn = (P)ONEM;
                   1474:        if ( !g ) {
                   1475:                *rp = 0; *dnp = dn; return;
                   1476:        }
                   1477:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1478:                wb = (int *)ALLOCA(n*sizeof(int));
                   1479:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1480:                wb[i] = QTOS((Q)BDY(l));
                   1481:        sugar = g->sugar;
                   1482:        for ( d = 0; g; ) {
                   1483:                for ( u = 0, i = 0; i < n; i++ ) {
                   1484:                        if ( dp_redble(g,p = ps[wb[i]]) ) {
                   1485:                                dp_red_mod(d,g,p,mod,&t,&u,&tdn);
                   1486:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1487:                                sugar = MAX(sugar,psugar);
                   1488:                                if ( !u ) {
                   1489:                                        if ( d )
                   1490:                                                d->sugar = sugar;
                   1491:                                        *rp = d; *dnp = dn; return;
                   1492:                                } else {
                   1493:                                        d = t;
                   1494:                                        mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
                   1495:                                }
                   1496:                                break;
                   1497:                        }
                   1498:                }
                   1499:                if ( u )
                   1500:                        g = u;
                   1501:                else if ( !full ) {
                   1502:                        if ( g ) {
                   1503:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                   1504:                        }
                   1505:                        *rp = g; *dnp = dn; return;
                   1506:                } else {
                   1507:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1508:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1509:                        addmd(CO,mod,d,t,&s); d = s;
                   1510:                        dp_rest(g,&t); g = t;
                   1511:                }
                   1512:        }
                   1513:        if ( d )
                   1514:                d->sugar = sugar;
                   1515:        *rp = d; *dnp = dn;
                   1516: }
                   1517:
1.20      noro     1518: void _dp_nf_mod_destructive(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5       noro     1519: {
1.20      noro     1520:        DP u,p,d;
1.7       noro     1521:        NODE l;
1.20      noro     1522:        MP m,mrd;
                   1523:        int sugar,psugar,n,h_reducible;
1.5       noro     1524:
1.7       noro     1525:        if ( !g ) {
                   1526:                *rp = 0; return;
1.5       noro     1527:        }
1.7       noro     1528:        sugar = g->sugar;
                   1529:        n = g->nv;
                   1530:        for ( d = 0; g; ) {
                   1531:                for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
                   1532:                        if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
                   1533:                                h_reducible = 1;
                   1534:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1535:                                _dp_red_mod_destructive(g,p,mod,&u); g = u;
                   1536:                                sugar = MAX(sugar,psugar);
                   1537:                                if ( !g ) {
                   1538:                                        if ( d )
                   1539:                                                d->sugar = sugar;
                   1540:                                        _dptodp(d,rp); _free_dp(d); return;
                   1541:                                }
                   1542:                                break;
                   1543:                        }
                   1544:                }
                   1545:                if ( !h_reducible ) {
                   1546:                        /* head term is not reducible */
                   1547:                        if ( !full ) {
                   1548:                                if ( g )
                   1549:                                        g->sugar = sugar;
                   1550:                                _dptodp(g,rp); _free_dp(g); return;
                   1551:                        } else {
                   1552:                                m = BDY(g);
                   1553:                                if ( NEXT(m) ) {
                   1554:                                        BDY(g) = NEXT(m); NEXT(m) = 0;
                   1555:                                } else {
                   1556:                                        _FREEDP(g); g = 0;
                   1557:                                }
                   1558:                                if ( d ) {
                   1559:                                        for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
                   1560:                                        NEXT(mrd) = m;
                   1561:                                } else {
                   1562:                                        _MKDP(n,m,d);
                   1563:                                }
                   1564:                        }
                   1565:                }
1.5       noro     1566:        }
1.7       noro     1567:        if ( d )
                   1568:                d->sugar = sugar;
                   1569:        _dptodp(d,rp); _free_dp(d);
1.5       noro     1570: }
1.13      noro     1571:
                   1572: /* reduction by linear base over a field */
                   1573:
1.20      noro     1574: void dp_lnf_f(DP p1,DP p2,NODE g,DP *r1p,DP *r2p)
1.13      noro     1575: {
                   1576:        DP r1,r2,b1,b2,t,s;
                   1577:        Obj c,c1,c2;
                   1578:        NODE l,b;
                   1579:        int n;
                   1580:
                   1581:        if ( !p1 ) {
                   1582:                *r1p = p1; *r2p = p2; return;
                   1583:        }
                   1584:        n = p1->nv;
                   1585:        for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
                   1586:                        if ( !r1 ) {
                   1587:                                *r1p = r1; *r2p = r2; return;
                   1588:                        }
                   1589:                        b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
                   1590:                        if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
                   1591:                                b2 = (DP)BDY(NEXT(b));
                   1592:                                divr(CO,(Obj)ONE,(Obj)BDY(b1)->c,&c1);
                   1593:                                mulr(CO,c1,(Obj)BDY(r1)->c,&c2); chsgnr(c2,&c);
                   1594:                                muldc(CO,b1,(P)c,&t); addd(CO,r1,t,&s); r1 = s;
                   1595:                                muldc(CO,b2,(P)c,&t); addd(CO,r2,t,&s); r2 = s;
                   1596:                        }
                   1597:        }
                   1598:        *r1p = r1; *r2p = r2;
                   1599: }
                   1600:
                   1601: /* reduction by linear base over GF(mod) */
1.5       noro     1602:
1.20      noro     1603: void dp_lnf_mod(DP p1,DP p2,NODE g,int mod,DP *r1p,DP *r2p)
1.5       noro     1604: {
1.7       noro     1605:        DP r1,r2,b1,b2,t,s;
                   1606:        P c;
                   1607:        MQ c1,c2;
                   1608:        NODE l,b;
                   1609:        int n;
                   1610:
                   1611:        if ( !p1 ) {
                   1612:                *r1p = p1; *r2p = p2; return;
                   1613:        }
                   1614:        n = p1->nv;
                   1615:        for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
                   1616:                        if ( !r1 ) {
                   1617:                                *r1p = r1; *r2p = r2; return;
                   1618:                        }
                   1619:                        b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
                   1620:                        if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
                   1621:                                b2 = (DP)BDY(NEXT(b));
                   1622:                                invmq(mod,(MQ)BDY(b1)->c,&c1);
                   1623:                                mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
                   1624:                                mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
                   1625:                                mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
                   1626:                        }
                   1627:        }
                   1628:        *r1p = r1; *r2p = r2;
1.5       noro     1629: }
                   1630:
1.20      noro     1631: void dp_nf_tab_mod(DP p,LIST *tab,int mod,DP *rp)
1.5       noro     1632: {
1.7       noro     1633:        DP s,t,u;
                   1634:        MP m;
                   1635:        DL h;
                   1636:        int i,n;
                   1637:
                   1638:        if ( !p ) {
                   1639:                *rp = p; return;
                   1640:        }
                   1641:        n = p->nv;
                   1642:        for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1643:                h = m->dl;
                   1644:                while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
                   1645:                        i++;
                   1646:                mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
                   1647:                addmd(CO,mod,s,t,&u); s = u;
1.24      noro     1648:        }
                   1649:        *rp = s;
                   1650: }
                   1651:
                   1652: void dp_nf_tab_f(DP p,LIST *tab,DP *rp)
                   1653: {
                   1654:        DP s,t,u;
                   1655:        MP m;
                   1656:        DL h;
                   1657:        int i,n;
                   1658:
                   1659:        if ( !p ) {
                   1660:                *rp = p; return;
                   1661:        }
                   1662:        n = p->nv;
                   1663:        for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1664:                h = m->dl;
                   1665:                while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
                   1666:                        i++;
                   1667:                muldc(CO,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
                   1668:                addd(CO,s,t,&u); s = u;
1.7       noro     1669:        }
                   1670:        *rp = s;
1.5       noro     1671: }
                   1672:
1.7       noro     1673: /*
                   1674:  * setting flags
1.30      noro     1675:  * call create_order_spec with vl=0 to set old type order.
1.7       noro     1676:  *
                   1677:  */
                   1678:
1.27      noro     1679: int create_order_spec(VL vl,Obj obj,struct order_spec **specp)
1.5       noro     1680: {
1.37      noro     1681:        int i,j,n,s,row,col,ret;
1.27      noro     1682:        struct order_spec *spec;
1.7       noro     1683:        struct order_pair *l;
                   1684:        NODE node,t,tn;
                   1685:        MAT m;
                   1686:        pointer **b;
                   1687:        int **w;
1.5       noro     1688:
1.37      noro     1689:        if ( vl && obj && OID(obj) == O_LIST ) {
                   1690:                ret = create_composite_order_spec(vl,(LIST)obj,specp);
                   1691:                if ( show_orderspec )
                   1692:                        print_composite_order_spec(*specp);
                   1693:                return ret;
                   1694:        }
1.27      noro     1695:
                   1696:        *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1.7       noro     1697:        if ( !obj || NUM(obj) ) {
                   1698:                spec->id = 0; spec->obj = obj;
                   1699:                spec->ord.simple = QTOS((Q)obj);
                   1700:                return 1;
                   1701:        } else if ( OID(obj) == O_LIST ) {
                   1702:                node = BDY((LIST)obj);
                   1703:                for ( n = 0, t = node; t; t = NEXT(t), n++ );
                   1704:                l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
                   1705:                for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
                   1706:                        tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
                   1707:                        tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
                   1708:                        s += l[i].length;
                   1709:                }
                   1710:                spec->id = 1; spec->obj = obj;
                   1711:                spec->ord.block.order_pair = l;
                   1712:                spec->ord.block.length = n; spec->nv = s;
                   1713:                return 1;
                   1714:        } else if ( OID(obj) == O_MAT ) {
                   1715:                m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
                   1716:                w = almat(row,col);
                   1717:                for ( i = 0; i < row; i++ )
                   1718:                        for ( j = 0; j < col; j++ )
                   1719:                                w[i][j] = QTOS((Q)b[i][j]);
                   1720:                spec->id = 2; spec->obj = obj;
                   1721:                spec->nv = col; spec->ord.matrix.row = row;
                   1722:                spec->ord.matrix.matrix = w;
                   1723:                return 1;
                   1724:        } else
1.5       noro     1725:                return 0;
                   1726: }
                   1727:
1.28      noro     1728: void print_composite_order_spec(struct order_spec *spec)
                   1729: {
                   1730:        int nv,n,len,i,j,k,start;
                   1731:        struct weight_or_block *worb;
                   1732:
                   1733:        nv = spec->nv;
                   1734:        n = spec->ord.composite.length;
                   1735:        worb = spec->ord.composite.w_or_b;
                   1736:        for ( i = 0; i < n; i++, worb++ ) {
                   1737:                len = worb->length;
                   1738:                printf("[ ");
                   1739:                switch ( worb->type ) {
                   1740:                        case IS_DENSE_WEIGHT:
                   1741:                                for ( j = 0; j < len; j++ )
                   1742:                                        printf("%d ",worb->body.dense_weight[j]);
                   1743:                                for ( ; j < nv; j++ )
                   1744:                                        printf("0 ");
                   1745:                                break;
                   1746:                        case IS_SPARSE_WEIGHT:
                   1747:                                for ( j = 0, k = 0; j < nv; j++ )
                   1748:                                        if ( j == worb->body.sparse_weight[k].pos )
                   1749:                                                printf("%d ",worb->body.sparse_weight[k++].value);
                   1750:                                        else
                   1751:                                                printf("0 ");
                   1752:                                break;
                   1753:                        case IS_BLOCK:
                   1754:                                start = worb->body.block.start;
                   1755:                                for ( j = 0; j < start; j++ ) printf("0 ");
                   1756:                                switch ( worb->body.block.order ) {
                   1757:                                        case 0:
                   1758:                                                for ( k = 0; k < len; k++, j++ ) printf("R ");
                   1759:                                                break;
                   1760:                                        case 1:
                   1761:                                                for ( k = 0; k < len; k++, j++ ) printf("G ");
                   1762:                                                break;
                   1763:                                        case 2:
                   1764:                                                for ( k = 0; k < len; k++, j++ ) printf("L ");
                   1765:                                                break;
                   1766:                                }
                   1767:                                for ( ; j < nv; j++ ) printf("0 ");
                   1768:                                break;
                   1769:                }
                   1770:                printf("]\n");
                   1771:        }
1.38      noro     1772: }
                   1773:
                   1774: struct order_spec *append_block(struct order_spec *spec,
                   1775:        int nv,int nalg,int ord)
                   1776: {
                   1777:        MAT m,mat;
                   1778:        int i,j,row,col,n;
                   1779:        Q **b,**wp;
                   1780:        int **w;
                   1781:        NODE t,s,s0;
                   1782:        struct order_pair *l,*l0;
                   1783:        int n0,nv0;
                   1784:        LIST list0,list1,list;
                   1785:        Q oq,nq;
                   1786:        struct order_spec *r;
                   1787:
                   1788:        r = (struct order_spec *)MALLOC(sizeof(struct order_spec));
                   1789:        switch ( spec->id ) {
                   1790:                case 0:
                   1791:                        STOQ(spec->ord.simple,oq); STOQ(nv,nq);
                   1792:                        t = mknode(2,oq,nq); MKLIST(list0,t);
                   1793:                        STOQ(ord,oq); STOQ(nalg,nq);
                   1794:                        t = mknode(2,oq,nq); MKLIST(list1,t);
                   1795:                        t = mknode(2,list0,list1); MKLIST(list,t);
                   1796:                        l = (struct order_pair *)MALLOC_ATOMIC(2*sizeof(struct order_pair));
                   1797:                        l[0].order = spec->ord.simple; l[0].length = nv;
                   1798:                        l[1].order = ord; l[1].length = nalg;
                   1799:                        r->id = 1;  r->obj = (Obj)list;
                   1800:                        r->ord.block.order_pair = l;
                   1801:                        r->ord.block.length = 2;
                   1802:                        r->nv = nv+nalg;
                   1803:                        break;
                   1804:                case 1:
                   1805:                        if ( spec->nv != nv )
                   1806:                                error("append_block : number of variables mismatch");
                   1807:                        l0 = spec->ord.block.order_pair;
                   1808:                        n0 = spec->ord.block.length;
                   1809:                        nv0 = spec->nv;
                   1810:                        list0 = (LIST)spec->obj;
                   1811:                        n = n0+1;
                   1812:                        l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
                   1813:                        for ( i = 0; i < n0; i++ )
                   1814:                                l[i] = l0[i];
                   1815:                        l[i].order = ord; l[i].length = nalg;
                   1816:                         for ( t = BDY(list0), s0 = 0; t; t = NEXT(t) ) {
                   1817:                                NEXTNODE(s0,s); BDY(s) = BDY(t);
                   1818:                        }
                   1819:                        STOQ(ord,oq); STOQ(nalg,nq);
                   1820:                        t = mknode(2,oq,nq); MKLIST(list,t);
                   1821:                        NEXTNODE(s0,s); BDY(s) = (pointer)list; NEXT(s) = 0;
                   1822:                        MKLIST(list,s0);
                   1823:                        r->id = 1;  r->obj = (Obj)list;
                   1824:                        r->ord.block.order_pair = l;
                   1825:                        r->ord.block.length = n;
                   1826:                        r->nv = nv+nalg;
                   1827:                        break;
                   1828:                case 2:
                   1829:                        if ( spec->nv != nv )
                   1830:                                error("append_block : number of variables mismatch");
                   1831:                        m = (MAT)spec->obj;
                   1832:                        row = m->row; col = m->col; b = (Q **)BDY(m);
                   1833:                        w = almat(row+nalg,col+nalg);
                   1834:                        MKMAT(mat,row+nalg,col+nalg); wp = (Q **)BDY(mat);
                   1835:                        for ( i = 0; i < row; i++ )
                   1836:                                for ( j = 0; j < col; j++ ) {
                   1837:                                        w[i][j] = QTOS(b[i][j]);
                   1838:                                        wp[i][j] = b[i][j];
                   1839:                                }
                   1840:                        for ( i = 0; i < nalg; i++ ) {
                   1841:                                w[i+row][i+col] = 1;
                   1842:                                wp[i+row][i+col] = ONE;
                   1843:                        }
                   1844:                        r->id = 2; r->obj = (Obj)mat;
                   1845:                        r->nv = col+nalg; r->ord.matrix.row = row+nalg;
                   1846:                        r->ord.matrix.matrix = w;
                   1847:                        break;
                   1848:                case 3:
                   1849:                default:
                   1850:                        /* XXX */
                   1851:                        error("append_block : not implemented yet");
                   1852:        }
                   1853:        return r;
1.28      noro     1854: }
                   1855:
1.37      noro     1856: int comp_sw(struct sparse_weight *a, struct sparse_weight *b)
                   1857: {
                   1858:        if ( a->pos > b->pos ) return 1;
                   1859:        else if ( a->pos < b->pos ) return -1;
                   1860:        else return 0;
                   1861: }
                   1862:
1.27      noro     1863: /* order = [w_or_b, w_or_b, ... ] */
                   1864: /* w_or_b = w or b                */
                   1865: /* w = [1,2,...] or [x,1,y,2,...] */
                   1866: /* b = [@lex,x,y,...,z] etc       */
                   1867:
                   1868: int create_composite_order_spec(VL vl,LIST order,struct order_spec **specp)
                   1869: {
                   1870:        NODE wb,t,p;
                   1871:        struct order_spec *spec;
                   1872:        VL tvl;
1.29      noro     1873:        int n,i,j,k,l,start,end,len,w;
1.27      noro     1874:        int *dw;
                   1875:        struct sparse_weight *sw;
                   1876:        struct weight_or_block *w_or_b;
                   1877:        Obj a0;
                   1878:        NODE a;
1.29      noro     1879:        V v,sv,ev;
                   1880:        SYMBOL sym;
                   1881:        int *top;
1.27      noro     1882:
                   1883:        /* l = number of vars in vl */
                   1884:        for ( l = 0, tvl = vl; tvl; tvl = NEXT(tvl), l++ );
                   1885:        /* n = number of primitives in order */
                   1886:        wb = BDY(order);
                   1887:        n = length(wb);
                   1888:        *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
                   1889:        spec->id = 3;
                   1890:        spec->obj = (Obj)order;
                   1891:        spec->nv = l;
                   1892:        spec->ord.composite.length = n;
1.28      noro     1893:        w_or_b = spec->ord.composite.w_or_b = (struct weight_or_block *)
1.29      noro     1894:                MALLOC(sizeof(struct weight_or_block)*(n+1));
                   1895:
                   1896:        /* top : register the top variable in each w_or_b specification */
                   1897:        top = (int *)ALLOCA(l*sizeof(int));
                   1898:        for ( i = 0; i < l; i++ ) top[i] = 0;
                   1899:
1.28      noro     1900:        for ( t = wb, i = 0; t; t = NEXT(t), i++ ) {
1.30      noro     1901:                if ( !BDY(t) || OID((Obj)BDY(t)) != O_LIST )
                   1902:                        error("a list of lists must be specified for the key \"order\"");
1.28      noro     1903:                a = BDY((LIST)BDY(t));
1.27      noro     1904:                len = length(a);
                   1905:                a0 = (Obj)BDY(a);
                   1906:                if ( !a0 || OID(a0) == O_N ) {
1.28      noro     1907:                        /* a is a dense weight vector */
1.27      noro     1908:                        dw = (int *)MALLOC(sizeof(int)*len);
1.30      noro     1909:                        for ( j = 0, p = a; j < len; p = NEXT(p), j++ ) {
                   1910:                                if ( !INT((Q)BDY(p)) )
                   1911:                                        error("a dense weight vector must be specified as a list of integers");
1.27      noro     1912:                                dw[j] = QTOS((Q)BDY(p));
1.30      noro     1913:                        }
1.27      noro     1914:                        w_or_b[i].type = IS_DENSE_WEIGHT;
                   1915:                        w_or_b[i].length = len;
                   1916:                        w_or_b[i].body.dense_weight = dw;
1.29      noro     1917:
                   1918:                        /* find the top */
                   1919:                        for ( k = 0; k < len && !dw[k]; k++ );
                   1920:                        if ( k < len ) top[k] = 1;
                   1921:
1.27      noro     1922:                } else if ( OID(a0) == O_P ) {
1.28      noro     1923:                        /* a is a sparse weight vector */
                   1924:                        len >>= 1;
1.27      noro     1925:                        sw = (struct sparse_weight *)
                   1926:                                MALLOC(sizeof(struct sparse_weight)*len);
                   1927:                        for ( j = 0, p = a; j < len; j++ ) {
1.30      noro     1928:                                if ( !BDY(p) || OID((P)BDY(p)) != O_P )
                   1929:                                        error("a sparse weight vector must be specified as [var1,weight1,...]");
1.28      noro     1930:                                v = VR((P)BDY(p)); p = NEXT(p);
1.27      noro     1931:                                for ( tvl = vl, k = 0; tvl && tvl->v != v;
                   1932:                                        k++, tvl = NEXT(tvl) );
                   1933:                                if ( !tvl )
1.30      noro     1934:                                        error("invalid variable name in a sparse weight vector");
1.27      noro     1935:                                sw[j].pos = k;
1.30      noro     1936:                                if ( !INT((Q)BDY(p)) )
                   1937:                                        error("a sparse weight vector must be specified as [var1,weight1,...]");
1.28      noro     1938:                                sw[j].value = QTOS((Q)BDY(p)); p = NEXT(p);
1.27      noro     1939:                        }
1.37      noro     1940:                        qsort(sw,len,sizeof(struct sparse_weight),
                   1941:                                (int (*)(const void *,const void *))comp_sw);
1.27      noro     1942:                        w_or_b[i].type = IS_SPARSE_WEIGHT;
                   1943:                        w_or_b[i].length = len;
                   1944:                        w_or_b[i].body.sparse_weight = sw;
1.29      noro     1945:
                   1946:                        /* find the top */
                   1947:                        for ( k = 0; k < len && !sw[k].value; k++ );
                   1948:                        if ( k < len ) top[sw[k].pos] = 1;
                   1949:                } else if ( OID(a0) == O_RANGE ) {
                   1950:                        /* [range(v1,v2),w] */
                   1951:                        sv = VR((P)(((RANGE)a0)->start));
                   1952:                        ev = VR((P)(((RANGE)a0)->end));
                   1953:                        for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
                   1954:                        if ( !tvl )
                   1955:                                error("invalid range");
                   1956:                        for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
                   1957:                        if ( !tvl )
                   1958:                                error("invalid range");
                   1959:                        len = end-start+1;
                   1960:                        sw = (struct sparse_weight *)
                   1961:                                MALLOC(sizeof(struct sparse_weight)*len);
                   1962:                        w = QTOS((Q)BDY(NEXT(a)));
                   1963:                        for ( tvl = vl, k = 0; k < start; k++, tvl = NEXT(tvl) );
                   1964:                        for ( j = 0 ; k <= end; k++, tvl = NEXT(tvl), j++ ) {
                   1965:                                sw[j].pos = k;
                   1966:                                sw[j].value = w;
                   1967:                        }
                   1968:                        w_or_b[i].type = IS_SPARSE_WEIGHT;
                   1969:                        w_or_b[i].length = len;
                   1970:                        w_or_b[i].body.sparse_weight = sw;
                   1971:
                   1972:                        /* register the top */
                   1973:                        if ( w ) top[start] = 1;
1.28      noro     1974:                } else if ( OID(a0) == O_SYMBOL ) {
                   1975:                        /* a is a block */
1.29      noro     1976:                        sym = (SYMBOL)a0; a = NEXT(a); len--;
                   1977:                        if ( OID((Obj)BDY(a)) == O_RANGE ) {
                   1978:                                sv = VR((P)(((RANGE)BDY(a))->start));
                   1979:                                ev = VR((P)(((RANGE)BDY(a))->end));
                   1980:                                for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
                   1981:                                if ( !tvl )
                   1982:                                        error("invalid range");
                   1983:                                for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
                   1984:                                if ( !tvl )
                   1985:                                        error("invalid range");
                   1986:                                len = end-start+1;
                   1987:                        } else {
                   1988:                                for ( start = 0, tvl = vl; tvl->v != VR((P)BDY(a));
1.28      noro     1989:                                tvl = NEXT(tvl), start++ );
1.29      noro     1990:                                for ( p = NEXT(a), tvl = NEXT(tvl); p;
1.30      noro     1991:                                        p = NEXT(p), tvl = NEXT(tvl) ) {
                   1992:                                        if ( !BDY(p) || OID((P)BDY(p)) != O_P )
                   1993:                                                error("a block must be specified as [ordsymbol,var1,var2,...]");
1.29      noro     1994:                                        if ( tvl->v != VR((P)BDY(p)) ) break;
1.30      noro     1995:                                }
1.29      noro     1996:                                if ( p )
1.30      noro     1997:                                        error("a block must be contiguous in the variable list");
1.29      noro     1998:                        }
1.28      noro     1999:                        w_or_b[i].type = IS_BLOCK;
                   2000:                        w_or_b[i].length = len;
                   2001:                        w_or_b[i].body.block.start = start;
                   2002:                        if ( !strcmp(sym->name,"@grlex") )
                   2003:                                w_or_b[i].body.block.order = 0;
                   2004:                        else if ( !strcmp(sym->name,"@glex") )
                   2005:                                w_or_b[i].body.block.order = 1;
                   2006:                        else if ( !strcmp(sym->name,"@lex") )
                   2007:                                w_or_b[i].body.block.order = 2;
                   2008:                        else
1.29      noro     2009:                                error("invalid ordername");
                   2010:                        /* register the tops */
                   2011:                        for ( j = 0, k = start; j < len; j++, k++ )
                   2012:                                top[k] = 1;
1.28      noro     2013:                }
1.29      noro     2014:        }
                   2015:        for ( k = 0; k < l && top[k]; k++ );
                   2016:        if ( k < l ) {
                   2017:                /* incomplete order specification; add @grlex */
                   2018:                w_or_b[n].type = IS_BLOCK;
                   2019:                w_or_b[n].length = l;
                   2020:                w_or_b[n].body.block.start = 0;
                   2021:                w_or_b[n].body.block.order = 0;
                   2022:                spec->ord.composite.length = n+1;
1.27      noro     2023:        }
                   2024: }
                   2025:
1.35      noro     2026: /* module order spec */
                   2027:
                   2028: void create_modorder_spec(int id,LIST shift,struct modorder_spec **s)
                   2029: {
                   2030:        struct modorder_spec *spec;
                   2031:        NODE n,t;
                   2032:        LIST list;
                   2033:        int *ds;
                   2034:        int i,l;
                   2035:        Q q;
                   2036:
                   2037:        *s = spec = (struct modorder_spec *)MALLOC(sizeof(struct modorder_spec));
                   2038:        spec->id = id;
                   2039:        if ( shift ) {
                   2040:                n = BDY(shift);
                   2041:                spec->len = l = length(n);
                   2042:                spec->degree_shift = ds = (int *)MALLOC_ATOMIC(l*sizeof(int));
                   2043:                for ( t = n, i = 0; t; t = NEXT(t), i++ )
                   2044:                        ds[i] = QTOS((Q)BDY(t));
                   2045:        } else {
                   2046:                spec->len = 0;
                   2047:                spec->degree_shift = 0;
                   2048:        }
                   2049:        STOQ(id,q);
                   2050:        n = mknode(2,q,shift);
                   2051:        MKLIST(list,n);
                   2052:        spec->obj = (Obj)list;
                   2053: }
                   2054:
1.7       noro     2055: /*
                   2056:  * converters
                   2057:  *
                   2058:  */
                   2059:
1.20      noro     2060: void dp_homo(DP p,DP *rp)
1.5       noro     2061: {
1.7       noro     2062:        MP m,mr,mr0;
                   2063:        int i,n,nv,td;
                   2064:        DL dl,dlh;
1.5       noro     2065:
1.7       noro     2066:        if ( !p )
                   2067:                *rp = 0;
                   2068:        else {
                   2069:                n = p->nv; nv = n + 1;
                   2070:                m = BDY(p); td = sugard(m);
                   2071:                for ( mr0 = 0; m; m = NEXT(m) ) {
                   2072:                        NEXTMP(mr0,mr); mr->c = m->c;
                   2073:                        dl = m->dl;
                   2074:                        mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
                   2075:                        dlh->td = td;
                   2076:                        for ( i = 0; i < n; i++ )
                   2077:                                dlh->d[i] = dl->d[i];
                   2078:                        dlh->d[n] = td - dl->td;
                   2079:                }
                   2080:                NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5       noro     2081:        }
                   2082: }
                   2083:
1.20      noro     2084: void dp_dehomo(DP p,DP *rp)
1.5       noro     2085: {
1.7       noro     2086:        MP m,mr,mr0;
                   2087:        int i,n,nv;
                   2088:        DL dl,dlh;
1.5       noro     2089:
1.7       noro     2090:        if ( !p )
                   2091:                *rp = 0;
                   2092:        else {
                   2093:                n = p->nv; nv = n - 1;
                   2094:                m = BDY(p);
                   2095:                for ( mr0 = 0; m; m = NEXT(m) ) {
                   2096:                        NEXTMP(mr0,mr); mr->c = m->c;
                   2097:                        dlh = m->dl;
                   2098:                        mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
                   2099:                        dl->td = dlh->td - dlh->d[nv];
                   2100:                        for ( i = 0; i < nv; i++ )
                   2101:                                dl->d[i] = dlh->d[i];
                   2102:                }
                   2103:                NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
                   2104:        }
1.5       noro     2105: }
                   2106:
1.20      noro     2107: void dp_mod(DP p,int mod,NODE subst,DP *rp)
1.5       noro     2108: {
1.7       noro     2109:        MP m,mr,mr0;
                   2110:        P t,s,s1;
                   2111:        V v;
                   2112:        NODE tn;
1.5       noro     2113:
1.7       noro     2114:        if ( !p )
                   2115:                *rp = 0;
                   2116:        else {
                   2117:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   2118:                        for ( tn = subst, s = m->c; tn; tn = NEXT(tn) ) {
                   2119:                                v = VR((P)BDY(tn)); tn = NEXT(tn);
                   2120:                                substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
                   2121:                        }
                   2122:                        ptomp(mod,s,&t);
                   2123:                        if ( t ) {
                   2124:                                NEXTMP(mr0,mr); mr->c = t; mr->dl = m->dl;
                   2125:                        }
                   2126:                }
                   2127:                if ( mr0 ) {
                   2128:                        NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                   2129:                } else
                   2130:                        *rp = 0;
                   2131:        }
1.5       noro     2132: }
                   2133:
1.20      noro     2134: void dp_rat(DP p,DP *rp)
1.5       noro     2135: {
1.7       noro     2136:        MP m,mr,mr0;
1.5       noro     2137:
1.7       noro     2138:        if ( !p )
                   2139:                *rp = 0;
                   2140:        else {
                   2141:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   2142:                        NEXTMP(mr0,mr); mptop(m->c,&mr->c); mr->dl = m->dl;
1.5       noro     2143:                }
1.7       noro     2144:                if ( mr0 ) {
                   2145:                        NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                   2146:                } else
                   2147:                        *rp = 0;
1.5       noro     2148:        }
                   2149: }
                   2150:
                   2151:
1.27      noro     2152: void homogenize_order(struct order_spec *old,int n,struct order_spec **newp)
1.5       noro     2153: {
1.7       noro     2154:        struct order_pair *l;
                   2155:        int length,nv,row,i,j;
                   2156:        int **newm,**oldm;
1.27      noro     2157:        struct order_spec *new;
1.31      noro     2158:        int onv,nnv,nlen,olen,owlen;
                   2159:        struct weight_or_block *owb,*nwb;
1.5       noro     2160:
1.27      noro     2161:        *newp = new = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1.7       noro     2162:        switch ( old->id ) {
                   2163:                case 0:
                   2164:                        switch ( old->ord.simple ) {
                   2165:                                case 0:
                   2166:                                        new->id = 0; new->ord.simple = 0; break;
                   2167:                                case 1:
                   2168:                                        l = (struct order_pair *)
                   2169:                                                MALLOC_ATOMIC(2*sizeof(struct order_pair));
                   2170:                                        l[0].length = n; l[0].order = 1;
                   2171:                                        l[1].length = 1; l[1].order = 2;
                   2172:                                        new->id = 1;
                   2173:                                        new->ord.block.order_pair = l;
                   2174:                                        new->ord.block.length = 2; new->nv = n+1;
                   2175:                                        break;
                   2176:                                case 2:
                   2177:                                        new->id = 0; new->ord.simple = 1; break;
                   2178:                                case 3: case 4: case 5:
                   2179:                                        new->id = 0; new->ord.simple = old->ord.simple+3;
                   2180:                                        dp_nelim = n-1; break;
                   2181:                                case 6: case 7: case 8: case 9:
                   2182:                                        new->id = 0; new->ord.simple = old->ord.simple; break;
                   2183:                                default:
                   2184:                                        error("homogenize_order : invalid input");
                   2185:                        }
                   2186:                        break;
                   2187:                case 1:
                   2188:                        length = old->ord.block.length;
                   2189:                        l = (struct order_pair *)
                   2190:                                MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
                   2191:                        bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
                   2192:                        l[length].order = 2; l[length].length = 1;
                   2193:                        new->id = 1; new->nv = n+1;
                   2194:                        new->ord.block.order_pair = l;
                   2195:                        new->ord.block.length = length+1;
                   2196:                        break;
                   2197:                case 2:
                   2198:                        nv = old->nv; row = old->ord.matrix.row;
                   2199:                        oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
                   2200:                        for ( i = 0; i <= nv; i++ )
                   2201:                                newm[0][i] = 1;
                   2202:                        for ( i = 0; i < row; i++ ) {
                   2203:                                for ( j = 0; j < nv; j++ )
                   2204:                                        newm[i+1][j] = oldm[i][j];
                   2205:                                newm[i+1][j] = 0;
                   2206:                        }
                   2207:                        new->id = 2; new->nv = nv+1;
                   2208:                        new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
1.31      noro     2209:                        break;
                   2210:                case 3:
                   2211:                        onv = old->nv;
                   2212:                        nnv = onv+1;
                   2213:                        olen = old->ord.composite.length;
                   2214:                        nlen = olen+1;
                   2215:                        owb = old->ord.composite.w_or_b;
                   2216:                        nwb = (struct weight_or_block *)
                   2217:                                MALLOC(nlen*sizeof(struct weight_or_block));
                   2218:                        for ( i = 0; i < olen; i++ ) {
                   2219:                                nwb[i].type = owb[i].type;
                   2220:                                switch ( owb[i].type ) {
                   2221:                                        case IS_DENSE_WEIGHT:
                   2222:                                                owlen = owb[i].length;
                   2223:                                                nwb[i].length = owlen+1;
                   2224:                                                nwb[i].body.dense_weight = (int *)MALLOC((owlen+1)*sizeof(int));
                   2225:                                                for ( j = 0; j < owlen; j++ )
                   2226:                                                        nwb[i].body.dense_weight[j] = owb[i].body.dense_weight[j];
                   2227:                                                nwb[i].body.dense_weight[owlen] = 0;
                   2228:                                                break;
                   2229:                                        case IS_SPARSE_WEIGHT:
                   2230:                                                nwb[i].length = owb[i].length;
                   2231:                                                nwb[i].body.sparse_weight = owb[i].body.sparse_weight;
                   2232:                                                break;
                   2233:                                        case IS_BLOCK:
                   2234:                                                nwb[i].length = owb[i].length;
                   2235:                                                nwb[i].body.block = owb[i].body.block;
                   2236:                                                break;
                   2237:                                }
                   2238:                        }
                   2239:                        nwb[i].type = IS_SPARSE_WEIGHT;
                   2240:                        nwb[i].body.sparse_weight =
                   2241:                                (struct sparse_weight *)MALLOC(sizeof(struct sparse_weight));
                   2242:                        nwb[i].body.sparse_weight[0].pos = onv;
                   2243:                        nwb[i].body.sparse_weight[0].value = 1;
                   2244:                        new->id = 3;
                   2245:                        new->nv = nnv;
                   2246:                        new->ord.composite.length = nlen;
                   2247:                        new->ord.composite.w_or_b = nwb;
                   2248:                        print_composite_order_spec(new);
1.7       noro     2249:                        break;
                   2250:                default:
                   2251:                        error("homogenize_order : invalid input");
1.5       noro     2252:        }
1.7       noro     2253: }
                   2254:
1.20      noro     2255: void qltozl(Q *w,int n,Q *dvr)
1.7       noro     2256: {
                   2257:        N nm,dn;
                   2258:        N g,l1,l2,l3;
                   2259:        Q c,d;
                   2260:        int i;
                   2261:        struct oVECT v;
1.5       noro     2262:
                   2263:        for ( i = 0; i < n; i++ )
1.7       noro     2264:                if ( w[i] && !INT(w[i]) )
                   2265:                        break;
                   2266:        if ( i == n ) {
                   2267:                v.id = O_VECT; v.len = n; v.body = (pointer *)w;
                   2268:                igcdv(&v,dvr); return;
                   2269:        }
                   2270:        c = w[0]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
                   2271:        for ( i = 1; i < n; i++ ) {
                   2272:                c = w[i]; l1 = INT(c) ? ONEN : DN(c);
                   2273:                gcdn(nm,NM(c),&g); nm = g;
                   2274:                gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
1.5       noro     2275:        }
1.7       noro     2276:        if ( UNIN(dn) )
                   2277:                NTOQ(nm,1,d);
                   2278:        else
                   2279:                NDTOQ(nm,dn,1,d);
                   2280:        *dvr = d;
                   2281: }
1.5       noro     2282:
1.20      noro     2283: int comp_nm(Q *a,Q *b)
1.7       noro     2284: {
                   2285:        return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
                   2286: }
                   2287:
1.20      noro     2288: void sortbynm(Q *w,int n)
1.7       noro     2289: {
                   2290:        qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
                   2291: }
1.5       noro     2292:
                   2293:
1.7       noro     2294: /*
                   2295:  * simple operations
                   2296:  *
                   2297:  */
1.5       noro     2298:
1.20      noro     2299: int dp_redble(DP p1,DP p2)
1.7       noro     2300: {
                   2301:        int i,n;
                   2302:        DL d1,d2;
1.5       noro     2303:
1.7       noro     2304:        d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                   2305:        if ( d1->td < d2->td )
                   2306:                return 0;
                   2307:        else {
                   2308:                for ( i = 0, n = p1->nv; i < n; i++ )
                   2309:                        if ( d1->d[i] < d2->d[i] )
                   2310:                                return 0;
                   2311:                return 1;
1.5       noro     2312:        }
                   2313: }
                   2314:
1.20      noro     2315: void dp_subd(DP p1,DP p2,DP *rp)
1.5       noro     2316: {
1.7       noro     2317:        int i,n;
1.5       noro     2318:        DL d1,d2,d;
                   2319:        MP m;
1.7       noro     2320:        DP s;
1.5       noro     2321:
                   2322:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.7       noro     2323:        NEWDL(d,n); d->td = d1->td - d2->td;
1.5       noro     2324:        for ( i = 0; i < n; i++ )
1.7       noro     2325:                d->d[i] = d1->d[i]-d2->d[i];
                   2326:        NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
                   2327:        MKDP(n,m,s); s->sugar = d->td;
                   2328:        *rp = s;
                   2329: }
                   2330:
1.20      noro     2331: void dltod(DL d,int n,DP *rp)
1.7       noro     2332: {
                   2333:        MP m;
                   2334:        DP s;
                   2335:
                   2336:        NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
                   2337:        MKDP(n,m,s); s->sugar = d->td;
                   2338:        *rp = s;
1.5       noro     2339: }
                   2340:
1.20      noro     2341: void dp_hm(DP p,DP *rp)
1.5       noro     2342: {
                   2343:        MP m,mr;
                   2344:
                   2345:        if ( !p )
                   2346:                *rp = 0;
                   2347:        else {
                   2348:                m = BDY(p);
                   2349:                NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
                   2350:                MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td;  /* XXX */
                   2351:        }
                   2352: }
                   2353:
1.35      noro     2354: void dp_ht(DP p,DP *rp)
                   2355: {
                   2356:        MP m,mr;
                   2357:
                   2358:        if ( !p )
                   2359:                *rp = 0;
                   2360:        else {
                   2361:                m = BDY(p);
                   2362:                NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
                   2363:                MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td;  /* XXX */
                   2364:        }
                   2365: }
                   2366:
1.20      noro     2367: void dp_rest(DP p,DP *rp)
1.5       noro     2368: {
                   2369:        MP m;
                   2370:
                   2371:        m = BDY(p);
                   2372:        if ( !NEXT(m) )
                   2373:                *rp = 0;
                   2374:        else {
                   2375:                MKDP(p->nv,NEXT(m),*rp);
                   2376:                if ( *rp )
                   2377:                        (*rp)->sugar = p->sugar;
                   2378:        }
                   2379: }
                   2380:
1.20      noro     2381: DL lcm_of_DL(int nv,DL dl1,DL dl2,DL dl)
1.5       noro     2382: {
1.21      noro     2383:        register int i, *d1, *d2, *d, td;
1.5       noro     2384:
                   2385:        if ( !dl ) NEWDL(dl,nv);
                   2386:        d = dl->d,  d1 = dl1->d,  d2 = dl2->d;
1.21      noro     2387:        for ( td = 0, i = 0; i < nv; d1++, d2++, d++, i++ ) {
                   2388:                *d = *d1 > *d2 ? *d1 : *d2;
                   2389:                td += MUL_WEIGHT(*d,i);
                   2390:        }
1.5       noro     2391:        dl->td = td;
                   2392:        return dl;
                   2393: }
                   2394:
1.20      noro     2395: int dl_equal(int nv,DL dl1,DL dl2)
1.5       noro     2396: {
                   2397:     register int *d1, *d2, n;
                   2398:
                   2399:     if ( dl1->td != dl2->td ) return 0;
                   2400:     for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
                   2401:         if ( *d1 != *d2 ) return 0;
                   2402:     return 1;
                   2403: }
                   2404:
1.20      noro     2405: int dp_nt(DP p)
1.5       noro     2406: {
                   2407:        int i;
                   2408:        MP m;
                   2409:
                   2410:        if ( !p )
                   2411:                return 0;
                   2412:        else {
                   2413:                for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
                   2414:                return i;
                   2415:        }
                   2416: }
                   2417:
1.20      noro     2418: int dp_homogeneous(DP p)
1.15      noro     2419: {
                   2420:        MP m;
                   2421:        int d;
                   2422:
                   2423:        if ( !p )
                   2424:                return 1;
                   2425:        else {
                   2426:                m = BDY(p);
                   2427:                d = m->dl->td;
                   2428:                m = NEXT(m);
                   2429:                for ( ; m; m = NEXT(m) ) {
                   2430:                        if ( m->dl->td != d )
                   2431:                                return 0;
                   2432:                }
                   2433:                return 1;
                   2434:        }
1.16      noro     2435: }
                   2436:
1.20      noro     2437: void _print_mp(int nv,MP m)
1.16      noro     2438: {
                   2439:        int i;
                   2440:
1.17      noro     2441:        if ( !m )
1.16      noro     2442:                return;
                   2443:        for ( ; m; m = NEXT(m) ) {
                   2444:                fprintf(stderr,"%d<",ITOS(C(m)));
                   2445:                for ( i = 0; i < nv; i++ ) {
                   2446:                        fprintf(stderr,"%d",m->dl->d[i]);
                   2447:                        if ( i != nv-1 )
                   2448:                                fprintf(stderr," ");
                   2449:                }
                   2450:                fprintf(stderr,">",C(m));
                   2451:        }
                   2452:        fprintf(stderr,"\n");
1.15      noro     2453: }
1.26      noro     2454:
                   2455: static int cmp_mp_nvar;
                   2456:
                   2457: int comp_mp(MP *a,MP *b)
                   2458: {
                   2459:        return -(*cmpdl)(cmp_mp_nvar,(*a)->dl,(*b)->dl);
                   2460: }
                   2461:
                   2462: void dp_sort(DP p,DP *rp)
                   2463: {
                   2464:        MP t,mp,mp0;
                   2465:        int i,n;
                   2466:        DP r;
                   2467:        MP *w;
                   2468:
                   2469:        if ( !p ) {
                   2470:                *rp = 0;
                   2471:                return;
                   2472:        }
                   2473:        for ( t = BDY(p), n = 0; t; t = NEXT(t), n++ );
                   2474:        w = (MP *)ALLOCA(n*sizeof(MP));
                   2475:        for ( t = BDY(p), i = 0; i < n; t = NEXT(t), i++ )
                   2476:                w[i] = t;
                   2477:        cmp_mp_nvar = NV(p);
                   2478:        qsort(w,n,sizeof(MP),(int (*)(const void *,const void *))comp_mp);
                   2479:        mp0 = 0;
                   2480:        for ( i = n-1; i >= 0; i-- ) {
                   2481:                NEWMP(mp); mp->dl = w[i]->dl; C(mp) = C(w[i]);
                   2482:                NEXT(mp) = mp0; mp0 = mp;
                   2483:        }
                   2484:        MKDP(p->nv,mp0,r);
                   2485:        r->sugar = p->sugar;
                   2486:        *rp = r;
                   2487: }
                   2488:
1.32      noro     2489: DP extract_initial_term_from_dp(DP p,int *weight,int n);
                   2490: LIST extract_initial_term(LIST f,int *weight,int n);
                   2491:
                   2492: DP extract_initial_term_from_dp(DP p,int *weight,int n)
                   2493: {
1.34      noro     2494:        int w,t,i,top;
1.32      noro     2495:        MP m,r0,r;
                   2496:        DP dp;
                   2497:
                   2498:        if ( !p ) return 0;
1.34      noro     2499:        top = 1;
1.32      noro     2500:        for ( m = BDY(p); m; m = NEXT(m) ) {
                   2501:                for ( i = 0, t = 0; i < n; i++ )
                   2502:                        t += weight[i]*m->dl->d[i];
1.34      noro     2503:                if ( top || t > w ) {
1.32      noro     2504:                        r0 = 0;
                   2505:                        w = t;
1.34      noro     2506:                        top = 0;
1.32      noro     2507:                }
                   2508:                if ( t == w ) {
                   2509:                        NEXTMP(r0,r);
                   2510:                        r->dl = m->dl;
                   2511:                        r->c = m->c;
                   2512:                }
                   2513:        }
                   2514:        NEXT(r) = 0;
                   2515:        MKDP(p->nv,r0,dp);
                   2516:        return dp;
                   2517: }
                   2518:
                   2519: LIST extract_initial_term(LIST f,int *weight,int n)
                   2520: {
                   2521:        NODE nd,r0,r;
                   2522:        Obj p;
                   2523:        LIST l;
                   2524:
                   2525:        nd = BDY(f);
                   2526:        for ( r0 = 0; nd; nd = NEXT(nd) ) {
                   2527:                NEXTNODE(r0,r);
                   2528:                p = (Obj)BDY(nd);
                   2529:                BDY(r) = (pointer)extract_initial_term_from_dp((DP)p,weight,n);
                   2530:        }
                   2531:        if ( r0 ) NEXT(r) = 0;
                   2532:        MKLIST(l,r0);
                   2533:        return l;
                   2534: }
                   2535:
                   2536: LIST dp_initial_term(LIST f,struct order_spec *ord)
                   2537: {
                   2538:        int n,l,i;
                   2539:        struct weight_or_block *worb;
                   2540:        int *weight;
                   2541:
                   2542:        switch ( ord->id ) {
                   2543:                case 2: /* matrix order */
                   2544:                        /* extract the first row */
                   2545:                        n = ord->nv;
                   2546:                        weight = ord->ord.matrix.matrix[0];
                   2547:                        return extract_initial_term(f,weight,n);
                   2548:                case 3: /* composite order */
                   2549:                        /* the first w_or_b */
                   2550:                        worb = ord->ord.composite.w_or_b;
                   2551:                        switch ( worb->type ) {
                   2552:                                case IS_DENSE_WEIGHT:
                   2553:                                        n = worb->length;
                   2554:                                        weight = worb->body.dense_weight;
                   2555:                                        return extract_initial_term(f,weight,n);
                   2556:                                case IS_SPARSE_WEIGHT:
                   2557:                                        n = ord->nv;
                   2558:                                        weight = (int *)ALLOCA(n*sizeof(int));
1.33      noro     2559:                                        for ( i = 0; i < n; i++ ) weight[i] = 0;
1.32      noro     2560:                                        l = worb->length;
                   2561:                                        for ( i = 0; i < l; i++ )
                   2562:                                                weight[worb->body.sparse_weight[i].pos]
                   2563:                                                        =  worb->body.sparse_weight[i].value;
                   2564:                                        return extract_initial_term(f,weight,n);
                   2565:                                default:
                   2566:                                        error("dp_initial_term : unsupported order");
                   2567:                        }
                   2568:                default:
                   2569:                        error("dp_initial_term : unsupported order");
                   2570:        }
                   2571: }
                   2572:
                   2573: int highest_order_dp(DP p,int *weight,int n);
                   2574: LIST highest_order(LIST f,int *weight,int n);
                   2575:
                   2576: int highest_order_dp(DP p,int *weight,int n)
                   2577: {
1.34      noro     2578:        int w,t,i,top;
1.32      noro     2579:        MP m;
                   2580:
                   2581:        if ( !p ) return -1;
1.34      noro     2582:        top = 1;
1.32      noro     2583:        for ( m = BDY(p); m; m = NEXT(m) ) {
                   2584:                for ( i = 0, t = 0; i < n; i++ )
                   2585:                        t += weight[i]*m->dl->d[i];
1.34      noro     2586:                if ( top || t > w ) {
1.32      noro     2587:                        w = t;
1.34      noro     2588:                        top = 0;
                   2589:                }
1.32      noro     2590:        }
                   2591:        return w;
                   2592: }
                   2593:
                   2594: LIST highest_order(LIST f,int *weight,int n)
                   2595: {
                   2596:        int h;
                   2597:        NODE nd,r0,r;
                   2598:        Obj p;
                   2599:        LIST l;
                   2600:        Q q;
                   2601:
                   2602:        nd = BDY(f);
                   2603:        for ( r0 = 0; nd; nd = NEXT(nd) ) {
                   2604:                NEXTNODE(r0,r);
                   2605:                p = (Obj)BDY(nd);
                   2606:                h = highest_order_dp((DP)p,weight,n);
                   2607:                STOQ(h,q);
                   2608:                BDY(r) = (pointer)q;
                   2609:        }
                   2610:        if ( r0 ) NEXT(r) = 0;
                   2611:        MKLIST(l,r0);
                   2612:        return l;
                   2613: }
                   2614:
                   2615: LIST dp_order(LIST f,struct order_spec *ord)
                   2616: {
                   2617:        int n,l,i;
                   2618:        struct weight_or_block *worb;
                   2619:        int *weight;
                   2620:
                   2621:        switch ( ord->id ) {
                   2622:                case 2: /* matrix order */
                   2623:                        /* extract the first row */
                   2624:                        n = ord->nv;
                   2625:                        weight = ord->ord.matrix.matrix[0];
                   2626:                        return highest_order(f,weight,n);
                   2627:                case 3: /* composite order */
                   2628:                        /* the first w_or_b */
                   2629:                        worb = ord->ord.composite.w_or_b;
                   2630:                        switch ( worb->type ) {
                   2631:                                case IS_DENSE_WEIGHT:
                   2632:                                        n = worb->length;
                   2633:                                        weight = worb->body.dense_weight;
                   2634:                                        return highest_order(f,weight,n);
                   2635:                                case IS_SPARSE_WEIGHT:
                   2636:                                        n = ord->nv;
                   2637:                                        weight = (int *)ALLOCA(n*sizeof(int));
1.33      noro     2638:                                        for ( i = 0; i < n; i++ ) weight[i] = 0;
1.32      noro     2639:                                        l = worb->length;
                   2640:                                        for ( i = 0; i < l; i++ )
                   2641:                                                weight[worb->body.sparse_weight[i].pos]
                   2642:                                                        =  worb->body.sparse_weight[i].value;
                   2643:                                        return highest_order(f,weight,n);
                   2644:                                default:
                   2645:                                        error("dp_initial_term : unsupported order");
                   2646:                        }
                   2647:                default:
                   2648:                        error("dp_initial_term : unsupported order");
1.35      noro     2649:        }
                   2650: }
                   2651:
                   2652: int dpv_ht(DPV p,DP *h)
                   2653: {
                   2654:        int len,max,maxi,i,t;
                   2655:        DP *e;
                   2656:        MP m,mr;
                   2657:
                   2658:        len = p->len;
                   2659:        e = p->body;
                   2660:        max = -1;
                   2661:        maxi = -1;
                   2662:        for ( i = 0; i < len; i++ )
                   2663:                if ( e[i] && (t = BDY(e[i])->dl->td) > max ) {
                   2664:                        max = t;
                   2665:                        maxi = i;
                   2666:                }
                   2667:        if ( max < 0 ) {
                   2668:                *h = 0;
                   2669:                return -1;
                   2670:        } else {
                   2671:                m = BDY(e[maxi]);
                   2672:                NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
                   2673:                MKDP(e[maxi]->nv,mr,*h); (*h)->sugar = mr->dl->td;  /* XXX */
                   2674:                return maxi;
1.32      noro     2675:        }
                   2676: }
1.42      noro     2677:
                   2678: /* return 1 if 0 <_w1 v && v <_w2 0 */
                   2679:
                   2680: int in_c12(int n,int *v,int row1,int **w1,int row2, int **w2)
                   2681: {
                   2682:        int t1,t2;
                   2683:
                   2684:        t1 = compare_zero(n,v,row1,w1);
                   2685:        t2 = compare_zero(n,v,row2,w2);
                   2686:        if ( t1 > 0 && t2 < 0 ) return 1;
                   2687:        else return 0;
                   2688: }
                   2689:
                   2690: /* 0 < u => 1, 0 > u => -1 */
                   2691:
                   2692: int compare_zero(int n,int *u,int row,int **w)
                   2693: {
                   2694:        int i,j,t;
                   2695:        int *wi;
                   2696:
                   2697:        for ( i = 0; i < row; i++ ) {
                   2698:                wi = w[i];
                   2699:                for ( j = 0, t = 0; j < n; j++ ) t += u[j]*wi[j];
                   2700:                if ( t > 0 ) return 1;
                   2701:                else if ( t < 0 ) return -1;
                   2702:        }
                   2703:        return 0;
                   2704: }
                   2705:
                   2706: /* functions for generic groebner walk */
                   2707: /* u=0 means u=-infty */
                   2708:
                   2709: int compare_facet_preorder(int n,int *u,int *v,
                   2710:        int row1,int **w1,int row2,int **w2)
                   2711: {
                   2712:        int i,j,s,t,tu,tv;
                   2713:        int *w2i,*uv;
                   2714:
                   2715:        if ( !u ) return 1;
                   2716:        uv = W_ALLOC(n);
                   2717:        for ( i = 0; i < row2; i++ ) {
                   2718:                w2i = w2[i];
                   2719:                for ( j = 0, tu = tv = 0; j < n; j++ )
                   2720:                        if ( s = w2i[j] ) {
                   2721:                                tu += s*u[j]; tv += s*v[j];
                   2722:                        }
                   2723:                for ( j = 0; j < n; j++ ) uv[j] = u[j]*tv-v[j]*tu;
                   2724:                t = compare_zero(n,uv,row1,w1);
                   2725:                if ( t > 0 ) return 1;
                   2726:                else if ( t < 0 ) return 0;
                   2727:        }
                   2728:        return 1;
                   2729: }
                   2730:
                   2731: /* return 0 if last_w = infty */
                   2732:
                   2733: NODE compute_last_w(NODE g,NODE gh,int n,int **w,
                   2734:        int row1,int **w1,int row2,int **w2)
                   2735: {
                   2736:        DP d;
                   2737:        MP f,m0,m;
                   2738:        int *wt,*v,*h;
                   2739:        NODE t,s,n0,tn,n1,r0,r;
                   2740:        int i;
                   2741:
                   2742:        wt = W_ALLOC(n);
                   2743:        n0 = 0;
                   2744:        for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) {
                   2745:                f = BDY((DP)BDY(t));
                   2746:                h = BDY((DP)BDY(s))->dl->d;
                   2747:                for ( ; f; f = NEXT(f) ) {
                   2748:                        for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
                   2749:                        for ( i = 0; i < n && !wt[i]; i++ );
                   2750:                        if ( i == n ) continue;
                   2751:
                   2752:                        if ( in_c12(n,wt,row1,w1,row2,w2) &&
                   2753:                                compare_facet_preorder(n,*w,wt,row1,w1,row2,w2) ) {
                   2754:                                v = (int *)MALLOC_ATOMIC(n*sizeof(int));
                   2755:                                for ( i = 0; i < n; i++ ) v[i] = wt[i];
                   2756:                                MKNODE(n1,v,n0); n0 = n1;
                   2757:                        }
                   2758:                }
                   2759:        }
                   2760:        if ( !n0 ) return 0;
                   2761:        for ( t = n0; t; t = NEXT(t) ) {
                   2762:                v = (int *)BDY(t);
                   2763:                for ( s = n0; s; s = NEXT(s) )
                   2764:                        if ( !compare_facet_preorder(n,v,(int *)BDY(s),row1,w1,row2,w2) )
                   2765:                                break;
                   2766:                if ( !s ) {
                   2767:                        *w = v;
                   2768:                        break;
                   2769:                }
                   2770:        }
                   2771:        if ( !t )
                   2772:                error("compute_last_w : cannot happen");
                   2773:        r0 = 0;
                   2774:        for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) {
                   2775:                f = BDY((DP)BDY(t));
                   2776:                h = BDY((DP)BDY(s))->dl->d;
                   2777:                for ( m0 = 0; f; f = NEXT(f) ) {
                   2778:                        for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
                   2779:                        for ( i = 0; i < n && !wt[i]; i++ );
                   2780:                        if ( i == n  ||
                   2781:                                (compare_facet_preorder(n,wt,*w,row1,w1,row2,w2)
                   2782:                                && compare_facet_preorder(n,*w,wt,row1,w1,row2,w2)) ) {
                   2783:                                NEXTMP(m0,m); m->c = f->c; m->dl = f->dl;
                   2784:                        }
                   2785:                }
1.43      noro     2786:                NEXT(m) = 0;
1.42      noro     2787:                MKDP(((DP)BDY(t))->nv,m0,d);  d->sugar = ((DP)BDY(t))->sugar;
                   2788:                NEXTNODE(r0,r); BDY(r) = (pointer)d;
                   2789:        }
                   2790:        NEXT(r) = 0;
                   2791:        return r0;
                   2792: }
1.44      noro     2793:
                   2794: /* compute a sufficient set of d(f)=u-v */
                   2795:
                   2796: NODE compute_essential_df(DP *g,DP *gh,int ng)
                   2797: {
1.45      noro     2798:        int nv,i,j,k,t,lj;
                   2799:        NODE r,r1,ri,rt,r0;
                   2800:        MP m;
                   2801:        MP *mj;
                   2802:        DL di,hj,dl,dlt;
                   2803:        int *d,*dt;
                   2804:        LIST l;
1.44      noro     2805:        Q q;
1.45      noro     2806:
                   2807:        nv = g[0]->nv;
                   2808:        r = 0;
                   2809:        for ( j = 0; j < ng; j++ ) {
                   2810:                for ( m = BDY(g[j]), lj = 0; m; m = NEXT(m), lj++ );
                   2811:                mj = (MP *)ALLOCA(lj*sizeof(MP));
                   2812:                for ( m = BDY(g[j]), k = 0; m; m = NEXT(m), k++ )
                   2813:                        mj[k] = m;
                   2814:                for ( i = 0; i < lj; i++ ) {
                   2815:                        for ( di = mj[i]->dl, k = i+1; k < lj; k++ )
                   2816:                                if ( _dl_redble(di,mj[k]->dl,nv) ) break;
                   2817:                        if ( k < lj ) mj[i] = 0;
                   2818:                }
                   2819:                hj = BDY(gh[j])->dl;
                   2820:                _NEWDL(dl,nv); d = dl->d;
                   2821:                r0 = r;
                   2822:                for ( i = 0; i < lj; i++ ) {
                   2823:                        if ( mj[i] && !dl_equal(nv,di=mj[i]->dl,hj) ) {
                   2824:                                for ( k = 0, t = 0; k < nv; k++ ) {
                   2825:                                        d[k] = hj->d[k]-di->d[k];
                   2826:                                        t += d[k];
                   2827:                                }
                   2828:                                dl->td = t;
                   2829: #if 1
                   2830:                                for ( rt = r0; rt; rt = NEXT(rt) ) {
                   2831:                                        dlt = (DL)BDY(rt);
                   2832:                                        if ( dlt->td != dl->td ) continue;
                   2833:                                        for ( dt = dlt->d, k = 0; k < nv; k++ )
                   2834:                                                if ( d[k] != dt[k] ) break;
                   2835:                                        if ( k == nv ) break;
                   2836:                                }
                   2837: #else
                   2838:                                rt = 0;
                   2839: #endif
                   2840:                                if ( !rt ) {
                   2841:                                        MKNODE(r1,dl,r); r = r1;
                   2842:                                        _NEWDL(dl,nv); d = dl->d;
                   2843:                                }
1.44      noro     2844:                        }
                   2845:                }
                   2846:        }
1.45      noro     2847:        for ( rt = r; rt; rt = NEXT(rt) ) {
                   2848:                dl = (DL)BDY(rt); d = dl->d;
                   2849:                ri = 0;
                   2850:                for ( k = nv-1; k >= 0; k-- ) {
                   2851:                        STOQ(d[k],q);
                   2852:                        MKNODE(r1,q,ri); ri = r1;
1.44      noro     2853:                }
1.45      noro     2854:                MKNODE(r1,0,ri); MKLIST(l,r1);
                   2855:                BDY(rt) = (pointer)l;
1.44      noro     2856:        }
                   2857:        return r;
                   2858: }

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