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

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.10    ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.9 2000/12/08 08:26:08 noro Exp $
1.2       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include "base.h"
                     52: #include "parse.h"
                     53: #include "ox.h"
                     54:
1.5       noro       55: #define HMAG(p) (p_mag(BDY(p)->c))
                     56:
1.1       noro       57: extern int (*cmpdl)();
1.5       noro       58: extern double pz_t_e,pz_t_d,pz_t_d1,pz_t_c;
                     59: extern int dp_nelim,dp_fcoeffs;
1.7       noro       60: extern int NoGCD;
                     61: extern int GenTrace;
                     62: extern NODE TraceList;
                     63:
                     64: /*
                     65:  * content reduction
                     66:  *
                     67:  */
                     68:
                     69: void dp_ptozp(p,rp)
                     70: DP p,*rp;
                     71: {
                     72:        MP m,mr,mr0;
                     73:        int i,n;
                     74:        Q *w;
                     75:        Q dvr;
                     76:        P t;
                     77:
                     78:        if ( !p )
                     79:                *rp = 0;
                     80:        else {
                     81:                for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ );
                     82:                w = (Q *)ALLOCA(n*sizeof(Q));
                     83:                for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
                     84:                        if ( NUM(m->c) )
                     85:                                w[i] = (Q)m->c;
                     86:                        else
                     87:                                ptozp(m->c,1,&w[i],&t);
                     88:                sortbynm(w,n);
                     89:                qltozl(w,n,&dvr);
                     90:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                     91:                        NEXTMP(mr0,mr); divsp(CO,m->c,(P)dvr,&mr->c); mr->dl = m->dl;
                     92:                }
                     93:                NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                     94:        }
                     95: }
                     96:
                     97: void dp_ptozp2(p0,p1,hp,rp)
                     98: DP p0,p1;
                     99: DP *hp,*rp;
                    100: {
                    101:        DP t,s,h,r;
                    102:        MP m,mr,mr0,m0;
                    103:
                    104:        addd(CO,p0,p1,&t); dp_ptozp(t,&s);
                    105:        if ( !p0 ) {
                    106:                h = 0; r = s;
                    107:        } else if ( !p1 ) {
                    108:                h = s; r = 0;
                    109:        } else {
                    110:                for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
                    111:                        m = NEXT(m), m0 = NEXT(m0) ) {
                    112:                        NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
                    113:                }
                    114:                NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
                    115:        }
                    116:        if ( h )
                    117:                h->sugar = p0->sugar;
                    118:        if ( r )
                    119:                r->sugar = p1->sugar;
                    120:        *hp = h; *rp = r;
                    121: }
1.1       noro      122:
                    123: void dp_idiv(p,c,rp)
                    124: DP p;
                    125: Q c;
                    126: DP *rp;
                    127: {
                    128:        Q t;
                    129:        N nm,q;
                    130:        int sgn,s;
                    131:        MP mr0,m,mr;
                    132:
                    133:        if ( !p )
                    134:                *rp = 0;
                    135:        else if ( MUNIQ((Q)c) )
                    136:                *rp = p;
                    137:        else if ( MUNIQ((Q)c) )
                    138:                chsgnd(p,rp);
                    139:        else {
                    140:                nm = NM(c); sgn = SGN(c);
                    141:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    142:                        NEXTMP(mr0,mr);
                    143:
                    144:                        divsn(NM((Q)(m->c)),nm,&q);
                    145:                        s = sgn*SGN((Q)(m->c));
                    146:                        NTOQ(q,s,t);
                    147:                        mr->c = (P)t;
                    148:                        mr->dl = m->dl;
                    149:                }
                    150:                NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
                    151:                if ( *rp )
                    152:                        (*rp)->sugar = p->sugar;
                    153:        }
                    154: }
                    155:
                    156: void dp_mbase(hlist,mbase)
                    157: NODE hlist;
                    158: NODE *mbase;
                    159: {
                    160:        DL *dl;
                    161:        DL d;
                    162:        int i,j,n,nvar,td;
                    163:
                    164:        n = length(hlist); nvar = ((DP)BDY(hlist))->nv;
                    165:        dl = (DL *)MALLOC(n*sizeof(DL));
                    166:        for ( i = 0; i < n; i++, hlist = NEXT(hlist) )
                    167:                dl[i] = BDY((DP)BDY(hlist))->dl;
                    168:        NEWDL(d,nvar); *mbase = 0;
                    169:        while ( 1 ) {
                    170:                insert_to_node(d,mbase,nvar);
                    171:                for ( i = nvar-1; i >= 0; ) {
                    172:                        d->d[i]++; d->td++;
                    173:                        for ( j = 0; j < n; j++ ) {
                    174:                                if ( _dl_redble(dl[j],d,nvar) )
                    175:                                        break;
                    176:                        }
                    177:                        if ( j < n ) {
                    178:                                for ( j = nvar-1; j >= i; j-- )
                    179:                                        d->d[j] = 0;
                    180:                                for ( j = 0, td = 0; j < i; j++ )
                    181:                                        td += d->d[j];
                    182:                                d->td = td;
                    183:                                i--;
                    184:                        } else
                    185:                                break;
                    186:                }
                    187:                if ( i < 0 )
                    188:                        break;
                    189:        }
                    190: }
                    191:
                    192: int _dl_redble(d1,d2,nvar)
                    193: DL d1,d2;
                    194: int nvar;
                    195: {
                    196:        int i;
                    197:
                    198:        if ( d1->td > d2->td )
                    199:                return 0;
                    200:        for ( i = 0; i < nvar; i++ )
                    201:                if ( d1->d[i] > d2->d[i] )
                    202:                        break;
                    203:        if ( i < nvar )
                    204:                return 0;
                    205:        else
                    206:                return 1;
                    207: }
                    208:
                    209: void insert_to_node(d,n,nvar)
                    210: DL d;
                    211: NODE *n;
                    212: int nvar;
                    213: {
                    214:        DL d1;
                    215:        MP m;
                    216:        DP dp;
                    217:        NODE n0,n1,n2;
                    218:
                    219:        NEWDL(d1,nvar); d1->td = d->td;
                    220:        bcopy((char *)d->d,(char *)d1->d,nvar*sizeof(int));
                    221:        NEWMP(m); m->dl = d1; m->c = (P)ONE; NEXT(m) = 0;
                    222:        MKDP(nvar,m,dp); dp->sugar = d->td;
                    223:        if ( !(*n) ) {
                    224:                MKNODE(n1,dp,0); *n = n1;
                    225:        } else {
                    226:                for ( n1 = *n, n0 = 0; n1; n0 = n1, n1 = NEXT(n1) )
                    227:                        if ( (*cmpdl)(nvar,d,BDY((DP)BDY(n1))->dl) > 0 ) {
                    228:                                MKNODE(n2,dp,n1);
                    229:                                if ( !n0 )
                    230:                                        *n = n2;
                    231:                                else
                    232:                                        NEXT(n0) = n2;
                    233:                                break;
                    234:                        }
                    235:                if ( !n1 ) {
                    236:                        MKNODE(n2,dp,0); NEXT(n0) = n2;
                    237:                }
                    238:        }
                    239: }
                    240:
                    241: void dp_vtod(c,p,rp)
                    242: Q *c;
                    243: DP p;
                    244: DP *rp;
                    245: {
                    246:        MP mr0,m,mr;
                    247:        int i;
                    248:
                    249:        if ( !p )
                    250:                *rp = 0;
                    251:        else {
                    252:                for ( mr0 = 0, m = BDY(p), i = 0; m; m = NEXT(m), i++ ) {
                    253:                        NEXTMP(mr0,mr); mr->c = (P)c[i]; mr->dl = m->dl;
                    254:                }
                    255:                NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
                    256:                (*rp)->sugar = p->sugar;
                    257:        }
                    258: }
                    259:
1.8       noro      260: extern int mpi_mag;
                    261: extern int PCoeffs;
                    262:
                    263: void dp_ptozp_d(p,rp)
1.1       noro      264: DP p,*rp;
                    265: {
                    266:        int i,j,k,l,n,nsep;
                    267:        MP m;
                    268:        NODE tn,n0,n1,n2,n3;
                    269:        struct oVECT v;
                    270:        VECT c,cs;
                    271:        VECT qi,ri;
                    272:        LIST *qr;
                    273:        int s,id;
                    274:        Obj dmy;
                    275:        Q d0,d1,gcd,a,u,u1;
                    276:        Q *q,*r;
                    277:        STRING iqr_v;
                    278:        pointer *b;
                    279:        N qn,gn;
                    280:        double get_rtime();
                    281:        int blen;
1.8       noro      282:        NODE dist;
                    283:        int ndist;
1.1       noro      284:        double t0;
                    285:        double t_e,t_d,t_d1,t_c;
1.8       noro      286:        extern int DP_NFStat;
                    287:        extern LIST Dist;
1.1       noro      288:
                    289:        if ( !p )
                    290:                *rp = 0;
                    291:        else {
1.8       noro      292:                if ( PCoeffs ) {
                    293:                        dp_ptozp(p,rp); return;
                    294:                }
1.9       noro      295:                if ( !Dist || p_mag(BDY(p)->c) <= mpi_mag ) {
1.8       noro      296:                        dist = 0; ndist = 0;
                    297:                        if ( DP_NFStat ) fprintf(asir_out,"L");
                    298:                } else {
                    299:                        dist = BDY(Dist); ndist = length(dist);
                    300:                        if ( DP_NFStat ) fprintf(asir_out,"D");
                    301:                }
1.1       noro      302:                for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
                    303:                nsep = ndist + 1;
                    304:                if ( n <= nsep ) {
                    305:                        dp_ptozp(p,rp); return;
                    306:                }
                    307:                t0 = get_rtime();
                    308:                dp_dtov(p,&c);
                    309:                igcdv_estimate(c,&d0);
                    310:                t_e = get_rtime()-t0;
                    311:                t0 = get_rtime();
                    312:                dp_dtov(p,&c);
                    313:                sepvect(c,nsep,&cs);
                    314:                MKSTR(iqr_v,"iqr");
                    315:                qr = (LIST *)CALLOC(nsep,sizeof(LIST));
                    316:                q = (Q *)CALLOC(n,sizeof(Q));
                    317:                r = (Q *)CALLOC(n,sizeof(Q));
                    318:                for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
                    319:                        MKNODE(n3,d0,0); MKNODE(n2,b[i],n3);
                    320:                        MKNODE(n1,iqr_v,n2); MKNODE(n0,BDY(tn),n1);
                    321:                        Pox_rpc(n0,&dmy);
                    322:                }
                    323:                iqrv(b[i],d0,&qr[i]);
                    324:                dp_dtov(p,&c);
                    325:                for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
                    326:                        Pox_pop_local(tn,&qr[i]);
                    327:                        if ( OID(qr[i]) == O_ERR ) {
                    328:                                printexpr(CO,(Obj)qr[i]);
                    329:                                error("dp_ptozp_d : aborted");
                    330:                        }
                    331:                }
                    332:                t_d = get_rtime()-t0;
                    333:                t_d1 = t_d/n;
                    334:                t0 = get_rtime();
                    335:                for ( i = j = 0; i < nsep; i++ ) {
                    336:                        tn = BDY(qr[i]); qi = (VECT)BDY(tn); ri = (VECT)BDY(NEXT(tn));
                    337:                        for ( k = 0, l = qi->len; k < l; k++, j++ ) {
                    338:                                q[j] = (Q)BDY(qi)[k]; r[j] = (Q)BDY(ri)[k];
                    339:                        }
                    340:                }
                    341:                v.id = O_VECT; v.len = n; v.body = (pointer *)r; igcdv(&v,&d1);
                    342:                if ( d1 ) {
                    343:                        gcdn(NM(d0),NM(d1),&gn); NTOQ(gn,1,gcd);
                    344:                        divsn(NM(d0),gn,&qn); NTOQ(qn,1,a);
                    345:                        for ( i = 0; i < n; i++ ) {
                    346:                                mulq(a,q[i],&u);
                    347:                                if ( r[i] ) {
                    348:                                        divsn(NM(r[i]),gn,&qn); NTOQ(qn,SGN(r[i]),u1);
                    349:                                        addq(u,u1,&q[i]);
                    350:                                } else
                    351:                                        q[i] = u;
                    352:                        }
                    353:                } else
                    354:                        gcd = d0;
                    355:                dp_vtod(q,p,rp);
                    356:                t_c = get_rtime()-t0;
                    357:                blen=p_mag((P)gcd);
                    358:                pz_t_e += t_e; pz_t_d += t_d; pz_t_d1 += t_d1; pz_t_c += t_c;
                    359:                if ( 0 )
                    360:                        fprintf(stderr,"(%d,%d)",p_mag((P)d0)-blen,blen);
                    361:        }
                    362: }
                    363:
1.8       noro      364: void dp_ptozp2_d(p0,p1,hp,rp)
1.1       noro      365: DP p0,p1;
                    366: DP *hp,*rp;
                    367: {
                    368:        DP t,s,h,r;
                    369:        MP m,mr,mr0,m0;
                    370:
1.8       noro      371:        addd(CO,p0,p1,&t); dp_ptozp_d(t,&s);
1.1       noro      372:        if ( !p0 ) {
                    373:                h = 0; r = s;
                    374:        } else if ( !p1 ) {
                    375:                h = s; r = 0;
                    376:        } else {
                    377:                for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
                    378:                        m = NEXT(m), m0 = NEXT(m0) ) {
                    379:                        NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
                    380:                }
                    381:                NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
                    382:        }
                    383:        if ( h )
                    384:                h->sugar = p0->sugar;
                    385:        if ( r )
                    386:                r->sugar = p1->sugar;
                    387:        *hp = h; *rp = r;
1.5       noro      388: }
                    389:
1.7       noro      390: void dp_prim(p,rp)
                    391: DP p,*rp;
