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

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.64    ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.63 2016/03/31 07:33:32 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:
1.58      noro      654:                node = mknode(4,ONE,NULLP,s1,ONE);
1.14      noro      655:                MKLIST(hist,node);
                    656:                MKNODE(TraceList,hist,0);
                    657:
1.58      noro      658:                node = mknode(4,ONE,NULLP,NULLP,ONE);
1.14      noro      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:
1.58      noro      707:                node = mknode(4,ONE,NULLP,s1,ONE);
1.7       noro      708:                MKLIST(hist,node);
                    709:                MKNODE(TraceList,hist,0);
                    710:
1.58      noro      711:                node = mknode(4,ONE,NULLP,NULLP,ONE);
1.7       noro      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.55      noro      890: void dp_red_marked_mod(DP p0,DP p1,DP p2,DP hp2,int mod,DP *head,DP *rest,P *dnp,DP *multp)
1.44      noro      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:        }
1.55      noro      908:        NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0;
                    909:        MKDP(n,m,s); s->sugar = d->td;
                    910:        *multp = s;
                    911:        mulmd(CO,mod,s,p2,&t);
1.44      noro      912:        if ( NUM(c2) ) {
1.55      noro      913:                submd(CO,mod,p1,t,&r); h = p0;
1.44      noro      914:        } else {
1.55      noro      915:                mulmdc(CO,mod,p1,c2,&s); submd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
1.44      noro      916:        }
                    917:        *head = h; *rest = r; *dnp = c2;
                    918: }
                    919:
1.13      noro      920: /* m-reduction over a field */
                    921:
1.20      noro      922: void dp_red_f(DP p1,DP p2,DP *rest)
1.13      noro      923: {
                    924:        int i,n;
                    925:        DL d1,d2,d;
                    926:        MP m;
1.20      noro      927:        DP t,s;
1.13      noro      928:        Obj a,b;
                    929:
                    930:        n = p1->nv;
                    931:        d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    932:
                    933:        NEWDL(d,n); d->td = d1->td - d2->td;
                    934:        for ( i = 0; i < n; i++ )
                    935:                d->d[i] = d1->d[i]-d2->d[i];
                    936:
                    937:        NEWMP(m); m->dl = d;
                    938:        divr(CO,(Obj)BDY(p1)->c,(Obj)BDY(p2)->c,&a); chsgnr(a,&b);
                    939:        C(m) = (P)b;
                    940:        NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
                    941:
                    942:        muld(CO,s,p2,&t); addd(CO,p1,t,rest);
                    943: }
                    944:
1.20      noro      945: void dp_red_mod(DP p0,DP p1,DP p2,int mod,DP *head,DP *rest,P *dnp)
1.7       noro      946: {
                    947:        int i,n;
                    948:        DL d1,d2,d;
                    949:        MP m;
                    950:        DP t,s,r,h;
                    951:        P c1,c2,g,u;
                    952:
                    953:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    954:        NEWDL(d,n); d->td = d1->td - d2->td;
                    955:        for ( i = 0; i < n; i++ )
                    956:                d->d[i] = d1->d[i]-d2->d[i];
                    957:        c1 = (P)BDY(p1)->c; c2 = (P)BDY(p2)->c;
                    958:        gcdprsmp(CO,mod,c1,c2,&g);
                    959:        divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
                    960:        if ( NUM(c2) ) {
                    961:                divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
                    962:        }
                    963:        NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0;
1.11      noro      964:        MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,s,p2,&t);
1.7       noro      965:        if ( NUM(c2) ) {
                    966:                addmd(CO,mod,p1,t,&r); h = p0;
                    967:        } else {
                    968:                mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
                    969:        }
                    970:        *head = h; *rest = r; *dnp = c2;
                    971: }
                    972:
1.10      noro      973: struct oEGT eg_red_mod;
                    974:
1.20      noro      975: void _dp_red_mod_destructive(DP p1,DP p2,int mod,DP *rp)
1.7       noro      976: {
                    977:        int i,n;
                    978:        DL d1,d2,d;
                    979:        MP m;
                    980:        DP t,s;
1.16      noro      981:        int c,c1,c2;
                    982:        extern int do_weyl;
1.7       noro      983:
                    984:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    985:        _NEWDL(d,n); d->td = d1->td - d2->td;
                    986:        for ( i = 0; i < n; i++ )
                    987:                d->d[i] = d1->d[i]-d2->d[i];
1.16      noro      988:        c = invm(ITOS(BDY(p2)->c),mod);
                    989:        c2 = ITOS(BDY(p1)->c);
                    990:        DMAR(c,c2,0,mod,c1);
1.7       noro      991:        _NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
1.16      noro      992: #if 0
1.7       noro      993:        _MKDP(n,m,s); s->sugar = d->td;
                    994:        _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16      noro      995: #else
                    996:        if ( do_weyl ) {
1.19      noro      997:                _MKDP(n,m,s); s->sugar = d->td;
                    998:                _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16      noro      999:        } else {
                   1000:                _mulmdm_dup(mod,p2,m,&t); _FREEMP(m);
                   1001:        }
                   1002: #endif
1.10      noro     1003: /* get_eg(&t0); */
1.7       noro     1004:        _addmd_destructive(mod,p1,t,rp);
1.10      noro     1005: /* get_eg(&t1); add_eg(&eg_red_mod,&t0,&t1); */
1.7       noro     1006: }
                   1007:
                   1008: /*
                   1009:  * normal form computation
                   1010:  *
                   1011:  */
1.5       noro     1012:
1.20      noro     1013: void dp_true_nf(NODE b,DP g,DP *ps,int full,DP *rp,P *dnp)
1.5       noro     1014: {
                   1015:        DP u,p,d,s,t,dmy;
                   1016:        NODE l;
                   1017:        MP m,mr;
                   1018:        int i,n;
                   1019:        int *wb;
                   1020:        int sugar,psugar;
                   1021:        P dn,tdn,tdn1;
                   1022:
                   1023:        dn = (P)ONE;
                   1024:        if ( !g ) {
                   1025:                *rp = 0; *dnp = dn; return;
                   1026:        }
                   1027:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1028:        wb = (int *)ALLOCA(n*sizeof(int));
                   1029:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1030:                wb[i] = QTOS((Q)BDY(l));
                   1031:        sugar = g->sugar;
                   1032:        for ( d = 0; g; ) {
                   1033:                for ( u = 0, i = 0; i < n; i++ ) {
                   1034:                        if ( dp_redble(g,p = ps[wb[i]]) ) {
                   1035:                                dp_red(d,g,p,&t,&u,&tdn,&dmy);
                   1036:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1037:                                sugar = MAX(sugar,psugar);
                   1038:                                if ( !u ) {
                   1039:                                        if ( d )
                   1040:                                                d->sugar = sugar;
                   1041:                                        *rp = d; *dnp = dn; return;
                   1042:                                } else {
                   1043:                                        d = t;
                   1044:                                        mulp(CO,dn,tdn,&tdn1); dn = tdn1;
                   1045:                                }
                   1046:                                break;
                   1047:                        }
                   1048:                }
                   1049:                if ( u )
                   1050:                        g = u;
                   1051:                else if ( !full ) {
                   1052:                        if ( g ) {
                   1053:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                   1054:                        }
                   1055:                        *rp = g; *dnp = dn; return;
                   1056:                } else {
                   1057:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1058:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1059:                        addd(CO,d,t,&s); d = s;
                   1060:                        dp_rest(g,&t); g = t;
                   1061:                }
                   1062:        }
                   1063:        if ( d )
                   1064:                d->sugar = sugar;
                   1065:        *rp = d; *dnp = dn;
                   1066: }
                   1067:
1.43      noro     1068: void dp_removecont2(DP p1,DP p2,DP *r1p,DP *r2p,Q *contp)
                   1069: {
                   1070:        struct oVECT v;
                   1071:        int i,n1,n2,n;
                   1072:        MP m,m0,t;
                   1073:        Q *w;
                   1074:        Q h;
                   1075:
                   1076:        if ( p1 ) {
                   1077:                for ( i = 0, m = BDY(p1); m; m = NEXT(m), i++ );
                   1078:                n1 = i;
                   1079:        } else
                   1080:                n1 = 0;
                   1081:        if ( p2 ) {
                   1082:                for ( i = 0, m = BDY(p2); m; m = NEXT(m), i++ );
                   1083:                n2 = i;
                   1084:        } else
                   1085:                n2 = 0;
                   1086:        n = n1+n2;
                   1087:        if ( !n ) {
                   1088:                *r1p = 0; *r2p = 0; *contp = ONE; return;
                   1089:        }
                   1090:        w = (Q *)ALLOCA(n*sizeof(Q));
                   1091:        v.len = n;
                   1092:        v.body = (pointer *)w;
                   1093:        i = 0;
                   1094:        if ( p1 )
                   1095:                for ( m = BDY(p1); i < n1; m = NEXT(m), i++ ) w[i] = (Q)m->c;
                   1096:        if ( p2 )
                   1097:                for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) w[i] = (Q)m->c;
                   1098:        h = w[0]; removecont_array((P *)w,n,1); divq(h,w[0],contp);
                   1099:        i = 0;
                   1100:        if ( p1 ) {
                   1101:                for ( m0 = 0, t = BDY(p1); i < n1; i++, t = NEXT(t) ) {
                   1102:                        NEXTMP(m0,m); m->c = (P)w[i]; m->dl = t->dl;
                   1103:                }
                   1104:                NEXT(m) = 0;
                   1105:                MKDP(p1->nv,m0,*r1p); (*r1p)->sugar = p1->sugar;
                   1106:        } else
                   1107:                *r1p = 0;
                   1108:        if ( p2 ) {
                   1109:                for ( m0 = 0, t = BDY(p2); i < n; i++, t = NEXT(t) ) {
                   1110:                        NEXTMP(m0,m); m->c = (P)w[i]; m->dl = t->dl;
                   1111:                }
                   1112:                NEXT(m) = 0;
                   1113:                MKDP(p2->nv,m0,*r2p); (*r2p)->sugar = p2->sugar;
                   1114:        } else
                   1115:                *r2p = 0;
                   1116: }
                   1117:
1.41      noro     1118: /* true nf by a marked GB */
                   1119:
1.43      noro     1120: void dp_true_nf_marked(NODE b,DP g,DP *ps,DP *hps,DP *rp,P *nmp,P *dnp)
1.41      noro     1121: {
                   1122:        DP u,p,d,s,t,dmy,hp;
                   1123:        NODE l;
                   1124:        MP m,mr;
1.43      noro     1125:        int i,n,hmag;
1.41      noro     1126:        int *wb;
1.43      noro     1127:        int sugar,psugar,multiple;
                   1128:        P nm,tnm1,dn,tdn,tdn1;
                   1129:        Q cont;
1.41      noro     1130:
1.43      noro     1131:        multiple = 0;
                   1132:        hmag = multiple*HMAG(g);
                   1133:        nm = (P)ONE;
1.41      noro     1134:        dn = (P)ONE;
                   1135:        if ( !g ) {
                   1136:                *rp = 0; *dnp = dn; return;
                   1137:        }
                   1138:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1139:        wb = (int *)ALLOCA(n*sizeof(int));
                   1140:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1141:                wb[i] = QTOS((Q)BDY(l));
                   1142:        sugar = g->sugar;
                   1143:        for ( d = 0; g; ) {
                   1144:                for ( u = 0, i = 0; i < n; i++ ) {
                   1145:                        if ( dp_redble(g,hp = hps[wb[i]]) ) {
                   1146:                                p = ps[wb[i]];
                   1147:                                dp_red_marked(d,g,p,hp,&t,&u,&tdn,&dmy);
                   1148:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1149:                                sugar = MAX(sugar,psugar);
                   1150:                                if ( !u ) {
1.43      noro     1151:                                        goto last;
1.41      noro     1152:                                } else {
                   1153:                                        d = t;
                   1154:                                        mulp(CO,dn,tdn,&tdn1); dn = tdn1;
                   1155:                                }
                   1156:                                break;
                   1157:                        }
                   1158:                }
1.43      noro     1159:                if ( u ) {
1.41      noro     1160:                        g = u;
1.43      noro     1161:                        if ( multiple && ((d && HMAG(d)>hmag) || (HMAG(g)>hmag)) ) {
                   1162:                                dp_removecont2(d,g,&t,&u,&cont); d = t; g = u;
                   1163:                                mulp(CO,nm,(P)cont,&tnm1); nm = tnm1;
                   1164:                                if ( d )
                   1165:                                        hmag = multiple*HMAG(d);
                   1166:                                else
                   1167:                                        hmag = multiple*HMAG(g);
                   1168:                        }
                   1169:                } else {
1.41      noro     1170:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1171:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1172:                        addd(CO,d,t,&s); d = s;
                   1173:                        dp_rest(g,&t); g = t;
                   1174:                }
                   1175:        }
1.43      noro     1176: last:
                   1177:        if ( d ) {
                   1178:                dp_removecont2(d,0,&t,&u,&cont); d = t;
                   1179:                mulp(CO,nm,(P)cont,&tnm1); nm = tnm1;
1.41      noro     1180:                d->sugar = sugar;
1.43      noro     1181:        }
                   1182:        *rp = d; *nmp = nm; *dnp = dn;
1.41      noro     1183: }
                   1184:
1.44      noro     1185: void dp_true_nf_marked_mod(NODE b,DP g,DP *ps,DP *hps,int mod,DP *rp,P *dnp)
                   1186: {
1.55      noro     1187:        DP hp,u,p,d,s,t,dmy;
1.44      noro     1188:        NODE l;
                   1189:        MP m,mr;
                   1190:        int i,n;
                   1191:        int *wb;
                   1192:        int sugar,psugar;
                   1193:        P dn,tdn,tdn1;
                   1194:
                   1195:        dn = (P)ONEM;
                   1196:        if ( !g ) {
                   1197:                *rp = 0; *dnp = dn; return;
                   1198:        }
                   1199:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1200:                wb = (int *)ALLOCA(n*sizeof(int));
                   1201:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1202:                wb[i] = QTOS((Q)BDY(l));
                   1203:        sugar = g->sugar;
                   1204:        for ( d = 0; g; ) {
                   1205:                for ( u = 0, i = 0; i < n; i++ ) {
                   1206:                        if ( dp_redble(g,hp = hps[wb[i]]) ) {
                   1207:                                p = ps[wb[i]];
1.55      noro     1208:                                dp_red_marked_mod(d,g,p,hp,mod,&t,&u,&tdn,&dmy);
1.44      noro     1209:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1210:                                sugar = MAX(sugar,psugar);
                   1211:                                if ( !u ) {
                   1212:                                        if ( d )
                   1213:                                                d->sugar = sugar;
                   1214:                                        *rp = d; *dnp = dn; return;
                   1215:                                } else {
                   1216:                                        d = t;
                   1217:                                        mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
                   1218:                                }
                   1219:                                break;
                   1220:                        }
                   1221:                }
                   1222:                if ( u )
                   1223:                        g = u;
                   1224:                else {
                   1225:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1226:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1227:                        addmd(CO,mod,d,t,&s); d = s;
                   1228:                        dp_rest(g,&t); g = t;
                   1229:                }
                   1230:        }
                   1231:        if ( d )
                   1232:                d->sugar = sugar;
                   1233:        *rp = d; *dnp = dn;
                   1234: }
                   1235:
1.47      noro     1236: /* true nf by a marked GB and collect quotients */
                   1237:
                   1238: DP *dp_true_nf_and_quotient_marked (NODE b,DP g,DP *ps,DP *hps,DP *rp,P *dnp)
                   1239: {
                   1240:        DP u,p,d,s,t,dmy,hp,mult;
                   1241:        DP *q;
                   1242:        NODE l;
                   1243:        MP m,mr;
                   1244:        int i,n,j;
                   1245:        int *wb;
                   1246:        int sugar,psugar,multiple;
                   1247:        P nm,tnm1,dn,tdn,tdn1;
                   1248:        Q cont;
                   1249:
                   1250:        dn = (P)ONE;
                   1251:        if ( !g ) {
1.59      noro     1252:                *rp = 0; *dnp = dn; return 0;
1.47      noro     1253:        }
                   1254:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1255:        wb = (int *)ALLOCA(n*sizeof(int));
                   1256:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1257:                wb[i] = QTOS((Q)BDY(l));
                   1258:        q = (DP *)MALLOC(n*sizeof(DP));
                   1259:        for ( i = 0; i < n; i++ ) q[i] = 0;
                   1260:        sugar = g->sugar;
                   1261:        for ( d = 0; g; ) {
                   1262:                for ( u = 0, i = 0; i < n; i++ ) {
                   1263:                        if ( dp_redble(g,hp = hps[wb[i]]) ) {
                   1264:                                p = ps[wb[i]];
                   1265:                                dp_red_marked(d,g,p,hp,&t,&u,&tdn,&mult);
                   1266:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1267:                                sugar = MAX(sugar,psugar);
                   1268:                                for ( j = 0; j < n; j++ ) {
                   1269:                                        muldc(CO,q[j],(P)tdn,&dmy); q[j] = dmy;
                   1270:                                }
                   1271:                                addd(CO,q[wb[i]],mult,&dmy); q[wb[i]] = dmy;
                   1272:                                mulp(CO,dn,tdn,&tdn1); dn = tdn1;
                   1273:                                d = t;
                   1274:                                if ( !u ) goto last;
                   1275:                                break;
                   1276:                        }
                   1277:                }
                   1278:                if ( u ) {
                   1279:                        g = u;
                   1280:                } else {
                   1281:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1282:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1283:                        addd(CO,d,t,&s); d = s;
                   1284:                        dp_rest(g,&t); g = t;
                   1285:                }
                   1286:        }
                   1287: last:
                   1288:        if ( d ) d->sugar = sugar;
                   1289:        *rp = d; *dnp = dn;
                   1290:        return q;
                   1291: }
                   1292:
1.55      noro     1293: DP *dp_true_nf_and_quotient_marked_mod(NODE b,DP g,DP *ps,DP *hps,int mod,DP *rp,P *dnp)
                   1294: {
                   1295:        DP u,p,d,s,t,dmy,hp,mult;
                   1296:        DP *q;
                   1297:        NODE l;
                   1298:        MP m,mr;
                   1299:        int i,n,j;
                   1300:        int *wb;
                   1301:        int sugar,psugar;
                   1302:        P dn,tdn,tdn1;
                   1303:
                   1304:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1305:        q = (DP *)MALLOC(n*sizeof(DP));
                   1306:        for ( i = 0; i < n; i++ ) q[i] = 0;
                   1307:        dn = (P)ONEM;
                   1308:        if ( !g ) {
1.59      noro     1309:                *rp = 0; *dnp = dn; return 0;
1.55      noro     1310:        }
                   1311:        wb = (int *)ALLOCA(n*sizeof(int));
                   1312:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1313:                wb[i] = QTOS((Q)BDY(l));
                   1314:        sugar = g->sugar;
                   1315:        for ( d = 0; g; ) {
                   1316:                for ( u = 0, i = 0; i < n; i++ ) {
                   1317:                        if ( dp_redble(g,hp = hps[wb[i]]) ) {
                   1318:                                p = ps[wb[i]];
                   1319:                                dp_red_marked_mod(d,g,p,hp,mod,&t,&u,&tdn,&mult);
                   1320:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1321:                                sugar = MAX(sugar,psugar);
                   1322:                                for ( j = 0; j < n; j++ ) {
                   1323:                                        mulmdc(CO,mod,q[j],(P)tdn,&dmy); q[j] = dmy;
                   1324:                                }
                   1325:                                addmd(CO,mod,q[wb[i]],mult,&dmy); q[wb[i]] = dmy;
                   1326:                                mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
                   1327:                                d = t;
                   1328:                                if ( !u ) goto last;
                   1329:                                break;
                   1330:                        }
                   1331:                }
                   1332:                if ( u )
                   1333:                        g = u;
                   1334:                else {
                   1335:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1336:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1337:                        addmd(CO,mod,d,t,&s); d = s;
                   1338:                        dp_rest(g,&t); g = t;
                   1339:                }
                   1340:        }
                   1341: last:
                   1342:        if ( d )
                   1343:                d->sugar = sugar;
                   1344:        *rp = d; *dnp = dn;
                   1345:        return q;
                   1346: }
                   1347:
1.13      noro     1348: /* nf computation over Z */
                   1349:
1.20      noro     1350: void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *rp)
1.5       noro     1351: {
                   1352:        DP u,p,d,s,t,dmy1;
                   1353:        P dmy;
                   1354:        NODE l;
                   1355:        MP m,mr;
                   1356:        int i,n;
                   1357:        int *wb;
                   1358:        int hmag;
                   1359:        int sugar,psugar;
                   1360:
                   1361:        if ( !g ) {
                   1362:                *rp = 0; return;
                   1363:        }
                   1364:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1365:        wb = (int *)ALLOCA(n*sizeof(int));
                   1366:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1367:                wb[i] = QTOS((Q)BDY(l));
1.12      noro     1368:
1.13      noro     1369:        hmag = multiple*HMAG(g);
1.5       noro     1370:        sugar = g->sugar;
1.12      noro     1371:
1.5       noro     1372:        for ( d = 0; g; ) {
                   1373:                for ( u = 0, i = 0; i < n; i++ ) {
                   1374:                        if ( dp_redble(g,p = ps[wb[i]]) ) {
                   1375:                                dp_red(d,g,p,&t,&u,&dmy,&dmy1);
                   1376:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1377:                                sugar = MAX(sugar,psugar);
                   1378:                                if ( !u ) {
                   1379:                                        if ( d )
                   1380:                                                d->sugar = sugar;
                   1381:                                        *rp = d; return;
                   1382:                                }
                   1383:                                d = t;
                   1384:                                break;
                   1385:                        }
                   1386:                }
                   1387:                if ( u ) {
                   1388:                        g = u;
                   1389:                        if ( d ) {
1.13      noro     1390:                                if ( multiple && HMAG(d) > hmag ) {
1.5       noro     1391:                                        dp_ptozp2(d,g,&t,&u); d = t; g = u;
                   1392:                                        hmag = multiple*HMAG(d);
                   1393:                                }
                   1394:                        } else {
1.13      noro     1395:                                if ( multiple && HMAG(g) > hmag ) {
1.5       noro     1396:                                        dp_ptozp(g,&t); g = t;
                   1397:                                        hmag = multiple*HMAG(g);
                   1398:                                }
                   1399:                        }
                   1400:                }
                   1401:                else if ( !full ) {
                   1402:                        if ( g ) {
                   1403:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                   1404:                        }
                   1405:                        *rp = g; return;
                   1406:                } else {
                   1407:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1408:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1409:                        addd(CO,d,t,&s); d = s;
                   1410:                        dp_rest(g,&t); g = t;
                   1411:
                   1412:                }
                   1413:        }
                   1414:        if ( d )
                   1415:                d->sugar = sugar;
                   1416:        *rp = d;
                   1417: }
                   1418:
1.13      noro     1419: /* nf computation over a field */
                   1420:
1.20      noro     1421: void dp_nf_f(NODE b,DP g,DP *ps,int full,DP *rp)
1.13      noro     1422: {
                   1423:        DP u,p,d,s,t;
                   1424:        NODE l;
                   1425:        MP m,mr;
                   1426:        int i,n;
                   1427:        int *wb;
                   1428:        int sugar,psugar;
                   1429:
                   1430:        if ( !g ) {
                   1431:                *rp = 0; return;
                   1432:        }
                   1433:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1434:        wb = (int *)ALLOCA(n*sizeof(int));
                   1435:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1436:                wb[i] = QTOS((Q)BDY(l));
                   1437:
                   1438:        sugar = g->sugar;
                   1439:        for ( d = 0; g; ) {
                   1440:                for ( u = 0, i = 0; i < n; i++ ) {
                   1441:                        if ( dp_redble(g,p = ps[wb[i]]) ) {
                   1442:                                dp_red_f(g,p,&u);
                   1443:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1444:                                sugar = MAX(sugar,psugar);
                   1445:                                if ( !u ) {
                   1446:                                        if ( d )
                   1447:                                                d->sugar = sugar;
                   1448:                                        *rp = d; return;
                   1449:                                }
                   1450:                                break;
                   1451:                        }
                   1452:                }
                   1453:                if ( u )
                   1454:                        g = u;
                   1455:                else if ( !full ) {
                   1456:                        if ( g ) {
                   1457:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                   1458:                        }
                   1459:                        *rp = g; return;
                   1460:                } else {
                   1461:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1462:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1463:                        addd(CO,d,t,&s); d = s;
                   1464:                        dp_rest(g,&t); g = t;
                   1465:                }
                   1466:        }
                   1467:        if ( d )
                   1468:                d->sugar = sugar;
                   1469:        *rp = d;
                   1470: }
                   1471:
                   1472: /* nf computation over GF(mod) (only for internal use) */
                   1473:
