![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to E.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to E.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/E.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #include "ca.h"
! 3:
! 4: void henmv(vl,vn,f,g0,h0,a0,b0,lg,lh,lg0,lh0,q,k,gp,hp)
! 5: VL vl;
! 6: VN vn;
! 7: P f,g0,h0,a0,b0,lg,lh,lg0,lh0;
! 8: Q q;
! 9: int k;
! 10: P *gp,*hp;
! 11: {
! 12: P g1,h1,a1,b1;
! 13: N qn;
! 14: Q q2;
! 15:
! 16: divin((NM(q)),2,&qn); NTOQ(qn,1,q2);
! 17: adjc(vl,g0,a0,lg0,q,&g1,&a1); adjc(vl,h0,b0,lh0,q,&h1,&b1);
! 18: henmvmain(vl,vn,f,g1,h1,b1,a1,lg,lh,q,q2,k,gp,hp);
! 19: }
! 20:
! 21: void henmvmain(vl,vn,f,fi0,fi1,gi0,gi1,l0,l1,mod,mod2,k,fr0,fr1)
! 22: VL vl;
! 23: VN vn;
! 24: P f,fi0,fi1,gi0,gi1,l0,l1;
! 25: Q mod,mod2;
! 26: int k;
! 27: P *fr0,*fr1;
! 28: {
! 29: V v;
! 30: int n,i,j;
! 31: int *md;
! 32: P x,m,m1,c,q,r,a,s,u,ff,f0,f1;
! 33: P w0,w1,cf,cfi,t,q1,dvr;
! 34: P *c0,*c1;
! 35: P *f0h,*f1h;
! 36:
! 37: v = VR(f); n = deg(v,f); MKV(v,x);
! 38: c0 = (P *)ALLOCA((n+1)*sizeof(P));
! 39: c1 = (P *)ALLOCA((n+1)*sizeof(P));
! 40: invl((Q)LC(fi1),mod,(Q *)&c); mulp(vl,fi1,c,&t); cmp(mod,t,&dvr);
! 41: cm2p(mod,mod2,gi0,&c0[0]); cm2p(mod,mod2,gi1,&c1[0]);
! 42: for ( i = 1; i <= n; i++ ) {
! 43: mulp(vl,x,c1[i-1],&m); divsrp(vl,m,dvr,&q,&r); mulp(vl,q,c,&q1);
! 44: cm2p(mod,mod2,r,&c1[i]);
! 45: mulp(vl,x,c0[i-1],&m); mulp(vl,q1,fi0,&m1); addp(vl,m,m1,&a);
! 46: cm2p(mod,mod2,a,&c0[i]);
! 47: }
! 48: affine(vl,f,vn,&t); cm2p(mod,mod2,t,&ff);
! 49: for ( i = 0; vn[i].v; i++ );
! 50: md = ( int *) ALLOCA((i+1)*sizeof(int));
! 51: for ( i = 0; vn[i].v; i++ )
! 52: md[i] = getdeg(vn[i].v,ff);
! 53: cm2p(mod,mod2,fi0,&f0); affine(vl,l0,vn,&t);
! 54: if ( NUM(f0) )
! 55: cm2p(mod,mod2,t,&f0);
! 56: else
! 57: cm2p(mod,mod2,t,&COEF(DC(f0)));
! 58: cm2p(mod,mod2,fi1,&f1); affine(vl,l1,vn,&t);
! 59: if ( NUM(f1) )
! 60: cm2p(mod,mod2,t,&f1);
! 61: else
! 62: cm2p(mod,mod2,t,&COEF(DC(f1)));
! 63: W_CALLOC(k,P,f0h); W_CALLOC(k,P,f1h);
! 64: for ( i = 0; i <= k; i++ ) {
! 65: exthpc(vl,v,f0,i,&f0h[i]); exthpc(vl,v,f1,i,&f1h[i]);
! 66: }
! 67: for ( j = 1; j <= k; j++ ) {
! 68: for ( i = 0; vn[i].v; i++ )
! 69: if ( getdeg(vn[i].v,f0)+getdeg(vn[i].v,f1) > md[i] )
! 70: goto END;
! 71: for ( i = 0, s = 0; i <= j; i++ ) {
! 72: mulp(vl,f0h[i],f1h[j-i],&t); addp(vl,s,t,&u); s = u;
! 73: }
! 74: cm2p(mod,mod2,s,&t);
! 75: exthpc(vl,v,ff,j,&u); subm2p(vl,mod,mod2,u,t,&cf);
! 76: for ( i = 0, w0 = 0, w1 = 0; i <= n; i++ ) {
! 77: if ( !cf )
! 78: cfi = 0;
! 79: else if ( VR(cf) == v )
! 80: coefp(cf,i,&cfi);
! 81: else if ( i == 0 )
! 82: cfi = cf;
! 83: else
! 84: cfi = 0;
! 85: if ( cfi ) {
! 86: mulp(vl,cfi,c0[i],&m); addp(vl,w0,m,&a); w0 = a;
! 87: mulp(vl,cfi,c1[i],&m); addp(vl,w1,m,&a); w1 = a;
! 88: }
! 89: }
! 90: cm2p(mod,mod2,w0,&t); addm2p(vl,mod,mod2,f0,t,&a);
! 91: addm2p(vl,mod,mod2,f0h[j],t,&s); f0h[j] = s; f0 = a;
! 92: cm2p(mod,mod2,w1,&u); addm2p(vl,mod,mod2,f1,u,&a);
! 93: addm2p(vl,mod,mod2,f1h[j],u,&s); f1h[j] = s; f1 = a;
! 94: if ( !t ) {
! 95: restore(vl,f0,vn,&s); cm2p(mod,mod2,s,&t);
! 96: if ( divtpz(vl,f,t,&s) ) {
! 97: *fr0 = t; *fr1 = s;
! 98: return;
! 99: }
! 100: }
! 101: if ( !u ) {
! 102: restore(vl,f1,vn,&s); cm2p(mod,mod2,s,&t);
! 103: if ( divtpz(vl,f,t,&s) ) {
! 104: *fr0 = s; *fr1 = t;
! 105: return;
! 106: }
! 107: }
! 108: }
! 109: END:
! 110: restore(vl,f0,vn,&t); cm2p(mod,mod2,t,fr0);
! 111: restore(vl,f1,vn,&t); cm2p(mod,mod2,t,fr1);
! 112: }
! 113:
! 114: /*
! 115: input : f, fi0, fi1, p, k; f = fi0 * fi1 mod p; ( p, k : integer )
! 116: output : gr1 * fr0 + gr0 * fr1 = 1 mod qr; ( qr = p^(2^(k-1)) )
! 117: */
! 118:
! 119: void henzq(f,i0,fi0,i1,fi1,p,k,fr0p,fr1p,gr0p,gr1p,qrp)
! 120: P f;
! 121: UM fi0,fi1;
! 122: int p,k;
! 123: P i0,i1;
! 124: P *fr0p,*fr1p,*gr0p,*gr1p;
! 125: Q *qrp;
! 126: {
! 127: N qn;
! 128: Q q,qq,q2;
! 129: int n,i;
! 130: UM wg0,wg1,wf0,wf1;
! 131: P f0,f1,g0,g1,m,m1,s,c,d,q1,r1,rm,rm1,a,a1,a2;
! 132:
! 133: n = UDEG(f);
! 134: wf0 = W_UMALLOC(n); wf1 = W_UMALLOC(n);
! 135: wg0 = W_UMALLOC(n); wg1 = W_UMALLOC(n);
! 136: cpyum(fi0,wf0); cpyum(fi1,wf1);
! 137: eucum(p,wf0,wf1,wg1,wg0);
! 138: umtop(VR(f),fi0,&f0); umtop(VR(f),fi1,&f1);
! 139: umtop(VR(f),wg0,&g0); umtop(VR(f),wg1,&g1);
! 140:
! 141: STOQ(p,q); divin(NM(q),2,&qn); NTOQ(qn,1,q2);
! 142: for ( i = 1; i < k; i++ ) {
! 143: #if 0
! 144: mulp(CO,i1,g0,&m); mulp(CO,i0,g1,&m1); addp(CO,m,m1,&a);
! 145: if ( NUM(a) ) {
! 146: for ( STOQ(p,q), j = 1; j < k; j++ ) {
! 147: mulq(q,q,&qq); q = qq;
! 148: }
! 149: f0 = i0; f1 = i1;
! 150: invl(a,q,&qq);
! 151: mulp(CO,g0,qq,&s); g0 = s; mulp(CO,g1,qq,&s); g1 = s;
! 152: break;
! 153: }
! 154: #endif
! 155: /* c = ((f - f0*f1)/q) mod q;
! 156: q1 = (c*g1) / f1;
! 157: r1 = (c*g1) mod f1;
! 158: f1 += (r1 mod q)*q;
! 159: f0 += ((c*g0 + q1*f0) mod q)*q;
! 160:
! 161: d = ((1 - (f1*g0 + f0*g1))/q) mod q;
! 162: q1 = (d*g0) / f1;
! 163: r1 = (d*g0) mod f1;
! 164: g1 += (r1 mod q)*q;
! 165: g0 += ((c*g0 + q1*f0) mod q)*q;
! 166: q = q^2;
! 167: */
! 168:
! 169: /* c = ((f - f0*f1)/q) mod q */
! 170: mulp(CO,f0,f1,&m); subp(CO,f,m,&s);
! 171: divcp(s,q,&m); cm2p(q,q2,m,&c);
! 172:
! 173: /* q1 = (c*g1) / f1; r1 = (c*g1) mod f1; */
! 174: mulp(CO,c,g1,&m); cm2p(q,q2,m,&s);
! 175: udivpwm(q,s,f1,&q1,&r1);
! 176:
! 177: /* f1 = f1 + (r1 mod q)*q; */
! 178: cm2p(q,q2,r1,&rm);
! 179: mulpq(rm,(P)q,&rm1); addp(CO,f1,rm1,&a);
! 180: f1 = a;
! 181:
! 182: /* a1 = (c*g0 + q1*f0) mod q; */
! 183: mulp(CO,c,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);
! 184: cm2p(q,q2,a,&a1);
! 185:
! 186: /* f0 = f0 + a1*q; */
! 187: mulpq(a1,(P)q,&a2);
! 188: addp(CO,f0,a2,&a);
! 189: f0 = a;
! 190:
! 191: /* d = ((1 - (f1*g0 + f0*g1))/q) mod q; */
! 192: mulp(CO,f1,g0,&m); mulp(CO,f0,g1,&m1); addp(CO,m,m1,&a);
! 193: subp(CO,(P)ONE,a,&s);
! 194: divcp(s,q,&m); cm2p(q,q2,m,&d);
! 195:
! 196: /* q1 = (d*g0) / f1; r1 = (d*g0) mod f1; */
! 197: mulp(CO,d,g1,&m); udivpwm(q,m,f1,&q1,&r1);
! 198:
! 199: /* g1 = g1 + (r1 mod q )*q; */
! 200: cm2p(q,q2,r1,&rm);
! 201: mulpq(rm,(P)q,&rm1); addp(CO,g1,rm1,&a);
! 202: g1 = a;
! 203:
! 204: /* g0 = g0 + (c*g0 + q1*f0 mod q)*q; */
! 205: mulp(CO,d,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);
! 206: cm2p(q,q2,a,&a1); mulpq(a1,(P)q,&a2);
! 207: addp(CO,g0,a2,&a);
! 208: g0 = a;
! 209:
! 210: /* q = q^2; */
! 211: mulq(q,q,&qq);
! 212: q = qq;
! 213: divin(NM(q),2,&qn); NTOQ(qn,1,q2);
! 214: }
! 215: *fr0p = f0; *fr1p = f1; *gr0p = g0; *gr1p = g1; *qrp = q;
! 216: }
! 217:
! 218: void henzq1(g,h,bound,gcp,hcp,qp)
! 219: P g,h;
! 220: Q bound;
! 221: P *gcp,*hcp;
! 222: Q *qp;
! 223: {
! 224: V v;
! 225: Q f,q,q1;
! 226: Q u,t,a,b,s;
! 227: P c,c1;
! 228: P tg,th,tr;
! 229: UM wg,wh,ws,wt,wm;
! 230: int n,m,modres,mod,index,i;
! 231: P gc0,hc0;
! 232: P z,zz,zzz;
! 233:
! 234:
! 235: v = VR(g); n=deg(v,g); m=deg(v,h);
! 236: norm(g,&a); norm(h,&b);
! 237: STOQ(m,u); pwrq(a,u,&t);
! 238: STOQ(n,u); pwrq(b,u,&s);
! 239: mulq(t,s,&u);
! 240:
! 241: factorial(n+m,&t); mulq(u,t,&s);
! 242: addq(s,s,&f);
! 243:
! 244: wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);
! 245: wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);
! 246: wm = W_UMALLOC(m+n);
! 247:
! 248: for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
! 249: mod = lprime[index++];
! 250: if ( !mod )
! 251: error("henzq1 : lprime[] exhausted.");
! 252: if ( !rem(NM((Q)LC(g)),mod) || !rem(NM((Q)LC(h)),mod) )
! 253: continue;
! 254: ptomp(mod,g,&tg); ptomp(mod,h,&th);
! 255: srchump(mod,tg,th,&tr);
! 256: if ( !tr )
! 257: continue;
! 258: else
! 259: modres = CONT((MQ)tr);
! 260:
! 261: mptoum(tg,wg); mptoum(th,wh);
! 262: eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
! 263: DEG(wm) = 0; COEF(wm)[0] = modres;
! 264: mulum(mod,ws,wm,wt);
! 265: for ( i = DEG(wt); i >= 0; i-- )
! 266: if ( ( COEF(wt)[i] * 2 ) > mod )
! 267: COEF(wt)[i] -= mod;
! 268: chnrem(mod,v,c,q,wt,&c1,&q1);
! 269: if ( !ucmpp(c,c1) ) {
! 270: mulp(CO,c,g,&z); divsrp(CO,z,h,&zz,&zzz);
! 271: if ( NUM(zzz) ) {
! 272: q = q1; c = c1;
! 273: break;
! 274: }
! 275: }
! 276: q = q1; c = c1;
! 277:
! 278: if ( cmpq(f,q) < 0 )
! 279: break;
! 280: }
! 281: ptozp(c,1,&s,&gc0);
! 282: /* gc0*g = zz*h + zzz -> c*g + (-zz)*h = zzz */
! 283: mulp(CO,gc0,g,&z);
! 284: divsrp(CO,z,h,&zz,&zzz);
! 285: ptozp(zz,1,&s,(P *)&t);
! 286: if ( INT((Q)s) )
! 287: chsgnp(zz,&hc0);
! 288: else {
! 289: NTOQ(DN((Q)s),1,q); mulq(q,(Q)zzz,&q1); zzz = (P)q1;
! 290: mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,zz,(P)q,&c); chsgnp(c,&hc0);
! 291: }
! 292: if ( !INT((Q)zzz) ) {
! 293: NTOQ(DN((Q)zzz),1,q); NTOQ(NM((Q)zzz),1,q1); zzz = (P)q1;
! 294: mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,hc0,(P)q,&c); hc0 = c;
! 295: }
! 296: for ( index = 0; ; ) {
! 297: mod = lprime[index++];
! 298: if ( !mod )
! 299: error("henzq1 : lprime[] exhausted.");
! 300: if ( !rem(NM((Q)zzz),mod) ||
! 301: !rem(NM((Q)LC(g)),mod) ||
! 302: !rem(NM((Q)LC(h)),mod) )
! 303: continue;
! 304: for ( STOQ(mod,q); cmpq(q,bound) < 0; ) {
! 305: mulq(q,q,&q1); q = q1;
! 306: }
! 307: *qp = q;
! 308: invl((Q)zzz,q,&q1);
! 309: mulp(CO,gc0,(P)q1,hcp); mulp(CO,hc0,(P)q1,gcp);
! 310: return;
! 311: }
! 312: }
! 313:
! 314: void addm2p(vl,mod,mod2,n1,n2,nr)
! 315: VL vl;
! 316: Q mod,mod2;
! 317: P n1,n2,*nr;
! 318: {
! 319: P t;
! 320:
! 321: addp(vl,n1,n2,&t);
! 322: if ( !t )
! 323: *nr = 0;
! 324: else
! 325: cm2p(mod,mod2,t,nr);
! 326: }
! 327:
! 328: void subm2p(vl,mod,mod2,n1,n2,nr)
! 329: VL vl;
! 330: Q mod,mod2;
! 331: P n1,n2,*nr;
! 332: {
! 333: P t;
! 334:
! 335: subp(vl,n1,n2,&t);
! 336: if ( !t )
! 337: *nr = 0;
! 338: else
! 339: cm2p(mod,mod2,t,nr);
! 340: }
! 341:
! 342: void mulm2p(vl,mod,mod2,n1,n2,nr)
! 343: VL vl;
! 344: Q mod,mod2;
! 345: P n1,n2,*nr;
! 346: {
! 347: P t;
! 348:
! 349: mulp(vl,n1,n2,&t);
! 350: if ( !t )
! 351: *nr = 0;
! 352: else
! 353: cm2p(mod,mod2,t,nr);
! 354: }
! 355:
! 356: void cmp(mod,p,pr)
! 357: Q mod;
! 358: P p,*pr;
! 359: {
! 360: P t;
! 361: DCP dc,dcr,dcr0;
! 362:
! 363: if ( !p )
! 364: *pr = 0;
! 365: else if ( NUM(p) )
! 366: remq((Q)p,mod,(Q *)pr);
! 367: else {
! 368: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 369: cmp(mod,COEF(dc),&t);
! 370: if ( t ) {
! 371: NEXTDC(dcr0,dcr);
! 372: DEG(dcr) = DEG(dc);
! 373: COEF(dcr) = t;
! 374: }
! 375: }
! 376: if ( !dcr0 )
! 377: *pr = 0;
! 378: else {
! 379: NEXT(dcr) = 0;
! 380: MKP(VR(p),dcr0,*pr);
! 381: }
! 382: }
! 383: }
! 384:
! 385: void cm2p(mod,m,p,pr)
! 386: Q mod,m;
! 387: P p,*pr;
! 388: {
! 389: P t;
! 390: DCP dc,dcr,dcr0;
! 391:
! 392: if ( !p )
! 393: *pr = 0;
! 394: else if ( NUM(p) )
! 395: rem2q((Q)p,mod,m,(Q *)pr);
! 396: else {
! 397: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 398: cm2p(mod,m,COEF(dc),&t);
! 399: if ( t ) {
! 400: NEXTDC(dcr0,dcr);
! 401: DEG(dcr) = DEG(dc);
! 402: COEF(dcr) = t;
! 403: }
! 404: }
! 405: if ( !dcr0 )
! 406: *pr = 0;
! 407: else {
! 408: NEXT(dcr) = 0;
! 409: MKP(VR(p),dcr0,*pr);
! 410: }
! 411: }
! 412: }
! 413:
! 414: void addm2q(mod,mod2,n1,n2,nr)
! 415: Q mod,mod2;
! 416: Q n1,n2,*nr;
! 417: {
! 418: Q t;
! 419:
! 420: addq(n1,n2,&t);
! 421: if ( !t )
! 422: *nr = 0;
! 423: else
! 424: rem2q(t,mod,mod2,nr);
! 425: }
! 426:
! 427: void subm2q(mod,mod2,n1,n2,nr)
! 428: Q mod,mod2;
! 429: Q n1,n2,*nr;
! 430: {
! 431: Q t;
! 432:
! 433: subq(n1,n2,&t);
! 434: if ( !t )
! 435: *nr = 0;
! 436: else
! 437: rem2q(t,mod,mod2,nr);
! 438: }
! 439:
! 440: void mulm2q(mod,mod2,n1,n2,nr)
! 441: Q mod,mod2;
! 442: Q n1,n2,*nr;
! 443: {
! 444: Q t;
! 445:
! 446: mulq(n1,n2,&t);
! 447: if ( !t )
! 448: *nr = 0;
! 449: else
! 450: rem2q(t,mod,mod2,nr);
! 451: }
! 452:
! 453: void rem2q(n,m,m2,nr)
! 454: Q n,m,m2,*nr;
! 455: {
! 456: N q,r,s;
! 457: int sgn;
! 458:
! 459: divn(NM(n),NM(m),&q,&r);
! 460: if ( !r )
! 461: *nr = 0;
! 462: else {
! 463: sgn = cmpn(r,NM(m2));
! 464: if ( sgn > 0 ) {
! 465: subn(NM(m),r,&s); NTOQ(s,-1*SGN(n),*nr);
! 466: } else
! 467: NTOQ(r,SGN(n),*nr);
! 468: }
! 469: }
! 470:
! 471: void exthp(vl,p,d,pr)
! 472: VL vl;
! 473: P p;
! 474: int d;
! 475: P *pr;
! 476: {
! 477: P t,t1,a,w,x,xd;
! 478: DCP dc;
! 479:
! 480: if ( d < 0 )
! 481: *pr = 0;
! 482: else if ( NUM(p) )
! 483: if ( d == 0 )
! 484: *pr = p;
! 485: else
! 486: *pr = 0;
! 487: else {
! 488: for ( MKV(VR(p),x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
! 489: exthp(vl,COEF(dc),d - QTOS(DEG(dc)),&t);
! 490: pwrp(vl,x,DEG(dc),&xd);
! 491: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
! 492: }
! 493: *pr = w;
! 494: }
! 495: }
! 496:
! 497: void exthpc(vl,v,p,d,pr)
! 498: VL vl;
! 499: V v;
! 500: P p;
! 501: int d;
! 502: P *pr;
! 503: {
! 504: P t,t1,a,w,x,xd;
! 505: DCP dc;
! 506:
! 507: if ( v != VR(p) )
! 508: exthp(vl,p,d,pr);
! 509: else if ( d < 0 )
! 510: *pr = 0;
! 511: else {
! 512: for ( MKV(v,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
! 513: exthp(vl,COEF(dc),d,&t);
! 514: pwrp(vl,x,DEG(dc),&xd);
! 515: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
! 516: }
! 517: *pr = w;
! 518: }
! 519: }
! 520:
! 521: void cbound(vl,p,b)
! 522: VL vl;
! 523: P p;
! 524: Q *b;
! 525: {
! 526: Q n,e,t,m;
! 527: int k;
! 528:
! 529: cmax(p,&n);
! 530: addq(n,n,&m);
! 531:
! 532: k = geldb(vl,p);
! 533: STOQ(3,t); STOQ(k,e);
! 534:
! 535: pwrq(t,e,&n);
! 536: mulq(m,n,b);
! 537: }
! 538:
! 539: int geldb(vl,p)
! 540: VL vl;
! 541: P p;
! 542: {
! 543: int m;
! 544:
! 545: for ( m = 0; vl; vl = NEXT(vl) )
! 546: m += getdeg(vl->v,p);
! 547: return ( m );
! 548: }
! 549:
! 550: int getdeg(v,p)
! 551: V v;
! 552: P p;
! 553: {
! 554: int m;
! 555: DCP dc;
! 556:
! 557: if ( !p || NUM(p) )
! 558: return ( 0 );
! 559: else if ( v == VR(p) )
! 560: return ( deg(v,p) );
! 561: else {
! 562: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) )
! 563: m = MAX(m,getdeg(v,COEF(dc)));
! 564: return ( m );
! 565: }
! 566: }
! 567:
! 568: void cmax(p,b)
! 569: P p;
! 570: Q *b;
! 571: {
! 572: DCP dc;
! 573: Q m,m1;
! 574: N tn;
! 575:
! 576: if ( NUM(p) ) {
! 577: tn = NM((Q)p);
! 578: NTOQ(tn,1,*b);
! 579: } else {
! 580: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
! 581: cmax(COEF(dc),&m1);
! 582: if ( cmpq(m1,m) > 0 )
! 583: m = m1;
! 584: }
! 585: *b = m;
! 586: }
! 587: }
! 588:
! 589: int nextbin(vn,n)
! 590: VN vn;
! 591: int n;
! 592: {
! 593: int tmp,i,carry;
! 594:
! 595: if ( n == 0 )
! 596: return ( 1 );
! 597:
! 598: for ( i = n - 1, carry = 1; i >= 0; i-- ) {
! 599: tmp = vn[i].n + carry;
! 600: vn[i].n = tmp % 2;
! 601: carry = tmp / 2;
! 602: }
! 603: return ( carry );
! 604: }
! 605:
! 606: void mulsgn(vn,vnt,n,vn1)
! 607: VN vn,vnt,vn1;
! 608: int n;
! 609: {
! 610: int i;
! 611:
! 612: for ( i = 0; vn[i].v; i++ )
! 613: vn1[i].n = vn[i].n;
! 614: for ( i = 0; i < n; i++ )
! 615: if ( vnt[i].n )
! 616: vn1[(int)vnt[i].v].n *= -1;
! 617: }
! 618:
! 619: void next(vn)
! 620: VN vn;
! 621: {
! 622: int i,m,n,tmp,carry;
! 623:
! 624: for ( m = 0, i = 0; vn[i].v; i++ )
! 625: m = MAX(m,ABS(vn[i].n));
! 626: if ( m == 0 ) {
! 627: vn[--i].n = 1;
! 628: return;
! 629: }
! 630: for ( carry = 1, m++, i--; ( i >= 0 ) && carry; i-- ) {
! 631: tmp = vn[i].n + carry;
! 632: vn[i].n = tmp % m;
! 633: carry = tmp / m;
! 634: }
! 635: if ( ( i == -1 ) && carry ) {
! 636: for ( i = 0; vn[i].v; i++ )
! 637: vn[i].n = 0;
! 638: vn[--i].n = m;
! 639: } else {
! 640: for ( n = 0, i = 0; vn[i].v; i++ )
! 641: n = MAX(n,ABS(vn[i].n));
! 642: if ( n < m - 1 )
! 643: vn[--i].n = m - 1;
! 644: }
! 645: }
! 646:
! 647: void clctv(vl,p,nvlp)
! 648: VL vl;
! 649: P p;
! 650: VL *nvlp;
! 651: {
! 652: int i,n;
! 653: VL tvl;
! 654: VN tvn;
! 655:
! 656: if ( !p || NUM(p) ) {
! 657: *nvlp = 0;
! 658: return;
! 659: }
! 660:
! 661: for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
! 662: tvn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
! 663: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) {
! 664: tvn[i].v = tvl->v;
! 665: tvn[i].n = 0;
! 666: }
! 667:
! 668: markv(tvn,n,p);
! 669: vntovl(tvn,n,nvlp);
! 670: }
! 671:
! 672: void markv(vn,n,p)
! 673: VN vn;
! 674: int n;
! 675: P p;
! 676: {
! 677: V v;
! 678: DCP dc;
! 679: int i;
! 680:
! 681: if ( NUM(p) )
! 682: return;
! 683: v = VR(p);
! 684: for ( i = 0, v = VR(p); i < n; i++ )
! 685: if ( v == vn[i].v ) {
! 686: vn[i].n = 1;
! 687: break;
! 688: }
! 689: for ( dc = DC(p); dc; dc = NEXT(dc) )
! 690: markv(vn,n,COEF(dc));
! 691: }
! 692:
! 693: void vntovl(vn,n,vlp)
! 694: VN vn;
! 695: int n;
! 696: VL *vlp;
! 697: {
! 698: int i;
! 699: VL tvl,tvl0;
! 700:
! 701: for ( i = 0, tvl0 = 0; ; ) {
! 702: while ( ( i < n ) && ( vn[i].n == 0 ) ) i++;
! 703: if ( i == n )
! 704: break;
! 705: else {
! 706: if ( !tvl0 ) {
! 707: NEWVL(tvl0);
! 708: tvl = tvl0;
! 709: } else {
! 710: NEWVL(NEXT(tvl));
! 711: tvl = NEXT(tvl);
! 712: }
! 713: tvl->v = vn[i++].v;
! 714: }
! 715: }
! 716: if ( tvl0 )
! 717: NEXT(tvl) = 0;
! 718: *vlp = tvl0;
! 719: }
! 720:
! 721: int dbound(v,f)
! 722: V v;
! 723: P f;
! 724: {
! 725: int m;
! 726: DCP dc;
! 727:
! 728: if ( !f )
! 729: return ( -1 );
! 730: else if ( v != VR(f) )
! 731: return homdeg(f);
! 732: else {
! 733: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) )
! 734: m = MAX(m,homdeg(COEF(dc)));
! 735: return ( m );
! 736: }
! 737: }
! 738:
! 739: int homdeg(f)
! 740: P f;
! 741: {
! 742: int m;
! 743: DCP dc;
! 744:
! 745: if ( !f )
! 746: return ( -1 );
! 747: else if ( NUM(f) )
! 748: return( 0 );
! 749: else {
! 750: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) )
! 751: m = MAX(m,QTOS(DEG(dc))+homdeg(COEF(dc)));
! 752: return ( m );
! 753: }
! 754: }
! 755:
! 756: int minhomdeg(f)
! 757: P f;
! 758: {
! 759: int m;
! 760: DCP dc;
! 761:
! 762: if ( !f )
! 763: return ( -1 );
! 764: else if ( NUM(f) )
! 765: return( 0 );
! 766: else {
! 767: for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) )
! 768: m = MIN(m,QTOS(DEG(dc))+minhomdeg(COEF(dc)));
! 769: return ( m );
! 770: }
! 771: }
! 772:
! 773: void adjc(vl,f,a,lc0,q,fr,ar)
! 774: VL vl;
! 775: P f,a,lc0;
! 776: Q q;
! 777: P *fr,*ar;
! 778: {
! 779: P m,m1;
! 780: Q t;
! 781:
! 782: invl((Q)LC(f),q,&t);
! 783: mulq((Q)lc0,t,(Q *)&m);
! 784: mulpq(f,m,&m1); cmp(q,m1,fr);
! 785: invl((Q)m,q,&t);
! 786: mulpq(a,(P)t,&m1);
! 787: cmp(q,m1,ar);
! 788: }
! 789:
! 790: #if 1
! 791: void affinemain(vl,p,v0,n,pl,pr)
! 792: VL vl;
! 793: V v0;
! 794: int n;
! 795: P *pl;
! 796: P p;
! 797: P *pr;
! 798: {
! 799: P x,t,m,c,s,a;
! 800: DCP dc;
! 801: Q d;
! 802:
! 803: if ( !p )
! 804: *pr = 0;
! 805: else if ( NUM(p) )
! 806: *pr = p;
! 807: else if ( VR(p) != v0 ) {
! 808: MKV(VR(p),x);
! 809: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 810: affinemain(vl,COEF(dc),v0,n,pl,&t);
! 811: if ( DEG(dc) ) {
! 812: pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
! 813: addp(vl,m,c,&a); c = a;
! 814: } else {
! 815: addp(vl,t,c,&a); c = a;
! 816: }
! 817: }
! 818: *pr = c;
! 819: } else {
! 820: dc = DC(p);
! 821: c = COEF(dc);
! 822: for ( d = DEG(dc), dc = NEXT(dc);
! 823: dc; d = DEG(dc), dc = NEXT(dc) ) {
! 824: mulp(vl,pl[QTOS(d)-QTOS(DEG(dc))],c,&m);
! 825: addp(vl,m,COEF(dc),&c);
! 826: }
! 827: if ( d ) {
! 828: mulp(vl,pl[QTOS(d)],c,&m); c = m;
! 829: }
! 830: *pr = c;
! 831: }
! 832: }
! 833: #endif
! 834:
! 835: #if 0
! 836: affine(vl,p,vn,r)
! 837: VL vl;
! 838: P p;
! 839: VN vn;
! 840: P *r;
! 841: {
! 842: int n,d,d1,i;
! 843: Q *q;
! 844: Q **bc;
! 845:
! 846: if ( !p || NUM(p) )
! 847: *r = p;
! 848: else {
! 849: for ( i = 0, d = 0; vn[i].v; i++ )
! 850: d1 = getdeg(vn[i].v,p), d = MAX(d,d1);
! 851: W_CALLOC(d+1,Q *,bc);
! 852: for ( i = 0; i <= d; i++ )
! 853: W_CALLOC(i+1,Q,q), mkbc(i,q), bc[i] = q;
! 854: afmain(vl,bc,p,vn,r);
! 855: }
! 856: }
! 857:
! 858: afmain(vl,bc,p,vn,r)
! 859: VL vl;
! 860: Q **bc;
! 861: P p;
! 862: VN vn;
! 863: P *r;
! 864: {
! 865: P t,s,u;
! 866: P *c,*rc;
! 867: Q *q;
! 868: DCP dc;
! 869: int n,i,j;
! 870:
! 871: if ( !p || NUM(p) || !vn[0].v )
! 872: *r = p;
! 873: else if ( vn[0].v != VR(p) ) {
! 874: for ( i = 0; vn[i].v && (vn[i].v != VR(p)); i++ );
! 875: if ( vn[i].v )
! 876: afmain(vl,bc,p,vn+i,r);
! 877: else {
! 878: n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);
! 879: for ( dc = DC(p); dc; dc = NEXT(dc) )
! 880: afmain(vl,bc,COEF(dc),vn,&c[QTOS(DEG(dc))]);
! 881: plisttop(c,VR(p),n,r);
! 882: }
! 883: } else {
! 884: n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);
! 885: W_CALLOC(n+1,P,rc); W_CALLOC(n+1,Q,q);
! 886: for ( dc = DC(p); dc; dc = NEXT(dc) )
! 887: afmain(vl,bc,COEF(dc),vn+1,&c[QTOS(DEG(dc))]);
! 888: if ( !vn[0].n )
! 889: bcopy(c,rc,sizeof(P)*(n+1));
! 890: else {
! 891: for ( q[0] = ONE, STOQ(vn[0].n,q[1]), i = 2; i <= n; i++ )
! 892: mulq(q[i-1],q[1],&q[i]);
! 893: for ( j = 0; j <= n; rc[j] = t, j++ )
! 894: for ( i = j, t = 0; i <= n; i++ )
! 895: if ( c[i] )
! 896: mulq(bc[i][j],q[i-j],&s), mulp(vl,s,c[i],&u),
! 897: addp(CO,u,t,&s), t = s;
! 898: }
! 899: plisttop(rc,VR(p),n,r);
! 900: }
! 901: }
! 902: #endif
! 903:
! 904: void restore(vl,f,vn,fr)
! 905: VL vl;
! 906: P f;
! 907: VN vn;
! 908: P *fr;
! 909: {
! 910: int i;
! 911: P vv,g,g1,t;
! 912: Q s;
! 913:
! 914: g = f;
! 915: for ( i = 0; vn[i].v; i++ ) {
! 916: MKV(vn[i].v,t);
! 917: if ( vn[i].n ) {
! 918: STOQ(-vn[i].n,s);
! 919: addp(vl,t,(P)s,&vv);
! 920: } else
! 921: vv = t;
! 922:
! 923: substp(vl,g,vn[i].v,vv,&g1); g = g1;
! 924: }
! 925: *fr = g;
! 926: }
! 927:
! 928: void mergev(vl,vl1,vl2,nvlp)
! 929: VL vl,vl1,vl2,*nvlp;
! 930: {
! 931: int i,n;
! 932: VL tvl;
! 933: VN vn;
! 934:
! 935: if ( !vl1 ) {
! 936: *nvlp = vl2; return;
! 937: } else if ( !vl2 ) {
! 938: *nvlp = vl1; return;
! 939: }
! 940: for ( i = 0, tvl = vl; tvl; i++, tvl = NEXT(tvl) );
! 941: n = i;
! 942: W_CALLOC(n,struct oVN,vn);
! 943: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ )
! 944: vn[i].v = tvl->v;
! 945: for ( i = 0, tvl = vl1; (i < n) && tvl; i++, tvl = NEXT(tvl)) {
! 946: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
! 947: i++;
! 948: if ( i == n )
! 949: break;
! 950: else
! 951: vn[i].n = 1;
! 952: }
! 953: for ( i = 0, tvl = vl2; (i < n) && tvl; i++, tvl = NEXT(tvl) ) {
! 954: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
! 955: i++;
! 956: if ( i == n )
! 957: break;
! 958: else
! 959: vn[i].n = 1;
! 960: }
! 961: vntovl(vn,n,nvlp);
! 962: }
! 963:
! 964: #if 0
! 965: void substvp(vl,f,vn,g)
! 966: VL vl;
! 967: P f;
! 968: VN vn;
! 969: P *g;
! 970: {
! 971: V v;
! 972: int i;
! 973: P h,h1;
! 974: Q t;
! 975:
! 976: h = f;
! 977: for ( i = 0; v = vn[i].v; i++ ) {
! 978: STOQ(vn[i].n,t);
! 979: substp(vl,h,v,(P)t,&h1); h = h1;
! 980: }
! 981: *g = h;
! 982: }
! 983:
! 984: void affine(vl,f,vn,fr)
! 985: VL vl;
! 986: P f;
! 987: VN vn;
! 988: P *fr;
! 989: {
! 990: int i,j,n;
! 991: P vv,g,g1,t,u;
! 992: Q s;
! 993: int *dlist;
! 994: P **plist;
! 995:
! 996: for ( n = 0; vn[n].v; n++);
! 997: dlist = (int *)ALLOCA((n+1)*sizeof(int));
! 998: plist = (P **)ALLOCA((n+1)*sizeof(P *));
! 999: for ( i = 0; vn[i].v; i++ ) {
! 1000: if ( !vn[i].n )
! 1001: continue;
! 1002: dlist[i] = getdeg(vn[i].v,f);
! 1003: plist[i] = (P *)ALLOCA((dlist[i]+1)*sizeof(P));
! 1004:
! 1005: MKV(vn[i].v,t);
! 1006: if ( vn[i].n ) {
! 1007: STOQ(vn[i].n,s); addp(vl,t,(P)s,&vv);
! 1008: } else
! 1009: vv = t;
! 1010:
! 1011: for ( j = 0, t = (P)ONE; j < dlist[i]; j++ ) {
! 1012: plist[i][j] = t;
! 1013: mulp(vl,t,vv,&u);
! 1014: t = u;
! 1015: }
! 1016: plist[i][j] = t;
! 1017: }
! 1018:
! 1019: g = f;
! 1020: for ( i = 0; vn[i].v; i++ ) {
! 1021: if ( !vn[i].n )
! 1022: continue;
! 1023: affinemain(vl,g,vn[i].v,dlist[i],plist[i],&g1); g = g1;
! 1024: }
! 1025: *fr = g;
! 1026: }
! 1027: #endif