1.5       noro      392: {
1.7       noro      393:        P t,g;
                    394:        DP p1;
                    395:        MP m,mr,mr0;
                    396:        int i,n;
                    397:        P *w;
                    398:        Q *c;
                    399:        Q dvr;
1.5       noro      400:
1.7       noro      401:        if ( !p )
                    402:                *rp = 0;
                    403:        else if ( dp_fcoeffs )
                    404:                *rp = p;
                    405:        else if ( NoGCD )
                    406:                dp_ptozp(p,rp);
                    407:        else {
                    408:                dp_ptozp(p,&p1); p = p1;
                    409:                for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
                    410:                if ( n == 1 ) {
                    411:                        m = BDY(p);
                    412:                        NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
                    413:                        MKDP(p->nv,mr,*rp); (*rp)->sugar = p->sugar;
                    414:                        return;
                    415:                }
                    416:                w = (P *)ALLOCA(n*sizeof(P));
                    417:                c = (Q *)ALLOCA(n*sizeof(Q));
                    418:                for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
                    419:                        if ( NUM(m->c) ) {
                    420:                                c[i] = (Q)m->c; w[i] = (P)ONE;
                    421:                        } else
                    422:                                ptozp(m->c,1,&c[i],&w[i]);
                    423:                qltozl(c,n,&dvr); heu_nezgcdnpz(CO,w,n,&t); mulp(CO,t,(P)dvr,&g);
                    424:                if ( NUM(g) )
                    425:                        *rp = p;
                    426:                else {
                    427:                        for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    428:                                NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl;
                    429:                        }
                    430:                        NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5       noro      431:                }
1.7       noro      432:        }
1.5       noro      433: }
                    434:
                    435: void heu_nezgcdnpz(vl,pl,m,pr)
                    436: VL vl;
                    437: P *pl,*pr;
                    438: int m;
                    439: {
                    440:        int i,r;
                    441:        P gcd,t,s1,s2,u;
                    442:        Q rq;
                    443:
                    444:        while ( 1 ) {
                    445:                for ( i = 0, s1 = 0; i < m; i++ ) {
                    446:                        r = random(); UTOQ(r,rq);
                    447:                        mulp(vl,pl[i],(P)rq,&t); addp(vl,s1,t,&u); s1 = u;
                    448:                }
                    449:                for ( i = 0, s2 = 0; i < m; i++ ) {
                    450:                        r = random(); UTOQ(r,rq);
                    451:                        mulp(vl,pl[i],(P)rq,&t); addp(vl,s2,t,&u); s2 = u;
                    452:                }
                    453:                ezgcdp(vl,s1,s2,&gcd);
                    454:                for ( i = 0; i < m; i++ ) {
                    455:                        if ( !divtpz(vl,pl[i],gcd,&t) )
                    456:                                break;
                    457:                }
                    458:                if ( i == m )
                    459:                        break;
                    460:        }
                    461:        *pr = gcd;
                    462: }
                    463:
                    464: void dp_prim_mod(p,mod,rp)
                    465: int mod;
                    466: DP p,*rp;
                    467: {
                    468:        P t,g;
                    469:        MP m,mr,mr0;
                    470:
                    471:        if ( !p )
                    472:                *rp = 0;
                    473:        else if ( NoGCD )
                    474:                *rp = p;
                    475:        else {
                    476:                for ( m = BDY(p), g = m->c, m = NEXT(m); m; m = NEXT(m) ) {
                    477:                        gcdprsmp(CO,mod,g,m->c,&t); g = t;
                    478:                }
                    479:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    480:                        NEXTMP(mr0,mr); divsmp(CO,mod,m->c,g,&mr->c); mr->dl = m->dl;
                    481:                }
                    482:                NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                    483:        }
                    484: }
                    485:
1.7       noro      486: void dp_cont(p,rp)
1.5       noro      487: DP p;
1.7       noro      488: Q *rp;
1.5       noro      489: {
1.7       noro      490:        VECT v;
1.5       noro      491:
1.7       noro      492:        dp_dtov(p,&v); igcdv(v,rp);
1.5       noro      493: }
                    494:
1.7       noro      495: void dp_dtov(dp,rp)
                    496: DP dp;
                    497: VECT *rp;
1.5       noro      498: {
1.7       noro      499:        MP m,t;
                    500:        int i,n;
                    501:        VECT v;
                    502:        pointer *p;
1.5       noro      503:
1.7       noro      504:        m = BDY(dp);
                    505:        for ( t = m, n = 0; t; t = NEXT(t), n++ );
                    506:        MKVECT(v,n);
                    507:        for ( i = 0, p = BDY(v), t = m; i < n; t = NEXT(t), i++ )
                    508:                p[i] = (pointer)(t->c);
                    509:        *rp = v;
1.5       noro      510: }
                    511:
1.7       noro      512: /*
                    513:  * s-poly computation
                    514:  *
                    515:  */
1.5       noro      516:
1.7       noro      517: void dp_sp(p1,p2,rp)
                    518: DP p1,p2;
