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

Diff for /OpenXM_contrib2/asir2000/engine/E.c between version 1.2 and 1.9

version 1.2, 2000/08/21 08:31:24 version 1.9, 2018/03/29 01:32:51
Line 23 
Line 23 
  * shall be made on your publication or presentation in any form of the   * shall be made on your publication or presentation in any form of the
  * results obtained by use of the SOFTWARE.   * results obtained by use of the SOFTWARE.
  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by   * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
  * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification   * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
  * for such modification or the source code of the modified part of the   * for such modification or the source code of the modified part of the
  * SOFTWARE.   * SOFTWARE.
  *   *
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/E.c,v 1.1.1.1 1999/12/03 07:39:07 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/E.c,v 1.8 2001/10/09 01:36:09 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
   
 void henmv(vl,vn,f,g0,h0,a0,b0,lg,lh,lg0,lh0,q,k,gp,hp)  void henmv(VL vl,VN vn,P f,P g0,P h0,P a0,P b0,P lg,P lh,P lg0,P lh0,Q q,int k,P *gp,P *hp)
 VL vl;  
 VN vn;  
 P f,g0,h0,a0,b0,lg,lh,lg0,lh0;  
 Q q;  
 int k;  
 P *gp,*hp;  
 {  {
         P g1,h1,a1,b1;    P g1,h1,a1,b1;
         N qn;    N qn;
         Q q2;    Q q2;
   
         divin((NM(q)),2,&qn); NTOQ(qn,1,q2);    divin((NM(q)),2,&qn); NTOQ(qn,1,q2);
         adjc(vl,g0,a0,lg0,q,&g1,&a1); adjc(vl,h0,b0,lh0,q,&h1,&b1);    adjc(vl,g0,a0,lg0,q,&g1,&a1); adjc(vl,h0,b0,lh0,q,&h1,&b1);
         henmvmain(vl,vn,f,g1,h1,b1,a1,lg,lh,q,q2,k,gp,hp);    henmvmain(vl,vn,f,g1,h1,b1,a1,lg,lh,q,q2,k,gp,hp);
 }  }
   
 void henmvmain(vl,vn,f,fi0,fi1,gi0,gi1,l0,l1,mod,mod2,k,fr0,fr1)  void henmvmain(VL vl,VN vn,P f,P fi0,P fi1,P gi0,P gi1,P l0,P l1,Q mod,Q mod2,int k,P *fr0,P *fr1)
 VL vl;  
 VN vn;  
 P f,fi0,fi1,gi0,gi1,l0,l1;  
 Q mod,mod2;  
 int k;  
 P *fr0,*fr1;  
 {  {
         V v;    V v;
         int n,i,j;    int n,i,j;
         int *md;    int *md;
         P x,m,m1,c,q,r,a,s,u,ff,f0,f1;    P x,m,m1,c,q,r,a,s,u,ff,f0,f1;
         P w0,w1,cf,cfi,t,q1,dvr;    P w0,w1,cf,cfi,t,q1,dvr;
         P *c0,*c1;    P *c0,*c1;
         P *f0h,*f1h;    P *f0h,*f1h;
   
         v = VR(f); n = deg(v,f); MKV(v,x);    v = VR(f); n = deg(v,f); MKV(v,x);
         c0 = (P *)ALLOCA((n+1)*sizeof(P));    c0 = (P *)ALLOCA((n+1)*sizeof(P));
         c1 = (P *)ALLOCA((n+1)*sizeof(P));    c1 = (P *)ALLOCA((n+1)*sizeof(P));
         invl((Q)LC(fi1),mod,(Q *)&c); mulp(vl,fi1,c,&t); cmp(mod,t,&dvr);    invl((Q)LC(fi1),mod,(Q *)&c); mulp(vl,fi1,c,&t); cmp(mod,t,&dvr);
         cm2p(mod,mod2,gi0,&c0[0]); cm2p(mod,mod2,gi1,&c1[0]);    cm2p(mod,mod2,gi0,&c0[0]); cm2p(mod,mod2,gi1,&c1[0]);
         for ( i = 1; i <= n; i++ ) {    for ( i = 1; i <= n; i++ ) {
                 mulp(vl,x,c1[i-1],&m); divsrp(vl,m,dvr,&q,&r); mulp(vl,q,c,&q1);      mulp(vl,x,c1[i-1],&m); divsrp(vl,m,dvr,&q,&r); mulp(vl,q,c,&q1);
                 cm2p(mod,mod2,r,&c1[i]);      cm2p(mod,mod2,r,&c1[i]);
                 mulp(vl,x,c0[i-1],&m); mulp(vl,q1,fi0,&m1); addp(vl,m,m1,&a);      mulp(vl,x,c0[i-1],&m); mulp(vl,q1,fi0,&m1); addp(vl,m,m1,&a);
                 cm2p(mod,mod2,a,&c0[i]);      cm2p(mod,mod2,a,&c0[i]);
         }    }
         affine(vl,f,vn,&t); cm2p(mod,mod2,t,&ff);    affine(vl,f,vn,&t); cm2p(mod,mod2,t,&ff);
         for ( i = 0; vn[i].v; i++ );    for ( i = 0; vn[i].v; i++ );
         md = ( int *) ALLOCA((i+1)*sizeof(int));    md = ( int *) ALLOCA((i+1)*sizeof(int));
         for ( i = 0; vn[i].v; i++ )    for ( i = 0; vn[i].v; i++ )
                 md[i] = getdeg(vn[i].v,ff);      md[i] = getdeg(vn[i].v,ff);
         cm2p(mod,mod2,fi0,&f0); affine(vl,l0,vn,&t);    cm2p(mod,mod2,fi0,&f0); affine(vl,l0,vn,&t);
         if ( NUM(f0) )    if ( NUM(f0) )
                 cm2p(mod,mod2,t,&f0);      cm2p(mod,mod2,t,&f0);
         else    else
                 cm2p(mod,mod2,t,&COEF(DC(f0)));      cm2p(mod,mod2,t,&COEF(DC(f0)));
         cm2p(mod,mod2,fi1,&f1); affine(vl,l1,vn,&t);    cm2p(mod,mod2,fi1,&f1); affine(vl,l1,vn,&t);
         if ( NUM(f1) )    if ( NUM(f1) )
                 cm2p(mod,mod2,t,&f1);      cm2p(mod,mod2,t,&f1);
         else    else
                 cm2p(mod,mod2,t,&COEF(DC(f1)));      cm2p(mod,mod2,t,&COEF(DC(f1)));
         W_CALLOC(k,P,f0h); W_CALLOC(k,P,f1h);    W_CALLOC(k,P,f0h); W_CALLOC(k,P,f1h);
         for ( i = 0; i <= k; i++ ) {    for ( i = 0; i <= k; i++ ) {
                 exthpc(vl,v,f0,i,&f0h[i]); exthpc(vl,v,f1,i,&f1h[i]);      exthpc(vl,v,f0,i,&f0h[i]); exthpc(vl,v,f1,i,&f1h[i]);
         }    }
         for ( j = 1; j <= k; j++ ) {    for ( j = 1; j <= k; j++ ) {
                 for ( i = 0; vn[i].v; i++ )      for ( i = 0; vn[i].v; i++ )
                         if ( getdeg(vn[i].v,f0)+getdeg(vn[i].v,f1) > md[i] )        if ( getdeg(vn[i].v,f0)+getdeg(vn[i].v,f1) > md[i] )
                                 goto END;          goto END;
                 for ( i = 0, s = 0; i <= j; i++ ) {      for ( i = 0, s = 0; i <= j; i++ ) {
                         mulp(vl,f0h[i],f1h[j-i],&t); addp(vl,s,t,&u); s = u;        mulp(vl,f0h[i],f1h[j-i],&t); addp(vl,s,t,&u); s = u;
                 }      }
                 cm2p(mod,mod2,s,&t);      cm2p(mod,mod2,s,&t);
                 exthpc(vl,v,ff,j,&u); subm2p(vl,mod,mod2,u,t,&cf);      exthpc(vl,v,ff,j,&u); subm2p(vl,mod,mod2,u,t,&cf);
                 for ( i = 0, w0 = 0, w1 = 0; i <= n; i++ ) {      for ( i = 0, w0 = 0, w1 = 0; i <= n; i++ ) {
                         if ( !cf )        if ( !cf )
                                 cfi = 0;          cfi = 0;
                         else if ( VR(cf) == v )        else if ( VR(cf) == v )
                                 coefp(cf,i,&cfi);          coefp(cf,i,&cfi);
                         else if ( i == 0 )        else if ( i == 0 )
                                 cfi = cf;          cfi = cf;
                         else        else
                                 cfi = 0;          cfi = 0;
                         if ( cfi ) {        if ( cfi ) {
                                 mulp(vl,cfi,c0[i],&m); addp(vl,w0,m,&a); w0 = a;          mulp(vl,cfi,c0[i],&m); addp(vl,w0,m,&a); w0 = a;
                                 mulp(vl,cfi,c1[i],&m); addp(vl,w1,m,&a); w1 = a;          mulp(vl,cfi,c1[i],&m); addp(vl,w1,m,&a); w1 = a;
                         }        }
                 }      }
                 cm2p(mod,mod2,w0,&t); addm2p(vl,mod,mod2,f0,t,&a);      cm2p(mod,mod2,w0,&t); addm2p(vl,mod,mod2,f0,t,&a);
                 addm2p(vl,mod,mod2,f0h[j],t,&s); f0h[j] = s; f0 = a;      addm2p(vl,mod,mod2,f0h[j],t,&s); f0h[j] = s; f0 = a;
                 cm2p(mod,mod2,w1,&u); addm2p(vl,mod,mod2,f1,u,&a);      cm2p(mod,mod2,w1,&u); addm2p(vl,mod,mod2,f1,u,&a);
                 addm2p(vl,mod,mod2,f1h[j],u,&s); f1h[j] = s; f1 = a;      addm2p(vl,mod,mod2,f1h[j],u,&s); f1h[j] = s; f1 = a;
                 if ( !t ) {      if ( !t ) {
                         restore(vl,f0,vn,&s); cm2p(mod,mod2,s,&t);        restore(vl,f0,vn,&s); cm2p(mod,mod2,s,&t);
                         if ( divtpz(vl,f,t,&s) ) {        if ( divtpz(vl,f,t,&s) ) {
                                 *fr0 = t; *fr1 = s;          *fr0 = t; *fr1 = s;
                                 return;          return;
                         }        }
                 }      }
                 if ( !u ) {      if ( !u ) {
                         restore(vl,f1,vn,&s); cm2p(mod,mod2,s,&t);        restore(vl,f1,vn,&s); cm2p(mod,mod2,s,&t);
                         if ( divtpz(vl,f,t,&s) ) {        if ( divtpz(vl,f,t,&s) ) {
                                 *fr0 = s; *fr1 = t;          *fr0 = s; *fr1 = t;
                                 return;          return;
                         }        }
                 }      }
         }    }
 END:  END:
         restore(vl,f0,vn,&t); cm2p(mod,mod2,t,fr0);    restore(vl,f0,vn,&t); cm2p(mod,mod2,t,fr0);
         restore(vl,f1,vn,&t); cm2p(mod,mod2,t,fr1);    restore(vl,f1,vn,&t); cm2p(mod,mod2,t,fr1);
 }  }
   
 /*  /*
         input : f, fi0, fi1, p, k; f = fi0 * fi1 mod p; ( p, k : integer )    input : f, fi0, fi1, p, k; f = fi0 * fi1 mod p; ( p, k : integer )
         output : gr1 * fr0 + gr0 * fr1 = 1 mod qr; ( qr = p^(2^(k-1)) )    output : gr1 * fr0 + gr0 * fr1 = 1 mod qr; ( qr = p^(2^(k-1)) )
 */  */
   
 void henzq(f,i0,fi0,i1,fi1,p,k,fr0p,fr1p,gr0p,gr1p,qrp)  void henzq(P f,P i0,UM fi0,P i1,UM fi1,int p,int k,P *fr0p,P *fr1p,P *gr0p,P *gr1p,Q *qrp)
 P f;  
 UM fi0,fi1;  
 int p,k;  
 P i0,i1;  
 P *fr0p,*fr1p,*gr0p,*gr1p;  
 Q *qrp;  
 {  {
         N qn;    N qn;
         Q q,qq,q2;    Q q,qq,q2;
         int n,i;    int n,i;
         UM wg0,wg1,wf0,wf1;    UM wg0,wg1,wf0,wf1;
         P f0,f1,g0,g1,m,m1,s,c,d,q1,r1,rm,rm1,a,a1,a2;    P f0,f1,g0,g1,m,m1,s,c,d,q1,r1,rm,rm1,a,a1,a2;
   
         n = UDEG(f);    n = UDEG(f);
         wf0 = W_UMALLOC(n); wf1 = W_UMALLOC(n);    wf0 = W_UMALLOC(n); wf1 = W_UMALLOC(n);
         wg0 = W_UMALLOC(n); wg1 = W_UMALLOC(n);    wg0 = W_UMALLOC(n); wg1 = W_UMALLOC(n);
         cpyum(fi0,wf0); cpyum(fi1,wf1);    cpyum(fi0,wf0); cpyum(fi1,wf1);
         eucum(p,wf0,wf1,wg1,wg0);    eucum(p,wf0,wf1,wg1,wg0);
         umtop(VR(f),fi0,&f0); umtop(VR(f),fi1,&f1);    umtop(VR(f),fi0,&f0); umtop(VR(f),fi1,&f1);
         umtop(VR(f),wg0,&g0); umtop(VR(f),wg1,&g1);    umtop(VR(f),wg0,&g0); umtop(VR(f),wg1,&g1);
   
         STOQ(p,q); divin(NM(q),2,&qn); NTOQ(qn,1,q2);    STOQ(p,q); divin(NM(q),2,&qn); NTOQ(qn,1,q2);
         for ( i = 1; i < k; i++ ) {    for ( i = 1; i < k; i++ ) {
 #if 0  #if 0
                 mulp(CO,i1,g0,&m); mulp(CO,i0,g1,&m1); addp(CO,m,m1,&a);      mulp(CO,i1,g0,&m); mulp(CO,i0,g1,&m1); addp(CO,m,m1,&a);
                 if ( NUM(a) ) {      if ( NUM(a) ) {
                         for ( STOQ(p,q), j = 1; j < k; j++ ) {        for ( STOQ(p,q), j = 1; j < k; j++ ) {
                                 mulq(q,q,&qq); q = qq;          mulq(q,q,&qq); q = qq;
                         }        }
                         f0 = i0; f1 = i1;        f0 = i0; f1 = i1;
                         invl(a,q,&qq);        invl(a,q,&qq);
                         mulp(CO,g0,qq,&s); g0 = s; mulp(CO,g1,qq,&s); g1 = s;        mulp(CO,g0,qq,&s); g0 = s; mulp(CO,g1,qq,&s); g1 = s;
                         break;        break;
                 }      }
 #endif  #endif
         /*      c = ((f - f0*f1)/q) mod q;    /*  c = ((f - f0*f1)/q) mod q;
                 q1 = (c*g1) / f1;      q1 = (c*g1) / f1;
                 r1 = (c*g1) mod f1;      r1 = (c*g1) mod f1;
                 f1 += (r1 mod q)*q;      f1 += (r1 mod q)*q;
                 f0 += ((c*g0 + q1*f0) mod q)*q;      f0 += ((c*g0 + q1*f0) mod q)*q;
   
                 d = ((1 - (f1*g0 + f0*g1))/q) mod q;      d = ((1 - (f1*g0 + f0*g1))/q) mod q;
                 q1 = (d*g0) / f1;      q1 = (d*g0) / f1;
                 r1 = (d*g0) mod f1;      r1 = (d*g0) mod f1;
                 g1 += (r1 mod q)*q;      g1 += (r1 mod q)*q;
                 g0 += ((c*g0 + q1*f0) mod q)*q;      g0 += ((c*g0 + q1*f0) mod q)*q;
                 q = q^2;      q = q^2;
         */    */
   
         /* c = ((f - f0*f1)/q) mod q */    /* c = ((f - f0*f1)/q) mod q */
                 mulp(CO,f0,f1,&m); subp(CO,f,m,&s);      mulp(CO,f0,f1,&m); subp(CO,f,m,&s);
                 divcp(s,q,&m); cm2p(q,q2,m,&c);      divcp(s,q,&m); cm2p(q,q2,m,&c);
   
         /* q1 = (c*g1) / f1; r1 = (c*g1) mod f1; */    /* q1 = (c*g1) / f1; r1 = (c*g1) mod f1; */
                 mulp(CO,c,g1,&m); cm2p(q,q2,m,&s);      mulp(CO,c,g1,&m); cm2p(q,q2,m,&s);
                 udivpwm(q,s,f1,&q1,&r1);      udivpwm(q,s,f1,&q1,&r1);
   
         /* f1 = f1 + (r1 mod q)*q; */    /* f1 = f1 + (r1 mod q)*q; */
                 cm2p(q,q2,r1,&rm);      cm2p(q,q2,r1,&rm);
                 mulpq(rm,(P)q,&rm1); addp(CO,f1,rm1,&a);      mulpq(rm,(P)q,&rm1); addp(CO,f1,rm1,&a);
                 f1 = a;      f1 = a;
   
         /* a1 = (c*g0 + q1*f0) mod q; */    /* a1 = (c*g0 + q1*f0) mod q; */
                 mulp(CO,c,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);      mulp(CO,c,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);
                 cm2p(q,q2,a,&a1);      cm2p(q,q2,a,&a1);
   
         /* f0 = f0 + a1*q; */    /* f0 = f0 + a1*q; */
                 mulpq(a1,(P)q,&a2);      mulpq(a1,(P)q,&a2);
                 addp(CO,f0,a2,&a);      addp(CO,f0,a2,&a);
                 f0 = a;      f0 = a;
   
         /* d = ((1 - (f1*g0 + f0*g1))/q) mod q; */    /* d = ((1 - (f1*g0 + f0*g1))/q) mod q; */
                 mulp(CO,f1,g0,&m); mulp(CO,f0,g1,&m1); addp(CO,m,m1,&a);      mulp(CO,f1,g0,&m); mulp(CO,f0,g1,&m1); addp(CO,m,m1,&a);
                 subp(CO,(P)ONE,a,&s);      subp(CO,(P)ONE,a,&s);
                 divcp(s,q,&m); cm2p(q,q2,m,&d);      divcp(s,q,&m); cm2p(q,q2,m,&d);
   
         /* q1 = (d*g0) / f1; r1 = (d*g0) mod f1; */    /* q1 = (d*g0) / f1; r1 = (d*g0) mod f1; */
                 mulp(CO,d,g1,&m); udivpwm(q,m,f1,&q1,&r1);      mulp(CO,d,g1,&m); udivpwm(q,m,f1,&q1,&r1);
   
         /* g1 = g1 + (r1 mod q )*q; */    /* g1 = g1 + (r1 mod q )*q; */
                 cm2p(q,q2,r1,&rm);      cm2p(q,q2,r1,&rm);
                 mulpq(rm,(P)q,&rm1); addp(CO,g1,rm1,&a);      mulpq(rm,(P)q,&rm1); addp(CO,g1,rm1,&a);
                 g1 = a;      g1 = a;
   
         /* g0 = g0 + (c*g0 + q1*f0 mod q)*q; */    /* g0 = g0 + (c*g0 + q1*f0 mod q)*q; */
                 mulp(CO,d,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);      mulp(CO,d,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);
                 cm2p(q,q2,a,&a1); mulpq(a1,(P)q,&a2);      cm2p(q,q2,a,&a1); mulpq(a1,(P)q,&a2);
                 addp(CO,g0,a2,&a);      addp(CO,g0,a2,&a);
                 g0 = a;      g0 = a;
   
         /* q = q^2; */    /* q = q^2; */
                 mulq(q,q,&qq);      mulq(q,q,&qq);
                 q = qq;      q = qq;
                 divin(NM(q),2,&qn); NTOQ(qn,1,q2);      divin(NM(q),2,&qn); NTOQ(qn,1,q2);
         }    }
         *fr0p = f0; *fr1p = f1; *gr0p = g0; *gr1p = g1; *qrp = q;    *fr0p = f0; *fr1p = f1; *gr0p = g0; *gr1p = g1; *qrp = q;
 }  }
   
 void henzq1(g,h,bound,gcp,hcp,qp)  void henzq1(P g,P h,Q bound,P *gcp,P *hcp,Q *qp)
 P g,h;  
 Q bound;  
 P *gcp,*hcp;  
 Q *qp;  
 {  {
         V v;    V v;
         Q f,q,q1;    Q f,q,q1;
         Q u,t,a,b,s;    Q u,t,a,b,s;
         P c,c1;    P c,c1;
         P tg,th,tr;    P tg,th,tr;
         UM wg,wh,ws,wt,wm;    UM wg,wh,ws,wt,wm;
         int n,m,modres,mod,index,i;    int n,m,modres,mod,index,i;
         P gc0,hc0;    P gc0,hc0;
         P z,zz,zzz;    P z,zz,zzz;
   
   
         v = VR(g); n=deg(v,g); m=deg(v,h);    v = VR(g); n=deg(v,g); m=deg(v,h);
         norm(g,&a); norm(h,&b);    norm(g,&a); norm(h,&b);
         STOQ(m,u); pwrq(a,u,&t);    STOQ(m,u); pwrq(a,u,&t);
         STOQ(n,u); pwrq(b,u,&s);    STOQ(n,u); pwrq(b,u,&s);
         mulq(t,s,&u);    mulq(t,s,&u);
   
         factorial(n+m,&t); mulq(u,t,&s);    factorial(n+m,&t); mulq(u,t,&s);
         addq(s,s,&f);    addq(s,s,&f);
   
         wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);    wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);
         wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);    wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);
         wm = W_UMALLOC(m+n);    wm = W_UMALLOC(m+n);
   
         for ( q = ONE, t = 0, c = 0, index = 0; ; ) {    for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
                 mod = lprime[index++];      mod = get_lprime(index++);
                 if ( !mod )      if ( !rem(NM((Q)LC(g)),mod) || !rem(NM((Q)LC(h)),mod) )
                         error("henzq1 : lprime[] exhausted.");        continue;
                 if ( !rem(NM((Q)LC(g)),mod) || !rem(NM((Q)LC(h)),mod) )      ptomp(mod,g,&tg); ptomp(mod,h,&th);
                         continue;      srchump(mod,tg,th,&tr);
                 ptomp(mod,g,&tg); ptomp(mod,h,&th);      if ( !tr )
                 srchump(mod,tg,th,&tr);        continue;
                 if ( !tr )      else
                         continue;        modres = CONT((MQ)tr);
                 else  
                         modres = CONT((MQ)tr);  
   
                 mptoum(tg,wg); mptoum(th,wh);      mptoum(tg,wg); mptoum(th,wh);
                 eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */      eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
                 DEG(wm) = 0; COEF(wm)[0] = modres;      DEG(wm) = 0; COEF(wm)[0] = modres;
                 mulum(mod,ws,wm,wt);      mulum(mod,ws,wm,wt);
                 for ( i = DEG(wt); i >= 0; i-- )      for ( i = DEG(wt); i >= 0; i-- )
                         if ( ( COEF(wt)[i] * 2 ) > mod )        if ( ( COEF(wt)[i] * 2 ) > mod )
                                 COEF(wt)[i] -= mod;          COEF(wt)[i] -= mod;
                 chnrem(mod,v,c,q,wt,&c1,&q1);      chnrem(mod,v,c,q,wt,&c1,&q1);
                 if ( !ucmpp(c,c1) ) {      if ( !ucmpp(c,c1) ) {
                         mulp(CO,c,g,&z); divsrp(CO,z,h,&zz,&zzz);        mulp(CO,c,g,&z); divsrp(CO,z,h,&zz,&zzz);
                         if ( NUM(zzz) ) {        if ( NUM(zzz) ) {
                                 q = q1; c = c1;          q = q1; c = c1;
                                 break;          break;
                         }        }
                 }      }
                 q = q1; c = c1;      q = q1; c = c1;
   
                 if ( cmpq(f,q) < 0 )      if ( cmpq(f,q) < 0 )
                         break;        break;
         }    }
         ptozp(c,1,&s,&gc0);    ptozp(c,1,&s,&gc0);
         /* gc0*g = zz*h + zzz -> c*g + (-zz)*h = zzz */    /* gc0*g = zz*h + zzz -> c*g + (-zz)*h = zzz */
         mulp(CO,gc0,g,&z);    mulp(CO,gc0,g,&z);
         divsrp(CO,z,h,&zz,&zzz);    divsrp(CO,z,h,&zz,&zzz);
         ptozp(zz,1,&s,(P *)&t);    ptozp(zz,1,&s,(P *)&t);
         if ( INT((Q)s) )    if ( INT((Q)s) )
                 chsgnp(zz,&hc0);      chsgnp(zz,&hc0);
         else {    else {
                 NTOQ(DN((Q)s),1,q); mulq(q,(Q)zzz,&q1); zzz = (P)q1;      NTOQ(DN((Q)s),1,q); mulq(q,(Q)zzz,&q1); zzz = (P)q1;
                 mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,zz,(P)q,&c); chsgnp(c,&hc0);      mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,zz,(P)q,&c); chsgnp(c,&hc0);
         }    }
         if ( !INT((Q)zzz) ) {    if ( !INT((Q)zzz) ) {
                 NTOQ(DN((Q)zzz),1,q); NTOQ(NM((Q)zzz),1,q1); zzz = (P)q1;      NTOQ(DN((Q)zzz),1,q); NTOQ(NM((Q)zzz),1,q1); zzz = (P)q1;
                 mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,hc0,(P)q,&c); hc0 = c;      mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,hc0,(P)q,&c); hc0 = c;
         }    }
         for ( index = 0; ; ) {    for ( index = 0; ; ) {
                 mod = lprime[index++];      mod = get_lprime(index++);
                 if ( !mod )      if ( !rem(NM((Q)zzz),mod) ||
                         error("henzq1 : lprime[] exhausted.");        !rem(NM((Q)LC(g)),mod) ||
                 if ( !rem(NM((Q)zzz),mod) ||        !rem(NM((Q)LC(h)),mod) )
                         !rem(NM((Q)LC(g)),mod) ||        continue;
                         !rem(NM((Q)LC(h)),mod) )      for ( STOQ(mod,q); cmpq(q,bound) < 0; ) {
                         continue;        mulq(q,q,&q1); q = q1;
                 for ( STOQ(mod,q); cmpq(q,bound) < 0; ) {      }
                         mulq(q,q,&q1); q = q1;      *qp = q;
                 }      invl((Q)zzz,q,&q1);
                 *qp = q;      mulp(CO,gc0,(P)q1,hcp); mulp(CO,hc0,(P)q1,gcp);
                 invl((Q)zzz,q,&q1);      return;
                 mulp(CO,gc0,(P)q1,hcp); mulp(CO,hc0,(P)q1,gcp);    }
                 return;  
         }  
 }  }
   
 void addm2p(vl,mod,mod2,n1,n2,nr)  void addm2p(VL vl,Q mod,Q mod2,P n1,P n2,P *nr)
 VL vl;  
 Q mod,mod2;  
 P n1,n2,*nr;  
 {  {
         P t;    P t;
   
         addp(vl,n1,n2,&t);    addp(vl,n1,n2,&t);
         if ( !t )    if ( !t )
                 *nr = 0;      *nr = 0;
         else    else
                 cm2p(mod,mod2,t,nr);      cm2p(mod,mod2,t,nr);
 }  }
   
 void subm2p(vl,mod,mod2,n1,n2,nr)  void subm2p(VL vl,Q mod,Q mod2,P n1,P n2,P *nr)
 VL vl;  
 Q mod,mod2;  
 P n1,n2,*nr;  
 {  {
         P t;    P t;
   
         subp(vl,n1,n2,&t);    subp(vl,n1,n2,&t);
         if ( !t )    if ( !t )
                 *nr = 0;      *nr = 0;
         else    else
                 cm2p(mod,mod2,t,nr);      cm2p(mod,mod2,t,nr);
 }  }
   
 void mulm2p(vl,mod,mod2,n1,n2,nr)  void mulm2p(VL vl,Q mod,Q mod2,P n1,P n2,P *nr)
 VL vl;  
 Q mod,mod2;  
 P n1,n2,*nr;  
 {  {
         P t;    P t;
   
         mulp(vl,n1,n2,&t);    mulp(vl,n1,n2,&t);
         if ( !t )    if ( !t )
                 *nr = 0;      *nr = 0;
         else    else
                 cm2p(mod,mod2,t,nr);      cm2p(mod,mod2,t,nr);
 }  }
   
 void cmp(mod,p,pr)  void cmp(Q mod,P p,P *pr)
 Q mod;  
 P p,*pr;  
 {  {
         P t;    P t;
         DCP dc,dcr,dcr0;    DCP dc,dcr,dcr0;
   
         if ( !p )    if ( !p )
                 *pr = 0;      *pr = 0;
         else if ( NUM(p) )    else if ( NUM(p) )
                 remq((Q)p,mod,(Q *)pr);      remq((Q)p,mod,(Q *)pr);
         else {    else {
                 for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {      for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
                         cmp(mod,COEF(dc),&t);        cmp(mod,COEF(dc),&t);
                         if ( t ) {        if ( t ) {
                                 NEXTDC(dcr0,dcr);          NEXTDC(dcr0,dcr);
                                 DEG(dcr) = DEG(dc);          DEG(dcr) = DEG(dc);
                                 COEF(dcr) = t;          COEF(dcr) = t;
                         }        }
                 }      }
                 if ( !dcr0 )      if ( !dcr0 )
                         *pr = 0;        *pr = 0;
                 else {      else {
                         NEXT(dcr) = 0;        NEXT(dcr) = 0;
                         MKP(VR(p),dcr0,*pr);        MKP(VR(p),dcr0,*pr);
                 }      }
         }    }
 }  }
   
 void cm2p(mod,m,p,pr)  void cm2p(Q mod,Q m,P p,P *pr)
 Q mod,m;  
 P p,*pr;  
 {  {
         P t;    P t;
         DCP dc,dcr,dcr0;    DCP dc,dcr,dcr0;
   
         if ( !p )    if ( !p )
                 *pr = 0;      *pr = 0;
         else if ( NUM(p) )    else if ( NUM(p) )
                 rem2q((Q)p,mod,m,(Q *)pr);      rem2q((Q)p,mod,m,(Q *)pr);
         else {    else {
                 for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {      for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
                         cm2p(mod,m,COEF(dc),&t);        cm2p(mod,m,COEF(dc),&t);
                         if ( t ) {        if ( t ) {
                                 NEXTDC(dcr0,dcr);          NEXTDC(dcr0,dcr);
                                 DEG(dcr) = DEG(dc);          DEG(dcr) = DEG(dc);
                                 COEF(dcr) = t;          COEF(dcr) = t;
                         }        }
                 }      }
                 if ( !dcr0 )      if ( !dcr0 )
                         *pr = 0;        *pr = 0;
                 else {      else {
                         NEXT(dcr) = 0;        NEXT(dcr) = 0;
                         MKP(VR(p),dcr0,*pr);        MKP(VR(p),dcr0,*pr);
                 }      }
         }    }
 }  }
   
 void addm2q(mod,mod2,n1,n2,nr)  void addm2q(Q mod,Q mod2,Q n1,Q n2,Q *nr)
 Q mod,mod2;  
 Q n1,n2,*nr;  
 {  {
         Q t;    Q t;
   
         addq(n1,n2,&t);    addq(n1,n2,&t);
         if ( !t )    if ( !t )
                 *nr = 0;      *nr = 0;
         else    else
                 rem2q(t,mod,mod2,nr);      rem2q(t,mod,mod2,nr);
 }  }
   
 void subm2q(mod,mod2,n1,n2,nr)  void subm2q(Q mod,Q mod2,Q n1,Q n2,Q *nr)
 Q mod,mod2;  
 Q n1,n2,*nr;  
 {  {
         Q t;    Q t;
   
         subq(n1,n2,&t);    subq(n1,n2,&t);
         if ( !t )    if ( !t )
                 *nr = 0;      *nr = 0;
         else    else
                 rem2q(t,mod,mod2,nr);      rem2q(t,mod,mod2,nr);
 }  }
   
 void mulm2q(mod,mod2,n1,n2,nr)  void mulm2q(Q mod,Q mod2,Q n1,Q n2,Q *nr)
 Q mod,mod2;  
 Q n1,n2,*nr;  
 {  {
         Q t;    Q t;
   
         mulq(n1,n2,&t);    mulq(n1,n2,&t);
         if ( !t )    if ( !t )
                 *nr = 0;      *nr = 0;
         else    else
                 rem2q(t,mod,mod2,nr);      rem2q(t,mod,mod2,nr);
 }  }
   
 void rem2q(n,m,m2,nr)  void rem2q(Q n,Q m,Q m2,Q *nr)
 Q n,m,m2,*nr;  
 {  {
         N q,r,s;    N q,r,s;
         int sgn;    int sgn;
   
         divn(NM(n),NM(m),&q,&r);    divn(NM(n),NM(m),&q,&r);
         if ( !r )    if ( !r )
                 *nr = 0;      *nr = 0;
         else {    else {
                 sgn = cmpn(r,NM(m2));      sgn = cmpn(r,NM(m2));
                 if ( sgn > 0 ) {      if ( sgn > 0 ) {
                         subn(NM(m),r,&s); NTOQ(s,-1*SGN(n),*nr);        subn(NM(m),r,&s); NTOQ(s,-1*SGN(n),*nr);
                 } else      } else
                         NTOQ(r,SGN(n),*nr);        NTOQ(r,SGN(n),*nr);
         }    }
 }  }
   
 void exthp(vl,p,d,pr)  /*
 VL vl;    extract d-homogeneous part with respect to vl - {v}
 P p;  */
 int d;  
 P *pr;  void exthpc_generic(VL vl,P p,int d,V v,P *pr)
 {  {
         P t,t1,a,w,x,xd;    P w,x,t,t1,a,xd;
         DCP dc;    V v0;
     DCP dc;
   
         if ( d < 0 )    if ( d < 0 || !p )
                 *pr = 0;      *pr = 0;
         else if ( NUM(p) )    else if ( NUM(p) )
                 if ( d == 0 )      if ( d == 0 )
                         *pr = p;        *pr = p;
                 else      else
                         *pr = 0;        *pr = 0;
         else {    else if ( v == VR(p) )
                 for ( MKV(VR(p),x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {      exthpc(vl,v,p,d,pr);
                         exthp(vl,COEF(dc),d - QTOS(DEG(dc)),&t);    else {
                         pwrp(vl,x,DEG(dc),&xd);      v0 = VR(p);
                         mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;      for ( MKV(v0,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
                 }        exthpc_generic(vl,COEF(dc),d-QTOS(DEG(dc)),v,&t);
                 *pr = w;        pwrp(vl,x,DEG(dc),&xd);
         }        mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
       }
       *pr = w;
     }
 }  }
   
 void exthpc(vl,v,p,d,pr)  void exthp(VL vl,P p,int d,P *pr)
 VL vl;  
 V v;  
 P p;  
 int d;  
 P *pr;  
 {  {
         P t,t1,a,w,x,xd;    P t,t1,a,w,x,xd;
         DCP dc;    DCP dc;
   
         if ( v != VR(p) )    if ( d < 0 )
                 exthp(vl,p,d,pr);      *pr = 0;
         else if ( d < 0 )    else if ( NUM(p) )
                 *pr = 0;      if ( d == 0 )
         else {        *pr = p;
                 for ( MKV(v,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {      else
                         exthp(vl,COEF(dc),d,&t);        *pr = 0;
                         pwrp(vl,x,DEG(dc),&xd);    else {
                         mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;      for ( MKV(VR(p),x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
                 }        exthp(vl,COEF(dc),d - QTOS(DEG(dc)),&t);
                 *pr = w;        pwrp(vl,x,DEG(dc),&xd);
         }        mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
       }
       *pr = w;
     }
 }  }
   
 void cbound(vl,p,b)  void exthpc(VL vl,V v,P p,int d,P *pr)
 VL vl;  
 P p;  
 Q *b;  
 {  {
         Q n,e,t,m;    P t,t1,a,w,x,xd;
         int k;    DCP dc;
   
         cmax(p,&n);    if ( v != VR(p) )
         addq(n,n,&m);      exthp(vl,p,d,pr);
     else if ( d < 0 )
       *pr = 0;
     else {
       for ( MKV(v,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
         exthp(vl,COEF(dc),d,&t);
         pwrp(vl,x,DEG(dc),&xd);
         mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
       }
       *pr = w;
     }
   }
   
         k = geldb(vl,p);  void cbound(VL vl,P p,Q *b)
         STOQ(3,t); STOQ(k,e);  {
     Q n,e,t,m;
     int k;
   
         pwrq(t,e,&n);    cmax(p,&n);
         mulq(m,n,b);    addq(n,n,&m);
   
     k = geldb(vl,p);
     STOQ(3,t); STOQ(k,e);
   
     pwrq(t,e,&n);
     mulq(m,n,b);
 }  }
   
 int geldb(vl,p)  int geldb(VL vl,P p)
 VL vl;  
 P p;  
 {  {
         int m;    int m;
   
         for ( m = 0; vl; vl = NEXT(vl) )    for ( m = 0; vl; vl = NEXT(vl) )
                 m += getdeg(vl->v,p);      m += getdeg(vl->v,p);
         return ( m );    return ( m );
 }  }
   
 int getdeg(v,p)  int getdeg(V v,P p)
 V v;  
 P p;  
 {  {
         int m;    int m,t;
         DCP dc;    DCP dc;
   
         if ( !p || NUM(p) )    if ( !p || NUM(p) )
                 return ( 0 );      return ( 0 );
         else if ( v == VR(p) )    else if ( v == VR(p) )
                 return ( deg(v,p) );      return ( deg(v,p) );
         else {    else {
                 for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) )      for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
                         m = MAX(m,getdeg(v,COEF(dc)));        t = getdeg(v,COEF(dc));
                 return ( m );        m = MAX(m,t);
         }      }
       return ( m );
     }
 }  }
   
 void cmax(p,b)  void cmax(P p,Q *b)
 P p;  
 Q *b;  
 {  {
         DCP dc;    DCP dc;
         Q m,m1;    Q m,m1;
         N tn;    N tn;
   
         if ( NUM(p) ) {    if ( NUM(p) ) {
                 tn = NM((Q)p);      tn = NM((Q)p);
                 NTOQ(tn,1,*b);      NTOQ(tn,1,*b);
         } else {    } else {
                 for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {      for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
                         cmax(COEF(dc),&m1);        cmax(COEF(dc),&m1);
                         if ( cmpq(m1,m) > 0 )        if ( cmpq(m1,m) > 0 )
                                 m = m1;          m = m1;
                 }      }
                 *b = m;      *b = m;
         }    }
 }  }
   
 int nextbin(vn,n)  int nextbin(VN vn,int n)
 VN vn;  
 int n;  
 {  {
         int tmp,i,carry;    int tmp,i,carry;
   
         if ( n == 0 )    if ( n == 0 )
                 return ( 1 );      return ( 1 );
   
         for ( i = n - 1, carry = 1; i >= 0; i-- ) {    for ( i = n - 1, carry = 1; i >= 0; i-- ) {
                 tmp =  vn[i].n + carry;      tmp =  vn[i].n + carry;
                 vn[i].n = tmp % 2;      vn[i].n = tmp % 2;
                 carry = tmp / 2;      carry = tmp / 2;
         }    }
         return ( carry );    return ( carry );
 }  }
   
 void mulsgn(vn,vnt,n,vn1)  void mulsgn(VN vn,VN vnt,int n,VN vn1)
 VN vn,vnt,vn1;  
 int n;  
 {  {
         int i;    int i;
   
         for ( i = 0; vn[i].v; i++ )    for ( i = 0; vn[i].v; i++ )
                 vn1[i].n = vn[i].n;      vn1[i].n = vn[i].n;
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 if ( vnt[i].n )      if ( vnt[i].n )
                         vn1[(int)vnt[i].v].n *= -1;        vn1[(int)vnt[i].v].n *= -1;
 }  }
   
 void next(vn)  void next(VN vn)
 VN vn;  
 {  {
         int i,m,n,tmp,carry;    int i,m,n,tmp,carry;
   
         for ( m = 0, i = 0; vn[i].v; i++ )    for ( m = 0, i = 0; vn[i].v; i++ )
                 m = MAX(m,ABS(vn[i].n));      m = MAX(m,ABS(vn[i].n));
         if ( m == 0 ) {    if ( m == 0 ) {
                 vn[--i].n = 1;      vn[--i].n = 1;
                 return;      return;
         }    }
         for ( carry = 1, m++, i--; ( i >= 0 ) && carry; i-- ) {    for ( carry = 1, m++, i--; ( i >= 0 ) && carry; i-- ) {
                 tmp = vn[i].n + carry;      tmp = vn[i].n + carry;
                 vn[i].n = tmp % m;      vn[i].n = tmp % m;
                 carry = tmp / m;      carry = tmp / m;
         }    }
         if ( ( i == -1 ) && carry ) {    if ( ( i == -1 ) && carry ) {
                 for ( i = 0; vn[i].v; i++ )      for ( i = 0; vn[i].v; i++ )
                         vn[i].n = 0;        vn[i].n = 0;
                 vn[--i].n = m;      vn[--i].n = m;
         } else {    } else {
                 for ( n = 0, i = 0; vn[i].v; i++ )      for ( n = 0, i = 0; vn[i].v; i++ )
                         n = MAX(n,ABS(vn[i].n));        n = MAX(n,ABS(vn[i].n));
                 if ( n < m - 1 )      if ( n < m - 1 )
                         vn[--i].n = m - 1;        vn[--i].n = m - 1;
         }    }
 }  }
   
 void clctv(vl,p,nvlp)  void clctv(VL vl,P p,VL *nvlp)
 VL vl;  
 P p;  
 VL *nvlp;  
 {  {
         int i,n;    int i,n;
         VL tvl;    VL tvl;
         VN tvn;    VN tvn;
   
         if ( !p || NUM(p) ) {    if ( !p || NUM(p) ) {
                 *nvlp = 0;      *nvlp = 0;
                 return;      return;
         }    }
   
         for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );    for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
         tvn = (VN) ALLOCA((n+1)*sizeof(struct oVN));    tvn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
         for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) {    for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) {
                 tvn[i].v = tvl->v;      tvn[i].v = tvl->v;
                 tvn[i].n = 0;      tvn[i].n = 0;
         }    }
   
         markv(tvn,n,p);    markv(tvn,n,p);
         vntovl(tvn,n,nvlp);    vntovl(tvn,n,nvlp);
 }  }
   
 void markv(vn,n,p)  void markv(VN vn,int n,P p)
 VN vn;  
 int n;  
 P p;  
 {  {
         V v;    V v;
         DCP dc;    DCP dc;
         int i;    int i;
   
         if ( NUM(p) )    if ( NUM(p) )
                 return;      return;
         v = VR(p);    v = VR(p);
         for ( i = 0, v = VR(p); i < n; i++ )    for ( i = 0, v = VR(p); i < n; i++ )
                 if ( v == vn[i].v ) {      if ( v == vn[i].v ) {
                         vn[i].n = 1;        vn[i].n = 1;
                         break;        break;
                 }      }
         for ( dc = DC(p); dc; dc = NEXT(dc) )    for ( dc = DC(p); dc; dc = NEXT(dc) )
                 markv(vn,n,COEF(dc));      markv(vn,n,COEF(dc));
 }  }
   
 void vntovl(vn,n,vlp)  void vntovl(VN vn,int n,VL *vlp)
 VN vn;  
 int n;  
 VL *vlp;  
 {  {
         int i;    int i;
         VL tvl,tvl0;    VL tvl,tvl0;
   
         for ( i = 0, tvl0 = 0; ; ) {    for ( i = 0, tvl0 = 0; ; ) {
                 while ( ( i < n ) && ( vn[i].n == 0 ) ) i++;      while ( ( i < n ) && ( vn[i].n == 0 ) ) i++;
                 if ( i == n )      if ( i == n )
                         break;        break;
                 else {      else {
                         if ( !tvl0 ) {        if ( !tvl0 ) {
                                 NEWVL(tvl0);          NEWVL(tvl0);
                                 tvl = tvl0;          tvl = tvl0;
                         } else {        } else {
                                 NEWVL(NEXT(tvl));          NEWVL(NEXT(tvl));
                                 tvl = NEXT(tvl);          tvl = NEXT(tvl);
                         }        }
                         tvl->v = vn[i++].v;        tvl->v = vn[i++].v;
                 }      }
         }    }
         if ( tvl0 )    if ( tvl0 )
                 NEXT(tvl) = 0;      NEXT(tvl) = 0;
         *vlp = tvl0;    *vlp = tvl0;
 }  }
   
 int dbound(v,f)  int dbound(V v,P f)
 V v;  
 P f;  
 {  {
         int m;    int m,t;
         DCP dc;    DCP dc;
   
         if ( !f )    if ( !f )
                 return ( -1 );      return ( -1 );
         else if ( v != VR(f) )    else if ( v != VR(f) )
                 return homdeg(f);      return homdeg(f);
         else {    else {
                 for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) )      for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
                         m = MAX(m,homdeg(COEF(dc)));        t = homdeg(COEF(dc));
                 return ( m );        m = MAX(m,t);
         }      }
       return ( m );
     }
 }  }
   
 int homdeg(f)  int homdeg(P f)
 P f;  
 {  {
         int m;    int m,t;
         DCP dc;    DCP dc;
   
         if ( !f )    if ( !f )
                 return ( -1 );      return ( -1 );
         else if ( NUM(f) )    else if ( NUM(f) )
                 return( 0 );      return( 0 );
         else {    else {
                 for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) )      for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
                         m = MAX(m,QTOS(DEG(dc))+homdeg(COEF(dc)));        t = QTOS(DEG(dc))+homdeg(COEF(dc));
                 return ( m );        m = MAX(m,t);
         }      }
       return ( m );
     }
 }  }
   
 int minhomdeg(f)  int minhomdeg(P f)
 P f;  
 {  {
         int m;    int m,t;
         DCP dc;    DCP dc;
   
         if ( !f )    if ( !f )
                 return ( -1 );      return ( -1 );
         else if ( NUM(f) )    else if ( NUM(f) )
                 return( 0 );      return( 0 );
         else {    else {
                 for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) )      for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) ) {
                         m = MIN(m,QTOS(DEG(dc))+minhomdeg(COEF(dc)));        t = QTOS(DEG(dc))+minhomdeg(COEF(dc));
                 return ( m );        m = MIN(m,t);
         }      }
       return ( m );
     }
 }  }
   
 void adjc(vl,f,a,lc0,q,fr,ar)  void adjc(VL vl,P f,P a,P lc0,Q q,P *fr,P *ar)
 VL vl;  
 P f,a,lc0;  
 Q q;  
 P *fr,*ar;  
 {  {
         P m,m1;    P m,m1;
         Q t;    Q t;
   
         invl((Q)LC(f),q,&t);    invl((Q)LC(f),q,&t);
         mulq((Q)lc0,t,(Q *)&m);    mulq((Q)lc0,t,(Q *)&m);
         mulpq(f,m,&m1); cmp(q,m1,fr);    mulpq(f,m,&m1); cmp(q,m1,fr);
         invl((Q)m,q,&t);    invl((Q)m,q,&t);
         mulpq(a,(P)t,&m1);    mulpq(a,(P)t,&m1);
         cmp(q,m1,ar);    cmp(q,m1,ar);
 }  }
   
 #if 1  #if 1
 void affinemain(vl,p,v0,n,pl,pr)  void affinemain(VL vl,P p,V v0,int n,P *pl,P *pr)
 VL vl;  
 V v0;  
 int n;  
 P *pl;  
 P p;  
 P *pr;  
 {  {
         P x,t,m,c,s,a;    P x,t,m,c,s,a;
         DCP dc;    DCP dc;
         Q d;    Q d;
   
         if ( !p )    if ( !p )
                 *pr = 0;      *pr = 0;
         else if ( NUM(p) )    else if ( NUM(p) )
                 *pr = p;      *pr = p;
         else if ( VR(p) != v0 ) {    else if ( VR(p) != v0 ) {
                 MKV(VR(p),x);      MKV(VR(p),x);
                 for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {      for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
                         affinemain(vl,COEF(dc),v0,n,pl,&t);        affinemain(vl,COEF(dc),v0,n,pl,&t);
                         if ( DEG(dc) ) {        if ( DEG(dc) ) {
                                 pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);          pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
                                 addp(vl,m,c,&a); c = a;          addp(vl,m,c,&a); c = a;
                         } else {        } else {
                                 addp(vl,t,c,&a); c = a;          addp(vl,t,c,&a); c = a;
                         }        }
                 }      }
                 *pr = c;      *pr = c;
         } else {    } else {
                 dc = DC(p);      dc = DC(p);
                 c = COEF(dc);      c = COEF(dc);
                 for ( d = DEG(dc), dc = NEXT(dc);      for ( d = DEG(dc), dc = NEXT(dc);
                         dc; d = DEG(dc), dc = NEXT(dc) ) {        dc; d = DEG(dc), dc = NEXT(dc) ) {
                                 mulp(vl,pl[QTOS(d)-QTOS(DEG(dc))],c,&m);          mulp(vl,pl[QTOS(d)-QTOS(DEG(dc))],c,&m);
                                 addp(vl,m,COEF(dc),&c);          addp(vl,m,COEF(dc),&c);
                 }      }
                 if ( d ) {      if ( d ) {
                         mulp(vl,pl[QTOS(d)],c,&m); c = m;        mulp(vl,pl[QTOS(d)],c,&m); c = m;
                 }      }
                 *pr = c;      *pr = c;
         }    }
 }  }
 #endif  #endif
   
 #if 0  #if 0
 affine(vl,p,vn,r)  affine(VL vl,P p,VN vn,P *r)
 VL vl;  
 P p;  
 VN vn;  
 P *r;  
 {  {
         int n,d,d1,i;    int n,d,d1,i;
         Q *q;    Q *q;
         Q **bc;    Q **bc;
   
         if ( !p || NUM(p) )    if ( !p || NUM(p) )
                 *r = p;      *r = p;
         else {    else {
                 for ( i = 0, d = 0; vn[i].v; i++ )      for ( i = 0, d = 0; vn[i].v; i++ )
                         d1 = getdeg(vn[i].v,p), d = MAX(d,d1);        d1 = getdeg(vn[i].v,p), d = MAX(d,d1);
                 W_CALLOC(d+1,Q *,bc);      W_CALLOC(d+1,Q *,bc);
                 for ( i = 0; i <= d; i++ )      for ( i = 0; i <= d; i++ )
                         W_CALLOC(i+1,Q,q), mkbc(i,q), bc[i] = q;        W_CALLOC(i+1,Q,q), mkbc(i,q), bc[i] = q;
                 afmain(vl,bc,p,vn,r);      afmain(vl,bc,p,vn,r);
         }    }
 }  }
   
 afmain(vl,bc,p,vn,r)  afmain(VL vl,Q **bc,P p,VN vn,P *r)
 VL vl;  
 Q **bc;  
 P p;  
 VN vn;  
 P *r;  
 {  {
         P t,s,u;    P t,s,u;
         P *c,*rc;    P *c,*rc;
         Q *q;    Q *q;
         DCP dc;    DCP dc;
         int n,i,j;    int n,i,j;
   
         if ( !p || NUM(p) || !vn[0].v )    if ( !p || NUM(p) || !vn[0].v )
                 *r = p;      *r = p;
         else if ( vn[0].v != VR(p) ) {    else if ( vn[0].v != VR(p) ) {
                 for ( i = 0; vn[i].v && (vn[i].v != VR(p)); i++ );      for ( i = 0; vn[i].v && (vn[i].v != VR(p)); i++ );
                 if ( vn[i].v )      if ( vn[i].v )
                         afmain(vl,bc,p,vn+i,r);        afmain(vl,bc,p,vn+i,r);
                 else {      else {
                         n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);        n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);
                         for ( dc = DC(p); dc; dc = NEXT(dc) )        for ( dc = DC(p); dc; dc = NEXT(dc) )
                                 afmain(vl,bc,COEF(dc),vn,&c[QTOS(DEG(dc))]);          afmain(vl,bc,COEF(dc),vn,&c[QTOS(DEG(dc))]);
                         plisttop(c,VR(p),n,r);        plisttop(c,VR(p),n,r);
                 }      }
         } else {    } else {
                 n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);      n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);
                 W_CALLOC(n+1,P,rc); W_CALLOC(n+1,Q,q);      W_CALLOC(n+1,P,rc); W_CALLOC(n+1,Q,q);
                 for ( dc = DC(p); dc; dc = NEXT(dc) )      for ( dc = DC(p); dc; dc = NEXT(dc) )
                         afmain(vl,bc,COEF(dc),vn+1,&c[QTOS(DEG(dc))]);        afmain(vl,bc,COEF(dc),vn+1,&c[QTOS(DEG(dc))]);
                 if ( !vn[0].n )      if ( !vn[0].n )
                         bcopy(c,rc,sizeof(P)*(n+1));        bcopy(c,rc,sizeof(P)*(n+1));
                 else {      else {
                         for ( q[0] = ONE, STOQ(vn[0].n,q[1]), i = 2; i <= n; i++ )        for ( q[0] = ONE, STOQ(vn[0].n,q[1]), i = 2; i <= n; i++ )
                                 mulq(q[i-1],q[1],&q[i]);          mulq(q[i-1],q[1],&q[i]);
                         for ( j = 0; j <= n; rc[j] = t, j++ )        for ( j = 0; j <= n; rc[j] = t, j++ )
                                 for ( i = j, t = 0; i <= n; i++ )          for ( i = j, t = 0; i <= n; i++ )
                                         if ( c[i] )            if ( c[i] )
                                                 mulq(bc[i][j],q[i-j],&s), mulp(vl,s,c[i],&u),              mulq(bc[i][j],q[i-j],&s), mulp(vl,s,c[i],&u),
                                                 addp(CO,u,t,&s), t = s;              addp(CO,u,t,&s), t = s;
                 }      }
                 plisttop(rc,VR(p),n,r);      plisttop(rc,VR(p),n,r);
         }    }
 }  }
 #endif  #endif
   
 void restore(vl,f,vn,fr)  void restore(VL vl,P f,VN vn,P *fr)
 VL vl;  
 P f;  
 VN vn;  
 P *fr;  
 {  {
         int i;    int i;
         P vv,g,g1,t;    P vv,g,g1,t;
         Q s;    Q s;
   
         g = f;    g = f;
         for ( i = 0; vn[i].v; i++ ) {    for ( i = 0; vn[i].v; i++ ) {
                 MKV(vn[i].v,t);      MKV(vn[i].v,t);
                 if ( vn[i].n ) {      if ( vn[i].n ) {
                         STOQ(-vn[i].n,s);        STOQ(-vn[i].n,s);
                         addp(vl,t,(P)s,&vv);        addp(vl,t,(P)s,&vv);
                 } else      } else
                         vv = t;        vv = t;
   
                 substp(vl,g,vn[i].v,vv,&g1); g = g1;      substp(vl,g,vn[i].v,vv,&g1); g = g1;
         }    }
         *fr = g;    *fr = g;
 }  }
   
 void mergev(vl,vl1,vl2,nvlp)  void mergev(VL vl,VL vl1,VL vl2,VL *nvlp)
 VL vl,vl1,vl2,*nvlp;  
 {  {
         int i,n;    int i,n;
         VL tvl;    VL tvl;
         VN vn;    VN vn;
   
         if ( !vl1 ) {    if ( !vl1 ) {
                 *nvlp = vl2; return;      *nvlp = vl2; return;
         } else if ( !vl2 ) {    } else if ( !vl2 ) {
                 *nvlp = vl1; return;      *nvlp = vl1; return;
         }    }
         for ( i = 0, tvl = vl; tvl; i++, tvl = NEXT(tvl) );    for ( i = 0, tvl = vl; tvl; i++, tvl = NEXT(tvl) );
         n = i;    n = i;
         W_CALLOC(n,struct oVN,vn);    W_CALLOC(n,struct oVN,vn);
         for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ )    for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ )
                 vn[i].v = tvl->v;      vn[i].v = tvl->v;
         for ( i = 0, tvl = vl1; (i < n) && tvl; i++, tvl = NEXT(tvl)) {    for ( i = 0, tvl = vl1; (i < n) && tvl; i++, tvl = NEXT(tvl)) {
                 while ( ( i < n ) && ( vn[i].v != tvl->v ) )      while ( ( i < n ) && ( vn[i].v != tvl->v ) )
                         i++;        i++;
                 if ( i == n )      if ( i == n )
                         break;        break;
                 else      else
                         vn[i].n = 1;        vn[i].n = 1;
         }    }
         for ( i = 0, tvl = vl2; (i < n) && tvl; i++, tvl = NEXT(tvl) ) {    for ( i = 0, tvl = vl2; (i < n) && tvl; i++, tvl = NEXT(tvl) ) {
                 while ( ( i < n ) && ( vn[i].v != tvl->v ) )      while ( ( i < n ) && ( vn[i].v != tvl->v ) )
                         i++;        i++;
                 if ( i == n )      if ( i == n )
                         break;        break;
                 else      else
                         vn[i].n = 1;        vn[i].n = 1;
         }    }
         vntovl(vn,n,nvlp);    vntovl(vn,n,nvlp);
 }  }
   
 #if 0  #if 0
 void substvp(vl,f,vn,g)  void substvp(VL vl,P f,VN vn,P *g)
 VL vl;  
 P f;  
 VN vn;  
 P *g;  
 {  {
         V v;    V v;
         int i;    int i;
         P h,h1;    P h,h1;
         Q t;    Q t;
   
         h = f;    h = f;
         for ( i = 0; v = vn[i].v; i++ ) {    for ( i = 0; v = vn[i].v; i++ ) {
                 STOQ(vn[i].n,t);      STOQ(vn[i].n,t);
                 substp(vl,h,v,(P)t,&h1); h = h1;      substp(vl,h,v,(P)t,&h1); h = h1;
         }    }
         *g = h;    *g = h;
 }  }
   
 void affine(vl,f,vn,fr)  void affine(VL vl,P f,VN vn,P *fr)
 VL vl;  
 P f;  
 VN vn;  
 P *fr;  
 {  {
         int i,j,n;    int i,j,n;
         P vv,g,g1,t,u;    P vv,g,g1,t,u;
         Q s;    Q s;
         int *dlist;    int *dlist;
         P **plist;    P **plist;
   
         for ( n = 0; vn[n].v; n++);    for ( n = 0; vn[n].v; n++);
         dlist = (int *)ALLOCA((n+1)*sizeof(int));    dlist = (int *)ALLOCA((n+1)*sizeof(int));
         plist = (P **)ALLOCA((n+1)*sizeof(P *));    plist = (P **)ALLOCA((n+1)*sizeof(P *));
         for ( i = 0; vn[i].v; i++ ) {    for ( i = 0; vn[i].v; i++ ) {
                 if ( !vn[i].n )      if ( !vn[i].n )
                         continue;        continue;
                 dlist[i] = getdeg(vn[i].v,f);      dlist[i] = getdeg(vn[i].v,f);
                 plist[i] = (P *)ALLOCA((dlist[i]+1)*sizeof(P));      plist[i] = (P *)ALLOCA((dlist[i]+1)*sizeof(P));
   
                 MKV(vn[i].v,t);      MKV(vn[i].v,t);
                 if ( vn[i].n ) {      if ( vn[i].n ) {
                         STOQ(vn[i].n,s); addp(vl,t,(P)s,&vv);        STOQ(vn[i].n,s); addp(vl,t,(P)s,&vv);
                 } else      } else
                         vv = t;        vv = t;
   
                 for ( j = 0, t = (P)ONE; j < dlist[i]; j++ ) {      for ( j = 0, t = (P)ONE; j < dlist[i]; j++ ) {
                         plist[i][j] = t;        plist[i][j] = t;
                         mulp(vl,t,vv,&u);        mulp(vl,t,vv,&u);
                         t = u;        t = u;
                 }      }
                 plist[i][j] = t;      plist[i][j] = t;
         }    }
   
         g = f;    g = f;
         for ( i = 0; vn[i].v; i++ ) {    for ( i = 0; vn[i].v; i++ ) {
                 if ( !vn[i].n )      if ( !vn[i].n )
                         continue;        continue;
                 affinemain(vl,g,vn[i].v,dlist[i],plist[i],&g1); g = g1;      affinemain(vl,g,vn[i].v,dlist[i],plist[i],&g1); g = g1;
         }    }
         *fr = g;    *fr = g;
 }  }
 #endif  #endif
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.9
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>