![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Fgfs.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 Fgfs.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.6 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/Fgfs.c,v 1.5 2002/10/23 07:54:58 noro Exp $ */
1.1 noro 2:
3: #include "ca.h"
4:
1.3 noro 5: void cont_pp_mv_sf(VL vl,VL rvl,P p,P *c,P *pp);
6: void gcdsf_main(VL vl,P *pa,int m,P *r);
7: void ugcdsf(P *pa,int m,P *r);
8: void head_monomial(V v,P p,P *coef,P *term);
1.4 noro 9: void sqfrsfmain(VL vl,P f,DCP *dcp);
10: void pthrootsf(P f,Q m,P *r);
11: void partial_sqfrsf(VL vl,V v,P f,P *r,DCP *dcp);
12: void gcdsf(VL vl,P *pa,int k,P *r);
1.5 noro 13: void mfctrsfmain(VL vl, P f, DCP *dcp);
1.6 ! noro 14: int next_evaluation_point(int *mev,int n);
! 15: void estimatelc_sf(VL vl,P c,DCP dc,int *mev,P *lcp);
! 16: void mfctrsf_hensel(VL vl,int *mev,P f,P pp0,P u0,P v0,P lcu,P lcv,P *up);
1.4 noro 17:
18: void lex_lc(P f,P *c)
19: {
20: if ( !f || NUM(f) )
21: *c = f;
22: else
23: lex_lc(COEF(DC(f)),c);
24: }
25:
26: DCP append_dc(DCP dc,DCP dct)
27: {
28: DCP dcs;
29:
30: if ( !dc )
31: return dct;
32: else {
33: for ( dcs = dc; NEXT(dcs); dcs = NEXT(dcs) );
34: NEXT (dcs) = dct;
35: return dc;
36: }
37: }
38:
39: void sqfrsf(VL vl, P f, DCP *dcp)
40: {
41: DCP dc,dct;
42: Obj obj;
43: P t,s,c;
44: VL tvl,nvl;
45:
46: simp_ff((Obj)f,&obj); f = (P)obj;
47: lex_lc(f,&c); divsp(vl,f,c,&t); f = t;
48: monomialfctr(vl,f,&t,&dc); f = t;
49: clctv(vl,f,&tvl); vl = tvl;
50: if ( !vl )
51: ;
52: else if ( !NEXT(vl) ) {
53: sfusqfr(f,&dct);
54: dc = append_dc(dc,NEXT(dct));
55: } else {
56: t = f;
57: for ( tvl = vl; tvl; tvl = NEXT(tvl) ) {
58: reordvar(vl,tvl->v,&nvl);
59: cont_pp_mv_sf(vl,NEXT(nvl),t,&c,&s); t = s;
60: sqfrsf(vl,c,&dct);
61: dc = append_dc(dc,NEXT(dct));
62: }
63: sqfrsfmain(vl,t,&dct);
64: dc = append_dc(dc,dct);
65: }
66: NEWDC(dct); DEG(dct) = ONE; COEF(dct) = (P)c; NEXT(dct) = dc;
67: *dcp = dct;
68: }
69:
70: void sqfrsfmain(VL vl,P f,DCP *dcp)
71: {
72: VL tvl;
73: DCP dc,dct,dcs;
74: P t,s;
75: Q m,m1;
76: V v;
77:
78: clctv(vl,f,&tvl); vl = tvl;
79: dc = 0;
80: t = f;
81: for ( tvl = vl; tvl; tvl = NEXT(tvl) ) {
82: v = tvl->v;
83: partial_sqfrsf(vl,v,t,&s,&dct); t = s;
84: dc = append_dc(dc,dct);
85: }
86: if ( !NUM(t) ) {
87: STOQ(characteristic_sf(),m);
88: pthrootsf(t,m,&s);
89: sqfrsfmain(vl,s,&dct);
90: for ( dcs = dct; dcs; dcs = NEXT(dcs) ) {
91: mulq(DEG(dcs),m,&m1); DEG(dcs) = m1;
92: }
93: dc = append_dc(dc,dct);
94: }
95: *dcp = dc;
96: }
97:
98: void pthrootsf(P f,Q m,P *r)
99: {
100: DCP dc,dc0,dct;
101: N qn,rn;
102:
103: if ( NUM(f) )
104: pthrootgfs(f,r);
105: else {
106: dc = DC(f);
107: dc0 = 0;
108: for ( dc0 = 0; dc; dc = NEXT(dc) ) {
109: NEXTDC(dc0,dct);
110: pthrootsf(COEF(dc),m,&COEF(dct));
111: if ( DEG(dc) ) {
112: divn(NM(DEG(dc)),NM(m),&qn,&rn);
113: if ( rn )
114: error("pthrootsf : cannot happen");
115: NTOQ(qn,1,DEG(dct));
116: } else
117: DEG(dct) = 0;
118: }
119: NEXT(dct) = 0;
120: MKP(VR(f),dc0,*r);
121: }
122: }
123:
124: void partial_sqfrsf(VL vl,V v,P f,P *r,DCP *dcp)
125: {
126: P ps[2];
127: DCP dc0,dc;
128: int m;
129: P t,flat,flat1,g,df,q;
130:
131: diffp(vl,f,v,&df);
132: if ( !df ) {
133: *dcp = 0;
134: *r = f;
135: return;
136: }
137: ps[0] = f; ps[1] = df;
138: gcdsf(vl,ps,2,&g);
139: divsp(vl,f,g,&flat);
140: m = 0;
141: t = f;
142: dc0 = 0;
143: while ( !NUM(flat) ) {
144: while ( divtp(vl,t,flat,&q) ) {
145: t = q; m++;
146: }
147: ps[0] = t; ps[1] = flat;
148: gcdsf(vl,ps,2,&flat1);
149: divsp(vl,flat,flat1,&g);
150: flat = flat1;
151: NEXTDC(dc0,dc);
152: COEF(dc) = g;
153: STOQ(m,DEG(dc));
154: }
155: NEXT(dc) = 0;
156: *dcp = dc0;
157: *r = t;
158: }
1.1 noro 159:
160: void gcdsf(VL vl,P *pa,int k,P *r)
161: {
1.3 noro 162: P *ps,*pl,*pm;
163: P **cp;
1.1 noro 164: int *cn;
165: DCP *ml;
166: Obj obj;
167: int i,j,l,m;
168: P mg,mgsf,t;
169: VL avl,nvl,tvl,svl;
170:
171: ps = (P *)ALLOCA(k*sizeof(P));
172: for ( i = 0, m = 0; i < k; i++ ) {
173: simp_ff((Obj)pa[i],&obj);
174: if ( obj )
1.3 noro 175: ps[m++] = (P)obj;
1.1 noro 176: }
177: if ( !m ) {
178: *r = 0;
179: return;
180: }
181: if ( m == 1 ) {
1.3 noro 182: *r = ps[0];
1.1 noro 183: return;
184: }
185: pl = (P *)ALLOCA(m*sizeof(P));
186: ml = (DCP *)ALLOCA(m*sizeof(DCP));
187: for ( i = 0; i < m; i++ )
188: monomialfctr(vl,ps[i],&pl[i],&ml[i]);
1.3 noro 189: gcdmonomial(vl,ml,m,&mg); simp_ff((Obj)mg,&obj); mgsf = (P)obj;
1.1 noro 190: for ( i = 0, nvl = vl, avl = 0; nvl && i < m; i++ ) {
191: clctv(vl,pl[i],&tvl);
192: intersectv(nvl,tvl,&svl); nvl = svl;
193: mergev(vl,avl,tvl,&svl); avl = svl;
194: }
195: if ( !nvl ) {
196: *r = mgsf;
197: return;
198: }
199: if ( !NEXT(avl) ) {
200: ugcdsf(pl,m,&t);
201: mulp(vl,mgsf,t,r);
202: return;
203: }
204: for ( tvl = nvl, i = 0; tvl; tvl = NEXT(tvl), i++ );
205: for ( tvl = avl, j = 0; tvl; tvl = NEXT(tvl), j++ );
206: if ( i == j ) {
207: /* all the pl[i]'s have the same variables */
208: gcdsf_main(avl,pl,m,&t);
209: } else {
210: cp = (P **)ALLOCA(m*sizeof(P *));
211: cn = (int *)ALLOCA(m*sizeof(int));
212: for ( i = 0; i < m; i++ ) {
213: cp[i] = (P *)ALLOCA(lengthp(pl[i])*sizeof(P));
214: cn[i] = pcoef(vl,nvl,pl[i],cp[i]);
215: }
216: for ( i = j = 0; i < m; i++ )
217: j += cn[i];
218: pm = (P *)ALLOCA(j*sizeof(P));
219: for ( i = l = 0; i < m; i++ )
220: for ( j = 0; j < cn[i]; j++ )
221: pm[l++] = cp[i][j];
222: gcdsf(vl,pm,l,&t);
223: }
224: mulp(vl,mgsf,t,r);
225: }
226:
227: /* univariate gcd */
228:
229: void ugcdsf(P *pa,int m,P *r)
230: {
1.3 noro 231: P *ps;
1.1 noro 232: int i;
233: UM w1,w2,w3,w;
234: int d;
235: V v;
236:
237: if ( m == 1 ) {
238: *r = pa[0];
239: return;
240: }
1.3 noro 241: for ( i = 0; i < m; i++ )
242: if ( NUM(pa[i]) ) {
243: itogfs(1,r);
244: return;
245: }
1.1 noro 246: ps = (P *)ALLOCA(m*sizeof(P));
247: sort_by_deg(m,pa,ps);
1.3 noro 248: v = VR(ps[m-1]);
249: d = getdeg(v,ps[m-1]);
1.1 noro 250: w1 = W_UMALLOC(d);
251: w2 = W_UMALLOC(d);
252: w3 = W_UMALLOC(d);
253: ptosfum(ps[0],w1);
254: for ( i = 1; i < m; i++ ) {
255: ptosfum(ps[i],w2);
256: gcdsfum(w1,w2,w3);
257: w = w1; w1 = w3; w3 = w;
258: if ( !DEG(w1) ) {
1.3 noro 259: itogfs(1,r);
1.1 noro 260: return;
261: }
262: }
263: sfumtop(v,w1,r);
264: }
265:
1.4 noro 266: /* deg(HT(p),v), where p is considered as distributed poly over F[v] */
267: int gethdeg(VL vl,V v,P p)
268: {
269: DCP dc;
270: Q dmax;
271: P cmax;
272:
273: if ( !p )
274: return -1;
275: else if ( NUM(p) )
276: return 0;
277: else if ( VR(p) != v )
278: /* HT(p) = HT(lc(p))*x^D */
279: return gethdeg(vl,v,COEF(DC(p)));
280: else {
281: /* VR(p) = v */
282: dc = DC(p); dmax = DEG(dc); cmax = COEF(dc);
283: for ( dc = NEXT(dc); dc; dc = NEXT(dc) )
284: if ( compp(vl,COEF(dc),cmax) > 0 ) {
285: dmax = DEG(dc); cmax = COEF(dc);
286: }
287: return QTOS(dmax);
288: }
289: }
1.1 noro 290:
291: /* all the pa[i]'s have the same variables (=vl) */
292:
293: void gcdsf_main(VL vl,P *pa,int m,P *r)
294: {
1.3 noro 295: int nv,i,i0,imin,d,d0,d1,d2,dmin,index;
296: V v,v0,vmin;
1.2 noro 297: VL tvl,nvl,rvl,nvl0,rvl0;
1.3 noro 298: P *pc, *ps, *ph,*lps;
299: P x,t,cont,hg,g,hm,mod,s;
300: P hge,ge,ce,he,u,cof1e,mode,mod1,adj,cof1,coadj,q;
301: GFS sf;
1.2 noro 302:
1.1 noro 303: for ( nv = 0, tvl = vl; tvl; tvl = NEXT(tvl), nv++);
304: if ( nv == 1 ) {
305: ugcdsf(pa,m,r);
306: return;
307: }
1.4 noro 308: /* find v s.t. min(deg(pa[i],v)+gethdeg(pa[i],v)) is minimal */
1.1 noro 309: tvl = vl;
310: do {
311: v = tvl->v;
312: i = 0;
313: do {
1.4 noro 314: d = getdeg(v,pa[i])+gethdeg(vl,v,pa[i]);
1.1 noro 315: if ( i == 0 || (d < d0) ) {
316: d0 = d; i0 = i; v0 = v;
317: }
318: } while ( ++i < m );
319: if ( tvl == vl || (d0 < dmin) ) {
320: dmin = d0; imin = i0; vmin = v0;
321: }
322: } while ( tvl = NEXT(tvl) );
323:
324: /* reorder variables so that vmin is the last variable */
325: for ( nvl0 = 0, rvl0 = 0, tvl = vl; tvl; tvl = NEXT(tvl) )
326: if ( tvl->v != vmin ) {
327: NEXTVL(nvl0,nvl); nvl->v = tvl->v;
328: NEXTVL(rvl0,rvl); rvl->v = tvl->v;
329: }
330: /* rvl = remaining variables */
1.3 noro 331: NEXT(rvl) = 0; rvl = rvl0;
1.1 noro 332: /* nvl = ...,vmin */
1.3 noro 333: NEXTVL(nvl0,nvl); nvl->v = vmin; NEXT(nvl) = 0; nvl = nvl0;
1.2 noro 334: MKV(vmin,x);
1.1 noro 335:
336: /* for content and primitive part */
337: pc = (P *)ALLOCA(m*sizeof(P));
338: ps = (P *)ALLOCA(m*sizeof(P));
339: ph = (P *)ALLOCA(m*sizeof(P));
340: /* separate the contents */
341: for ( i = 0; i < m; i++ ) {
1.3 noro 342: reorderp(nvl,vl,pa[i],&t);
1.1 noro 343: cont_pp_mv_sf(nvl,rvl,t,&pc[i],&ps[i]);
344: head_monomial(vmin,ps[i],&ph[i],&t);
345: }
346: ugcdsf(pc,m,&cont);
347: ugcdsf(ph,m,&hg);
348:
349: /* for hg*pp (used in check phase) */
350: lps = (P *)ALLOCA(m*sizeof(P));
351: for ( i = 0; i < m; i++ )
352: mulp(nvl,hg,ps[i],&lps[i]);
353:
354: while ( 1 ) {
355: g = 0;
1.3 noro 356: cof1 = 0;
1.1 noro 357: hm = 0;
1.3 noro 358: itogfs(1,&mod);
1.1 noro 359: index = 0;
1.3 noro 360: for ( index = 0; getdeg(vmin,mod) <= d+1; index++ ) {
1.1 noro 361: /* evaluation pt */
1.3 noro 362: indextogfs(index,&s);
1.1 noro 363: substp(nvl,hg,vmin,s,&hge);
364: if ( !hge )
365: continue;
366: for ( i = 0; i < m; i++ )
367: substp(nvl,ps[i],vmin,s,&ph[i]);
368: /* ge = GCD(ps[0]|x=s,...,ps[m-1]|x=s) */
369: gcdsf(nvl,ph,m,&ge);
370: head_monomial(vmin,ge,&ce,&he);
1.3 noro 371: if ( NUM(he) ) {
1.1 noro 372: *r = cont;
373: return;
374: }
1.3 noro 375: divgfs((GFS)hge,(GFS)ce,&sf); t = (P)sf;
376: mulp(nvl,t,ge,&u); ge = u;
1.1 noro 377: divsp(nvl,ph[imin],ge,&t); mulp(nvl,hge,t,&cof1e);
1.2 noro 378: /* hm=0 : reset; he==hm : lucky */
1.3 noro 379: if ( !hm || !compp(nvl,he,hm) ) {
1.2 noro 380: substp(nvl,mod,vmin,s,&mode); divsp(nvl,mod,mode,&mod1);
381: /* adj = mod/(mod|x=s)*(ge-g|x=s) */
382: substp(nvl,g,vmin,s,&t);
383: subp(nvl,ge,t,&u); mulp(nvl,mod1,u,&adj);
384: /* coadj = mod/(mod|vmin=s)*(cof1e-cof1e|vmin=s) */
385: substp(nvl,cof1,vmin,s,&t);
1.3 noro 386: subp(nvl,cof1e,t,&u); mulp(nvl,mod1,u,&coadj);
1.2 noro 387: if ( !adj ) {
388: /* adj == gcd ? */
389: for ( i = 0; i < m; i++ )
1.3 noro 390: if ( !divtp(nvl,lps[i],g,&t) )
1.2 noro 391: break;
392: if ( i == m ) {
1.3 noro 393: cont_pp_mv_sf(nvl,rvl,g,&t,&u);
1.2 noro 394: mulp(nvl,cont,u,&t);
1.3 noro 395: reorderp(vl,nvl,t,r);
1.2 noro 396: return;
397: }
398: } else if ( !coadj ) {
1.3 noro 399: /* ps[imin]/coadj == gcd ? */
400: if ( divtp(nvl,lps[imin],cof1,&q) ) {
1.2 noro 401: for ( i = 0; i < m; i++ )
402: if ( !divtp(nvl,lps[i],q,&t) )
403: break;
404: if ( i == m ) {
405: cont_pp_mv_sf(nvl,rvl,q,&t,&u);
406: mulp(nvl,cont,u,&t);
1.3 noro 407: reorderp(vl,nvl,t,r);
1.2 noro 408: return;
409: }
410: }
411: }
412: addp(nvl,g,adj,&t); g = t;
413: addp(nvl,cof1,coadj,&t); cof1 = t;
414: subp(nvl,x,s,&t); mulp(nvl,mod,t,&u); mod = u;
415: hm = he;
416: } else {
417: d1 = homdeg(hm); d2 = homdeg(he);
418: if ( d1 < d2 ) /* we use current hm */
419: continue;
420: else if ( d1 > d2 ) {
421: /* use he */
422: g = ge;
423: cof1 = cof1e;
424: hm = he;
425: subp(nvl,x,s,&mod);
426: } else {
427: /* d1==d2, but hm!=he => both are unlucky */
428: g = 0;
429: cof1 = 0;
1.3 noro 430: itogfs(1,&mod);
1.2 noro 431: }
1.1 noro 432: }
433: }
434: }
435: }
436:
437: void head_monomial(V v,P p,P *coef,P *term)
438: {
439: P t,s,u;
440: DCP dc;
441: GFS one;
1.3 noro 442: VL vl;
1.1 noro 443:
1.3 noro 444: itogfs(1,&one);
445: t = (P)one;
1.1 noro 446: while ( 1 ) {
447: if ( NUM(p) || VR(p) == v ) {
448: *coef = p;
449: *term = t;
450: return;
451: } else {
1.3 noro 452: NEWDC(dc);
453: COEF(dc) = (P)one; DEG(dc) = DEG(DC(p));
1.1 noro 454: MKP(VR(p),dc,s);
455: mulp(vl,t,s,&u); t = u;
456: p = COEF(DC(p));
457: }
458: }
459: }
460:
461: void cont_pp_mv_sf(VL vl,VL rvl,P p,P *c,P *pp)
462: {
463: DP dp;
464: MP t;
465: int i,m;
466: P *ps;
467:
468: ptod(vl,rvl,p,&dp);
469: for ( t = BDY(dp), m = 0; t; t = NEXT(t), m++ );
470: ps = (P *)ALLOCA(m*sizeof(P));
1.3 noro 471: for ( t = BDY(dp), i = 0; t; t = NEXT(t), i++ )
1.1 noro 472: ps[i] = C(t);
473: ugcdsf(ps,m,c);
1.3 noro 474: divsp(vl,p,*c,pp);
1.5 noro 475: }
476:
477: void mfctrsf(VL vl, P f, DCP *dcp)
478: {
479: DCP dc0,dc,dct,dcs,dcr;
480: Obj obj;
481:
482: simp_ff((Obj)f,&obj); f = (P)obj;
483: sqfrsf(vl,f,&dct);
484: dc = dc0 = dct; dct = NEXT(dct); NEXT(dc) = 0;
485: for ( ; dct; dct = NEXT(dct) ) {
486: mfctrsfmain(vl,COEF(dct),&dcs);
487: for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
488: DEG(dcr) = DEG(dct);
489: for ( ; NEXT(dc); dc = NEXT(dc) );
490: NEXT(dc) = dcs;
491: }
492: *dcp = dc0;
493: }
494:
495: /* f : sqfr, non const */
496:
497: void mfctrsfmain(VL vl, P f, DCP *dcp)
498: {
1.6 ! noro 499: VL tvl,nvl,rvl;
! 500: DCP dc,dc0,dc1,dc2,dct,lcfdc;
! 501: int imin,inext,i,n,k,np;
1.5 noro 502: int *da;
503: V vx,vy;
504: V *va;
505: P gcd,g,df,dfmin;
506: P pa[2];
1.6 ! noro 507: P g0,pp0,spp0,c,c0,x,y,u,v,lcf,lcu,lcv,u0,v0,t,s;
! 508: GFS ev,evy;
! 509: P *fp0;
! 510: int *mev,*win;
1.5 noro 511:
512: clctv(vl,f,&tvl); vl = tvl;
513: if ( !vl )
514: error("mfctrsfmain : cannot happen");
515: if ( !NEXT(vl) ) {
516: /* univariate */
517: ufctrsf(f,&dc);
518: /* remove lc */
519: *dcp = NEXT(dc);
520: return;
521: }
522: for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
523: va = (V *)ALLOCA(n*sizeof(int));
524: da = (int *)ALLOCA(n*sizeof(int));
525: /* find v s.t. diff(f,v) is nonzero and deg(f,v) is minimal */
526: imin = -1;
527: for ( i = 0, tvl = vl; i < n; tvl = NEXT(tvl), i++ ) {
528: va[i] = tvl->v;
529: da[i] = getdeg(va[i],f);
530: diffp(vl,f,va[i],&df);
531: if ( !df )
532: continue;
533: if ( imin < 0 || da[i] < da[imin] ) {
534: dfmin = df;
535: imin = i;
536: }
537: }
538: /* find v1 neq v s.t. deg(f,v) is minimal */
539: inext = -1;
540: for ( i = 0; i < n; i++ ) {
541: if ( i == imin )
542: continue;
543: if ( inext < 0 || da[i] < da[inext] )
544: inext = i;
545: }
546: pa[0] = f;
547: pa[1] = dfmin;
548: gcdsf_main(vl,pa,2,&gcd);
549: if ( !NUM(gcd) ) {
550: /* f = gcd * f/gcd */
551: mfctrsfmain(vl,gcd,&dc1);
552: divsp(vl,f,gcd,&g);
553: mfctrsfmain(vl,g,&dc2);
554: for ( dct = dc1; NEXT(dct); dct = NEXT(dct) );
555: NEXT(dct) = dc2;
556: *dcp = dc1;
557: return;
558: }
559: /* create vl s.t. vl[0] = va[imin], vl[1] = va[inext] */
560: nvl = 0;
561: NEXTVL(nvl,tvl); tvl->v = va[imin];
562: NEXTVL(nvl,tvl); tvl->v = va[inext];
563: for ( i = 0; i < n; i++ ) {
564: if ( i == imin || i == inext )
565: continue;
566: NEXTVL(nvl,tvl); tvl->v = va[i];
567: }
568: NEXT(tvl) = 0;
569:
570: reorderp(nvl,vl,f,&g);
571: vx = nvl->v;
572: vy = NEXT(nvl)->v;
1.6 ! noro 573: MKV(vx,x);
! 574: MKV(vy,y);
! 575: /* remaining variables */
! 576: rvl = NEXT(NEXT(nvl));
! 577: if ( !rvl ) {
1.5 noro 578: /* bivariate */
579: sfbfctr(g,vx,vy,getdeg(vx,g),&dc1);
580: for ( dc0 = 0; dc1; dc1 = NEXT(dc1) ) {
581: NEXTDC(dc0,dc);
582: DEG(dc) = ONE;
583: reorderp(vl,nvl,COEF(dc1),&COEF(dc));
584: }
585: NEXT(dc) = 0;
586: *dcp = dc0;
587: return;
588: }
1.6 ! noro 589: /* n >= 3; nvl = (vx,vy,X) */
! 590: /* find good evaluation pt for X */
! 591: mev = (int *)CALLOC(n-2,sizeof(int));
! 592: while ( 1 ) {
! 593: for ( g0 = g, tvl = rvl, i = 0; tvl; tvl = NEXT(tvl), i++ ) {
! 594: indextogfs(mev[i],&ev);
! 595: substp(nvl,g0,tvl->v,(P)ev,&t); g0 = t;
! 596: }
! 597: pa[0] = g0;
! 598: diffp(nvl,g0,vx,&pa[1]);
! 599: if ( pa[1] ) {
! 600: gcdsf(nvl,pa,2,&gcd);
! 601: /* XXX maybe we have to accept the case where gcd is a poly of y */
! 602: if ( NUM(gcd) )
! 603: break;
! 604: }
! 605: if ( next_evaluation_point(mev,n-2) )
! 606: error("mfctrsfhmain : short of evaluation points");
! 607: }
! 608: /* g0 = g(x,y,mev) */
! 609: /* separate content; g0 may have the content wrt x */
! 610: cont_pp_sfp(nvl,g0,&c0,&pp0);
! 611:
! 612: /* factorize pp0; spp0 = pp0(x,y+evy) = prod dc */
! 613: sfbfctr_shift(pp0,vx,vy,getdeg(vx,pp0),&evy,&spp0,&dc);
! 614:
! 615: if ( !NEXT(dc) ) {
! 616: /* f is irreducible */
! 617: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0;
! 618: *dcp = dc;
! 619: return;
! 620: }
! 621: /* shift c0; c0 <- c0(y+evy) */
! 622: addp(nvl,y,(P)evy,&t);
! 623: substp(nvl,c0,vy,t,&s);
! 624: c0 = s;
! 625:
! 626: /* now f(x,y+ev,mev) = c0 * prod dc */
! 627: /* factorize lc_x(f) */
! 628: lcf = COEF(DC(f));
! 629: mfctrsf(nvl,lcf,&lcfdc); lcfdc = NEXT(lcfdc);
! 630:
! 631: /* np = number of bivariate factors */
! 632: for ( np = 0, dct = dc; dct; dct = NEXT(dct), np++ );
! 633: fp0 = (P *)ALLOCA((np+1)*sizeof(P));
! 634: for ( i = 0, dct = dc; i < np; dct = NEXT(dct), i++ )
! 635: fp0[i] = COEF(dct);
! 636: fp0[np] = 0;
! 637: win = W_ALLOC(np+1);
! 638: for ( k = 1, win[0] = 1, --np; ; ) {
! 639: itogfs(1,&u0);
! 640: /* u0 = product of selected factors */
! 641: for ( i = 0; i < k; i++ ) {
! 642: mulp(nvl,u0,fp0[win[i]],&t); u0 = t;
! 643: }
! 644: /* we have to consider the content */
! 645: /* g0(y+yev) = c0*u0*v0 */
! 646: mulp(nvl,LC(u0),c0,&c); estimatelc_sf(nvl,c,dc,mev,&lcu);
! 647: divsp(nvl,pp0,u0,&v0);
! 648: mulp(nvl,LC(v0),c0,&c); estimatelc_sf(nvl,c,dc,mev,&lcv);
! 649: mfctrsf_hensel(nvl,mev,f,pp0,u0,v0,lcu,lcv,&u);
! 650: }
! 651: }
! 652:
! 653: int next_evaluation_point(int *mev,int n)
! 654: {
! 655: }
! 656:
! 657: void estimatelc_sf(VL vl,P c,DCP dc,int *mev,P *lcp)
! 658: {
! 659: }
! 660:
! 661: void mfctrsf_hensel(VL vl,int *mev,P f,P pp0,P u0,P v0,P lcu,P lcv,P *up)
! 662: {
1.1 noro 663: }