1.5       noro      519: DP *rp;
                    520: {
1.7       noro      521:        int i,n,td;
                    522:        int *w;
                    523:        DL d1,d2,d;
                    524:        MP m;
                    525:        DP t,s1,s2,u;
                    526:        Q c,c1,c2;
                    527:        N gn,tn;
1.5       noro      528:
1.7       noro      529:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    530:        w = (int *)ALLOCA(n*sizeof(int));
                    531:        for ( i = 0, td = 0; i < n; i++ ) {
                    532:                w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
1.5       noro      533:        }
1.7       noro      534:
                    535:        NEWDL(d,n); d->td = td - d1->td;
                    536:        for ( i = 0; i < n; i++ )
                    537:                d->d[i] = w[i] - d1->d[i];
                    538:        c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
                    539:        if ( INT(c1) && INT(c2) ) {
                    540:                gcdn(NM(c1),NM(c2),&gn);
                    541:                if ( !UNIN(gn) ) {
                    542:                        divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
                    543:                        divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
1.5       noro      544:                }
                    545:        }
1.7       noro      546:
                    547:        NEWMP(m); m->dl = d; m->c = (P)c2; NEXT(m) = 0;
                    548:        MKDP(n,m,s1); s1->sugar = d->td; muld(CO,s1,p1,&t);
                    549:
                    550:        NEWDL(d,n); d->td = td - d2->td;
                    551:        for ( i = 0; i < n; i++ )
                    552:                d->d[i] = w[i] - d2->d[i];
                    553:        NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0;
                    554:        MKDP(n,m,s2); s2->sugar = d->td; muld(CO,s2,p2,&u);
                    555:
                    556:        subd(CO,t,u,rp);
                    557:        if ( GenTrace ) {
                    558:                LIST hist;
                    559:                NODE node;
                    560:
                    561:                node = mknode(4,ONE,0,s1,ONE);
                    562:                MKLIST(hist,node);
                    563:                MKNODE(TraceList,hist,0);
                    564:
                    565:                node = mknode(4,ONE,0,0,ONE);
                    566:                chsgnd(s2,(DP *)&ARG2(node));
                    567:                MKLIST(hist,node);
                    568:                MKNODE(node,hist,TraceList); TraceList = node;
                    569:        }
                    570: }
                    571:
                    572: void dp_sp_mod(p1,p2,mod,rp)
                    573: DP p1,p2;
                    574: int mod;
                    575: DP *rp;
                    576: {
                    577:        int i,n,td;
                    578:        int *w;
                    579:        DL d1,d2,d;
                    580:        MP m;
                    581:        DP t,s,u;
                    582:
                    583:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    584:        w = (int *)ALLOCA(n*sizeof(int));
                    585:        for ( i = 0, td = 0; i < n; i++ ) {
                    586:                w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
                    587:        }
                    588:        NEWDL(d,n); d->td = td - d1->td;
                    589:        for ( i = 0; i < n; i++ )
                    590:                d->d[i] = w[i] - d1->d[i];
                    591:        NEWMP(m); m->dl = d; m->c = (P)BDY(p2)->c; NEXT(m) = 0;
                    592:        MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p1,s,&t);
                    593:        NEWDL(d,n); d->td = td - d2->td;
                    594:        for ( i = 0; i < n; i++ )
                    595:                d->d[i] = w[i] - d2->d[i];
                    596:        NEWMP(m); m->dl = d; m->c = (P)BDY(p1)->c; NEXT(m) = 0;
                    597:        MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p2,s,&u);
                    598:        submd(CO,mod,t,u,rp);
                    599: }
                    600:
                    601: void _dp_sp_mod_dup(p1,p2,mod,rp)
                    602: DP p1,p2;
                    603: int mod;
                    604: DP *rp;
                    605: {
                    606:        int i,n,td;
                    607:        int *w;
                    608:        DL d1,d2,d;
                    609:        MP m;
                    610:        DP t,s,u;
                    611:
                    612:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    613:        w = (int *)ALLOCA(n*sizeof(int));
                    614:        for ( i = 0, td = 0; i < n; i++ ) {
                    615:                w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
                    616:        }
                    617:        _NEWDL(d,n); d->td = td - d1->td;
                    618:        for ( i = 0; i < n; i++ )
                    619:                d->d[i] = w[i] - d1->d[i];
                    620:        _NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
                    621:        _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p1,&t); _free_dp(s);
                    622:        _NEWDL(d,n); d->td = td - d2->td;
                    623:        for ( i = 0; i < n; i++ )
                    624:                d->d[i] = w[i] - d2->d[i];
                    625:        _NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
                    626:        _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p2,&u); _free_dp(s);
                    627:        _addmd_destructive(mod,t,u,rp);
                    628: }
                    629:
                    630: void _dp_sp_mod(p1,p2,mod,rp)
                    631: DP p1,p2;
                    632: int mod;
                    633: DP *rp;
                    634: {
                    635:        int i,n,td;
                    636:        int *w;
                    637:        DL d1,d2,d;
                    638:        MP m;
                    639:        DP t,s,u;
                    640:
                    641:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    642:        w = (int *)ALLOCA(n*sizeof(int));
                    643:        for ( i = 0, td = 0; i < n; i++ ) {
                    644:                w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
                    645:        }
                    646:        NEWDL(d,n); d->td = td - d1->td;
                    647:        for ( i = 0; i < n; i++ )
                    648:                d->d[i] = w[i] - d1->d[i];
                    649:        NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
                    650:        MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p1,&t);
                    651:        NEWDL(d,n); d->td = td - d2->td;
                    652:        for ( i = 0; i < n; i++ )
                    653:                d->d[i] = w[i] - d2->d[i];
                    654:        NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
                    655:        MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p2,&u);
                    656:        addmd_destructive(mod,t,u,rp);
                    657: }
                    658:
                    659: /*
                    660:  * m-reduction
                    661:  *
                    662:  */
                    663:
                    664: void dp_red(p0,p1,p2,head,rest,dnp,multp)
                    665: DP p0,p1,p2;
                    666: DP *head,*rest;
                    667: P *dnp;
                    668: DP *multp;
                    669: {
                    670:        int i,n;
                    671:        DL d1,d2,d;
                    672:        MP m;
                    673:        DP t,s,r,h;
                    674:        Q c,c1,c2;
                    675:        N gn,tn;
                    676:        P g,a;
                    677:
                    678:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    679:        NEWDL(d,n); d->td = d1->td - d2->td;
                    680:        for ( i = 0; i < n; i++ )
                    681:                d->d[i] = d1->d[i]-d2->d[i];
                    682:        c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
                    683:        if ( dp_fcoeffs ) {
                    684:                /* do nothing */
                    685:        } else 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:        } else {
                    692:                ezgcdpz(CO,(P)c1,(P)c2,&g);
                    693:                divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
                    694:        }
                    695:        NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
                    696:        *multp = s;
                    697:        muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); addd(CO,s,t,&r);
                    698:        muldc(CO,p0,(P)c2,&h);
                    699:        *head = h; *rest = r; *dnp = (P)c2;
                    700: }
                    701:
                    702: void dp_red_mod(p0,p1,p2,mod,head,rest,dnp)
                    703: DP p0,p1,p2;
                    704: int mod;
                    705: DP *head,*rest;
                    706: P *dnp;
                    707: {
                    708:        int i,n;
                    709:        DL d1,d2,d;
                    710:        MP m;
                    711:        DP t,s,r,h;
                    712:        P c1,c2,g,u;
                    713:
                    714:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    715:        NEWDL(d,n); d->td = d1->td - d2->td;
                    716:        for ( i = 0; i < n; i++ )
                    717:                d->d[i] = d1->d[i]-d2->d[i];
                    718:        c1 = (P)BDY(p1)->c; c2 = (P)BDY(p2)->c;
                    719:        gcdprsmp(CO,mod,c1,c2,&g);
                    720:        divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
                    721:        if ( NUM(c2) ) {
                    722:                divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
                    723:        }
                    724:        NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0;
                    725:        MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p2,s,&t);
                    726:        if ( NUM(c2) ) {
                    727:                addmd(CO,mod,p1,t,&r); h = p0;
                    728:        } else {
                    729:                mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
                    730:        }
                    731:        *head = h; *rest = r; *dnp = c2;
                    732: }
                    733:
1.10    ! noro      734: struct oEGT eg_red_mod;
        !           735:
1.7       noro      736: void _dp_red_mod_destructive(p1,p2,mod,rp)
                    737: DP p1,p2;
                    738: int mod;
                    739: DP *rp;
                    740: {
                    741:        int i,n;
                    742:        DL d1,d2,d;
                    743:        MP m;
                    744:        DP t,s;
                    745:        int c,c1;
1.10    ! noro      746: struct oEGT t0,t1;
1.7       noro      747:
                    748:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    749:        _NEWDL(d,n); d->td = d1->td - d2->td;
                    750:        for ( i = 0; i < n; i++ )
                    751:                d->d[i] = d1->d[i]-d2->d[i];
                    752:        c = invm(ITOS(BDY(p2)->c),mod); c1 = dmar(c,ITOS(BDY(p1)->c),0,mod);
                    753:        _NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
                    754:        _MKDP(n,m,s); s->sugar = d->td;
                    755:        _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.10    ! noro      756: /* get_eg(&t0); */
1.7       noro      757:        _addmd_destructive(mod,p1,t,rp);
1.10    ! noro      758: /* get_eg(&t1); add_eg(&eg_red_mod,&t0,&t1); */
1.7       noro      759: }
                    760:
                    761: /*
                    762:  * normal form computation
                    763:  *
                    764:  */
