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

Diff for /OpenXM_contrib2/asir2000/engine/EZ.c between version 1.3 and 1.4

version 1.3, 2000/08/22 05:04:03 version 1.4, 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/EZ.c,v 1.2 2000/08/21 08:31:24 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/EZ.c,v 1.3 2000/08/22 05:04:03 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
   
Line 54  VL vl;
Line 54  VL vl;
 P p1,p2;  P p1,p2;
 P *pr;  P *pr;
 {  {
         P f1,f2;    P f1,f2;
         Q c;    Q c;
   
         if ( !p1 )    if ( !p1 )
                 if ( !p2 )      if ( !p2 )
                         *pr = 0;        *pr = 0;
                 else      else
                         *pr = p2;        *pr = p2;
         else if ( !p2 )    else if ( !p2 )
                 *pr = p1;      *pr = p1;
         else {    else {
                 if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )      if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
                         error("ezgcdp : invalid argument");        error("ezgcdp : invalid argument");
                 ptozp(p1,1,&c,&f1);      ptozp(p1,1,&c,&f1);
                 ptozp(p2,1,&c,&f2);      ptozp(p2,1,&c,&f2);
                 ezgcdpz(vl,f1,f2,pr);      ezgcdpz(vl,f1,f2,pr);
         }    }
 }  }
   
 /*  /*
Line 82  void ezgcdpz(vl,p1,p2,pr)
Line 82  void ezgcdpz(vl,p1,p2,pr)
 VL vl;  VL vl;
 P p1,p2,*pr;  P p1,p2,*pr;
 {  {
         P f1,f2,t,s,g,gcd;    P f1,f2,t,s,g,gcd;
         P fc1,fc2,cont;    P fc1,fc2,cont;
         VL tvl,svl;    VL tvl,svl;
         V v1,v2;    V v1,v2;
         DCP dc,dc1,dc2;    DCP dc,dc1,dc2;
         N cn;    N cn;
         Q cq;    Q cq;
         Q c1,c2,c;    Q c1,c2,c;
         P g1,g2;    P g1,g2;
         P mgcd;    P mgcd;
         extern int nez;    extern int nez;
         P pm[2];    P pm[2];
   
         if ( nez ) {    if ( nez ) {
                 pm[0] = p1; pm[1] = p2; nezgcdnpz(vl,pm,2,pr); return;      pm[0] = p1; pm[1] = p2; nezgcdnpz(vl,pm,2,pr); return;
         }    }
         monomialfctr(vl,p1,&f1,&dc1); p1 = f1;    monomialfctr(vl,p1,&f1,&dc1); p1 = f1;
         monomialfctr(vl,p2,&f2,&dc2); p2 = f2;    monomialfctr(vl,p2,&f2,&dc2); p2 = f2;
         for ( mgcd = (P)ONE; dc1; dc1 = NEXT(dc1) )    for ( mgcd = (P)ONE; dc1; dc1 = NEXT(dc1) )
                 for ( v1 = VR(dc1->c), dc = dc2; dc; dc = NEXT(dc) )      for ( v1 = VR(dc1->c), dc = dc2; dc; dc = NEXT(dc) )
                         if ( v1 == VR(dc->c) ) {        if ( v1 == VR(dc->c) ) {
                                 pwrp(vl,dc->c,cmpq(dc1->d,dc->d)>=0?dc->d:dc1->d,&t);          pwrp(vl,dc->c,cmpq(dc1->d,dc->d)>=0?dc->d:dc1->d,&t);
                                 mulp(vl,mgcd,t,&s); mgcd = s; break;          mulp(vl,mgcd,t,&s); mgcd = s; break;
                         }        }
   
         if ( NUM(p1) ) {    if ( NUM(p1) ) {
                 if ( NUM(p2) ) {      if ( NUM(p2) ) {
                         gcdn(NM((Q)p1),NM((Q)p2),&cn); NTOQ(cn,1,c);        gcdn(NM((Q)p1),NM((Q)p2),&cn); NTOQ(cn,1,c);
                 } else {      } else {
                         ptozp(p2,1,&c2,&t);        ptozp(p2,1,&c2,&t);
                         gcdn(NM((Q)p1),NM(c2),&cn); NTOQ(cn,1,c);        gcdn(NM((Q)p1),NM(c2),&cn); NTOQ(cn,1,c);
                 }      }
                 mulp(vl,(P)c,mgcd,pr);      mulp(vl,(P)c,mgcd,pr);
                 return;      return;
         } else if ( NUM(p2) ) {    } else if ( NUM(p2) ) {
                 ptozp(p1,1,&c1,&t);      ptozp(p1,1,&c1,&t);
                 gcdn(NM(c1),NM((Q)p2),&cn); NTOQ(cn,1,c);      gcdn(NM(c1),NM((Q)p2),&cn); NTOQ(cn,1,c);
                 mulp(vl,(P)c,mgcd,pr);      mulp(vl,(P)c,mgcd,pr);
                 return;      return;
         }    }
   
         ptozp(p1,1,&c1,&g1); ptozp(p2,1,&c2,&g2);    ptozp(p1,1,&c1,&g1); ptozp(p2,1,&c2,&g2);
         gcdn(NM(c1),NM(c2),&cn); NTOQ(cn,1,cq);    gcdn(NM(c1),NM(c2),&cn); NTOQ(cn,1,cq);
   
         if ( ( v1 = VR(f1) ) != ( v2 = VR(f2) ) ) {    if ( ( v1 = VR(f1) ) != ( v2 = VR(f2) ) ) {
                 while ( v1 != VR(vl) && v2 != VR(vl) )      while ( v1 != VR(vl) && v2 != VR(vl) )
                         vl = NEXT(vl);        vl = NEXT(vl);
                 if ( v1 == VR(vl) ) {      if ( v1 == VR(vl) ) {
                         pcp(vl,f1,&g,&fc1);        pcp(vl,f1,&g,&fc1);
                         ezgcdpz(vl,fc1,f2,&t);        ezgcdpz(vl,fc1,f2,&t);
                 } else {      } else {
                         pcp(vl,f2,&g,&fc2);        pcp(vl,f2,&g,&fc2);
                         ezgcdpz(vl,fc2,f1,&t);        ezgcdpz(vl,fc2,f1,&t);
                 }      }
                 mulp(vl,t,(P)cq,&s); mulp(vl,s,mgcd,pr);      mulp(vl,t,(P)cq,&s); mulp(vl,s,mgcd,pr);
                 return;      return;
         }    }
   
         pcp(vl,f1,&g1,&fc1); pcp(vl,f2,&g2,&fc2);    pcp(vl,f1,&g1,&fc1); pcp(vl,f2,&g2,&fc2);
         ezgcdpz(vl,fc1,fc2,&cont);    ezgcdpz(vl,fc1,fc2,&cont);
         clctv(vl,g1,&tvl); clctv(vl,g2,&svl);    clctv(vl,g1,&tvl); clctv(vl,g2,&svl);
         if ( !NEXT(tvl) && !NEXT(svl) ) {    if ( !NEXT(tvl) && !NEXT(svl) ) {
                 uezgcdpz(vl,g1,g2,&t);      uezgcdpz(vl,g1,g2,&t);
                 mulp(vl,t,cont,&s); mulp(vl,s,(P)cq,&t); mulp(vl,t,mgcd,pr);      mulp(vl,t,cont,&s); mulp(vl,s,(P)cq,&t); mulp(vl,t,mgcd,pr);
                 return;      return;
         }    }
   
         if ( homdeg(g1) > homdeg(g2) ) {    if ( homdeg(g1) > homdeg(g2) ) {
                 t = g1; g1 = g2; g2 = t;      t = g1; g1 = g2; g2 = t;
         }    }
         sqfrp(vl,g1,&dc);    sqfrp(vl,g1,&dc);
         for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) {    for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) {
                 if ( NUM(COEF(dc)) )      if ( NUM(COEF(dc)) )
                         continue;        continue;
                 ezgcdpp(vl,dc,g,&t);      ezgcdpp(vl,dc,g,&t);
                 if ( NUM(t) )      if ( NUM(t) )
                         continue;        continue;
                 mulp(vl,gcd,t,&s); gcd = s;      mulp(vl,gcd,t,&s); gcd = s;
                 divsp(vl,g,t,&s); g = s;      divsp(vl,g,t,&s); g = s;
         }    }
         mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,&s); mulp(vl,s,mgcd,pr);    mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,&s); mulp(vl,s,mgcd,pr);
 }  }
   
 void uezgcdpz(vl,p1,p2,pr)  void uezgcdpz(vl,p1,p2,pr)
 VL vl;  VL vl;
 P p1,p2,*pr;  P p1,p2,*pr;
 {  {
         P g1,g2,gcd,t,s,g;    P g1,g2,gcd,t,s,g;
         DCP dc;    DCP dc;
         N cn;    N cn;
         Q c1,c2,cq;    Q c1,c2,cq;
         P f1,f2;    P f1,f2;
   
         if ( NUM(p1) ) {    if ( NUM(p1) ) {
                 if ( NUM(p2) ) {      if ( NUM(p2) ) {
                         gcdn(NM((Q)p1),NM((Q)p2),&cn); NTOQ(cn,1,cq),*pr = (P)cq;        gcdn(NM((Q)p1),NM((Q)p2),&cn); NTOQ(cn,1,cq),*pr = (P)cq;
                         return;        return;
                 } else {      } else {
                         ptozp(p2,1,&c2,&t);        ptozp(p2,1,&c2,&t);
                         gcdn(NM((Q)p1),NM(c2),&cn); NTOQ(cn,1,cq),*pr = (P)cq;        gcdn(NM((Q)p1),NM(c2),&cn); NTOQ(cn,1,cq),*pr = (P)cq;
                         return;        return;
                 }      }
         } else if ( NUM(p2) ) {    } else if ( NUM(p2) ) {
                 ptozp(p1,1,&c1,&t);      ptozp(p1,1,&c1,&t);
                 gcdn(NM(c1),NM((Q)p2),&cn); NTOQ(cn,1,cq),*pr = (P)cq;      gcdn(NM(c1),NM((Q)p2),&cn); NTOQ(cn,1,cq),*pr = (P)cq;
                 return;      return;
         }    }
   
         ptozp(p1,1,&c1,&f1); ptozp(p2,1,&c2,&f2);    ptozp(p1,1,&c1,&f1); ptozp(p2,1,&c2,&f2);
         gcdn(NM(c1),NM(c2),&cn); NTOQ(cn,1,cq);    gcdn(NM(c1),NM(c2),&cn); NTOQ(cn,1,cq);
         if ( UDEG(f1) > UDEG(f2) ) {    if ( UDEG(f1) > UDEG(f2) ) {
                 g1 = f2; g2 = f1;      g1 = f2; g2 = f1;
         } else {    } else {
                 g1 = f1; g2 = f2;      g1 = f1; g2 = f2;
         }    }
         usqp(g1,&dc);    usqp(g1,&dc);
         for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) {    for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) {
                 if ( NUM(COEF(dc)) )      if ( NUM(COEF(dc)) )
                         continue;        continue;
                 uezgcdpp(dc,g,&t);      uezgcdpp(dc,g,&t);
                 if ( NUM(t) )      if ( NUM(t) )
                         continue;        continue;
                 mulp(CO,gcd,t,&s); gcd = s;      mulp(CO,gcd,t,&s); gcd = s;
                 divsp(CO,g,t,&s); g = s;      divsp(CO,g,t,&s); g = s;
         }    }
         mulp(vl,gcd,(P)cq,pr);    mulp(vl,gcd,(P)cq,pr);
 }  }
   
 /*  /*
         dc : dc pair c^d    dc : dc pair c^d
         r = gcd(c^d,f)    r = gcd(c^d,f)
 */  */
   
 void ezgcdpp(vl,dc,f,r)  void ezgcdpp(vl,dc,f,r)
