[BACK]Return to EZ.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Annotation of OpenXM_contrib2/asir2000/engine/EZ.c, Revision 1.1.1.1

1.1       noro        1: /* $OpenXM: OpenXM/src/asir99/engine/EZ.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
                      2: #include "ca.h"
                      3:
                      4: void ezgcdp(vl,p1,p2,pr)
                      5: VL vl;
                      6: P p1,p2;
                      7: P *pr;
                      8: {
                      9:        P f1,f2;
                     10:        Q c;
                     11:
                     12:        if ( !p1 )
                     13:                if ( !p2 )
                     14:                        *pr = 0;
                     15:                else
                     16:                        *pr = p2;
                     17:        else if ( !p2 )
                     18:                *pr = p1;
                     19:        else {
                     20:                if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
                     21:                        error("ezgcdp : invalid argument");
                     22:                ptozp(p1,1,&c,&f1);
                     23:                ptozp(p2,1,&c,&f2);
                     24:                ezgcdpz(vl,f1,f2,pr);
                     25:        }
                     26: }
                     27:
                     28: /*
                     29:  * p1, p2 : integer coefficient
                     30:  * returns gcd(pp(p1),pp(p2)) * gcd(cont(p1),cont(p2))
                     31:  */
                     32:
                     33: void ezgcdpz(vl,p1,p2,pr)
                     34: VL vl;
                     35: P p1,p2,*pr;
                     36: {
                     37:        P f1,f2,t,s,g,gcd;
                     38:        P fc1,fc2,cont;
                     39:        VL tvl,svl;
                     40:        V v1,v2;
                     41:        DCP dc,dc1,dc2;
                     42:        N cn;
                     43:        Q cq;
                     44:        Q c1,c2,c;
                     45:        P g1,g2;
                     46:        P mgcd;
                     47:        extern int nez;
                     48:        P pm[2];
                     49:
                     50:        if ( nez ) {
                     51:                pm[0] = p1; pm[1] = p2; nezgcdnpz(vl,pm,2,pr); return;
                     52:        }
                     53:        monomialfctr(vl,p1,&f1,&dc1); p1 = f1;
                     54:        monomialfctr(vl,p2,&f2,&dc2); p2 = f2;
                     55:        for ( mgcd = (P)ONE; dc1; dc1 = NEXT(dc1) )
                     56:                for ( v1 = VR(dc1->c), dc = dc2; dc; dc = NEXT(dc) )
                     57:                        if ( v1 == VR(dc->c) ) {
                     58:                                pwrp(vl,dc->c,cmpq(dc1->d,dc->d)>=0?dc->d:dc1->d,&t);
                     59:                                mulp(vl,mgcd,t,&s); mgcd = s; break;
                     60:                        }
                     61:
                     62:        if ( NUM(p1) ) {
                     63:                if ( NUM(p2) ) {
                     64:                        gcdn(NM((Q)p1),NM((Q)p2),&cn); NTOQ(cn,1,c);
                     65:                } else {
                     66:                        ptozp(p2,1,&c2,&t);
                     67:                        gcdn(NM((Q)p1),NM(c2),&cn); NTOQ(cn,1,c);
                     68:                }
                     69:                mulp(vl,(P)c,mgcd,pr);
                     70:                return;
                     71:        } else if ( NUM(p2) ) {
                     72:                ptozp(p1,1,&c1,&t);
                     73:                gcdn(NM(c1),NM((Q)p2),&cn); NTOQ(cn,1,c);
                     74:                mulp(vl,(P)c,mgcd,pr);
                     75:                return;
                     76:        }
                     77:
                     78:        ptozp(p1,1,&c1,&g1); ptozp(p2,1,&c2,&g2);
                     79:        gcdn(NM(c1),NM(c2),&cn); NTOQ(cn,1,cq);
                     80:
                     81:        if ( ( v1 = VR(f1) ) != ( v2 = VR(f2) ) ) {
                     82:                while ( v1 != VR(vl) && v2 != VR(vl) )
                     83:                        vl = NEXT(vl);
                     84:                if ( v1 == VR(vl) ) {
                     85:                        pcp(vl,f1,&g,&fc1);
                     86:                        ezgcdpz(vl,fc1,f2,&t);
                     87:                } else {
                     88:                        pcp(vl,f2,&g,&fc2);
                     89:                        ezgcdpz(vl,fc2,f1,&t);
                     90:                }
                     91:                mulp(vl,t,(P)cq,&s); mulp(vl,s,mgcd,pr);
                     92:                return;
                     93:        }
                     94:
                     95:        pcp(vl,f1,&g1,&fc1); pcp(vl,f2,&g2,&fc2);
                     96:        ezgcdpz(vl,fc1,fc2,&cont);
                     97:        clctv(vl,g1,&tvl); clctv(vl,g2,&svl);
                     98:        if ( !NEXT(tvl) && !NEXT(svl) ) {
                     99:                uezgcdpz(vl,g1,g2,&t);
                    100:                mulp(vl,t,cont,&s); mulp(vl,s,(P)cq,&t); mulp(vl,t,mgcd,pr);
                    101:                return;
                    102:        }
                    103:
                    104:        if ( homdeg(g1) > homdeg(g2) ) {
                    105:                t = g1; g1 = g2; g2 = t;
                    106:        }
                    107:        sqfrp(vl,g1,&dc);
                    108:        for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) {
                    109:                if ( NUM(COEF(dc)) )
                    110:                        continue;
                    111:                ezgcdpp(vl,dc,g,&t);
                    112:                if ( NUM(t) )
                    113:                        continue;
                    114:                mulp(vl,gcd,t,&s); gcd = s;
                    115:                divsp(vl,g,t,&s); g = s;
                    116:        }
                    117:        mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,&s); mulp(vl,s,mgcd,pr);
                    118: }
                    119:
                    120: void uezgcdpz(vl,p1,p2,pr)
                    121: VL vl;
                    122: P p1,p2,*pr;
                    123: {
                    124:        P g1,g2,gcd,t,s,g;
                    125:        DCP dc;
                    126:        N cn;
                    127:        Q c1,c2,cq;
                    128:        P f1,f2;
                    129:
                    130:        if ( NUM(p1) ) {
                    131:                if ( NUM(p2) ) {
                    132:                        gcdn(NM((Q)p1),NM((Q)p2),&cn); NTOQ(cn,1,cq),*pr = (P)cq;
                    133:                        return;
                    134:                } else {
                    135:                        ptozp(p2,1,&c2,&t);
                    136:                        gcdn(NM((Q)p1),NM(c2),&cn); NTOQ(cn,1,cq),*pr = (P)cq;
                    137:                        return;
                    138:                }
                    139:        } else if ( NUM(p2) ) {
                    140:                ptozp(p1,1,&c1,&t);
                    141:                gcdn(NM(c1),NM((Q)p2),&cn); NTOQ(cn,1,cq),*pr = (P)cq;
                    142:                return;
                    143:        }
                    144:
                    145:        ptozp(p1,1,&c1,&f1); ptozp(p2,1,&c2,&f2);
                    146:        gcdn(NM(c1),NM(c2),&cn); NTOQ(cn,1,cq);
                    147:        if ( UDEG(f1) > UDEG(f2) ) {
                    148:                g1 = f2; g2 = f1;
                    149:        } else {
                    150:                g1 = f1; g2 = f2;
                    151:        }
                    152:        usqp(g1,&dc);
                    153:        for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) {
                    154:                if ( NUM(COEF(dc)) )
                    155:                        continue;
                    156:                uezgcdpp(dc,g,&t);
                    157:                if ( NUM(t) )
                    158:                        continue;
                    159:                mulp(CO,gcd,t,&s); gcd = s;
                    160:                divsp(CO,g,t,&s); g = s;
                    161:        }
                    162:        mulp(vl,gcd,(P)cq,pr);
                    163: }
                    164:
                    165: /*
                    166:        dc : dc pair c^d
                    167:        r = gcd(c^d,f)
                    168: */
                    169:
                    170: void ezgcdpp(vl,dc,f,r)
                    171: VL vl;
                    172: DCP dc;
                    173: P f;
                    174: P *r;
                    175: {
                    176:        int k;
                    177:        P g,h,t,s,gcd;
                    178:
                    179:        if ( NUM(f) ) {
                    180:                *r = (P)ONE;
                    181:                return;
                    182:        }
                    183:        k = QTOS(DEG(dc)) - 1;
                    184: /*     ezgcd1p(vl,COEF(dc),f,&gcd); */
                    185:        ezgcdnp(vl,COEF(dc),&f,1,&gcd);
                    186:        g = gcd; divsp(vl,f,gcd,&h);
                    187:        for ( ; k > 0; k-- ) {
                    188: /*             ezgcd1p(vl,g,h,&t); */
                    189:                ezgcdnp(vl,g,&h,1,&t);
                    190:                if ( NUM(t) )
                    191:                        break;
                    192:                mulp(vl,gcd,t,&s); gcd = s;
                    193:                divsp(vl,h,t,&s); h = s;
                    194:        }
                    195:        *r = gcd;
                    196: }
                    197:
                    198: void ezgcd1p(vl,p0,p1,gcdp)
                    199: VL vl;
                    200: P p0,p1;
                    201: P *gcdp;
                    202: {
                    203:        VL nvl,tvl,vl0,vl1,vl2;
                    204:        P f0,f1,t;
                    205:
                    206:        clctv(vl,p0,&vl0); clctv(vl,p1,&vl1); mergev(vl,vl0,vl1,&vl2);
                    207:        minchdegp(vl,p0,&tvl,&t); reordvar(vl2,tvl->v,&nvl);
                    208:        reorderp(nvl,vl,p0,&f0); reorderp(nvl,vl,p1,&f1);
                    209:        ezgcdnp(nvl,f0,&f1,1,&t);
                    210:        reorderp(vl,nvl,t,gcdp);
                    211: }
                    212:
                    213: void uezgcdpp(dc,f,r)
                    214: DCP dc;
                    215: P f;
                    216: P *r;
                    217: {
                    218:        int k;
                    219:        P g,h,t,s,gcd;
                    220:
                    221:        if ( NUM(f) ) {
                    222:                *r = (P)ONE;
                    223:                return;
                    224:        }
                    225:
                    226:        k = QTOS(DEG(dc)) - 1;
                    227:        uezgcd1p(COEF(dc),f,&gcd);
                    228:        g = gcd; divsp(CO,f,gcd,&h);
                    229:        for ( ; k > 0; k-- ) {
                    230:                uezgcd1p(g,h,&t);
                    231:                if ( NUM(t) )
                    232:                        break;
                    233:                mulp(CO,gcd,t,&s); gcd = s;
                    234:                divsp(CO,h,t,&s); h = s;
                    235:        }
                    236:        *r = gcd;
                    237: }
                    238:
                    239: /*
                    240:  * pr = gcd(p0,ps[0],...,ps[m-1])
                    241:  *
                    242:  * p0 is already square-free and primitive.
                    243:  * ps[i] is at least primitive.
                    244:  *
                    245:  */
                    246:
                    247: void ezgcdnp(vl,p0,ps,m,pr)
                    248: VL vl;
                    249: int m;
                    250: P p0,*ps,*pr;
                    251: {
                    252:        /* variables */
                    253:        P p00,g,h,g0,h0,a0,b0;
                    254:        P lp0,lgp0,lp00,lg,lg0,cbd,tq,t;
                    255:        P *tps;
                    256:        VL nvl1,nvl2,nvl,tvl;
                    257:        V v;
                    258:        int i,j,k,d0,dg,dg0,dmin;
                    259:        VN vn0,vn1,vnt;
                    260:        int nv,flag;
                    261:
                    262:        /* set-up */
                    263:        if ( NUM(p0) ) {
                    264:                *pr = (P) ONE;
                    265:                return;
                    266:        }
                    267:        for ( v = VR(p0), i = 0; i < m; i++ )
                    268:                if ( NUM(ps[i]) || (v != VR(ps[i])) ) {
                    269:                        *pr = (P)ONE;
                    270:                        return;
                    271:                }
                    272:        tps = (P *) ALLOCA(m*sizeof(P));
                    273:        for ( i = 0; i < m; i++ )
                    274:                tps[i] = ps[i];
                    275:        sortplist(tps,m);
                    276:        /* deg(tps[0]) <= deg(tps[1]) <= ... */
                    277:
                    278:        if ( !cmpq(DEG(DC(p0)),ONE) ) {
                    279:                if ( divcheck(vl,tps,m,(P)ONE,p0) )
                    280:                        *pr = p0;
                    281:                else
                    282:                        *pr = (P)ONE;
                    283:                return;
                    284:        }
                    285:
                    286:        lp0 = LC(p0); lg = lp0; dmin = d0 = deg(v,p0);
                    287:        clctv(vl,p0,&nvl);
                    288:        for ( i = 0; i < m; i++ ) {
                    289:                clctv(vl,tps[i],&nvl1); mergev(vl,nvl,nvl1,&nvl2);
                    290:                nvl = nvl2;
                    291:                ezgcdpz(nvl,lg,LC(tps[i]),&t); lg = t;
                    292:        }
                    293:
                    294:        mulp(nvl,p0,lg,&lgp0);
                    295:        k = dbound(v,lgp0) + 1;
                    296:        cbound(nvl,lgp0,(Q *)&cbd);
                    297:        for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++ );
                    298:        W_CALLOC(nv,struct oVN,vn0); W_CALLOC(nv,struct oVN,vnt);
                    299:        W_CALLOC(nv,struct oVN,vn1);
                    300:        for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
                    301:                vn1[i].v = vn0[i].v = tvl->v;
                    302:
                    303:        /* main loop */
                    304:        for ( dg = deg(v,tps[0]) + 1; ; next(vn0) )
                    305:                do {
                    306:                        for ( i = 0, j = 0; vn0[i].v; i++ )
                    307:                                if ( vn0[i].n ) vnt[j++].v = (V)i;
                    308:                        vnt[j].n = 0;
                    309:
                    310:                        /* find b s.t. LC(p0)(b), LC(tps[i])(b) != 0 */
                    311:                        mulsgn(vn0,vnt,j,vn1);
                    312:                        substvp(nvl,p0,vn1,&p00);
                    313:                        flag = (!zerovpchk(nvl,lp0,vn1) && sqfrchk(p00));
                    314:                        for ( i = 0; flag && (i < m); i++ )
                    315:                                flag &= (!zerovpchk(nvl,LC(tps[i]),vn1));
                    316:                        if ( !flag )
                    317:                                continue;
                    318:
                    319:                        /* substitute y -> b */
                    320:                        substvp(nvl,lg,vn1,&lg0);
                    321:                        lp00 = LC(p00);
                    322:                        /* extended-gcd in 1-variable */
                    323:                        uexgcdnp(nvl,p00,tps,m,vn1,(Q)cbd,&g0,&h0,&a0,&b0,(Q *)&tq);
                    324:                        if ( NUM(g0) ) {
                    325:                                *pr = (P)ONE;
                    326:                                return;
                    327:                        } else if ( dg > ( dg0 = deg(v,g0) ) ) {
                    328:                                dg = dg0;
                    329:                                if ( dg == d0 ) {
                    330:                                        if ( divcheck(nvl,tps,m,lp0,p0) ) {
                    331:                                                *pr = p0;
                    332:                                                return;
                    333:                                        }
                    334:                                } else if ( dg == deg(v,tps[0]) ) {
                    335:                                        if ( divtpz(nvl,p0,tps[0],&t) &&
                    336:                                                divcheck(nvl,tps+1,m-1,LC(tps[0]),tps[0]) ) {
                    337:                                                *pr = tps[0];
                    338:                                                return;
                    339:                                        } else
                    340:                                                break;
                    341:                                } else {
                    342:                                        henmv(nvl,vn1,lgp0,g0,h0,a0,b0,
                    343:                                                lg,lp0,lg0,lp00,(Q)tq,k,&g,&h);
                    344:                                        if ( divtpz(nvl,lgp0,g,&t) &&
                    345:                                                divcheck(nvl,tps,m,lg,g) ) {
                    346:                                                pcp(nvl,g,pr,&t);
                    347:                                                return;
                    348:                                        }
                    349:                                }
                    350:                        }
                    351:                } while ( !nextbin(vnt,j) );
                    352: }
                    353:
                    354: /*
                    355:        f : sqfr
                    356: */
                    357:
                    358: void uezgcd1p(f,g,r)
                    359: P f,g;
                    360: P *r;
                    361: {
                    362:        UM wf,wf1,wf2,wg,wg1,wgcd,wcof;
                    363:        int d1,d2,dg,m,index,n;
                    364:        ML list;
                    365:        DCP dc;
                    366:        P t;
                    367:        Q c;
                    368:
                    369:        if ( NUM(f) || NUM(g) ) {
                    370:                *r = (P)ONE;
                    371:                return;
                    372:        }
                    373:        ptozp(f,1,&c,&t); f = t; ptozp(g,1,&c,&t); g = t;
                    374:        d1 = UDEG(f); d2 = UDEG(g);
                    375:        n = MAX(d1,d2)+1;
                    376:        wf = W_UMALLOC(n); wf1 = W_UMALLOC(n);
                    377:        wf2 = W_UMALLOC(n); wg = W_UMALLOC(n);
                    378:        wg1 = W_UMALLOC(n); wgcd = W_UMALLOC(n);
                    379:        wcof = W_UMALLOC(n);
                    380:
                    381:        for ( index = INDEX, dg = MIN(d1,d2)+1; ; ) {
                    382:                m = sprime[index++];
                    383:                if ( !rem(NM((Q)UCOEF(f)),m) || !rem(NM((Q)UCOEF(g)),m))
                    384:                        continue;
                    385:                ptoum(m,f,wf); cpyum(wf,wf1);
                    386:                diffum(m,wf1,wf2); gcdum(m,wf1,wf2,wgcd);
                    387:                if ( DEG(wgcd) > 0 )
                    388:                        continue;
                    389:                ptoum(m,g,wg); cpyum(wg,wg1); cpyum(wf,wf1);
                    390:                gcdum(m,wf1,wg1,wgcd);
                    391:                if ( dg <= DEG(wgcd) )
                    392:                        continue;
                    393:                else
                    394:                        dg = DEG(wgcd);
                    395:                if ( dg == 0 ) {
                    396:                        *r = (P)ONE;
                    397:                        return;
                    398:                } else if ( dg == d1 ) {
                    399:                        if ( divtpz(CO,g,f,&t) ) {
                    400:                                *r = f;
                    401:                                return;
                    402:                        }
                    403:                } else if ( dg == d2 ) {
                    404:                        if ( divtpz(CO,f,g,&t) ) {
                    405:                                *r = g;
                    406:                                return;
                    407:                        }
                    408:                } else {
                    409:                        divum(m,wf,wgcd,wcof);
                    410:                        ezgcdhensel(f,m,wcof,wgcd,&list);
                    411:                        dtest(f,list,1,&dc);
                    412:                        if ( NEXT(dc) ) {
                    413:                                if ( divtpz(CO,g,COEF(dc),&t) ) {
                    414:                                        *r = COEF(dc);
                    415:                                        return;
                    416:                                }
                    417:                        }
                    418:                }
                    419:        }
                    420: }
                    421:
                    422: void uexgcdnp(vl,p1,ps,n,vn,cbd,g,h,a,b,q)
                    423: VL vl;
                    424: P p1;
                    425: P *ps;
                    426: int n;
                    427: VN vn;
                    428: Q cbd;
                    429: P *g,*h,*a,*b;
                    430: Q *q;
                    431: {
                    432:        UM wg,wh,wg1,wh1,wgcd,wt;
                    433:        P t,s,gcd,cof,gk,hk,ak,bk;
                    434:        Q c,tq,bb,tq1,qk;
                    435:        int mod,k,i,index,d;
                    436:
                    437:        ptozp(p1,1,&c,&gcd);
                    438:        for ( i = 0; i < n; i++ ) {
                    439:                substvp(vl,ps[i],vn,&t);
                    440:                uezgcd1p(gcd,t,&s);
                    441:                if ( NUM(s) ) {
                    442:                        *g = (P)ONE; *h = p1; *a = 0; *b = (P)ONE;
                    443:                        return;
                    444:                } else
                    445:                        gcd = s;
                    446:        }
                    447:
                    448:        if ( !dcomp(p1,gcd) ) {
                    449:                *g = p1; *h = (P)ONE; *a = 0; *b = (P)ONE;
                    450:                return;
                    451:        }
                    452:
                    453:        divsp(CO,p1,gcd,&cof);
                    454:        d = UDEG(p1)+1;
                    455:        wg = W_UMALLOC(d); wh = W_UMALLOC(d);
                    456:        wg1 = W_UMALLOC(d); wh1 = W_UMALLOC(d);
                    457:        wgcd = W_UMALLOC(d); wt = W_UMALLOC(d);
                    458:        for ( index = INDEX; ; ) {
                    459:                mod = sprime[index++];
                    460:                if ( !rem(NM((Q)LC(p1)),mod) )
                    461:                        continue;
                    462:                ptoum(mod,gcd,wg); ptoum(mod,cof,wh);
                    463:                cpyum(wg,wg1); cpyum(wh,wh1);
                    464:                gcdum(mod,wg,wh,wt); cpyum(wt,wgcd);
                    465:                if ( DEG(wgcd) > 0 )
                    466:                        continue;
                    467:                STOQ(mod,tq); addq(cbd,cbd,&bb);
                    468:                for ( k = 1; cmpq(tq,bb) < 0; k++ ) {
                    469:                        mulq(tq,tq,&tq1); tq = tq1;
                    470:                }
                    471:                henzq(p1,gcd,wg1,cof,wh1,mod,k,&gk,&hk,&bk,&ak,&qk);
                    472:                break;
                    473:        }
                    474:        *g = gk; *h = hk; *a = ak; *b = bk; *q = tq;
                    475: }
                    476:
                    477: void ezgcdhensel(f,mod,cof,gcd,listp)
                    478: P f;
                    479: int mod;
                    480: UM gcd,cof;
                    481: ML *listp;
                    482: {
                    483:        register int i,j;
                    484:        int q,n,bound;
                    485:        int *p;
                    486:        int **pp;
                    487:        ML blist,clist,bqlist,cqlist,rlist;
                    488:        UM *b;
                    489:        LUM fl,tl;
                    490:        LUM *l;
                    491:        LUM LUMALLOC();
                    492:
                    493:        blist = W_MLALLOC(2);
                    494:        blist->n = 2; blist->mod = mod;
                    495:        blist->c[0] = (pointer)cof; blist->c[1] = (pointer)gcd;
                    496:        gcdgen(f,blist,&clist);
                    497:        henprep(f,blist,clist,&bqlist,&cqlist);
                    498:        n = bqlist->n; q = bqlist->mod;
                    499:        bqlist->bound = cqlist->bound = bound = mignotte(q,f);
                    500:
                    501:        if ( bound == 1 ) {
                    502:                *listp = rlist = MLALLOC(n);
                    503:                rlist->n = n;
                    504:                rlist->mod = q;
                    505:                rlist->bound = bound;
                    506:
                    507:                for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c;
                    508:                        i < n; i++ ) {
                    509:                        tl = LUMALLOC(DEG(b[i]),1);
                    510:                        l[i] = tl;
                    511:                        p = COEF(b[i]);
                    512:
                    513:                        for ( j = 0, pp = COEF(tl);
                    514:                                  j <= DEG(tl); j++ )
                    515:                                        pp[j][0] = p[j];
                    516:                }
                    517:        } else {
                    518:                W_LUMALLOC((int)UDEG(f),bound,fl);
                    519:                ptolum(q,bound,f,fl);
                    520:                henmain(fl,bqlist,cqlist,listp);
                    521:        }
                    522: }
                    523:
                    524: void ezgcdnpz(vl,ps,m,pr)
                    525: VL vl;
                    526: P *ps,*pr;
                    527: int m;
                    528: {
                    529:        P t,s,gcd;
                    530:        P cont;
                    531:        VL tvl,svl,nvl;
                    532:        int i;
                    533:        DCP dc;
                    534:        N cn;
                    535:        Q cq;
                    536:        P *pl,*ppl,*pcl;
                    537:        Q *cl;
                    538:        V v;
                    539:
                    540:        pl = (P *)ALLOCA((m+1)*sizeof(P));
                    541:        cl = (Q *)ALLOCA((m+1)*sizeof(Q));
                    542:        for ( i = 0; i < m; i++ )
                    543:                ptozp(ps[i],1,&cl[i],&pl[i]);
                    544:        for ( i = 1, cq = cl[0]; i < m; i++ ) {
                    545:                gcdn(NM(cl[i]),NM(cq),&cn); NTOQ(cn,1,cq);
                    546:        }
                    547:        for ( i = 0; i < m; i++ )
                    548:                if ( NUM(pl[i]) ) {
                    549:                        *pr = (P)cq;
                    550:                        return;
                    551:                }
                    552:
                    553:        for ( i = 0, nvl = 0; i < m; i++ ) {
                    554:                clctv(vl,ps[i],&tvl); mergev(vl,nvl,tvl,&svl); nvl = svl;
                    555:        }
                    556:
                    557:        ppl = (P *)ALLOCA((m+1)*sizeof(P));
                    558:        pcl = (P *)ALLOCA((m+1)*sizeof(P));
                    559:        for ( i = 0, v = nvl->v; i < m; i++ )
                    560:                if ( v == VR(pl[i]) )
                    561:                        pcp(nvl,pl[i],&ppl[i],&pcl[i]);
                    562:                else {
                    563:                        ppl[i] = (P)ONE; pcl[i] = pl[i];
                    564:                }
                    565:        ezgcdnpz(nvl,pcl,m,&cont);
                    566:
                    567:        for ( i = 0; i < m; i++ )
                    568:                if ( NUM(ppl[i]) ) {
                    569:                        mulp(nvl,cont,(P)cq,pr);
                    570:                        return;
                    571:                }
                    572:        sortplist(ppl,m);
                    573:        sqfrp(nvl,ppl[0],&dc);
                    574:        for ( dc = NEXT(dc), gcd = (P)ONE; dc; dc = NEXT(dc) ) {
                    575:                if ( NUM(COEF(dc)) )
                    576:                        continue;
                    577:                ezgcdnpp(vl,dc,ppl+1,m-1,&t);
                    578:                if ( NUM(t) )
                    579:                        continue;
                    580:                mulp(vl,gcd,t,&s); gcd = s;
                    581:        }
                    582:        mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,pr);
                    583: }
                    584:
                    585: void ezgcdnpp(vl,dc,pl,m,r)
                    586: VL vl;
                    587: DCP dc;
                    588: P *pl;
                    589: int m;
                    590: P *r;
                    591: {
                    592:        int i,k;
                    593:        P g,t,s,gcd;
                    594:        P *pm;
                    595:
                    596:        ezgcdnp(vl,COEF(dc),pl,m,&gcd);
                    597:        if ( NUM(gcd) ) {
                    598:                *r = (P)ONE;
                    599:                return;
                    600:        }
                    601:        pm = (P *) ALLOCA((m+1)*sizeof(P));
                    602:        for ( i = 0; i < m; i++ ) {
                    603:                divsp(vl,pl[i],gcd,&pm[i]);
                    604:                if ( NUM(pm[i]) ) {
                    605:                        *r = gcd;
                    606:                        return;
                    607:                }
                    608:        }
                    609:        for ( g = gcd, k = QTOS(DEG(dc)) - 1; k > 0; k-- ) {
                    610:                ezgcdnp(vl,g,pm,m,&t);
                    611:                if ( NUM(t) )
                    612:                        break;
                    613:                mulp(vl,gcd,t,&s); gcd = s;
                    614:                for ( i = 0; i < m; i++ ) {
                    615:                        divsp(vl,pm[i],t,&s);
                    616:                        if ( NUM(s) ) {
                    617:                                *r = gcd;
                    618:                                return;
                    619:                        }
                    620:                        pm[i] = s;
                    621:                }
                    622:        }
                    623:        *r = gcd;
                    624: }
                    625:
                    626: void pcp(vl,p,pp,c)
                    627: VL vl;
                    628: P p;
                    629: P *pp,*c;
                    630: {
                    631:        P f,g,n;
                    632:        DCP dc;
                    633:        P *pl;
                    634:        int i,m;
                    635:        extern int nez;
                    636:
                    637:        if ( NUM(p) ) {
                    638:                *c = p;
                    639:                *pp = (P)ONE;
                    640:        } else {
                    641:                for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
                    642:                if ( m == 1 ) {
                    643:                        *c = COEF(DC(p));
                    644:                        divsp(vl,p,*c,pp);
                    645:                        return;
                    646:                }
                    647:                ptozp(p,1,(Q *)&n,&f);
                    648:                pl = (P *)ALLOCA((m+1)*sizeof(P));
                    649:                for ( i = 0, dc = DC(f); i < m; dc = NEXT(dc), i++ )
                    650:                        pl[i] = COEF(dc);
                    651:                if ( nez )
                    652:                        nezgcdnpz(vl,pl,m,&g);
                    653:                else
                    654:                        ezgcdnpz(vl,pl,m,&g);
                    655:                mulp(vl,g,n,c); divsp(vl,f,g,pp);
                    656:        }
                    657: }
                    658:
                    659: #if 0
                    660: pcp(vl,p,pp,c)
                    661: VL vl;
                    662: P p;
                    663: P *pp,*c;
                    664: {
                    665:        P f,g,n,t;
                    666:        DCP dc;
                    667:
                    668:        if ( NUM(p) ) {
                    669:                *c = p;
                    670:                *pp = (P)ONE;
                    671:        } else {
                    672:                ptozp(p,1,&n,&f);
                    673:                for ( g = LC(f), dc = NEXT(DC(f)); dc; dc = NEXT(dc) ) {
                    674:                        ezgcdpz(vl,g,COEF(dc),&t); g = t;
                    675:                }
                    676:                mulp(vl,g,n,c);
                    677:                divsp(vl,f,g,pp);
                    678:        }
                    679: }
                    680: #endif
                    681:
                    682: int divcheck(vl,ps,m,l,p)
                    683: VL vl;
                    684: int m;
                    685: P *ps,l,p;
                    686: {
                    687:        int i;
                    688:        P r,t;
                    689:
                    690:        for ( i = 0; i < m; i++ ) {
                    691:                mulp(vl,ps[i],l,&t);
                    692:                if ( !divtpz(vl,t,p,&r) )
                    693:                        return ( 0 );
                    694:        }
                    695:        return ( 1 );
                    696: }
                    697:
                    698: void sortplist(plist,n)
                    699: P *plist;
                    700: int n;
                    701: {
                    702:        int i,j;
                    703:        P t;
                    704:
                    705:        for ( i = 0; i < n; i++ )
                    706:                for ( j = i + 1; j < n; j++ )
                    707:                        if ( cmpq(DEG(DC(plist[i])),DEG(DC(plist[j]))) > 0 ) {
                    708:                                t = plist[i]; plist[i] = plist[j]; plist[j] = t;
                    709:                        }
                    710: }
                    711:
                    712: int dcomp(p1,p2)
                    713: P p1,p2;
                    714: {
                    715:        return ( cmpq(DEG(DC(p1)),DEG(DC(p2))) );
                    716: }

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