1.20      noro     1474: void dp_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5       noro     1475: {
                   1476:        DP u,p,d,s,t;
                   1477:        P dmy;
                   1478:        NODE l;
                   1479:        MP m,mr;
                   1480:        int sugar,psugar;
                   1481:
                   1482:        if ( !g ) {
                   1483:                *rp = 0; return;
                   1484:        }
                   1485:        sugar = g->sugar;
                   1486:        for ( d = 0; g; ) {
                   1487:                for ( u = 0, l = b; l; l = NEXT(l) ) {
                   1488:                        if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
                   1489:                                dp_red_mod(d,g,p,mod,&t,&u,&dmy);
                   1490:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1491:                                sugar = MAX(sugar,psugar);
                   1492:                                if ( !u ) {
                   1493:                                        if ( d )
                   1494:                                                d->sugar = sugar;
                   1495:                                        *rp = d; return;
                   1496:                                }
                   1497:                                d = t;
                   1498:                                break;
                   1499:                        }
                   1500:                }
                   1501:                if ( u )
                   1502:                        g = u;
                   1503:                else if ( !full ) {
                   1504:                        if ( g ) {
                   1505:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                   1506:                        }
                   1507:                        *rp = g; return;
                   1508:                } else {
                   1509:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1510:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1511:                        addmd(CO,mod,d,t,&s); d = s;
                   1512:                        dp_rest(g,&t); g = t;
                   1513:                }
                   1514:        }
                   1515:        if ( d )
                   1516:                d->sugar = sugar;
                   1517:        *rp = d;
                   1518: }
                   1519:
1.20      noro     1520: void dp_true_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp,P *dnp)
1.5       noro     1521: {
                   1522:        DP u,p,d,s,t;
                   1523:        NODE l;
                   1524:        MP m,mr;
                   1525:        int i,n;
                   1526:        int *wb;
                   1527:        int sugar,psugar;
                   1528:        P dn,tdn,tdn1;
                   1529:
                   1530:        dn = (P)ONEM;
                   1531:        if ( !g ) {
                   1532:                *rp = 0; *dnp = dn; return;
                   1533:        }
                   1534:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                   1535:                wb = (int *)ALLOCA(n*sizeof(int));
                   1536:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                   1537:                wb[i] = QTOS((Q)BDY(l));
                   1538:        sugar = g->sugar;
                   1539:        for ( d = 0; g; ) {
                   1540:                for ( u = 0, i = 0; i < n; i++ ) {
                   1541:                        if ( dp_redble(g,p = ps[wb[i]]) ) {
                   1542:                                dp_red_mod(d,g,p,mod,&t,&u,&tdn);
                   1543:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1544:                                sugar = MAX(sugar,psugar);
                   1545:                                if ( !u ) {
                   1546:                                        if ( d )
                   1547:                                                d->sugar = sugar;
                   1548:                                        *rp = d; *dnp = dn; return;
                   1549:                                } else {
                   1550:                                        d = t;
                   1551:                                        mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
                   1552:                                }
                   1553:                                break;
                   1554:                        }
                   1555:                }
                   1556:                if ( u )
                   1557:                        g = u;
                   1558:                else if ( !full ) {
                   1559:                        if ( g ) {
                   1560:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                   1561:                        }
                   1562:                        *rp = g; *dnp = dn; return;
                   1563:                } else {
                   1564:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1565:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1566:                        addmd(CO,mod,d,t,&s); d = s;
                   1567:                        dp_rest(g,&t); g = t;
                   1568:                }
                   1569:        }
                   1570:        if ( d )
                   1571:                d->sugar = sugar;
                   1572:        *rp = d; *dnp = dn;
                   1573: }
                   1574:
1.20      noro     1575: void _dp_nf_mod_destructive(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5       noro     1576: {
1.20      noro     1577:        DP u,p,d;
1.7       noro     1578:        NODE l;
1.20      noro     1579:        MP m,mrd;
                   1580:        int sugar,psugar,n,h_reducible;
1.5       noro     1581:
1.7       noro     1582:        if ( !g ) {
                   1583:                *rp = 0; return;
1.5       noro     1584:        }
1.7       noro     1585:        sugar = g->sugar;
                   1586:        n = g->nv;
                   1587:        for ( d = 0; g; ) {
                   1588:                for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
                   1589:                        if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
                   1590:                                h_reducible = 1;
                   1591:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1592:                                _dp_red_mod_destructive(g,p,mod,&u); g = u;
                   1593:                                sugar = MAX(sugar,psugar);
                   1594:                                if ( !g ) {
                   1595:                                        if ( d )
                   1596:                                                d->sugar = sugar;
                   1597:                                        _dptodp(d,rp); _free_dp(d); return;
                   1598:                                }
                   1599:                                break;
                   1600:                        }
                   1601:                }
                   1602:                if ( !h_reducible ) {
                   1603:                        /* head term is not reducible */
                   1604:                        if ( !full ) {
                   1605:                                if ( g )
                   1606:                                        g->sugar = sugar;
                   1607:                                _dptodp(g,rp); _free_dp(g); return;
                   1608:                        } else {
                   1609:                                m = BDY(g);
                   1610:                                if ( NEXT(m) ) {
                   1611:                                        BDY(g) = NEXT(m); NEXT(m) = 0;
                   1612:                                } else {
                   1613:                                        _FREEDP(g); g = 0;
                   1614:                                }
                   1615:                                if ( d ) {
                   1616:                                        for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
                   1617:                                        NEXT(mrd) = m;
                   1618:                                } else {
                   1619:                                        _MKDP(n,m,d);
                   1620:                                }
                   1621:                        }
                   1622:                }
1.5       noro     1623:        }
1.7       noro     1624:        if ( d )
                   1625:                d->sugar = sugar;
                   1626:        _dptodp(d,rp); _free_dp(d);
1.5       noro     1627: }
1.13      noro     1628:
                   1629: /* reduction by linear base over a field */
                   1630:
1.20      noro     1631: void dp_lnf_f(DP p1,DP p2,NODE g,DP *r1p,DP *r2p)
1.13      noro     1632: {
                   1633:        DP r1,r2,b1,b2,t,s;
                   1634:        Obj c,c1,c2;
                   1635:        NODE l,b;
                   1636:        int n;
                   1637:
                   1638:        if ( !p1 ) {
                   1639:                *r1p = p1; *r2p = p2; return;
                   1640:        }
                   1641:        n = p1->nv;
                   1642:        for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
                   1643:                        if ( !r1 ) {
                   1644:                                *r1p = r1; *r2p = r2; return;
                   1645:                        }
                   1646:                        b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
                   1647:                        if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
                   1648:                                b2 = (DP)BDY(NEXT(b));
                   1649:                                divr(CO,(Obj)ONE,(Obj)BDY(b1)->c,&c1);
                   1650:                                mulr(CO,c1,(Obj)BDY(r1)->c,&c2); chsgnr(c2,&c);
                   1651:                                muldc(CO,b1,(P)c,&t); addd(CO,r1,t,&s); r1 = s;
                   1652:                                muldc(CO,b2,(P)c,&t); addd(CO,r2,t,&s); r2 = s;
                   1653:                        }
                   1654:        }
                   1655:        *r1p = r1; *r2p = r2;
                   1656: }
                   1657:
                   1658: /* reduction by linear base over GF(mod) */
1.5       noro     1659:
1.20      noro     1660: void dp_lnf_mod(DP p1,DP p2,NODE g,int mod,DP *r1p,DP *r2p)
1.5       noro     1661: {
1.7       noro     1662:        DP r1,r2,b1,b2,t,s;
                   1663:        P c;
                   1664:        MQ c1,c2;
                   1665:        NODE l,b;
                   1666:        int n;
                   1667:
                   1668:        if ( !p1 ) {
                   1669:                *r1p = p1; *r2p = p2; return;
                   1670:        }
                   1671:        n = p1->nv;
                   1672:        for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
                   1673:                        if ( !r1 ) {
                   1674:                                *r1p = r1; *r2p = r2; return;
                   1675:                        }
                   1676:                        b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
                   1677:                        if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
                   1678:                                b2 = (DP)BDY(NEXT(b));
                   1679:                                invmq(mod,(MQ)BDY(b1)->c,&c1);
                   1680:                                mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
                   1681:                                mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
                   1682:                                mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
                   1683:                        }
                   1684:        }
                   1685:        *r1p = r1; *r2p = r2;
1.5       noro     1686: }
                   1687:
1.20      noro     1688: void dp_nf_tab_mod(DP p,LIST *tab,int mod,DP *rp)
1.5       noro     1689: {
1.7       noro     1690:        DP s,t,u;
                   1691:        MP m;
                   1692:        DL h;
                   1693:        int i,n;
                   1694:
                   1695:        if ( !p ) {
                   1696:                *rp = p; return;
                   1697:        }
                   1698:        n = p->nv;
                   1699:        for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1700:                h = m->dl;
                   1701:                while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
                   1702:                        i++;
                   1703:                mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
                   1704:                addmd(CO,mod,s,t,&u); s = u;
1.24      noro     1705:        }
                   1706:        *rp = s;
                   1707: }
                   1708:
                   1709: void dp_nf_tab_f(DP p,LIST *tab,DP *rp)
                   1710: {
                   1711:        DP s,t,u;
                   1712:        MP m;
                   1713:        DL h;
                   1714:        int i,n;
                   1715:
                   1716:        if ( !p ) {
                   1717:                *rp = p; return;
                   1718:        }
                   1719:        n = p->nv;
                   1720:        for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1721:                h = m->dl;
                   1722:                while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
                   1723:                        i++;
                   1724:                muldc(CO,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
                   1725:                addd(CO,s,t,&u); s = u;
1.7       noro     1726:        }
                   1727:        *rp = s;
1.5       noro     1728: }
                   1729:
1.7       noro     1730: /*
                   1731:  * setting flags
1.30      noro     1732:  * call create_order_spec with vl=0 to set old type order.
1.7       noro     1733:  *
                   1734:  */
                   1735:
1.27      noro     1736: int create_order_spec(VL vl,Obj obj,struct order_spec **specp)
1.5       noro     1737: {
1.62      noro     1738:        int i,j,n,s,row,col,ret,wlen;
1.27      noro     1739:        struct order_spec *spec;
1.7       noro     1740:        struct order_pair *l;
1.62      noro     1741:   Obj wp,wm;
                   1742:        NODE node,t,tn,wpair;
1.7       noro     1743:        MAT m;
1.49      noro     1744:        VECT v;
                   1745:        pointer **b,*bv;
1.7       noro     1746:        int **w;
1.5       noro     1747:
1.37      noro     1748:        if ( vl && obj && OID(obj) == O_LIST ) {
                   1749:                ret = create_composite_order_spec(vl,(LIST)obj,specp);
                   1750:                if ( show_orderspec )
                   1751:                        print_composite_order_spec(*specp);
                   1752:                return ret;
                   1753:        }
1.27      noro     1754:
                   1755:        *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1.7       noro     1756:        if ( !obj || NUM(obj) ) {
                   1757:                spec->id = 0; spec->obj = obj;
                   1758:                spec->ord.simple = QTOS((Q)obj);
                   1759:                return 1;
                   1760:        } else if ( OID(obj) == O_LIST ) {
1.62      noro     1761:     /* module order; obj = [0|1,w,ord] or [0|1,ord] */
1.7       noro     1762:                node = BDY((LIST)obj);
1.53      noro     1763:                if ( !BDY(node) || NUM(BDY(node)) ) {
1.62      noro     1764:       switch ( length(node) ) {
                   1765:       case 2:
                   1766:                          create_order_spec(0,(Obj)BDY(NEXT(node)),&spec);
                   1767:                          spec->id += 256; spec->obj = obj;
                   1768:         spec->top_weight = 0;
                   1769:         spec->module_rank = 0;
                   1770:         spec->module_top_weight = 0;
                   1771:                          spec->ispot = (BDY(node)!=0);
                   1772:                          if ( spec->ispot ) {
                   1773:                                n = QTOS((Q)BDY(node));
                   1774:                                if ( n < 0 )
                   1775:                                        spec->pot_nelim = -n;
                   1776:                                else
                   1777:                                        spec->pot_nelim = 0;
                   1778:                          }
                   1779:         break;
                   1780:
                   1781:       case 3:
                   1782:                          create_order_spec(0,(Obj)BDY(NEXT(NEXT(node))),&spec);
                   1783:                          spec->id += 256; spec->obj = obj;
                   1784:                          spec->ispot = (BDY(node)!=0);
                   1785:         node = NEXT(node);
                   1786:         if ( !BDY(node) || OID(BDY(node)) != O_LIST )
                   1787:                                  error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight");
                   1788:         wpair = BDY((LIST)BDY(node));
                   1789:         if ( length(wpair) != 2 )
                   1790:                                  error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight");
                   1791:
                   1792:         wp = BDY(wpair);
                   1793:         wm = BDY(NEXT(wpair));
                   1794:         if ( !wp || OID(wp) != O_LIST || !wm || OID(wm) != O_LIST )
                   1795:                                  error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight");
                   1796:         spec->nv = length(BDY((LIST)wp));
                   1797:         spec->top_weight = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
                   1798:                    for ( i = 0, t = BDY((LIST)wp); i < spec->nv; t = NEXT(t), i++ )
                   1799:           spec->top_weight[i] = QTOS((Q)BDY(t));
                   1800:
                   1801:         spec->module_rank = length(BDY((LIST)wm));
                   1802:         spec->module_top_weight = (int *)MALLOC_ATOMIC(spec->module_rank*sizeof(int));
                   1803:                    for ( i = 0, t = BDY((LIST)wm); i < spec->module_rank; t = NEXT(t), i++ )
                   1804:           spec->module_top_weight[i] = QTOS((Q)BDY(t));
                   1805:         break;
                   1806:       default:
                   1807:                                error("create_order_spec : invalid arguments for module order");
                   1808:       }
                   1809:
1.53      noro     1810:                        *specp = spec;
                   1811:                        return 1;
1.62      noro     1812:                } else {
                   1813:       /* block order in polynomial ring */
                   1814:                  for ( n = 0, t = node; t; t = NEXT(t), n++ );
                   1815:                  l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
                   1816:                  for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
                   1817:                          tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
                   1818:                          tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
                   1819:                          s += l[i].length;
                   1820:                  }
                   1821:                  spec->id = 1; spec->obj = obj;
                   1822:                  spec->ord.block.order_pair = l;
                   1823:                  spec->ord.block.length = n; spec->nv = s;
                   1824:                  return 1;
                   1825:     }
1.7       noro     1826:        } else if ( OID(obj) == O_MAT ) {
                   1827:                m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
                   1828:                w = almat(row,col);
                   1829:                for ( i = 0; i < row; i++ )
                   1830:                        for ( j = 0; j < col; j++ )
                   1831:                                w[i][j] = QTOS((Q)b[i][j]);
                   1832:                spec->id = 2; spec->obj = obj;
                   1833:                spec->nv = col; spec->ord.matrix.row = row;
                   1834:                spec->ord.matrix.matrix = w;
                   1835:                return 1;
                   1836:        } else