Line 221  DCP dc;
Line 221  DCP dc;
 P f;  P f;
 P *r;  P *r;
 {  {
         int k;    int k;
         P g,h,t,s,gcd;    P g,h,t,s,gcd;
   
         if ( NUM(f) ) {    if ( NUM(f) ) {
                 *r = (P)ONE;      *r = (P)ONE;
                 return;      return;
         }    }
         k = QTOS(DEG(dc)) - 1;    k = QTOS(DEG(dc)) - 1;
 /*      ezgcd1p(vl,COEF(dc),f,&gcd); */  /*  ezgcd1p(vl,COEF(dc),f,&gcd); */
         ezgcdnp(vl,COEF(dc),&f,1,&gcd);    ezgcdnp(vl,COEF(dc),&f,1,&gcd);
         g = gcd; divsp(vl,f,gcd,&h);    g = gcd; divsp(vl,f,gcd,&h);
         for ( ; k > 0; k-- ) {    for ( ; k > 0; k-- ) {
 /*              ezgcd1p(vl,g,h,&t); */  /*    ezgcd1p(vl,g,h,&t); */
                 ezgcdnp(vl,g,&h,1,&t);      ezgcdnp(vl,g,&h,1,&t);
                 if ( NUM(t) )      if ( NUM(t) )
                         break;        break;
                 mulp(vl,gcd,t,&s); gcd = s;      mulp(vl,gcd,t,&s); gcd = s;
                 divsp(vl,h,t,&s); h = s;      divsp(vl,h,t,&s); h = s;
         }    }
         *r = gcd;    *r = gcd;
 }  }
   
 void ezgcd1p(vl,p0,p1,gcdp)  void ezgcd1p(vl,p0,p1,gcdp)
