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

Annotation of OpenXM_contrib2/asir2000/engine/F.c, Revision 1.7

1.3       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.4       noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.3       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.7     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/engine/F.c,v 1.6 2001/04/20 02:34:21 noro Exp $
1.3       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include <math.h>
                     52:
                     53: void homfctr();
                     54:
                     55: void fctrp(vl,f,dcp)
                     56: VL vl;
                     57: P f;
                     58: DCP *dcp;
                     59: {
                     60:        VL nvl;
                     61:        DCP dc;
                     62:
                     63:        if ( !f || NUM(f) ) {
                     64:                NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
                     65:                NEXT(dc) = 0; *dcp = dc;
                     66:                return;
                     67:        } else if ( !qpcheck((Obj)f) )
                     68:                error("fctrp : invalid argument");
                     69:        else {
                     70:                clctv(vl,f,&nvl);
                     71:                if ( !NEXT(nvl) )
                     72:                        ufctr(f,1,dcp);
                     73:                else
                     74:                        mfctr(nvl,f,dcp);
                     75:        }
                     76: }
                     77:
                     78: void homfctr(vl,g,dcp)
                     79: VL vl;
                     80: P g;
                     81: DCP *dcp;
                     82: {
                     83:        P s,t,u,f,h,x;
                     84:        Q e;
                     85:        int d,d0;
                     86:        DCP dc,dct;
                     87:
                     88:        substp(vl,g,vl->v,(P)ONE,&t); fctrp(vl,t,&dc); MKV(vl->v,x);
                     89:        for ( dct = dc; dct; dct = NEXT(dct) )
                     90:                if ( !NUM(dc->c) ) {
                     91:                        for ( s = 0, f = dc->c, d = d0 = homdeg(f); f; d = homdeg(f) ) {
                     92:                                exthp(vl,f,d,&h); STOQ(d0-d,e); pwrp(vl,x,e,&t);
                     93:                                mulp(vl,t,h,&u); addp(vl,s,u,&t); s = t;
                     94:                                subp(vl,f,h,&u); f = u;
                     95:                        }
                     96:                        dc->c = s;
                     97:                }
                     98:        *dcp = dc;
                     99: }
                    100:
                    101: void mfctr(vl,f,dcp)
                    102: VL vl;
                    103: P f;
                    104: DCP *dcp;
                    105: {
                    106:        DCP dc,dc0,dct,dcs,dcr;
                    107:        P p,pmin,ppmin,cmin,t;
                    108:        VL mvl;
                    109:        Q c;
                    110:
                    111:        ptozp(f,1,&c,&p);
                    112:        NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
                    113:        msqfr(vl,p,&dct);
                    114:        for ( ; dct; dct = NEXT(dct) ) {
                    115:                mindegp(vl,COEF(dct),&mvl,&pmin);
                    116: #if 0
                    117:                minlcdegp(vl,COEF(dct),&mvl,&pmin);
                    118:                min_common_vars_in_coefp(vl,COEF(dct),&mvl,&pmin);
                    119: #endif
                    120:                pcp(mvl,pmin,&ppmin,&cmin);
                    121:                if ( !NUM(cmin) ) {
                    122:                        mfctrmain(mvl,cmin,&dcs);
                    123:                        for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
                    124:                                DEG(dcr) = DEG(dct);
                    125:                                reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
                    126:                        }
                    127:                        for ( ; NEXT(dc); dc = NEXT(dc) );
                    128:                        NEXT(dc) = dcs;
                    129:                }
                    130:                mfctrmain(mvl,ppmin,&dcs);
                    131:                for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
                    132:                        DEG(dcr) = DEG(dct);
                    133:                        reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
                    134:                }
                    135:                for ( ; NEXT(dc); dc = NEXT(dc) );
                    136:                NEXT(dc) = dcs;
                    137:        }
                    138:        adjsgn(f,dc0); *dcp = dc0;
                    139: }
                    140:
1.2       noro      141: #if 0
1.1       noro      142: void adjsgn(p,dc)
                    143: P p;
                    144: DCP dc;
                    145: {
                    146:        int sgn;
                    147:        DCP dct;
                    148:        P c;
                    149:
                    150:        for ( dct = dc, sgn = headsgn(p); dct; dct = NEXT(dct) )
                    151:                if ( !EVENN(NM(dct->d)) )
                    152:                        sgn *= headsgn(COEF(dct));
                    153:        if ( sgn < 0 ) {
                    154:                chsgnp(COEF(dc),&c); COEF(dc) = c;
                    155:        }
                    156: }
1.2       noro      157: #else
                    158: void adjsgn(p,dc)
                    159: P p;
                    160: DCP dc;
                    161: {
                    162:        int sgn;
                    163:        DCP dct;
                    164:        P c;
                    165:
                    166:        if ( headsgn(COEF(dc)) != headsgn(p) ) {
                    167:                chsgnp(COEF(dc),&c); COEF(dc) = c;
                    168:        }
                    169:        for ( dct = NEXT(dc); dct; dct = NEXT(dct) )
                    170:                if ( headsgn(COEF(dct)) < 0 ) {
                    171:                        chsgnp(COEF(dct),&c); COEF(dct) = c;
                    172:                }
                    173: }
                    174: #endif
