![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to F.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 F.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/F.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #include "ca.h"
! 3: #include <math.h>
! 4:
! 5: void homfctr();
! 6:
! 7: void fctrp(vl,f,dcp)
! 8: VL vl;
! 9: P f;
! 10: DCP *dcp;
! 11: {
! 12: VL nvl;
! 13: DCP dc;
! 14:
! 15: if ( !f || NUM(f) ) {
! 16: NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
! 17: NEXT(dc) = 0; *dcp = dc;
! 18: return;
! 19: } else if ( !qpcheck((Obj)f) )
! 20: error("fctrp : invalid argument");
! 21: else {
! 22: clctv(vl,f,&nvl);
! 23: if ( !NEXT(nvl) )
! 24: ufctr(f,1,dcp);
! 25: else
! 26: mfctr(nvl,f,dcp);
! 27: }
! 28: }
! 29:
! 30: void homfctr(vl,g,dcp)
! 31: VL vl;
! 32: P g;
! 33: DCP *dcp;
! 34: {
! 35: P s,t,u,f,h,x;
! 36: Q e;
! 37: int d,d0;
! 38: DCP dc,dct;
! 39:
! 40: substp(vl,g,vl->v,(P)ONE,&t); fctrp(vl,t,&dc); MKV(vl->v,x);
! 41: for ( dct = dc; dct; dct = NEXT(dct) )
! 42: if ( !NUM(dc->c) ) {
! 43: for ( s = 0, f = dc->c, d = d0 = homdeg(f); f; d = homdeg(f) ) {
! 44: exthp(vl,f,d,&h); STOQ(d0-d,e); pwrp(vl,x,e,&t);
! 45: mulp(vl,t,h,&u); addp(vl,s,u,&t); s = t;
! 46: subp(vl,f,h,&u); f = u;
! 47: }
! 48: dc->c = s;
! 49: }
! 50: *dcp = dc;
! 51: }
! 52:
! 53: void mfctr(vl,f,dcp)
! 54: VL vl;
! 55: P f;
! 56: DCP *dcp;
! 57: {
! 58: DCP dc,dc0,dct,dcs,dcr;
! 59: P p,pmin,ppmin,cmin,t;
! 60: VL mvl;
! 61: Q c;
! 62:
! 63: ptozp(f,1,&c,&p);
! 64: NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
! 65: msqfr(vl,p,&dct);
! 66: for ( ; dct; dct = NEXT(dct) ) {
! 67: mindegp(vl,COEF(dct),&mvl,&pmin);
! 68: #if 0
! 69: minlcdegp(vl,COEF(dct),&mvl,&pmin);
! 70: min_common_vars_in_coefp(vl,COEF(dct),&mvl,&pmin);
! 71: #endif
! 72: pcp(mvl,pmin,&ppmin,&cmin);
! 73: if ( !NUM(cmin) ) {
! 74: mfctrmain(mvl,cmin,&dcs);
! 75: for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
! 76: DEG(dcr) = DEG(dct);
! 77: reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
! 78: }
! 79: for ( ; NEXT(dc); dc = NEXT(dc) );
! 80: NEXT(dc) = dcs;
! 81: }
! 82: mfctrmain(mvl,ppmin,&dcs);
! 83: for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
! 84: DEG(dcr) = DEG(dct);
! 85: reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
! 86: }
! 87: for ( ; NEXT(dc); dc = NEXT(dc) );
! 88: NEXT(dc) = dcs;
! 89: }
! 90: adjsgn(f,dc0); *dcp = dc0;
! 91: }
! 92:
! 93: void adjsgn(p,dc)
! 94: P p;
! 95: DCP dc;
! 96: {
! 97: int sgn;
! 98: DCP dct;
! 99: P c;
! 100:
! 101: for ( dct = dc, sgn = headsgn(p); dct; dct = NEXT(dct) )
! 102: if ( !EVENN(NM(dct->d)) )
! 103: sgn *= headsgn(COEF(dct));
! 104: if ( sgn < 0 ) {
! 105: chsgnp(COEF(dc),&c); COEF(dc) = c;
! 106: }
! 107: }
! 108:
! 109: int headsgn(p)
! 110: P p;
! 111: {
! 112: if ( !p )
! 113: return 0;
! 114: else {
! 115: while ( !NUM(p) )
! 116: p = COEF(DC(p));
! 117: return SGN((Q)p);
! 118: }
! 119: }
! 120:
! 121: void fctrwithmvp(vl,f,v,dcp)
! 122: VL vl;
! 123: P f;
! 124: V v;
! 125: DCP *dcp;
! 126: {
! 127: VL nvl;
! 128: DCP dc;
! 129:
! 130: if ( NUM(f) ) {
! 131: NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
! 132: NEXT(dc) = 0; *dcp = dc;
! 133: return;
! 134: }
! 135:
! 136: clctv(vl,f,&nvl);
! 137: if ( !NEXT(nvl) )
! 138: ufctr(f,1,dcp);
! 139: else
! 140: mfctrwithmv(nvl,f,v,dcp);
! 141: }
! 142:
! 143: void mfctrwithmv(vl,f,v,dcp)
! 144: VL vl;
! 145: P f;
! 146: V v;
! 147: DCP *dcp;
! 148: {
! 149: DCP dc,dc0,dct,dcs,dcr;
! 150: P p,pmin,ppmin,cmin,t;
! 151: VL mvl;
! 152: Q c;
! 153:
! 154: ptozp(f,1,&c,&p);
! 155: NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
! 156: msqfr(vl,p,&dct);
! 157: for ( ; dct; dct = NEXT(dct) ) {
! 158: if ( !v )
! 159: mindegp(vl,COEF(dct),&mvl,&pmin);
! 160: else {
! 161: reordvar(vl,v,&mvl); reorderp(mvl,vl,COEF(dct),&pmin);
! 162: }
! 163: pcp(mvl,pmin,&ppmin,&cmin);
! 164: if ( !NUM(cmin) ) {
! 165: mfctrmain(mvl,cmin,&dcs);
! 166: for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
! 167: DEG(dcr) = DEG(dct);
! 168: reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
! 169: }
! 170: for ( ; NEXT(dc); dc = NEXT(dc) );
! 171: NEXT(dc) = dcs;
! 172: }
! 173: mfctrmain(mvl,ppmin,&dcs);
! 174: for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
! 175: DEG(dcr) = DEG(dct);
! 176: reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
! 177: }
! 178: for ( ; NEXT(dc); dc = NEXT(dc) );
! 179: NEXT(dc) = dcs;
! 180: }
! 181: *dcp = dc0;
! 182: }
! 183:
! 184: void ufctr(f,hint,dcp)
! 185: P f;
! 186: int hint;
! 187: DCP *dcp;
! 188: {
! 189: P p,c;
! 190: DCP dc,dct,dcs,dcr;
! 191:
! 192: ptozp(f,SGN((Q)UCOEF(f)),(Q *)&c,&p);
! 193: NEWDC(dc); COEF(dc) = c; DEG(dc) = ONE; NEXT(dc) = 0;
! 194: usqp(p,&dct);
! 195: for ( *dcp = dc; dct; dct = NEXT(dct) ) {
! 196: ufctrmain(COEF(dct),hint,&dcs);
! 197: for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
! 198: DEG(dcr) = DEG(dct);
! 199: for ( ; NEXT(dc); dc = NEXT(dc) );
! 200: NEXT(dc) = dcs;
! 201: }
! 202: }
! 203:
! 204: void mfctrmain(vl,p,dcp)
! 205: VL vl;
! 206: P p;
! 207: DCP *dcp;
! 208: {
! 209: int i,j,k,*win,np;
! 210: VL nvl,tvl;
! 211: VN vn,vnt,vn1;
! 212: P p0,f,g,f0,g0,s,t,lp,m;
! 213: P *fp0,*fpt,*l,*tl;
! 214: DCP dc,dc0,dcl;
! 215: int count,nv;
! 216:
! 217: if ( !cmpq(DEG(DC(p)),ONE) ) {
! 218: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 219: *dcp = dc; return;
! 220: }
! 221: lp = LC(p); fctrp(vl,lp,&dcl); clctv(vl,p,&nvl);
! 222: for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
! 223: W_CALLOC(nv,struct oVN,vn); W_CALLOC(nv,struct oVN,vnt);
! 224: W_CALLOC(nv,struct oVN,vn1);
! 225: for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
! 226: vn1[i].v = vn[i].v = tvl->v;
! 227: vn1[i].v = vn[i].v = 0;
! 228: count = 0;
! 229: while ( 1 ) {
! 230: while ( 1 ) {
! 231: count++;
! 232: for ( i = 0, j = 0; vn[i].v; i++ )
! 233: if ( vn[i].n )
! 234: vnt[j++].v = (V)i;
! 235: vnt[j].n = 0; mulsgn(vn,vnt,j,vn1);
! 236: for ( i = 0; vn1[i].v; i++ )
! 237: if ( vn1[i].n ) {
! 238: if ( vn1[i].n > 0 )
! 239: vn1[i].n = sprime[vn1[i].n];
! 240: else
! 241: vn1[i].n = -sprime[-vn1[i].n];
! 242: }
! 243: if ( valideval(nvl,dcl,vn1) ) {
! 244: substvp(nvl,p,vn1,&p0);
! 245: if ( sqfrchk(p0) ) {
! 246: ufctr(p0,1,&dc0);
! 247: if ( NEXT(NEXT(dc0)) == 0 ) {
! 248: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 249: *dcp = dc;
! 250: return;
! 251: } else
! 252: goto MAIN;
! 253: }
! 254: }
! 255: if ( nextbin(vnt,j) )
! 256: break;
! 257: }
! 258: next(vn);
! 259: }
! 260: MAIN :
! 261: #if 0
! 262: for ( i = 0; vn1[i].v; i++ )
! 263: fprintf(stderr,"%d ",vn1[i].n);
! 264: fprintf(stderr,"\n");
! 265: #endif
! 266: for ( np = 0, dc = NEXT(dc0); dc; dc = NEXT(dc), np++ );
! 267: fp0 = (P *)ALLOCA((np + 1)*sizeof(P));
! 268: fpt = (P *)ALLOCA((np + 1)*sizeof(P));
! 269: l = tl = (P *)ALLOCA((np + 1)*sizeof(P));
! 270: for ( i = 1, dc = NEXT(dc0); dc; i++, dc = NEXT(dc) )
! 271: fp0[i-1] = COEF(dc);
! 272: #if 0
! 273: sort_by_deg(np,fp0,fpt);
! 274: sort_by_deg_rev(np-1,fpt+1,fp0+1);
! 275: #endif
! 276: #if 0
! 277: sort_by_deg_rev(np,fp0,fpt); fp0[0] = fpt[0];
! 278: sort_by_deg(np-1,fpt+1,fp0+1); fp0[np] = 0;
! 279: #endif
! 280: fp0[np] = 0;
! 281: win = W_ALLOC(np + 1); f = p; divsp(vl,p0,COEF(dc0),&f0);
! 282: for ( k = 1, win[0] = 1, --np; ; ) {
! 283: P h0,lcg,lch;
! 284: Q c;
! 285:
! 286: g0 = fp0[win[0]];
! 287: for ( i = 1; i < k; i++ ) {
! 288: mulp(vl,g0,fp0[win[i]],&m); g0 = m;
! 289: }
! 290: mulq((Q)LC(g0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lcg);
! 291: divsp(nvl,f0,g0,&h0);
! 292: mulq((Q)LC(h0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lch);
! 293: mfctrhen2(nvl,vn1,f,f0,g0,h0,lcg,lch,&g);
! 294: if ( g ) {
! 295: *tl++ = g; divsp(vl,f,g,&t);
! 296: f = t; divsp(vl,f0,g0,&t); ptozp(t,1,(Q *)&s,&f0);
! 297: for ( i = 0; i < k - 1; i++ )
! 298: for ( j = win[i] + 1; j < win[i + 1]; j++ )
! 299: fp0[j - i - 1] = fp0[j];
! 300: for ( j = win[k - 1] + 1; j <= np; j++ )
! 301: fp0[j - k] = fp0[j];
! 302: if ( ( np -= k ) < k ) break;
! 303: if ( np - win[0] + 1 < k )
! 304: if ( ++k <= np ) {
! 305: for ( i = 0; i < k; i++ ) win[i] = i + 1;
! 306: continue;
! 307: } else
! 308: break;
! 309: else for ( i = 1; i < k; i++ ) win[i] = win[0] + i;
! 310: } else {
! 311: if ( ncombi(1,np,k,win) == 0 )
! 312: if ( k == np ) break;
! 313: else for ( i = 0, ++k; i < k; i++ ) win[i] = i + 1;
! 314: }
! 315: }
! 316: *tl++ = f; *tl = 0;
! 317: for ( dc0 = 0; *l; l++ ) {
! 318: NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = *l;
! 319: }
! 320: NEXT(dc) = 0; *dcp = dc0;
! 321: }
! 322:
! 323: void ufctrmain(p,hint,dcp)
! 324: P p;
! 325: int hint;
! 326: DCP *dcp;
! 327: {
! 328: ML list;
! 329: DCP dc;
! 330:
! 331: if ( NUM(p) || (UDEG(p) == 1) ) {
! 332: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 333: *dcp = dc;
! 334: } else if ( iscycm(p) )
! 335: cycm(VR(p),UDEG(p),dcp);
! 336: else if ( iscycp(p) )
! 337: cycp(VR(p),UDEG(p),dcp);
! 338: else {
! 339: hensel(5,5,p,&list);
! 340: if ( list->n == 1 ) {
! 341: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 342: *dcp = dc;
! 343: } else
! 344: dtest(p,list,hint,dcp);
! 345: }
! 346: }
! 347:
! 348: struct oMF {
! 349: int m;
! 350: P f;
! 351: };
! 352:
! 353: void calcphi();
! 354:
! 355: void cycm(v,n,dcp)
! 356: V v;
! 357: register int n;
! 358: DCP *dcp;
! 359: {
! 360: register int i,j;
! 361: struct oMF *mfp;
! 362: DCP dc,dc0;
! 363:
! 364: if ( n == 1 ) {
! 365: NEWDC(dc); mkssum(v,1,1,-1,&COEF(dc)); DEG(dc) = ONE;
! 366: } else {
! 367: mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
! 368: bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
! 369: for ( i = 1, j = 0; i <= n; i++ )
! 370: if ( !(n%i) ) mfp[j++].m = i;
! 371: mkssum(v,1,1,-1,&mfp[0].f);
! 372: for ( i = 1; i < j; i++ )
! 373: calcphi(v,i,mfp);
! 374: for ( dc0 = 0, i = 0; i < j; i++ ) {
! 375: NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
! 376: }
! 377: }
! 378: NEXT(dc) = 0; *dcp = dc0;
! 379: }
! 380:
! 381: void cycp(v,n,dcp)
! 382: V v;
! 383: register int n;
! 384: DCP *dcp;
! 385: {
! 386: register int i,j;
! 387: int n0;
! 388: struct oMF *mfp;
! 389: DCP dc,dc0;
! 390:
! 391: if ( n == 1 ) {
! 392: NEWDC(dc); mkssum(v,1,1,1,&COEF(dc)); DEG(dc) = ONE;
! 393: } else {
! 394: n0 = n; n *= 2;
! 395: mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
! 396: bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
! 397: for ( i = 1, j = 0; i <= n; i++ )
! 398: if ( !(n%i) ) mfp[j++].m = i;
! 399: mkssum(v,1,1,-1,&mfp[0].f);
! 400: for ( i = 1; i < j; i++ )
! 401: calcphi(v,i,mfp);
! 402: for ( dc0 = 0, i = 0; i < j; i++ )
! 403: if ( n0 % mfp[i].m ) {
! 404: NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
! 405: }
! 406: }
! 407: NEXT(dc) = 0; *dcp = dc0;
! 408: }
! 409:
! 410: void calcphi(v,n,mfp)
! 411: V v;
! 412: int n;
! 413: register struct oMF *mfp;
! 414: {
! 415: register int i,m;
! 416: P t,s,tmp;
! 417:
! 418: for ( t = (P)ONE, i = 0, m = mfp[n].m; i < n; i++ )
! 419: if ( !(m%mfp[i].m) ) {
! 420: mulp(CO,t,mfp[i].f,&s); t = s;
! 421: }
! 422: mkssum(v,m,1,-1,&s); udivpz(s,t,&mfp[n].f,&tmp);
! 423: if ( tmp )
! 424: error("calcphi: cannot happen");
! 425: }
! 426:
! 427: void mkssum(v,e,s,sgn,r)
! 428: V v;
! 429: int e,s,sgn;
! 430: P *r;
! 431: {
! 432: register int i,sgnt;
! 433: DCP dc,dc0;
! 434: Q q;
! 435:
! 436: for ( dc0 = 0, i = s, sgnt = 1; i >= 0; i--, sgnt *= sgn ) {
! 437: if ( !dc0 ) {
! 438: NEWDC(dc0); dc = dc0;
! 439: } else {
! 440: NEWDC(NEXT(dc)); dc = NEXT(dc);
! 441: }
! 442: STOQ(i*e,DEG(dc)); STOQ(sgnt,q),COEF(dc) = (P)q;
! 443: }
! 444: NEXT(dc) = 0; MKP(v,dc0,*r);
! 445: }
! 446:
! 447: int iscycp(f)
! 448: P f;
! 449: {
! 450: DCP dc;
! 451: dc = DC(f);
! 452:
! 453: if ( !UNIQ((Q)COEF(dc)) )
! 454: return ( 0 );
! 455: dc = NEXT(dc);
! 456: if ( NEXT(dc) || DEG(dc) || !UNIQ((Q)COEF(dc)) )
! 457: return ( 0 );
! 458: return ( 1 );
! 459: }
! 460:
! 461: int iscycm(f)
! 462: P f;
! 463: {
! 464: DCP dc;
! 465:
! 466: dc = DC(f);
! 467: if ( !UNIQ((Q)COEF(dc)) )
! 468: return ( 0 );
! 469: dc = NEXT(dc);
! 470: if ( NEXT(dc) || DEG(dc) || !MUNIQ((Q)COEF(dc)) )
! 471: return ( 0 );
! 472: return ( 1 );
! 473: }
! 474:
! 475: void sortfs(dcp)
! 476: DCP *dcp;
! 477: {
! 478: int i,k,n,k0,d;
! 479: DCP dc,dct,t;
! 480: DCP *a;
! 481:
! 482: dc = *dcp;
! 483: for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
! 484: a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
! 485: for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ )
! 486: a[i] = dct;
! 487: a[n] = 0;
! 488:
! 489: for ( i = 0; i < n; i++ ) {
! 490: for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ )
! 491: if ( (int)UDEG(COEF(a[k])) < d ) {
! 492: k0 = k;
! 493: d = UDEG(COEF(a[k]));
! 494: }
! 495: if ( i != k0 ) {
! 496: t = a[i]; a[i] = a[k0]; a[k0] = t;
! 497: }
! 498: }
! 499: for ( *dcp = dct = a[0], i = 1; i < n; i++ )
! 500: dct = NEXT(dct) = a[i];
! 501: NEXT(dct) = 0;
! 502: }
! 503:
! 504: void sortfsrev(dcp)
! 505: DCP *dcp;
! 506: {
! 507: int i,k,n,k0,d;
! 508: DCP dc,dct,t;
! 509: DCP *a;
! 510:
! 511: dc = *dcp;
! 512: for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
! 513: a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
! 514: for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ )
! 515: a[i] = dct;
! 516: a[n] = 0;
! 517:
! 518: for ( i = 0; i < n; i++ ) {
! 519: for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ )
! 520: if ( (int)UDEG(COEF(a[k])) > d ) {
! 521: k0 = k;
! 522: d = UDEG(COEF(a[k]));
! 523: }
! 524: if ( i != k0 ) {
! 525: t = a[i]; a[i] = a[k0]; a[k0] = t;
! 526: }
! 527: }
! 528: for ( *dcp = dct = a[0], i = 1; i < n; i++ )
! 529: dct = NEXT(dct) = a[i];
! 530: NEXT(dct) = 0;
! 531: }
! 532:
! 533: void nthrootchk(f,dc,fp,dcp)
! 534: P f;
! 535: struct oDUM *dc;
! 536: ML fp;
! 537: DCP *dcp;
! 538: {
! 539: register int i,k;
! 540: int e,n,dr,tmp,t;
! 541: int *tmpp,**tmppp;
! 542: int **pp,**wpp;
! 543: LUM fpa,tpa,f0l;
! 544: UM wt,wq,ws,dvr,f0,fe;
! 545: N lc;
! 546: int lcm;
! 547: int m,b;
! 548:
! 549: m = fp->mod; b = fp->bound; fe = *((UM *)fp->c);
! 550: e = dc->n; f0 = dc->f; nthrootn(NM((Q)COEF(DC(f))),e,&lc);
! 551: if ( !lc ) {
! 552: *dcp = 0;
! 553: return;
! 554: }
! 555: lcm = rem(lc,m); W_LUMALLOC(DEG(f0),b,f0l);
! 556: for ( i = DEG(f0), tmppp = COEF(f0l), tmpp = COEF(f0);
! 557: i >= 0; i-- )
! 558: *(tmppp[i]) = dmb(m,tmpp[i],lcm,&tmp);
! 559: dtestroot(m,1,f,f0l,dc,dcp);
! 560: if ( *dcp )
! 561: return;
! 562: n = UDEG(f); W_LUMALLOC(n,b,fpa); W_LUMALLOC(n,b,tpa);
! 563: ptolum(m,b,f,fpa);
! 564: dvr = W_UMALLOC(n); wq = W_UMALLOC(n);
! 565: wt = W_UMALLOC(n); ws = W_UMALLOC(n);
! 566: cpyum(fe,dvr); divum(m,dvr,f0,wq);
! 567: t = dmb(m,pwrm(m,lcm,e-1),e,&tmp); DEG(dvr) = DEG(wq);
! 568: for ( k = 0; k <= DEG(wq); k++ )
! 569: COEF(dvr)[k] = dmb(m,COEF(wq)[k],t,&tmp);
! 570: for ( i = 1; i < b; i++ ) {
! 571: pwrlum(m,i+1,f0l,e,tpa);
! 572: for ( k = 0, pp = COEF(fpa), wpp = COEF(tpa);
! 573: k <= n; k++ )
! 574: COEF(wt)[k] = (pp[k][i]-wpp[k][i]+m)%m;
! 575: degum(wt,n); dr = divum(m,wt,dvr,ws);
! 576: if ( dr >= 0 ) {
! 577: *dcp = 0;
! 578: return;
! 579: }
! 580: for ( k = 0, pp = COEF(f0l); k <= DEG(ws); k++ )
! 581: pp[k][i] = COEF(ws)[k];
! 582: dtestroot(m,i+1,f,f0l,dc,dcp);
! 583: if ( *dcp )
! 584: return;
! 585: }
! 586: }
! 587:
! 588: void sqfrp(vl,f,dcp)
! 589: VL vl;
! 590: P f;
! 591: DCP *dcp;
! 592: {
! 593: P c,p;
! 594: DCP dc,dc0;
! 595:
! 596: if ( !f || NUM(f) ) {
! 597: NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
! 598: COEF(dc0) = f; NEXT(dc0) = 0;
! 599: } else if ( !qpcheck((Obj)f) )
! 600: error("sqfrp : invalid argument");
! 601: else {
! 602: NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
! 603: ptozp(f,1,(Q *)&c,&p); msqfr(vl,p,&dc);
! 604: COEF(dc0) = c; NEXT(dc0) = dc;
! 605: adjsgn(f,dc0);
! 606: }
! 607: }
! 608:
! 609: /*
! 610: * f : must be a poly with int coef, ignore int content
! 611: */
! 612: void msqfr(vl,f,dcp)
! 613: VL vl;
! 614: P f;
! 615: DCP *dcp;
! 616: {
! 617: DCP dc,dct,dcs;
! 618: P c,p,t,s,pc;
! 619: VL mvl;
! 620:
! 621: ptozp(f,1,(Q *)&c,&t); monomialfctr(vl,t,&p,&dc);
! 622: if ( NUM(p) ) {
! 623: *dcp = dc;
! 624: return;
! 625: }
! 626: mindegp(vl,p,&mvl,&s);
! 627: #if 0
! 628: minlcdegp(vl,p,&mvl,&s);
! 629: min_common_vars_in_coefp(vl,p,&mvl,&s);
! 630: #endif
! 631: pcp(mvl,s,&p,&pc);
! 632: if ( !NUM(pc) ) {
! 633: msqfr(mvl,pc,&dcs);
! 634: if ( !dc )
! 635: dc = dcs;
! 636: else {
! 637: for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
! 638: NEXT(dct) = dcs;
! 639: }
! 640: }
! 641: msqfrmain(mvl,p,&dcs);
! 642: for ( dct = dcs; dct; dct = NEXT(dct) ) {
! 643: reorderp(vl,mvl,COEF(dct),&t); COEF(dct) = t;
! 644: }
! 645: if ( !dc )
! 646: *dcp = dcs;
! 647: else {
! 648: for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
! 649: NEXT(dct) = dcs; *dcp = dc;
! 650: }
! 651: }
! 652:
! 653: void usqp(f,dcp)
! 654: P f;
! 655: DCP *dcp;
! 656: {
! 657: int index,nindex;
! 658: P g,c,h;
! 659: DCP dc;
! 660:
! 661: ptozp(f,1,(Q *)&c,&h);
! 662: if ( SGN((Q)LC(h)) < 0 )
! 663: chsgnp(h,&g);
! 664: else
! 665: g = h;
! 666: for ( index = 0, dc = 0; !dc; index = nindex )
! 667: hsq(index,5,g,&nindex,&dc);
! 668: *dcp = dc;
! 669: }
! 670:
! 671: void msqfrmain(vl,p,dcp)
! 672: VL vl;
! 673: P p;
! 674: DCP *dcp;
! 675: {
! 676: int i,j;
! 677: VL nvl,tvl;
! 678: VN vn,vnt,vn1;
! 679: P gg,tmp,p0,g;
! 680: DCP dc,dc0,dcr,dcr0;
! 681: int nv,d,d1;
! 682: int found;
! 683: VN svn1;
! 684: P sp0;
! 685: DCP sdc0;
! 686:
! 687: if ( deg(VR(p),p) == 1 ) {
! 688: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 689: *dcp = dc;
! 690: return;
! 691: }
! 692: clctv(vl,p,&nvl);
! 693: for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
! 694: if ( nv == 1 ) {
! 695: usqp(p,dcp);
! 696: return;
! 697: }
! 698: #if 0
! 699: if ( heusqfrchk(nvl,p) ) {
! 700: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 701: *dcp = dc;
! 702: return;
! 703: }
! 704: #endif
! 705: W_CALLOC(nv,struct oVN,vn);
! 706: W_CALLOC(nv,struct oVN,vnt);
! 707: W_CALLOC(nv,struct oVN,vn1);
! 708: W_CALLOC(nv,struct oVN,svn1);
! 709: for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
! 710: vn1[i].v = vn[i].v = tvl->v;
! 711: vn1[i].v = vn[i].v = 0;
! 712:
! 713: for ( dcr0 = 0, g = p, d = deg(VR(g),g), found = 0; ; ) {
! 714: while ( 1 ) {
! 715: for ( i = 0, j = 0; vn[i].v; i++ )
! 716: if ( vn[i].n ) vnt[j++].v = (V)i;
! 717: vnt[j].n = 0;
! 718:
! 719: mulsgn(vn,vnt,j,vn1);
! 720: substvp(nvl,LC(g),vn1,&tmp);
! 721: if ( tmp ) {
! 722: substvp(nvl,g,vn1,&p0);
! 723: usqp(p0,&dc0);
! 724: for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )
! 725: if ( DEG(dc) )
! 726: d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));
! 727:
! 728: if ( d1 == 0 ) {
! 729: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = g; NEXT(dc) = 0;
! 730: if ( !dcr0 )
! 731: dcr0 = dc;
! 732: else
! 733: NEXT(dcr) = dc;
! 734: *dcp = dcr0;
! 735: return;
! 736: }
! 737:
! 738: if ( d < d1 )
! 739: goto END;
! 740: if ( d > d1 ) {
! 741: d = d1;
! 742: found = 1;
! 743: bcopy((char *)vn1,(char *)svn1,(int)(sizeof(struct oVN)*nv));
! 744: sp0 = p0; sdc0 = dc0;
! 745: goto END;
! 746: }
! 747: /* d1 == d */
! 748: if ( found ) {
! 749: found = 0;
! 750: #if 0
! 751: if ( d > d1 ) {
! 752: d = d1;
! 753: /*} */
! 754: #endif
! 755: msqfrmainmain(nvl,g,svn1,sp0,sdc0,&dc,&gg); g = gg;
! 756: if ( dc ) {
! 757: if ( !dcr0 )
! 758: dcr0 = dc;
! 759: else
! 760: NEXT(dcr) = dc;
! 761: for ( dcr = dc; NEXT(dcr); dcr = NEXT(dcr) );
! 762: if ( NUM(g) ) {
! 763: NEXT(dcr) = 0; *dcp = dcr0;
! 764: return;
! 765: }
! 766: d = deg(VR(g),g);
! 767: }
! 768: }
! 769: }
! 770: END:
! 771: if ( nextbin(vnt,j) )
! 772: break;
! 773: }
! 774: next(vn);
! 775: }
! 776: }
! 777:
! 778: void msqfrmainmain(vl,p,vn,p0,dc0,dcp,pp)
! 779: VL vl;
! 780: P p;
! 781: VN vn;
! 782: P p0;
! 783: DCP dc0;
! 784: DCP *dcp;
! 785: P *pp;
! 786: {
! 787: int i,j,k,np;
! 788: DCP *a;
! 789: DCP dc,dcr,dcr0,dct;
! 790: P g,t,s,u,t0,f,f0,d,d0,g0,h0,x,xx;
! 791: Q q;
! 792: V v;
! 793:
! 794: for ( np = 0, dc = dc0; dc; dc = NEXT(dc), np++ );
! 795: a = (DCP *)ALLOCA((np + 1)*sizeof(DCP));
! 796: for ( i = 0, dc = dc0; dc; i++, dc = NEXT(dc) )
! 797: a[np-i-1] = dc;
! 798:
! 799: for ( i = 0, dcr0 = 0, f = p, f0 = p0, v = VR(p);
! 800: i < np; i++ ) {
! 801: if ( (i == (np-1))&&UNIQ(DEG(a[i])) ) {
! 802: NEXTDC(dcr0,dcr);
! 803: DEG(dcr) = DEG(a[i]);
! 804: COEF(dcr) = f;
! 805: f = (P)ONE;
! 806: } else if ( (i == (np-1))&&UNIQ(DEG(DC(COEF(a[i])))) ) {
! 807: diffp(vl,f,v,&s); pcp(vl,s,&t,&u);
! 808: if ( divtpz(vl,f,t,&s) ) {
! 809: NEXTDC(dcr0,dcr);
! 810: DEG(dcr) = DEG(a[i]);
! 811: COEF(dcr) = s;
! 812: f = (P)ONE;
! 813: } else
! 814: break;
! 815: } else {
! 816: for ( t = f, t0 = f0,
! 817: j = 0, k = QTOS(DEG(a[i]))-1; j < k; j++ ) {
! 818: diffp(vl,t,v,&s); t = s;
! 819: diffp(vl,t0,v,&s); t0 = s;
! 820: }
! 821: factorial(k,&q);
! 822: divsp(vl,t,(P)q,&s);
! 823: for ( dct = DC(s); NEXT(dct); dct = NEXT(dct) );
! 824: if ( DEG(dct) ) {
! 825: MKV(VR(s),x); pwrp(vl,x,DEG(dct),&xx);
! 826: divsp(vl,s,xx,&d);
! 827: } else {
! 828: xx = (P)ONE; d = s;
! 829: }
! 830: divsp(vl,t0,xx,&t);
! 831: divsp(vl,t,(P)q,&s);
! 832: ptozp(s,1,(Q *)&t,&d0);
! 833:
! 834: for ( dct = DC(COEF(a[i])); NEXT(dct); dct = NEXT(dct) );
! 835: if ( DEG(dct) )
! 836: divsp(vl,COEF(a[i]),xx,&g0);
! 837: else {
! 838: xx = (P)ONE; g0 = COEF(a[i]);
! 839: }
! 840:
! 841: pcp(vl,d,&t,&u); d = t;
! 842: ptozp(g0,1,(Q *)&u,&t); g0 = t;
! 843:
! 844: {
! 845: DCP dca,dcb;
! 846:
! 847: fctrp(vl,LC(d),&dca);
! 848: for ( dcb = dca, u = (P)ONE; dcb; dcb = NEXT(dcb) ) {
! 849: if ( NUM(COEF(dcb)) ) {
! 850: mulp(vl,u,COEF(dcb),&t); u = t;
! 851: } else {
! 852: Q qq;
! 853: P tt;
! 854:
! 855: pwrp(vl,COEF(dcb),DEG(a[i]),&s);
! 856: for ( t = LC(f), j = 0; divtpz(vl,t,s,&tt); j++, t = tt );
! 857: STOQ(j,qq);
! 858: if ( cmpq(qq,DEG(dcb)) > 0 )
! 859: qq = DEG(dcb);
! 860: pwrp(vl,COEF(dcb),qq,&t); mulp(vl,u,t,&s); u = s;
! 861: }
! 862: }
! 863: divsp(vl,d0,g0,&h0);
! 864: }
! 865: mfctrhen2(vl,vn,d,d0,g0,h0,u,LC(d),&s);
! 866: if ( s ) {
! 867: mulp(vl,s,xx,&g);
! 868: pwrp(vl,g,DEG(a[i]),&t);
! 869: if ( divtpz(vl,f,t,&s) ) {
! 870: NEXTDC(dcr0,dcr);
! 871: DEG(dcr) = DEG(a[i]); COEF(dcr) = g;
! 872: f = s; substvp(vl,f,vn,&f0);
! 873: } else
! 874: break;
! 875: } else
! 876: break;
! 877: }
! 878: }
! 879: *pp = f;
! 880: if ( dcr0 )
! 881: NEXT(dcr) = 0;
! 882: *dcp = dcr0;
! 883: }
! 884:
! 885: void mfctrhen2(vl,vn,f,f0,g0,h0,lcg,lch,gp)
! 886: VL vl;
! 887: VN vn;
! 888: P f;
! 889: P f0,g0,h0,lcg,lch;
! 890: P *gp;
! 891: {
! 892: V v;
! 893: P f1,lc,lc0,lcg0,lch0;
! 894: P m,ff0,gg0,hh0,gk,hk,ggg,gggr,hhh,ak,bk,tmp;
! 895: Q bb,qk,s;
! 896: Q cbd;
! 897: int dbd;
! 898: int d;
! 899:
! 900: if ( NUM(g0) ) {
! 901: *gp = (P)ONE;
! 902: return;
! 903: }
! 904:
! 905: v = VR(f); d = deg(v,f);
! 906: if ( d == deg(v,g0) ) {
! 907: pcp(vl,f,gp,&tmp);
! 908: return;
! 909: }
! 910:
! 911: mulp(vl,lcg,lch,&lc);
! 912: if ( !divtpz(vl,lc,LC(f),(P *)&s) ) {
! 913: *gp = 0;
! 914: return;
! 915: }
! 916: mulp(vl,(P)s,f,&f1);
! 917: dbd = dbound(VR(f1),f1) + 1; cbound(vl,f1,&cbd);
! 918:
! 919: substvp(vl,lc,vn,&lc0);
! 920: divq((Q)lc0,(Q)LC(f0),(Q *)&m); mulp(vl,f0,m,&ff0);
! 921: substvp(vl,lcg,vn,&lcg0);
! 922: divq((Q)lcg0,(Q)LC(g0),(Q *)&m); mulp(vl,g0,m,&gg0);
! 923: substvp(vl,lch,vn,&lch0);
! 924: divq((Q)lch0,(Q)LC(h0),(Q *)&m); mulp(vl,h0,m,&hh0);
! 925: addq(cbd,cbd,&bb);
! 926: henzq1(gg0,hh0,bb,&bk,&ak,&qk); gk = gg0; hk = hh0;
! 927: henmv(vl,vn,f1,gk,hk,ak,bk,
! 928: lcg,lch,lcg0,lch0,qk,dbd,&ggg,&hhh);
! 929:
! 930: if ( divtpz(vl,f1,ggg,&gggr) )
! 931: pcp(vl,ggg,gp,&tmp);
! 932: else
! 933: *gp = 0;
! 934: }
! 935:
! 936: int sqfrchk(p)
! 937: P p;
! 938: {
! 939: Q c;
! 940: P f;
! 941: DCP dc;
! 942:
! 943: ptozp(p,SGN((Q)UCOEF(p)),&c,&f); usqp(f,&dc);
! 944: if ( NEXT(dc) || !UNIQ(DEG(dc)) )
! 945: return ( 0 );
! 946: else
! 947: return ( 1 );
! 948: }
! 949:
! 950: int cycchk(p)
! 951: P p;
! 952: {
! 953: Q c;
! 954: P f;
! 955:
! 956: ptozp(p,SGN((Q)UCOEF(p)),&c,&f);
! 957: if ( iscycp(f) || iscycm(f) )
! 958: return 0;
! 959: else
! 960: return 1;
! 961: }
! 962:
! 963: int zerovpchk(vl,p,vn)
! 964: VL vl;
! 965: P p;
! 966: VN vn;
! 967: {
! 968: P t;
! 969:
! 970: substvp(vl,p,vn,&t);
! 971: if ( t )
! 972: return ( 0 );
! 973: else
! 974: return ( 1 );
! 975: }
! 976:
! 977: int valideval(vl,dc,vn)
! 978: VL vl;
! 979: DCP dc;
! 980: VN vn;
! 981: {
! 982: DCP dct;
! 983: Q *a;
! 984: int i,j,n;
! 985: N q,r;
! 986:
! 987: for ( dct = NEXT(dc), n = 0; dct; dct = NEXT(dct), n++ );
! 988: W_CALLOC(n,Q,a);
! 989: for ( dct = NEXT(dc), i = 0; dct; dct = NEXT(dct), i++ ) {
! 990: substvp(vl,COEF(dct),vn,(P *)&a[i]);
! 991: if ( !a[i] )
! 992: return ( 0 );
! 993:
! 994: for ( j = 0; j < i; j++ ) {
! 995: divn(NM(a[j]),NM(a[i]),&q,&r);
! 996: if ( !r )
! 997: return ( 0 );
! 998: divn(NM(a[i]),NM(a[j]),&q,&r);
! 999: if ( !r )
! 1000: return ( 0 );
! 1001: }
! 1002: }
! 1003: return ( 1 );
! 1004: }
! 1005:
! 1006: void estimatelc(vl,c,dc,vn,lcp)
! 1007: VL vl;
! 1008: Q c;
! 1009: DCP dc;
! 1010: VN vn;
! 1011: P *lcp;
! 1012: {
! 1013: int i;
! 1014: DCP dct;
! 1015: P r,s,t;
! 1016: Q c0,c1,c2;
! 1017:
! 1018: for ( dct = dc, r = (P)ONE; dct; dct = NEXT(dct) ) {
! 1019: if ( NUM(COEF(dct)) ) {
! 1020: mulp(vl,r,COEF(dct),&s); r = s;
! 1021: } else {
! 1022: substvp(vl,COEF(dct),vn,(P *)&c0);
! 1023: for ( i = 0, c1 = c; i < (int)QTOS(DEG(dct)); i++ ) {
! 1024: divq(c1,c0,&c2);
! 1025: if ( !INT(c2) )
! 1026: break;
! 1027: else
! 1028: c1 = c2;
! 1029: }
! 1030: if ( i ) {
! 1031: STOQ(i,c1);
! 1032: pwrp(vl,COEF(dct),c1,&s); mulp(vl,r,s,&t); r = t;
! 1033: }
! 1034: }
! 1035: }
! 1036: *lcp = r;
! 1037: }
! 1038:
! 1039: void monomialfctr(vl,p,pr,dcp)
! 1040: VL vl;
! 1041: P p;
! 1042: P *pr;
! 1043: DCP *dcp;
! 1044: {
! 1045: VL nvl,avl;
! 1046: Q d;
! 1047: P f,t,s;
! 1048: DCP dc0,dc;
! 1049:
! 1050: clctv(vl,p,&nvl);
! 1051: for ( dc0 = 0, avl = nvl, f = p; avl; avl = NEXT(avl) ) {
! 1052: getmindeg(avl->v,f,&d);
! 1053: if ( d ) {
! 1054: MKV(avl->v,t); NEXTDC(dc0,dc); DEG(dc) = d; COEF(dc) = t;
! 1055: pwrp(vl,t,d,&s); divsp(vl,f,s,&t); f = t;
! 1056: }
! 1057: }
! 1058: if ( dc0 )
! 1059: NEXT(dc) = 0;
! 1060: *pr = f; *dcp = dc0;
! 1061: }
! 1062:
! 1063: void afctr(vl,p0,p,dcp)
! 1064: VL vl;
! 1065: P p,p0;
! 1066: DCP *dcp;
! 1067: {
! 1068: DCP dc,dc0,dcr,dct,dcs;
! 1069: P t;
! 1070: VL nvl;
! 1071:
! 1072: if ( VR(p) == VR(p0) ) {
! 1073: NEWDC(dc);
! 1074: DEG(dc) = ONE;
! 1075: COEF(dc) = p;
! 1076: NEXT(dc) = 0;
! 1077: *dcp = dc;
! 1078: return;
! 1079: }
! 1080:
! 1081: clctv(vl,p,&nvl);
! 1082: if ( !NEXT(nvl) )
! 1083: ufctr(p,1,&dc);
! 1084: else {
! 1085: sqa(vl,p0,p,&dc);
! 1086: for ( dct = dc; dct; dct = NEXT(dct) ) {
! 1087: pmonic(vl,p0,COEF(dct),&t); COEF(dct) = t;
! 1088: }
! 1089: }
! 1090: if ( NUM(COEF(dc)) )
! 1091: dcr = NEXT(dc);
! 1092: else
! 1093: dcr = dc;
! 1094: for ( dc0 = 0; dcr; dcr = NEXT(dcr) ) {
! 1095: afctrmain(vl,p0,COEF(dcr),1,&dcs);
! 1096:
! 1097: for ( dct = dcs; dct; dct = NEXT(dct) )
! 1098: DEG(dct) = DEG(dcr);
! 1099: if ( !dc0 )
! 1100: dc0 = dcs;
! 1101: else {
! 1102: for ( dct = dc0; NEXT(dct); dct = NEXT(dct) );
! 1103: NEXT(dct) = dcs;
! 1104: }
! 1105: }
! 1106: *dcp = dc0;
! 1107: }
! 1108:
! 1109: void afctrmain(vl,p0,p,init,dcp)
! 1110: VL vl;
! 1111: P p,p0;
! 1112: int init;
! 1113: DCP *dcp;
! 1114: {
! 1115: P x,y,s,m,a,t,u,pt,pt1,res,g;
! 1116: Q q;
! 1117: DCP dc,dc0,dcsq0,dcr0,dcr,dct,dcs;
! 1118: V v,v0;
! 1119:
! 1120: if ( !cmpq(DEG(DC(p)),ONE) ) {
! 1121: NEWDC(dc); DEG(dc) = ONE; NEXT(dc) = 0;
! 1122: pmonic(vl,p0,p,&COEF(dc)); *dcp = dc;
! 1123: return;
! 1124: }
! 1125:
! 1126: v = VR(p); MKV(v,x);
! 1127: v0 = VR(p0); MKV(v0,y);
! 1128: STOQ(init,q),s = (P)q;
! 1129: mulp(vl,s,y,&m); subp(vl,x,m,&t); addp(vl,x,m,&a);
! 1130: substp(vl,p,v,t,&pt);
! 1131: remsdcp(vl,pt,p0,&pt1);
! 1132:
! 1133: /*
! 1134: if ( ( deg(v0,p0) <= 3 ) || ( TCPQ(p0) <= 2 ) )
! 1135: resultp(vl,v0,p0,pt1,&res);
! 1136: else
! 1137: srcrnorm(vl,v0,pt1,p0,&res);
! 1138: */
! 1139: #if 0
! 1140: srcr(vl,v0,pt1,p0,&res);
! 1141: #endif
! 1142: resultp(vl,v0,p0,pt1,&res);
! 1143: usqp(res,&dcsq0);
! 1144: for ( dc0 = 0, dct = dcsq0; dct; dct = NEXT(dct) ) {
! 1145: if ( UNIQ(DEG(dct)) )
! 1146: ufctr(COEF(dct),deg(v0,p0),&dcs);
! 1147: else
! 1148: ufctr(COEF(dct),1,&dcs);
! 1149: for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
! 1150: DEG(dcr) = DEG(dct);
! 1151: if ( !dc0 ) {
! 1152: dc0 = NEXT(dcs);
! 1153: dc = dc0;
! 1154: } else {
! 1155: for ( ; NEXT(dc); dc = NEXT(dc) );
! 1156: NEXT(dc) = NEXT(dcs);
! 1157: }
! 1158: }
! 1159: sortfs(&dc0);
! 1160:
! 1161: for ( g = p, dcr = dcr0 = 0, dc = dc0; dc; dc = NEXT(dc) ) {
! 1162: if ( !UNIQ(DEG(dc)) ) {
! 1163: substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s);
! 1164: gcda(vl,p0,s,g,&u);
! 1165: if ( !NUM(u) && (VR(u) == v)) {
! 1166: afctrmain(vl,p0,u,init+1,&dct);
! 1167: for ( dcs = dct, t = (P)ONE; dcs; dcs = NEXT(dcs) ) {
! 1168: mulp(vl,t,COEF(dcs),&s); remsdcp(vl,s,p0,&t);
! 1169: }
! 1170: pdiva(vl,p0,g,t,&s); g = s;
! 1171: if ( !dcr0 )
! 1172: dcr0 = dct;
! 1173: else
! 1174: NEXT(dcr) = dct;
! 1175: for ( dcr = dct; NEXT(dcr); dcr = NEXT(dcr) );
! 1176: }
! 1177: } else {
! 1178: substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s);
! 1179: gcda(vl,p0,s,g,&u);
! 1180: if ( !NUM(u) && (VR(u) == v)) {
! 1181: NEXTDC(dcr0,dcr);
! 1182: DEG(dcr) = ONE;
! 1183: COEF(dcr) = u;
! 1184: pdiva(vl,p0,g,u,&t); g = t;
! 1185: }
! 1186: }
! 1187: }
! 1188: if ( dcr0 )
! 1189: NEXT(dcr) = 0;
! 1190: *dcp = dcr0;
! 1191: }