1.5       noro      765:
                    766: void dp_true_nf(b,g,ps,full,rp,dnp)
                    767: NODE b;
                    768: DP g;
                    769: DP *ps;
                    770: int full;
                    771: DP *rp;
                    772: P *dnp;
                    773: {
                    774:        DP u,p,d,s,t,dmy;
                    775:        NODE l;
                    776:        MP m,mr;
                    777:        int i,n;
                    778:        int *wb;
                    779:        int sugar,psugar;
                    780:        P dn,tdn,tdn1;
                    781:
                    782:        dn = (P)ONE;
                    783:        if ( !g ) {
                    784:                *rp = 0; *dnp = dn; return;
                    785:        }
                    786:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                    787:        wb = (int *)ALLOCA(n*sizeof(int));
                    788:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                    789:                wb[i] = QTOS((Q)BDY(l));
                    790:        sugar = g->sugar;
                    791:        for ( d = 0; g; ) {
                    792:                for ( u = 0, i = 0; i < n; i++ ) {
                    793:                        if ( dp_redble(g,p = ps[wb[i]]) ) {
                    794:                                dp_red(d,g,p,&t,&u,&tdn,&dmy);
                    795:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                    796:                                sugar = MAX(sugar,psugar);
                    797:                                if ( !u ) {
                    798:                                        if ( d )
                    799:                                                d->sugar = sugar;
                    800:                                        *rp = d; *dnp = dn; return;
                    801:                                } else {
                    802:                                        d = t;
                    803:                                        mulp(CO,dn,tdn,&tdn1); dn = tdn1;
                    804:                                }
                    805:                                break;
                    806:                        }
                    807:                }
                    808:                if ( u )
                    809:                        g = u;
                    810:                else if ( !full ) {
                    811:                        if ( g ) {
                    812:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                    813:                        }
                    814:                        *rp = g; *dnp = dn; return;
                    815:                } else {
                    816:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                    817:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                    818:                        addd(CO,d,t,&s); d = s;
                    819:                        dp_rest(g,&t); g = t;
                    820:                }
                    821:        }
                    822:        if ( d )
                    823:                d->sugar = sugar;
                    824:        *rp = d; *dnp = dn;
                    825: }
                    826:
                    827: void dp_nf_ptozp(b,g,ps,full,multiple,rp)
                    828: NODE b;
                    829: DP g;
                    830: DP *ps;
                    831: int full,multiple;
                    832: DP *rp;
                    833: {
                    834:        DP u,p,d,s,t,dmy1;
                    835:        P dmy;
                    836:        NODE l;
                    837:        MP m,mr;
                    838:        int i,n;
                    839:        int *wb;
                    840:        int hmag;
                    841:        int sugar,psugar;
                    842:
                    843:        if ( !g ) {
                    844:                *rp = 0; return;
                    845:        }
                    846:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                    847:        wb = (int *)ALLOCA(n*sizeof(int));
                    848:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                    849:                wb[i] = QTOS((Q)BDY(l));
                    850:        hmag = multiple*HMAG(g);
                    851:        sugar = g->sugar;
                    852:        for ( d = 0; g; ) {
                    853:                for ( u = 0, i = 0; i < n; i++ ) {
                    854:                        if ( dp_redble(g,p = ps[wb[i]]) ) {
                    855:                                dp_red(d,g,p,&t,&u,&dmy,&dmy1);
                    856:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                    857:                                sugar = MAX(sugar,psugar);
                    858:                                if ( !u ) {
                    859:                                        if ( d )
                    860:                                                d->sugar = sugar;
                    861:                                        *rp = d; return;
                    862:                                }
                    863:                                d = t;
                    864:                                break;
                    865:                        }
                    866:                }
                    867:                if ( u ) {
                    868:                        g = u;
                    869:                        if ( d ) {
1.7       noro      870:                                if ( multiple && HMAG(d) > hmag ) {
1.5       noro      871:                                        dp_ptozp2(d,g,&t,&u); d = t; g = u;
                    872:                                        hmag = multiple*HMAG(d);
                    873:                                }
                    874:                        } else {
1.7       noro      875:                                if ( multiple && HMAG(g) > hmag ) {
1.5       noro      876:                                        dp_ptozp(g,&t); g = t;
                    877:                                        hmag = multiple*HMAG(g);
                    878:                                }
                    879:                        }
                    880:                }
                    881:                else if ( !full ) {
                    882:                        if ( g ) {
                    883:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                    884:                        }
                    885:                        *rp = g; return;
                    886:                } else {
                    887:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                    888:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                    889:                        addd(CO,d,t,&s); d = s;
                    890:                        dp_rest(g,&t); g = t;
                    891:
                    892:                }
                    893:        }
                    894:        if ( d )
                    895:                d->sugar = sugar;
                    896:        *rp = d;
                    897: }
                    898:
                    899: void dp_nf_mod(b,g,ps,mod,full,rp)
                    900: NODE b;
                    901: DP g;
                    902: DP *ps;
                    903: int mod,full;
                    904: DP *rp;
                    905: {
                    906:        DP u,p,d,s,t;
                    907:        P dmy;
                    908:        NODE l;
                    909:        MP m,mr;
                    910:        int sugar,psugar;
                    911:
                    912:        if ( !g ) {
                    913:                *rp = 0; return;
                    914:        }
                    915:        sugar = g->sugar;
                    916:        for ( d = 0; g; ) {
                    917:                for ( u = 0, l = b; l; l = NEXT(l) ) {
                    918:                        if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
                    919:                                dp_red_mod(d,g,p,mod,&t,&u,&dmy);
                    920:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                    921:                                sugar = MAX(sugar,psugar);
                    922:                                if ( !u ) {
                    923:                                        if ( d )
                    924:                                                d->sugar = sugar;
                    925:                                        *rp = d; return;
                    926:                                }
                    927:                                d = t;
                    928:                                break;
                    929:                        }
                    930:                }
                    931:                if ( u )
                    932:                        g = u;
                    933:                else if ( !full ) {
                    934:                        if ( g ) {
                    935:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                    936:                        }
                    937:                        *rp = g; return;
                    938:                } else {
                    939:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                    940:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                    941:                        addmd(CO,mod,d,t,&s); d = s;
                    942:                        dp_rest(g,&t); g = t;
                    943:                }
                    944:        }
                    945:        if ( d )
                    946:                d->sugar = sugar;
                    947:        *rp = d;
                    948: }
                    949:
                    950: void dp_true_nf_mod(b,g,ps,mod,full,rp,dnp)
                    951: NODE b;
                    952: DP g;
                    953: DP *ps;
                    954: int mod,full;
                    955: DP *rp;
                    956: P *dnp;
                    957: {
                    958:        DP u,p,d,s,t;
                    959:        NODE l;
                    960:        MP m,mr;
                    961:        int i,n;
                    962:        int *wb;
                    963:        int sugar,psugar;
                    964:        P dn,tdn,tdn1;
                    965:
                    966:        dn = (P)ONEM;
                    967:        if ( !g ) {
                    968:                *rp = 0; *dnp = dn; return;
                    969:        }
                    970:        for ( n = 0, l = b; l; l = NEXT(l), n++ );
                    971:                wb = (int *)ALLOCA(n*sizeof(int));
                    972:        for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
                    973:                wb[i] = QTOS((Q)BDY(l));
                    974:        sugar = g->sugar;
                    975:        for ( d = 0; g; ) {
                    976:                for ( u = 0, i = 0; i < n; i++ ) {
                    977:                        if ( dp_redble(g,p = ps[wb[i]]) ) {
                    978:                                dp_red_mod(d,g,p,mod,&t,&u,&tdn);
                    979:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                    980:                                sugar = MAX(sugar,psugar);
                    981:                                if ( !u ) {
                    982:                                        if ( d )
                    983:                                                d->sugar = sugar;
                    984:                                        *rp = d; *dnp = dn; return;
                    985:                                } else {
                    986:                                        d = t;
                    987:                                        mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
                    988:                                }
                    989:                                break;
                    990:                        }
                    991:                }
                    992:                if ( u )
                    993:                        g = u;
                    994:                else if ( !full ) {
                    995:                        if ( g ) {
                    996:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                    997:                        }
                    998:                        *rp = g; *dnp = dn; return;
                    999:                } else {
                   1000:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   1001:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   1002:                        addmd(CO,mod,d,t,&s); d = s;
                   1003:                        dp_rest(g,&t); g = t;
                   1004:                }
                   1005:        }
                   1006:        if ( d )
                   1007:                d->sugar = sugar;
                   1008:        *rp = d; *dnp = dn;
                   1009: }
                   1010:
1.7       noro     1011: void _dp_nf_mod_destructive(b,g,ps,mod,full,rp)
                   1012: NODE b;
                   1013: DP g;
                   1014: DP *ps;
                   1015: int mod,full;
                   1016: DP *rp;
1.5       noro     1017: {
1.7       noro     1018:        DP u,p,d,s,t;
                   1019:        NODE l;
                   1020:        MP m,mr,mrd;
                   1021:        int sugar,psugar,n,h_reducible,i;
1.5       noro     1022:
1.7       noro     1023:        if ( !g ) {
                   1024:                *rp = 0; return;
1.5       noro     1025:        }
1.7       noro     1026:        sugar = g->sugar;
                   1027:        n = g->nv;
                   1028:        for ( d = 0; g; ) {
                   1029:                for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
                   1030:                        if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
                   1031:                                h_reducible = 1;
                   1032:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   1033:                                _dp_red_mod_destructive(g,p,mod,&u); g = u;
                   1034:                                sugar = MAX(sugar,psugar);
                   1035:                                if ( !g ) {
                   1036:                                        if ( d )
                   1037:                                                d->sugar = sugar;
                   1038:                                        _dptodp(d,rp); _free_dp(d); return;
                   1039:                                }
                   1040:                                break;
                   1041:                        }
                   1042:                }
                   1043:                if ( !h_reducible ) {
                   1044:                        /* head term is not reducible */
                   1045:                        if ( !full ) {
                   1046:                                if ( g )
                   1047:                                        g->sugar = sugar;
                   1048:                                _dptodp(g,rp); _free_dp(g); return;
                   1049:                        } else {
                   1050:                                m = BDY(g);
                   1051:                                if ( NEXT(m) ) {
                   1052:                                        BDY(g) = NEXT(m); NEXT(m) = 0;
                   1053:                                } else {
                   1054:                                        _FREEDP(g); g = 0;
                   1055:                                }
                   1056:                                if ( d ) {
                   1057:                                        for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
                   1058:                                        NEXT(mrd) = m;
                   1059:                                } else {
                   1060:                                        _MKDP(n,m,d);
                   1061:                                }
                   1062:                        }
                   1063:                }
1.5       noro     1064:        }
1.7       noro     1065:        if ( d )
                   1066:                d->sugar = sugar;
                   1067:        _dptodp(d,rp); _free_dp(d);
1.5       noro     1068: }
                   1069:
1.7       noro     1070: void dp_lnf_mod(p1,p2,g,mod,r1p,r2p)
                   1071: DP p1,p2;
                   1072: NODE g;
                   1073: int mod;
                   1074: DP *r1p,*r2p;
1.5       noro     1075: {
1.7       noro     1076:        DP r1,r2,b1,b2,t,s;
                   1077:        P c;
                   1078:        MQ c1,c2;
                   1079:        NODE l,b;
                   1080:        int n;
                   1081:
                   1082:        if ( !p1 ) {
                   1083:                *r1p = p1; *r2p = p2; return;
                   1084:        }
                   1085:        n = p1->nv;
                   1086:        for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
                   1087:                        if ( !r1 ) {
                   1088:                                *r1p = r1; *r2p = r2; return;
                   1089:                        }
                   1090:                        b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
                   1091:                        if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
                   1092:                                b2 = (DP)BDY(NEXT(b));
                   1093:                                invmq(mod,(MQ)BDY(b1)->c,&c1);
                   1094:                                mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
                   1095:                                mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
                   1096:                                mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
                   1097:                        }
                   1098:        }
                   1099:        *r1p = r1; *r2p = r2;
1.5       noro     1100: }
                   1101:
1.7       noro     1102: void dp_nf_tab_mod(p,tab,mod,rp)
                   1103: DP p;
                   1104: LIST *tab;
                   1105: int mod;
                   1106: DP *rp;
1.5       noro     1107: {
1.7       noro     1108:        DP s,t,u;
                   1109:        MP m;
                   1110:        DL h;
                   1111:        int i,n;
                   1112:
                   1113:        if ( !p ) {
                   1114:                *rp = p; return;
                   1115:        }
                   1116:        n = p->nv;
                   1117:        for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1118:                h = m->dl;
                   1119:                while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
                   1120:                        i++;
                   1121:                mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
                   1122:                addmd(CO,mod,s,t,&u); s = u;
                   1123:        }
                   1124:        *rp = s;
1.5       noro     1125: }
                   1126:
1.7       noro     1127: /*
                   1128:  * setting flags
                   1129:  *
                   1130:  */
                   1131:
                   1132: int create_order_spec(obj,spec)
                   1133: Obj obj;
                   1134: struct order_spec *spec;