Line 248  VL vl;
Line 248  VL vl;
 P p0,p1;  P p0,p1;
 P *gcdp;  P *gcdp;
 {  {
         VL nvl,tvl,vl0,vl1,vl2;    VL nvl,tvl,vl0,vl1,vl2;
         P f0,f1,t;    P f0,f1,t;
   
         clctv(vl,p0,&vl0); clctv(vl,p1,&vl1); mergev(vl,vl0,vl1,&vl2);    clctv(vl,p0,&vl0); clctv(vl,p1,&vl1); mergev(vl,vl0,vl1,&vl2);
         minchdegp(vl,p0,&tvl,&t); reordvar(vl2,tvl->v,&nvl);    minchdegp(vl,p0,&tvl,&t); reordvar(vl2,tvl->v,&nvl);
         reorderp(nvl,vl,p0,&f0); reorderp(nvl,vl,p1,&f1);    reorderp(nvl,vl,p0,&f0); reorderp(nvl,vl,p1,&f1);
         ezgcdnp(nvl,f0,&f1,1,&t);    ezgcdnp(nvl,f0,&f1,1,&t);
         reorderp(vl,nvl,t,gcdp);    reorderp(vl,nvl,t,gcdp);
 }  }
   
 void uezgcdpp(dc,f,r)  void uezgcdpp(dc,f,r)
