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