![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to EZ.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 EZ.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/EZ.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #include "ca.h"
! 3:
! 4: void ezgcdp(vl,p1,p2,pr)
! 5: VL vl;
! 6: P p1,p2;
! 7: P *pr;
! 8: {
! 9: P f1,f2;
! 10: Q c;
! 11:
! 12: if ( !p1 )
! 13: if ( !p2 )
! 14: *pr = 0;
! 15: else
! 16: *pr = p2;
! 17: else if ( !p2 )
! 18: *pr = p1;
! 19: else {
! 20: if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
! 21: error("ezgcdp : invalid argument");
! 22: ptozp(p1,1,&c,&f1);
! 23: ptozp(p2,1,&c,&f2);
! 24: ezgcdpz(vl,f1,f2,pr);
! 25: }
! 26: }
! 27:
! 28: /*
! 29: * p1, p2 : integer coefficient
! 30: * returns gcd(pp(p1),pp(p2)) * gcd(cont(p1),cont(p2))
! 31: */
! 32:
! 33: void ezgcdpz(vl,p1,p2,pr)
! 34: VL vl;
! 35: P p1,p2,*pr;
! 36: {
! 37: P f1,f2,t,s,g,gcd;
! 38: P fc1,fc2,cont;
! 39: VL tvl,svl;
! 40: V v1,v2;
! 41: DCP dc,dc1,dc2;
! 42: N cn;
! 43: Q cq;
! 44: Q c1,c2,c;
! 45: P g1,g2;
! 46: P mgcd;
! 47: extern int nez;
! 48: P pm[2];
! 49:
! 50: if ( nez ) {
! 51: pm[0] = p1; pm[1] = p2; nezgcdnpz(vl,pm,2,pr); return;
! 52: }
! 53: monomialfctr(vl,p1,&f1,&dc1); p1 = f1;
! 54: monomialfctr(vl,p2,&f2,&dc2); p2 = f2;
! 55: for ( mgcd = (P)ONE; dc1; dc1 = NEXT(dc1) )
! 56: for ( v1 = VR(dc1->c), dc = dc2; dc; dc = NEXT(dc) )
! 57: if ( v1 == VR(dc->c) ) {
! 58: pwrp(vl,dc->c,cmpq(dc1->d,dc->d)>=0?dc->d:dc1->d,&t);
! 59: mulp(vl,mgcd,t,&s); mgcd = s; break;
! 60: }
! 61:
! 62: if ( NUM(p1) ) {
! 63: if ( NUM(p2) ) {
! 64: gcdn(NM((Q)p1),NM((Q)p2),&cn); NTOQ(cn,1,c);
! 65: } else {
! 66: ptozp(p2,1,&c2,&t);
! 67: gcdn(NM((Q)p1),NM(c2),&cn); NTOQ(cn,1,c);
! 68: }
! 69: mulp(vl,(P)c,mgcd,pr);
! 70: return;
! 71: } else if ( NUM(p2) ) {
! 72: ptozp(p1,1,&c1,&t);
! 73: gcdn(NM(c1),NM((Q)p2),&cn); NTOQ(cn,1,c);
! 74: mulp(vl,(P)c,mgcd,pr);
! 75: return;
! 76: }
! 77:
! 78: ptozp(p1,1,&c1,&g1); ptozp(p2,1,&c2,&g2);
! 79: gcdn(NM(c1),NM(c2),&cn); NTOQ(cn,1,cq);
! 80:
! 81: if ( ( v1 = VR(f1) ) != ( v2 = VR(f2) ) ) {
! 82: while ( v1 != VR(vl) && v2 != VR(vl) )
! 83: vl = NEXT(vl);
! 84: if ( v1 == VR(vl) ) {
! 85: pcp(vl,f1,&g,&fc1);
! 86: ezgcdpz(vl,fc1,f2,&t);
! 87: } else {
! 88: pcp(vl,f2,&g,&fc2);
! 89: ezgcdpz(vl,fc2,f1,&t);
! 90: }
! 91: mulp(vl,t,(P)cq,&s); mulp(vl,s,mgcd,pr);
! 92: return;
! 93: }
! 94:
! 95: pcp(vl,f1,&g1,&fc1); pcp(vl,f2,&g2,&fc2);
! 96: ezgcdpz(vl,fc1,fc2,&cont);
! 97: clctv(vl,g1,&tvl); clctv(vl,g2,&svl);
! 98: if ( !NEXT(tvl) && !NEXT(svl) ) {
! 99: uezgcdpz(vl,g1,g2,&t);
! 100: mulp(vl,t,cont,&s); mulp(vl,s,(P)cq,&t); mulp(vl,t,mgcd,pr);
! 101: return;
! 102: }
! 103:
! 104: if ( homdeg(g1) > homdeg(g2) ) {
! 105: t = g1; g1 = g2; g2 = t;
! 106: }
! 107: sqfrp(vl,g1,&dc);
! 108: for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) {
! 109: if ( NUM(COEF(dc)) )
! 110: continue;
! 111: ezgcdpp(vl,dc,g,&t);
! 112: if ( NUM(t) )
! 113: continue;
! 114: mulp(vl,gcd,t,&s); gcd = s;
! 115: divsp(vl,g,t,&s); g = s;
! 116: }
! 117: mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,&s); mulp(vl,s,mgcd,pr);
! 118: }
! 119:
! 120: void uezgcdpz(vl,p1,p2,pr)
! 121: VL vl;
! 122: P p1,p2,*pr;
! 123: {
! 124: P g1,g2,gcd,t,s,g;
! 125: DCP dc;
! 126: N cn;
! 127: Q c1,c2,cq;
! 128: P f1,f2;
! 129:
! 130: if ( NUM(p1) ) {
! 131: if ( NUM(p2) ) {
! 132: gcdn(NM((Q)p1),NM((Q)p2),&cn); NTOQ(cn,1,cq),*pr = (P)cq;
! 133: return;
! 134: } else {
! 135: ptozp(p2,1,&c2,&t);
! 136: gcdn(NM((Q)p1),NM(c2),&cn); NTOQ(cn,1,cq),*pr = (P)cq;
! 137: return;
! 138: }
! 139: } else if ( NUM(p2) ) {
! 140: ptozp(p1,1,&c1,&t);
! 141: gcdn(NM(c1),NM((Q)p2),&cn); NTOQ(cn,1,cq),*pr = (P)cq;
! 142: return;
! 143: }
! 144:
! 145: ptozp(p1,1,&c1,&f1); ptozp(p2,1,&c2,&f2);
! 146: gcdn(NM(c1),NM(c2),&cn); NTOQ(cn,1,cq);
! 147: if ( UDEG(f1) > UDEG(f2) ) {
! 148: g1 = f2; g2 = f1;
! 149: } else {
! 150: g1 = f1; g2 = f2;
! 151: }
! 152: usqp(g1,&dc);
! 153: for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) {
! 154: if ( NUM(COEF(dc)) )
! 155: continue;
! 156: uezgcdpp(dc,g,&t);
! 157: if ( NUM(t) )
! 158: continue;
! 159: mulp(CO,gcd,t,&s); gcd = s;
! 160: divsp(CO,g,t,&s); g = s;
! 161: }
! 162: mulp(vl,gcd,(P)cq,pr);
! 163: }
! 164:
! 165: /*
! 166: dc : dc pair c^d
! 167: r = gcd(c^d,f)
! 168: */
! 169:
! 170: void ezgcdpp(vl,dc,f,r)
! 171: VL vl;
! 172: DCP dc;
! 173: P f;
! 174: P *r;
! 175: {
! 176: int k;
! 177: P g,h,t,s,gcd;
! 178:
! 179: if ( NUM(f) ) {
! 180: *r = (P)ONE;
! 181: return;
! 182: }
! 183: k = QTOS(DEG(dc)) - 1;
! 184: /* ezgcd1p(vl,COEF(dc),f,&gcd); */
! 185: ezgcdnp(vl,COEF(dc),&f,1,&gcd);
! 186: g = gcd; divsp(vl,f,gcd,&h);
! 187: for ( ; k > 0; k-- ) {
! 188: /* ezgcd1p(vl,g,h,&t); */
! 189: ezgcdnp(vl,g,&h,1,&t);
! 190: if ( NUM(t) )
! 191: break;
! 192: mulp(vl,gcd,t,&s); gcd = s;
! 193: divsp(vl,h,t,&s); h = s;
! 194: }
! 195: *r = gcd;
! 196: }
! 197:
! 198: void ezgcd1p(vl,p0,p1,gcdp)
! 199: VL vl;
! 200: P p0,p1;
! 201: P *gcdp;
! 202: {
! 203: VL nvl,tvl,vl0,vl1,vl2;
! 204: P f0,f1,t;
! 205:
! 206: clctv(vl,p0,&vl0); clctv(vl,p1,&vl1); mergev(vl,vl0,vl1,&vl2);
! 207: minchdegp(vl,p0,&tvl,&t); reordvar(vl2,tvl->v,&nvl);
! 208: reorderp(nvl,vl,p0,&f0); reorderp(nvl,vl,p1,&f1);
! 209: ezgcdnp(nvl,f0,&f1,1,&t);
! 210: reorderp(vl,nvl,t,gcdp);
! 211: }
! 212:
! 213: void uezgcdpp(dc,f,r)
! 214: DCP dc;
! 215: P f;
! 216: P *r;
! 217: {
! 218: int k;
! 219: P g,h,t,s,gcd;
! 220:
! 221: if ( NUM(f) ) {
! 222: *r = (P)ONE;
! 223: return;
! 224: }
! 225:
! 226: k = QTOS(DEG(dc)) - 1;
! 227: uezgcd1p(COEF(dc),f,&gcd);
! 228: g = gcd; divsp(CO,f,gcd,&h);
! 229: for ( ; k > 0; k-- ) {
! 230: uezgcd1p(g,h,&t);
! 231: if ( NUM(t) )
! 232: break;
! 233: mulp(CO,gcd,t,&s); gcd = s;
! 234: divsp(CO,h,t,&s); h = s;
! 235: }
! 236: *r = gcd;
! 237: }
! 238:
! 239: /*
! 240: * pr = gcd(p0,ps[0],...,ps[m-1])
! 241: *
! 242: * p0 is already square-free and primitive.
! 243: * ps[i] is at least primitive.
! 244: *
! 245: */
! 246:
! 247: void ezgcdnp(vl,p0,ps,m,pr)
! 248: VL vl;
! 249: int m;
! 250: P p0,*ps,*pr;
! 251: {
! 252: /* variables */
! 253: P p00,g,h,g0,h0,a0,b0;
! 254: P lp0,lgp0,lp00,lg,lg0,cbd,tq,t;
! 255: P *tps;
! 256: VL nvl1,nvl2,nvl,tvl;
! 257: V v;
! 258: int i,j,k,d0,dg,dg0,dmin;
! 259: VN vn0,vn1,vnt;
! 260: int nv,flag;
! 261:
! 262: /* set-up */
! 263: if ( NUM(p0) ) {
! 264: *pr = (P) ONE;
! 265: return;
! 266: }
! 267: for ( v = VR(p0), i = 0; i < m; i++ )
! 268: if ( NUM(ps[i]) || (v != VR(ps[i])) ) {
! 269: *pr = (P)ONE;
! 270: return;
! 271: }
! 272: tps = (P *) ALLOCA(m*sizeof(P));
! 273: for ( i = 0; i < m; i++ )
! 274: tps[i] = ps[i];
! 275: sortplist(tps,m);
! 276: /* deg(tps[0]) <= deg(tps[1]) <= ... */
! 277:
! 278: if ( !cmpq(DEG(DC(p0)),ONE) ) {
! 279: if ( divcheck(vl,tps,m,(P)ONE,p0) )
! 280: *pr = p0;
! 281: else
! 282: *pr = (P)ONE;
! 283: return;
! 284: }
! 285:
! 286: lp0 = LC(p0); lg = lp0; dmin = d0 = deg(v,p0);
! 287: clctv(vl,p0,&nvl);
! 288: for ( i = 0; i < m; i++ ) {
! 289: clctv(vl,tps[i],&nvl1); mergev(vl,nvl,nvl1,&nvl2);
! 290: nvl = nvl2;
! 291: ezgcdpz(nvl,lg,LC(tps[i]),&t); lg = t;
! 292: }
! 293:
! 294: mulp(nvl,p0,lg,&lgp0);
! 295: k = dbound(v,lgp0) + 1;
! 296: cbound(nvl,lgp0,(Q *)&cbd);
! 297: for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++ );
! 298: W_CALLOC(nv,struct oVN,vn0); W_CALLOC(nv,struct oVN,vnt);
! 299: W_CALLOC(nv,struct oVN,vn1);
! 300: for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
! 301: vn1[i].v = vn0[i].v = tvl->v;
! 302:
! 303: /* main loop */
! 304: for ( dg = deg(v,tps[0]) + 1; ; next(vn0) )
! 305: do {
! 306: for ( i = 0, j = 0; vn0[i].v; i++ )
! 307: if ( vn0[i].n ) vnt[j++].v = (V)i;
! 308: vnt[j].n = 0;
! 309:
! 310: /* find b s.t. LC(p0)(b), LC(tps[i])(b) != 0 */
! 311: mulsgn(vn0,vnt,j,vn1);
! 312: substvp(nvl,p0,vn1,&p00);
! 313: flag = (!zerovpchk(nvl,lp0,vn1) && sqfrchk(p00));
! 314: for ( i = 0; flag && (i < m); i++ )
! 315: flag &= (!zerovpchk(nvl,LC(tps[i]),vn1));
! 316: if ( !flag )
! 317: continue;
! 318:
! 319: /* substitute y -> b */
! 320: substvp(nvl,lg,vn1,&lg0);
! 321: lp00 = LC(p00);
! 322: /* extended-gcd in 1-variable */
! 323: uexgcdnp(nvl,p00,tps,m,vn1,(Q)cbd,&g0,&h0,&a0,&b0,(Q *)&tq);
! 324: if ( NUM(g0) ) {
! 325: *pr = (P)ONE;
! 326: return;
! 327: } else if ( dg > ( dg0 = deg(v,g0) ) ) {
! 328: dg = dg0;
! 329: if ( dg == d0 ) {
! 330: if ( divcheck(nvl,tps,m,lp0,p0) ) {
! 331: *pr = p0;
! 332: return;
! 333: }
! 334: } else if ( dg == deg(v,tps[0]) ) {
! 335: if ( divtpz(nvl,p0,tps[0],&t) &&
! 336: divcheck(nvl,tps+1,m-1,LC(tps[0]),tps[0]) ) {
! 337: *pr = tps[0];
! 338: return;
! 339: } else
! 340: break;
! 341: } else {
! 342: henmv(nvl,vn1,lgp0,g0,h0,a0,b0,
! 343: lg,lp0,lg0,lp00,(Q)tq,k,&g,&h);
! 344: if ( divtpz(nvl,lgp0,g,&t) &&
! 345: divcheck(nvl,tps,m,lg,g) ) {
! 346: pcp(nvl,g,pr,&t);
! 347: return;
! 348: }
! 349: }
! 350: }
! 351: } while ( !nextbin(vnt,j) );
! 352: }
! 353:
! 354: /*
! 355: f : sqfr
! 356: */
! 357:
! 358: void uezgcd1p(f,g,r)
! 359: P f,g;
! 360: P *r;
! 361: {
! 362: UM wf,wf1,wf2,wg,wg1,wgcd,wcof;
! 363: int d1,d2,dg,m,index,n;
! 364: ML list;
! 365: DCP dc;
! 366: P t;
! 367: Q c;
! 368:
! 369: if ( NUM(f) || NUM(g) ) {
! 370: *r = (P)ONE;
! 371: return;
! 372: }
! 373: ptozp(f,1,&c,&t); f = t; ptozp(g,1,&c,&t); g = t;
! 374: d1 = UDEG(f); d2 = UDEG(g);
! 375: n = MAX(d1,d2)+1;
! 376: wf = W_UMALLOC(n); wf1 = W_UMALLOC(n);
! 377: wf2 = W_UMALLOC(n); wg = W_UMALLOC(n);
! 378: wg1 = W_UMALLOC(n); wgcd = W_UMALLOC(n);
! 379: wcof = W_UMALLOC(n);
! 380:
! 381: for ( index = INDEX, dg = MIN(d1,d2)+1; ; ) {
! 382: m = sprime[index++];
! 383: if ( !rem(NM((Q)UCOEF(f)),m) || !rem(NM((Q)UCOEF(g)),m))
! 384: continue;
! 385: ptoum(m,f,wf); cpyum(wf,wf1);
! 386: diffum(m,wf1,wf2); gcdum(m,wf1,wf2,wgcd);
! 387: if ( DEG(wgcd) > 0 )
! 388: continue;
! 389: ptoum(m,g,wg); cpyum(wg,wg1); cpyum(wf,wf1);
! 390: gcdum(m,wf1,wg1,wgcd);
! 391: if ( dg <= DEG(wgcd) )
! 392: continue;
! 393: else
! 394: dg = DEG(wgcd);
! 395: if ( dg == 0 ) {
! 396: *r = (P)ONE;
! 397: return;
! 398: } else if ( dg == d1 ) {
! 399: if ( divtpz(CO,g,f,&t) ) {
! 400: *r = f;
! 401: return;
! 402: }
! 403: } else if ( dg == d2 ) {
! 404: if ( divtpz(CO,f,g,&t) ) {
! 405: *r = g;
! 406: return;
! 407: }
! 408: } else {
! 409: divum(m,wf,wgcd,wcof);
! 410: ezgcdhensel(f,m,wcof,wgcd,&list);
! 411: dtest(f,list,1,&dc);
! 412: if ( NEXT(dc) ) {
! 413: if ( divtpz(CO,g,COEF(dc),&t) ) {
! 414: *r = COEF(dc);
! 415: return;
! 416: }
! 417: }
! 418: }
! 419: }
! 420: }
! 421:
! 422: void uexgcdnp(vl,p1,ps,n,vn,cbd,g,h,a,b,q)
! 423: VL vl;
! 424: P p1;
! 425: P *ps;
! 426: int n;
! 427: VN vn;
! 428: Q cbd;
! 429: P *g,*h,*a,*b;
! 430: Q *q;
! 431: {
! 432: UM wg,wh,wg1,wh1,wgcd,wt;
! 433: P t,s,gcd,cof,gk,hk,ak,bk;
! 434: Q c,tq,bb,tq1,qk;
! 435: int mod,k,i,index,d;
! 436:
! 437: ptozp(p1,1,&c,&gcd);
! 438: for ( i = 0; i < n; i++ ) {
! 439: substvp(vl,ps[i],vn,&t);
! 440: uezgcd1p(gcd,t,&s);
! 441: if ( NUM(s) ) {
! 442: *g = (P)ONE; *h = p1; *a = 0; *b = (P)ONE;
! 443: return;
! 444: } else
! 445: gcd = s;
! 446: }
! 447:
! 448: if ( !dcomp(p1,gcd) ) {
! 449: *g = p1; *h = (P)ONE; *a = 0; *b = (P)ONE;
! 450: return;
! 451: }
! 452:
! 453: divsp(CO,p1,gcd,&cof);
! 454: d = UDEG(p1)+1;
! 455: wg = W_UMALLOC(d); wh = W_UMALLOC(d);
! 456: wg1 = W_UMALLOC(d); wh1 = W_UMALLOC(d);
! 457: wgcd = W_UMALLOC(d); wt = W_UMALLOC(d);
! 458: for ( index = INDEX; ; ) {
! 459: mod = sprime[index++];
! 460: if ( !rem(NM((Q)LC(p1)),mod) )
! 461: continue;
! 462: ptoum(mod,gcd,wg); ptoum(mod,cof,wh);
! 463: cpyum(wg,wg1); cpyum(wh,wh1);
! 464: gcdum(mod,wg,wh,wt); cpyum(wt,wgcd);
! 465: if ( DEG(wgcd) > 0 )
! 466: continue;
! 467: STOQ(mod,tq); addq(cbd,cbd,&bb);
! 468: for ( k = 1; cmpq(tq,bb) < 0; k++ ) {
! 469: mulq(tq,tq,&tq1); tq = tq1;
! 470: }
! 471: henzq(p1,gcd,wg1,cof,wh1,mod,k,&gk,&hk,&bk,&ak,&qk);
! 472: break;
! 473: }
! 474: *g = gk; *h = hk; *a = ak; *b = bk; *q = tq;
! 475: }
! 476:
! 477: void ezgcdhensel(f,mod,cof,gcd,listp)
! 478: P f;
! 479: int mod;
! 480: UM gcd,cof;
! 481: ML *listp;
! 482: {
! 483: register int i,j;
! 484: int q,n,bound;
! 485: int *p;
! 486: int **pp;
! 487: ML blist,clist,bqlist,cqlist,rlist;
! 488: UM *b;
! 489: LUM fl,tl;
! 490: LUM *l;
! 491: LUM LUMALLOC();
! 492:
! 493: blist = W_MLALLOC(2);
! 494: blist->n = 2; blist->mod = mod;
! 495: blist->c[0] = (pointer)cof; blist->c[1] = (pointer)gcd;
! 496: gcdgen(f,blist,&clist);
! 497: henprep(f,blist,clist,&bqlist,&cqlist);
! 498: n = bqlist->n; q = bqlist->mod;
! 499: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 500:
! 501: if ( bound == 1 ) {
! 502: *listp = rlist = MLALLOC(n);
! 503: rlist->n = n;
! 504: rlist->mod = q;
! 505: rlist->bound = bound;
! 506:
! 507: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c;
! 508: i < n; i++ ) {
! 509: tl = LUMALLOC(DEG(b[i]),1);
! 510: l[i] = tl;
! 511: p = COEF(b[i]);
! 512:
! 513: for ( j = 0, pp = COEF(tl);
! 514: j <= DEG(tl); j++ )
! 515: pp[j][0] = p[j];
! 516: }
! 517: } else {
! 518: W_LUMALLOC((int)UDEG(f),bound,fl);
! 519: ptolum(q,bound,f,fl);
! 520: henmain(fl,bqlist,cqlist,listp);
! 521: }
! 522: }
! 523:
! 524: void ezgcdnpz(vl,ps,m,pr)
! 525: VL vl;
! 526: P *ps,*pr;
! 527: int m;
! 528: {
! 529: P t,s,gcd;
! 530: P cont;
! 531: VL tvl,svl,nvl;
! 532: int i;
! 533: DCP dc;
! 534: N cn;
! 535: Q cq;
! 536: P *pl,*ppl,*pcl;
! 537: Q *cl;
! 538: V v;
! 539:
! 540: pl = (P *)ALLOCA((m+1)*sizeof(P));
! 541: cl = (Q *)ALLOCA((m+1)*sizeof(Q));
! 542: for ( i = 0; i < m; i++ )
! 543: ptozp(ps[i],1,&cl[i],&pl[i]);
! 544: for ( i = 1, cq = cl[0]; i < m; i++ ) {
! 545: gcdn(NM(cl[i]),NM(cq),&cn); NTOQ(cn,1,cq);
! 546: }
! 547: for ( i = 0; i < m; i++ )
! 548: if ( NUM(pl[i]) ) {
! 549: *pr = (P)cq;
! 550: return;
! 551: }
! 552:
! 553: for ( i = 0, nvl = 0; i < m; i++ ) {
! 554: clctv(vl,ps[i],&tvl); mergev(vl,nvl,tvl,&svl); nvl = svl;
! 555: }
! 556:
! 557: ppl = (P *)ALLOCA((m+1)*sizeof(P));
! 558: pcl = (P *)ALLOCA((m+1)*sizeof(P));
! 559: for ( i = 0, v = nvl->v; i < m; i++ )
! 560: if ( v == VR(pl[i]) )
! 561: pcp(nvl,pl[i],&ppl[i],&pcl[i]);
! 562: else {
! 563: ppl[i] = (P)ONE; pcl[i] = pl[i];
! 564: }
! 565: ezgcdnpz(nvl,pcl,m,&cont);
! 566:
! 567: for ( i = 0; i < m; i++ )
! 568: if ( NUM(ppl[i]) ) {
! 569: mulp(nvl,cont,(P)cq,pr);
! 570: return;
! 571: }
! 572: sortplist(ppl,m);
! 573: sqfrp(nvl,ppl[0],&dc);
! 574: for ( dc = NEXT(dc), gcd = (P)ONE; dc; dc = NEXT(dc) ) {
! 575: if ( NUM(COEF(dc)) )
! 576: continue;
! 577: ezgcdnpp(vl,dc,ppl+1,m-1,&t);
! 578: if ( NUM(t) )
! 579: continue;
! 580: mulp(vl,gcd,t,&s); gcd = s;
! 581: }
! 582: mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,pr);
! 583: }
! 584:
! 585: void ezgcdnpp(vl,dc,pl,m,r)
! 586: VL vl;
! 587: DCP dc;
! 588: P *pl;
! 589: int m;
! 590: P *r;
! 591: {
! 592: int i,k;
! 593: P g,t,s,gcd;
! 594: P *pm;
! 595:
! 596: ezgcdnp(vl,COEF(dc),pl,m,&gcd);
! 597: if ( NUM(gcd) ) {
! 598: *r = (P)ONE;
! 599: return;
! 600: }
! 601: pm = (P *) ALLOCA((m+1)*sizeof(P));
! 602: for ( i = 0; i < m; i++ ) {
! 603: divsp(vl,pl[i],gcd,&pm[i]);
! 604: if ( NUM(pm[i]) ) {
! 605: *r = gcd;
! 606: return;
! 607: }
! 608: }
! 609: for ( g = gcd, k = QTOS(DEG(dc)) - 1; k > 0; k-- ) {
! 610: ezgcdnp(vl,g,pm,m,&t);
! 611: if ( NUM(t) )
! 612: break;
! 613: mulp(vl,gcd,t,&s); gcd = s;
! 614: for ( i = 0; i < m; i++ ) {
! 615: divsp(vl,pm[i],t,&s);
! 616: if ( NUM(s) ) {
! 617: *r = gcd;
! 618: return;
! 619: }
! 620: pm[i] = s;
! 621: }
! 622: }
! 623: *r = gcd;
! 624: }
! 625:
! 626: void pcp(vl,p,pp,c)
! 627: VL vl;
! 628: P p;
! 629: P *pp,*c;
! 630: {
! 631: P f,g,n;
! 632: DCP dc;
! 633: P *pl;
! 634: int i,m;
! 635: extern int nez;
! 636:
! 637: if ( NUM(p) ) {
! 638: *c = p;
! 639: *pp = (P)ONE;
! 640: } else {
! 641: for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
! 642: if ( m == 1 ) {
! 643: *c = COEF(DC(p));
! 644: divsp(vl,p,*c,pp);
! 645: return;
! 646: }
! 647: ptozp(p,1,(Q *)&n,&f);
! 648: pl = (P *)ALLOCA((m+1)*sizeof(P));
! 649: for ( i = 0, dc = DC(f); i < m; dc = NEXT(dc), i++ )
! 650: pl[i] = COEF(dc);
! 651: if ( nez )
! 652: nezgcdnpz(vl,pl,m,&g);
! 653: else
! 654: ezgcdnpz(vl,pl,m,&g);
! 655: mulp(vl,g,n,c); divsp(vl,f,g,pp);
! 656: }
! 657: }
! 658:
! 659: #if 0
! 660: pcp(vl,p,pp,c)
! 661: VL vl;
! 662: P p;
! 663: P *pp,*c;
! 664: {
! 665: P f,g,n,t;
! 666: DCP dc;
! 667:
! 668: if ( NUM(p) ) {
! 669: *c = p;
! 670: *pp = (P)ONE;
! 671: } else {
! 672: ptozp(p,1,&n,&f);
! 673: for ( g = LC(f), dc = NEXT(DC(f)); dc; dc = NEXT(dc) ) {
! 674: ezgcdpz(vl,g,COEF(dc),&t); g = t;
! 675: }
! 676: mulp(vl,g,n,c);
! 677: divsp(vl,f,g,pp);
! 678: }
! 679: }
! 680: #endif
! 681:
! 682: int divcheck(vl,ps,m,l,p)
! 683: VL vl;
! 684: int m;
! 685: P *ps,l,p;
! 686: {
! 687: int i;
! 688: P r,t;
! 689:
! 690: for ( i = 0; i < m; i++ ) {
! 691: mulp(vl,ps[i],l,&t);
! 692: if ( !divtpz(vl,t,p,&r) )
! 693: return ( 0 );
! 694: }
! 695: return ( 1 );
! 696: }
! 697:
! 698: void sortplist(plist,n)
! 699: P *plist;
! 700: int n;
! 701: {
! 702: int i,j;
! 703: P t;
! 704:
! 705: for ( i = 0; i < n; i++ )
! 706: for ( j = i + 1; j < n; j++ )
! 707: if ( cmpq(DEG(DC(plist[i])),DEG(DC(plist[j]))) > 0 ) {
! 708: t = plist[i]; plist[i] = plist[j]; plist[j] = t;
! 709: }
! 710: }
! 711:
! 712: int dcomp(p1,p2)
! 713: P p1,p2;
! 714: {
! 715: return ( cmpq(DEG(DC(p1)),DEG(DC(p2))) );
! 716: }