1.5       noro     1837:                return 0;
                   1838: }
                   1839:
1.28      noro     1840: void print_composite_order_spec(struct order_spec *spec)
                   1841: {
                   1842:        int nv,n,len,i,j,k,start;
                   1843:        struct weight_or_block *worb;
                   1844:
                   1845:        nv = spec->nv;
                   1846:        n = spec->ord.composite.length;
                   1847:        worb = spec->ord.composite.w_or_b;
                   1848:        for ( i = 0; i < n; i++, worb++ ) {
                   1849:                len = worb->length;
                   1850:                printf("[ ");
                   1851:                switch ( worb->type ) {
                   1852:                        case IS_DENSE_WEIGHT:
                   1853:                                for ( j = 0; j < len; j++ )
                   1854:                                        printf("%d ",worb->body.dense_weight[j]);
                   1855:                                for ( ; j < nv; j++ )
                   1856:                                        printf("0 ");
                   1857:                                break;
                   1858:                        case IS_SPARSE_WEIGHT:
                   1859:                                for ( j = 0, k = 0; j < nv; j++ )
                   1860:                                        if ( j == worb->body.sparse_weight[k].pos )
                   1861:                                                printf("%d ",worb->body.sparse_weight[k++].value);
                   1862:                                        else
                   1863:                                                printf("0 ");
                   1864:                                break;
                   1865:                        case IS_BLOCK:
                   1866:                                start = worb->body.block.start;
                   1867:                                for ( j = 0; j < start; j++ ) printf("0 ");
                   1868:                                switch ( worb->body.block.order ) {
                   1869:                                        case 0:
                   1870:                                                for ( k = 0; k < len; k++, j++ ) printf("R ");
                   1871:                                                break;
                   1872:                                        case 1:
                   1873:                                                for ( k = 0; k < len; k++, j++ ) printf("G ");
                   1874:                                                break;
                   1875:                                        case 2:
                   1876:                                                for ( k = 0; k < len; k++, j++ ) printf("L ");
                   1877:                                                break;
                   1878:                                }
                   1879:                                for ( ; j < nv; j++ ) printf("0 ");
                   1880:                                break;
                   1881:                }
                   1882:                printf("]\n");
                   1883:        }
1.38      noro     1884: }
                   1885:
                   1886: struct order_spec *append_block(struct order_spec *spec,
                   1887:        int nv,int nalg,int ord)
                   1888: {
                   1889:        MAT m,mat;
                   1890:        int i,j,row,col,n;
                   1891:        Q **b,**wp;
                   1892:        int **w;
                   1893:        NODE t,s,s0;
                   1894:        struct order_pair *l,*l0;
                   1895:        int n0,nv0;
                   1896:        LIST list0,list1,list;
                   1897:        Q oq,nq;
                   1898:        struct order_spec *r;
                   1899:
                   1900:        r = (struct order_spec *)MALLOC(sizeof(struct order_spec));
                   1901:        switch ( spec->id ) {
                   1902:                case 0:
                   1903:                        STOQ(spec->ord.simple,oq); STOQ(nv,nq);
                   1904:                        t = mknode(2,oq,nq); MKLIST(list0,t);
                   1905:                        STOQ(ord,oq); STOQ(nalg,nq);
                   1906:                        t = mknode(2,oq,nq); MKLIST(list1,t);
                   1907:                        t = mknode(2,list0,list1); MKLIST(list,t);
                   1908:                        l = (struct order_pair *)MALLOC_ATOMIC(2*sizeof(struct order_pair));
                   1909:                        l[0].order = spec->ord.simple; l[0].length = nv;
                   1910:                        l[1].order = ord; l[1].length = nalg;
                   1911:                        r->id = 1;  r->obj = (Obj)list;
                   1912:                        r->ord.block.order_pair = l;
                   1913:                        r->ord.block.length = 2;
                   1914:                        r->nv = nv+nalg;
                   1915:                        break;
                   1916:                case 1:
                   1917:                        if ( spec->nv != nv )
                   1918:                                error("append_block : number of variables mismatch");
                   1919:                        l0 = spec->ord.block.order_pair;
                   1920:                        n0 = spec->ord.block.length;
                   1921:                        nv0 = spec->nv;
                   1922:                        list0 = (LIST)spec->obj;
                   1923:                        n = n0+1;
                   1924:                        l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
                   1925:                        for ( i = 0; i < n0; i++ )
                   1926:                                l[i] = l0[i];
                   1927:                        l[i].order = ord; l[i].length = nalg;
                   1928:                         for ( t = BDY(list0), s0 = 0; t; t = NEXT(t) ) {
                   1929:                                NEXTNODE(s0,s); BDY(s) = BDY(t);
                   1930:                        }
                   1931:                        STOQ(ord,oq); STOQ(nalg,nq);
                   1932:                        t = mknode(2,oq,nq); MKLIST(list,t);
                   1933:                        NEXTNODE(s0,s); BDY(s) = (pointer)list; NEXT(s) = 0;
                   1934:                        MKLIST(list,s0);
                   1935:                        r->id = 1;  r->obj = (Obj)list;
                   1936:                        r->ord.block.order_pair = l;
                   1937:                        r->ord.block.length = n;
                   1938:                        r->nv = nv+nalg;
                   1939:                        break;
                   1940:                case 2:
                   1941:                        if ( spec->nv != nv )
                   1942:                                error("append_block : number of variables mismatch");
                   1943:                        m = (MAT)spec->obj;
                   1944:                        row = m->row; col = m->col; b = (Q **)BDY(m);
                   1945:                        w = almat(row+nalg,col+nalg);
                   1946:                        MKMAT(mat,row+nalg,col+nalg); wp = (Q **)BDY(mat);
                   1947:                        for ( i = 0; i < row; i++ )
                   1948:                                for ( j = 0; j < col; j++ ) {
                   1949:                                        w[i][j] = QTOS(b[i][j]);
                   1950:                                        wp[i][j] = b[i][j];
                   1951:                                }
                   1952:                        for ( i = 0; i < nalg; i++ ) {
                   1953:                                w[i+row][i+col] = 1;
                   1954:                                wp[i+row][i+col] = ONE;
                   1955:                        }
                   1956:                        r->id = 2; r->obj = (Obj)mat;
                   1957:                        r->nv = col+nalg; r->ord.matrix.row = row+nalg;
                   1958:                        r->ord.matrix.matrix = w;
                   1959:                        break;
                   1960:                case 3:
                   1961:                default:
                   1962:                        /* XXX */
                   1963:                        error("append_block : not implemented yet");
                   1964:        }
                   1965:        return r;
1.28      noro     1966: }
                   1967:
1.37      noro     1968: int comp_sw(struct sparse_weight *a, struct sparse_weight *b)
                   1969: {
                   1970:        if ( a->pos > b->pos ) return 1;
                   1971:        else if ( a->pos < b->pos ) return -1;
                   1972:        else return 0;
                   1973: }
                   1974:
1.27      noro     1975: /* order = [w_or_b, w_or_b, ... ] */
                   1976: /* w_or_b = w or b                */
                   1977: /* w = [1,2,...] or [x,1,y,2,...] */
                   1978: /* b = [@lex,x,y,...,z] etc       */
                   1979:
                   1980: int create_composite_order_spec(VL vl,LIST order,struct order_spec **specp)
                   1981: {
                   1982:        NODE wb,t,p;
                   1983:        struct order_spec *spec;
                   1984:        VL tvl;
1.29      noro     1985:        int n,i,j,k,l,start,end,len,w;
1.27      noro     1986:        int *dw;
                   1987:        struct sparse_weight *sw;
                   1988:        struct weight_or_block *w_or_b;
                   1989:        Obj a0;
                   1990:        NODE a;
1.29      noro     1991:        V v,sv,ev;
                   1992:        SYMBOL sym;
                   1993:        int *top;
1.27      noro     1994:
                   1995:        /* l = number of vars in vl */
                   1996:        for ( l = 0, tvl = vl; tvl; tvl = NEXT(tvl), l++ );
                   1997:        /* n = number of primitives in order */
                   1998:        wb = BDY(order);
                   1999:        n = length(wb);
                   2000:        *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
                   2001:        spec->id = 3;
                   2002:        spec->obj = (Obj)order;
                   2003:        spec->nv = l;
                   2004:        spec->ord.composite.length = n;
1.28      noro     2005:        w_or_b = spec->ord.composite.w_or_b = (struct weight_or_block *)
1.29      noro     2006:                MALLOC(sizeof(struct weight_or_block)*(n+1));
                   2007:
                   2008:        /* top : register the top variable in each w_or_b specification */
                   2009:        top = (int *)ALLOCA(l*sizeof(int));
                   2010:        for ( i = 0; i < l; i++ ) top[i] = 0;
                   2011:
1.28      noro     2012:        for ( t = wb, i = 0; t; t = NEXT(t), i++ ) {
1.30      noro     2013:                if ( !BDY(t) || OID((Obj)BDY(t)) != O_LIST )
                   2014:                        error("a list of lists must be specified for the key \"order\"");
1.28      noro     2015:                a = BDY((LIST)BDY(t));
1.27      noro     2016:                len = length(a);
                   2017:                a0 = (Obj)BDY(a);
                   2018:                if ( !a0 || OID(a0) == O_N ) {
1.28      noro     2019:                        /* a is a dense weight vector */
1.27      noro     2020:                        dw = (int *)MALLOC(sizeof(int)*len);
1.30      noro     2021:                        for ( j = 0, p = a; j < len; p = NEXT(p), j++ ) {
                   2022:                                if ( !INT((Q)BDY(p)) )
                   2023:                                        error("a dense weight vector must be specified as a list of integers");
1.27      noro     2024:                                dw[j] = QTOS((Q)BDY(p));
1.30      noro     2025:                        }
1.27      noro     2026:                        w_or_b[i].type = IS_DENSE_WEIGHT;
                   2027:                        w_or_b[i].length = len;
                   2028:                        w_or_b[i].body.dense_weight = dw;
1.29      noro     2029:
                   2030:                        /* find the top */
                   2031:                        for ( k = 0; k < len && !dw[k]; k++ );
                   2032:                        if ( k < len ) top[k] = 1;
                   2033:
1.27      noro     2034:                } else if ( OID(a0) == O_P ) {
1.28      noro     2035:                        /* a is a sparse weight vector */
                   2036:                        len >>= 1;
1.27      noro     2037:                        sw = (struct sparse_weight *)
                   2038:                                MALLOC(sizeof(struct sparse_weight)*len);
                   2039:                        for ( j = 0, p = a; j < len; j++ ) {
1.30      noro     2040:                                if ( !BDY(p) || OID((P)BDY(p)) != O_P )
                   2041:                                        error("a sparse weight vector must be specified as [var1,weight1,...]");
1.28      noro     2042:                                v = VR((P)BDY(p)); p = NEXT(p);
1.27      noro     2043:                                for ( tvl = vl, k = 0; tvl && tvl->v != v;
                   2044:                                        k++, tvl = NEXT(tvl) );
                   2045:                                if ( !tvl )
1.30      noro     2046:                                        error("invalid variable name in a sparse weight vector");
1.27      noro     2047:                                sw[j].pos = k;
1.30      noro     2048:                                if ( !INT((Q)BDY(p)) )
                   2049:                                        error("a sparse weight vector must be specified as [var1,weight1,...]");
1.28      noro     2050:                                sw[j].value = QTOS((Q)BDY(p)); p = NEXT(p);
1.27      noro     2051:                        }
1.37      noro     2052:                        qsort(sw,len,sizeof(struct sparse_weight),
                   2053:                                (int (*)(const void *,const void *))comp_sw);
1.27      noro     2054:                        w_or_b[i].type = IS_SPARSE_WEIGHT;
                   2055:                        w_or_b[i].length = len;
                   2056:                        w_or_b[i].body.sparse_weight = sw;
1.29      noro     2057:
                   2058:                        /* find the top */
                   2059:                        for ( k = 0; k < len && !sw[k].value; k++ );
                   2060:                        if ( k < len ) top[sw[k].pos] = 1;
                   2061:                } else if ( OID(a0) == O_RANGE ) {
                   2062:                        /* [range(v1,v2),w] */
                   2063:                        sv = VR((P)(((RANGE)a0)->start));
                   2064:                        ev = VR((P)(((RANGE)a0)->end));
                   2065:                        for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
                   2066:                        if ( !tvl )
                   2067:                                error("invalid range");
                   2068:                        for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
                   2069:                        if ( !tvl )
                   2070:                                error("invalid range");
                   2071:                        len = end-start+1;
                   2072:                        sw = (struct sparse_weight *)
                   2073:                                MALLOC(sizeof(struct sparse_weight)*len);
                   2074:                        w = QTOS((Q)BDY(NEXT(a)));
                   2075:                        for ( tvl = vl, k = 0; k < start; k++, tvl = NEXT(tvl) );
                   2076:                        for ( j = 0 ; k <= end; k++, tvl = NEXT(tvl), j++ ) {
                   2077:                                sw[j].pos = k;
                   2078:                                sw[j].value = w;
                   2079:                        }
                   2080:                        w_or_b[i].type = IS_SPARSE_WEIGHT;
                   2081:                        w_or_b[i].length = len;
                   2082:                        w_or_b[i].body.sparse_weight = sw;
                   2083:
                   2084:                        /* register the top */
                   2085:                        if ( w ) top[start] = 1;
1.28      noro     2086:                } else if ( OID(a0) == O_SYMBOL ) {
                   2087:                        /* a is a block */
1.29      noro     2088:                        sym = (SYMBOL)a0; a = NEXT(a); len--;
                   2089:                        if ( OID((Obj)BDY(a)) == O_RANGE ) {
                   2090:                                sv = VR((P)(((RANGE)BDY(a))->start));
                   2091:                                ev = VR((P)(((RANGE)BDY(a))->end));
                   2092:                                for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
                   2093:                                if ( !tvl )
                   2094:                                        error("invalid range");
                   2095:                                for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
                   2096:                                if ( !tvl )
                   2097:                                        error("invalid range");
                   2098:                                len = end-start+1;
                   2099:                        } else {
                   2100:                                for ( start = 0, tvl = vl; tvl->v != VR((P)BDY(a));
1.28      noro     2101:                                tvl = NEXT(tvl), start++ );
1.29      noro     2102:                                for ( p = NEXT(a), tvl = NEXT(tvl); p;
1.30      noro     2103:                                        p = NEXT(p), tvl = NEXT(tvl) ) {
                   2104:                                        if ( !BDY(p) || OID((P)BDY(p)) != O_P )
                   2105:                                                error("a block must be specified as [ordsymbol,var1,var2,...]");
1.29      noro     2106:                                        if ( tvl->v != VR((P)BDY(p)) ) break;
1.30      noro     2107:                                }
1.29      noro     2108:                                if ( p )
1.30      noro     2109:                                        error("a block must be contiguous in the variable list");
1.29      noro     2110:                        }
1.28      noro     2111:                        w_or_b[i].type = IS_BLOCK;
                   2112:                        w_or_b[i].length = len;
                   2113:                        w_or_b[i].body.block.start = start;
                   2114:                        if ( !strcmp(sym->name,"@grlex") )
                   2115:                                w_or_b[i].body.block.order = 0;
                   2116:                        else if ( !strcmp(sym->name,"@glex") )
                   2117:                                w_or_b[i].body.block.order = 1;
                   2118:                        else if ( !strcmp(sym->name,"@lex") )
                   2119:                                w_or_b[i].body.block.order = 2;
                   2120:                        else
1.29      noro     2121:                                error("invalid ordername");
                   2122:                        /* register the tops */
                   2123:                        for ( j = 0, k = start; j < len; j++, k++ )
                   2124:                                top[k] = 1;
1.28      noro     2125:                }
1.29      noro     2126:        }
                   2127:        for ( k = 0; k < l && top[k]; k++ );
                   2128:        if ( k < l ) {
                   2129:                /* incomplete order specification; add @grlex */
                   2130:                w_or_b[n].type = IS_BLOCK;
                   2131:                w_or_b[n].length = l;
                   2132:                w_or_b[n].body.block.start = 0;
                   2133:                w_or_b[n].body.block.order = 0;
                   2134:                spec->ord.composite.length = n+1;
1.27      noro     2135:        }
                   2136: }
                   2137:
1.35      noro     2138: /* module order spec */
                   2139:
                   2140: void create_modorder_spec(int id,LIST shift,struct modorder_spec **s)
                   2141: {
                   2142:        struct modorder_spec *spec;
                   2143:        NODE n,t;
                   2144:        LIST list;
                   2145:        int *ds;
                   2146:        int i,l;
                   2147:        Q q;
                   2148:
                   2149:        *s = spec = (struct modorder_spec *)MALLOC(sizeof(struct modorder_spec));
                   2150:        spec->id = id;
                   2151:        if ( shift ) {
                   2152:                n = BDY(shift);
                   2153:                spec->len = l = length(n);
                   2154:                spec->degree_shift = ds = (int *)MALLOC_ATOMIC(l*sizeof(int));
                   2155:                for ( t = n, i = 0; t; t = NEXT(t), i++ )
                   2156:                        ds[i] = QTOS((Q)BDY(t));
                   2157:        } else {
                   2158:                spec->len = 0;
                   2159:                spec->degree_shift = 0;
                   2160:        }
                   2161:        STOQ(id,q);
                   2162:        n = mknode(2,q,shift);
                   2163:        MKLIST(list,n);
                   2164:        spec->obj = (Obj)list;
                   2165: }
                   2166:
1.7       noro     2167: /*
                   2168:  * converters
                   2169:  *
                   2170:  */
                   2171:
1.20      noro     2172: void dp_homo(DP p,DP *rp)
1.5       noro     2173: {
1.7       noro     2174:        MP m,mr,mr0;
                   2175:        int i,n,nv,td;
                   2176:        DL dl,dlh;
1.5       noro     2177:
1.7       noro     2178:        if ( !p )
                   2179:                *rp = 0;
                   2180:        else {
                   2181:                n = p->nv; nv = n + 1;
                   2182:                m = BDY(p); td = sugard(m);
                   2183:                for ( mr0 = 0; m; m = NEXT(m) ) {
                   2184:                        NEXTMP(mr0,mr); mr->c = m->c;
                   2185:                        dl = m->dl;
                   2186:                        mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
                   2187:                        dlh->td = td;
                   2188:                        for ( i = 0; i < n; i++ )
                   2189:                                dlh->d[i] = dl->d[i];
                   2190:                        dlh->d[n] = td - dl->td;
                   2191:                }
                   2192:                NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5       noro     2193:        }
                   2194: }
                   2195:
1.20      noro     2196: void dp_dehomo(DP p,DP *rp)
1.5       noro     2197: {
1.7       noro     2198:        MP m,mr,mr0;
                   2199:        int i,n,nv;
                   2200:        DL dl,dlh;
1.5       noro     2201:
1.7       noro     2202:        if ( !p )
                   2203:                *rp = 0;
                   2204:        else {
                   2205:                n = p->nv; nv = n - 1;
                   2206:                m = BDY(p);
                   2207:                for ( mr0 = 0; m; m = NEXT(m) ) {
                   2208:                        NEXTMP(mr0,mr); mr->c = m->c;
                   2209:                        dlh = m->dl;
                   2210:                        mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
                   2211:                        dl->td = dlh->td - dlh->d[nv];
                   2212:                        for ( i = 0; i < nv; i++ )
                   2213:                                dl->d[i] = dlh->d[i];
                   2214:                }
                   2215:                NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
                   2216:        }
1.5       noro     2217: }
                   2218:
1.20      noro     2219: void dp_mod(DP p,int mod,NODE subst,DP *rp)
1.5       noro     2220: {
1.7       noro     2221:        MP m,mr,mr0;
                   2222:        P t,s,s1;
                   2223:        V v;
                   2224:        NODE tn;
1.5       noro     2225:
1.7       noro     2226:        if ( !p )
                   2227:                *rp = 0;
                   2228:        else {
                   2229:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   2230:                        for ( tn = subst, s = m->c; tn; tn = NEXT(tn) ) {
                   2231:                                v = VR((P)BDY(tn)); tn = NEXT(tn);
                   2232:                                substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
                   2233:                        }
                   2234:                        ptomp(mod,s,&t);
                   2235:                        if ( t ) {
                   2236:                                NEXTMP(mr0,mr); mr->c = t; mr->dl = m->dl;
                   2237:                        }
                   2238:                }
                   2239:                if ( mr0 ) {
                   2240:                        NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                   2241:                } else
                   2242:                        *rp = 0;
                   2243:        }
1.5       noro     2244: }
                   2245:
1.20      noro     2246: void dp_rat(DP p,DP *rp)
1.5       noro     2247: {
1.7       noro     2248:        MP m,mr,mr0;
1.5       noro     2249:
1.7       noro     2250:        if ( !p )
                   2251:                *rp = 0;
                   2252:        else {
                   2253:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   2254:                        NEXTMP(mr0,mr); mptop(m->c,&mr->c); mr->dl = m->dl;
1.5       noro     2255:                }
1.7       noro     2256:                if ( mr0 ) {
                   2257:                        NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                   2258:                } else
                   2259:                        *rp = 0;
1.5       noro     2260:        }
                   2261: }
                   2262:
                   2263:
