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