Line 263  DCP dc;
Line 263  DCP dc;
 P f;  P f;
 P *r;  P *r;
 {  {
         int k;    int k;
         P g,h,t,s,gcd;    P g,h,t,s,gcd;
   
         if ( NUM(f) ) {    if ( NUM(f) ) {
                 *r = (P)ONE;      *r = (P)ONE;
                 return;      return;
         }    }
   
         k = QTOS(DEG(dc)) - 1;    k = QTOS(DEG(dc)) - 1;
         uezgcd1p(COEF(dc),f,&gcd);    uezgcd1p(COEF(dc),f,&gcd);
         g = gcd; divsp(CO,f,gcd,&h);    g = gcd; divsp(CO,f,gcd,&h);
         for ( ; k > 0; k-- ) {    for ( ; k > 0; k-- ) {
                 uezgcd1p(g,h,&t);      uezgcd1p(g,h,&t);
                 if ( NUM(t) )      if ( NUM(t) )
                         break;        break;
                 mulp(CO,gcd,t,&s); gcd = s;      mulp(CO,gcd,t,&s); gcd = s;
                 divsp(CO,h,t,&s); h = s;      divsp(CO,h,t,&s); h = s;
         }    }
         *r = gcd;    *r = gcd;
 }  }
   
 /*  /*
Line 297  VL vl;
Line 297  VL vl;
 int m;  int m;
 P p0,*ps,*pr;  P p0,*ps,*pr;
 {  {
         /* variables */    /* variables */
         P p00,g,h,g0,h0,a0,b0;    P p00,g,h,g0,h0,a0,b0;
         P lp0,lgp0,lp00,lg,lg0,cbd,tq,t;    P lp0,lgp0,lp00,lg,lg0,cbd,tq,t;
         P *tps;    P *tps;
         VL nvl1,nvl2,nvl,tvl;    VL nvl1,nvl2,nvl,tvl;
         V v;    V v;
         int i,j,k,d0,dg,dg0,dmin;    int i,j,k,d0,dg,dg0,dmin;
         VN vn0,vn1,vnt;    VN vn0,vn1,vnt;
         int nv,flag;    int nv,flag;
   
         /* set-up */    /* set-up */
         if ( NUM(p0) ) {    if ( NUM(p0) ) {
                 *pr = (P) ONE;      *pr = (P) ONE;
                 return;      return;
         }    }
         for ( v = VR(p0), i = 0; i < m; i++ )    for ( v = VR(p0), i = 0; i < m; i++ )
                 if ( NUM(ps[i]) || (v != VR(ps[i])) ) {      if ( NUM(ps[i]) || (v != VR(ps[i])) ) {
                         *pr = (P)ONE;        *pr = (P)ONE;
                         return;        return;
                 }      }
         tps = (P *) ALLOCA(m*sizeof(P));    tps = (P *) ALLOCA(m*sizeof(P));
         for ( i = 0; i < m; i++ )    for ( i = 0; i < m; i++ )
                 tps[i] = ps[i];      tps[i] = ps[i];
         sortplist(tps,m);    sortplist(tps,m);
         /* deg(tps[0]) <= deg(tps[1]) <= ... */    /* deg(tps[0]) <= deg(tps[1]) <= ... */
   
         if ( !cmpq(DEG(DC(p0)),ONE) ) {    if ( !cmpq(DEG(DC(p0)),ONE) ) {
                 if ( divcheck(vl,tps,m,(P)ONE,p0) )      if ( divcheck(vl,tps,m,(P)ONE,p0) )
                         *pr = p0;        *pr = p0;
                 else      else
                         *pr = (P)ONE;        *pr = (P)ONE;
                 return;      return;
         }    }
   
         lp0 = LC(p0); lg = lp0; dmin = d0 = deg(v,p0);    lp0 = LC(p0); lg = lp0; dmin = d0 = deg(v,p0);
         clctv(vl,p0,&nvl);    clctv(vl,p0,&nvl);
         for ( i = 0; i < m; i++ ) {    for ( i = 0; i < m; i++ ) {
                 clctv(vl,tps[i],&nvl1); mergev(vl,nvl,nvl1,&nvl2);      clctv(vl,tps[i],&nvl1); mergev(vl,nvl,nvl1,&nvl2);
                 nvl = nvl2;      nvl = nvl2;
                 ezgcdpz(nvl,lg,LC(tps[i]),&t); lg = t;      ezgcdpz(nvl,lg,LC(tps[i]),&t); lg = t;
         }    }
   
         mulp(nvl,p0,lg,&lgp0);    mulp(nvl,p0,lg,&lgp0);
         k = dbound(v,lgp0) + 1;    k = dbound(v,lgp0) + 1;
         cbound(nvl,lgp0,(Q *)&cbd);    cbound(nvl,lgp0,(Q *)&cbd);
         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,vn0); W_CALLOC(nv,struct oVN,vnt);    W_CALLOC(nv,struct oVN,vn0); W_CALLOC(nv,struct oVN,vnt);
         W_CALLOC(nv,struct oVN,vn1);    W_CALLOC(nv,struct oVN,vn1);
         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 = vn0[i].v = tvl->v;      vn1[i].v = vn0[i].v = tvl->v;
   
         /* main loop */    /* main loop */
         for ( dg = deg(v,tps[0]) + 1; ; next(vn0) )    for ( dg = deg(v,tps[0]) + 1; ; next(vn0) )
                 do {      do {
                         for ( i = 0, j = 0; vn0[i].v; i++ )        for ( i = 0, j = 0; vn0[i].v; i++ )
                                 if ( vn0[i].n ) vnt[j++].v = (V)i;          if ( vn0[i].n ) vnt[j++].v = (V)i;
                         vnt[j].n = 0;        vnt[j].n = 0;
   
                         /* find b s.t. LC(p0)(b), LC(tps[i])(b) != 0 */        /* find b s.t. LC(p0)(b), LC(tps[i])(b) != 0 */
                         mulsgn(vn0,vnt,j,vn1);        mulsgn(vn0,vnt,j,vn1);
                         substvp(nvl,p0,vn1,&p00);        substvp(nvl,p0,vn1,&p00);
                         flag = (!zerovpchk(nvl,lp0,vn1) && sqfrchk(p00));        flag = (!zerovpchk(nvl,lp0,vn1) && sqfrchk(p00));
                         for ( i = 0; flag && (i < m); i++ )        for ( i = 0; flag && (i < m); i++ )
                                 flag &= (!zerovpchk(nvl,LC(tps[i]),vn1));          flag &= (!zerovpchk(nvl,LC(tps[i]),vn1));
                         if ( !flag )        if ( !flag )
                                 continue;          continue;
   
                         /* substitute y -> b */        /* substitute y -> b */
                         substvp(nvl,lg,vn1,&lg0);        substvp(nvl,lg,vn1,&lg0);
                         lp00 = LC(p00);        lp00 = LC(p00);
                         /* extended-gcd in 1-variable */        /* extended-gcd in 1-variable */
                         uexgcdnp(nvl,p00,tps,m,vn1,(Q)cbd,&g0,&h0,&a0,&b0,(Q *)&tq);        uexgcdnp(nvl,p00,tps,m,vn1,(Q)cbd,&g0,&h0,&a0,&b0,(Q *)&tq);
                         if ( NUM(g0) ) {        if ( NUM(g0) ) {
                                 *pr = (P)ONE;          *pr = (P)ONE;
                                 return;          return;
                         } else if ( dg > ( dg0 = deg(v,g0) ) ) {        } else if ( dg > ( dg0 = deg(v,g0) ) ) {
                                 dg = dg0;          dg = dg0;
                                 if ( dg == d0 ) {          if ( dg == d0 ) {
                                         if ( divcheck(nvl,tps,m,lp0,p0) ) {            if ( divcheck(nvl,tps,m,lp0,p0) ) {
                                                 *pr = p0;              *pr = p0;
                                                 return;               return;
                                         }            }
                                 } else if ( dg == deg(v,tps[0]) ) {          } else if ( dg == deg(v,tps[0]) ) {
                                         if ( divtpz(nvl,p0,tps[0],&t) &&            if ( divtpz(nvl,p0,tps[0],&t) &&
                                                 divcheck(nvl,tps+1,m-1,LC(tps[0]),tps[0]) ) {              divcheck(nvl,tps+1,m-1,LC(tps[0]),tps[0]) ) {
                                                 *pr = tps[0];              *pr = tps[0];
                                                 return;               return;
                                         } else            } else
                                                 break;              break;
                                 } else {          } else {
                                         henmv(nvl,vn1,lgp0,g0,h0,a0,b0,            henmv(nvl,vn1,lgp0,g0,h0,a0,b0,
                                                 lg,lp0,lg0,lp00,(Q)tq,k,&g,&h);              lg,lp0,lg0,lp00,(Q)tq,k,&g,&h);
                                         if ( divtpz(nvl,lgp0,g,&t) &&            if ( divtpz(nvl,lgp0,g,&t) &&
                                                 divcheck(nvl,tps,m,lg,g) ) {              divcheck(nvl,tps,m,lg,g) ) {
                                                 pcp(nvl,g,pr,&t);              pcp(nvl,g,pr,&t);
                                                 return;              return;
                                         }            }
                                 }          }
                         }        }
                 } while ( !nextbin(vnt,j) );      } while ( !nextbin(vnt,j) );
 }  }
   
 /*  /*
         f : sqfr    f : sqfr
 */  */
   
 void uezgcd1p(f,g,r)  void uezgcd1p(f,g,r)
 P f,g;  P f,g;
 P *r;  P *r;
 {  {
         UM wf,wf1,wf2,wg,wg1,wgcd,wcof;    UM wf,wf1,wf2,wg,wg1,wgcd,wcof;
         int d1,d2,dg,m,index,n;    int d1,d2,dg,m,index,n;
         ML list;    ML list;
         DCP dc;    DCP dc;
         P t;    P t;
         Q c;    Q c;
   
         if ( NUM(f) || NUM(g) ) {    if ( NUM(f) || NUM(g) ) {
                 *r = (P)ONE;      *r = (P)ONE;
                 return;      return;
         }    }
         ptozp(f,1,&c,&t); f = t; ptozp(g,1,&c,&t); g = t;    ptozp(f,1,&c,&t); f = t; ptozp(g,1,&c,&t); g = t;
         d1 = UDEG(f); d2 = UDEG(g);    d1 = UDEG(f); d2 = UDEG(g);
         n = MAX(d1,d2)+1;    n = MAX(d1,d2)+1;
         wf = W_UMALLOC(n); wf1 = W_UMALLOC(n);    wf = W_UMALLOC(n); wf1 = W_UMALLOC(n);
         wf2 = W_UMALLOC(n); wg = W_UMALLOC(n);    wf2 = W_UMALLOC(n); wg = W_UMALLOC(n);
         wg1 = W_UMALLOC(n); wgcd = W_UMALLOC(n);    wg1 = W_UMALLOC(n); wgcd = W_UMALLOC(n);
         wcof = W_UMALLOC(n);    wcof = W_UMALLOC(n);
   
         for ( index = INDEX, dg = MIN(d1,d2)+1; ; ) {    for ( index = INDEX, dg = MIN(d1,d2)+1; ; ) {
                 m = sprime[index++];      m = sprime[index++];
                 if ( !rem(NM((Q)UCOEF(f)),m) || !rem(NM((Q)UCOEF(g)),m))      if ( !rem(NM((Q)UCOEF(f)),m) || !rem(NM((Q)UCOEF(g)),m))
                         continue;        continue;
                 ptoum(m,f,wf); cpyum(wf,wf1);      ptoum(m,f,wf); cpyum(wf,wf1);
                 diffum(m,wf1,wf2); gcdum(m,wf1,wf2,wgcd);      diffum(m,wf1,wf2); gcdum(m,wf1,wf2,wgcd);
                 if ( DEG(wgcd) > 0 )      if ( DEG(wgcd) > 0 )
                         continue;        continue;
                 ptoum(m,g,wg); cpyum(wg,wg1); cpyum(wf,wf1);      ptoum(m,g,wg); cpyum(wg,wg1); cpyum(wf,wf1);
                 gcdum(m,wf1,wg1,wgcd);      gcdum(m,wf1,wg1,wgcd);
                 if ( dg <= DEG(wgcd) )      if ( dg <= DEG(wgcd) )
                         continue;        continue;
                 else      else
                         dg = DEG(wgcd);        dg = DEG(wgcd);
                 if ( dg == 0 ) {      if ( dg == 0 ) {
                         *r = (P)ONE;        *r = (P)ONE;
                         return;        return;
                 } else if ( dg == d1 ) {      } else if ( dg == d1 ) {
                         if ( divtpz(CO,g,f,&t) ) {        if ( divtpz(CO,g,f,&t) ) {
                                 *r = f;          *r = f;
                                 return;          return;
                         }        }
                 } else if ( dg == d2 ) {      } else if ( dg == d2 ) {
                         if ( divtpz(CO,f,g,&t) ) {        if ( divtpz(CO,f,g,&t) ) {
                                 *r = g;          *r = g;
                                 return;          return;
                         }        }
                 } else {      } else {
                         divum(m,wf,wgcd,wcof);        divum(m,wf,wgcd,wcof);
                         ezgcdhensel(f,m,wcof,wgcd,&list);        ezgcdhensel(f,m,wcof,wgcd,&list);
                         dtest(f,list,1,&dc);        dtest(f,list,1,&dc);
                         if ( NEXT(dc) ) {        if ( NEXT(dc) ) {
                                 if ( divtpz(CO,g,COEF(dc),&t) ) {          if ( divtpz(CO,g,COEF(dc),&t) ) {
                                         *r = COEF(dc);            *r = COEF(dc);
                                         return;            return;
                                 }          }
                         }        }
                 }      }
         }    }
 }  }
   
 void uexgcdnp(vl,p1,ps,n,vn,cbd,g,h,a,b,q)  void uexgcdnp(vl,p1,ps,n,vn,cbd,g,h,a,b,q)
 VL vl;  VL vl;
 P p1;  P p1;