1.27      noro     2264: void homogenize_order(struct order_spec *old,int n,struct order_spec **newp)
1.5       noro     2265: {
1.7       noro     2266:        struct order_pair *l;
                   2267:        int length,nv,row,i,j;
                   2268:        int **newm,**oldm;
1.27      noro     2269:        struct order_spec *new;
1.31      noro     2270:        int onv,nnv,nlen,olen,owlen;
                   2271:        struct weight_or_block *owb,*nwb;
1.5       noro     2272:
1.27      noro     2273:        *newp = new = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1.7       noro     2274:        switch ( old->id ) {
                   2275:                case 0:
                   2276:                        switch ( old->ord.simple ) {
                   2277:                                case 0:
                   2278:                                        new->id = 0; new->ord.simple = 0; break;
                   2279:                                case 1:
                   2280:                                        l = (struct order_pair *)
                   2281:                                                MALLOC_ATOMIC(2*sizeof(struct order_pair));
                   2282:                                        l[0].length = n; l[0].order = 1;
                   2283:                                        l[1].length = 1; l[1].order = 2;
                   2284:                                        new->id = 1;
                   2285:                                        new->ord.block.order_pair = l;
                   2286:                                        new->ord.block.length = 2; new->nv = n+1;
                   2287:                                        break;
                   2288:                                case 2:
                   2289:                                        new->id = 0; new->ord.simple = 1; break;
                   2290:                                case 3: case 4: case 5:
                   2291:                                        new->id = 0; new->ord.simple = old->ord.simple+3;
                   2292:                                        dp_nelim = n-1; break;
                   2293:                                case 6: case 7: case 8: case 9:
                   2294:                                        new->id = 0; new->ord.simple = old->ord.simple; break;
                   2295:                                default:
                   2296:                                        error("homogenize_order : invalid input");
                   2297:                        }
                   2298:                        break;
1.50      noro     2299:                case 1: case 257:
1.7       noro     2300:                        length = old->ord.block.length;
                   2301:                        l = (struct order_pair *)
                   2302:                                MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
                   2303:                        bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
                   2304:                        l[length].order = 2; l[length].length = 1;
1.50      noro     2305:                        new->id = old->id; new->nv = n+1;
1.7       noro     2306:                        new->ord.block.order_pair = l;
                   2307:                        new->ord.block.length = length+1;
1.51      noro     2308:                        new->ispot = old->ispot;
1.7       noro     2309:                        break;
1.50      noro     2310:                case 2: case 258:
1.7       noro     2311:                        nv = old->nv; row = old->ord.matrix.row;
                   2312:                        oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
                   2313:                        for ( i = 0; i <= nv; i++ )
                   2314:                                newm[0][i] = 1;
                   2315:                        for ( i = 0; i < row; i++ ) {
                   2316:                                for ( j = 0; j < nv; j++ )
                   2317:                                        newm[i+1][j] = oldm[i][j];
                   2318:                                newm[i+1][j] = 0;
                   2319:                        }
1.50      noro     2320:                        new->id = old->id; new->nv = nv+1;
1.7       noro     2321:                        new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
1.51      noro     2322:                        new->ispot = old->ispot;
1.31      noro     2323:                        break;
1.50      noro     2324:                case 3: case 259:
1.31      noro     2325:                        onv = old->nv;
                   2326:                        nnv = onv+1;
                   2327:                        olen = old->ord.composite.length;
                   2328:                        nlen = olen+1;
                   2329:                        owb = old->ord.composite.w_or_b;
                   2330:                        nwb = (struct weight_or_block *)
                   2331:                                MALLOC(nlen*sizeof(struct weight_or_block));
                   2332:                        for ( i = 0; i < olen; i++ ) {
                   2333:                                nwb[i].type = owb[i].type;
                   2334:                                switch ( owb[i].type ) {
                   2335:                                        case IS_DENSE_WEIGHT:
                   2336:                                                owlen = owb[i].length;
                   2337:                                                nwb[i].length = owlen+1;
                   2338:                                                nwb[i].body.dense_weight = (int *)MALLOC((owlen+1)*sizeof(int));
                   2339:                                                for ( j = 0; j < owlen; j++ )
                   2340:                                                        nwb[i].body.dense_weight[j] = owb[i].body.dense_weight[j];
                   2341:                                                nwb[i].body.dense_weight[owlen] = 0;
                   2342:                                                break;
                   2343:                                        case IS_SPARSE_WEIGHT:
                   2344:                                                nwb[i].length = owb[i].length;
                   2345:                                                nwb[i].body.sparse_weight = owb[i].body.sparse_weight;
                   2346:                                                break;
                   2347:                                        case IS_BLOCK:
                   2348:                                                nwb[i].length = owb[i].length;
                   2349:                                                nwb[i].body.block = owb[i].body.block;
                   2350:                                                break;
                   2351:                                }
                   2352:                        }
                   2353:                        nwb[i].type = IS_SPARSE_WEIGHT;
                   2354:                        nwb[i].body.sparse_weight =
                   2355:                                (struct sparse_weight *)MALLOC(sizeof(struct sparse_weight));
                   2356:                        nwb[i].body.sparse_weight[0].pos = onv;
                   2357:                        nwb[i].body.sparse_weight[0].value = 1;
1.50      noro     2358:                        new->id = old->id;
1.31      noro     2359:                        new->nv = nnv;
                   2360:                        new->ord.composite.length = nlen;
                   2361:                        new->ord.composite.w_or_b = nwb;
1.51      noro     2362:                        new->ispot = old->ispot;
1.31      noro     2363:                        print_composite_order_spec(new);
1.7       noro     2364:                        break;
1.50      noro     2365:                case 256: /* simple module order */
                   2366:                        switch ( old->ord.simple ) {
                   2367:                                case 0:
                   2368:                                        new->id = 256; new->ord.simple = 0; break;
                   2369:                                case 1:
                   2370:                                        l = (struct order_pair *)
                   2371:                                                MALLOC_ATOMIC(2*sizeof(struct order_pair));
                   2372:                                        l[0].length = n; l[0].order = old->ord.simple;
                   2373:                                        l[1].length = 1; l[1].order = 2;
                   2374:                                        new->id = 257;
                   2375:                                        new->ord.block.order_pair = l;
                   2376:                                        new->ord.block.length = 2; new->nv = n+1;
                   2377:                                        break;
                   2378:                                case 2:
                   2379:                                        new->id = 256; new->ord.simple = 1; break;
                   2380:                                default:
                   2381:                                        error("homogenize_order : invalid input");
                   2382:                        }
1.51      noro     2383:                        new->ispot = old->ispot;
1.50      noro     2384:                        break;
1.7       noro     2385:                default:
                   2386:                        error("homogenize_order : invalid input");
1.5       noro     2387:        }
1.7       noro     2388: }
                   2389:
1.20      noro     2390: void qltozl(Q *w,int n,Q *dvr)
1.7       noro     2391: {
                   2392:        N nm,dn;
                   2393:        N g,l1,l2,l3;
                   2394:        Q c,d;
                   2395:        int i;
                   2396:        struct oVECT v;
1.5       noro     2397:
                   2398:        for ( i = 0; i < n; i++ )
1.7       noro     2399:                if ( w[i] && !INT(w[i]) )
                   2400:                        break;
                   2401:        if ( i == n ) {
                   2402:                v.id = O_VECT; v.len = n; v.body = (pointer *)w;
                   2403:                igcdv(&v,dvr); return;
                   2404:        }
1.56      noro     2405:        for ( i = 0; !w[i]; i++ );
                   2406:        c = w[i]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
                   2407:        for ( i++; i < n; i++ ) {
                   2408:                c = w[i];
                   2409:                if ( !c ) continue;
                   2410:                l1 = INT(c) ? ONEN : DN(c);
1.7       noro     2411:                gcdn(nm,NM(c),&g); nm = g;
                   2412:                gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
1.5       noro     2413:        }
1.7       noro     2414:        if ( UNIN(dn) )
                   2415:                NTOQ(nm,1,d);
                   2416:        else
                   2417:                NDTOQ(nm,dn,1,d);
                   2418:        *dvr = d;
                   2419: }
1.5       noro     2420:
1.20      noro     2421: int comp_nm(Q *a,Q *b)
1.7       noro     2422: {
                   2423:        return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
                   2424: }
                   2425:
1.20      noro     2426: void sortbynm(Q *w,int n)
1.7       noro     2427: {
                   2428:        qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
                   2429: }
1.5       noro     2430:
                   2431:
1.7       noro     2432: /*
                   2433:  * simple operations
                   2434:  *
                   2435:  */
1.5       noro     2436:
1.20      noro     2437: int dp_redble(DP p1,DP p2)
1.7       noro     2438: {
                   2439:        int i,n;
                   2440:        DL d1,d2;
1.5       noro     2441:
1.7       noro     2442:        d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                   2443:        if ( d1->td < d2->td )
                   2444:                return 0;
                   2445:        else {
                   2446:                for ( i = 0, n = p1->nv; i < n; i++ )
                   2447:                        if ( d1->d[i] < d2->d[i] )
                   2448:                                return 0;
                   2449:                return 1;
1.5       noro     2450:        }
                   2451: }
                   2452:
1.20      noro     2453: void dp_subd(DP p1,DP p2,DP *rp)
1.5       noro     2454: {
1.7       noro     2455:        int i,n;
1.5       noro     2456:        DL d1,d2,d;
                   2457:        MP m;
1.7       noro     2458:        DP s;
1.5       noro     2459:
                   2460:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.7       noro     2461:        NEWDL(d,n); d->td = d1->td - d2->td;
1.5       noro     2462:        for ( i = 0; i < n; i++ )
1.7       noro     2463:                d->d[i] = d1->d[i]-d2->d[i];
                   2464:        NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
                   2465:        MKDP(n,m,s); s->sugar = d->td;
                   2466:        *rp = s;
                   2467: }
                   2468:
1.20      noro     2469: void dltod(DL d,int n,DP *rp)
1.7       noro     2470: {
                   2471:        MP m;
                   2472:        DP s;
                   2473:
                   2474:        NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
                   2475:        MKDP(n,m,s); s->sugar = d->td;
                   2476:        *rp = s;
1.5       noro     2477: }
                   2478:
1.20      noro     2479: void dp_hm(DP p,DP *rp)
1.5       noro     2480: {
                   2481:        MP m,mr;
                   2482:
                   2483:        if ( !p )
                   2484:                *rp = 0;
                   2485:        else {
                   2486:                m = BDY(p);
                   2487:                NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
                   2488:                MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td;  /* XXX */
                   2489:        }
                   2490: }
                   2491:
1.35      noro     2492: void dp_ht(DP p,DP *rp)
                   2493: {
                   2494:        MP m,mr;
                   2495:
                   2496:        if ( !p )
                   2497:                *rp = 0;
                   2498:        else {
                   2499:                m = BDY(p);
                   2500:                NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
                   2501:                MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td;  /* XXX */
                   2502:        }
                   2503: }
                   2504:
1.20      noro     2505: void dp_rest(DP p,DP *rp)
1.5       noro     2506: {
                   2507:        MP m;
                   2508:
                   2509:        m = BDY(p);
                   2510:        if ( !NEXT(m) )
                   2511:                *rp = 0;
                   2512:        else {
                   2513:                MKDP(p->nv,NEXT(m),*rp);
                   2514:                if ( *rp )
                   2515:                        (*rp)->sugar = p->sugar;
                   2516:        }
                   2517: }
                   2518:
1.20      noro     2519: DL lcm_of_DL(int nv,DL dl1,DL dl2,DL dl)
1.5       noro     2520: {
1.21      noro     2521:        register int i, *d1, *d2, *d, td;
1.5       noro     2522:
                   2523:        if ( !dl ) NEWDL(dl,nv);
                   2524:        d = dl->d,  d1 = dl1->d,  d2 = dl2->d;
1.21      noro     2525:        for ( td = 0, i = 0; i < nv; d1++, d2++, d++, i++ ) {
                   2526:                *d = *d1 > *d2 ? *d1 : *d2;
                   2527:                td += MUL_WEIGHT(*d,i);
                   2528:        }
1.5       noro     2529:        dl->td = td;
                   2530:        return dl;
                   2531: }
                   2532:
1.20      noro     2533: int dl_equal(int nv,DL dl1,DL dl2)
1.5       noro     2534: {
                   2535:     register int *d1, *d2, n;
                   2536:
                   2537:     if ( dl1->td != dl2->td ) return 0;
                   2538:     for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
                   2539:         if ( *d1 != *d2 ) return 0;
                   2540:     return 1;
                   2541: }
                   2542:
1.20      noro     2543: int dp_nt(DP p)
1.5       noro     2544: {
                   2545:        int i;
                   2546:        MP m;
                   2547:
                   2548:        if ( !p )
                   2549:                return 0;
                   2550:        else {
                   2551:                for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
                   2552:                return i;
                   2553:        }
                   2554: }
                   2555:
1.20      noro     2556: int dp_homogeneous(DP p)
1.15      noro     2557: {
                   2558:        MP m;
                   2559:        int d;
                   2560:
                   2561:        if ( !p )
                   2562:                return 1;
                   2563:        else {
                   2564:                m = BDY(p);
                   2565:                d = m->dl->td;
                   2566:                m = NEXT(m);
                   2567:                for ( ; m; m = NEXT(m) ) {
                   2568:                        if ( m->dl->td != d )
                   2569:                                return 0;
                   2570:                }
                   2571:                return 1;
                   2572:        }
1.16      noro     2573: }
                   2574:
1.20      noro     2575: void _print_mp(int nv,MP m)
1.16      noro     2576: {
                   2577:        int i;
                   2578:
1.17      noro     2579:        if ( !m )
1.16      noro     2580:                return;
                   2581:        for ( ; m; m = NEXT(m) ) {
                   2582:                fprintf(stderr,"%d<",ITOS(C(m)));
                   2583:                for ( i = 0; i < nv; i++ ) {
                   2584:                        fprintf(stderr,"%d",m->dl->d[i]);
                   2585:                        if ( i != nv-1 )
                   2586:                                fprintf(stderr," ");
                   2587:                }
                   2588:                fprintf(stderr,">",C(m));
                   2589:        }
                   2590:        fprintf(stderr,"\n");
1.15      noro     2591: }
1.26      noro     2592:
                   2593: static int cmp_mp_nvar;
                   2594:
                   2595: int comp_mp(MP *a,MP *b)
                   2596: {
                   2597:        return -(*cmpdl)(cmp_mp_nvar,(*a)->dl,(*b)->dl);
                   2598: }
                   2599:
                   2600: void dp_sort(DP p,DP *rp)
                   2601: {
                   2602:        MP t,mp,mp0;
                   2603:        int i,n;
                   2604:        DP r;
                   2605:        MP *w;
                   2606:
                   2607:        if ( !p ) {
                   2608:                *rp = 0;
                   2609:                return;
                   2610:        }
                   2611:        for ( t = BDY(p), n = 0; t; t = NEXT(t), n++ );
                   2612:        w = (MP *)ALLOCA(n*sizeof(MP));
                   2613:        for ( t = BDY(p), i = 0; i < n; t = NEXT(t), i++ )
                   2614:                w[i] = t;
                   2615:        cmp_mp_nvar = NV(p);
                   2616:        qsort(w,n,sizeof(MP),(int (*)(const void *,const void *))comp_mp);
                   2617:        mp0 = 0;
                   2618:        for ( i = n-1; i >= 0; i-- ) {
                   2619:                NEWMP(mp); mp->dl = w[i]->dl; C(mp) = C(w[i]);
                   2620:                NEXT(mp) = mp0; mp0 = mp;
                   2621:        }
                   2622:        MKDP(p->nv,mp0,r);
                   2623:        r->sugar = p->sugar;
                   2624:        *rp = r;
                   2625: }
                   2626:
1.32      noro     2627: DP extract_initial_term_from_dp(DP p,int *weight,int n);
                   2628: LIST extract_initial_term(LIST f,int *weight,int n);
                   2629:
                   2630: DP extract_initial_term_from_dp(DP p,int *weight,int n)
                   2631: {
1.34      noro     2632:        int w,t,i,top;
1.32      noro     2633:        MP m,r0,r;
                   2634:        DP dp;
                   2635:
                   2636:        if ( !p ) return 0;
1.34      noro     2637:        top = 1;
1.32      noro     2638:        for ( m = BDY(p); m; m = NEXT(m) ) {
                   2639:                for ( i = 0, t = 0; i < n; i++ )
                   2640:                        t += weight[i]*m->dl->d[i];
1.34      noro     2641:                if ( top || t > w ) {
1.32      noro     2642:                        r0 = 0;
                   2643:                        w = t;
1.34      noro     2644:                        top = 0;
1.32      noro     2645:                }
                   2646:                if ( t == w ) {
                   2647:                        NEXTMP(r0,r);
                   2648:                        r->dl = m->dl;
                   2649:                        r->c = m->c;
                   2650:                }
                   2651:        }
                   2652:        NEXT(r) = 0;
                   2653:        MKDP(p->nv,r0,dp);
                   2654:        return dp;
                   2655: }
                   2656:
                   2657: LIST extract_initial_term(LIST f,int *weight,int n)
                   2658: {
                   2659:        NODE nd,r0,r;
                   2660:        Obj p;
                   2661:        LIST l;
                   2662:
                   2663:        nd = BDY(f);
                   2664:        for ( r0 = 0; nd; nd = NEXT(nd) ) {
                   2665:                NEXTNODE(r0,r);
                   2666:                p = (Obj)BDY(nd);
                   2667:                BDY(r) = (pointer)extract_initial_term_from_dp((DP)p,weight,n);
                   2668:        }
                   2669:        if ( r0 ) NEXT(r) = 0;
                   2670:        MKLIST(l,r0);
                   2671:        return l;
                   2672: }
                   2673:
                   2674: LIST dp_initial_term(LIST f,struct order_spec *ord)
                   2675: {
                   2676:        int n,l,i;
                   2677:        struct weight_or_block *worb;
                   2678:        int *weight;
                   2679:
                   2680:        switch ( ord->id ) {
                   2681:                case 2: /* matrix order */
                   2682:                        /* extract the first row */
                   2683:                        n = ord->nv;
                   2684:                        weight = ord->ord.matrix.matrix[0];
                   2685:                        return extract_initial_term(f,weight,n);
                   2686:                case 3: /* composite order */
                   2687:                        /* the first w_or_b */
                   2688:                        worb = ord->ord.composite.w_or_b;
                   2689:                        switch ( worb->type ) {
                   2690:                                case IS_DENSE_WEIGHT:
                   2691:                                        n = worb->length;
                   2692:                                        weight = worb->body.dense_weight;
                   2693:                                        return extract_initial_term(f,weight,n);
                   2694:                                case IS_SPARSE_WEIGHT:
                   2695:                                        n = ord->nv;
                   2696:                                        weight = (int *)ALLOCA(n*sizeof(int));
1.33      noro     2697:                                        for ( i = 0; i < n; i++ ) weight[i] = 0;
1.32      noro     2698:                                        l = worb->length;
                   2699:                                        for ( i = 0; i < l; i++ )
                   2700:                                                weight[worb->body.sparse_weight[i].pos]
                   2701:                                                        =  worb->body.sparse_weight[i].value;
                   2702:                                        return extract_initial_term(f,weight,n);
                   2703:                                default:
                   2704:                                        error("dp_initial_term : unsupported order");
                   2705:                        }
                   2706:                default:
                   2707:                        error("dp_initial_term : unsupported order");
                   2708:        }
                   2709: }
                   2710:
                   2711: int highest_order_dp(DP p,int *weight,int n);
                   2712: LIST highest_order(LIST f,int *weight,int n);
                   2713:
                   2714: int highest_order_dp(DP p,int *weight,int n)
                   2715: {
1.34      noro     2716:        int w,t,i,top;
1.32      noro     2717:        MP m;
                   2718:
                   2719:        if ( !p ) return -1;
1.34      noro     2720:        top = 1;
1.32      noro     2721:        for ( m = BDY(p); m; m = NEXT(m) ) {
                   2722:                for ( i = 0, t = 0; i < n; i++ )
                   2723:                        t += weight[i]*m->dl->d[i];
1.34      noro     2724:                if ( top || t > w ) {
1.32      noro     2725:                        w = t;
1.34      noro     2726:                        top = 0;
                   2727:                }
1.32      noro     2728:        }
                   2729:        return w;
                   2730: }
                   2731:
                   2732: LIST highest_order(LIST f,int *weight,int n)
                   2733: {
                   2734:        int h;
                   2735:        NODE nd,r0,r;
                   2736:        Obj p;
                   2737:        LIST l;
                   2738:        Q q;
                   2739:
                   2740:        nd = BDY(f);
                   2741:        for ( r0 = 0; nd; nd = NEXT(nd) ) {
                   2742:                NEXTNODE(r0,r);
                   2743:                p = (Obj)BDY(nd);
                   2744:                h = highest_order_dp((DP)p,weight,n);
                   2745:                STOQ(h,q);
                   2746:                BDY(r) = (pointer)q;
                   2747:        }
                   2748:        if ( r0 ) NEXT(r) = 0;
                   2749:        MKLIST(l,r0);
                   2750:        return l;
                   2751: }
                   2752:
                   2753: LIST dp_order(LIST f,struct order_spec *ord)
                   2754: {
                   2755:        int n,l,i;
                   2756:        struct weight_or_block *worb;
                   2757:        int *weight;
                   2758:
                   2759:        switch ( ord->id ) {
                   2760:                case 2: /* matrix order */
                   2761:                        /* extract the first row */
                   2762:                        n = ord->nv;
                   2763:                        weight = ord->ord.matrix.matrix[0];
                   2764:                        return highest_order(f,weight,n);
                   2765:                case 3: /* composite order */
                   2766:                        /* the first w_or_b */
                   2767:                        worb = ord->ord.composite.w_or_b;
                   2768:                        switch ( worb->type ) {
                   2769:                                case IS_DENSE_WEIGHT:
                   2770:                                        n = worb->length;
                   2771:                                        weight = worb->body.dense_weight;
                   2772:                                        return highest_order(f,weight,n);
                   2773:                                case IS_SPARSE_WEIGHT:
                   2774:                                        n = ord->nv;
                   2775:                                        weight = (int *)ALLOCA(n*sizeof(int));
1.33      noro     2776:                                        for ( i = 0; i < n; i++ ) weight[i] = 0;
1.32      noro     2777:                                        l = worb->length;
                   2778:                                        for ( i = 0; i < l; i++ )
                   2779:                                                weight[worb->body.sparse_weight[i].pos]
                   2780:                                                        =  worb->body.sparse_weight[i].value;
                   2781:                                        return highest_order(f,weight,n);
                   2782:                                default:
                   2783:                                        error("dp_initial_term : unsupported order");
                   2784:                        }
                   2785:                default:
                   2786:                        error("dp_initial_term : unsupported order");
1.35      noro     2787:        }
                   2788: }
                   2789:
                   2790: int dpv_ht(DPV p,DP *h)
                   2791: {
                   2792:        int len,max,maxi,i,t;
                   2793:        DP *e;
                   2794:        MP m,mr;
                   2795:
                   2796:        len = p->len;
                   2797:        e = p->body;
                   2798:        max = -1;
                   2799:        maxi = -1;
                   2800:        for ( i = 0; i < len; i++ )
                   2801:                if ( e[i] && (t = BDY(e[i])->dl->td) > max ) {
                   2802:                        max = t;
                   2803:                        maxi = i;
                   2804:                }
                   2805:        if ( max < 0 ) {
                   2806:                *h = 0;
                   2807:                return -1;
                   2808:        } else {
                   2809:                m = BDY(e[maxi]);
                   2810:                NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
                   2811:                MKDP(e[maxi]->nv,mr,*h); (*h)->sugar = mr->dl->td;  /* XXX */
                   2812:                return maxi;
1.32      noro     2813:        }
                   2814: }
1.42      noro     2815:
                   2816: /* return 1 if 0 <_w1 v && v <_w2 0 */
                   2817:
                   2818: int in_c12(int n,int *v,int row1,int **w1,int row2, int **w2)
                   2819: {
                   2820:        int t1,t2;
                   2821:
                   2822:        t1 = compare_zero(n,v,row1,w1);
                   2823:        t2 = compare_zero(n,v,row2,w2);
                   2824:        if ( t1 > 0 && t2 < 0 ) return 1;
                   2825:        else return 0;
                   2826: }
                   2827:
                   2828: /* 0 < u => 1, 0 > u => -1 */
                   2829:
                   2830: int compare_zero(int n,int *u,int row,int **w)
                   2831: {
                   2832:        int i,j,t;
                   2833:        int *wi;
                   2834:
                   2835:        for ( i = 0; i < row; i++ ) {
                   2836:                wi = w[i];
                   2837:                for ( j = 0, t = 0; j < n; j++ ) t += u[j]*wi[j];
                   2838:                if ( t > 0 ) return 1;
                   2839:                else if ( t < 0 ) return -1;
                   2840:        }
                   2841:        return 0;
                   2842: }
                   2843:
                   2844: /* functions for generic groebner walk */
                   2845: /* u=0 means u=-infty */
                   2846:
                   2847: int compare_facet_preorder(int n,int *u,int *v,
                   2848:        int row1,int **w1,int row2,int **w2)
                   2849: {
                   2850:        int i,j,s,t,tu,tv;
                   2851:        int *w2i,*uv;
                   2852:
                   2853:        if ( !u ) return 1;
                   2854:        uv = W_ALLOC(n);
                   2855:        for ( i = 0; i < row2; i++ ) {
                   2856:                w2i = w2[i];
                   2857:                for ( j = 0, tu = tv = 0; j < n; j++ )
                   2858:                        if ( s = w2i[j] ) {
                   2859:                                tu += s*u[j]; tv += s*v[j];
                   2860:                        }
                   2861:                for ( j = 0; j < n; j++ ) uv[j] = u[j]*tv-v[j]*tu;
                   2862:                t = compare_zero(n,uv,row1,w1);
                   2863:                if ( t > 0 ) return 1;
                   2864:                else if ( t < 0 ) return 0;
                   2865:        }
                   2866:        return 1;
                   2867: }
                   2868:
1.48      noro     2869: Q inner_product_with_small_vector(VECT w,int *v)
                   2870: {
                   2871:        int n,i;
                   2872:        Q q,s,t,u;
                   2873:
                   2874:        n = w->len;
                   2875:        s = 0;
                   2876:        for ( i = 0; i < n; i++ ) {
                   2877:                STOQ(v[i],q); mulq((Q)w->body[i],q,&t); addq(t,s,&u); s = u;
                   2878:        }
                   2879:        return s;
                   2880: }
                   2881:
                   2882: Q compute_last_t(NODE g,NODE gh,Q t,VECT w1,VECT w2,NODE *homo,VECT *wp)
                   2883: {
                   2884:        int n,i;
                   2885:        int *wt;
                   2886:        Q last,d1,d2,dn,nm,s,t1;
                   2887:        VECT wd,wt1,wt2,w;
                   2888:        NODE tg,tgh;
                   2889:        MP f;
                   2890:        int *h;
                   2891:        NODE r0,r;
                   2892:        MP m0,m;
                   2893:        DP d;
                   2894:
                   2895:        n = w1->len;
                   2896:        wt = W_ALLOC(n);
                   2897:        last = ONE;
                   2898:        /* t1 = 1-t */
                   2899:        for ( tg = g, tgh = gh; tg; tg = NEXT(tg), tgh = NEXT(tgh ) ) {
                   2900:                f = BDY((DP)BDY(tg));
                   2901:                h = BDY((DP)BDY(tgh))->dl->d;
                   2902:                for ( ; f; f = NEXT(f) ) {
                   2903:                        for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
                   2904:                        for ( i = 0; i < n && !wt[i]; i++ );
                   2905:                        if ( i == n ) continue;
                   2906:                        d1 = inner_product_with_small_vector(w1,wt);
                   2907:                        d2 = inner_product_with_small_vector(w2,wt);
                   2908:                        nm = d1; subq(d1,d2,&dn);
                   2909:                        /* if d1=d2 then nothing happens */
                   2910:                        if ( !dn ) continue;
                   2911:                        /* s satisfies ds = 0*/
                   2912:                        divq(nm,dn,&s);
                   2913:
                   2914:                        if ( cmpq(s,t) > 0 && cmpq(s,last) < 0 )
                   2915:                                last = s;
                   2916:                        else if ( !cmpq(s,t) ) {
                   2917:                                if ( cmpq(d2,0) < 0 ) {
                   2918:                                        last = t;
                   2919:                                        break;
                   2920:                                }
                   2921:                        }
                   2922:                }
                   2923:        }
                   2924:        if ( !last ) {
                   2925:                dn = ONE; nm = 0;
                   2926:        } else {
                   2927:                NTOQ(NM(last),1,nm);
                   2928:                if ( INT(last) ) dn = ONE;
                   2929:                else {
                   2930:                        NTOQ(DN(last),1,dn);
                   2931:                }
                   2932:        }
                   2933:        /* (1-n/d)*w1+n/d*w2 -> w=(d-n)*w1+n*w2 */
                   2934:        subq(dn,nm,&t1); mulvect(CO,(Obj)w1,(Obj)t1,(Obj *)&wt1);
                   2935:        mulvect(CO,(Obj)w2,(Obj)nm,(Obj *)&wt2); addvect(CO,wt1,wt2,&w);
                   2936:
                   2937:        r0 = 0;
                   2938:        for ( tg = g, tgh = gh; tg; tg = NEXT(tg), tgh = NEXT(tgh ) ) {
                   2939:                f = BDY((DP)BDY(tg));
                   2940:                h = BDY((DP)BDY(tgh))->dl->d;
                   2941:                for ( m0 = 0; f; f = NEXT(f) ) {
                   2942:                        for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
                   2943:                        for ( i = 0; i < n && !wt[i]; i++ );
                   2944:                        if ( !inner_product_with_small_vector(w,wt) ) {
                   2945:                                NEXTMP(m0,m); m->c = f->c; m->dl = f->dl;
                   2946:                        }
                   2947:                }
                   2948:                NEXT(m) = 0;
                   2949:                MKDP(((DP)BDY(tg))->nv,m0,d);  d->sugar = ((DP)BDY(tg))->sugar;
                   2950:                NEXTNODE(r0,r); BDY(r) = (pointer)d;
                   2951:        }
                   2952:        NEXT(r) = 0;
                   2953:        *homo = r0;
                   2954:        *wp = w;
                   2955:        return last;
                   2956: }
                   2957:
1.42      noro     2958: /* return 0 if last_w = infty */
                   2959:
                   2960: NODE compute_last_w(NODE g,NODE gh,int n,int **w,
                   2961:        int row1,int **w1,int row2,int **w2)
                   2962: {
                   2963:        DP d;
                   2964:        MP f,m0,m;
                   2965:        int *wt,*v,*h;
                   2966:        NODE t,s,n0,tn,n1,r0,r;
                   2967:        int i;
                   2968:
                   2969:        wt = W_ALLOC(n);
                   2970:        n0 = 0;
                   2971:        for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) {
                   2972:                f = BDY((DP)BDY(t));
                   2973:                h = BDY((DP)BDY(s))->dl->d;
                   2974:                for ( ; f; f = NEXT(f) ) {
                   2975:                        for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
                   2976:                        for ( i = 0; i < n && !wt[i]; i++ );
                   2977:                        if ( i == n ) continue;
                   2978:
                   2979:                        if ( in_c12(n,wt,row1,w1,row2,w2) &&
                   2980:                                compare_facet_preorder(n,*w,wt,row1,w1,row2,w2) ) {
                   2981:                                v = (int *)MALLOC_ATOMIC(n*sizeof(int));
                   2982:                                for ( i = 0; i < n; i++ ) v[i] = wt[i];
                   2983:                                MKNODE(n1,v,n0); n0 = n1;
                   2984:                        }
                   2985:                }
                   2986:        }
                   2987:        if ( !n0 ) return 0;
                   2988:        for ( t = n0; t; t = NEXT(t) ) {
                   2989:                v = (int *)BDY(t);
                   2990:                for ( s = n0; s; s = NEXT(s) )
                   2991:                        if ( !compare_facet_preorder(n,v,(int *)BDY(s),row1,w1,row2,w2) )
                   2992:                                break;
                   2993:                if ( !s ) {
                   2994:                        *w = v;
                   2995:                        break;
                   2996:                }
                   2997:        }
                   2998:        if ( !t )
                   2999:                error("compute_last_w : cannot happen");
                   3000:        r0 = 0;
                   3001:        for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) {
                   3002:                f = BDY((DP)BDY(t));
                   3003:                h = BDY((DP)BDY(s))->dl->d;
                   3004:                for ( m0 = 0; f; f = NEXT(f) ) {
                   3005:                        for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
                   3006:                        for ( i = 0; i < n && !wt[i]; i++ );
                   3007:                        if ( i == n  ||
                   3008:                                (compare_facet_preorder(n,wt,*w,row1,w1,row2,w2)
                   3009:                                && compare_facet_preorder(n,*w,wt,row1,w1,row2,w2)) ) {
                   3010:                                NEXTMP(m0,m); m->c = f->c; m->dl = f->dl;
                   3011:                        }
                   3012:                }
1.43      noro     3013:                NEXT(m) = 0;
1.42      noro     3014:                MKDP(((DP)BDY(t))->nv,m0,d);  d->sugar = ((DP)BDY(t))->sugar;
                   3015:                NEXTNODE(r0,r); BDY(r) = (pointer)d;
                   3016:        }
                   3017:        NEXT(r) = 0;
                   3018:        return r0;
                   3019: }
1.44      noro     3020:
                   3021: /* compute a sufficient set of d(f)=u-v */
                   3022:
                   3023: NODE compute_essential_df(DP *g,DP *gh,int ng)
                   3024: {
1.45      noro     3025:        int nv,i,j,k,t,lj;
                   3026:        NODE r,r1,ri,rt,r0;
                   3027:        MP m;
                   3028:        MP *mj;
                   3029:        DL di,hj,dl,dlt;
                   3030:        int *d,*dt;
                   3031:        LIST l;
1.44      noro     3032:        Q q;
1.45      noro     3033:
                   3034:        nv = g[0]->nv;
                   3035:        r = 0;
                   3036:        for ( j = 0; j < ng; j++ ) {
                   3037:                for ( m = BDY(g[j]), lj = 0; m; m = NEXT(m), lj++ );
                   3038:                mj = (MP *)ALLOCA(lj*sizeof(MP));
                   3039:                for ( m = BDY(g[j]), k = 0; m; m = NEXT(m), k++ )
                   3040:                        mj[k] = m;
                   3041:                for ( i = 0; i < lj; i++ ) {
                   3042:                        for ( di = mj[i]->dl, k = i+1; k < lj; k++ )
                   3043:                                if ( _dl_redble(di,mj[k]->dl,nv) ) break;
                   3044:                        if ( k < lj ) mj[i] = 0;
                   3045:                }
                   3046:                hj = BDY(gh[j])->dl;
                   3047:                _NEWDL(dl,nv); d = dl->d;
                   3048:                r0 = r;
                   3049:                for ( i = 0; i < lj; i++ ) {
                   3050:                        if ( mj[i] && !dl_equal(nv,di=mj[i]->dl,hj) ) {
                   3051:                                for ( k = 0, t = 0; k < nv; k++ ) {
                   3052:                                        d[k] = hj->d[k]-di->d[k];
                   3053:                                        t += d[k];
                   3054:                                }
                   3055:                                dl->td = t;
                   3056: #if 1
                   3057:                                for ( rt = r0; rt; rt = NEXT(rt) ) {
                   3058:                                        dlt = (DL)BDY(rt);
                   3059:                                        if ( dlt->td != dl->td ) continue;
                   3060:                                        for ( dt = dlt->d, k = 0; k < nv; k++ )
                   3061:                                                if ( d[k] != dt[k] ) break;
                   3062:                                        if ( k == nv ) break;
                   3063:                                }
                   3064: #else
                   3065:                                rt = 0;
                   3066: #endif
                   3067:                                if ( !rt ) {
                   3068:                                        MKNODE(r1,dl,r); r = r1;
                   3069:                                        _NEWDL(dl,nv); d = dl->d;
                   3070:                                }
1.44      noro     3071:                        }
                   3072:                }
                   3073:        }
1.45      noro     3074:        for ( rt = r; rt; rt = NEXT(rt) ) {
                   3075:                dl = (DL)BDY(rt); d = dl->d;
                   3076:                ri = 0;
                   3077:                for ( k = nv-1; k >= 0; k-- ) {
                   3078:                        STOQ(d[k],q);
                   3079:                        MKNODE(r1,q,ri); ri = r1;
1.44      noro     3080:                }
1.45      noro     3081:                MKNODE(r1,0,ri); MKLIST(l,r1);
                   3082:                BDY(rt) = (pointer)l;
1.44      noro     3083:        }
                   3084:        return r;
                   3085: }
1.57      noro     3086:
                   3087: int comp_bits_divisible(int *a,int *b,int n)
                   3088: {
                   3089:        int bpi,i,wi,bi;
                   3090:
                   3091:        bpi = (sizeof(int)/sizeof(char))*8;
                   3092:        for ( i = 0; i < n; i++ ) {
                   3093:                wi = i/bpi; bi = i%bpi;
                   3094:                if ( !(a[wi]&(1<<bi)) && (b[wi]&(1<<bi)) ) return 0;
                   3095:        }
                   3096:        return 1;
                   3097: }
                   3098:
                   3099: int comp_bits_lex(int *a,int *b,int n)
                   3100: {
                   3101:        int bpi,i,wi,ba,bb,bi;
                   3102:
                   3103:        bpi = (sizeof(int)/sizeof(char))*8;
                   3104:        for ( i = 0; i < n; i++ ) {
                   3105:                wi = i/bpi; bi = i%bpi;
                   3106:                ba = (a[wi]&(1<<bi))?1:0;
                   3107:                bb = (b[wi]&(1<<bi))?1:0;
                   3108:                if ( ba > bb ) return 1;
                   3109:                else if ( ba < bb ) return -1;
                   3110:        }
                   3111:        return 0;
                   3112: }
                   3113:
                   3114: NODE mono_raddec(NODE ideal)
                   3115: {
                   3116:        DP p;
                   3117:        int nv,w,i,bpi,di,c,len;
                   3118:        int *d,*s,*u,*new;
                   3119:        NODE t,t1,v,r,rem,prev;
                   3120:
                   3121:        if( !ideal ) return 0;
                   3122:        p = (DP)BDY(ideal);
                   3123:        nv = NV(p);
                   3124:        bpi = (sizeof(int)/sizeof(char))*8;
                   3125:        w = (nv+(bpi-1))/bpi;
                   3126:        d = p->body->dl->d;
                   3127:        if ( !NEXT(ideal) )     {
                   3128:                for ( t = 0, i = nv-1; i >= 0; i-- ) {
                   3129:                        if ( d[i] ) {
                   3130:                                s = (int *)CALLOC(w,sizeof(int));
                   3131:                                s[i/bpi] |= 1<<(i%bpi);
                   3132:                                MKNODE(t1,s,t);
                   3133:                                t = t1;
                   3134:                        }
                   3135:                }
                   3136:                return t;
                   3137:        }
                   3138:        rem = mono_raddec(NEXT(ideal));
                   3139:        r = 0;
                   3140:        len = w*sizeof(int);
                   3141:        u = (int *)CALLOC(w,sizeof(int));
                   3142:        for ( i = nv-1; i >= 0; i-- ) {
                   3143:                if ( d[i] ) {
                   3144:                        for ( t = rem; t; t = NEXT(t) ) {
                   3145:                                bcopy((char *)BDY(t),(char *)u,len);
                   3146:                                u[i/bpi] |= 1<<(i%bpi);
                   3147:                                for ( v = r; v; v = NEXT(v) ) {
                   3148:                                        if ( comp_bits_divisible(u,(int *)BDY(v),nv) ) break;
                   3149:                                }
                   3150:                                if ( v ) continue;
                   3151:                                for ( v = r, prev = 0; v; v = NEXT(v) ) {
                   3152:                                        if ( comp_bits_divisible((int *)BDY(v),u,nv) ) {
                   3153:                                                if ( prev ) NEXT(prev) = NEXT(v);
                   3154:                                                else r = NEXT(r);
                   3155:                                        } else prev =v;
                   3156:                                }
                   3157:                                for ( v = r, prev = 0; v; prev = v, v = NEXT(v) ) {
                   3158:                                        if ( comp_bits_lex(u,(int *)BDY(v),nv) < 0 ) break;
                   3159:                                }
                   3160:                                new = (int *)CALLOC(w,sizeof(int));
                   3161:                                bcopy((char *)u,(char *)new,len);
                   3162:                                MKNODE(t1,new,v);
                   3163:                                if ( prev ) NEXT(prev) = t1;
                   3164:                                else r = t1;
                   3165:                        }
                   3166:                }
                   3167:        }
                   3168:        return r;
                   3169: }

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