1.5       noro     1135: {
1.7       noro     1136:        int i,j,n,s,row,col;
                   1137:        struct order_pair *l;
                   1138:        NODE node,t,tn;
                   1139:        MAT m;
                   1140:        pointer **b;
                   1141:        int **w;
1.5       noro     1142:
1.7       noro     1143:        if ( !obj || NUM(obj) ) {
                   1144:                spec->id = 0; spec->obj = obj;
                   1145:                spec->ord.simple = QTOS((Q)obj);
                   1146:                return 1;
                   1147:        } else if ( OID(obj) == O_LIST ) {
                   1148:                node = BDY((LIST)obj);
                   1149:                for ( n = 0, t = node; t; t = NEXT(t), n++ );
                   1150:                l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
                   1151:                for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
                   1152:                        tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
                   1153:                        tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
                   1154:                        s += l[i].length;
                   1155:                }
                   1156:                spec->id = 1; spec->obj = obj;
                   1157:                spec->ord.block.order_pair = l;
                   1158:                spec->ord.block.length = n; spec->nv = s;
                   1159:                return 1;
                   1160:        } else if ( OID(obj) == O_MAT ) {
                   1161:                m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
                   1162:                w = almat(row,col);
                   1163:                for ( i = 0; i < row; i++ )
                   1164:                        for ( j = 0; j < col; j++ )
                   1165:                                w[i][j] = QTOS((Q)b[i][j]);
                   1166:                spec->id = 2; spec->obj = obj;
                   1167:                spec->nv = col; spec->ord.matrix.row = row;
                   1168:                spec->ord.matrix.matrix = w;
                   1169:                return 1;
                   1170:        } else
1.5       noro     1171:                return 0;
                   1172: }
                   1173:
1.7       noro     1174: /*
                   1175:  * converters
                   1176:  *
                   1177:  */
                   1178:
                   1179: void dp_homo(p,rp)
                   1180: DP p;
                   1181: DP *rp;
1.5       noro     1182: {
1.7       noro     1183:        MP m,mr,mr0;
                   1184:        int i,n,nv,td;
                   1185:        DL dl,dlh;
1.5       noro     1186:
1.7       noro     1187:        if ( !p )
                   1188:                *rp = 0;
                   1189:        else {
                   1190:                n = p->nv; nv = n + 1;
                   1191:                m = BDY(p); td = sugard(m);
                   1192:                for ( mr0 = 0; m; m = NEXT(m) ) {
                   1193:                        NEXTMP(mr0,mr); mr->c = m->c;
                   1194:                        dl = m->dl;
                   1195:                        mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
                   1196:                        dlh->td = td;
                   1197:                        for ( i = 0; i < n; i++ )
                   1198:                                dlh->d[i] = dl->d[i];
                   1199:                        dlh->d[n] = td - dl->td;
                   1200:                }
                   1201:                NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5       noro     1202:        }
                   1203: }
                   1204:
1.7       noro     1205: void dp_dehomo(p,rp)
                   1206: DP p;
1.5       noro     1207: DP *rp;
                   1208: {
1.7       noro     1209:        MP m,mr,mr0;
                   1210:        int i,n,nv;
                   1211:        DL dl,dlh;
1.5       noro     1212:
1.7       noro     1213:        if ( !p )
                   1214:                *rp = 0;
                   1215:        else {
                   1216:                n = p->nv; nv = n - 1;
                   1217:                m = BDY(p);
                   1218:                for ( mr0 = 0; m; m = NEXT(m) ) {
                   1219:                        NEXTMP(mr0,mr); mr->c = m->c;
                   1220:                        dlh = m->dl;
                   1221:                        mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
                   1222:                        dl->td = dlh->td - dlh->d[nv];
                   1223:                        for ( i = 0; i < nv; i++ )
                   1224:                                dl->d[i] = dlh->d[i];
                   1225:                }
                   1226:                NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
                   1227:        }
1.5       noro     1228: }
                   1229:
1.7       noro     1230: void dp_mod(p,mod,subst,rp)
                   1231: DP p;
                   1232: int mod;
                   1233: NODE subst;
1.5       noro     1234: DP *rp;
                   1235: {
1.7       noro     1236:        MP m,mr,mr0;
                   1237:        P t,s,s1;
                   1238:        V v;
                   1239:        NODE tn;
1.5       noro     1240:
1.7       noro     1241:        if ( !p )
                   1242:                *rp = 0;
                   1243:        else {
                   1244:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1245:                        for ( tn = subst, s = m->c; tn; tn = NEXT(tn) ) {
                   1246:                                v = VR((P)BDY(tn)); tn = NEXT(tn);
                   1247:                                substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
                   1248:                        }
                   1249:                        ptomp(mod,s,&t);
                   1250:                        if ( t ) {
                   1251:                                NEXTMP(mr0,mr); mr->c = t; mr->dl = m->dl;
                   1252:                        }
                   1253:                }
                   1254:                if ( mr0 ) {
                   1255:                        NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                   1256:                } else
                   1257:                        *rp = 0;
                   1258:        }
1.5       noro     1259: }
                   1260:
1.7       noro     1261: void dp_rat(p,rp)
                   1262: DP p;
                   1263: DP *rp;
1.5       noro     1264: {
1.7       noro     1265:        MP m,mr,mr0;
1.5       noro     1266:
1.7       noro     1267:        if ( !p )
                   1268:                *rp = 0;
                   1269:        else {
                   1270:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   1271:                        NEXTMP(mr0,mr); mptop(m->c,&mr->c); mr->dl = m->dl;
1.5       noro     1272:                }
1.7       noro     1273:                if ( mr0 ) {
                   1274:                        NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                   1275:                } else
                   1276:                        *rp = 0;
1.5       noro     1277:        }
                   1278: }
                   1279:
                   1280:
1.7       noro     1281: void homogenize_order(old,n,new)
                   1282: struct order_spec *old,*new;
                   1283: int n;
