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

Diff for /OpenXM_contrib2/asir2000/engine/F.c between version 1.12 and 1.13

version 1.12, 2013/01/08 07:25:58 version 1.13, 2018/03/29 01:32:51
Line 45 
Line 45 
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *   *
  * $OpenXM: OpenXM_contrib2/asir2000/engine/F.c,v 1.11 2002/01/15 01:09:55 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/F.c,v 1.12 2013/01/08 07:25:58 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include <math.h>  #include <math.h>
Line 54  int use_new_hensel;
Line 54  int use_new_hensel;
   
 void fctrp(VL vl,P f,DCP *dcp)  void fctrp(VL vl,P f,DCP *dcp)
 {  {
         VL nvl;    VL nvl;
         DCP dc;    DCP dc;
   
         if ( !f || NUM(f) ) {    if ( !f || NUM(f) ) {
                 NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;      NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
                 NEXT(dc) = 0; *dcp = dc;      NEXT(dc) = 0; *dcp = dc;
                 return;      return;
         } else if ( !qpcheck((Obj)f) )    } else if ( !qpcheck((Obj)f) )
                 error("fctrp : invalid argument");      error("fctrp : invalid argument");
         else {    else {
                 clctv(vl,f,&nvl);      clctv(vl,f,&nvl);
                 if ( !NEXT(nvl) )      if ( !NEXT(nvl) )
                         ufctr(f,1,dcp);        ufctr(f,1,dcp);
                 else      else
                         mfctr(nvl,f,dcp);        mfctr(nvl,f,dcp);
         }    }
 }  }
   
 void fctr_wrt_v_p(VL vl,P f,V v,DCP *dcp)  void fctr_wrt_v_p(VL vl,P f,V v,DCP *dcp)
 {  {
         VL nvl;    VL nvl;
         DCP dc;    DCP dc;
   
         if ( !f || NUM(f) ) {    if ( !f || NUM(f) ) {
                 NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;      NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
                 NEXT(dc) = 0; *dcp = dc;      NEXT(dc) = 0; *dcp = dc;
                 return;      return;
         } else if ( !qpcheck((Obj)f) )    } else if ( !qpcheck((Obj)f) )
                 error("fctrp : invalid argument");      error("fctrp : invalid argument");
         else {    else {
                 clctv(vl,f,&nvl);      clctv(vl,f,&nvl);
                 if ( !NEXT(nvl) )      if ( !NEXT(nvl) )
                         ufctr(f,1,dcp);        ufctr(f,1,dcp);
                 else      else
                         mfctr_wrt_v(nvl,f,v,dcp);        mfctr_wrt_v(nvl,f,v,dcp);
         }    }
 }  }
   
 void homfctr(VL vl,P g,DCP *dcp)  void homfctr(VL vl,P g,DCP *dcp)
 {  {
         P s,t,u,f,h,x;    P s,t,u,f,h,x;
         Q e;    Q e;
         int d,d0;    int d,d0;
         DCP dc,dct;    DCP dc,dct;
   
         substp(vl,g,vl->v,(P)ONE,&t); fctrp(vl,t,&dc); MKV(vl->v,x);    substp(vl,g,vl->v,(P)ONE,&t); fctrp(vl,t,&dc); MKV(vl->v,x);
         for ( dct = dc; dct; dct = NEXT(dct) )    for ( dct = dc; dct; dct = NEXT(dct) )
                 if ( !NUM(dc->c) ) {      if ( !NUM(dc->c) ) {
                         for ( s = 0, f = dc->c, d = d0 = homdeg(f); f; d = homdeg(f) ) {        for ( s = 0, f = dc->c, d = d0 = homdeg(f); f; d = homdeg(f) ) {
                                 exthp(vl,f,d,&h); STOQ(d0-d,e); pwrp(vl,x,e,&t);          exthp(vl,f,d,&h); STOQ(d0-d,e); pwrp(vl,x,e,&t);
                                 mulp(vl,t,h,&u); addp(vl,s,u,&t); s = t;          mulp(vl,t,h,&u); addp(vl,s,u,&t); s = t;
                                 subp(vl,f,h,&u); f = u;          subp(vl,f,h,&u); f = u;
                         }        }
                         dc->c = s;        dc->c = s;
                 }      }
         *dcp = dc;    *dcp = dc;
 }  }
   
 void mfctr(VL vl,P f,DCP *dcp)  void mfctr(VL vl,P f,DCP *dcp)
 {  {
         DCP dc,dc0,dct,dcs,dcr;    DCP dc,dc0,dct,dcs,dcr;
         P p,pmin,ppmin,cmin,t;    P p,pmin,ppmin,cmin,t;
         VL mvl;    VL mvl;
         Q c;    Q c;
   
         ptozp(f,1,&c,&p);    ptozp(f,1,&c,&p);
         NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;    NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
         msqfr(vl,p,&dct);    msqfr(vl,p,&dct);
         for ( ; dct; dct = NEXT(dct) ) {    for ( ; dct; dct = NEXT(dct) ) {
                 mindegp(vl,COEF(dct),&mvl,&pmin);      mindegp(vl,COEF(dct),&mvl,&pmin);
 #if 0  #if 0
                 minlcdegp(vl,COEF(dct),&mvl,&pmin);      minlcdegp(vl,COEF(dct),&mvl,&pmin);
                 min_common_vars_in_coefp(vl,COEF(dct),&mvl,&pmin);      min_common_vars_in_coefp(vl,COEF(dct),&mvl,&pmin);
 #endif  #endif
                 pcp(mvl,pmin,&ppmin,&cmin);      pcp(mvl,pmin,&ppmin,&cmin);
                 if ( !NUM(cmin) ) {      if ( !NUM(cmin) ) {
                         mfctr(mvl,cmin,&dcs);        mfctr(mvl,cmin,&dcs);
                         for ( dcr = NEXT(dcs); dcr; dcr = NEXT(dcr) ) {        for ( dcr = NEXT(dcs); dcr; dcr = NEXT(dcr) ) {
                                 DEG(dcr) = DEG(dct);          DEG(dcr) = DEG(dct);
                                 reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;          reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
                         }        }
                         for ( ; NEXT(dc); dc = NEXT(dc) );        for ( ; NEXT(dc); dc = NEXT(dc) );
                         NEXT(dc) = NEXT(dcs);        NEXT(dc) = NEXT(dcs);
                 }      }
                 mfctrmain(mvl,ppmin,&dcs);      mfctrmain(mvl,ppmin,&dcs);
                 for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {      for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
                         DEG(dcr) = DEG(dct);        DEG(dcr) = DEG(dct);
                         reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;        reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
                 }      }
                 for ( ; NEXT(dc); dc = NEXT(dc) );      for ( ; NEXT(dc); dc = NEXT(dc) );
                 NEXT(dc) = dcs;      NEXT(dc) = dcs;
         }    }
         adjsgn(f,dc0); *dcp = dc0;    adjsgn(f,dc0); *dcp = dc0;
 }  }
   
 void mfctr_wrt_v(VL vl,P f,V v,DCP *dcp)  void mfctr_wrt_v(VL vl,P f,V v,DCP *dcp)
 {  {
         DCP dc,dc0,dct,dcs,dcr;    DCP dc,dc0,dct,dcs,dcr;
         P p,pmin,ppmin,cmin,t;    P p,pmin,ppmin,cmin,t;
         VL nvl,mvl;    VL nvl,mvl;
         Q c;    Q c;
   
         ptozp(f,1,&c,&p);    ptozp(f,1,&c,&p);
         NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;    NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
         msqfr(vl,p,&dct);    msqfr(vl,p,&dct);
         for ( ; dct; dct = NEXT(dct) ) {    for ( ; dct; dct = NEXT(dct) ) {
                 clctv(vl,COEF(dct),&nvl);      clctv(vl,COEF(dct),&nvl);
                 reordvar(nvl,v,&mvl);      reordvar(nvl,v,&mvl);
                 reorderp(mvl,vl,COEF(dct),&pmin);      reorderp(mvl,vl,COEF(dct),&pmin);
                 pcp(mvl,pmin,&ppmin,&cmin);      pcp(mvl,pmin,&ppmin,&cmin);
                 if ( !NUM(cmin) ) {      if ( !NUM(cmin) ) {
                         mfctrmain(mvl,cmin,&dcs);        mfctrmain(mvl,cmin,&dcs);
                         for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {        for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
                                 DEG(dcr) = DEG(dct);          DEG(dcr) = DEG(dct);
                                 reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;          reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
                         }        }
                         for ( ; NEXT(dc); dc = NEXT(dc) );        for ( ; NEXT(dc); dc = NEXT(dc) );
                         NEXT(dc) = dcs;        NEXT(dc) = dcs;
                 }      }
                 mfctrmain(mvl,ppmin,&dcs);      mfctrmain(mvl,ppmin,&dcs);
                 for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {      for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
                         DEG(dcr) = DEG(dct);        DEG(dcr) = DEG(dct);
                         reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;        reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
                 }      }
                 for ( ; NEXT(dc); dc = NEXT(dc) );      for ( ; NEXT(dc); dc = NEXT(dc) );
                 NEXT(dc) = dcs;      NEXT(dc) = dcs;
         }    }
         adjsgn(f,dc0); *dcp = dc0;    adjsgn(f,dc0); *dcp = dc0;
 }  }
   
 #if 0  #if 0
 void adjsgn(P p,DCP dc)  void adjsgn(P p,DCP dc)
 {  {
         int sgn;    int sgn;
         DCP dct;    DCP dct;
         P c;    P c;
   
         for ( dct = dc, sgn = headsgn(p); dct; dct = NEXT(dct) )    for ( dct = dc, sgn = headsgn(p); dct; dct = NEXT(dct) )
                 if ( !EVENN(NM(dct->d)) )      if ( !EVENN(NM(dct->d)) )
                         sgn *= headsgn(COEF(dct));        sgn *= headsgn(COEF(dct));
         if ( sgn < 0 ) {    if ( sgn < 0 ) {
                 chsgnp(COEF(dc),&c); COEF(dc) = c;      chsgnp(COEF(dc),&c); COEF(dc) = c;
         }    }
 }  }
 #else  #else
 void adjsgn(P p,DCP dc)  void adjsgn(P p,DCP dc)
 {  {
         DCP dct;    DCP dct;
         P c;    P c;
   
         if ( headsgn(COEF(dc)) != headsgn(p) ) {    if ( headsgn(COEF(dc)) != headsgn(p) ) {
                 chsgnp(COEF(dc),&c); COEF(dc) = c;      chsgnp(COEF(dc),&c); COEF(dc) = c;
         }    }
         for ( dct = NEXT(dc); dct; dct = NEXT(dct) )    for ( dct = NEXT(dc); dct; dct = NEXT(dct) )
                 if ( headsgn(COEF(dct)) < 0 ) {      if ( headsgn(COEF(dct)) < 0 ) {
                         chsgnp(COEF(dct),&c); COEF(dct) = c;        chsgnp(COEF(dct),&c); COEF(dct) = c;
                 }      }
 }  }
 #endif  #endif
   
 int headsgn(P p)  int headsgn(P p)
 {  {
         if ( !p )    if ( !p )
                 return 0;      return 0;
         else {    else {
                 while ( !NUM(p) )      while ( !NUM(p) )
                         p = COEF(DC(p));        p = COEF(DC(p));
                 return SGN((Q)p);      return SGN((Q)p);
         }    }
 }  }
   
 void fctrwithmvp(VL vl,P f,V v,DCP *dcp)  void fctrwithmvp(VL vl,P f,V v,DCP *dcp)
 {  {
         VL nvl;    VL nvl;
         DCP dc;    DCP dc;
   
         if ( NUM(f) ) {    if ( NUM(f) ) {
                 NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;      NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
                 NEXT(dc) = 0; *dcp = dc;      NEXT(dc) = 0; *dcp = dc;
                 return;      return;
         }    }
   
         clctv(vl,f,&nvl);    clctv(vl,f,&nvl);
         if ( !NEXT(nvl) )    if ( !NEXT(nvl) )
                 ufctr(f,1,dcp);      ufctr(f,1,dcp);
         else    else
                 mfctrwithmv(nvl,f,v,dcp);      mfctrwithmv(nvl,f,v,dcp);
 }  }
   
 void mfctrwithmv(VL vl,P f,V v,DCP *dcp)  void mfctrwithmv(VL vl,P f,V v,DCP *dcp)
 {  {
         DCP dc,dc0,dct,dcs,dcr;    DCP dc,dc0,dct,dcs,dcr;
         P p,pmin,ppmin,cmin,t;    P p,pmin,ppmin,cmin,t;
         VL mvl;    VL mvl;
         Q c;    Q c;
   
         ptozp(f,1,&c,&p);    ptozp(f,1,&c,&p);
         NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;    NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
         msqfr(vl,p,&dct);    msqfr(vl,p,&dct);
         for ( ; dct; dct = NEXT(dct) ) {    for ( ; dct; dct = NEXT(dct) ) {
                 if ( !v )      if ( !v )
                         mindegp(vl,COEF(dct),&mvl,&pmin);        mindegp(vl,COEF(dct),&mvl,&pmin);
                 else {      else {
                         reordvar(vl,v,&mvl); reorderp(mvl,vl,COEF(dct),&pmin);        reordvar(vl,v,&mvl); reorderp(mvl,vl,COEF(dct),&pmin);
                 }      }
                 pcp(mvl,pmin,&ppmin,&cmin);      pcp(mvl,pmin,&ppmin,&cmin);
                 if ( !NUM(cmin) ) {      if ( !NUM(cmin) ) {
                         mfctrmain(mvl,cmin,&dcs);        mfctrmain(mvl,cmin,&dcs);
                         for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {        for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
                                 DEG(dcr) = DEG(dct);          DEG(dcr) = DEG(dct);
                                 reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;          reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
                         }        }
                         for ( ; NEXT(dc); dc = NEXT(dc) );        for ( ; NEXT(dc); dc = NEXT(dc) );
                         NEXT(dc) = dcs;        NEXT(dc) = dcs;
                 }      }
                 mfctrmain(mvl,ppmin,&dcs);      mfctrmain(mvl,ppmin,&dcs);
                 for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {      for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
                         DEG(dcr) = DEG(dct);        DEG(dcr) = DEG(dct);
                         reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;        reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
                 }      }
                 for ( ; NEXT(dc); dc = NEXT(dc) );      for ( ; NEXT(dc); dc = NEXT(dc) );
                 NEXT(dc) = dcs;      NEXT(dc) = dcs;
         }    }
         *dcp = dc0;    *dcp = dc0;
 }  }
   
 void ufctr(P f,int hint,DCP *dcp)  void ufctr(P f,int hint,DCP *dcp)
 {  {
         P p,c;    P p,c;
         DCP dc,dct,dcs,dcr;    DCP dc,dct,dcs,dcr;
   
         ptozp(f,SGN((Q)UCOEF(f)),(Q *)&c,&p);    ptozp(f,SGN((Q)UCOEF(f)),(Q *)&c,&p);
         NEWDC(dc); COEF(dc) = c; DEG(dc) = ONE; NEXT(dc) = 0;    NEWDC(dc); COEF(dc) = c; DEG(dc) = ONE; NEXT(dc) = 0;
         usqp(p,&dct);    usqp(p,&dct);
         for ( *dcp = dc; dct; dct = NEXT(dct) ) {    for ( *dcp = dc; dct; dct = NEXT(dct) ) {
                 ufctrmain(COEF(dct),hint,&dcs);      ufctrmain(COEF(dct),hint,&dcs);
                 for ( dcr = dcs; dcr; dcr = NEXT(dcr) )      for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
                         DEG(dcr) = DEG(dct);        DEG(dcr) = DEG(dct);
                 for ( ; NEXT(dc); dc = NEXT(dc) );      for ( ; NEXT(dc); dc = NEXT(dc) );
                 NEXT(dc) = dcs;      NEXT(dc) = dcs;
         }    }
 }  }
   
 void mfctrmain(VL vl,P p,DCP *dcp)  void mfctrmain(VL vl,P p,DCP *dcp)
 {  {
         int i,j,k,*win,np,x;    int i,j,k,*win,np,x;
         VL nvl,tvl;    VL nvl,tvl;
         VN vn,vnt,vn1;    VN vn,vnt,vn1;
         P p0,f,g,f0,g0,s,t,lp,m;    P p0,f,g,f0,g0,s,t,lp,m;
         P *fp0,*fpt,*l,*tl;    P *fp0,*fpt,*l,*tl;
         DCP dc,dc0,dcl;    DCP dc,dc0,dcl;
         int count,nv;    int count,nv;
         int *nonzero;    int *nonzero;
   
         if ( !cmpq(DEG(DC(p)),ONE) ) {    if ( !cmpq(DEG(DC(p)),ONE) ) {
                 NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;      NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                 *dcp = dc; return;      *dcp = dc; return;
         }    }
         lp = LC(p); fctrp(vl,lp,&dcl); clctv(vl,p,&nvl);    lp = LC(p); fctrp(vl,lp,&dcl); clctv(vl,p,&nvl);
         for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);    for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
         W_CALLOC(nv,struct oVN,vn); W_CALLOC(nv,struct oVN,vnt);    W_CALLOC(nv,struct oVN,vn); W_CALLOC(nv,struct oVN,vnt);
         W_CALLOC(nv,struct oVN,vn1);    W_CALLOC(nv,struct oVN,vn1);
         W_CALLOC(nv,int,nonzero);    W_CALLOC(nv,int,nonzero);
   
         for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )    for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
                 vn1[i].v = vn[i].v = tvl->v;      vn1[i].v = vn[i].v = tvl->v;
         vn1[i].v = vn[i].v = 0;    vn1[i].v = vn[i].v = 0;
   
         /* determine a heuristic bound of deg(GCD(p,p')) */    /* determine a heuristic bound of deg(GCD(p,p')) */
         while ( 1 ) {    while ( 1 ) {
                 for ( i = 0; vn1[i].v; i++ )      for ( i = 0; vn1[i].v; i++ )
                         vn1[i].n = ((unsigned int)random())%256+1;        vn1[i].n = ((unsigned int)random())%256+1;
                 substvp(nvl,LC(p),vn1,&p0);      substvp(nvl,LC(p),vn1,&p0);
                 if ( p0 ) {      if ( p0 ) {
                         substvp(nvl,p,vn1,&p0);        substvp(nvl,p,vn1,&p0);
                         if ( sqfrchk(p0) ) {        if ( sqfrchk(p0) ) {
                                 ufctr(p0,1,&dc0);          ufctr(p0,1,&dc0);
                                 if ( NEXT(NEXT(dc0)) == 0 ) {          if ( NEXT(NEXT(dc0)) == 0 ) {
                                         NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;            NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                                         *dcp = dc;            *dcp = dc;
                                         return;            return;
                                 } else {          } else {
                                         for ( dc0 = NEXT(dc0), np = 0; dc0; dc0 = NEXT(dc0), np++ );            for ( dc0 = NEXT(dc0), np = 0; dc0; dc0 = NEXT(dc0), np++ );
                                         break;            break;
                                 }          }
                         }        }
                 }      }
         }    }
   
         /* determine the position of variables which is not allowed to    /* determine the position of variables which is not allowed to
            be set to 0 */       be set to 0 */
         for ( i = 0; vn1[i].v; i++ ) {    for ( i = 0; vn1[i].v; i++ ) {
                 x = vn1[i].n; vn1[i].n = 0;      x = vn1[i].n; vn1[i].n = 0;
                 substvp(nvl,LC(p),vn1,&p0);      substvp(nvl,LC(p),vn1,&p0);
                 if ( !p0 )      if ( !p0 )
                         vn1[i].n = x;        vn1[i].n = x;
                 else {      else {
                         substvp(nvl,p,vn1,&p0);        substvp(nvl,p,vn1,&p0);
                         if ( !sqfrchk(p0) )        if ( !sqfrchk(p0) )
                                 vn1[i].n = x;          vn1[i].n = x;
                         else {        else {
                                 ufctr(p0,1,&dc0);          ufctr(p0,1,&dc0);
                                 for ( dc0 = NEXT(dc0), j = 0; dc0; dc0 = NEXT(dc0), j++ );          for ( dc0 = NEXT(dc0), j = 0; dc0; dc0 = NEXT(dc0), j++ );
                                 if ( j > np )          if ( j > np )
                                         vn1[i].n = x;            vn1[i].n = x;
                         }        }
                 }      }
         }    }
         for ( i = 0; vn1[i].v; i++ )    for ( i = 0; vn1[i].v; i++ )
                 if (vn1[i].n )      if (vn1[i].n )
                         nonzero[i] = 1;        nonzero[i] = 1;
   
         count = 0;    count = 0;
         while  ( 1 ) {    while  ( 1 ) {
                 while ( 1 ) {      while ( 1 ) {
                         count++;        count++;
                         for ( i = 0, j = 0; vn[i].v; i++ )        for ( i = 0, j = 0; vn[i].v; i++ )
                                 if ( vn[i].n )          if ( vn[i].n )
                                         vnt[j++].v = (V)i;            vnt[j++].v = (V)i;
                         vnt[j].n = 0; mulsgn(vn,vnt,j,vn1);        vnt[j].n = 0; mulsgn(vn,vnt,j,vn1);
                         for ( i = 0; vn1[i].v; i++ )        for ( i = 0; vn1[i].v; i++ )
                                 if ( vn1[i].n ) {          if ( vn1[i].n ) {
                                         if ( vn1[i].n > 0 )            if ( vn1[i].n > 0 )
                                                 vn1[i].n = sprime[vn1[i].n];              vn1[i].n = sprime[vn1[i].n];
                                         else            else
                                                 vn1[i].n = -sprime[-vn1[i].n];              vn1[i].n = -sprime[-vn1[i].n];
                                 }          }
                         if ( valideval(nvl,dcl,vn1) ) {        if ( valideval(nvl,dcl,vn1) ) {
                                 substvp(nvl,p,vn1,&p0);          substvp(nvl,p,vn1,&p0);
                                 if ( sqfrchk(p0) ) {          if ( sqfrchk(p0) ) {
                                         ufctr(p0,1,&dc0);            ufctr(p0,1,&dc0);
                                         if ( NEXT(NEXT(dc0)) == 0 ) {            if ( NEXT(NEXT(dc0)) == 0 ) {
                                                 NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;              NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                                                 *dcp = dc;              *dcp = dc;
                                                 return;              return;
                                         } else {            } else {
                                                 for ( dc = NEXT(dc0), i = 0; dc; dc = NEXT(dc), i++ );              for ( dc = NEXT(dc0), i = 0; dc; dc = NEXT(dc), i++ );
                                                 if ( i <= np )              if ( i <= np )
                                                         goto MAIN;                goto MAIN;
                                                 if ( i < np )              if ( i < np )
                                                         np = i;                np = i;
                                         }            }
                                 }          }
                         }        }
                         if ( nextbin(vnt,j) )        if ( nextbin(vnt,j) )
                                 break;          break;
                 }      }
                 while ( 1 ) {      while ( 1 ) {
                         next(vn);        next(vn);
                         for ( i = 0; vn[i].v; i++ )        for ( i = 0; vn[i].v; i++ )
                                 if ( nonzero[i] && !vn[i].n )          if ( nonzero[i] && !vn[i].n )
                                         break;            break;
                         if ( !vn[i].v )        if ( !vn[i].v )
                                 break;          break;
                 }      }
         }    }
 MAIN :  MAIN :
 #if 0  #if 0
         for ( i = 0; vn1[i].v; i++ )    for ( i = 0; vn1[i].v; i++ )
                 fprintf(stderr,"%d ",vn1[i].n);      fprintf(stderr,"%d ",vn1[i].n);
         fprintf(stderr,"\n");    fprintf(stderr,"\n");
 #endif  #endif
         for ( np = 0, dc = NEXT(dc0); dc; dc = NEXT(dc), np++ );    for ( np = 0, dc = NEXT(dc0); dc; dc = NEXT(dc), np++ );
         fp0 = (P *)ALLOCA((np + 1)*sizeof(P));    fp0 = (P *)ALLOCA((np + 1)*sizeof(P));
         fpt = (P *)ALLOCA((np + 1)*sizeof(P));    fpt = (P *)ALLOCA((np + 1)*sizeof(P));
         l = tl = (P *)ALLOCA((np + 1)*sizeof(P));    l = tl = (P *)ALLOCA((np + 1)*sizeof(P));
         for ( i = 1, dc = NEXT(dc0); dc; i++, dc = NEXT(dc) )    for ( i = 1, dc = NEXT(dc0); dc; i++, dc = NEXT(dc) )
                 fp0[i-1] = COEF(dc);      fp0[i-1] = COEF(dc);
 #if 0  #if 0
         sort_by_deg(np,fp0,fpt);    sort_by_deg(np,fp0,fpt);
         sort_by_deg_rev(np-1,fpt+1,fp0+1);    sort_by_deg_rev(np-1,fpt+1,fp0+1);
 #endif  #endif
 #if 0  #if 0
         sort_by_deg_rev(np,fp0,fpt); fp0[0] = fpt[0];    sort_by_deg_rev(np,fp0,fpt); fp0[0] = fpt[0];
         sort_by_deg(np-1,fpt+1,fp0+1); fp0[np] = 0;    sort_by_deg(np-1,fpt+1,fp0+1); fp0[np] = 0;
 #endif  #endif
         fp0[np] = 0;    fp0[np] = 0;
         win = W_ALLOC(np + 1); f = p; divsp(vl,p0,COEF(dc0),&f0);    win = W_ALLOC(np + 1); f = p; divsp(vl,p0,COEF(dc0),&f0);
         for ( k = 1, win[0] = 1, --np; ; ) {    for ( k = 1, win[0] = 1, --np; ; ) {
                 P h0,lcg,lch;      P h0,lcg,lch;
                 Q c;      Q c;
   
                 g0 = fp0[win[0]];      g0 = fp0[win[0]];
                 for ( i = 1; i < k; i++ ) {      for ( i = 1; i < k; i++ ) {
                         mulp(vl,g0,fp0[win[i]],&m); g0 = m;        mulp(vl,g0,fp0[win[i]],&m); g0 = m;
                 }      }
                 mulq((Q)LC(g0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lcg);      mulq((Q)LC(g0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lcg);
                 divsp(nvl,f0,g0,&h0);      divsp(nvl,f0,g0,&h0);
                 mulq((Q)LC(h0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lch);      mulq((Q)LC(h0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lch);
                 mfctrhen2(nvl,vn1,f,f0,g0,h0,lcg,lch,&g);      mfctrhen2(nvl,vn1,f,f0,g0,h0,lcg,lch,&g);
                 if ( g ) {      if ( g ) {
                         *tl++ = g; divsp(vl,f,g,&t);        *tl++ = g; divsp(vl,f,g,&t);
                         f = t; divsp(vl,f0,g0,&t); ptozp(t,1,(Q *)&s,&f0);        f = t; divsp(vl,f0,g0,&t); ptozp(t,1,(Q *)&s,&f0);
                         for ( i = 0; i < k - 1; i++ )        for ( i = 0; i < k - 1; i++ )
                                 for ( j = win[i] + 1; j < win[i + 1]; j++ )          for ( j = win[i] + 1; j < win[i + 1]; j++ )
                                         fp0[j - i - 1] = fp0[j];            fp0[j - i - 1] = fp0[j];
                         for ( j = win[k - 1] + 1; j <= np; j++ )        for ( j = win[k - 1] + 1; j <= np; j++ )
                                         fp0[j - k] = fp0[j];            fp0[j - k] = fp0[j];
                         if ( ( np -= k ) < k ) break;        if ( ( np -= k ) < k ) break;
                         if ( np - win[0] + 1 < k )        if ( np - win[0] + 1 < k )
                                 if ( ++k <= np ) {          if ( ++k <= np ) {
                                         for ( i = 0; i < k; i++ ) win[i] = i + 1;            for ( i = 0; i < k; i++ ) win[i] = i + 1;
                                         continue;            continue;
                                 } else          } else
                                         break;            break;
                         else for ( i = 1; i < k; i++ ) win[i] = win[0] + i;        else for ( i = 1; i < k; i++ ) win[i] = win[0] + i;
                 } else {      } else {
                         if ( ncombi(1,np,k,win) == 0 )        if ( ncombi(1,np,k,win) == 0 )
                                 if ( k == np ) break;          if ( k == np ) break;
                                 else for ( i = 0, ++k; i < k; i++ ) win[i] = i + 1;          else for ( i = 0, ++k; i < k; i++ ) win[i] = i + 1;
                 }      }
         }    }
         *tl++ = f; *tl = 0;    *tl++ = f; *tl = 0;
         for ( dc0 = 0; *l; l++ ) {    for ( dc0 = 0; *l; l++ ) {
                 NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = *l;      NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = *l;
         }    }
         NEXT(dc) = 0; *dcp = dc0;    NEXT(dc) = 0; *dcp = dc0;
 }  }
   
 void ufctrmain(P p,int hint,DCP *dcp)  void ufctrmain(P p,int hint,DCP *dcp)
 {  {
         ML list;    ML list;
         DCP dc;    DCP dc;
   
         if ( NUM(p) || (UDEG(p) == 1) ) {    if ( NUM(p) || (UDEG(p) == 1) ) {
                 NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;      NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                 *dcp = dc;      *dcp = dc;
         } else if ( iscycm(p) )    } else if ( iscycm(p) )
                 cycm(VR(p),UDEG(p),dcp);      cycm(VR(p),UDEG(p),dcp);
         else if ( iscycp(p) )    else if ( iscycp(p) )
                 cycp(VR(p),UDEG(p),dcp);      cycp(VR(p),UDEG(p),dcp);
         else {    else {
                 if ( use_new_hensel )      if ( use_new_hensel )
                         hensel2(5,5,p,&list);        hensel2(5,5,p,&list);
                 else      else
                         hensel(5,5,p,&list);        hensel(5,5,p,&list);
                 if ( list->n == 1 ) {      if ( list->n == 1 ) {
                         NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;        NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                         *dcp = dc;        *dcp = dc;
                 } else      } else
                         dtest(p,list,hint,dcp);        dtest(p,list,hint,dcp);
         }    }
 }  }
   
 void cycm(V v,int n,DCP *dcp)  void cycm(V v,int n,DCP *dcp)
 {  {
         register int i,j;    register int i,j;
         struct oMF *mfp;    struct oMF *mfp;
         DCP dc,dc0;    DCP dc,dc0;
   
         if ( n == 1 ) {    if ( n == 1 ) {
                 NEWDC(dc); mkssum(v,1,1,-1,&COEF(dc)); DEG(dc) = ONE;      NEWDC(dc); mkssum(v,1,1,-1,&COEF(dc)); DEG(dc) = ONE;
         } else {    } else {
                 mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));      mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
                 bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));      bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
                 for ( i = 1, j = 0; i <= n; i++ )      for ( i = 1, j = 0; i <= n; i++ )
                         if ( !(n%i) ) mfp[j++].m = i;        if ( !(n%i) ) mfp[j++].m = i;
                 mkssum(v,1,1,-1,&mfp[0].f);      mkssum(v,1,1,-1,&mfp[0].f);
                 for ( i = 1; i < j; i++ )      for ( i = 1; i < j; i++ )
                         calcphi(v,i,mfp);        calcphi(v,i,mfp);
                 for ( dc0 = 0, i = 0; i < j; i++ ) {      for ( dc0 = 0, i = 0; i < j; i++ ) {
                         NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;        NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
                 }      }
         }    }
         NEXT(dc) = 0; *dcp = dc0;    NEXT(dc) = 0; *dcp = dc0;
 }  }
   
 void cycp(V v,int n,DCP *dcp)  void cycp(V v,int n,DCP *dcp)
 {  {
         register int i,j;    register int i,j;
         int n0;    int n0;
         struct oMF *mfp;    struct oMF *mfp;
         DCP dc,dc0;    DCP dc,dc0;
   
         if ( n == 1 ) {    if ( n == 1 ) {
                 NEWDC(dc); mkssum(v,1,1,1,&COEF(dc)); DEG(dc) = ONE;      NEWDC(dc); mkssum(v,1,1,1,&COEF(dc)); DEG(dc) = ONE;
         } else {    } else {
                 n0 = n; n *= 2;      n0 = n; n *= 2;
                 mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));      mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
                 bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));      bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
                 for ( i = 1, j = 0; i <= n; i++ )      for ( i = 1, j = 0; i <= n; i++ )
                         if ( !(n%i) ) mfp[j++].m = i;        if ( !(n%i) ) mfp[j++].m = i;
                 mkssum(v,1,1,-1,&mfp[0].f);      mkssum(v,1,1,-1,&mfp[0].f);
                 for ( i = 1; i < j; i++ )      for ( i = 1; i < j; i++ )
                         calcphi(v,i,mfp);        calcphi(v,i,mfp);
                 for ( dc0 = 0, i = 0; i < j; i++ )      for ( dc0 = 0, i = 0; i < j; i++ )
                         if ( n0 % mfp[i].m ) {        if ( n0 % mfp[i].m ) {
                                 NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;          NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
                         }        }
         }    }
         NEXT(dc) = 0; *dcp = dc0;    NEXT(dc) = 0; *dcp = dc0;
 }  }
   
 void calcphi(V v,int n,struct oMF *mfp)  void calcphi(V v,int n,struct oMF *mfp)
 {  {
         register int i,m;    register int i,m;
         P t,s,tmp;    P t,s,tmp;
   
         for ( t = (P)ONE, i = 0, m = mfp[n].m; i < n; i++ )    for ( t = (P)ONE, i = 0, m = mfp[n].m; i < n; i++ )
                 if ( !(m%mfp[i].m) ) {      if ( !(m%mfp[i].m) ) {
                         mulp(CO,t,mfp[i].f,&s); t = s;        mulp(CO,t,mfp[i].f,&s); t = s;
                 }      }
         mkssum(v,m,1,-1,&s); udivpz(s,t,&mfp[n].f,&tmp);    mkssum(v,m,1,-1,&s); udivpz(s,t,&mfp[n].f,&tmp);
         if ( tmp )    if ( tmp )
                 error("calcphi: cannot happen");      error("calcphi: cannot happen");
 }  }
   
 void mkssum(V v,int e,int s,int sgn,P *r)  void mkssum(V v,int e,int s,int sgn,P *r)
 {  {
         register int i,sgnt;    register int i,sgnt;
         DCP dc,dc0;    DCP dc,dc0;
         Q q;    Q q;
   
         for ( dc0 = 0, i = s, sgnt = 1; i >= 0; i--, sgnt *= sgn ) {    for ( dc0 = 0, i = s, sgnt = 1; i >= 0; i--, sgnt *= sgn ) {
                 if ( !dc0 ) {      if ( !dc0 ) {
                         NEWDC(dc0); dc = dc0;        NEWDC(dc0); dc = dc0;
                 } else {      } else {
                         NEWDC(NEXT(dc)); dc = NEXT(dc);        NEWDC(NEXT(dc)); dc = NEXT(dc);
                 }      }
                 STOQ(i*e,DEG(dc)); STOQ(sgnt,q),COEF(dc) = (P)q;      STOQ(i*e,DEG(dc)); STOQ(sgnt,q),COEF(dc) = (P)q;
         }    }
         NEXT(dc) = 0; MKP(v,dc0,*r);    NEXT(dc) = 0; MKP(v,dc0,*r);
 }  }
   
 int iscycp(P f)  int iscycp(P f)
 {  {
         DCP dc;    DCP dc;
         dc = DC(f);    dc = DC(f);
   
         if ( !UNIQ((Q)COEF(dc)) )    if ( !UNIQ((Q)COEF(dc)) )
                 return ( 0 );      return ( 0 );
         dc = NEXT(dc);    dc = NEXT(dc);
         if ( NEXT(dc) || DEG(dc) || !UNIQ((Q)COEF(dc)) )    if ( NEXT(dc) || DEG(dc) || !UNIQ((Q)COEF(dc)) )
                 return ( 0 );      return ( 0 );
         return ( 1 );    return ( 1 );
 }  }
   
 int iscycm(P f)  int iscycm(P f)
 {  {
         DCP dc;    DCP dc;
   
         dc = DC(f);    dc = DC(f);
         if ( !UNIQ((Q)COEF(dc)) )    if ( !UNIQ((Q)COEF(dc)) )
                 return ( 0 );      return ( 0 );
         dc = NEXT(dc);    dc = NEXT(dc);
         if ( NEXT(dc) || DEG(dc) || !MUNIQ((Q)COEF(dc)) )    if ( NEXT(dc) || DEG(dc) || !MUNIQ((Q)COEF(dc)) )
                 return ( 0 );      return ( 0 );
         return ( 1 );    return ( 1 );
 }  }
   
 void sortfs(DCP *dcp)  void sortfs(DCP *dcp)
 {  {
         int i,k,n,k0,d;    int i,k,n,k0,d;
         DCP dc,dct,t;    DCP dc,dct,t;
         DCP *a;    DCP *a;
   
         dc = *dcp;    dc = *dcp;
         for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );    for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
         a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));    a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
         for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ )    for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ )
                 a[i] = dct;      a[i] = dct;
         a[n] = 0;    a[n] = 0;
   
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ )      for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ )
                         if ( (int)UDEG(COEF(a[k])) < d ) {        if ( (int)UDEG(COEF(a[k])) < d ) {
                                 k0 = k;          k0 = k;
                                 d = UDEG(COEF(a[k]));          d = UDEG(COEF(a[k]));
                         }        }
                 if ( i != k0 ) {      if ( i != k0 ) {
                         t = a[i]; a[i] = a[k0]; a[k0] = t;        t = a[i]; a[i] = a[k0]; a[k0] = t;
                 }      }
         }    }
         for ( *dcp = dct = a[0], i = 1; i < n; i++ )    for ( *dcp = dct = a[0], i = 1; i < n; i++ )
                 dct = NEXT(dct) = a[i];      dct = NEXT(dct) = a[i];
         NEXT(dct) = 0;    NEXT(dct) = 0;
 }  }
   
 void sortfsrev(DCP *dcp)  void sortfsrev(DCP *dcp)
 {  {
         int i,k,n,k0,d;    int i,k,n,k0,d;
         DCP dc,dct,t;    DCP dc,dct,t;
         DCP *a;    DCP *a;
   
         dc = *dcp;    dc = *dcp;
         for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );    for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
         a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));    a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
         for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ )    for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ )
                 a[i] = dct;      a[i] = dct;
         a[n] = 0;    a[n] = 0;
   
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ )      for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ )
                         if ( (int)UDEG(COEF(a[k])) > d ) {        if ( (int)UDEG(COEF(a[k])) > d ) {
                                 k0 = k;          k0 = k;
                                 d = UDEG(COEF(a[k]));          d = UDEG(COEF(a[k]));
                         }        }
                 if ( i != k0 ) {      if ( i != k0 ) {
                         t = a[i]; a[i] = a[k0]; a[k0] = t;        t = a[i]; a[i] = a[k0]; a[k0] = t;
                 }      }
         }    }
         for ( *dcp = dct = a[0], i = 1; i < n; i++ )    for ( *dcp = dct = a[0], i = 1; i < n; i++ )
                 dct = NEXT(dct) = a[i];      dct = NEXT(dct) = a[i];
         NEXT(dct) = 0;    NEXT(dct) = 0;
 }  }
   
 void nthrootchk(P f,struct oDUM *dc,ML fp,DCP *dcp)  void nthrootchk(P f,struct oDUM *dc,ML fp,DCP *dcp)
 {  {
         register int i,k;    register int i,k;
         int e,n,dr,tmp,t;    int e,n,dr,tmp,t;
         int *tmpp,**tmppp;    int *tmpp,**tmppp;
         int **pp,**wpp;    int **pp,**wpp;
         LUM fpa,tpa,f0l;    LUM fpa,tpa,f0l;
         UM wt,wq,ws,dvr,f0,fe;    UM wt,wq,ws,dvr,f0,fe;
         N lc;    N lc;
         int lcm;    int lcm;
         int m,b;    int m,b;
   
         m = fp->mod; b = fp->bound; fe = *((UM *)fp->c);    m = fp->mod; b = fp->bound; fe = *((UM *)fp->c);
         e = dc->n; f0 = dc->f; nthrootn(NM((Q)COEF(DC(f))),e,&lc);    e = dc->n; f0 = dc->f; nthrootn(NM((Q)COEF(DC(f))),e,&lc);
         if ( !lc ) {    if ( !lc ) {
                 *dcp = 0;      *dcp = 0;
                 return;      return;
         }    }
         lcm = rem(lc,m); W_LUMALLOC(DEG(f0),b,f0l);    lcm = rem(lc,m); W_LUMALLOC(DEG(f0),b,f0l);
         for ( i = DEG(f0), tmppp = COEF(f0l), tmpp = COEF(f0);    for ( i = DEG(f0), tmppp = COEF(f0l), tmpp = COEF(f0);
                 i >= 0; i-- )      i >= 0; i-- )
                 *(tmppp[i]) = dmb(m,tmpp[i],lcm,&tmp);      *(tmppp[i]) = dmb(m,tmpp[i],lcm,&tmp);
         dtestroot(m,1,f,f0l,dc,dcp);    dtestroot(m,1,f,f0l,dc,dcp);
         if ( *dcp )    if ( *dcp )
                 return;      return;
         n = UDEG(f); W_LUMALLOC(n,b,fpa); W_LUMALLOC(n,b,tpa);    n = UDEG(f); W_LUMALLOC(n,b,fpa); W_LUMALLOC(n,b,tpa);
         ptolum(m,b,f,fpa);    ptolum(m,b,f,fpa);
         dvr = W_UMALLOC(n); wq = W_UMALLOC(n);    dvr = W_UMALLOC(n); wq = W_UMALLOC(n);
         wt = W_UMALLOC(n); ws = W_UMALLOC(n);    wt = W_UMALLOC(n); ws = W_UMALLOC(n);
         cpyum(fe,dvr); divum(m,dvr,f0,wq);    cpyum(fe,dvr); divum(m,dvr,f0,wq);
         t = dmb(m,pwrm(m,lcm,e-1),e,&tmp); DEG(dvr) = DEG(wq);    t = dmb(m,pwrm(m,lcm,e-1),e,&tmp); DEG(dvr) = DEG(wq);
         for ( k = 0; k <= DEG(wq); k++ )    for ( k = 0; k <= DEG(wq); k++ )
                 COEF(dvr)[k] = dmb(m,COEF(wq)[k],t,&tmp);      COEF(dvr)[k] = dmb(m,COEF(wq)[k],t,&tmp);
         for ( i = 1; i < b; i++ ) {    for ( i = 1; i < b; i++ ) {
                 pwrlum(m,i+1,f0l,e,tpa);      pwrlum(m,i+1,f0l,e,tpa);
                 for ( k = 0, pp = COEF(fpa), wpp = COEF(tpa);      for ( k = 0, pp = COEF(fpa), wpp = COEF(tpa);
                         k <= n; k++ )        k <= n; k++ )
                         COEF(wt)[k] = (pp[k][i]-wpp[k][i]+m)%m;        COEF(wt)[k] = (pp[k][i]-wpp[k][i]+m)%m;
                 degum(wt,n); dr = divum(m,wt,dvr,ws);      degum(wt,n); dr = divum(m,wt,dvr,ws);
                 if ( dr >= 0 ) {      if ( dr >= 0 ) {
                         *dcp = 0;        *dcp = 0;
                         return;        return;
                 }      }
                 for ( k = 0, pp = COEF(f0l); k <= DEG(ws); k++ )      for ( k = 0, pp = COEF(f0l); k <= DEG(ws); k++ )
                         pp[k][i] = COEF(ws)[k];        pp[k][i] = COEF(ws)[k];
                 dtestroot(m,i+1,f,f0l,dc,dcp);      dtestroot(m,i+1,f,f0l,dc,dcp);
                 if ( *dcp )      if ( *dcp )
                         return;        return;
         }    }
 }  }
   
 void sqfrp(VL vl,P f,DCP *dcp)  void sqfrp(VL vl,P f,DCP *dcp)
 {  {
         P c,p;    P c,p;
         DCP dc,dc0;    DCP dc,dc0;
   
         if ( !f || NUM(f) ) {    if ( !f || NUM(f) ) {
                 NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;      NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
                 COEF(dc0) = f; NEXT(dc0) = 0;      COEF(dc0) = f; NEXT(dc0) = 0;
         } else if ( !qpcheck((Obj)f) )    } else if ( !qpcheck((Obj)f) )
                 error("sqfrp : invalid argument");      error("sqfrp : invalid argument");
         else {    else {
                 NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;      NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
                 ptozp(f,1,(Q *)&c,&p); msqfr(vl,p,&dc);      ptozp(f,1,(Q *)&c,&p); msqfr(vl,p,&dc);
                 COEF(dc0) = c; NEXT(dc0) = dc;      COEF(dc0) = c; NEXT(dc0) = dc;
                 adjsgn(f,dc0);      adjsgn(f,dc0);
         }    }
 }  }
   
 /*  /*
Line 734  void sqfrp(VL vl,P f,DCP *dcp)
Line 734  void sqfrp(VL vl,P f,DCP *dcp)
  */   */
 void msqfr(VL vl,P f,DCP *dcp)  void msqfr(VL vl,P f,DCP *dcp)
 {  {
         DCP dc,dct,dcs;    DCP dc,dct,dcs;
         P c,p,t,s,pc;    P c,p,t,s,pc;
         VL mvl;    VL mvl;
   
         ptozp(f,1,(Q *)&c,&t); monomialfctr(vl,t,&p,&dc);    ptozp(f,1,(Q *)&c,&t); monomialfctr(vl,t,&p,&dc);
         if ( NUM(p) ) {    if ( NUM(p) ) {
                 *dcp = dc;      *dcp = dc;
                 return;      return;
         }    }
         mindegp(vl,p,&mvl,&s);    mindegp(vl,p,&mvl,&s);
 #if 0  #if 0
         minlcdegp(vl,p,&mvl,&s);    minlcdegp(vl,p,&mvl,&s);
         min_common_vars_in_coefp(vl,p,&mvl,&s);    min_common_vars_in_coefp(vl,p,&mvl,&s);
 #endif  #endif
         pcp(mvl,s,&p,&pc);    pcp(mvl,s,&p,&pc);
         if ( !NUM(pc) ) {    if ( !NUM(pc) ) {
                 msqfr(mvl,pc,&dcs);      msqfr(mvl,pc,&dcs);
                 if ( !dc )      if ( !dc )
                         dc = dcs;        dc = dcs;
                 else {      else {
                         for ( dct = dc; NEXT(dct); dct = NEXT(dct) );        for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
                         NEXT(dct) = dcs;        NEXT(dct) = dcs;
                 }      }
         }    }
         msqfrmain(mvl,p,&dcs);    msqfrmain(mvl,p,&dcs);
         for ( dct = dcs; dct; dct = NEXT(dct) ) {    for ( dct = dcs; dct; dct = NEXT(dct) ) {
                 reorderp(vl,mvl,COEF(dct),&t); COEF(dct) = t;      reorderp(vl,mvl,COEF(dct),&t); COEF(dct) = t;
         }    }
         if ( !dc )    if ( !dc )
                 *dcp = dcs;      *dcp = dcs;
         else {    else {
                 for ( dct = dc; NEXT(dct); dct = NEXT(dct) );      for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
                 NEXT(dct) = dcs; *dcp = dc;      NEXT(dct) = dcs; *dcp = dc;
         }    }
 }  }
   
 void usqp(P f,DCP *dcp)  void usqp(P f,DCP *dcp)
 {  {
         int index,nindex;    int index,nindex;
         P g,c,h;    P g,c,h;
         DCP dc;    DCP dc;
   
         ptozp(f,1,(Q *)&c,&h);    ptozp(f,1,(Q *)&c,&h);
         if ( SGN((Q)LC(h)) < 0 )    if ( SGN((Q)LC(h)) < 0 )
                 chsgnp(h,&g);      chsgnp(h,&g);
         else    else
                 g = h;      g = h;
         for ( index = 0, dc = 0; !dc; index = nindex )    for ( index = 0, dc = 0; !dc; index = nindex )
                 hsq(index,5,g,&nindex,&dc);      hsq(index,5,g,&nindex,&dc);
         *dcp = dc;    *dcp = dc;
 }  }
   
 void msqfrmain(VL vl,P p,DCP *dcp)  void msqfrmain(VL vl,P p,DCP *dcp)
 {  {
         int i,j;    int i,j;
         VL nvl,tvl;    VL nvl,tvl;
         VN vn,vnt,vn1;    VN vn,vnt,vn1;
         P gg,tmp,p0,g;    P gg,tmp,p0,g;
         DCP dc,dc0,dcr,dcr0;    DCP dc,dc0,dcr,dcr0;
         int nv,d,d1;    int nv,d,d1;
         int found;    int found;
         VN svn1;    VN svn1;
         P sp0;    P sp0;
         DCP sdc0;    DCP sdc0;
   
         if ( deg(VR(p),p) == 1 ) {    if ( deg(VR(p),p) == 1 ) {
                 NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;      NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                 *dcp = dc;      *dcp = dc;
                 return;      return;
         }    }
         clctv(vl,p,&nvl);    clctv(vl,p,&nvl);
         for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);    for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
         if ( nv == 1 ) {    if ( nv == 1 ) {
                 usqp(p,dcp);      usqp(p,dcp);
                 return;      return;
         }    }
 #if 0  #if 0
         if ( heusqfrchk(nvl,p) ) {    if ( heusqfrchk(nvl,p) ) {
                 NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;      NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                 *dcp = dc;      *dcp = dc;
                 return;      return;
         }    }
 #endif  #endif
         W_CALLOC(nv,struct oVN,vn);    W_CALLOC(nv,struct oVN,vn);
         W_CALLOC(nv,struct oVN,vnt);    W_CALLOC(nv,struct oVN,vnt);
         W_CALLOC(nv,struct oVN,vn1);    W_CALLOC(nv,struct oVN,vn1);
         W_CALLOC(nv,struct oVN,svn1);    W_CALLOC(nv,struct oVN,svn1);
         for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )    for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
                 vn1[i].v = vn[i].v = tvl->v;      vn1[i].v = vn[i].v = tvl->v;
         vn1[i].v = vn[i].v = 0;    vn1[i].v = vn[i].v = 0;
   
         /* determine a heuristic bound of deg(GCD(p,p')) */    /* determine a heuristic bound of deg(GCD(p,p')) */
         while ( 1 ) {    while ( 1 ) {
                 for ( i = 0; vn1[i].v; i++ )      for ( i = 0; vn1[i].v; i++ )
                         vn1[i].n = ((unsigned int)random())%256+1;        vn1[i].n = ((unsigned int)random())%256+1;
                 substvp(nvl,LC(p),vn1,&tmp);      substvp(nvl,LC(p),vn1,&tmp);
                 if ( tmp ) {      if ( tmp ) {
                         substvp(nvl,p,vn1,&p0);        substvp(nvl,p,vn1,&p0);
                         usqp(p0,&dc0);        usqp(p0,&dc0);
                         for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )        for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )
                                 if ( DEG(dc) )          if ( DEG(dc) )
                                         d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));            d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));
                         if ( d1 == 0 ) {        if ( d1 == 0 ) {
                                 /* p is squarefree */          /* p is squarefree */
                                 NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;          NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
                                 *dcp = dc;          *dcp = dc;
                                 return;          return;
                         } else {        } else {
                                 d = d1+1; /* XXX : try searching better evaluation */          d = d1+1; /* XXX : try searching better evaluation */
                                 found = 0;          found = 0;
                                 break;          break;
                         }        }
                 }      }
         }    }
   
         for  ( dcr0 = 0, g = p; ; ) {    for  ( dcr0 = 0, g = p; ; ) {
                 while ( 1 ) {      while ( 1 ) {
                         for ( i = 0, j = 0; vn[i].v; i++ )        for ( i = 0, j = 0; vn[i].v; i++ )
                                 if ( vn[i].n ) vnt[j++].v = (V)i;          if ( vn[i].n ) vnt[j++].v = (V)i;
                         vnt[j].n = 0;        vnt[j].n = 0;
   
                         mulsgn(vn,vnt,j,vn1);        mulsgn(vn,vnt,j,vn1);
                         substvp(nvl,LC(g),vn1,&tmp);        substvp(nvl,LC(g),vn1,&tmp);
                         if ( tmp ) {        if ( tmp ) {
                                 substvp(nvl,g,vn1,&p0);          substvp(nvl,g,vn1,&p0);
                                 usqp(p0,&dc0);          usqp(p0,&dc0);
                                 for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )          for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )
                                         if ( DEG(dc) )            if ( DEG(dc) )
                                                 d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));              d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));
   
                                 if ( d1 == 0 ) {          if ( d1 == 0 ) {
                                         NEWDC(dc); DEG(dc) = ONE; COEF(dc) = g; NEXT(dc) = 0;            NEWDC(dc); DEG(dc) = ONE; COEF(dc) = g; NEXT(dc) = 0;
                                         if ( !dcr0 )            if ( !dcr0 )
                                                 dcr0 = dc;              dcr0 = dc;
                                         else            else
                                                 NEXT(dcr) = dc;              NEXT(dcr) = dc;
                                         *dcp = dcr0;            *dcp = dcr0;
                                         return;            return;
                                 }          }
   
                                 if ( d < d1 )          if ( d < d1 )
                                         goto END;            goto END;
                                 if ( d > d1 ) {          if ( d > d1 ) {
                                         d = d1;            d = d1;
                                         found = 1;            found = 1;
                                         bcopy((char *)vn1,(char *)svn1,(int)(sizeof(struct oVN)*nv));            bcopy((char *)vn1,(char *)svn1,(int)(sizeof(struct oVN)*nv));
                                         sp0 = p0; sdc0 = dc0;            sp0 = p0; sdc0 = dc0;
                                         goto END;            goto END;
                                 }          }
                                 /* d1 == d */          /* d1 == d */
                                 if ( found ) {          if ( found ) {
                                         found = 0;            found = 0;
 #if 0  #if 0
                                 if ( d > d1 ) {          if ( d > d1 ) {
                                         d = d1;            d = d1;
                                                         /*} */                /*} */
 #endif  #endif
                                         msqfrmainmain(nvl,g,svn1,sp0,sdc0,&dc,&gg); g = gg;            msqfrmainmain(nvl,g,svn1,sp0,sdc0,&dc,&gg); g = gg;
                                         if ( dc ) {            if ( dc ) {
                                                 if ( !dcr0 )              if ( !dcr0 )
                                                         dcr0 = dc;                dcr0 = dc;
                                                 else              else
                                                         NEXT(dcr) = dc;                NEXT(dcr) = dc;
                                                 for ( dcr = dc; NEXT(dcr); dcr = NEXT(dcr) );              for ( dcr = dc; NEXT(dcr); dcr = NEXT(dcr) );
                                                 if ( NUM(g) ) {              if ( NUM(g) ) {
                                                         NEXT(dcr) = 0; *dcp = dcr0;                NEXT(dcr) = 0; *dcp = dcr0;
                                                         return;                return;
                                                 }              }
                                                 d = deg(VR(g),g);              d = deg(VR(g),g);
                                         }            }
                                 }          }
                         }        }
 END:  END:
                         if ( nextbin(vnt,j) )        if ( nextbin(vnt,j) )
                                 break;          break;
                 }      }
                 next(vn);      next(vn);
         }    }
 }  }
   
 void msqfrmainmain(VL vl,P p,VN vn,P p0,DCP dc0,DCP *dcp,P *pp)  void msqfrmainmain(VL vl,P p,VN vn,P p0,DCP dc0,DCP *dcp,P *pp)
 {  {
         int i,j,k,np;    int i,j,k,np;
         DCP *a;    DCP *a;
         DCP dc,dcr,dcr0,dct;    DCP dc,dcr,dcr0,dct;
         P g,t,s,u,t0,f,f0,d,d0,g0,h0,x,xx;    P g,t,s,u,t0,f,f0,d,d0,g0,h0,x,xx;
         Q q;    Q q;
         V v;    V v;
   
         for ( np = 0, dc = dc0; dc; dc = NEXT(dc), np++ );    for ( np = 0, dc = dc0; dc; dc = NEXT(dc), np++ );
         a = (DCP *)ALLOCA((np + 1)*sizeof(DCP));    a = (DCP *)ALLOCA((np + 1)*sizeof(DCP));
         for ( i = 0, dc = dc0; dc; i++, dc = NEXT(dc) )    for ( i = 0, dc = dc0; dc; i++, dc = NEXT(dc) )
                 a[np-i-1] = dc;      a[np-i-1] = dc;
   
         for ( i = 0, dcr0 = 0, f = p, f0 = p0, v = VR(p);    for ( i = 0, dcr0 = 0, f = p, f0 = p0, v = VR(p);
                 i < np; i++ ) {      i < np; i++ ) {
                 if ( (i == (np-1))&&UNIQ(DEG(a[i])) ) {      if ( (i == (np-1))&&UNIQ(DEG(a[i])) ) {
                         NEXTDC(dcr0,dcr);        NEXTDC(dcr0,dcr);
                         DEG(dcr) = DEG(a[i]);        DEG(dcr) = DEG(a[i]);
                         COEF(dcr) = f;        COEF(dcr) = f;
                         f = (P)ONE;        f = (P)ONE;
                 } else if ( (i == (np-1))&&UNIQ(DEG(DC(COEF(a[i])))) ) {      } else if ( (i == (np-1))&&UNIQ(DEG(DC(COEF(a[i])))) ) {
                         diffp(vl,f,v,&s); pcp(vl,s,&t,&u);        diffp(vl,f,v,&s); pcp(vl,s,&t,&u);
                         if ( divtpz(vl,f,t,&s) ) {        if ( divtpz(vl,f,t,&s) ) {
                                 NEXTDC(dcr0,dcr);          NEXTDC(dcr0,dcr);
                                 DEG(dcr) = DEG(a[i]);          DEG(dcr) = DEG(a[i]);
                                 COEF(dcr) = s;          COEF(dcr) = s;
                                 f = (P)ONE;          f = (P)ONE;
                         } else        } else
                                 break;          break;
                 } else {      } else {
                         for ( t = f, t0 = f0,        for ( t = f, t0 = f0,
                                 j = 0, k = QTOS(DEG(a[i]))-1; j < k; j++ ) {          j = 0, k = QTOS(DEG(a[i]))-1; j < k; j++ ) {
                                 diffp(vl,t,v,&s); t = s;          diffp(vl,t,v,&s); t = s;
                                 diffp(vl,t0,v,&s); t0 = s;          diffp(vl,t0,v,&s); t0 = s;
                         }        }
                         factorial(k,&q);        factorial(k,&q);
                         divsp(vl,t,(P)q,&s);        divsp(vl,t,(P)q,&s);
                         for ( dct = DC(s); NEXT(dct); dct = NEXT(dct) );        for ( dct = DC(s); NEXT(dct); dct = NEXT(dct) );
                         if ( DEG(dct) ) {        if ( DEG(dct) ) {
                                 MKV(VR(s),x); pwrp(vl,x,DEG(dct),&xx);          MKV(VR(s),x); pwrp(vl,x,DEG(dct),&xx);
                                 divsp(vl,s,xx,&d);          divsp(vl,s,xx,&d);
                         } else {        } else {
                                 xx = (P)ONE; d = s;          xx = (P)ONE; d = s;
                         }        }
                         divsp(vl,t0,xx,&t);        divsp(vl,t0,xx,&t);
                         divsp(vl,t,(P)q,&s);        divsp(vl,t,(P)q,&s);
                         ptozp(s,1,(Q *)&t,&d0);        ptozp(s,1,(Q *)&t,&d0);
   
                         for ( dct = DC(COEF(a[i])); NEXT(dct); dct = NEXT(dct) );        for ( dct = DC(COEF(a[i])); NEXT(dct); dct = NEXT(dct) );
                         if ( DEG(dct) )        if ( DEG(dct) )
                                 divsp(vl,COEF(a[i]),xx,&g0);          divsp(vl,COEF(a[i]),xx,&g0);
                         else {        else {
                                 xx = (P)ONE; g0 = COEF(a[i]);          xx = (P)ONE; g0 = COEF(a[i]);
                         }        }
   
                         pcp(vl,d,&t,&u); d = t;        pcp(vl,d,&t,&u); d = t;
                         ptozp(g0,1,(Q *)&u,&t); g0 = t;        ptozp(g0,1,(Q *)&u,&t); g0 = t;
   
                 {      {
                         DCP dca,dcb;        DCP dca,dcb;
   
                         fctrp(vl,LC(d),&dca);        fctrp(vl,LC(d),&dca);
                         for ( dcb = dca, u = (P)ONE; dcb; dcb = NEXT(dcb) ) {        for ( dcb = dca, u = (P)ONE; dcb; dcb = NEXT(dcb) ) {
                                 if ( NUM(COEF(dcb)) ) {          if ( NUM(COEF(dcb)) ) {
                                         mulp(vl,u,COEF(dcb),&t); u = t;            mulp(vl,u,COEF(dcb),&t); u = t;
                                 } else {          } else {
                                         Q qq;            Q qq;
                                         P tt;            P tt;
   
                                         pwrp(vl,COEF(dcb),DEG(a[i]),&s);            pwrp(vl,COEF(dcb),DEG(a[i]),&s);
                                         for ( t = LC(f), j = 0; divtpz(vl,t,s,&tt); j++, t = tt );            for ( t = LC(f), j = 0; divtpz(vl,t,s,&tt); j++, t = tt );
                                         STOQ(j,qq);            STOQ(j,qq);
                                         if ( cmpq(qq,DEG(dcb)) > 0 )            if ( cmpq(qq,DEG(dcb)) > 0 )
                                                 qq = DEG(dcb);              qq = DEG(dcb);
                                         pwrp(vl,COEF(dcb),qq,&t); mulp(vl,u,t,&s); u = s;            pwrp(vl,COEF(dcb),qq,&t); mulp(vl,u,t,&s); u = s;
                                 }          }
                         }        }
                         divsp(vl,d0,g0,&h0);        divsp(vl,d0,g0,&h0);
                 }      }
                         mfctrhen2(vl,vn,d,d0,g0,h0,u,LC(d),&s);        mfctrhen2(vl,vn,d,d0,g0,h0,u,LC(d),&s);
                         if ( s ) {        if ( s ) {
                                 mulp(vl,s,xx,&g);          mulp(vl,s,xx,&g);
                                 pwrp(vl,g,DEG(a[i]),&t);          pwrp(vl,g,DEG(a[i]),&t);
                                 if ( divtpz(vl,f,t,&s) ) {          if ( divtpz(vl,f,t,&s) ) {
                                         NEXTDC(dcr0,dcr);            NEXTDC(dcr0,dcr);
                                         DEG(dcr) = DEG(a[i]); COEF(dcr) = g;            DEG(dcr) = DEG(a[i]); COEF(dcr) = g;
                                         f = s; substvp(vl,f,vn,&f0);            f = s; substvp(vl,f,vn,&f0);
                                 } else          } else
                                         break;            break;
                         } else        } else
                                 break;          break;
                 }      }
         }    }
         *pp = f;    *pp = f;
         if ( dcr0 )    if ( dcr0 )
                 NEXT(dcr) = 0;      NEXT(dcr) = 0;
         *dcp = dcr0;    *dcp = dcr0;
 }  }
   
 void mfctrhen2(VL vl,VN vn,P f,P f0,P g0,P h0,P lcg,P lch,P *gp)  void mfctrhen2(VL vl,VN vn,P f,P f0,P g0,P h0,P lcg,P lch,P *gp)
 {  {
         V v;    V v;
         P f1,lc,lc0,lcg0,lch0;    P f1,lc,lc0,lcg0,lch0;
         P m,ff0,gg0,hh0,gk,hk,ggg,gggr,hhh,ak,bk,tmp;    P m,ff0,gg0,hh0,gk,hk,ggg,gggr,hhh,ak,bk,tmp;
         Q bb,qk,s;    Q bb,qk,s;
         Q cbd;    Q cbd;
         int dbd;    int dbd;
         int d;    int d;
   
         if ( NUM(g0) ) {    if ( NUM(g0) ) {
                 *gp = (P)ONE;      *gp = (P)ONE;
                 return;      return;
         }    }
   
         v = VR(f); d = deg(v,f);    v = VR(f); d = deg(v,f);
         if ( d == deg(v,g0) ) {    if ( d == deg(v,g0) ) {
                 pcp(vl,f,gp,&tmp);      pcp(vl,f,gp,&tmp);
                 return;      return;
         }    }
   
         mulp(vl,lcg,lch,&lc);    mulp(vl,lcg,lch,&lc);
         if ( !divtpz(vl,lc,LC(f),(P *)&s) ) {    if ( !divtpz(vl,lc,LC(f),(P *)&s) ) {
                 *gp = 0;      *gp = 0;
                 return;      return;
         }    }
         mulp(vl,(P)s,f,&f1);    mulp(vl,(P)s,f,&f1);
         dbd = dbound(VR(f1),f1) + 1; cbound(vl,f1,&cbd);    dbd = dbound(VR(f1),f1) + 1; cbound(vl,f1,&cbd);
   
         substvp(vl,lc,vn,&lc0);    substvp(vl,lc,vn,&lc0);
         divq((Q)lc0,(Q)LC(f0),(Q *)&m); mulp(vl,f0,m,&ff0);    divq((Q)lc0,(Q)LC(f0),(Q *)&m); mulp(vl,f0,m,&ff0);
         substvp(vl,lcg,vn,&lcg0);    substvp(vl,lcg,vn,&lcg0);
         divq((Q)lcg0,(Q)LC(g0),(Q *)&m); mulp(vl,g0,m,&gg0);    divq((Q)lcg0,(Q)LC(g0),(Q *)&m); mulp(vl,g0,m,&gg0);
         substvp(vl,lch,vn,&lch0);    substvp(vl,lch,vn,&lch0);
         divq((Q)lch0,(Q)LC(h0),(Q *)&m); mulp(vl,h0,m,&hh0);    divq((Q)lch0,(Q)LC(h0),(Q *)&m); mulp(vl,h0,m,&hh0);
         addq(cbd,cbd,&bb);    addq(cbd,cbd,&bb);
         henzq1(gg0,hh0,bb,&bk,&ak,&qk); gk = gg0; hk = hh0;    henzq1(gg0,hh0,bb,&bk,&ak,&qk); gk = gg0; hk = hh0;
         henmv(vl,vn,f1,gk,hk,ak,bk,    henmv(vl,vn,f1,gk,hk,ak,bk,
                 lcg,lch,lcg0,lch0,qk,dbd,&ggg,&hhh);      lcg,lch,lcg0,lch0,qk,dbd,&ggg,&hhh);
   
         if ( divtpz(vl,f1,ggg,&gggr) )    if ( divtpz(vl,f1,ggg,&gggr) )
                 pcp(vl,ggg,gp,&tmp);      pcp(vl,ggg,gp,&tmp);
         else    else
                 *gp = 0;      *gp = 0;
 }  }
   
 int sqfrchk(P p)  int sqfrchk(P p)
 {  {
         Q c;    Q c;
         P f;    P f;
         DCP dc;    DCP dc;
   
         ptozp(p,SGN((Q)UCOEF(p)),&c,&f); usqp(f,&dc);    ptozp(p,SGN((Q)UCOEF(p)),&c,&f); usqp(f,&dc);
         if ( NEXT(dc) || !UNIQ(DEG(dc)) )    if ( NEXT(dc) || !UNIQ(DEG(dc)) )
                 return ( 0 );      return ( 0 );
         else    else
                 return ( 1 );      return ( 1 );
 }  }
   
 int cycchk(P p)  int cycchk(P p)
 {  {
         Q c;    Q c;
         P f;    P f;
   
         ptozp(p,SGN((Q)UCOEF(p)),&c,&f);    ptozp(p,SGN((Q)UCOEF(p)),&c,&f);
         if ( iscycp(f) || iscycm(f) )    if ( iscycp(f) || iscycm(f) )
                 return 0;      return 0;
         else    else
                 return 1;      return 1;
 }  }
   
 int zerovpchk(VL vl,P p,VN vn)  int zerovpchk(VL vl,P p,VN vn)
 {  {
         P t;    P t;
   
         substvp(vl,p,vn,&t);    substvp(vl,p,vn,&t);
         if ( t )    if ( t )
                 return ( 0 );      return ( 0 );
         else    else
                 return ( 1 );      return ( 1 );
 }  }
   
 int valideval(VL vl,DCP dc,VN vn)  int valideval(VL vl,DCP dc,VN vn)
 {  {
         DCP dct;    DCP dct;
         Q *a;    Q *a;
         int i,j,n;    int i,j,n;
         N q,r;    N q,r;
   
         for ( dct = NEXT(dc), n = 0; dct; dct = NEXT(dct), n++ );    for ( dct = NEXT(dc), n = 0; dct; dct = NEXT(dct), n++ );
         W_CALLOC(n,Q,a);    W_CALLOC(n,Q,a);
         for ( dct = NEXT(dc), i = 0; dct; dct = NEXT(dct), i++ ) {    for ( dct = NEXT(dc), i = 0; dct; dct = NEXT(dct), i++ ) {
                 substvp(vl,COEF(dct),vn,(P *)&a[i]);      substvp(vl,COEF(dct),vn,(P *)&a[i]);
                 if ( !a[i] )      if ( !a[i] )
                         return ( 0 );        return ( 0 );
   
                 for ( j = 0; j < i; j++ ) {      for ( j = 0; j < i; j++ ) {
                         divn(NM(a[j]),NM(a[i]),&q,&r);        divn(NM(a[j]),NM(a[i]),&q,&r);
                         if ( !r )        if ( !r )
                                 return ( 0 );          return ( 0 );
                         divn(NM(a[i]),NM(a[j]),&q,&r);        divn(NM(a[i]),NM(a[j]),&q,&r);
                         if ( !r )        if ( !r )
                                 return ( 0 );          return ( 0 );
                 }      }
         }    }
         return ( 1 );    return ( 1 );
 }  }
   
 void estimatelc(VL vl,Q c,DCP dc,VN vn,P *lcp)  void estimatelc(VL vl,Q c,DCP dc,VN vn,P *lcp)
 {  {
         int i;    int i;
         DCP dct;    DCP dct;
         P r,s,t;    P r,s,t;
         Q c0,c1,c2;    Q c0,c1,c2;
   
         for ( dct = dc, r = (P)ONE; dct; dct = NEXT(dct) ) {    for ( dct = dc, r = (P)ONE; dct; dct = NEXT(dct) ) {
                 if ( NUM(COEF(dct)) ) {      if ( NUM(COEF(dct)) ) {
                         mulp(vl,r,COEF(dct),&s); r = s;        mulp(vl,r,COEF(dct),&s); r = s;
                 } else {      } else {
                         substvp(vl,COEF(dct),vn,(P *)&c0);        substvp(vl,COEF(dct),vn,(P *)&c0);
                         for ( i = 0, c1 = c; i < (int)QTOS(DEG(dct)); i++ ) {        for ( i = 0, c1 = c; i < (int)QTOS(DEG(dct)); i++ ) {
                                 divq(c1,c0,&c2);          divq(c1,c0,&c2);
                                 if ( !INT(c2) )          if ( !INT(c2) )
                                         break;            break;
                                 else          else
                                         c1 = c2;            c1 = c2;
                         }        }
                         if ( i ) {        if ( i ) {
                                 STOQ(i,c1);          STOQ(i,c1);
                                 pwrp(vl,COEF(dct),c1,&s); mulp(vl,r,s,&t); r = t;          pwrp(vl,COEF(dct),c1,&s); mulp(vl,r,s,&t); r = t;
                         }        }
                 }      }
         }    }
         *lcp = r;    *lcp = r;
 }  }
   
 void monomialfctr(VL vl,P p,P *pr,DCP *dcp)  void monomialfctr(VL vl,P p,P *pr,DCP *dcp)
 {  {
         VL nvl,avl;    VL nvl,avl;
         Q d;    Q d;
         P f,t,s;    P f,t,s;
         DCP dc0,dc;    DCP dc0,dc;
   
         clctv(vl,p,&nvl);    clctv(vl,p,&nvl);
         for ( dc0 = 0, avl = nvl, f = p; avl; avl = NEXT(avl) ) {    for ( dc0 = 0, avl = nvl, f = p; avl; avl = NEXT(avl) ) {
                 getmindeg(avl->v,f,&d);      getmindeg(avl->v,f,&d);
                 if ( d ) {      if ( d ) {
                         MKV(avl->v,t); NEXTDC(dc0,dc); DEG(dc) = d; COEF(dc) = t;        MKV(avl->v,t); NEXTDC(dc0,dc); DEG(dc) = d; COEF(dc) = t;
                         pwrp(vl,t,d,&s); divsp(vl,f,s,&t); f = t;        pwrp(vl,t,d,&s); divsp(vl,f,s,&t); f = t;
                 }      }
         }    }
         if ( dc0 )    if ( dc0 )
                 NEXT(dc) = 0;      NEXT(dc) = 0;
         *pr = f; *dcp = dc0;    *pr = f; *dcp = dc0;
 }  }
   
 void afctr(VL vl,P p0,P p,DCP *dcp)  void afctr(VL vl,P p0,P p,DCP *dcp)
 {  {
         DCP dc,dc0,dcr,dct,dcs;    DCP dc,dc0,dcr,dct,dcs;
         P t;    P t;
         VL nvl;    VL nvl;
   
         if ( VR(p) == VR(p0) ) {    if ( VR(p) == VR(p0) ) {
                 NEWDC(dc);      NEWDC(dc);
                 DEG(dc) = ONE;      DEG(dc) = ONE;
                 COEF(dc) = p;      COEF(dc) = p;
                 NEXT(dc) = 0;      NEXT(dc) = 0;
                 *dcp = dc;      *dcp = dc;
                 return;      return;
         }    }
   
         clctv(vl,p,&nvl);    clctv(vl,p,&nvl);
         if ( !NEXT(nvl) )    if ( !NEXT(nvl) )
                 ufctr(p,1,&dc);      ufctr(p,1,&dc);
         else {    else {
                 sqa(vl,p0,p,&dc);      sqa(vl,p0,p,&dc);
                 for ( dct = dc; dct; dct = NEXT(dct) ) {      for ( dct = dc; dct; dct = NEXT(dct) ) {
                         pmonic(vl,p0,COEF(dct),&t); COEF(dct) = t;        pmonic(vl,p0,COEF(dct),&t); COEF(dct) = t;
                 }      }
         }    }
         if ( NUM(COEF(dc)) )    if ( NUM(COEF(dc)) )
                 dcr = NEXT(dc);      dcr = NEXT(dc);
         else    else
                 dcr = dc;      dcr = dc;
         for ( dc0 = 0; dcr; dcr = NEXT(dcr) ) {    for ( dc0 = 0; dcr; dcr = NEXT(dcr) ) {
                 afctrmain(vl,p0,COEF(dcr),1,&dcs);      afctrmain(vl,p0,COEF(dcr),1,&dcs);
   
                 for ( dct = dcs; dct; dct = NEXT(dct) )      for ( dct = dcs; dct; dct = NEXT(dct) )
                         DEG(dct) = DEG(dcr);        DEG(dct) = DEG(dcr);
                 if ( !dc0 )      if ( !dc0 )
                         dc0 = dcs;        dc0 = dcs;
                 else {      else {
                         for ( dct = dc0; NEXT(dct); dct = NEXT(dct) );        for ( dct = dc0; NEXT(dct); dct = NEXT(dct) );
                         NEXT(dct) = dcs;        NEXT(dct) = dcs;
                 }      }
         }    }
         *dcp = dc0;    *dcp = dc0;
 }  }
   
 void afctrmain(VL vl,P p0,P p,int init,DCP *dcp)  void afctrmain(VL vl,P p0,P p,int init,DCP *dcp)
 {  {
         P x,y,s,m,a,t,u,pt,pt1,res,g;    P x,y,s,m,a,t,u,pt,pt1,res,g;
         Q q;    Q q;
         DCP dc,dc0,dcsq0,dcr0,dcr,dct,dcs;    DCP dc,dc0,dcsq0,dcr0,dcr,dct,dcs;
         V v,v0;    V v,v0;
   
         if ( !cmpq(DEG(DC(p)),ONE) ) {    if ( !cmpq(DEG(DC(p)),ONE) ) {
                 NEWDC(dc); DEG(dc) = ONE; NEXT(dc) = 0;      NEWDC(dc); DEG(dc) = ONE; NEXT(dc) = 0;
                 pmonic(vl,p0,p,&COEF(dc)); *dcp = dc;      pmonic(vl,p0,p,&COEF(dc)); *dcp = dc;
                 return;      return;
         }    }
   
         v = VR(p); MKV(v,x);    v = VR(p); MKV(v,x);
         v0 = VR(p0); MKV(v0,y);    v0 = VR(p0); MKV(v0,y);
         STOQ(init,q),s = (P)q;    STOQ(init,q),s = (P)q;
         mulp(vl,s,y,&m); subp(vl,x,m,&t); addp(vl,x,m,&a);    mulp(vl,s,y,&m); subp(vl,x,m,&t); addp(vl,x,m,&a);
         substp(vl,p,v,t,&pt);    substp(vl,p,v,t,&pt);
         remsdcp(vl,pt,p0,&pt1);    remsdcp(vl,pt,p0,&pt1);
   
 /*  /*
         if ( ( deg(v0,p0) <= 3 ) || ( TCPQ(p0) <= 2 ) )    if ( ( deg(v0,p0) <= 3 ) || ( TCPQ(p0) <= 2 ) )
                 resultp(vl,v0,p0,pt1,&res);      resultp(vl,v0,p0,pt1,&res);
         else    else
                 srcrnorm(vl,v0,pt1,p0,&res);      srcrnorm(vl,v0,pt1,p0,&res);
 */  */
 #if 0  #if 0
         srcr(vl,v0,pt1,p0,&res);    srcr(vl,v0,pt1,p0,&res);
 #endif  #endif
         resultp(vl,v0,p0,pt1,&res);    resultp(vl,v0,p0,pt1,&res);
         usqp(res,&dcsq0);    usqp(res,&dcsq0);
         for ( dc0 = 0, dct = dcsq0; dct; dct = NEXT(dct) ) {    for ( dc0 = 0, dct = dcsq0; dct; dct = NEXT(dct) ) {
                 if ( UNIQ(DEG(dct)) )      if ( UNIQ(DEG(dct)) )
                         ufctr(COEF(dct),deg(v0,p0),&dcs);        ufctr(COEF(dct),deg(v0,p0),&dcs);
                 else      else
                         ufctr(COEF(dct),1,&dcs);        ufctr(COEF(dct),1,&dcs);
                 for ( dcr = dcs; dcr; dcr = NEXT(dcr) )      for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
                         DEG(dcr) = DEG(dct);        DEG(dcr) = DEG(dct);
                 if ( !dc0 ) {      if ( !dc0 ) {
                         dc0 = NEXT(dcs);        dc0 = NEXT(dcs);
                         dc = dc0;        dc = dc0;
                 } else {      } else {
                         for ( ; NEXT(dc); dc = NEXT(dc) );        for ( ; NEXT(dc); dc = NEXT(dc) );
                         NEXT(dc) = NEXT(dcs);        NEXT(dc) = NEXT(dcs);
                 }      }
         }    }
         sortfs(&dc0);    sortfs(&dc0);
   
         for ( g = p, dcr = dcr0 = 0, dc = dc0; dc; dc = NEXT(dc) ) {    for ( g = p, dcr = dcr0 = 0, dc = dc0; dc; dc = NEXT(dc) ) {
                 if ( !UNIQ(DEG(dc)) ) {      if ( !UNIQ(DEG(dc)) ) {
                         substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s);        substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s);
                         gcda(vl,p0,s,g,&u);        gcda(vl,p0,s,g,&u);
                         if ( !NUM(u) && (VR(u) == v)) {        if ( !NUM(u) && (VR(u) == v)) {
                                 afctrmain(vl,p0,u,init+1,&dct);          afctrmain(vl,p0,u,init+1,&dct);
                                 for ( dcs = dct, t = (P)ONE; dcs; dcs = NEXT(dcs) ) {          for ( dcs = dct, t = (P)ONE; dcs; dcs = NEXT(dcs) ) {
                                         mulp(vl,t,COEF(dcs),&s); remsdcp(vl,s,p0,&t);            mulp(vl,t,COEF(dcs),&s); remsdcp(vl,s,p0,&t);
                                 }          }
                                 pdiva(vl,p0,g,t,&s); g = s;          pdiva(vl,p0,g,t,&s); g = s;
                                 if ( !dcr0 )          if ( !dcr0 )
                                         dcr0 = dct;            dcr0 = dct;
                                 else          else
                                         NEXT(dcr) = dct;            NEXT(dcr) = dct;
                                 for ( dcr = dct; NEXT(dcr); dcr = NEXT(dcr) );          for ( dcr = dct; NEXT(dcr); dcr = NEXT(dcr) );
                         }        }
                 } else {      } else {
                         substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s);        substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s);
                         gcda(vl,p0,s,g,&u);        gcda(vl,p0,s,g,&u);
                         if ( !NUM(u) && (VR(u) == v)) {        if ( !NUM(u) && (VR(u) == v)) {
                                 NEXTDC(dcr0,dcr);          NEXTDC(dcr0,dcr);
                                 DEG(dcr) = ONE;          DEG(dcr) = ONE;
                                 COEF(dcr) = u;          COEF(dcr) = u;
                                 pdiva(vl,p0,g,u,&t); g = t;          pdiva(vl,p0,g,u,&t); g = t;
                         }        }
                 }      }
         }    }
         if ( dcr0 )    if ( dcr0 )
                 NEXT(dcr) = 0;      NEXT(dcr) = 0;
         *dcp = dcr0;    *dcp = dcr0;
 }  }
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.13
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>