1.1       noro      175:
                    176: int headsgn(p)
                    177: P p;
                    178: {
                    179:        if ( !p )
                    180:                return 0;
                    181:        else {
                    182:                while ( !NUM(p) )
                    183:                        p = COEF(DC(p));
                    184:                return SGN((Q)p);
                    185:        }
                    186: }
                    187:
                    188: void fctrwithmvp(vl,f,v,dcp)
                    189: VL vl;
                    190: P f;
                    191: V v;
                    192: DCP *dcp;
                    193: {
                    194:        VL nvl;
                    195:        DCP dc;
                    196:
                    197:        if ( NUM(f) ) {
                    198:                NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
                    199:                NEXT(dc) = 0; *dcp = dc;
                    200:                return;
                    201:        }
                    202:
                    203:        clctv(vl,f,&nvl);
                    204:        if ( !NEXT(nvl) )
                    205:                ufctr(f,1,dcp);
                    206:        else
                    207:                mfctrwithmv(nvl,f,v,dcp);
                    208: }
                    209:
                    210: void mfctrwithmv(vl,f,v,dcp)
                    211: VL vl;
                    212: P f;
                    213: V v;
                    214: DCP *dcp;
                    215: {
                    216:        DCP dc,dc0,dct,dcs,dcr;
                    217:        P p,pmin,ppmin,cmin,t;
                    218:        VL mvl;
                    219:        Q c;
                    220:
                    221:        ptozp(f,1,&c,&p);
                    222:        NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
                    223:        msqfr(vl,p,&dct);
                    224:        for ( ; dct; dct = NEXT(dct) ) {
                    225:                if ( !v )
                    226:                        mindegp(vl,COEF(dct),&mvl,&pmin);
                    227:                else {
                    228:                        reordvar(vl,v,&mvl); reorderp(mvl,vl,COEF(dct),&pmin);
                    229:                }
                    230:                pcp(mvl,pmin,&ppmin,&cmin);
                    231:                if ( !NUM(cmin) ) {
                    232:                        mfctrmain(mvl,cmin,&dcs);
                    233:                        for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
                    234:                                DEG(dcr) = DEG(dct);
                    235:                                reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
                    236:                        }
                    237:                        for ( ; NEXT(dc); dc = NEXT(dc) );
                    238:                        NEXT(dc) = dcs;
                    239:                }
                    240:                mfctrmain(mvl,ppmin,&dcs);
                    241:                for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
                    242:                        DEG(dcr) = DEG(dct);
                    243:                        reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
                    244:                }
                    245:                for ( ; NEXT(dc); dc = NEXT(dc) );
                    246:                NEXT(dc) = dcs;
                    247:        }
                    248:        *dcp = dc0;
                    249: }
                    250:
                    251: void ufctr(f,hint,dcp)
                    252: P f;
                    253: int hint;
                    254: DCP *dcp;
                    255: {
                    256:        P p,c;
                    257:        DCP dc,dct,dcs,dcr;
                    258:
                    259:        ptozp(f,SGN((Q)UCOEF(f)),(Q *)&c,&p);
                    260:        NEWDC(dc); COEF(dc) = c; DEG(dc) = ONE; NEXT(dc) = 0;
                    261:        usqp(p,&dct);
                    262:        for ( *dcp = dc; dct; dct = NEXT(dct) ) {
                    263:                ufctrmain(COEF(dct),hint,&dcs);
                    264:                for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
                    265:                        DEG(dcr) = DEG(dct);
                    266:                for ( ; NEXT(dc); dc = NEXT(dc) );
                    267:                NEXT(dc) = dcs;
                    268:        }
                    269: }
                    270:
                    271: void mfctrmain(vl,p,dcp)
                    272: VL vl;
                    273: P p;
                    274: DCP *dcp;
                    275: {
1.5       noro      276:        int i,j,k,*win,np,x;
1.1       noro      277:        VL nvl,tvl;
                    278:        VN vn,vnt,vn1;
                    279:        P p0,f,g,f0,g0,s,t,lp,m;
                    280:        P *fp0,*fpt,*l,*tl;
                    281:        DCP dc,dc0,dcl;
                    282:        int count,nv;
1.5       noro      283:        int *nonzero;
1.1       noro      284:
                    285:        if ( !cmpq(DEG(DC(p)),ONE) ) {
                    286:                NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                    287:                *dcp = dc; return;
                    288:        }
                    289:        lp = LC(p); fctrp(vl,lp,&dcl); clctv(vl,p,&nvl);
                    290:        for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
                    291:        W_CALLOC(nv,struct oVN,vn); W_CALLOC(nv,struct oVN,vnt);
                    292:        W_CALLOC(nv,struct oVN,vn1);
1.5       noro      293:        W_CALLOC(nv,int,nonzero);
                    294:
1.1       noro      295:        for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
                    296:                vn1[i].v = vn[i].v = tvl->v;
                    297:        vn1[i].v = vn[i].v = 0;
1.5       noro      298:
                    299:        /* determine a heuristic bound of deg(GCD(p,p')) */
                    300:        while ( 1 ) {
                    301:                for ( i = 0; vn1[i].v; i++ )
                    302:                        vn1[i].n = ((unsigned int)random())%256+1;
                    303:                substvp(nvl,LC(p),vn1,&p0);
                    304:                if ( p0 ) {
                    305:                        substvp(nvl,p,vn1,&p0);
                    306:                        if ( sqfrchk(p0) ) {
                    307:                                ufctr(p0,1,&dc0);
                    308:                                if ( NEXT(NEXT(dc0)) == 0 ) {
                    309:                                        NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                    310:                                        *dcp = dc;
                    311:                                        return;
                    312:                                } else {
                    313:                                        for ( dc0 = NEXT(dc0), np = 0; dc0; dc0 = NEXT(dc0), np++ );
                    314:                                        break;
                    315:                                }
                    316:                        }
                    317:                }
                    318:        }
                    319:
1.6       noro      320:        /* determine the position of variables which is not allowed to
                    321:           be set to 0 */
1.5       noro      322:        for ( i = 0; vn1[i].v; i++ ) {
                    323:                x = vn1[i].n; vn1[i].n = 0;
                    324:                substvp(nvl,LC(p),vn1,&p0);
                    325:                if ( !p0 )
                    326:                        vn1[i].n = x;
                    327:                else {
                    328:                        substvp(nvl,p,vn1,&p0);
                    329:                        if ( !sqfrchk(p0) )
                    330:                                vn1[i].n = x;
                    331:                        else {
                    332:                                ufctr(p0,1,&dc0);
                    333:                                for ( dc0 = NEXT(dc0), j = 0; dc0; dc0 = NEXT(dc0), j++ );
                    334:                                if ( j > np )
                    335:                                        vn1[i].n = x;
                    336:                        }
                    337:                }
                    338:        }
                    339:        for ( i = 0; vn1[i].v; i++ )
                    340:                if (vn1[i].n )
                    341:                        nonzero[i] = 1;
                    342:
1.1       noro      343:        count = 0;
                    344:        while  ( 1 ) {
                    345:                while ( 1 ) {
                    346:                        count++;
                    347:                        for ( i = 0, j = 0; vn[i].v; i++ )
                    348:                                if ( vn[i].n )
                    349:                                        vnt[j++].v = (V)i;
                    350:                        vnt[j].n = 0; mulsgn(vn,vnt,j,vn1);
                    351:                        for ( i = 0; vn1[i].v; i++ )
                    352:                                if ( vn1[i].n ) {
                    353:                                        if ( vn1[i].n > 0 )
                    354:                                                vn1[i].n = sprime[vn1[i].n];
                    355:                                        else
                    356:                                                vn1[i].n = -sprime[-vn1[i].n];
                    357:                                }
                    358:                        if ( valideval(nvl,dcl,vn1) ) {
                    359:                                substvp(nvl,p,vn1,&p0);
                    360:                                if ( sqfrchk(p0) ) {
                    361:                                        ufctr(p0,1,&dc0);
                    362:                                        if ( NEXT(NEXT(dc0)) == 0 ) {
                    363:                                                NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                    364:                                                *dcp = dc;
                    365:                                                return;
1.5       noro      366:                                        } else {
                    367:                                                for ( dc = NEXT(dc0), i = 0; dc; dc = NEXT(dc), i++ );
                    368:                                                if ( i <= np )
                    369:                                                        goto MAIN;
                    370:                                                if ( i < np )
                    371:                                                        np = i;
                    372:                                        }
                    373:                                }
1.1       noro      374:                        }
                    375:                        if ( nextbin(vnt,j) )
                    376:                                break;
                    377:                }
1.5       noro      378:                while ( 1 ) {
                    379:                        next(vn);
                    380:                        for ( i = 0; vn[i].v; i++ )
                    381:                                if ( nonzero[i] && !vn[i].n )
                    382:                                        break;
                    383:                        if ( !vn[i].v )
                    384:                                break;
                    385:                }
1.1       noro      386:        }
                    387: MAIN :
                    388: #if 0
                    389:        for ( i = 0; vn1[i].v; i++ )
                    390:                fprintf(stderr,"%d ",vn1[i].n);
                    391:        fprintf(stderr,"\n");
                    392: #endif
                    393:        for ( np = 0, dc = NEXT(dc0); dc; dc = NEXT(dc), np++ );
                    394:        fp0 = (P *)ALLOCA((np + 1)*sizeof(P));
                    395:        fpt = (P *)ALLOCA((np + 1)*sizeof(P));
                    396:        l = tl = (P *)ALLOCA((np + 1)*sizeof(P));
                    397:        for ( i = 1, dc = NEXT(dc0); dc; i++, dc = NEXT(dc) )
                    398:                fp0[i-1] = COEF(dc);
                    399: #if 0
                    400:        sort_by_deg(np,fp0,fpt);
                    401:        sort_by_deg_rev(np-1,fpt+1,fp0+1);
                    402: #endif
                    403: #if 0
                    404:        sort_by_deg_rev(np,fp0,fpt); fp0[0] = fpt[0];
                    405:        sort_by_deg(np-1,fpt+1,fp0+1); fp0[np] = 0;
                    406: #endif
                    407:        fp0[np] = 0;
                    408:        win = W_ALLOC(np + 1); f = p; divsp(vl,p0,COEF(dc0),&f0);
                    409:        for ( k = 1, win[0] = 1, --np; ; ) {
                    410:                P h0,lcg,lch;
                    411:                Q c;
                    412:
                    413:                g0 = fp0[win[0]];
                    414:                for ( i = 1; i < k; i++ ) {
                    415:                        mulp(vl,g0,fp0[win[i]],&m); g0 = m;
                    416:                }
                    417:                mulq((Q)LC(g0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lcg);
                    418:                divsp(nvl,f0,g0,&h0);
                    419:                mulq((Q)LC(h0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lch);
                    420:                mfctrhen2(nvl,vn1,f,f0,g0,h0,lcg,lch,&g);
                    421:                if ( g ) {
                    422:                        *tl++ = g; divsp(vl,f,g,&t);
                    423:                        f = t; divsp(vl,f0,g0,&t); ptozp(t,1,(Q *)&s,&f0);
                    424:                        for ( i = 0; i < k - 1; i++ )
                    425:                                for ( j = win[i] + 1; j < win[i + 1]; j++ )
                    426:                                        fp0[j - i - 1] = fp0[j];
                    427:                        for ( j = win[k - 1] + 1; j <= np; j++ )
                    428:                                        fp0[j - k] = fp0[j];
                    429:                        if ( ( np -= k ) < k ) break;
                    430:                        if ( np - win[0] + 1 < k )
                    431:                                if ( ++k <= np ) {
                    432:                                        for ( i = 0; i < k; i++ ) win[i] = i + 1;
                    433:                                        continue;
                    434:                                } else
                    435:                                        break;
                    436:                        else for ( i = 1; i < k; i++ ) win[i] = win[0] + i;
                    437:                } else {
                    438:                        if ( ncombi(1,np,k,win) == 0 )
                    439:                                if ( k == np ) break;
                    440:                                else for ( i = 0, ++k; i < k; i++ ) win[i] = i + 1;
                    441:                }
                    442:        }
                    443:        *tl++ = f; *tl = 0;
                    444:        for ( dc0 = 0; *l; l++ ) {
                    445:                NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = *l;
                    446:        }
                    447:        NEXT(dc) = 0; *dcp = dc0;
                    448: }
                    449:
                    450: void ufctrmain(p,hint,dcp)
                    451: P p;
                    452: int hint;
                    453: DCP *dcp;
                    454: {
                    455:        ML list;
                    456:        DCP dc;
                    457:
                    458:        if ( NUM(p) || (UDEG(p) == 1) ) {
                    459:                NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                    460:                *dcp = dc;
                    461:        } else if ( iscycm(p) )
                    462:                cycm(VR(p),UDEG(p),dcp);
                    463:        else if ( iscycp(p) )
                    464:                cycp(VR(p),UDEG(p),dcp);
                    465:        else {
                    466:                hensel(5,5,p,&list);
                    467:                if ( list->n == 1 ) {
                    468:                        NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                    469:                        *dcp = dc;
                    470:                } else
                    471:                        dtest(p,list,hint,dcp);
                    472:        }
                    473: }
                    474:
                    475: struct oMF {
                    476:        int m;
                    477:        P f;
                    478: };
                    479:
                    480: void calcphi();
                    481:
                    482: void cycm(v,n,dcp)
                    483: V v;
                    484: register int n;
                    485: DCP *dcp;
                    486: {
                    487:        register int i,j;
                    488:        struct oMF *mfp;
                    489:        DCP dc,dc0;
                    490:
                    491:        if ( n == 1 ) {
                    492:                NEWDC(dc); mkssum(v,1,1,-1,&COEF(dc)); DEG(dc) = ONE;
                    493:        } else {
                    494:                mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
                    495:                bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
                    496:                for ( i = 1, j = 0; i <= n; i++ )
                    497:                        if ( !(n%i) ) mfp[j++].m = i;
                    498:                mkssum(v,1,1,-1,&mfp[0].f);
                    499:                for ( i = 1; i < j; i++ )
                    500:                        calcphi(v,i,mfp);
                    501:                for ( dc0 = 0, i = 0; i < j; i++ ) {
                    502:                        NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
                    503:                }
                    504:        }
                    505:        NEXT(dc) = 0; *dcp = dc0;
                    506: }
                    507:
                    508: void cycp(v,n,dcp)
                    509: V v;
                    510: register int n;
                    511: DCP *dcp;
                    512: {
                    513:        register int i,j;
                    514:        int n0;
                    515:        struct oMF *mfp;
                    516:        DCP dc,dc0;
                    517:
                    518:        if ( n == 1 ) {
                    519:                NEWDC(dc); mkssum(v,1,1,1,&COEF(dc)); DEG(dc) = ONE;
                    520:        } else {
                    521:                n0 = n; n *= 2;
                    522:                mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
                    523:                bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
                    524:                for ( i = 1, j = 0; i <= n; i++ )
                    525:                        if ( !(n%i) ) mfp[j++].m = i;
                    526:                mkssum(v,1,1,-1,&mfp[0].f);
                    527:                for ( i = 1; i < j; i++ )
                    528:                        calcphi(v,i,mfp);
                    529:                for ( dc0 = 0, i = 0; i < j; i++ )
                    530:                        if ( n0 % mfp[i].m ) {
                    531:                                NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
                    532:                        }
                    533:        }
                    534:        NEXT(dc) = 0; *dcp = dc0;
                    535: }
                    536:
                    537: void calcphi(v,n,mfp)
                    538: V v;
                    539: int n;
                    540: register struct oMF *mfp;
                    541: {
                    542:        register int i,m;
                    543:        P t,s,tmp;
                    544:
                    545:        for ( t = (P)ONE, i = 0, m = mfp[n].m; i < n; i++ )
                    546:                if ( !(m%mfp[i].m) ) {
                    547:                        mulp(CO,t,mfp[i].f,&s); t = s;
                    548:                }
                    549:        mkssum(v,m,1,-1,&s); udivpz(s,t,&mfp[n].f,&tmp);
                    550:        if ( tmp )
                    551:                error("calcphi: cannot happen");
                    552: }
                    553:
                    554: void mkssum(v,e,s,sgn,r)
                    555: V v;
                    556: int e,s,sgn;
                    557: P *r;
                    558: {
                    559:        register int i,sgnt;
                    560:        DCP dc,dc0;
                    561:        Q q;
                    562:
                    563:        for ( dc0 = 0, i = s, sgnt = 1; i >= 0; i--, sgnt *= sgn ) {
                    564:                if ( !dc0 ) {
                    565:                        NEWDC(dc0); dc = dc0;
                    566:                } else {
                    567:                        NEWDC(NEXT(dc)); dc = NEXT(dc);
                    568:                }
                    569:                STOQ(i*e,DEG(dc)); STOQ(sgnt,q),COEF(dc) = (P)q;
                    570:        }
                    571:        NEXT(dc) = 0; MKP(v,dc0,*r);
                    572: }
                    573:
                    574: int iscycp(f)
                    575: P f;
                    576: {
                    577:        DCP dc;
                    578:        dc = DC(f);
                    579:
                    580:        if ( !UNIQ((Q)COEF(dc)) )
                    581:                return ( 0 );
                    582:        dc = NEXT(dc);
                    583:        if ( NEXT(dc) || DEG(dc) || !UNIQ((Q)COEF(dc)) )
                    584:                return ( 0 );
                    585:        return ( 1 );
                    586: }
                    587:
                    588: int iscycm(f)
                    589: P f;
                    590: {
                    591:        DCP dc;
                    592:
                    593:        dc = DC(f);
                    594:        if ( !UNIQ((Q)COEF(dc)) )
                    595:                return ( 0 );
                    596:        dc = NEXT(dc);
                    597:        if ( NEXT(dc) || DEG(dc) || !MUNIQ((Q)COEF(dc)) )
                    598:                return ( 0 );
                    599:        return ( 1 );
                    600: }
                    601:
                    602: void sortfs(dcp)
                    603: DCP *dcp;
                    604: {
                    605:        int i,k,n,k0,d;
                    606:        DCP dc,dct,t;
                    607:        DCP *a;
                    608:
                    609:        dc = *dcp;
                    610:        for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
                    611:        a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
                    612:        for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ )
                    613:                a[i] = dct;
                    614:        a[n] = 0;
                    615:
                    616:        for ( i = 0; i < n; i++ ) {
                    617:                for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ )
                    618:                        if ( (int)UDEG(COEF(a[k])) < d ) {
                    619:                                k0 = k;
                    620:                                d = UDEG(COEF(a[k]));
                    621:                        }
                    622:                if ( i != k0 ) {
                    623:                        t = a[i]; a[i] = a[k0]; a[k0] = t;
                    624:                }
                    625:        }
                    626:        for ( *dcp = dct = a[0], i = 1; i < n; i++ )
                    627:                dct = NEXT(dct) = a[i];
                    628:        NEXT(dct) = 0;
                    629: }
                    630:
                    631: void sortfsrev(dcp)
                    632: DCP *dcp;
                    633: {
                    634:        int i,k,n,k0,d;
                    635:        DCP dc,dct,t;
                    636:        DCP *a;
                    637:
                    638:        dc = *dcp;
                    639:        for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
                    640:        a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
                    641:        for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ )
                    642:                a[i] = dct;
                    643:        a[n] = 0;
                    644:
                    645:        for ( i = 0; i < n; i++ ) {
                    646:                for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ )
                    647:                        if ( (int)UDEG(COEF(a[k])) > d ) {
                    648:                                k0 = k;
                    649:                                d = UDEG(COEF(a[k]));
                    650:                        }
                    651:                if ( i != k0 ) {
                    652:                        t = a[i]; a[i] = a[k0]; a[k0] = t;
                    653:                }
                    654:        }
                    655:        for ( *dcp = dct = a[0], i = 1; i < n; i++ )
                    656:                dct = NEXT(dct) = a[i];
                    657:        NEXT(dct) = 0;
                    658: }
                    659:
                    660: void nthrootchk(f,dc,fp,dcp)
                    661: P f;
                    662: struct oDUM *dc;
                    663: ML fp;
                    664: DCP *dcp;
                    665: {
                    666:        register int i,k;
                    667:        int e,n,dr,tmp,t;
                    668:        int *tmpp,**tmppp;
                    669:        int **pp,**wpp;
                    670:        LUM fpa,tpa,f0l;
                    671:        UM wt,wq,ws,dvr,f0,fe;
                    672:        N lc;
                    673:        int lcm;
                    674:        int m,b;
                    675:
                    676:        m = fp->mod; b = fp->bound; fe = *((UM *)fp->c);
                    677:        e = dc->n; f0 = dc->f; nthrootn(NM((Q)COEF(DC(f))),e,&lc);
                    678:        if ( !lc ) {
                    679:                *dcp = 0;
                    680:                return;
                    681:        }
                    682:        lcm = rem(lc,m); W_LUMALLOC(DEG(f0),b,f0l);
                    683:        for ( i = DEG(f0), tmppp = COEF(f0l), tmpp = COEF(f0);
                    684:                i >= 0; i-- )
                    685:                *(tmppp[i]) = dmb(m,tmpp[i],lcm,&tmp);
                    686:        dtestroot(m,1,f,f0l,dc,dcp);
                    687:        if ( *dcp )
                    688:                return;
                    689:        n = UDEG(f); W_LUMALLOC(n,b,fpa); W_LUMALLOC(n,b,tpa);
                    690:        ptolum(m,b,f,fpa);
                    691:        dvr = W_UMALLOC(n); wq = W_UMALLOC(n);
                    692:        wt = W_UMALLOC(n); ws = W_UMALLOC(n);
                    693:        cpyum(fe,dvr); divum(m,dvr,f0,wq);
                    694:        t = dmb(m,pwrm(m,lcm,e-1),e,&tmp); DEG(dvr) = DEG(wq);
                    695:        for ( k = 0; k <= DEG(wq); k++ )
                    696:                COEF(dvr)[k] = dmb(m,COEF(wq)[k],t,&tmp);
                    697:        for ( i = 1; i < b; i++ ) {
                    698:                pwrlum(m,i+1,f0l,e,tpa);
                    699:                for ( k = 0, pp = COEF(fpa), wpp = COEF(tpa);
                    700:                        k <= n; k++ )
                    701:                        COEF(wt)[k] = (pp[k][i]-wpp[k][i]+m)%m;
                    702:                degum(wt,n); dr = divum(m,wt,dvr,ws);
                    703:                if ( dr >= 0 ) {
                    704:                        *dcp = 0;
                    705:                        return;
                    706:                }
                    707:                for ( k = 0, pp = COEF(f0l); k <= DEG(ws); k++ )
                    708:                        pp[k][i] = COEF(ws)[k];
                    709:                dtestroot(m,i+1,f,f0l,dc,dcp);
                    710:                if ( *dcp )
                    711:                        return;
                    712:        }
                    713: }
                    714:
                    715: void sqfrp(vl,f,dcp)
                    716: VL vl;
                    717: P f;
                    718: DCP *dcp;
                    719: {
                    720:        P c,p;
                    721:        DCP dc,dc0;
                    722:
                    723:        if ( !f || NUM(f) ) {
                    724:                NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
                    725:                COEF(dc0) = f; NEXT(dc0) = 0;
                    726:        } else if ( !qpcheck((Obj)f) )
                    727:                error("sqfrp : invalid argument");
                    728:        else {
                    729:                NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
                    730:                ptozp(f,1,(Q *)&c,&p); msqfr(vl,p,&dc);
                    731:                COEF(dc0) = c; NEXT(dc0) = dc;
                    732:                adjsgn(f,dc0);
                    733:        }
                    734: }
                    735:
                    736: /*
                    737:  * f : must be a poly with int coef, ignore int content
                    738:  */
                    739: void msqfr(vl,f,dcp)
                    740: VL vl;
                    741: P f;
                    742: DCP *dcp;
                    743: {
                    744:        DCP dc,dct,dcs;
                    745:        P c,p,t,s,pc;
                    746:        VL mvl;
                    747:
                    748:        ptozp(f,1,(Q *)&c,&t); monomialfctr(vl,t,&p,&dc);
                    749:        if ( NUM(p) ) {
                    750:                *dcp = dc;
                    751:                return;
                    752:        }
                    753:        mindegp(vl,p,&mvl,&s);
                    754: #if 0
                    755:        minlcdegp(vl,p,&mvl,&s);
                    756:        min_common_vars_in_coefp(vl,p,&mvl,&s);
                    757: #endif
                    758:        pcp(mvl,s,&p,&pc);
                    759:        if ( !NUM(pc) ) {
                    760:                msqfr(mvl,pc,&dcs);
                    761:                if ( !dc )
                    762:                        dc = dcs;
                    763:                else {
                    764:                        for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
                    765:                        NEXT(dct) = dcs;
                    766:                }
                    767:        }
                    768:        msqfrmain(mvl,p,&dcs);
                    769:        for ( dct = dcs; dct; dct = NEXT(dct) ) {
                    770:                reorderp(vl,mvl,COEF(dct),&t); COEF(dct) = t;
                    771:        }
                    772:        if ( !dc )
                    773:                *dcp = dcs;
                    774:        else {
                    775:                for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
                    776:                NEXT(dct) = dcs; *dcp = dc;
                    777:        }
                    778: }
                    779:
                    780: void usqp(f,dcp)
                    781: P f;
                    782: DCP *dcp;
                    783: {
                    784:        int index,nindex;
                    785:        P g,c,h;
                    786:        DCP dc;
                    787:
                    788:        ptozp(f,1,(Q *)&c,&h);
                    789:        if ( SGN((Q)LC(h)) < 0 )
                    790:                chsgnp(h,&g);
                    791:        else
                    792:                g = h;
                    793:        for ( index = 0, dc = 0; !dc; index = nindex )
                    794:                hsq(index,5,g,&nindex,&dc);
                    795:        *dcp = dc;
                    796: }
                    797:
                    798: void msqfrmain(vl,p,dcp)
                    799: VL vl;
                    800: P p;
                    801: DCP *dcp;
                    802: {
                    803:        int i,j;
                    804:        VL nvl,tvl;
                    805:        VN vn,vnt,vn1;
                    806:        P gg,tmp,p0,g;
                    807:        DCP dc,dc0,dcr,dcr0;
                    808:        int nv,d,d1;
                    809:        int found;
                    810:        VN svn1;
                    811:        P sp0;
                    812:        DCP sdc0;
                    813:
                    814:        if ( deg(VR(p),p) == 1 ) {
                    815:                NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                    816:                *dcp = dc;
                    817:                return;
                    818:        }
                    819:        clctv(vl,p,&nvl);
                    820:        for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
                    821:        if ( nv == 1 ) {
                    822:                usqp(p,dcp);
                    823:                return;
                    824:        }
                    825: #if 0
                    826:        if ( heusqfrchk(nvl,p) ) {
                    827:                NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                    828:                *dcp = dc;
                    829:                return;
                    830:        }
                    831: #endif
                    832:        W_CALLOC(nv,struct oVN,vn);
                    833:        W_CALLOC(nv,struct oVN,vnt);
                    834:        W_CALLOC(nv,struct oVN,vn1);
                    835:        W_CALLOC(nv,struct oVN,svn1);
                    836:        for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
                    837:                vn1[i].v = vn[i].v = tvl->v;
                    838:        vn1[i].v = vn[i].v = 0;
1.5       noro      839:
                    840:        /* determine a heuristic bound of deg(GCD(p,p')) */
                    841:        while ( 1 ) {
                    842:                for ( i = 0; vn1[i].v; i++ )
                    843:                        vn1[i].n = ((unsigned int)random())%256+1;
                    844:                substvp(nvl,LC(p),vn1,&tmp);
                    845:                if ( tmp ) {
                    846:                        substvp(nvl,p,vn1,&p0);
                    847:                        usqp(p0,&dc0);
                    848:                        for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )
                    849:                                if ( DEG(dc) )
                    850:                                        d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));
                    851:                        if ( d1 == 0 ) {
                    852:                                /* p is squarefree */
                    853:                                NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                    854:                                *dcp = dc;
                    855:                                return;
                    856:                        } else {
1.7     ! noro      857:                                d = d1+1; /* XXX : try searching better evaluation */
        !           858:                                found = 0;
1.5       noro      859:                                break;
                    860:                        }
                    861:                }
                    862:        }