1.5       noro     1284: {
1.7       noro     1285:        struct order_pair *l;
                   1286:        int length,nv,row,i,j;
                   1287:        int **newm,**oldm;
1.5       noro     1288:
1.7       noro     1289:        switch ( old->id ) {
                   1290:                case 0:
                   1291:                        switch ( old->ord.simple ) {
                   1292:                                case 0:
                   1293:                                        new->id = 0; new->ord.simple = 0; break;
                   1294:                                case 1:
                   1295:                                        l = (struct order_pair *)
                   1296:                                                MALLOC_ATOMIC(2*sizeof(struct order_pair));
                   1297:                                        l[0].length = n; l[0].order = 1;
                   1298:                                        l[1].length = 1; l[1].order = 2;
                   1299:                                        new->id = 1;
                   1300:                                        new->ord.block.order_pair = l;
                   1301:                                        new->ord.block.length = 2; new->nv = n+1;
                   1302:                                        break;
                   1303:                                case 2:
                   1304:                                        new->id = 0; new->ord.simple = 1; break;
                   1305:                                case 3: case 4: case 5:
                   1306:                                        new->id = 0; new->ord.simple = old->ord.simple+3;
                   1307:                                        dp_nelim = n-1; break;
                   1308:                                case 6: case 7: case 8: case 9:
                   1309:                                        new->id = 0; new->ord.simple = old->ord.simple; break;
                   1310:                                default:
                   1311:                                        error("homogenize_order : invalid input");
                   1312:                        }
                   1313:                        break;
                   1314:                case 1:
                   1315:                        length = old->ord.block.length;
                   1316:                        l = (struct order_pair *)
                   1317:                                MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
                   1318:                        bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
                   1319:                        l[length].order = 2; l[length].length = 1;
                   1320:                        new->id = 1; new->nv = n+1;
                   1321:                        new->ord.block.order_pair = l;
                   1322:                        new->ord.block.length = length+1;
                   1323:                        break;
                   1324:                case 2:
                   1325:                        nv = old->nv; row = old->ord.matrix.row;
                   1326:                        oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
                   1327:                        for ( i = 0; i <= nv; i++ )
                   1328:                                newm[0][i] = 1;
                   1329:                        for ( i = 0; i < row; i++ ) {
                   1330:                                for ( j = 0; j < nv; j++ )
                   1331:                                        newm[i+1][j] = oldm[i][j];
                   1332:                                newm[i+1][j] = 0;
                   1333:                        }
                   1334:                        new->id = 2; new->nv = nv+1;
                   1335:                        new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
                   1336:                        break;
                   1337:                default:
                   1338:                        error("homogenize_order : invalid input");
1.5       noro     1339:        }
1.7       noro     1340: }
                   1341:
                   1342: void qltozl(w,n,dvr)
                   1343: Q *w,*dvr;
                   1344: int n;
                   1345: {
                   1346:        N nm,dn;
                   1347:        N g,l1,l2,l3;
                   1348:        Q c,d;
                   1349:        int i;
                   1350:        struct oVECT v;
1.5       noro     1351:
                   1352:        for ( i = 0; i < n; i++ )
1.7       noro     1353:                if ( w[i] && !INT(w[i]) )
                   1354:                        break;
                   1355:        if ( i == n ) {
                   1356:                v.id = O_VECT; v.len = n; v.body = (pointer *)w;
                   1357:                igcdv(&v,dvr); return;
                   1358:        }
                   1359:        c = w[0]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
                   1360:        for ( i = 1; i < n; i++ ) {
                   1361:                c = w[i]; l1 = INT(c) ? ONEN : DN(c);
                   1362:                gcdn(nm,NM(c),&g); nm = g;
                   1363:                gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
1.5       noro     1364:        }
1.7       noro     1365:        if ( UNIN(dn) )
                   1366:                NTOQ(nm,1,d);
                   1367:        else
                   1368:                NDTOQ(nm,dn,1,d);
                   1369:        *dvr = d;
                   1370: }
1.5       noro     1371:
1.7       noro     1372: int comp_nm(a,b)
                   1373: Q *a,*b;
                   1374: {
                   1375:        return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
                   1376: }
                   1377:
                   1378: void sortbynm(w,n)
                   1379: Q *w;
                   1380: int n;
                   1381: {
                   1382:        qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
                   1383: }
1.5       noro     1384:
                   1385:
1.7       noro     1386: /*
                   1387:  * simple operations
                   1388:  *
                   1389:  */
1.5       noro     1390:
1.7       noro     1391: int dp_redble(p1,p2)
                   1392: DP p1,p2;
                   1393: {
                   1394:        int i,n;
                   1395:        DL d1,d2;
1.5       noro     1396:
1.7       noro     1397:        d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                   1398:        if ( d1->td < d2->td )
                   1399:                return 0;
                   1400:        else {
                   1401:                for ( i = 0, n = p1->nv; i < n; i++ )
                   1402:                        if ( d1->d[i] < d2->d[i] )
                   1403:                                return 0;
                   1404:                return 1;
1.5       noro     1405:        }
                   1406: }
                   1407:
1.7       noro     1408: void dp_subd(p1,p2,rp)
1.5       noro     1409: DP p1,p2;
                   1410: DP *rp;
                   1411: {
1.7       noro     1412:        int i,n;
1.5       noro     1413:        DL d1,d2,d;
                   1414:        MP m;
1.7       noro     1415:        DP s;
1.5       noro     1416:
                   1417:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.7       noro     1418:        NEWDL(d,n); d->td = d1->td - d2->td;
1.5       noro     1419:        for ( i = 0; i < n; i++ )
1.7       noro     1420:                d->d[i] = d1->d[i]-d2->d[i];
                   1421:        NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
                   1422:        MKDP(n,m,s); s->sugar = d->td;
                   1423:        *rp = s;
                   1424: }
                   1425:
                   1426: void dltod(d,n,rp)
                   1427: DL d;
                   1428: int n;
                   1429: DP *rp;
                   1430: {
                   1431:        MP m;
                   1432:        DP s;
                   1433:
                   1434:        NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
                   1435:        MKDP(n,m,s); s->sugar = d->td;
                   1436:        *rp = s;
1.5       noro     1437: }
                   1438:
                   1439: void dp_hm(p,rp)
                   1440: DP p;
                   1441: DP *rp;
                   1442: {
                   1443:        MP m,mr;
                   1444:
                   1445:        if ( !p )
                   1446:                *rp = 0;
                   1447:        else {
                   1448:                m = BDY(p);
                   1449:                NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
                   1450:                MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td;  /* XXX */
                   1451:        }
                   1452: }
                   1453:
                   1454: void dp_rest(p,rp)
                   1455: DP p,*rp;
                   1456: {
                   1457:        MP m;
                   1458:
                   1459:        m = BDY(p);
                   1460:        if ( !NEXT(m) )
                   1461:                *rp = 0;
                   1462:        else {
                   1463:                MKDP(p->nv,NEXT(m),*rp);
                   1464:                if ( *rp )
                   1465:                        (*rp)->sugar = p->sugar;
                   1466:        }
                   1467: }
                   1468:
                   1469: DL lcm_of_DL(nv,dl1,dl2,dl)
                   1470: int nv;
                   1471: DL dl1,dl2;
                   1472: register DL dl;
                   1473: {
                   1474:        register int n, *d1, *d2, *d, td;
                   1475:
                   1476:        if ( !dl ) NEWDL(dl,nv);
                   1477:        d = dl->d,  d1 = dl1->d,  d2 = dl2->d;
                   1478:        for ( td = 0, n = nv; --n >= 0; d1++, d2++, d++ )
                   1479:                td += (*d = *d1 > *d2 ? *d1 : *d2 );
                   1480:        dl->td = td;
                   1481:        return dl;
                   1482: }
                   1483:
                   1484: int dl_equal(nv,dl1,dl2)
                   1485: int nv;
                   1486: DL dl1, dl2;
                   1487: {
                   1488:     register int *d1, *d2, n;
                   1489:
                   1490:     if ( dl1->td != dl2->td ) return 0;
                   1491:     for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
                   1492:         if ( *d1 != *d2 ) return 0;
                   1493:     return 1;
                   1494: }
                   1495:
                   1496: int dp_nt(p)
                   1497: DP p;
                   1498: {
                   1499:        int i;
                   1500:        MP m;
                   1501:
                   1502:        if ( !p )
                   1503:                return 0;
                   1504:        else {
                   1505:                for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
                   1506:                return i;
                   1507:        }
                   1508: }
                   1509:

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