=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/E.c,v retrieving revision 1.8 retrieving revision 1.9 diff -u -p -r1.8 -r1.9 --- OpenXM_contrib2/asir2000/engine/E.c 2001/10/09 01:36:09 1.8 +++ OpenXM_contrib2/asir2000/engine/E.c 2018/03/29 01:32:51 1.9 @@ -45,955 +45,955 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM: OpenXM_contrib2/asir2000/engine/E.c,v 1.7 2001/06/07 04:54:40 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/engine/E.c,v 1.8 2001/10/09 01:36:09 noro Exp $ */ #include "ca.h" 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) { - P g1,h1,a1,b1; - N qn; - Q q2; + P g1,h1,a1,b1; + N qn; + Q 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); - henmvmain(vl,vn,f,g1,h1,b1,a1,lg,lh,q,q2,k,gp,hp); + 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); + henmvmain(vl,vn,f,g1,h1,b1,a1,lg,lh,q,q2,k,gp,hp); } 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) { - V v; - int n,i,j; - int *md; - P x,m,m1,c,q,r,a,s,u,ff,f0,f1; - P w0,w1,cf,cfi,t,q1,dvr; - P *c0,*c1; - P *f0h,*f1h; + V v; + int n,i,j; + int *md; + P x,m,m1,c,q,r,a,s,u,ff,f0,f1; + P w0,w1,cf,cfi,t,q1,dvr; + P *c0,*c1; + P *f0h,*f1h; - v = VR(f); n = deg(v,f); MKV(v,x); - c0 = (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); - cm2p(mod,mod2,gi0,&c0[0]); cm2p(mod,mod2,gi1,&c1[0]); - for ( i = 1; i <= n; i++ ) { - mulp(vl,x,c1[i-1],&m); divsrp(vl,m,dvr,&q,&r); mulp(vl,q,c,&q1); - cm2p(mod,mod2,r,&c1[i]); - mulp(vl,x,c0[i-1],&m); mulp(vl,q1,fi0,&m1); addp(vl,m,m1,&a); - cm2p(mod,mod2,a,&c0[i]); - } - affine(vl,f,vn,&t); cm2p(mod,mod2,t,&ff); - for ( i = 0; vn[i].v; i++ ); - md = ( int *) ALLOCA((i+1)*sizeof(int)); - for ( i = 0; vn[i].v; i++ ) - md[i] = getdeg(vn[i].v,ff); - cm2p(mod,mod2,fi0,&f0); affine(vl,l0,vn,&t); - if ( NUM(f0) ) - cm2p(mod,mod2,t,&f0); - else - cm2p(mod,mod2,t,&COEF(DC(f0))); - cm2p(mod,mod2,fi1,&f1); affine(vl,l1,vn,&t); - if ( NUM(f1) ) - cm2p(mod,mod2,t,&f1); - else - cm2p(mod,mod2,t,&COEF(DC(f1))); - W_CALLOC(k,P,f0h); W_CALLOC(k,P,f1h); - for ( i = 0; i <= k; i++ ) { - exthpc(vl,v,f0,i,&f0h[i]); exthpc(vl,v,f1,i,&f1h[i]); - } - for ( j = 1; j <= k; j++ ) { - for ( i = 0; vn[i].v; i++ ) - if ( getdeg(vn[i].v,f0)+getdeg(vn[i].v,f1) > md[i] ) - goto END; - for ( i = 0, s = 0; i <= j; i++ ) { - mulp(vl,f0h[i],f1h[j-i],&t); addp(vl,s,t,&u); s = u; - } - cm2p(mod,mod2,s,&t); - exthpc(vl,v,ff,j,&u); subm2p(vl,mod,mod2,u,t,&cf); - for ( i = 0, w0 = 0, w1 = 0; i <= n; i++ ) { - if ( !cf ) - cfi = 0; - else if ( VR(cf) == v ) - coefp(cf,i,&cfi); - else if ( i == 0 ) - cfi = cf; - else - cfi = 0; - if ( cfi ) { - 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; - } - } - 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; - 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; - if ( !t ) { - restore(vl,f0,vn,&s); cm2p(mod,mod2,s,&t); - if ( divtpz(vl,f,t,&s) ) { - *fr0 = t; *fr1 = s; - return; - } - } - if ( !u ) { - restore(vl,f1,vn,&s); cm2p(mod,mod2,s,&t); - if ( divtpz(vl,f,t,&s) ) { - *fr0 = s; *fr1 = t; - return; - } - } - } + v = VR(f); n = deg(v,f); MKV(v,x); + c0 = (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); + cm2p(mod,mod2,gi0,&c0[0]); cm2p(mod,mod2,gi1,&c1[0]); + for ( i = 1; i <= n; i++ ) { + mulp(vl,x,c1[i-1],&m); divsrp(vl,m,dvr,&q,&r); mulp(vl,q,c,&q1); + cm2p(mod,mod2,r,&c1[i]); + mulp(vl,x,c0[i-1],&m); mulp(vl,q1,fi0,&m1); addp(vl,m,m1,&a); + cm2p(mod,mod2,a,&c0[i]); + } + affine(vl,f,vn,&t); cm2p(mod,mod2,t,&ff); + for ( i = 0; vn[i].v; i++ ); + md = ( int *) ALLOCA((i+1)*sizeof(int)); + for ( i = 0; vn[i].v; i++ ) + md[i] = getdeg(vn[i].v,ff); + cm2p(mod,mod2,fi0,&f0); affine(vl,l0,vn,&t); + if ( NUM(f0) ) + cm2p(mod,mod2,t,&f0); + else + cm2p(mod,mod2,t,&COEF(DC(f0))); + cm2p(mod,mod2,fi1,&f1); affine(vl,l1,vn,&t); + if ( NUM(f1) ) + cm2p(mod,mod2,t,&f1); + else + cm2p(mod,mod2,t,&COEF(DC(f1))); + W_CALLOC(k,P,f0h); W_CALLOC(k,P,f1h); + for ( i = 0; i <= k; i++ ) { + exthpc(vl,v,f0,i,&f0h[i]); exthpc(vl,v,f1,i,&f1h[i]); + } + for ( j = 1; j <= k; j++ ) { + for ( i = 0; vn[i].v; i++ ) + if ( getdeg(vn[i].v,f0)+getdeg(vn[i].v,f1) > md[i] ) + goto END; + for ( i = 0, s = 0; i <= j; i++ ) { + mulp(vl,f0h[i],f1h[j-i],&t); addp(vl,s,t,&u); s = u; + } + cm2p(mod,mod2,s,&t); + exthpc(vl,v,ff,j,&u); subm2p(vl,mod,mod2,u,t,&cf); + for ( i = 0, w0 = 0, w1 = 0; i <= n; i++ ) { + if ( !cf ) + cfi = 0; + else if ( VR(cf) == v ) + coefp(cf,i,&cfi); + else if ( i == 0 ) + cfi = cf; + else + cfi = 0; + if ( cfi ) { + 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; + } + } + 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; + 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; + if ( !t ) { + restore(vl,f0,vn,&s); cm2p(mod,mod2,s,&t); + if ( divtpz(vl,f,t,&s) ) { + *fr0 = t; *fr1 = s; + return; + } + } + if ( !u ) { + restore(vl,f1,vn,&s); cm2p(mod,mod2,s,&t); + if ( divtpz(vl,f,t,&s) ) { + *fr0 = s; *fr1 = t; + return; + } + } + } END: - restore(vl,f0,vn,&t); cm2p(mod,mod2,t,fr0); - restore(vl,f1,vn,&t); cm2p(mod,mod2,t,fr1); + restore(vl,f0,vn,&t); cm2p(mod,mod2,t,fr0); + restore(vl,f1,vn,&t); cm2p(mod,mod2,t,fr1); } /* - 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)) ) + 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)) ) */ 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) { - N qn; - Q q,qq,q2; - int n,i; - UM wg0,wg1,wf0,wf1; - P f0,f1,g0,g1,m,m1,s,c,d,q1,r1,rm,rm1,a,a1,a2; + N qn; + Q q,qq,q2; + int n,i; + UM wg0,wg1,wf0,wf1; + P f0,f1,g0,g1,m,m1,s,c,d,q1,r1,rm,rm1,a,a1,a2; - n = UDEG(f); - wf0 = W_UMALLOC(n); wf1 = W_UMALLOC(n); - wg0 = W_UMALLOC(n); wg1 = W_UMALLOC(n); - cpyum(fi0,wf0); cpyum(fi1,wf1); - eucum(p,wf0,wf1,wg1,wg0); - umtop(VR(f),fi0,&f0); umtop(VR(f),fi1,&f1); - umtop(VR(f),wg0,&g0); umtop(VR(f),wg1,&g1); + n = UDEG(f); + wf0 = W_UMALLOC(n); wf1 = W_UMALLOC(n); + wg0 = W_UMALLOC(n); wg1 = W_UMALLOC(n); + cpyum(fi0,wf0); cpyum(fi1,wf1); + eucum(p,wf0,wf1,wg1,wg0); + umtop(VR(f),fi0,&f0); umtop(VR(f),fi1,&f1); + umtop(VR(f),wg0,&g0); umtop(VR(f),wg1,&g1); - STOQ(p,q); divin(NM(q),2,&qn); NTOQ(qn,1,q2); - for ( i = 1; i < k; i++ ) { + STOQ(p,q); divin(NM(q),2,&qn); NTOQ(qn,1,q2); + for ( i = 1; i < k; i++ ) { #if 0 - mulp(CO,i1,g0,&m); mulp(CO,i0,g1,&m1); addp(CO,m,m1,&a); - if ( NUM(a) ) { - for ( STOQ(p,q), j = 1; j < k; j++ ) { - mulq(q,q,&qq); q = qq; - } - f0 = i0; f1 = i1; - invl(a,q,&qq); - mulp(CO,g0,qq,&s); g0 = s; mulp(CO,g1,qq,&s); g1 = s; - break; - } + mulp(CO,i1,g0,&m); mulp(CO,i0,g1,&m1); addp(CO,m,m1,&a); + if ( NUM(a) ) { + for ( STOQ(p,q), j = 1; j < k; j++ ) { + mulq(q,q,&qq); q = qq; + } + f0 = i0; f1 = i1; + invl(a,q,&qq); + mulp(CO,g0,qq,&s); g0 = s; mulp(CO,g1,qq,&s); g1 = s; + break; + } #endif - /* c = ((f - f0*f1)/q) mod q; - q1 = (c*g1) / f1; - r1 = (c*g1) mod f1; - f1 += (r1 mod q)*q; - f0 += ((c*g0 + q1*f0) mod q)*q; + /* c = ((f - f0*f1)/q) mod q; + q1 = (c*g1) / f1; + r1 = (c*g1) mod f1; + f1 += (r1 mod q)*q; + f0 += ((c*g0 + q1*f0) mod q)*q; - d = ((1 - (f1*g0 + f0*g1))/q) mod q; - q1 = (d*g0) / f1; - r1 = (d*g0) mod f1; - g1 += (r1 mod q)*q; - g0 += ((c*g0 + q1*f0) mod q)*q; - q = q^2; - */ + d = ((1 - (f1*g0 + f0*g1))/q) mod q; + q1 = (d*g0) / f1; + r1 = (d*g0) mod f1; + g1 += (r1 mod q)*q; + g0 += ((c*g0 + q1*f0) mod q)*q; + q = q^2; + */ - /* c = ((f - f0*f1)/q) mod q */ - mulp(CO,f0,f1,&m); subp(CO,f,m,&s); - divcp(s,q,&m); cm2p(q,q2,m,&c); + /* c = ((f - f0*f1)/q) mod q */ + mulp(CO,f0,f1,&m); subp(CO,f,m,&s); + divcp(s,q,&m); cm2p(q,q2,m,&c); - /* q1 = (c*g1) / f1; r1 = (c*g1) mod f1; */ - mulp(CO,c,g1,&m); cm2p(q,q2,m,&s); - udivpwm(q,s,f1,&q1,&r1); + /* q1 = (c*g1) / f1; r1 = (c*g1) mod f1; */ + mulp(CO,c,g1,&m); cm2p(q,q2,m,&s); + udivpwm(q,s,f1,&q1,&r1); - /* f1 = f1 + (r1 mod q)*q; */ - cm2p(q,q2,r1,&rm); - mulpq(rm,(P)q,&rm1); addp(CO,f1,rm1,&a); - f1 = a; + /* f1 = f1 + (r1 mod q)*q; */ + cm2p(q,q2,r1,&rm); + mulpq(rm,(P)q,&rm1); addp(CO,f1,rm1,&a); + f1 = a; - /* a1 = (c*g0 + q1*f0) mod q; */ - mulp(CO,c,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a); - cm2p(q,q2,a,&a1); + /* a1 = (c*g0 + q1*f0) mod q; */ + mulp(CO,c,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a); + cm2p(q,q2,a,&a1); - /* f0 = f0 + a1*q; */ - mulpq(a1,(P)q,&a2); - addp(CO,f0,a2,&a); - f0 = a; + /* f0 = f0 + a1*q; */ + mulpq(a1,(P)q,&a2); + addp(CO,f0,a2,&a); + f0 = a; - /* d = ((1 - (f1*g0 + f0*g1))/q) mod q; */ - mulp(CO,f1,g0,&m); mulp(CO,f0,g1,&m1); addp(CO,m,m1,&a); - subp(CO,(P)ONE,a,&s); - divcp(s,q,&m); cm2p(q,q2,m,&d); + /* d = ((1 - (f1*g0 + f0*g1))/q) mod q; */ + mulp(CO,f1,g0,&m); mulp(CO,f0,g1,&m1); addp(CO,m,m1,&a); + subp(CO,(P)ONE,a,&s); + divcp(s,q,&m); cm2p(q,q2,m,&d); - /* q1 = (d*g0) / f1; r1 = (d*g0) mod f1; */ - mulp(CO,d,g1,&m); udivpwm(q,m,f1,&q1,&r1); + /* q1 = (d*g0) / f1; r1 = (d*g0) mod f1; */ + mulp(CO,d,g1,&m); udivpwm(q,m,f1,&q1,&r1); - /* g1 = g1 + (r1 mod q )*q; */ - cm2p(q,q2,r1,&rm); - mulpq(rm,(P)q,&rm1); addp(CO,g1,rm1,&a); - g1 = a; + /* g1 = g1 + (r1 mod q )*q; */ + cm2p(q,q2,r1,&rm); + mulpq(rm,(P)q,&rm1); addp(CO,g1,rm1,&a); + g1 = a; - /* g0 = g0 + (c*g0 + q1*f0 mod q)*q; */ - 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); - addp(CO,g0,a2,&a); - g0 = a; + /* g0 = g0 + (c*g0 + q1*f0 mod q)*q; */ + 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); + addp(CO,g0,a2,&a); + g0 = a; - /* q = q^2; */ - mulq(q,q,&qq); - q = qq; - divin(NM(q),2,&qn); NTOQ(qn,1,q2); - } - *fr0p = f0; *fr1p = f1; *gr0p = g0; *gr1p = g1; *qrp = q; + /* q = q^2; */ + mulq(q,q,&qq); + q = qq; + divin(NM(q),2,&qn); NTOQ(qn,1,q2); + } + *fr0p = f0; *fr1p = f1; *gr0p = g0; *gr1p = g1; *qrp = q; } void henzq1(P g,P h,Q bound,P *gcp,P *hcp,Q *qp) { - V v; - Q f,q,q1; - Q u,t,a,b,s; - P c,c1; - P tg,th,tr; - UM wg,wh,ws,wt,wm; - int n,m,modres,mod,index,i; - P gc0,hc0; - P z,zz,zzz; + V v; + Q f,q,q1; + Q u,t,a,b,s; + P c,c1; + P tg,th,tr; + UM wg,wh,ws,wt,wm; + int n,m,modres,mod,index,i; + P gc0,hc0; + P z,zz,zzz; - v = VR(g); n=deg(v,g); m=deg(v,h); - norm(g,&a); norm(h,&b); - STOQ(m,u); pwrq(a,u,&t); - STOQ(n,u); pwrq(b,u,&s); - mulq(t,s,&u); + v = VR(g); n=deg(v,g); m=deg(v,h); + norm(g,&a); norm(h,&b); + STOQ(m,u); pwrq(a,u,&t); + STOQ(n,u); pwrq(b,u,&s); + mulq(t,s,&u); - factorial(n+m,&t); mulq(u,t,&s); - addq(s,s,&f); + factorial(n+m,&t); mulq(u,t,&s); + addq(s,s,&f); - wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n); - wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n); - wm = W_UMALLOC(m+n); + wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n); + wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n); + wm = W_UMALLOC(m+n); - for ( q = ONE, t = 0, c = 0, index = 0; ; ) { - mod = get_lprime(index++); - if ( !rem(NM((Q)LC(g)),mod) || !rem(NM((Q)LC(h)),mod) ) - continue; - ptomp(mod,g,&tg); ptomp(mod,h,&th); - srchump(mod,tg,th,&tr); - if ( !tr ) - continue; - else - modres = CONT((MQ)tr); + for ( q = ONE, t = 0, c = 0, index = 0; ; ) { + mod = get_lprime(index++); + if ( !rem(NM((Q)LC(g)),mod) || !rem(NM((Q)LC(h)),mod) ) + continue; + ptomp(mod,g,&tg); ptomp(mod,h,&th); + srchump(mod,tg,th,&tr); + if ( !tr ) + continue; + else + modres = CONT((MQ)tr); - mptoum(tg,wg); mptoum(th,wh); - eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */ - DEG(wm) = 0; COEF(wm)[0] = modres; - mulum(mod,ws,wm,wt); - for ( i = DEG(wt); i >= 0; i-- ) - if ( ( COEF(wt)[i] * 2 ) > mod ) - COEF(wt)[i] -= mod; - chnrem(mod,v,c,q,wt,&c1,&q1); - if ( !ucmpp(c,c1) ) { - mulp(CO,c,g,&z); divsrp(CO,z,h,&zz,&zzz); - if ( NUM(zzz) ) { - q = q1; c = c1; - break; - } - } - q = q1; c = c1; + mptoum(tg,wg); mptoum(th,wh); + eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */ + DEG(wm) = 0; COEF(wm)[0] = modres; + mulum(mod,ws,wm,wt); + for ( i = DEG(wt); i >= 0; i-- ) + if ( ( COEF(wt)[i] * 2 ) > mod ) + COEF(wt)[i] -= mod; + chnrem(mod,v,c,q,wt,&c1,&q1); + if ( !ucmpp(c,c1) ) { + mulp(CO,c,g,&z); divsrp(CO,z,h,&zz,&zzz); + if ( NUM(zzz) ) { + q = q1; c = c1; + break; + } + } + q = q1; c = c1; - if ( cmpq(f,q) < 0 ) - break; - } - ptozp(c,1,&s,&gc0); - /* gc0*g = zz*h + zzz -> c*g + (-zz)*h = zzz */ - mulp(CO,gc0,g,&z); - divsrp(CO,z,h,&zz,&zzz); - ptozp(zz,1,&s,(P *)&t); - if ( INT((Q)s) ) - chsgnp(zz,&hc0); - else { - 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); - } - if ( !INT((Q)zzz) ) { - 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; - } - for ( index = 0; ; ) { - mod = get_lprime(index++); - if ( !rem(NM((Q)zzz),mod) || - !rem(NM((Q)LC(g)),mod) || - !rem(NM((Q)LC(h)),mod) ) - continue; - for ( STOQ(mod,q); cmpq(q,bound) < 0; ) { - mulq(q,q,&q1); q = q1; - } - *qp = q; - invl((Q)zzz,q,&q1); - mulp(CO,gc0,(P)q1,hcp); mulp(CO,hc0,(P)q1,gcp); - return; - } + if ( cmpq(f,q) < 0 ) + break; + } + ptozp(c,1,&s,&gc0); + /* gc0*g = zz*h + zzz -> c*g + (-zz)*h = zzz */ + mulp(CO,gc0,g,&z); + divsrp(CO,z,h,&zz,&zzz); + ptozp(zz,1,&s,(P *)&t); + if ( INT((Q)s) ) + chsgnp(zz,&hc0); + else { + 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); + } + if ( !INT((Q)zzz) ) { + 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; + } + for ( index = 0; ; ) { + mod = get_lprime(index++); + if ( !rem(NM((Q)zzz),mod) || + !rem(NM((Q)LC(g)),mod) || + !rem(NM((Q)LC(h)),mod) ) + continue; + for ( STOQ(mod,q); cmpq(q,bound) < 0; ) { + mulq(q,q,&q1); q = q1; + } + *qp = q; + invl((Q)zzz,q,&q1); + mulp(CO,gc0,(P)q1,hcp); mulp(CO,hc0,(P)q1,gcp); + return; + } } void addm2p(VL vl,Q mod,Q mod2,P n1,P n2,P *nr) { - P t; + P t; - addp(vl,n1,n2,&t); - if ( !t ) - *nr = 0; - else - cm2p(mod,mod2,t,nr); + addp(vl,n1,n2,&t); + if ( !t ) + *nr = 0; + else + cm2p(mod,mod2,t,nr); } void subm2p(VL vl,Q mod,Q mod2,P n1,P n2,P *nr) { - P t; + P t; - subp(vl,n1,n2,&t); - if ( !t ) - *nr = 0; - else - cm2p(mod,mod2,t,nr); + subp(vl,n1,n2,&t); + if ( !t ) + *nr = 0; + else + cm2p(mod,mod2,t,nr); } void mulm2p(VL vl,Q mod,Q mod2,P n1,P n2,P *nr) { - P t; + P t; - mulp(vl,n1,n2,&t); - if ( !t ) - *nr = 0; - else - cm2p(mod,mod2,t,nr); + mulp(vl,n1,n2,&t); + if ( !t ) + *nr = 0; + else + cm2p(mod,mod2,t,nr); } void cmp(Q mod,P p,P *pr) { - P t; - DCP dc,dcr,dcr0; + P t; + DCP dc,dcr,dcr0; - if ( !p ) - *pr = 0; - else if ( NUM(p) ) - remq((Q)p,mod,(Q *)pr); - else { - for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) { - cmp(mod,COEF(dc),&t); - if ( t ) { - NEXTDC(dcr0,dcr); - DEG(dcr) = DEG(dc); - COEF(dcr) = t; - } - } - if ( !dcr0 ) - *pr = 0; - else { - NEXT(dcr) = 0; - MKP(VR(p),dcr0,*pr); - } - } + if ( !p ) + *pr = 0; + else if ( NUM(p) ) + remq((Q)p,mod,(Q *)pr); + else { + for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) { + cmp(mod,COEF(dc),&t); + if ( t ) { + NEXTDC(dcr0,dcr); + DEG(dcr) = DEG(dc); + COEF(dcr) = t; + } + } + if ( !dcr0 ) + *pr = 0; + else { + NEXT(dcr) = 0; + MKP(VR(p),dcr0,*pr); + } + } } void cm2p(Q mod,Q m,P p,P *pr) { - P t; - DCP dc,dcr,dcr0; + P t; + DCP dc,dcr,dcr0; - if ( !p ) - *pr = 0; - else if ( NUM(p) ) - rem2q((Q)p,mod,m,(Q *)pr); - else { - for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) { - cm2p(mod,m,COEF(dc),&t); - if ( t ) { - NEXTDC(dcr0,dcr); - DEG(dcr) = DEG(dc); - COEF(dcr) = t; - } - } - if ( !dcr0 ) - *pr = 0; - else { - NEXT(dcr) = 0; - MKP(VR(p),dcr0,*pr); - } - } + if ( !p ) + *pr = 0; + else if ( NUM(p) ) + rem2q((Q)p,mod,m,(Q *)pr); + else { + for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) { + cm2p(mod,m,COEF(dc),&t); + if ( t ) { + NEXTDC(dcr0,dcr); + DEG(dcr) = DEG(dc); + COEF(dcr) = t; + } + } + if ( !dcr0 ) + *pr = 0; + else { + NEXT(dcr) = 0; + MKP(VR(p),dcr0,*pr); + } + } } void addm2q(Q mod,Q mod2,Q n1,Q n2,Q *nr) { - Q t; + Q t; - addq(n1,n2,&t); - if ( !t ) - *nr = 0; - else - rem2q(t,mod,mod2,nr); + addq(n1,n2,&t); + if ( !t ) + *nr = 0; + else + rem2q(t,mod,mod2,nr); } void subm2q(Q mod,Q mod2,Q n1,Q n2,Q *nr) { - Q t; + Q t; - subq(n1,n2,&t); - if ( !t ) - *nr = 0; - else - rem2q(t,mod,mod2,nr); + subq(n1,n2,&t); + if ( !t ) + *nr = 0; + else + rem2q(t,mod,mod2,nr); } void mulm2q(Q mod,Q mod2,Q n1,Q n2,Q *nr) { - Q t; + Q t; - mulq(n1,n2,&t); - if ( !t ) - *nr = 0; - else - rem2q(t,mod,mod2,nr); + mulq(n1,n2,&t); + if ( !t ) + *nr = 0; + else + rem2q(t,mod,mod2,nr); } void rem2q(Q n,Q m,Q m2,Q *nr) { - N q,r,s; - int sgn; + N q,r,s; + int sgn; - divn(NM(n),NM(m),&q,&r); - if ( !r ) - *nr = 0; - else { - sgn = cmpn(r,NM(m2)); - if ( sgn > 0 ) { - subn(NM(m),r,&s); NTOQ(s,-1*SGN(n),*nr); - } else - NTOQ(r,SGN(n),*nr); - } + divn(NM(n),NM(m),&q,&r); + if ( !r ) + *nr = 0; + else { + sgn = cmpn(r,NM(m2)); + if ( sgn > 0 ) { + subn(NM(m),r,&s); NTOQ(s,-1*SGN(n),*nr); + } else + NTOQ(r,SGN(n),*nr); + } } /* - extract d-homogeneous part with respect to vl - {v} + extract d-homogeneous part with respect to vl - {v} */ void exthpc_generic(VL vl,P p,int d,V v,P *pr) { - P w,x,t,t1,a,xd; - V v0; - DCP dc; + P w,x,t,t1,a,xd; + V v0; + DCP dc; - if ( d < 0 || !p ) - *pr = 0; - else if ( NUM(p) ) - if ( d == 0 ) - *pr = p; - else - *pr = 0; - else if ( v == VR(p) ) - exthpc(vl,v,p,d,pr); - else { - v0 = VR(p); - for ( MKV(v0,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) { - exthpc_generic(vl,COEF(dc),d-QTOS(DEG(dc)),v,&t); - pwrp(vl,x,DEG(dc),&xd); - mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a; - } - *pr = w; - } + if ( d < 0 || !p ) + *pr = 0; + else if ( NUM(p) ) + if ( d == 0 ) + *pr = p; + else + *pr = 0; + else if ( v == VR(p) ) + exthpc(vl,v,p,d,pr); + else { + v0 = VR(p); + for ( MKV(v0,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) { + exthpc_generic(vl,COEF(dc),d-QTOS(DEG(dc)),v,&t); + pwrp(vl,x,DEG(dc),&xd); + mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a; + } + *pr = w; + } } void exthp(VL vl,P p,int d,P *pr) { - P t,t1,a,w,x,xd; - DCP dc; + P t,t1,a,w,x,xd; + DCP dc; - if ( d < 0 ) - *pr = 0; - else if ( NUM(p) ) - if ( d == 0 ) - *pr = p; - else - *pr = 0; - else { - for ( MKV(VR(p),x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) { - exthp(vl,COEF(dc),d - QTOS(DEG(dc)),&t); - pwrp(vl,x,DEG(dc),&xd); - mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a; - } - *pr = w; - } + if ( d < 0 ) + *pr = 0; + else if ( NUM(p) ) + if ( d == 0 ) + *pr = p; + else + *pr = 0; + else { + for ( MKV(VR(p),x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) { + exthp(vl,COEF(dc),d - QTOS(DEG(dc)),&t); + pwrp(vl,x,DEG(dc),&xd); + mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a; + } + *pr = w; + } } void exthpc(VL vl,V v,P p,int d,P *pr) { - P t,t1,a,w,x,xd; - DCP dc; + P t,t1,a,w,x,xd; + DCP dc; - if ( v != VR(p) ) - 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; - } + if ( v != VR(p) ) + 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; + } } void cbound(VL vl,P p,Q *b) { - Q n,e,t,m; - int k; + Q n,e,t,m; + int k; - cmax(p,&n); - addq(n,n,&m); + cmax(p,&n); + addq(n,n,&m); - k = geldb(vl,p); - STOQ(3,t); STOQ(k,e); + k = geldb(vl,p); + STOQ(3,t); STOQ(k,e); - pwrq(t,e,&n); - mulq(m,n,b); + pwrq(t,e,&n); + mulq(m,n,b); } int geldb(VL vl,P p) { - int m; + int m; - for ( m = 0; vl; vl = NEXT(vl) ) - m += getdeg(vl->v,p); - return ( m ); + for ( m = 0; vl; vl = NEXT(vl) ) + m += getdeg(vl->v,p); + return ( m ); } int getdeg(V v,P p) { - int m,t; - DCP dc; - - if ( !p || NUM(p) ) - return ( 0 ); - else if ( v == VR(p) ) - return ( deg(v,p) ); - else { - for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) { - t = getdeg(v,COEF(dc)); - m = MAX(m,t); - } - return ( m ); - } + int m,t; + DCP dc; + + if ( !p || NUM(p) ) + return ( 0 ); + else if ( v == VR(p) ) + return ( deg(v,p) ); + else { + for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) { + t = getdeg(v,COEF(dc)); + m = MAX(m,t); + } + return ( m ); + } } void cmax(P p,Q *b) { - DCP dc; - Q m,m1; - N tn; + DCP dc; + Q m,m1; + N tn; - if ( NUM(p) ) { - tn = NM((Q)p); - NTOQ(tn,1,*b); - } else { - for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) { - cmax(COEF(dc),&m1); - if ( cmpq(m1,m) > 0 ) - m = m1; - } - *b = m; - } + if ( NUM(p) ) { + tn = NM((Q)p); + NTOQ(tn,1,*b); + } else { + for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) { + cmax(COEF(dc),&m1); + if ( cmpq(m1,m) > 0 ) + m = m1; + } + *b = m; + } } int nextbin(VN vn,int n) { - int tmp,i,carry; + int tmp,i,carry; - if ( n == 0 ) - return ( 1 ); + if ( n == 0 ) + return ( 1 ); - for ( i = n - 1, carry = 1; i >= 0; i-- ) { - tmp = vn[i].n + carry; - vn[i].n = tmp % 2; - carry = tmp / 2; - } - return ( carry ); + for ( i = n - 1, carry = 1; i >= 0; i-- ) { + tmp = vn[i].n + carry; + vn[i].n = tmp % 2; + carry = tmp / 2; + } + return ( carry ); } void mulsgn(VN vn,VN vnt,int n,VN vn1) { - int i; + int i; - for ( i = 0; vn[i].v; i++ ) - vn1[i].n = vn[i].n; - for ( i = 0; i < n; i++ ) - if ( vnt[i].n ) - vn1[(int)vnt[i].v].n *= -1; + for ( i = 0; vn[i].v; i++ ) + vn1[i].n = vn[i].n; + for ( i = 0; i < n; i++ ) + if ( vnt[i].n ) + vn1[(int)vnt[i].v].n *= -1; } void next(VN vn) { - int i,m,n,tmp,carry; + int i,m,n,tmp,carry; - for ( m = 0, i = 0; vn[i].v; i++ ) - m = MAX(m,ABS(vn[i].n)); - if ( m == 0 ) { - vn[--i].n = 1; - return; - } - for ( carry = 1, m++, i--; ( i >= 0 ) && carry; i-- ) { - tmp = vn[i].n + carry; - vn[i].n = tmp % m; - carry = tmp / m; - } - if ( ( i == -1 ) && carry ) { - for ( i = 0; vn[i].v; i++ ) - vn[i].n = 0; - vn[--i].n = m; - } else { - for ( n = 0, i = 0; vn[i].v; i++ ) - n = MAX(n,ABS(vn[i].n)); - if ( n < m - 1 ) - vn[--i].n = m - 1; - } + for ( m = 0, i = 0; vn[i].v; i++ ) + m = MAX(m,ABS(vn[i].n)); + if ( m == 0 ) { + vn[--i].n = 1; + return; + } + for ( carry = 1, m++, i--; ( i >= 0 ) && carry; i-- ) { + tmp = vn[i].n + carry; + vn[i].n = tmp % m; + carry = tmp / m; + } + if ( ( i == -1 ) && carry ) { + for ( i = 0; vn[i].v; i++ ) + vn[i].n = 0; + vn[--i].n = m; + } else { + for ( n = 0, i = 0; vn[i].v; i++ ) + n = MAX(n,ABS(vn[i].n)); + if ( n < m - 1 ) + vn[--i].n = m - 1; + } } - + void clctv(VL vl,P p,VL *nvlp) { - int i,n; - VL tvl; - VN tvn; + int i,n; + VL tvl; + VN tvn; - if ( !p || NUM(p) ) { - *nvlp = 0; - return; - } + if ( !p || NUM(p) ) { + *nvlp = 0; + return; + } - for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ ); - tvn = (VN) ALLOCA((n+1)*sizeof(struct oVN)); - for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) { - tvn[i].v = tvl->v; - tvn[i].n = 0; - } + for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ ); + tvn = (VN) ALLOCA((n+1)*sizeof(struct oVN)); + for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) { + tvn[i].v = tvl->v; + tvn[i].n = 0; + } - markv(tvn,n,p); - vntovl(tvn,n,nvlp); + markv(tvn,n,p); + vntovl(tvn,n,nvlp); } void markv(VN vn,int n,P p) { - V v; - DCP dc; - int i; + V v; + DCP dc; + int i; - if ( NUM(p) ) - return; - v = VR(p); - for ( i = 0, v = VR(p); i < n; i++ ) - if ( v == vn[i].v ) { - vn[i].n = 1; - break; - } - for ( dc = DC(p); dc; dc = NEXT(dc) ) - markv(vn,n,COEF(dc)); + if ( NUM(p) ) + return; + v = VR(p); + for ( i = 0, v = VR(p); i < n; i++ ) + if ( v == vn[i].v ) { + vn[i].n = 1; + break; + } + for ( dc = DC(p); dc; dc = NEXT(dc) ) + markv(vn,n,COEF(dc)); } - + void vntovl(VN vn,int n,VL *vlp) { - int i; - VL tvl,tvl0; + int i; + VL tvl,tvl0; - for ( i = 0, tvl0 = 0; ; ) { - while ( ( i < n ) && ( vn[i].n == 0 ) ) i++; - if ( i == n ) - break; - else { - if ( !tvl0 ) { - NEWVL(tvl0); - tvl = tvl0; - } else { - NEWVL(NEXT(tvl)); - tvl = NEXT(tvl); - } - tvl->v = vn[i++].v; - } - } - if ( tvl0 ) - NEXT(tvl) = 0; - *vlp = tvl0; + for ( i = 0, tvl0 = 0; ; ) { + while ( ( i < n ) && ( vn[i].n == 0 ) ) i++; + if ( i == n ) + break; + else { + if ( !tvl0 ) { + NEWVL(tvl0); + tvl = tvl0; + } else { + NEWVL(NEXT(tvl)); + tvl = NEXT(tvl); + } + tvl->v = vn[i++].v; + } + } + if ( tvl0 ) + NEXT(tvl) = 0; + *vlp = tvl0; } int dbound(V v,P f) { - int m,t; - DCP dc; + int m,t; + DCP dc; - if ( !f ) - return ( -1 ); - else if ( v != VR(f) ) - return homdeg(f); - else { - for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) { - t = homdeg(COEF(dc)); - m = MAX(m,t); - } - return ( m ); - } + if ( !f ) + return ( -1 ); + else if ( v != VR(f) ) + return homdeg(f); + else { + for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) { + t = homdeg(COEF(dc)); + m = MAX(m,t); + } + return ( m ); + } } int homdeg(P f) { - int m,t; - DCP dc; + int m,t; + DCP dc; - if ( !f ) - return ( -1 ); - else if ( NUM(f) ) - return( 0 ); - else { - for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) { - t = QTOS(DEG(dc))+homdeg(COEF(dc)); - m = MAX(m,t); - } - return ( m ); - } + if ( !f ) + return ( -1 ); + else if ( NUM(f) ) + return( 0 ); + else { + for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) { + t = QTOS(DEG(dc))+homdeg(COEF(dc)); + m = MAX(m,t); + } + return ( m ); + } } int minhomdeg(P f) { - int m,t; - DCP dc; + int m,t; + DCP dc; - if ( !f ) - return ( -1 ); - else if ( NUM(f) ) - return( 0 ); - else { - for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) ) { - t = QTOS(DEG(dc))+minhomdeg(COEF(dc)); - m = MIN(m,t); - } - return ( m ); - } + if ( !f ) + return ( -1 ); + else if ( NUM(f) ) + return( 0 ); + else { + for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) ) { + t = QTOS(DEG(dc))+minhomdeg(COEF(dc)); + m = MIN(m,t); + } + return ( m ); + } } void adjc(VL vl,P f,P a,P lc0,Q q,P *fr,P *ar) { - P m,m1; - Q t; + P m,m1; + Q t; - invl((Q)LC(f),q,&t); - mulq((Q)lc0,t,(Q *)&m); - mulpq(f,m,&m1); cmp(q,m1,fr); - invl((Q)m,q,&t); - mulpq(a,(P)t,&m1); - cmp(q,m1,ar); + invl((Q)LC(f),q,&t); + mulq((Q)lc0,t,(Q *)&m); + mulpq(f,m,&m1); cmp(q,m1,fr); + invl((Q)m,q,&t); + mulpq(a,(P)t,&m1); + cmp(q,m1,ar); } #if 1 void affinemain(VL vl,P p,V v0,int n,P *pl,P *pr) { - P x,t,m,c,s,a; - DCP dc; - Q d; + P x,t,m,c,s,a; + DCP dc; + Q d; - if ( !p ) - *pr = 0; - else if ( NUM(p) ) - *pr = p; - else if ( VR(p) != v0 ) { - MKV(VR(p),x); - for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) { - affinemain(vl,COEF(dc),v0,n,pl,&t); - if ( DEG(dc) ) { - pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m); - addp(vl,m,c,&a); c = a; - } else { - addp(vl,t,c,&a); c = a; - } - } - *pr = c; - } else { - dc = DC(p); - c = COEF(dc); - for ( d = DEG(dc), dc = NEXT(dc); - dc; d = DEG(dc), dc = NEXT(dc) ) { - mulp(vl,pl[QTOS(d)-QTOS(DEG(dc))],c,&m); - addp(vl,m,COEF(dc),&c); - } - if ( d ) { - mulp(vl,pl[QTOS(d)],c,&m); c = m; - } - *pr = c; - } + if ( !p ) + *pr = 0; + else if ( NUM(p) ) + *pr = p; + else if ( VR(p) != v0 ) { + MKV(VR(p),x); + for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) { + affinemain(vl,COEF(dc),v0,n,pl,&t); + if ( DEG(dc) ) { + pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m); + addp(vl,m,c,&a); c = a; + } else { + addp(vl,t,c,&a); c = a; + } + } + *pr = c; + } else { + dc = DC(p); + c = COEF(dc); + for ( d = DEG(dc), dc = NEXT(dc); + dc; d = DEG(dc), dc = NEXT(dc) ) { + mulp(vl,pl[QTOS(d)-QTOS(DEG(dc))],c,&m); + addp(vl,m,COEF(dc),&c); + } + if ( d ) { + mulp(vl,pl[QTOS(d)],c,&m); c = m; + } + *pr = c; + } } #endif #if 0 affine(VL vl,P p,VN vn,P *r) { - int n,d,d1,i; - Q *q; - Q **bc; + int n,d,d1,i; + Q *q; + Q **bc; - if ( !p || NUM(p) ) - *r = p; - else { - for ( i = 0, d = 0; vn[i].v; i++ ) - d1 = getdeg(vn[i].v,p), d = MAX(d,d1); - W_CALLOC(d+1,Q *,bc); - for ( i = 0; i <= d; i++ ) - W_CALLOC(i+1,Q,q), mkbc(i,q), bc[i] = q; - afmain(vl,bc,p,vn,r); - } + if ( !p || NUM(p) ) + *r = p; + else { + for ( i = 0, d = 0; vn[i].v; i++ ) + d1 = getdeg(vn[i].v,p), d = MAX(d,d1); + W_CALLOC(d+1,Q *,bc); + for ( i = 0; i <= d; i++ ) + W_CALLOC(i+1,Q,q), mkbc(i,q), bc[i] = q; + afmain(vl,bc,p,vn,r); + } } afmain(VL vl,Q **bc,P p,VN vn,P *r) { - P t,s,u; - P *c,*rc; - Q *q; - DCP dc; - int n,i,j; - - if ( !p || NUM(p) || !vn[0].v ) - *r = p; - else if ( vn[0].v != VR(p) ) { - for ( i = 0; vn[i].v && (vn[i].v != VR(p)); i++ ); - if ( vn[i].v ) - afmain(vl,bc,p,vn+i,r); - else { - n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c); - for ( dc = DC(p); dc; dc = NEXT(dc) ) - afmain(vl,bc,COEF(dc),vn,&c[QTOS(DEG(dc))]); - plisttop(c,VR(p),n,r); - } - } else { - n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c); - W_CALLOC(n+1,P,rc); W_CALLOC(n+1,Q,q); - for ( dc = DC(p); dc; dc = NEXT(dc) ) - afmain(vl,bc,COEF(dc),vn+1,&c[QTOS(DEG(dc))]); - if ( !vn[0].n ) - bcopy(c,rc,sizeof(P)*(n+1)); - else { - for ( q[0] = ONE, STOQ(vn[0].n,q[1]), i = 2; i <= n; i++ ) - mulq(q[i-1],q[1],&q[i]); - for ( j = 0; j <= n; rc[j] = t, j++ ) - for ( i = j, t = 0; i <= n; i++ ) - if ( c[i] ) - mulq(bc[i][j],q[i-j],&s), mulp(vl,s,c[i],&u), - addp(CO,u,t,&s), t = s; - } - plisttop(rc,VR(p),n,r); - } + P t,s,u; + P *c,*rc; + Q *q; + DCP dc; + int n,i,j; + + if ( !p || NUM(p) || !vn[0].v ) + *r = p; + else if ( vn[0].v != VR(p) ) { + for ( i = 0; vn[i].v && (vn[i].v != VR(p)); i++ ); + if ( vn[i].v ) + afmain(vl,bc,p,vn+i,r); + else { + n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c); + for ( dc = DC(p); dc; dc = NEXT(dc) ) + afmain(vl,bc,COEF(dc),vn,&c[QTOS(DEG(dc))]); + plisttop(c,VR(p),n,r); + } + } else { + n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c); + W_CALLOC(n+1,P,rc); W_CALLOC(n+1,Q,q); + for ( dc = DC(p); dc; dc = NEXT(dc) ) + afmain(vl,bc,COEF(dc),vn+1,&c[QTOS(DEG(dc))]); + if ( !vn[0].n ) + bcopy(c,rc,sizeof(P)*(n+1)); + else { + for ( q[0] = ONE, STOQ(vn[0].n,q[1]), i = 2; i <= n; i++ ) + mulq(q[i-1],q[1],&q[i]); + for ( j = 0; j <= n; rc[j] = t, j++ ) + for ( i = j, t = 0; i <= n; i++ ) + if ( c[i] ) + mulq(bc[i][j],q[i-j],&s), mulp(vl,s,c[i],&u), + addp(CO,u,t,&s), t = s; + } + plisttop(rc,VR(p),n,r); + } } #endif void restore(VL vl,P f,VN vn,P *fr) { - int i; - P vv,g,g1,t; - Q s; + int i; + P vv,g,g1,t; + Q s; - g = f; - for ( i = 0; vn[i].v; i++ ) { - MKV(vn[i].v,t); - if ( vn[i].n ) { - STOQ(-vn[i].n,s); - addp(vl,t,(P)s,&vv); - } else - vv = t; + g = f; + for ( i = 0; vn[i].v; i++ ) { + MKV(vn[i].v,t); + if ( vn[i].n ) { + STOQ(-vn[i].n,s); + addp(vl,t,(P)s,&vv); + } else + vv = t; - substp(vl,g,vn[i].v,vv,&g1); g = g1; - } - *fr = g; + substp(vl,g,vn[i].v,vv,&g1); g = g1; + } + *fr = g; } void mergev(VL vl,VL vl1,VL vl2,VL *nvlp) { - int i,n; - VL tvl; - VN vn; + int i,n; + VL tvl; + VN vn; - if ( !vl1 ) { - *nvlp = vl2; return; - } else if ( !vl2 ) { - *nvlp = vl1; return; - } - for ( i = 0, tvl = vl; tvl; i++, tvl = NEXT(tvl) ); - n = i; - W_CALLOC(n,struct oVN,vn); - for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) - vn[i].v = tvl->v; - for ( i = 0, tvl = vl1; (i < n) && tvl; i++, tvl = NEXT(tvl)) { - while ( ( i < n ) && ( vn[i].v != tvl->v ) ) - i++; - if ( i == n ) - break; - else - vn[i].n = 1; - } - for ( i = 0, tvl = vl2; (i < n) && tvl; i++, tvl = NEXT(tvl) ) { - while ( ( i < n ) && ( vn[i].v != tvl->v ) ) - i++; - if ( i == n ) - break; - else - vn[i].n = 1; - } - vntovl(vn,n,nvlp); + if ( !vl1 ) { + *nvlp = vl2; return; + } else if ( !vl2 ) { + *nvlp = vl1; return; + } + for ( i = 0, tvl = vl; tvl; i++, tvl = NEXT(tvl) ); + n = i; + W_CALLOC(n,struct oVN,vn); + for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) + vn[i].v = tvl->v; + for ( i = 0, tvl = vl1; (i < n) && tvl; i++, tvl = NEXT(tvl)) { + while ( ( i < n ) && ( vn[i].v != tvl->v ) ) + i++; + if ( i == n ) + break; + else + vn[i].n = 1; + } + for ( i = 0, tvl = vl2; (i < n) && tvl; i++, tvl = NEXT(tvl) ) { + while ( ( i < n ) && ( vn[i].v != tvl->v ) ) + i++; + if ( i == n ) + break; + else + vn[i].n = 1; + } + vntovl(vn,n,nvlp); } #if 0 void substvp(VL vl,P f,VN vn,P *g) { - V v; - int i; - P h,h1; - Q t; + V v; + int i; + P h,h1; + Q t; - h = f; - for ( i = 0; v = vn[i].v; i++ ) { - STOQ(vn[i].n,t); - substp(vl,h,v,(P)t,&h1); h = h1; - } - *g = h; + h = f; + for ( i = 0; v = vn[i].v; i++ ) { + STOQ(vn[i].n,t); + substp(vl,h,v,(P)t,&h1); h = h1; + } + *g = h; } void affine(VL vl,P f,VN vn,P *fr) { - int i,j,n; - P vv,g,g1,t,u; - Q s; - int *dlist; - P **plist; + int i,j,n; + P vv,g,g1,t,u; + Q s; + int *dlist; + P **plist; - for ( n = 0; vn[n].v; n++); - dlist = (int *)ALLOCA((n+1)*sizeof(int)); - plist = (P **)ALLOCA((n+1)*sizeof(P *)); - for ( i = 0; vn[i].v; i++ ) { - if ( !vn[i].n ) - continue; - dlist[i] = getdeg(vn[i].v,f); - plist[i] = (P *)ALLOCA((dlist[i]+1)*sizeof(P)); + for ( n = 0; vn[n].v; n++); + dlist = (int *)ALLOCA((n+1)*sizeof(int)); + plist = (P **)ALLOCA((n+1)*sizeof(P *)); + for ( i = 0; vn[i].v; i++ ) { + if ( !vn[i].n ) + continue; + dlist[i] = getdeg(vn[i].v,f); + plist[i] = (P *)ALLOCA((dlist[i]+1)*sizeof(P)); - MKV(vn[i].v,t); - if ( vn[i].n ) { - STOQ(vn[i].n,s); addp(vl,t,(P)s,&vv); - } else - vv = t; + MKV(vn[i].v,t); + if ( vn[i].n ) { + STOQ(vn[i].n,s); addp(vl,t,(P)s,&vv); + } else + vv = t; - for ( j = 0, t = (P)ONE; j < dlist[i]; j++ ) { - plist[i][j] = t; - mulp(vl,t,vv,&u); - t = u; - } - plist[i][j] = t; - } + for ( j = 0, t = (P)ONE; j < dlist[i]; j++ ) { + plist[i][j] = t; + mulp(vl,t,vv,&u); + t = u; + } + plist[i][j] = t; + } - g = f; - for ( i = 0; vn[i].v; i++ ) { - if ( !vn[i].n ) - continue; - affinemain(vl,g,vn[i].v,dlist[i],plist[i],&g1); g = g1; - } - *fr = g; + g = f; + for ( i = 0; vn[i].v; i++ ) { + if ( !vn[i].n ) + continue; + affinemain(vl,g,vn[i].v,dlist[i],plist[i],&g1); g = g1; + } + *fr = g; } #endif