1.1       noro      863:
1.6       noro      864:        for  ( dcr0 = 0, g = p; ; ) {
1.1       noro      865:                while ( 1 ) {
                    866:                        for ( i = 0, j = 0; vn[i].v; i++ )
                    867:                                if ( vn[i].n ) vnt[j++].v = (V)i;
                    868:                        vnt[j].n = 0;
                    869:
                    870:                        mulsgn(vn,vnt,j,vn1);
                    871:                        substvp(nvl,LC(g),vn1,&tmp);
                    872:                        if ( tmp ) {
                    873:                                substvp(nvl,g,vn1,&p0);
                    874:                                usqp(p0,&dc0);
                    875:                                for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )
                    876:                                        if ( DEG(dc) )
                    877:                                                d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));
                    878:
                    879:                                if ( d1 == 0 ) {
                    880:                                        NEWDC(dc); DEG(dc) = ONE; COEF(dc) = g; NEXT(dc) = 0;
                    881:                                        if ( !dcr0 )
                    882:                                                dcr0 = dc;
                    883:                                        else
                    884:                                                NEXT(dcr) = dc;
                    885:                                        *dcp = dcr0;
                    886:                                        return;
                    887:                                }
                    888:
                    889:                                if ( d < d1 )
                    890:                                        goto END;
                    891:                                if ( d > d1 ) {
                    892:                                        d = d1;
                    893:                                        found = 1;
                    894:                                        bcopy((char *)vn1,(char *)svn1,(int)(sizeof(struct oVN)*nv));
                    895:                                        sp0 = p0; sdc0 = dc0;
                    896:                                        goto END;
                    897:                                }
                    898:                                /* d1 == d */
                    899:                                if ( found ) {
                    900:                                        found = 0;
                    901: #if 0
                    902:                                if ( d > d1 ) {
                    903:                                        d = d1;
                    904:                                                        /*} */
                    905: #endif
                    906:                                        msqfrmainmain(nvl,g,svn1,sp0,sdc0,&dc,&gg); g = gg;
                    907:                                        if ( dc ) {
                    908:                                                if ( !dcr0 )
                    909:                                                        dcr0 = dc;
                    910:                                                else
                    911:                                                        NEXT(dcr) = dc;
                    912:                                                for ( dcr = dc; NEXT(dcr); dcr = NEXT(dcr) );
                    913:                                                if ( NUM(g) ) {
                    914:                                                        NEXT(dcr) = 0; *dcp = dcr0;
                    915:                                                        return;
                    916:                                                }
                    917:                                                d = deg(VR(g),g);
                    918:                                        }
                    919:                                }
                    920:                        }
                    921: END:
                    922:                        if ( nextbin(vnt,j) )
                    923:                                break;
                    924:                }
                    925:                next(vn);
                    926:        }
                    927: }
                    928:
                    929: void msqfrmainmain(vl,p,vn,p0,dc0,dcp,pp)
                    930: VL vl;
                    931: P p;
                    932: VN vn;
                    933: P p0;
                    934: DCP dc0;
                    935: DCP *dcp;
                    936: P *pp;
                    937: {
                    938:        int i,j,k,np;
                    939:        DCP *a;
                    940:        DCP dc,dcr,dcr0,dct;
                    941:        P g,t,s,u,t0,f,f0,d,d0,g0,h0,x,xx;
                    942:        Q q;
                    943:        V v;
                    944:
                    945:        for ( np = 0, dc = dc0; dc; dc = NEXT(dc), np++ );
                    946:        a = (DCP *)ALLOCA((np + 1)*sizeof(DCP));
                    947:        for ( i = 0, dc = dc0; dc; i++, dc = NEXT(dc) )
                    948:                a[np-i-1] = dc;
                    949:
                    950:        for ( i = 0, dcr0 = 0, f = p, f0 = p0, v = VR(p);
                    951:                i < np; i++ ) {
                    952:                if ( (i == (np-1))&&UNIQ(DEG(a[i])) ) {
                    953:                        NEXTDC(dcr0,dcr);
                    954:                        DEG(dcr) = DEG(a[i]);
                    955:                        COEF(dcr) = f;
                    956:                        f = (P)ONE;
                    957:                } else if ( (i == (np-1))&&UNIQ(DEG(DC(COEF(a[i])))) ) {
                    958:                        diffp(vl,f,v,&s); pcp(vl,s,&t,&u);
                    959:                        if ( divtpz(vl,f,t,&s) ) {
                    960:                                NEXTDC(dcr0,dcr);
                    961:                                DEG(dcr) = DEG(a[i]);
                    962:                                COEF(dcr) = s;
                    963:                                f = (P)ONE;
                    964:                        } else
                    965:                                break;
                    966:                } else {
                    967:                        for ( t = f, t0 = f0,
                    968:                                j = 0, k = QTOS(DEG(a[i]))-1; j < k; j++ ) {
                    969:                                diffp(vl,t,v,&s); t = s;
                    970:                                diffp(vl,t0,v,&s); t0 = s;
                    971:                        }
                    972:                        factorial(k,&q);
                    973:                        divsp(vl,t,(P)q,&s);
                    974:                        for ( dct = DC(s); NEXT(dct); dct = NEXT(dct) );
                    975:                        if ( DEG(dct) ) {
                    976:                                MKV(VR(s),x); pwrp(vl,x,DEG(dct),&xx);
                    977:                                divsp(vl,s,xx,&d);
                    978:                        } else {
                    979:                                xx = (P)ONE; d = s;
                    980:                        }
                    981:                        divsp(vl,t0,xx,&t);
                    982:                        divsp(vl,t,(P)q,&s);
                    983:                        ptozp(s,1,(Q *)&t,&d0);
                    984:
                    985:                        for ( dct = DC(COEF(a[i])); NEXT(dct); dct = NEXT(dct) );
                    986:                        if ( DEG(dct) )
                    987:                                divsp(vl,COEF(a[i]),xx,&g0);
                    988:                        else {
                    989:                                xx = (P)ONE; g0 = COEF(a[i]);
                    990:                        }
                    991:
                    992:                        pcp(vl,d,&t,&u); d = t;
                    993:                        ptozp(g0,1,(Q *)&u,&t); g0 = t;
                    994:
                    995:                {
                    996:                        DCP dca,dcb;
                    997:
                    998:                        fctrp(vl,LC(d),&dca);
                    999:                        for ( dcb = dca, u = (P)ONE; dcb; dcb = NEXT(dcb) ) {
                   1000:                                if ( NUM(COEF(dcb)) ) {
                   1001:                                        mulp(vl,u,COEF(dcb),&t); u = t;
                   1002:                                } else {
                   1003:                                        Q qq;
                   1004:                                        P tt;
                   1005:
                   1006:                                        pwrp(vl,COEF(dcb),DEG(a[i]),&s);
                   1007:                                        for ( t = LC(f), j = 0; divtpz(vl,t,s,&tt); j++, t = tt );
                   1008:                                        STOQ(j,qq);
                   1009:                                        if ( cmpq(qq,DEG(dcb)) > 0 )
                   1010:                                                qq = DEG(dcb);
                   1011:                                        pwrp(vl,COEF(dcb),qq,&t); mulp(vl,u,t,&s); u = s;
                   1012:                                }
                   1013:                        }
                   1014:                        divsp(vl,d0,g0,&h0);
                   1015:                }
                   1016:                        mfctrhen2(vl,vn,d,d0,g0,h0,u,LC(d),&s);
                   1017:                        if ( s ) {
                   1018:                                mulp(vl,s,xx,&g);
                   1019:                                pwrp(vl,g,DEG(a[i]),&t);
                   1020:                                if ( divtpz(vl,f,t,&s) ) {
                   1021:                                        NEXTDC(dcr0,dcr);
                   1022:                                        DEG(dcr) = DEG(a[i]); COEF(dcr) = g;
                   1023:                                        f = s; substvp(vl,f,vn,&f0);
                   1024:                                } else
                   1025:                                        break;
                   1026:                        } else
                   1027:                                break;
                   1028:                }
                   1029:        }
                   1030:        *pp = f;
                   1031:        if ( dcr0 )
                   1032:                NEXT(dcr) = 0;
                   1033:        *dcp = dcr0;
                   1034: }
                   1035:
                   1036: void mfctrhen2(vl,vn,f,f0,g0,h0,lcg,lch,gp)
                   1037: VL vl;
                   1038: VN vn;
                   1039: P f;
                   1040: P f0,g0,h0,lcg,lch;
                   1041: P *gp;
                   1042: {
                   1043:        V v;
                   1044:        P f1,lc,lc0,lcg0,lch0;
                   1045:        P m,ff0,gg0,hh0,gk,hk,ggg,gggr,hhh,ak,bk,tmp;
                   1046:        Q bb,qk,s;
                   1047:        Q cbd;
                   1048:        int dbd;
                   1049:        int d;
                   1050:
                   1051:        if ( NUM(g0) ) {
                   1052:                *gp = (P)ONE;
                   1053:                return;
                   1054:        }
                   1055:
                   1056:        v = VR(f); d = deg(v,f);
                   1057:        if ( d == deg(v,g0) ) {
                   1058:                pcp(vl,f,gp,&tmp);
                   1059:                return;
                   1060:        }
                   1061:
                   1062:        mulp(vl,lcg,lch,&lc);
                   1063:        if ( !divtpz(vl,lc,LC(f),(P *)&s) ) {
                   1064:                *gp = 0;
                   1065:                return;
                   1066:        }
                   1067:        mulp(vl,(P)s,f,&f1);
                   1068:        dbd = dbound(VR(f1),f1) + 1; cbound(vl,f1,&cbd);
                   1069:
                   1070:        substvp(vl,lc,vn,&lc0);
                   1071:        divq((Q)lc0,(Q)LC(f0),(Q *)&m); mulp(vl,f0,m,&ff0);
                   1072:        substvp(vl,lcg,vn,&lcg0);
                   1073:        divq((Q)lcg0,(Q)LC(g0),(Q *)&m); mulp(vl,g0,m,&gg0);
                   1074:        substvp(vl,lch,vn,&lch0);
                   1075:        divq((Q)lch0,(Q)LC(h0),(Q *)&m); mulp(vl,h0,m,&hh0);
                   1076:        addq(cbd,cbd,&bb);
                   1077:        henzq1(gg0,hh0,bb,&bk,&ak,&qk); gk = gg0; hk = hh0;
                   1078:        henmv(vl,vn,f1,gk,hk,ak,bk,
                   1079:                lcg,lch,lcg0,lch0,qk,dbd,&ggg,&hhh);
                   1080:
                   1081:        if ( divtpz(vl,f1,ggg,&gggr) )
                   1082:                pcp(vl,ggg,gp,&tmp);
                   1083:        else
                   1084:                *gp = 0;
                   1085: }
                   1086:
                   1087: int sqfrchk(p)
                   1088: P p;
                   1089: {
                   1090:        Q c;
                   1091:        P f;
                   1092:        DCP dc;
                   1093:
                   1094:        ptozp(p,SGN((Q)UCOEF(p)),&c,&f); usqp(f,&dc);
                   1095:        if ( NEXT(dc) || !UNIQ(DEG(dc)) )
                   1096:                return ( 0 );
                   1097:        else
                   1098:                return ( 1 );
                   1099: }
                   1100:
                   1101: int cycchk(p)
                   1102: P p;
                   1103: {
                   1104:        Q c;
                   1105:        P f;
                   1106:
                   1107:        ptozp(p,SGN((Q)UCOEF(p)),&c,&f);
                   1108:        if ( iscycp(f) || iscycm(f) )
                   1109:                return 0;
                   1110:        else
                   1111:                return 1;
                   1112: }
                   1113:
                   1114: int zerovpchk(vl,p,vn)
                   1115: VL vl;
                   1116: P p;
                   1117: VN vn;
                   1118: {
                   1119:        P t;
                   1120:
                   1121:        substvp(vl,p,vn,&t);
                   1122:        if ( t )
                   1123:                return ( 0 );
                   1124:        else
                   1125:                return ( 1 );
                   1126: }
                   1127:
                   1128: int valideval(vl,dc,vn)
                   1129: VL vl;
                   1130: DCP dc;
                   1131: VN vn;
                   1132: {
                   1133:        DCP dct;
                   1134:        Q *a;
                   1135:        int i,j,n;
                   1136:        N q,r;
                   1137:
                   1138:        for ( dct = NEXT(dc), n = 0; dct; dct = NEXT(dct), n++ );
                   1139:        W_CALLOC(n,Q,a);
                   1140:        for ( dct = NEXT(dc), i = 0; dct; dct = NEXT(dct), i++ ) {
                   1141:                substvp(vl,COEF(dct),vn,(P *)&a[i]);
                   1142:                if ( !a[i] )
                   1143:                        return ( 0 );
                   1144:
                   1145:                for ( j = 0; j < i; j++ ) {
                   1146:                        divn(NM(a[j]),NM(a[i]),&q,&r);
                   1147:                        if ( !r )
                   1148:                                return ( 0 );
                   1149:                        divn(NM(a[i]),NM(a[j]),&q,&r);
                   1150:                        if ( !r )
                   1151:                                return ( 0 );
                   1152:                }
                   1153:        }
                   1154:        return ( 1 );
                   1155: }
                   1156:
                   1157: void estimatelc(vl,c,dc,vn,lcp)
                   1158: VL vl;
                   1159: Q c;
                   1160: DCP dc;
                   1161: VN vn;
                   1162: P *lcp;
                   1163: {
                   1164:        int i;
                   1165:        DCP dct;
                   1166:        P r,s,t;
                   1167:        Q c0,c1,c2;
                   1168:
                   1169:        for ( dct = dc, r = (P)ONE; dct; dct = NEXT(dct) ) {
                   1170:                if ( NUM(COEF(dct)) ) {
                   1171:                        mulp(vl,r,COEF(dct),&s); r = s;
                   1172:                } else {
                   1173:                        substvp(vl,COEF(dct),vn,(P *)&c0);
                   1174:                        for ( i = 0, c1 = c; i < (int)QTOS(DEG(dct)); i++ ) {
                   1175:                                divq(c1,c0,&c2);
                   1176:                                if ( !INT(c2) )
                   1177:                                        break;
                   1178:                                else
                   1179:                                        c1 = c2;
                   1180:                        }
                   1181:                        if ( i ) {
                   1182:                                STOQ(i,c1);
                   1183:                                pwrp(vl,COEF(dct),c1,&s); mulp(vl,r,s,&t); r = t;
                   1184:                        }
                   1185:                }
                   1186:        }
                   1187:        *lcp = r;
                   1188: }
                   1189:
                   1190: void monomialfctr(vl,p,pr,dcp)
                   1191: VL vl;
                   1192: P p;
                   1193: P *pr;
                   1194: DCP *dcp;
                   1195: {
                   1196:        VL nvl,avl;
                   1197:        Q d;
                   1198:        P f,t,s;
                   1199:        DCP dc0,dc;
                   1200:
                   1201:        clctv(vl,p,&nvl);
                   1202:        for ( dc0 = 0, avl = nvl, f = p; avl; avl = NEXT(avl) ) {
                   1203:                getmindeg(avl->v,f,&d);
                   1204:                if ( d ) {
                   1205:                        MKV(avl->v,t); NEXTDC(dc0,dc); DEG(dc) = d; COEF(dc) = t;
                   1206:                        pwrp(vl,t,d,&s); divsp(vl,f,s,&t); f = t;
                   1207:                }
                   1208:        }
                   1209:        if ( dc0 )
                   1210:                NEXT(dc) = 0;
                   1211:        *pr = f; *dcp = dc0;
                   1212: }
                   1213:
                   1214: void afctr(vl,p0,p,dcp)
                   1215: VL vl;
                   1216: P p,p0;
                   1217: DCP *dcp;
                   1218: {
                   1219:        DCP dc,dc0,dcr,dct,dcs;
                   1220:        P t;
                   1221:        VL nvl;
                   1222:
                   1223:        if ( VR(p) == VR(p0) ) {
                   1224:                NEWDC(dc);
                   1225:                DEG(dc) = ONE;
                   1226:                COEF(dc) = p;
                   1227:                NEXT(dc) = 0;
                   1228:                *dcp = dc;
                   1229:                return;
                   1230:        }
                   1231:
                   1232:        clctv(vl,p,&nvl);
                   1233:        if ( !NEXT(nvl) )
                   1234:                ufctr(p,1,&dc);
                   1235:        else {
                   1236:                sqa(vl,p0,p,&dc);
                   1237:                for ( dct = dc; dct; dct = NEXT(dct) ) {
                   1238:                        pmonic(vl,p0,COEF(dct),&t); COEF(dct) = t;
                   1239:                }
                   1240:        }
                   1241:        if ( NUM(COEF(dc)) )
                   1242:                dcr = NEXT(dc);
                   1243:        else
                   1244:                dcr = dc;
                   1245:        for ( dc0 = 0; dcr; dcr = NEXT(dcr) ) {
                   1246:                afctrmain(vl,p0,COEF(dcr),1,&dcs);
                   1247:
                   1248:                for ( dct = dcs; dct; dct = NEXT(dct) )
                   1249:                        DEG(dct) = DEG(dcr);
                   1250:                if ( !dc0 )
                   1251:                        dc0 = dcs;
                   1252:                else {
                   1253:                        for ( dct = dc0; NEXT(dct); dct = NEXT(dct) );
                   1254:                        NEXT(dct) = dcs;
                   1255:                }
                   1256:        }
                   1257:        *dcp = dc0;
                   1258: }
                   1259:
                   1260: void afctrmain(vl,p0,p,init,dcp)
                   1261: VL vl;
                   1262: P p,p0;
                   1263: int init;
                   1264: DCP *dcp;
                   1265: {
                   1266:        P x,y,s,m,a,t,u,pt,pt1,res,g;
                   1267:        Q q;
                   1268:        DCP dc,dc0,dcsq0,dcr0,dcr,dct,dcs;
                   1269:        V v,v0;
                   1270:
                   1271:        if ( !cmpq(DEG(DC(p)),ONE) ) {
                   1272:                NEWDC(dc); DEG(dc) = ONE; NEXT(dc) = 0;
                   1273:                pmonic(vl,p0,p,&COEF(dc)); *dcp = dc;
                   1274:                return;
                   1275:        }
                   1276:
                   1277:        v = VR(p); MKV(v,x);
                   1278:        v0 = VR(p0); MKV(v0,y);
                   1279:        STOQ(init,q),s = (P)q;
                   1280:        mulp(vl,s,y,&m); subp(vl,x,m,&t); addp(vl,x,m,&a);
                   1281:        substp(vl,p,v,t,&pt);
                   1282:        remsdcp(vl,pt,p0,&pt1);
                   1283:
                   1284: /*
                   1285:        if ( ( deg(v0,p0) <= 3 ) || ( TCPQ(p0) <= 2 ) )
                   1286:                resultp(vl,v0,p0,pt1,&res);
                   1287:        else
                   1288:                srcrnorm(vl,v0,pt1,p0,&res);
                   1289: */
                   1290: #if 0
                   1291:        srcr(vl,v0,pt1,p0,&res);
                   1292: #endif
                   1293:        resultp(vl,v0,p0,pt1,&res);
                   1294:        usqp(res,&dcsq0);
                   1295:        for ( dc0 = 0, dct = dcsq0; dct; dct = NEXT(dct) ) {
                   1296:                if ( UNIQ(DEG(dct)) )
                   1297:                        ufctr(COEF(dct),deg(v0,p0),&dcs);
                   1298:                else
                   1299:                        ufctr(COEF(dct),1,&dcs);
                   1300:                for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
                   1301:                        DEG(dcr) = DEG(dct);
                   1302:                if ( !dc0 ) {
                   1303:                        dc0 = NEXT(dcs);
                   1304:                        dc = dc0;
                   1305:                } else {
                   1306:                        for ( ; NEXT(dc); dc = NEXT(dc) );
                   1307:                        NEXT(dc) = NEXT(dcs);
                   1308:                }
                   1309:        }
                   1310:        sortfs(&dc0);
                   1311:
                   1312:        for ( g = p, dcr = dcr0 = 0, dc = dc0; dc; dc = NEXT(dc) ) {
                   1313:                if ( !UNIQ(DEG(dc)) ) {
                   1314:                        substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s);
                   1315:                        gcda(vl,p0,s,g,&u);
                   1316:                        if ( !NUM(u) && (VR(u) == v)) {
                   1317:                                afctrmain(vl,p0,u,init+1,&dct);
                   1318:                                for ( dcs = dct, t = (P)ONE; dcs; dcs = NEXT(dcs) ) {
                   1319:                                        mulp(vl,t,COEF(dcs),&s); remsdcp(vl,s,p0,&t);
                   1320:                                }
                   1321:                                pdiva(vl,p0,g,t,&s); g = s;
                   1322:                                if ( !dcr0 )
                   1323:                                        dcr0 = dct;
                   1324:                                else
                   1325:                                        NEXT(dcr) = dct;
                   1326:                                for ( dcr = dct; NEXT(dcr); dcr = NEXT(dcr) );
                   1327:                        }
                   1328:                } else {
                   1329:                        substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s);
                   1330:                        gcda(vl,p0,s,g,&u);
                   1331:                        if ( !NUM(u) && (VR(u) == v)) {
                   1332:                                NEXTDC(dcr0,dcr);
                   1333:                                DEG(dcr) = ONE;
                   1334:                                COEF(dcr) = u;
                   1335:                                pdiva(vl,p0,g,u,&t); g = t;
                   1336:                        }
                   1337:                }
                   1338:        }
                   1339:        if ( dcr0 )
                   1340:                NEXT(dcr) = 0;
                   1341:        *dcp = dcr0;
                   1342: }

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