Line 477  Q cbd;
Line 477  Q cbd;
 P *g,*h,*a,*b;  P *g,*h,*a,*b;
 Q *q;  Q *q;
 {  {
         UM wg,wh,wg1,wh1,wgcd,wt;    UM wg,wh,wg1,wh1,wgcd,wt;
         P t,s,gcd,cof,gk,hk,ak,bk;    P t,s,gcd,cof,gk,hk,ak,bk;
         Q c,tq,bb,tq1,qk;    Q c,tq,bb,tq1,qk;
         int mod,k,i,index,d;    int mod,k,i,index,d;
   
         ptozp(p1,1,&c,&gcd);    ptozp(p1,1,&c,&gcd);
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 substvp(vl,ps[i],vn,&t);      substvp(vl,ps[i],vn,&t);
                 uezgcd1p(gcd,t,&s);      uezgcd1p(gcd,t,&s);
                 if ( NUM(s) ) {      if ( NUM(s) ) {
                         *g = (P)ONE; *h = p1; *a = 0; *b = (P)ONE;        *g = (P)ONE; *h = p1; *a = 0; *b = (P)ONE;
                         return;        return;
                 } else      } else
                         gcd = s;        gcd = s;
         }    }
   
         if ( !dcomp(p1,gcd) ) {    if ( !dcomp(p1,gcd) ) {
                 *g = p1; *h = (P)ONE; *a = 0; *b = (P)ONE;      *g = p1; *h = (P)ONE; *a = 0; *b = (P)ONE;
                 return;      return;
         }    }
   
         divsp(CO,p1,gcd,&cof);    divsp(CO,p1,gcd,&cof);
         d = UDEG(p1)+1;    d = UDEG(p1)+1;
         wg = W_UMALLOC(d); wh = W_UMALLOC(d);    wg = W_UMALLOC(d); wh = W_UMALLOC(d);
         wg1 = W_UMALLOC(d); wh1 = W_UMALLOC(d);    wg1 = W_UMALLOC(d); wh1 = W_UMALLOC(d);
         wgcd = W_UMALLOC(d); wt = W_UMALLOC(d);    wgcd = W_UMALLOC(d); wt = W_UMALLOC(d);
         for ( index = INDEX; ; ) {    for ( index = INDEX; ; ) {
                 mod = sprime[index++];      mod = sprime[index++];
                 if ( !rem(NM((Q)LC(p1)),mod) )      if ( !rem(NM((Q)LC(p1)),mod) )
                         continue;        continue;
                 ptoum(mod,gcd,wg); ptoum(mod,cof,wh);      ptoum(mod,gcd,wg); ptoum(mod,cof,wh);
                 cpyum(wg,wg1); cpyum(wh,wh1);      cpyum(wg,wg1); cpyum(wh,wh1);
                 gcdum(mod,wg,wh,wt); cpyum(wt,wgcd);      gcdum(mod,wg,wh,wt); cpyum(wt,wgcd);
                 if ( DEG(wgcd) > 0 )      if ( DEG(wgcd) > 0 )
                         continue;        continue;
                 STOQ(mod,tq); addq(cbd,cbd,&bb);      STOQ(mod,tq); addq(cbd,cbd,&bb);
                 for ( k = 1; cmpq(tq,bb) < 0; k++ ) {      for ( k = 1; cmpq(tq,bb) < 0; k++ ) {
                         mulq(tq,tq,&tq1); tq = tq1;        mulq(tq,tq,&tq1); tq = tq1;
                 }      }
                 henzq(p1,gcd,wg1,cof,wh1,mod,k,&gk,&hk,&bk,&ak,&qk);      henzq(p1,gcd,wg1,cof,wh1,mod,k,&gk,&hk,&bk,&ak,&qk);
                 break;      break;
         }    }
         *g = gk; *h = hk; *a = ak; *b = bk; *q = tq;    *g = gk; *h = hk; *a = ak; *b = bk; *q = tq;
 }  }
   
 void ezgcdhensel(f,mod,cof,gcd,listp)  void ezgcdhensel(f,mod,cof,gcd,listp)
