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