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