Line 528  int mod;
Line 528  int mod;
 UM gcd,cof;  UM gcd,cof;
 ML *listp;  ML *listp;
 {  {
         register int i,j;    register int i,j;
         int q,n,bound;    int q,n,bound;
         int *p;    int *p;
         int **pp;    int **pp;
         ML blist,clist,bqlist,cqlist,rlist;    ML blist,clist,bqlist,cqlist,rlist;
         UM *b;    UM *b;
         LUM fl,tl;    LUM fl,tl;
         LUM *l;    LUM *l;
         LUM LUMALLOC();    LUM LUMALLOC();
   
         blist = W_MLALLOC(2);    blist = W_MLALLOC(2);
         blist->n = 2; blist->mod = mod;    blist->n = 2; blist->mod = mod;
         blist->c[0] = (pointer)cof; blist->c[1] = (pointer)gcd;    blist->c[0] = (pointer)cof; blist->c[1] = (pointer)gcd;
         gcdgen(f,blist,&clist);    gcdgen(f,blist,&clist);
         henprep(f,blist,clist,&bqlist,&cqlist);    henprep(f,blist,clist,&bqlist,&cqlist);
         n = bqlist->n; q = bqlist->mod;    n = bqlist->n; q = bqlist->mod;
         bqlist->bound = cqlist->bound = bound = mignotte(q,f);    bqlist->bound = cqlist->bound = bound = mignotte(q,f);
   
         if ( bound == 1 ) {    if ( bound == 1 ) {
                 *listp = rlist = MLALLOC(n);      *listp = rlist = MLALLOC(n);
                 rlist->n = n;      rlist->n = n;
                 rlist->mod = q;      rlist->mod = q;
                 rlist->bound = bound;      rlist->bound = bound;
   
                 for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c;      for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c;
                         i < n; i++ ) {        i < n; i++ ) {
                         tl = LUMALLOC(DEG(b[i]),1);        tl = LUMALLOC(DEG(b[i]),1);
                         l[i] = tl;        l[i] = tl;
                         p = COEF(b[i]);        p = COEF(b[i]);
   
                         for ( j = 0, pp = COEF(tl);        for ( j = 0, pp = COEF(tl);
                                   j <= DEG(tl); j++ )            j <= DEG(tl); j++ )
                                         pp[j][0] = p[j];            pp[j][0] = p[j];
                 }      }
         } else {    } else {
                 W_LUMALLOC((int)UDEG(f),bound,fl);      W_LUMALLOC((int)UDEG(f),bound,fl);
                 ptolum(q,bound,f,fl);      ptolum(q,bound,f,fl);
                 henmain(fl,bqlist,cqlist,listp);      henmain(fl,bqlist,cqlist,listp);
         }    }
 }  }
   
 void ezgcdnpz(vl,ps,m,pr)  void ezgcdnpz(vl,ps,m,pr)
