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