![[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
1.4 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.3 noro 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: *
1.5 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/F.c,v 1.4 2000/08/22 05:04:04 noro Exp $
1.3 noro 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: {
1.5 ! noro 276: int i,j,k,*win,np,x;
1.1 noro 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;
1.5 ! noro 283: int *nonzero;
1.1 noro 284:
285: if ( !cmpq(DEG(DC(p)),ONE) ) {
286: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
287: *dcp = dc; return;
288: }
289: lp = LC(p); fctrp(vl,lp,&dcl); clctv(vl,p,&nvl);
290: for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
291: W_CALLOC(nv,struct oVN,vn); W_CALLOC(nv,struct oVN,vnt);
292: W_CALLOC(nv,struct oVN,vn1);
1.5 ! noro 293: W_CALLOC(nv,int,nonzero);
! 294:
1.1 noro 295: for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
296: vn1[i].v = vn[i].v = tvl->v;
297: vn1[i].v = vn[i].v = 0;
1.5 ! noro 298:
! 299: /* determine a heuristic bound of deg(GCD(p,p')) */
! 300: while ( 1 ) {
! 301: for ( i = 0; vn1[i].v; i++ )
! 302: vn1[i].n = ((unsigned int)random())%256+1;
! 303: substvp(nvl,LC(p),vn1,&p0);
! 304: if ( p0 ) {
! 305: substvp(nvl,p,vn1,&p0);
! 306: if ( sqfrchk(p0) ) {
! 307: ufctr(p0,1,&dc0);
! 308: if ( NEXT(NEXT(dc0)) == 0 ) {
! 309: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 310: *dcp = dc;
! 311: return;
! 312: } else {
! 313: for ( dc0 = NEXT(dc0), np = 0; dc0; dc0 = NEXT(dc0), np++ );
! 314: break;
! 315: }
! 316: }
! 317: }
! 318: }
! 319:
! 320: for ( i = 0; vn1[i].v; i++ ) {
! 321: x = vn1[i].n; vn1[i].n = 0;
! 322: substvp(nvl,LC(p),vn1,&p0);
! 323: if ( !p0 )
! 324: vn1[i].n = x;
! 325: else {
! 326: substvp(nvl,p,vn1,&p0);
! 327: if ( !sqfrchk(p0) )
! 328: vn1[i].n = x;
! 329: else {
! 330: ufctr(p0,1,&dc0);
! 331: for ( dc0 = NEXT(dc0), j = 0; dc0; dc0 = NEXT(dc0), j++ );
! 332: if ( j > np )
! 333: vn1[i].n = x;
! 334: }
! 335: }
! 336: }
! 337: for ( i = 0; vn1[i].v; i++ )
! 338: if (vn1[i].n )
! 339: nonzero[i] = 1;
! 340:
1.1 noro 341: count = 0;
342: while ( 1 ) {
343: while ( 1 ) {
344: count++;
345: for ( i = 0, j = 0; vn[i].v; i++ )
346: if ( vn[i].n )
347: vnt[j++].v = (V)i;
348: vnt[j].n = 0; mulsgn(vn,vnt,j,vn1);
349: for ( i = 0; vn1[i].v; i++ )
350: if ( vn1[i].n ) {
351: if ( vn1[i].n > 0 )
352: vn1[i].n = sprime[vn1[i].n];
353: else
354: vn1[i].n = -sprime[-vn1[i].n];
355: }
356: if ( valideval(nvl,dcl,vn1) ) {
357: substvp(nvl,p,vn1,&p0);
358: if ( sqfrchk(p0) ) {
359: ufctr(p0,1,&dc0);
360: if ( NEXT(NEXT(dc0)) == 0 ) {
361: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
362: *dcp = dc;
363: return;
1.5 ! noro 364: } else {
! 365: for ( dc = NEXT(dc0), i = 0; dc; dc = NEXT(dc), i++ );
! 366: if ( i <= np )
! 367: goto MAIN;
! 368: if ( i < np )
! 369: np = i;
! 370: }
! 371: }
1.1 noro 372: }
373: if ( nextbin(vnt,j) )
374: break;
375: }
1.5 ! noro 376: while ( 1 ) {
! 377: next(vn);
! 378: for ( i = 0; vn[i].v; i++ )
! 379: if ( nonzero[i] && !vn[i].n )
! 380: break;
! 381: if ( !vn[i].v )
! 382: break;
! 383: }
1.1 noro 384: }
385: MAIN :
386: #if 0
387: for ( i = 0; vn1[i].v; i++ )
388: fprintf(stderr,"%d ",vn1[i].n);
389: fprintf(stderr,"\n");
390: #endif
391: for ( np = 0, dc = NEXT(dc0); dc; dc = NEXT(dc), np++ );
392: fp0 = (P *)ALLOCA((np + 1)*sizeof(P));
393: fpt = (P *)ALLOCA((np + 1)*sizeof(P));
394: l = tl = (P *)ALLOCA((np + 1)*sizeof(P));
395: for ( i = 1, dc = NEXT(dc0); dc; i++, dc = NEXT(dc) )
396: fp0[i-1] = COEF(dc);
397: #if 0
398: sort_by_deg(np,fp0,fpt);
399: sort_by_deg_rev(np-1,fpt+1,fp0+1);
400: #endif
401: #if 0
402: sort_by_deg_rev(np,fp0,fpt); fp0[0] = fpt[0];
403: sort_by_deg(np-1,fpt+1,fp0+1); fp0[np] = 0;
404: #endif
405: fp0[np] = 0;
406: win = W_ALLOC(np + 1); f = p; divsp(vl,p0,COEF(dc0),&f0);
407: for ( k = 1, win[0] = 1, --np; ; ) {
408: P h0,lcg,lch;
409: Q c;
410:
411: g0 = fp0[win[0]];
412: for ( i = 1; i < k; i++ ) {
413: mulp(vl,g0,fp0[win[i]],&m); g0 = m;
414: }
415: mulq((Q)LC(g0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lcg);
416: divsp(nvl,f0,g0,&h0);
417: mulq((Q)LC(h0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lch);
418: mfctrhen2(nvl,vn1,f,f0,g0,h0,lcg,lch,&g);
419: if ( g ) {
420: *tl++ = g; divsp(vl,f,g,&t);
421: f = t; divsp(vl,f0,g0,&t); ptozp(t,1,(Q *)&s,&f0);
422: for ( i = 0; i < k - 1; i++ )
423: for ( j = win[i] + 1; j < win[i + 1]; j++ )
424: fp0[j - i - 1] = fp0[j];
425: for ( j = win[k - 1] + 1; j <= np; j++ )
426: fp0[j - k] = fp0[j];
427: if ( ( np -= k ) < k ) break;
428: if ( np - win[0] + 1 < k )
429: if ( ++k <= np ) {
430: for ( i = 0; i < k; i++ ) win[i] = i + 1;
431: continue;
432: } else
433: break;
434: else for ( i = 1; i < k; i++ ) win[i] = win[0] + i;
435: } else {
436: if ( ncombi(1,np,k,win) == 0 )
437: if ( k == np ) break;
438: else for ( i = 0, ++k; i < k; i++ ) win[i] = i + 1;
439: }
440: }
441: *tl++ = f; *tl = 0;
442: for ( dc0 = 0; *l; l++ ) {
443: NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = *l;
444: }
445: NEXT(dc) = 0; *dcp = dc0;
446: }
447:
448: void ufctrmain(p,hint,dcp)
449: P p;
450: int hint;
451: DCP *dcp;
452: {
453: ML list;
454: DCP dc;
455:
456: if ( NUM(p) || (UDEG(p) == 1) ) {
457: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
458: *dcp = dc;
459: } else if ( iscycm(p) )
460: cycm(VR(p),UDEG(p),dcp);
461: else if ( iscycp(p) )
462: cycp(VR(p),UDEG(p),dcp);
463: else {
464: hensel(5,5,p,&list);
465: if ( list->n == 1 ) {
466: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
467: *dcp = dc;
468: } else
469: dtest(p,list,hint,dcp);
470: }
471: }
472:
473: struct oMF {
474: int m;
475: P f;
476: };
477:
478: void calcphi();
479:
480: void cycm(v,n,dcp)
481: V v;
482: register int n;
483: DCP *dcp;
484: {
485: register int i,j;
486: struct oMF *mfp;
487: DCP dc,dc0;
488:
489: if ( n == 1 ) {
490: NEWDC(dc); mkssum(v,1,1,-1,&COEF(dc)); DEG(dc) = ONE;
491: } else {
492: mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
493: bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
494: for ( i = 1, j = 0; i <= n; i++ )
495: if ( !(n%i) ) mfp[j++].m = i;
496: mkssum(v,1,1,-1,&mfp[0].f);
497: for ( i = 1; i < j; i++ )
498: calcphi(v,i,mfp);
499: for ( dc0 = 0, i = 0; i < j; i++ ) {
500: NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
501: }
502: }
503: NEXT(dc) = 0; *dcp = dc0;
504: }
505:
506: void cycp(v,n,dcp)
507: V v;
508: register int n;
509: DCP *dcp;
510: {
511: register int i,j;
512: int n0;
513: struct oMF *mfp;
514: DCP dc,dc0;
515:
516: if ( n == 1 ) {
517: NEWDC(dc); mkssum(v,1,1,1,&COEF(dc)); DEG(dc) = ONE;
518: } else {
519: n0 = n; n *= 2;
520: mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
521: bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
522: for ( i = 1, j = 0; i <= n; i++ )
523: if ( !(n%i) ) mfp[j++].m = i;
524: mkssum(v,1,1,-1,&mfp[0].f);
525: for ( i = 1; i < j; i++ )
526: calcphi(v,i,mfp);
527: for ( dc0 = 0, i = 0; i < j; i++ )
528: if ( n0 % mfp[i].m ) {
529: NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
530: }
531: }
532: NEXT(dc) = 0; *dcp = dc0;
533: }
534:
535: void calcphi(v,n,mfp)
536: V v;
537: int n;
538: register struct oMF *mfp;
539: {
540: register int i,m;
541: P t,s,tmp;
542:
543: for ( t = (P)ONE, i = 0, m = mfp[n].m; i < n; i++ )
544: if ( !(m%mfp[i].m) ) {
545: mulp(CO,t,mfp[i].f,&s); t = s;
546: }
547: mkssum(v,m,1,-1,&s); udivpz(s,t,&mfp[n].f,&tmp);
548: if ( tmp )
549: error("calcphi: cannot happen");
550: }
551:
552: void mkssum(v,e,s,sgn,r)
553: V v;
554: int e,s,sgn;
555: P *r;
556: {
557: register int i,sgnt;
558: DCP dc,dc0;
559: Q q;
560:
561: for ( dc0 = 0, i = s, sgnt = 1; i >= 0; i--, sgnt *= sgn ) {
562: if ( !dc0 ) {
563: NEWDC(dc0); dc = dc0;
564: } else {
565: NEWDC(NEXT(dc)); dc = NEXT(dc);
566: }
567: STOQ(i*e,DEG(dc)); STOQ(sgnt,q),COEF(dc) = (P)q;
568: }
569: NEXT(dc) = 0; MKP(v,dc0,*r);
570: }
571:
572: int iscycp(f)
573: P f;
574: {
575: DCP dc;
576: dc = DC(f);
577:
578: if ( !UNIQ((Q)COEF(dc)) )
579: return ( 0 );
580: dc = NEXT(dc);
581: if ( NEXT(dc) || DEG(dc) || !UNIQ((Q)COEF(dc)) )
582: return ( 0 );
583: return ( 1 );
584: }
585:
586: int iscycm(f)
587: P f;
588: {
589: DCP dc;
590:
591: dc = DC(f);
592: if ( !UNIQ((Q)COEF(dc)) )
593: return ( 0 );
594: dc = NEXT(dc);
595: if ( NEXT(dc) || DEG(dc) || !MUNIQ((Q)COEF(dc)) )
596: return ( 0 );
597: return ( 1 );
598: }
599:
600: void sortfs(dcp)
601: DCP *dcp;
602: {
603: int i,k,n,k0,d;
604: DCP dc,dct,t;
605: DCP *a;
606:
607: dc = *dcp;
608: for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
609: a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
610: for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ )
611: a[i] = dct;
612: a[n] = 0;
613:
614: for ( i = 0; i < n; i++ ) {
615: for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ )
616: if ( (int)UDEG(COEF(a[k])) < d ) {
617: k0 = k;
618: d = UDEG(COEF(a[k]));
619: }
620: if ( i != k0 ) {
621: t = a[i]; a[i] = a[k0]; a[k0] = t;
622: }
623: }
624: for ( *dcp = dct = a[0], i = 1; i < n; i++ )
625: dct = NEXT(dct) = a[i];
626: NEXT(dct) = 0;
627: }
628:
629: void sortfsrev(dcp)
630: DCP *dcp;
631: {
632: int i,k,n,k0,d;
633: DCP dc,dct,t;
634: DCP *a;
635:
636: dc = *dcp;
637: for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
638: a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
639: for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ )
640: a[i] = dct;
641: a[n] = 0;
642:
643: for ( i = 0; i < n; i++ ) {
644: for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ )
645: if ( (int)UDEG(COEF(a[k])) > d ) {
646: k0 = k;
647: d = UDEG(COEF(a[k]));
648: }
649: if ( i != k0 ) {
650: t = a[i]; a[i] = a[k0]; a[k0] = t;
651: }
652: }
653: for ( *dcp = dct = a[0], i = 1; i < n; i++ )
654: dct = NEXT(dct) = a[i];
655: NEXT(dct) = 0;
656: }
657:
658: void nthrootchk(f,dc,fp,dcp)
659: P f;
660: struct oDUM *dc;
661: ML fp;
662: DCP *dcp;
663: {
664: register int i,k;
665: int e,n,dr,tmp,t;
666: int *tmpp,**tmppp;
667: int **pp,**wpp;
668: LUM fpa,tpa,f0l;
669: UM wt,wq,ws,dvr,f0,fe;
670: N lc;
671: int lcm;
672: int m,b;
673:
674: m = fp->mod; b = fp->bound; fe = *((UM *)fp->c);
675: e = dc->n; f0 = dc->f; nthrootn(NM((Q)COEF(DC(f))),e,&lc);
676: if ( !lc ) {
677: *dcp = 0;
678: return;
679: }
680: lcm = rem(lc,m); W_LUMALLOC(DEG(f0),b,f0l);
681: for ( i = DEG(f0), tmppp = COEF(f0l), tmpp = COEF(f0);
682: i >= 0; i-- )
683: *(tmppp[i]) = dmb(m,tmpp[i],lcm,&tmp);
684: dtestroot(m,1,f,f0l,dc,dcp);
685: if ( *dcp )
686: return;
687: n = UDEG(f); W_LUMALLOC(n,b,fpa); W_LUMALLOC(n,b,tpa);
688: ptolum(m,b,f,fpa);
689: dvr = W_UMALLOC(n); wq = W_UMALLOC(n);
690: wt = W_UMALLOC(n); ws = W_UMALLOC(n);
691: cpyum(fe,dvr); divum(m,dvr,f0,wq);
692: t = dmb(m,pwrm(m,lcm,e-1),e,&tmp); DEG(dvr) = DEG(wq);
693: for ( k = 0; k <= DEG(wq); k++ )
694: COEF(dvr)[k] = dmb(m,COEF(wq)[k],t,&tmp);
695: for ( i = 1; i < b; i++ ) {
696: pwrlum(m,i+1,f0l,e,tpa);
697: for ( k = 0, pp = COEF(fpa), wpp = COEF(tpa);
698: k <= n; k++ )
699: COEF(wt)[k] = (pp[k][i]-wpp[k][i]+m)%m;
700: degum(wt,n); dr = divum(m,wt,dvr,ws);
701: if ( dr >= 0 ) {
702: *dcp = 0;
703: return;
704: }
705: for ( k = 0, pp = COEF(f0l); k <= DEG(ws); k++ )
706: pp[k][i] = COEF(ws)[k];
707: dtestroot(m,i+1,f,f0l,dc,dcp);
708: if ( *dcp )
709: return;
710: }
711: }
712:
713: void sqfrp(vl,f,dcp)
714: VL vl;
715: P f;
716: DCP *dcp;
717: {
718: P c,p;
719: DCP dc,dc0;
720:
721: if ( !f || NUM(f) ) {
722: NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
723: COEF(dc0) = f; NEXT(dc0) = 0;
724: } else if ( !qpcheck((Obj)f) )
725: error("sqfrp : invalid argument");
726: else {
727: NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
728: ptozp(f,1,(Q *)&c,&p); msqfr(vl,p,&dc);
729: COEF(dc0) = c; NEXT(dc0) = dc;
730: adjsgn(f,dc0);
731: }
732: }
733:
734: /*
735: * f : must be a poly with int coef, ignore int content
736: */
737: void msqfr(vl,f,dcp)
738: VL vl;
739: P f;
740: DCP *dcp;
741: {
742: DCP dc,dct,dcs;
743: P c,p,t,s,pc;
744: VL mvl;
745:
746: ptozp(f,1,(Q *)&c,&t); monomialfctr(vl,t,&p,&dc);
747: if ( NUM(p) ) {
748: *dcp = dc;
749: return;
750: }
751: mindegp(vl,p,&mvl,&s);
752: #if 0
753: minlcdegp(vl,p,&mvl,&s);
754: min_common_vars_in_coefp(vl,p,&mvl,&s);
755: #endif
756: pcp(mvl,s,&p,&pc);
757: if ( !NUM(pc) ) {
758: msqfr(mvl,pc,&dcs);
759: if ( !dc )
760: dc = dcs;
761: else {
762: for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
763: NEXT(dct) = dcs;
764: }
765: }
766: msqfrmain(mvl,p,&dcs);
767: for ( dct = dcs; dct; dct = NEXT(dct) ) {
768: reorderp(vl,mvl,COEF(dct),&t); COEF(dct) = t;
769: }
770: if ( !dc )
771: *dcp = dcs;
772: else {
773: for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
774: NEXT(dct) = dcs; *dcp = dc;
775: }
776: }
777:
778: void usqp(f,dcp)
779: P f;
780: DCP *dcp;
781: {
782: int index,nindex;
783: P g,c,h;
784: DCP dc;
785:
786: ptozp(f,1,(Q *)&c,&h);
787: if ( SGN((Q)LC(h)) < 0 )
788: chsgnp(h,&g);
789: else
790: g = h;
791: for ( index = 0, dc = 0; !dc; index = nindex )
792: hsq(index,5,g,&nindex,&dc);
793: *dcp = dc;
794: }
795:
796: void msqfrmain(vl,p,dcp)
797: VL vl;
798: P p;
799: DCP *dcp;
800: {
801: int i,j;
802: VL nvl,tvl;
803: VN vn,vnt,vn1;
804: P gg,tmp,p0,g;
805: DCP dc,dc0,dcr,dcr0;
806: int nv,d,d1;
807: int found;
808: VN svn1;
809: P sp0;
810: DCP sdc0;
811:
812: if ( deg(VR(p),p) == 1 ) {
813: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
814: *dcp = dc;
815: return;
816: }
817: clctv(vl,p,&nvl);
818: for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
819: if ( nv == 1 ) {
820: usqp(p,dcp);
821: return;
822: }
823: #if 0
824: if ( heusqfrchk(nvl,p) ) {
825: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
826: *dcp = dc;
827: return;
828: }
829: #endif
830: W_CALLOC(nv,struct oVN,vn);
831: W_CALLOC(nv,struct oVN,vnt);
832: W_CALLOC(nv,struct oVN,vn1);
833: W_CALLOC(nv,struct oVN,svn1);
834: for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
835: vn1[i].v = vn[i].v = tvl->v;
836: vn1[i].v = vn[i].v = 0;
1.5 ! noro 837:
! 838: /* determine a heuristic bound of deg(GCD(p,p')) */
! 839: while ( 1 ) {
! 840: for ( i = 0; vn1[i].v; i++ )
! 841: vn1[i].n = ((unsigned int)random())%256+1;
! 842: substvp(nvl,LC(p),vn1,&tmp);
! 843: if ( tmp ) {
! 844: substvp(nvl,p,vn1,&p0);
! 845: usqp(p0,&dc0);
! 846: for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )
! 847: if ( DEG(dc) )
! 848: d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));
! 849: if ( d1 == 0 ) {
! 850: /* p is squarefree */
! 851: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 852: *dcp = dc;
! 853: return;
! 854: } else {
! 855: d = d1;
! 856: break;
! 857: }
! 858: }
! 859: }
1.1 noro 860:
861: for ( dcr0 = 0, g = p, d = deg(VR(g),g), found = 0; ; ) {
862: while ( 1 ) {
863: for ( i = 0, j = 0; vn[i].v; i++ )
864: if ( vn[i].n ) vnt[j++].v = (V)i;
865: vnt[j].n = 0;
866:
867: mulsgn(vn,vnt,j,vn1);
868: substvp(nvl,LC(g),vn1,&tmp);
869: if ( tmp ) {
870: substvp(nvl,g,vn1,&p0);
871: usqp(p0,&dc0);
872: for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )
873: if ( DEG(dc) )
874: d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));
875:
876: if ( d1 == 0 ) {
877: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = g; NEXT(dc) = 0;
878: if ( !dcr0 )
879: dcr0 = dc;
880: else
881: NEXT(dcr) = dc;
882: *dcp = dcr0;
883: return;
884: }
885:
886: if ( d < d1 )
887: goto END;
888: if ( d > d1 ) {
889: d = d1;
890: found = 1;
891: bcopy((char *)vn1,(char *)svn1,(int)(sizeof(struct oVN)*nv));
892: sp0 = p0; sdc0 = dc0;
893: goto END;
894: }
895: /* d1 == d */
896: if ( found ) {
897: found = 0;
898: #if 0
899: if ( d > d1 ) {
900: d = d1;
901: /*} */
902: #endif
903: msqfrmainmain(nvl,g,svn1,sp0,sdc0,&dc,&gg); g = gg;
904: if ( dc ) {
905: if ( !dcr0 )
906: dcr0 = dc;
907: else
908: NEXT(dcr) = dc;
909: for ( dcr = dc; NEXT(dcr); dcr = NEXT(dcr) );
910: if ( NUM(g) ) {
911: NEXT(dcr) = 0; *dcp = dcr0;
912: return;
913: }
914: d = deg(VR(g),g);
915: }
916: }
917: }
918: END:
919: if ( nextbin(vnt,j) )
920: break;
921: }
922: next(vn);
923: }
924: }
925:
926: void msqfrmainmain(vl,p,vn,p0,dc0,dcp,pp)
927: VL vl;
928: P p;
929: VN vn;
930: P p0;
931: DCP dc0;
932: DCP *dcp;
933: P *pp;
934: {
935: int i,j,k,np;
936: DCP *a;
937: DCP dc,dcr,dcr0,dct;
938: P g,t,s,u,t0,f,f0,d,d0,g0,h0,x,xx;
939: Q q;
940: V v;
941:
942: for ( np = 0, dc = dc0; dc; dc = NEXT(dc), np++ );
943: a = (DCP *)ALLOCA((np + 1)*sizeof(DCP));
944: for ( i = 0, dc = dc0; dc; i++, dc = NEXT(dc) )
945: a[np-i-1] = dc;
946:
947: for ( i = 0, dcr0 = 0, f = p, f0 = p0, v = VR(p);
948: i < np; i++ ) {
949: if ( (i == (np-1))&&UNIQ(DEG(a[i])) ) {
950: NEXTDC(dcr0,dcr);
951: DEG(dcr) = DEG(a[i]);
952: COEF(dcr) = f;
953: f = (P)ONE;
954: } else if ( (i == (np-1))&&UNIQ(DEG(DC(COEF(a[i])))) ) {
955: diffp(vl,f,v,&s); pcp(vl,s,&t,&u);
956: if ( divtpz(vl,f,t,&s) ) {
957: NEXTDC(dcr0,dcr);
958: DEG(dcr) = DEG(a[i]);
959: COEF(dcr) = s;
960: f = (P)ONE;
961: } else
962: break;
963: } else {
964: for ( t = f, t0 = f0,
965: j = 0, k = QTOS(DEG(a[i]))-1; j < k; j++ ) {
966: diffp(vl,t,v,&s); t = s;
967: diffp(vl,t0,v,&s); t0 = s;
968: }
969: factorial(k,&q);
970: divsp(vl,t,(P)q,&s);
971: for ( dct = DC(s); NEXT(dct); dct = NEXT(dct) );
972: if ( DEG(dct) ) {
973: MKV(VR(s),x); pwrp(vl,x,DEG(dct),&xx);
974: divsp(vl,s,xx,&d);
975: } else {
976: xx = (P)ONE; d = s;
977: }
978: divsp(vl,t0,xx,&t);
979: divsp(vl,t,(P)q,&s);
980: ptozp(s,1,(Q *)&t,&d0);
981:
982: for ( dct = DC(COEF(a[i])); NEXT(dct); dct = NEXT(dct) );
983: if ( DEG(dct) )
984: divsp(vl,COEF(a[i]),xx,&g0);
985: else {
986: xx = (P)ONE; g0 = COEF(a[i]);
987: }
988:
989: pcp(vl,d,&t,&u); d = t;
990: ptozp(g0,1,(Q *)&u,&t); g0 = t;
991:
992: {
993: DCP dca,dcb;
994:
995: fctrp(vl,LC(d),&dca);
996: for ( dcb = dca, u = (P)ONE; dcb; dcb = NEXT(dcb) ) {
997: if ( NUM(COEF(dcb)) ) {
998: mulp(vl,u,COEF(dcb),&t); u = t;
999: } else {
1000: Q qq;
1001: P tt;
1002:
1003: pwrp(vl,COEF(dcb),DEG(a[i]),&s);
1004: for ( t = LC(f), j = 0; divtpz(vl,t,s,&tt); j++, t = tt );
1005: STOQ(j,qq);
1006: if ( cmpq(qq,DEG(dcb)) > 0 )
1007: qq = DEG(dcb);
1008: pwrp(vl,COEF(dcb),qq,&t); mulp(vl,u,t,&s); u = s;
1009: }
1010: }
1011: divsp(vl,d0,g0,&h0);
1012: }
1013: mfctrhen2(vl,vn,d,d0,g0,h0,u,LC(d),&s);
1014: if ( s ) {
1015: mulp(vl,s,xx,&g);
1016: pwrp(vl,g,DEG(a[i]),&t);
1017: if ( divtpz(vl,f,t,&s) ) {
1018: NEXTDC(dcr0,dcr);
1019: DEG(dcr) = DEG(a[i]); COEF(dcr) = g;
1020: f = s; substvp(vl,f,vn,&f0);
1021: } else
1022: break;
1023: } else
1024: break;
1025: }
1026: }
1027: *pp = f;
1028: if ( dcr0 )
1029: NEXT(dcr) = 0;
1030: *dcp = dcr0;
1031: }
1032:
1033: void mfctrhen2(vl,vn,f,f0,g0,h0,lcg,lch,gp)
1034: VL vl;
1035: VN vn;
1036: P f;
1037: P f0,g0,h0,lcg,lch;
1038: P *gp;
1039: {
1040: V v;
1041: P f1,lc,lc0,lcg0,lch0;
1042: P m,ff0,gg0,hh0,gk,hk,ggg,gggr,hhh,ak,bk,tmp;
1043: Q bb,qk,s;
1044: Q cbd;
1045: int dbd;
1046: int d;
1047:
1048: if ( NUM(g0) ) {
1049: *gp = (P)ONE;
1050: return;
1051: }
1052:
1053: v = VR(f); d = deg(v,f);
1054: if ( d == deg(v,g0) ) {
1055: pcp(vl,f,gp,&tmp);
1056: return;
1057: }
1058:
1059: mulp(vl,lcg,lch,&lc);
1060: if ( !divtpz(vl,lc,LC(f),(P *)&s) ) {
1061: *gp = 0;
1062: return;
1063: }
1064: mulp(vl,(P)s,f,&f1);
1065: dbd = dbound(VR(f1),f1) + 1; cbound(vl,f1,&cbd);
1066:
1067: substvp(vl,lc,vn,&lc0);
1068: divq((Q)lc0,(Q)LC(f0),(Q *)&m); mulp(vl,f0,m,&ff0);
1069: substvp(vl,lcg,vn,&lcg0);
1070: divq((Q)lcg0,(Q)LC(g0),(Q *)&m); mulp(vl,g0,m,&gg0);
1071: substvp(vl,lch,vn,&lch0);
1072: divq((Q)lch0,(Q)LC(h0),(Q *)&m); mulp(vl,h0,m,&hh0);
1073: addq(cbd,cbd,&bb);
1074: henzq1(gg0,hh0,bb,&bk,&ak,&qk); gk = gg0; hk = hh0;
1075: henmv(vl,vn,f1,gk,hk,ak,bk,
1076: lcg,lch,lcg0,lch0,qk,dbd,&ggg,&hhh);
1077:
1078: if ( divtpz(vl,f1,ggg,&gggr) )
1079: pcp(vl,ggg,gp,&tmp);
1080: else
1081: *gp = 0;
1082: }
1083:
1084: int sqfrchk(p)
1085: P p;
1086: {
1087: Q c;
1088: P f;
1089: DCP dc;
1090:
1091: ptozp(p,SGN((Q)UCOEF(p)),&c,&f); usqp(f,&dc);
1092: if ( NEXT(dc) || !UNIQ(DEG(dc)) )
1093: return ( 0 );
1094: else
1095: return ( 1 );
1096: }
1097:
1098: int cycchk(p)
1099: P p;
1100: {
1101: Q c;
1102: P f;
1103:
1104: ptozp(p,SGN((Q)UCOEF(p)),&c,&f);
1105: if ( iscycp(f) || iscycm(f) )
1106: return 0;
1107: else
1108: return 1;
1109: }
1110:
1111: int zerovpchk(vl,p,vn)
1112: VL vl;
1113: P p;
1114: VN vn;
1115: {
1116: P t;
1117:
1118: substvp(vl,p,vn,&t);
1119: if ( t )
1120: return ( 0 );
1121: else
1122: return ( 1 );
1123: }
1124:
1125: int valideval(vl,dc,vn)
1126: VL vl;
1127: DCP dc;
1128: VN vn;
1129: {
1130: DCP dct;
1131: Q *a;
1132: int i,j,n;
1133: N q,r;
1134:
1135: for ( dct = NEXT(dc), n = 0; dct; dct = NEXT(dct), n++ );
1136: W_CALLOC(n,Q,a);
1137: for ( dct = NEXT(dc), i = 0; dct; dct = NEXT(dct), i++ ) {
1138: substvp(vl,COEF(dct),vn,(P *)&a[i]);
1139: if ( !a[i] )
1140: return ( 0 );
1141:
1142: for ( j = 0; j < i; j++ ) {
1143: divn(NM(a[j]),NM(a[i]),&q,&r);
1144: if ( !r )
1145: return ( 0 );
1146: divn(NM(a[i]),NM(a[j]),&q,&r);
1147: if ( !r )
1148: return ( 0 );
1149: }
1150: }
1151: return ( 1 );
1152: }
1153:
1154: void estimatelc(vl,c,dc,vn,lcp)
1155: VL vl;
1156: Q c;
1157: DCP dc;
1158: VN vn;
1159: P *lcp;
1160: {
1161: int i;
1162: DCP dct;
1163: P r,s,t;
1164: Q c0,c1,c2;
1165:
1166: for ( dct = dc, r = (P)ONE; dct; dct = NEXT(dct) ) {
1167: if ( NUM(COEF(dct)) ) {
1168: mulp(vl,r,COEF(dct),&s); r = s;
1169: } else {
1170: substvp(vl,COEF(dct),vn,(P *)&c0);
1171: for ( i = 0, c1 = c; i < (int)QTOS(DEG(dct)); i++ ) {
1172: divq(c1,c0,&c2);
1173: if ( !INT(c2) )
1174: break;
1175: else
1176: c1 = c2;
1177: }
1178: if ( i ) {
1179: STOQ(i,c1);
1180: pwrp(vl,COEF(dct),c1,&s); mulp(vl,r,s,&t); r = t;
1181: }
1182: }
1183: }
1184: *lcp = r;
1185: }
1186:
1187: void monomialfctr(vl,p,pr,dcp)
1188: VL vl;
1189: P p;
1190: P *pr;
1191: DCP *dcp;
1192: {
1193: VL nvl,avl;
1194: Q d;
1195: P f,t,s;
1196: DCP dc0,dc;
1197:
1198: clctv(vl,p,&nvl);
1199: for ( dc0 = 0, avl = nvl, f = p; avl; avl = NEXT(avl) ) {
1200: getmindeg(avl->v,f,&d);
1201: if ( d ) {
1202: MKV(avl->v,t); NEXTDC(dc0,dc); DEG(dc) = d; COEF(dc) = t;
1203: pwrp(vl,t,d,&s); divsp(vl,f,s,&t); f = t;
1204: }
1205: }
1206: if ( dc0 )
1207: NEXT(dc) = 0;
1208: *pr = f; *dcp = dc0;
1209: }
1210:
1211: void afctr(vl,p0,p,dcp)
1212: VL vl;
1213: P p,p0;
1214: DCP *dcp;
1215: {
1216: DCP dc,dc0,dcr,dct,dcs;
1217: P t;
1218: VL nvl;
1219:
1220: if ( VR(p) == VR(p0) ) {
1221: NEWDC(dc);
1222: DEG(dc) = ONE;
1223: COEF(dc) = p;
1224: NEXT(dc) = 0;
1225: *dcp = dc;
1226: return;
1227: }
1228:
1229: clctv(vl,p,&nvl);
1230: if ( !NEXT(nvl) )
1231: ufctr(p,1,&dc);
1232: else {
1233: sqa(vl,p0,p,&dc);
1234: for ( dct = dc; dct; dct = NEXT(dct) ) {
1235: pmonic(vl,p0,COEF(dct),&t); COEF(dct) = t;
1236: }
1237: }
1238: if ( NUM(COEF(dc)) )
1239: dcr = NEXT(dc);
1240: else
1241: dcr = dc;
1242: for ( dc0 = 0; dcr; dcr = NEXT(dcr) ) {
1243: afctrmain(vl,p0,COEF(dcr),1,&dcs);
1244:
1245: for ( dct = dcs; dct; dct = NEXT(dct) )
1246: DEG(dct) = DEG(dcr);
1247: if ( !dc0 )
1248: dc0 = dcs;
1249: else {
1250: for ( dct = dc0; NEXT(dct); dct = NEXT(dct) );
1251: NEXT(dct) = dcs;
1252: }
1253: }
1254: *dcp = dc0;
1255: }
1256:
1257: void afctrmain(vl,p0,p,init,dcp)
1258: VL vl;
1259: P p,p0;
1260: int init;
1261: DCP *dcp;
1262: {
1263: P x,y,s,m,a,t,u,pt,pt1,res,g;
1264: Q q;
1265: DCP dc,dc0,dcsq0,dcr0,dcr,dct,dcs;
1266: V v,v0;
1267:
1268: if ( !cmpq(DEG(DC(p)),ONE) ) {
1269: NEWDC(dc); DEG(dc) = ONE; NEXT(dc) = 0;
1270: pmonic(vl,p0,p,&COEF(dc)); *dcp = dc;
1271: return;
1272: }
1273:
1274: v = VR(p); MKV(v,x);
1275: v0 = VR(p0); MKV(v0,y);
1276: STOQ(init,q),s = (P)q;
1277: mulp(vl,s,y,&m); subp(vl,x,m,&t); addp(vl,x,m,&a);
1278: substp(vl,p,v,t,&pt);
1279: remsdcp(vl,pt,p0,&pt1);
1280:
1281: /*
1282: if ( ( deg(v0,p0) <= 3 ) || ( TCPQ(p0) <= 2 ) )
1283: resultp(vl,v0,p0,pt1,&res);
1284: else
1285: srcrnorm(vl,v0,pt1,p0,&res);
1286: */
1287: #if 0
1288: srcr(vl,v0,pt1,p0,&res);
1289: #endif
1290: resultp(vl,v0,p0,pt1,&res);
1291: usqp(res,&dcsq0);
1292: for ( dc0 = 0, dct = dcsq0; dct; dct = NEXT(dct) ) {
1293: if ( UNIQ(DEG(dct)) )
1294: ufctr(COEF(dct),deg(v0,p0),&dcs);
1295: else
1296: ufctr(COEF(dct),1,&dcs);
1297: for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
1298: DEG(dcr) = DEG(dct);
1299: if ( !dc0 ) {
1300: dc0 = NEXT(dcs);
1301: dc = dc0;
1302: } else {
1303: for ( ; NEXT(dc); dc = NEXT(dc) );
1304: NEXT(dc) = NEXT(dcs);
1305: }
1306: }
1307: sortfs(&dc0);
1308:
1309: for ( g = p, dcr = dcr0 = 0, dc = dc0; dc; dc = NEXT(dc) ) {
1310: if ( !UNIQ(DEG(dc)) ) {
1311: substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s);
1312: gcda(vl,p0,s,g,&u);
1313: if ( !NUM(u) && (VR(u) == v)) {
1314: afctrmain(vl,p0,u,init+1,&dct);
1315: for ( dcs = dct, t = (P)ONE; dcs; dcs = NEXT(dcs) ) {
1316: mulp(vl,t,COEF(dcs),&s); remsdcp(vl,s,p0,&t);
1317: }
1318: pdiva(vl,p0,g,t,&s); g = s;
1319: if ( !dcr0 )
1320: dcr0 = dct;
1321: else
1322: NEXT(dcr) = dct;
1323: for ( dcr = dct; NEXT(dcr); dcr = NEXT(dcr) );
1324: }
1325: } else {
1326: substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s);
1327: gcda(vl,p0,s,g,&u);
1328: if ( !NUM(u) && (VR(u) == v)) {
1329: NEXTDC(dcr0,dcr);
1330: DEG(dcr) = ONE;
1331: COEF(dcr) = u;
1332: pdiva(vl,p0,g,u,&t); g = t;
1333: }
1334: }
1335: }
1336: if ( dcr0 )
1337: NEXT(dcr) = 0;
1338: *dcp = dcr0;
1339: }