Line 574  VL vl;
Line 574  VL vl;
 P *ps,*pr;  P *ps,*pr;
 int m;  int m;
 {  {
         P t,s,gcd;    P t,s,gcd;
         P cont;    P cont;
         VL tvl,svl,nvl;    VL tvl,svl,nvl;
         int i;    int i;
         DCP dc;    DCP dc;
         N cn;    N cn;
         Q cq;    Q cq;
         P *pl,*ppl,*pcl;    P *pl,*ppl,*pcl;
         Q *cl;    Q *cl;
         V v;    V v;
   
         pl = (P *)ALLOCA((m+1)*sizeof(P));    pl = (P *)ALLOCA((m+1)*sizeof(P));
         cl = (Q *)ALLOCA((m+1)*sizeof(Q));    cl = (Q *)ALLOCA((m+1)*sizeof(Q));
         for ( i = 0; i < m; i++ )    for ( i = 0; i < m; i++ )
                 ptozp(ps[i],1,&cl[i],&pl[i]);      ptozp(ps[i],1,&cl[i],&pl[i]);
         for ( i = 1, cq = cl[0]; i < m; i++ ) {    for ( i = 1, cq = cl[0]; i < m; i++ ) {
                 gcdn(NM(cl[i]),NM(cq),&cn); NTOQ(cn,1,cq);      gcdn(NM(cl[i]),NM(cq),&cn); NTOQ(cn,1,cq);
         }    }
         for ( i = 0; i < m; i++ )    for ( i = 0; i < m; i++ )
                 if ( NUM(pl[i]) ) {      if ( NUM(pl[i]) ) {
                         *pr = (P)cq;        *pr = (P)cq;
                         return;        return;
                 }      }
   
         for ( i = 0, nvl = 0; i < m; i++ ) {    for ( i = 0, nvl = 0; i < m; i++ ) {
                 clctv(vl,ps[i],&tvl); mergev(vl,nvl,tvl,&svl); nvl = svl;      clctv(vl,ps[i],&tvl); mergev(vl,nvl,tvl,&svl); nvl = svl;
         }    }
   
         ppl = (P *)ALLOCA((m+1)*sizeof(P));    ppl = (P *)ALLOCA((m+1)*sizeof(P));
         pcl = (P *)ALLOCA((m+1)*sizeof(P));    pcl = (P *)ALLOCA((m+1)*sizeof(P));
         for ( i = 0, v = nvl->v; i < m; i++ )    for ( i = 0, v = nvl->v; i < m; i++ )
                 if ( v == VR(pl[i]) )      if ( v == VR(pl[i]) )
                         pcp(nvl,pl[i],&ppl[i],&pcl[i]);        pcp(nvl,pl[i],&ppl[i],&pcl[i]);
                 else {      else {
                         ppl[i] = (P)ONE; pcl[i] = pl[i];        ppl[i] = (P)ONE; pcl[i] = pl[i];
                 }      }
         ezgcdnpz(nvl,pcl,m,&cont);    ezgcdnpz(nvl,pcl,m,&cont);
   
         for ( i = 0; i < m; i++ )    for ( i = 0; i < m; i++ )
                 if ( NUM(ppl[i]) ) {      if ( NUM(ppl[i]) ) {
                         mulp(nvl,cont,(P)cq,pr);        mulp(nvl,cont,(P)cq,pr);
                         return;        return;
                 }      }
         sortplist(ppl,m);    sortplist(ppl,m);
         sqfrp(nvl,ppl[0],&dc);    sqfrp(nvl,ppl[0],&dc);
         for ( dc = NEXT(dc), gcd = (P)ONE; dc; dc = NEXT(dc) ) {    for ( dc = NEXT(dc), gcd = (P)ONE; dc; dc = NEXT(dc) ) {
                 if ( NUM(COEF(dc)) )      if ( NUM(COEF(dc)) )
                         continue;        continue;
                 ezgcdnpp(vl,dc,ppl+1,m-1,&t);      ezgcdnpp(vl,dc,ppl+1,m-1,&t);
                 if ( NUM(t) )      if ( NUM(t) )
                         continue;        continue;
                 mulp(vl,gcd,t,&s); gcd = s;      mulp(vl,gcd,t,&s); gcd = s;
         }    }
         mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,pr);    mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,pr);
 }  }
   
 void ezgcdnpp(vl,dc,pl,m,r)  void ezgcdnpp(vl,dc,pl,m,r)
