![[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: }