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

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

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