Line 637  P *pl;
Line 637  P *pl;
 int m;  int m;
 P *r;  P *r;
 {  {
         int i,k;    int i,k;
         P g,t,s,gcd;    P g,t,s,gcd;
         P *pm;    P *pm;
   
         ezgcdnp(vl,COEF(dc),pl,m,&gcd);    ezgcdnp(vl,COEF(dc),pl,m,&gcd);
         if ( NUM(gcd) ) {    if ( NUM(gcd) ) {
                 *r = (P)ONE;      *r = (P)ONE;
                 return;      return;
         }    }
         pm = (P *) ALLOCA((m+1)*sizeof(P));    pm = (P *) ALLOCA((m+1)*sizeof(P));
         for ( i = 0; i < m; i++ ) {    for ( i = 0; i < m; i++ ) {
                 divsp(vl,pl[i],gcd,&pm[i]);      divsp(vl,pl[i],gcd,&pm[i]);
                 if ( NUM(pm[i]) ) {      if ( NUM(pm[i]) ) {
                         *r = gcd;        *r = gcd;
                         return;        return;
                 }      }
         }    }
         for ( g = gcd, k = QTOS(DEG(dc)) - 1; k > 0; k-- ) {    for ( g = gcd, k = QTOS(DEG(dc)) - 1; k > 0; k-- ) {
                 ezgcdnp(vl,g,pm,m,&t);      ezgcdnp(vl,g,pm,m,&t);
                 if ( NUM(t) )      if ( NUM(t) )
                         break;        break;
                 mulp(vl,gcd,t,&s); gcd = s;      mulp(vl,gcd,t,&s); gcd = s;
                 for ( i = 0; i < m; i++ ) {      for ( i = 0; i < m; i++ ) {
                         divsp(vl,pm[i],t,&s);        divsp(vl,pm[i],t,&s);
                         if ( NUM(s) ) {        if ( NUM(s) ) {
                                 *r = gcd;          *r = gcd;
                                 return;          return;
                         }        }
                         pm[i] = s;        pm[i] = s;
                 }      }
         }    }
         *r = gcd;    *r = gcd;
 }  }
   
 void pcp(vl,p,pp,c)  void pcp(vl,p,pp,c)
Line 676  VL vl;
Line 676  VL vl;
 P p;  P p;
 P *pp,*c;  P *pp,*c;
 {  {
         P f,g,n;    P f,g,n;
         DCP dc;    DCP dc;
         P *pl;    P *pl;
         int i,m;    int i,m;
         extern int nez;    extern int nez;
   
         if ( NUM(p) ) {    if ( NUM(p) ) {
                 *c = p;      *c = p;
                 *pp = (P)ONE;      *pp = (P)ONE;
         } else {    } else {
                 for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );      for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
                 if ( m == 1 ) {      if ( m == 1 ) {
                         *c = COEF(DC(p));        *c = COEF(DC(p));
                         divsp(vl,p,*c,pp);        divsp(vl,p,*c,pp);
                         return;        return;
                 }      }
                 ptozp(p,1,(Q *)&n,&f);      ptozp(p,1,(Q *)&n,&f);
                 pl = (P *)ALLOCA((m+1)*sizeof(P));      pl = (P *)ALLOCA((m+1)*sizeof(P));
                 for ( i = 0, dc = DC(f); i < m; dc = NEXT(dc), i++ )      for ( i = 0, dc = DC(f); i < m; dc = NEXT(dc), i++ )
                         pl[i] = COEF(dc);        pl[i] = COEF(dc);
                 if ( nez )      if ( nez )
                         nezgcdnpz(vl,pl,m,&g);        nezgcdnpz(vl,pl,m,&g);
                 else      else
                         ezgcdnpz(vl,pl,m,&g);        ezgcdnpz(vl,pl,m,&g);
                 mulp(vl,g,n,c); divsp(vl,f,g,pp);      mulp(vl,g,n,c); divsp(vl,f,g,pp);
         }    }
 }  }
   
 #if 0  #if 0
Line 710  VL vl;
Line 710  VL vl;
 P p;  P p;
 P *pp,*c;  P *pp,*c;
 {  {
         P f,g,n,t;    P f,g,n,t;
         DCP dc;    DCP dc;
   
         if ( NUM(p) ) {    if ( NUM(p) ) {
                 *c = p;      *c = p;
                 *pp = (P)ONE;      *pp = (P)ONE;
         } else {    } else {
                 ptozp(p,1,&n,&f);      ptozp(p,1,&n,&f);
                 for ( g = LC(f), dc = NEXT(DC(f)); dc; dc = NEXT(dc) ) {      for ( g = LC(f), dc = NEXT(DC(f)); dc; dc = NEXT(dc) ) {
                         ezgcdpz(vl,g,COEF(dc),&t); g = t;        ezgcdpz(vl,g,COEF(dc),&t); g = t;
                 }      }
                 mulp(vl,g,n,c);      mulp(vl,g,n,c);
                 divsp(vl,f,g,pp);      divsp(vl,f,g,pp);
         }    }
 }  }
 #endif  #endif
   
Line 732  VL vl;
Line 732  VL vl;
 int m;  int m;
 P *ps,l,p;  P *ps,l,p;
 {  {
         int i;    int i;
         P r,t;    P r,t;
   
         for ( i = 0; i < m; i++ ) {    for ( i = 0; i < m; i++ ) {
                 mulp(vl,ps[i],l,&t);      mulp(vl,ps[i],l,&t);
                 if ( !divtpz(vl,t,p,&r) )      if ( !divtpz(vl,t,p,&r) )
                         return ( 0 );        return ( 0 );
         }    }
         return ( 1 );    return ( 1 );
 }  }
   
 void sortplist(plist,n)  void sortplist(plist,n)
 P *plist;  P *plist;
 int n;  int n;
 {  {
         int i,j;    int i,j;
         P t;    P t;
   
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 for ( j = i + 1; j < n; j++ )      for ( j = i + 1; j < n; j++ )
                         if ( cmpq(DEG(DC(plist[i])),DEG(DC(plist[j]))) > 0 ) {        if ( cmpq(DEG(DC(plist[i])),DEG(DC(plist[j]))) > 0 ) {
                                 t = plist[i]; plist[i] = plist[j]; plist[j] = t;          t = plist[i]; plist[i] = plist[j]; plist[j] = t;
                         }        }
 }  }
   
 int dcomp(p1,p2)  int dcomp(p1,p2)
 P p1,p2;  P p1,p2;
 {  {
         return ( cmpq(DEG(DC(p1)),DEG(DC(p2))) );    return ( cmpq(DEG(DC(p1)),DEG(DC(p2))) );
 }  }
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>