![[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.7 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/Fgfs.c,v 1.6 2002/10/25 02:43:40 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.7 ! noro 14: void next_evaluation_point(int *mev,int n);
! 15: void estimatelc_sf(VL vl,VL rvl,P c,DCP dc,P *lcp);
! 16: void mfctrsf_hensel(VL vl,VL rvl,P f,P pp0,P u0,P v0,P lcu,P lcv,P *up);
! 17: void substvp_sf(VL vl,VL rvl,P f,int *mev,P *r);
! 18: void shift_sf(VL vl, VL rvl, P f, int *mev, int sgn, P *r);
! 19: void adjust_coef_sf(VL vl,VL rvl,P lcu,P u0,P *r);
! 20: void extended_gcd_modyk(P u0,P v0,P *cu,P *cv);
! 21: void poly_to_gfsn_poly(VL vl,P f,V v,P *r);
! 22: void gfsn_poly_to_poly(VL vl,P f,V v,P *r);
1.4 noro 23:
24: void lex_lc(P f,P *c)
25: {
26: if ( !f || NUM(f) )
27: *c = f;
28: else
29: lex_lc(COEF(DC(f)),c);
30: }
31:
32: DCP append_dc(DCP dc,DCP dct)
33: {
34: DCP dcs;
35:
36: if ( !dc )
37: return dct;
38: else {
39: for ( dcs = dc; NEXT(dcs); dcs = NEXT(dcs) );
40: NEXT (dcs) = dct;
41: return dc;
42: }
43: }
44:
45: void sqfrsf(VL vl, P f, DCP *dcp)
46: {
47: DCP dc,dct;
48: Obj obj;
49: P t,s,c;
50: VL tvl,nvl;
51:
52: simp_ff((Obj)f,&obj); f = (P)obj;
53: lex_lc(f,&c); divsp(vl,f,c,&t); f = t;
54: monomialfctr(vl,f,&t,&dc); f = t;
55: clctv(vl,f,&tvl); vl = tvl;
56: if ( !vl )
57: ;
58: else if ( !NEXT(vl) ) {
59: sfusqfr(f,&dct);
60: dc = append_dc(dc,NEXT(dct));
61: } else {
62: t = f;
63: for ( tvl = vl; tvl; tvl = NEXT(tvl) ) {
64: reordvar(vl,tvl->v,&nvl);
65: cont_pp_mv_sf(vl,NEXT(nvl),t,&c,&s); t = s;
66: sqfrsf(vl,c,&dct);
67: dc = append_dc(dc,NEXT(dct));
68: }
69: sqfrsfmain(vl,t,&dct);
70: dc = append_dc(dc,dct);
71: }
72: NEWDC(dct); DEG(dct) = ONE; COEF(dct) = (P)c; NEXT(dct) = dc;
73: *dcp = dct;
74: }
75:
76: void sqfrsfmain(VL vl,P f,DCP *dcp)
77: {
78: VL tvl;
79: DCP dc,dct,dcs;
80: P t,s;
81: Q m,m1;
82: V v;
83:
84: clctv(vl,f,&tvl); vl = tvl;
85: dc = 0;
86: t = f;
87: for ( tvl = vl; tvl; tvl = NEXT(tvl) ) {
88: v = tvl->v;
89: partial_sqfrsf(vl,v,t,&s,&dct); t = s;
90: dc = append_dc(dc,dct);
91: }
92: if ( !NUM(t) ) {
93: STOQ(characteristic_sf(),m);
94: pthrootsf(t,m,&s);
95: sqfrsfmain(vl,s,&dct);
96: for ( dcs = dct; dcs; dcs = NEXT(dcs) ) {
97: mulq(DEG(dcs),m,&m1); DEG(dcs) = m1;
98: }
99: dc = append_dc(dc,dct);
100: }
101: *dcp = dc;
102: }
103:
104: void pthrootsf(P f,Q m,P *r)
105: {
106: DCP dc,dc0,dct;
107: N qn,rn;
108:
109: if ( NUM(f) )
110: pthrootgfs(f,r);
111: else {
112: dc = DC(f);
113: dc0 = 0;
114: for ( dc0 = 0; dc; dc = NEXT(dc) ) {
115: NEXTDC(dc0,dct);
116: pthrootsf(COEF(dc),m,&COEF(dct));
117: if ( DEG(dc) ) {
118: divn(NM(DEG(dc)),NM(m),&qn,&rn);
119: if ( rn )
120: error("pthrootsf : cannot happen");
121: NTOQ(qn,1,DEG(dct));
122: } else
123: DEG(dct) = 0;
124: }
125: NEXT(dct) = 0;
126: MKP(VR(f),dc0,*r);
127: }
128: }
129:
130: void partial_sqfrsf(VL vl,V v,P f,P *r,DCP *dcp)
131: {
132: P ps[2];
133: DCP dc0,dc;
134: int m;
135: P t,flat,flat1,g,df,q;
136:
137: diffp(vl,f,v,&df);
138: if ( !df ) {
139: *dcp = 0;
140: *r = f;
141: return;
142: }
143: ps[0] = f; ps[1] = df;
144: gcdsf(vl,ps,2,&g);
145: divsp(vl,f,g,&flat);
146: m = 0;
147: t = f;
148: dc0 = 0;
149: while ( !NUM(flat) ) {
150: while ( divtp(vl,t,flat,&q) ) {
151: t = q; m++;
152: }
153: ps[0] = t; ps[1] = flat;
154: gcdsf(vl,ps,2,&flat1);
155: divsp(vl,flat,flat1,&g);
156: flat = flat1;
157: NEXTDC(dc0,dc);
158: COEF(dc) = g;
159: STOQ(m,DEG(dc));
160: }
161: NEXT(dc) = 0;
162: *dcp = dc0;
163: *r = t;
164: }
1.1 noro 165:
166: void gcdsf(VL vl,P *pa,int k,P *r)
167: {
1.3 noro 168: P *ps,*pl,*pm;
169: P **cp;
1.1 noro 170: int *cn;
171: DCP *ml;
172: Obj obj;
173: int i,j,l,m;
174: P mg,mgsf,t;
175: VL avl,nvl,tvl,svl;
176:
177: ps = (P *)ALLOCA(k*sizeof(P));
178: for ( i = 0, m = 0; i < k; i++ ) {
179: simp_ff((Obj)pa[i],&obj);
180: if ( obj )
1.3 noro 181: ps[m++] = (P)obj;
1.1 noro 182: }
183: if ( !m ) {
184: *r = 0;
185: return;
186: }
187: if ( m == 1 ) {
1.3 noro 188: *r = ps[0];
1.1 noro 189: return;
190: }
191: pl = (P *)ALLOCA(m*sizeof(P));
192: ml = (DCP *)ALLOCA(m*sizeof(DCP));
193: for ( i = 0; i < m; i++ )
194: monomialfctr(vl,ps[i],&pl[i],&ml[i]);
1.3 noro 195: gcdmonomial(vl,ml,m,&mg); simp_ff((Obj)mg,&obj); mgsf = (P)obj;
1.1 noro 196: for ( i = 0, nvl = vl, avl = 0; nvl && i < m; i++ ) {
197: clctv(vl,pl[i],&tvl);
198: intersectv(nvl,tvl,&svl); nvl = svl;
199: mergev(vl,avl,tvl,&svl); avl = svl;
200: }
201: if ( !nvl ) {
202: *r = mgsf;
203: return;
204: }
205: if ( !NEXT(avl) ) {
206: ugcdsf(pl,m,&t);
207: mulp(vl,mgsf,t,r);
208: return;
209: }
210: for ( tvl = nvl, i = 0; tvl; tvl = NEXT(tvl), i++ );
211: for ( tvl = avl, j = 0; tvl; tvl = NEXT(tvl), j++ );
212: if ( i == j ) {
213: /* all the pl[i]'s have the same variables */
214: gcdsf_main(avl,pl,m,&t);
215: } else {
216: cp = (P **)ALLOCA(m*sizeof(P *));
217: cn = (int *)ALLOCA(m*sizeof(int));
218: for ( i = 0; i < m; i++ ) {
219: cp[i] = (P *)ALLOCA(lengthp(pl[i])*sizeof(P));
220: cn[i] = pcoef(vl,nvl,pl[i],cp[i]);
221: }
222: for ( i = j = 0; i < m; i++ )
223: j += cn[i];
224: pm = (P *)ALLOCA(j*sizeof(P));
225: for ( i = l = 0; i < m; i++ )
226: for ( j = 0; j < cn[i]; j++ )
227: pm[l++] = cp[i][j];
228: gcdsf(vl,pm,l,&t);
229: }
230: mulp(vl,mgsf,t,r);
231: }
232:
233: /* univariate gcd */
234:
235: void ugcdsf(P *pa,int m,P *r)
236: {
1.3 noro 237: P *ps;
1.1 noro 238: int i;
239: UM w1,w2,w3,w;
240: int d;
241: V v;
242:
243: if ( m == 1 ) {
244: *r = pa[0];
245: return;
246: }
1.3 noro 247: for ( i = 0; i < m; i++ )
248: if ( NUM(pa[i]) ) {
249: itogfs(1,r);
250: return;
251: }
1.1 noro 252: ps = (P *)ALLOCA(m*sizeof(P));
253: sort_by_deg(m,pa,ps);
1.3 noro 254: v = VR(ps[m-1]);
255: d = getdeg(v,ps[m-1]);
1.1 noro 256: w1 = W_UMALLOC(d);
257: w2 = W_UMALLOC(d);
258: w3 = W_UMALLOC(d);
259: ptosfum(ps[0],w1);
260: for ( i = 1; i < m; i++ ) {
261: ptosfum(ps[i],w2);
262: gcdsfum(w1,w2,w3);
263: w = w1; w1 = w3; w3 = w;
264: if ( !DEG(w1) ) {
1.3 noro 265: itogfs(1,r);
1.1 noro 266: return;
267: }
268: }
269: sfumtop(v,w1,r);
270: }
271:
1.4 noro 272: /* deg(HT(p),v), where p is considered as distributed poly over F[v] */
273: int gethdeg(VL vl,V v,P p)
274: {
275: DCP dc;
276: Q dmax;
277: P cmax;
278:
279: if ( !p )
280: return -1;
281: else if ( NUM(p) )
282: return 0;
283: else if ( VR(p) != v )
284: /* HT(p) = HT(lc(p))*x^D */
285: return gethdeg(vl,v,COEF(DC(p)));
286: else {
287: /* VR(p) = v */
288: dc = DC(p); dmax = DEG(dc); cmax = COEF(dc);
289: for ( dc = NEXT(dc); dc; dc = NEXT(dc) )
290: if ( compp(vl,COEF(dc),cmax) > 0 ) {
291: dmax = DEG(dc); cmax = COEF(dc);
292: }
293: return QTOS(dmax);
294: }
295: }
1.1 noro 296:
297: /* all the pa[i]'s have the same variables (=vl) */
298:
299: void gcdsf_main(VL vl,P *pa,int m,P *r)
300: {
1.3 noro 301: int nv,i,i0,imin,d,d0,d1,d2,dmin,index;
302: V v,v0,vmin;
1.2 noro 303: VL tvl,nvl,rvl,nvl0,rvl0;
1.3 noro 304: P *pc, *ps, *ph,*lps;
305: P x,t,cont,hg,g,hm,mod,s;
306: P hge,ge,ce,he,u,cof1e,mode,mod1,adj,cof1,coadj,q;
307: GFS sf;
1.2 noro 308:
1.1 noro 309: for ( nv = 0, tvl = vl; tvl; tvl = NEXT(tvl), nv++);
310: if ( nv == 1 ) {
311: ugcdsf(pa,m,r);
312: return;
313: }
1.4 noro 314: /* find v s.t. min(deg(pa[i],v)+gethdeg(pa[i],v)) is minimal */
1.1 noro 315: tvl = vl;
316: do {
317: v = tvl->v;
318: i = 0;
319: do {
1.4 noro 320: d = getdeg(v,pa[i])+gethdeg(vl,v,pa[i]);
1.1 noro 321: if ( i == 0 || (d < d0) ) {
322: d0 = d; i0 = i; v0 = v;
323: }
324: } while ( ++i < m );
325: if ( tvl == vl || (d0 < dmin) ) {
326: dmin = d0; imin = i0; vmin = v0;
327: }
328: } while ( tvl = NEXT(tvl) );
329:
330: /* reorder variables so that vmin is the last variable */
331: for ( nvl0 = 0, rvl0 = 0, tvl = vl; tvl; tvl = NEXT(tvl) )
332: if ( tvl->v != vmin ) {
333: NEXTVL(nvl0,nvl); nvl->v = tvl->v;
334: NEXTVL(rvl0,rvl); rvl->v = tvl->v;
335: }
336: /* rvl = remaining variables */
1.3 noro 337: NEXT(rvl) = 0; rvl = rvl0;
1.1 noro 338: /* nvl = ...,vmin */
1.3 noro 339: NEXTVL(nvl0,nvl); nvl->v = vmin; NEXT(nvl) = 0; nvl = nvl0;
1.2 noro 340: MKV(vmin,x);
1.1 noro 341:
342: /* for content and primitive part */
343: pc = (P *)ALLOCA(m*sizeof(P));
344: ps = (P *)ALLOCA(m*sizeof(P));
345: ph = (P *)ALLOCA(m*sizeof(P));
346: /* separate the contents */
347: for ( i = 0; i < m; i++ ) {
1.3 noro 348: reorderp(nvl,vl,pa[i],&t);
1.1 noro 349: cont_pp_mv_sf(nvl,rvl,t,&pc[i],&ps[i]);
350: head_monomial(vmin,ps[i],&ph[i],&t);
351: }
352: ugcdsf(pc,m,&cont);
353: ugcdsf(ph,m,&hg);
354:
355: /* for hg*pp (used in check phase) */
356: lps = (P *)ALLOCA(m*sizeof(P));
357: for ( i = 0; i < m; i++ )
358: mulp(nvl,hg,ps[i],&lps[i]);
359:
360: while ( 1 ) {
361: g = 0;
1.3 noro 362: cof1 = 0;
1.1 noro 363: hm = 0;
1.3 noro 364: itogfs(1,&mod);
1.1 noro 365: index = 0;
1.3 noro 366: for ( index = 0; getdeg(vmin,mod) <= d+1; index++ ) {
1.1 noro 367: /* evaluation pt */
1.3 noro 368: indextogfs(index,&s);
1.1 noro 369: substp(nvl,hg,vmin,s,&hge);
370: if ( !hge )
371: continue;
372: for ( i = 0; i < m; i++ )
373: substp(nvl,ps[i],vmin,s,&ph[i]);
374: /* ge = GCD(ps[0]|x=s,...,ps[m-1]|x=s) */
375: gcdsf(nvl,ph,m,&ge);
376: head_monomial(vmin,ge,&ce,&he);
1.3 noro 377: if ( NUM(he) ) {
1.1 noro 378: *r = cont;
379: return;
380: }
1.3 noro 381: divgfs((GFS)hge,(GFS)ce,&sf); t = (P)sf;
382: mulp(nvl,t,ge,&u); ge = u;
1.1 noro 383: divsp(nvl,ph[imin],ge,&t); mulp(nvl,hge,t,&cof1e);
1.2 noro 384: /* hm=0 : reset; he==hm : lucky */
1.3 noro 385: if ( !hm || !compp(nvl,he,hm) ) {
1.2 noro 386: substp(nvl,mod,vmin,s,&mode); divsp(nvl,mod,mode,&mod1);
387: /* adj = mod/(mod|x=s)*(ge-g|x=s) */
388: substp(nvl,g,vmin,s,&t);
389: subp(nvl,ge,t,&u); mulp(nvl,mod1,u,&adj);
390: /* coadj = mod/(mod|vmin=s)*(cof1e-cof1e|vmin=s) */
391: substp(nvl,cof1,vmin,s,&t);
1.3 noro 392: subp(nvl,cof1e,t,&u); mulp(nvl,mod1,u,&coadj);
1.2 noro 393: if ( !adj ) {
394: /* adj == gcd ? */
395: for ( i = 0; i < m; i++ )
1.3 noro 396: if ( !divtp(nvl,lps[i],g,&t) )
1.2 noro 397: break;
398: if ( i == m ) {
1.3 noro 399: cont_pp_mv_sf(nvl,rvl,g,&t,&u);
1.2 noro 400: mulp(nvl,cont,u,&t);
1.3 noro 401: reorderp(vl,nvl,t,r);
1.2 noro 402: return;
403: }
404: } else if ( !coadj ) {
1.3 noro 405: /* ps[imin]/coadj == gcd ? */
406: if ( divtp(nvl,lps[imin],cof1,&q) ) {
1.2 noro 407: for ( i = 0; i < m; i++ )
408: if ( !divtp(nvl,lps[i],q,&t) )
409: break;
410: if ( i == m ) {
411: cont_pp_mv_sf(nvl,rvl,q,&t,&u);
412: mulp(nvl,cont,u,&t);
1.3 noro 413: reorderp(vl,nvl,t,r);
1.2 noro 414: return;
415: }
416: }
417: }
418: addp(nvl,g,adj,&t); g = t;
419: addp(nvl,cof1,coadj,&t); cof1 = t;
420: subp(nvl,x,s,&t); mulp(nvl,mod,t,&u); mod = u;
421: hm = he;
422: } else {
423: d1 = homdeg(hm); d2 = homdeg(he);
424: if ( d1 < d2 ) /* we use current hm */
425: continue;
426: else if ( d1 > d2 ) {
427: /* use he */
428: g = ge;
429: cof1 = cof1e;
430: hm = he;
431: subp(nvl,x,s,&mod);
432: } else {
433: /* d1==d2, but hm!=he => both are unlucky */
434: g = 0;
435: cof1 = 0;
1.3 noro 436: itogfs(1,&mod);
1.2 noro 437: }
1.1 noro 438: }
439: }
440: }
441: }
442:
443: void head_monomial(V v,P p,P *coef,P *term)
444: {
445: P t,s,u;
446: DCP dc;
447: GFS one;
1.3 noro 448: VL vl;
1.1 noro 449:
1.3 noro 450: itogfs(1,&one);
451: t = (P)one;
1.1 noro 452: while ( 1 ) {
453: if ( NUM(p) || VR(p) == v ) {
454: *coef = p;
455: *term = t;
456: return;
457: } else {
1.3 noro 458: NEWDC(dc);
459: COEF(dc) = (P)one; DEG(dc) = DEG(DC(p));
1.1 noro 460: MKP(VR(p),dc,s);
461: mulp(vl,t,s,&u); t = u;
462: p = COEF(DC(p));
463: }
464: }
465: }
466:
467: void cont_pp_mv_sf(VL vl,VL rvl,P p,P *c,P *pp)
468: {
469: DP dp;
470: MP t;
471: int i,m;
472: P *ps;
473:
474: ptod(vl,rvl,p,&dp);
475: for ( t = BDY(dp), m = 0; t; t = NEXT(t), m++ );
476: ps = (P *)ALLOCA(m*sizeof(P));
1.3 noro 477: for ( t = BDY(dp), i = 0; t; t = NEXT(t), i++ )
1.1 noro 478: ps[i] = C(t);
479: ugcdsf(ps,m,c);
1.3 noro 480: divsp(vl,p,*c,pp);
1.5 noro 481: }
482:
483: void mfctrsf(VL vl, P f, DCP *dcp)
484: {
485: DCP dc0,dc,dct,dcs,dcr;
486: Obj obj;
487:
488: simp_ff((Obj)f,&obj); f = (P)obj;
489: sqfrsf(vl,f,&dct);
490: dc = dc0 = dct; dct = NEXT(dct); NEXT(dc) = 0;
491: for ( ; dct; dct = NEXT(dct) ) {
492: mfctrsfmain(vl,COEF(dct),&dcs);
493: for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
494: DEG(dcr) = DEG(dct);
495: for ( ; NEXT(dc); dc = NEXT(dc) );
496: NEXT(dc) = dcs;
497: }
498: *dcp = dc0;
499: }
500:
501: /* f : sqfr, non const */
502:
503: void mfctrsfmain(VL vl, P f, DCP *dcp)
504: {
1.6 noro 505: VL tvl,nvl,rvl;
1.7 ! noro 506: DCP dc,dc0,dc1,dc2,dct,lcfdc,dcs;
! 507: int imin,inext,i,j,n,k,np;
1.5 noro 508: int *da;
509: V vx,vy;
510: V *va;
1.7 ! noro 511: P *l,*tl;
1.5 noro 512: P gcd,g,df,dfmin;
513: P pa[2];
1.6 noro 514: P g0,pp0,spp0,c,c0,x,y,u,v,lcf,lcu,lcv,u0,v0,t,s;
1.7 ! noro 515: P ype,yme;
1.6 noro 516: GFS ev,evy;
517: P *fp0;
518: int *mev,*win;
1.5 noro 519:
520: clctv(vl,f,&tvl); vl = tvl;
521: if ( !vl )
522: error("mfctrsfmain : cannot happen");
523: if ( !NEXT(vl) ) {
524: /* univariate */
525: ufctrsf(f,&dc);
526: /* remove lc */
527: *dcp = NEXT(dc);
528: return;
529: }
530: for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
531: va = (V *)ALLOCA(n*sizeof(int));
532: da = (int *)ALLOCA(n*sizeof(int));
533: /* find v s.t. diff(f,v) is nonzero and deg(f,v) is minimal */
534: imin = -1;
535: for ( i = 0, tvl = vl; i < n; tvl = NEXT(tvl), i++ ) {
536: va[i] = tvl->v;
537: da[i] = getdeg(va[i],f);
538: diffp(vl,f,va[i],&df);
539: if ( !df )
540: continue;
541: if ( imin < 0 || da[i] < da[imin] ) {
542: dfmin = df;
543: imin = i;
544: }
545: }
546: /* find v1 neq v s.t. deg(f,v) is minimal */
547: inext = -1;
548: for ( i = 0; i < n; i++ ) {
549: if ( i == imin )
550: continue;
551: if ( inext < 0 || da[i] < da[inext] )
552: inext = i;
553: }
554: pa[0] = f;
555: pa[1] = dfmin;
556: gcdsf_main(vl,pa,2,&gcd);
557: if ( !NUM(gcd) ) {
558: /* f = gcd * f/gcd */
559: mfctrsfmain(vl,gcd,&dc1);
560: divsp(vl,f,gcd,&g);
561: mfctrsfmain(vl,g,&dc2);
562: for ( dct = dc1; NEXT(dct); dct = NEXT(dct) );
563: NEXT(dct) = dc2;
564: *dcp = dc1;
565: return;
566: }
567: /* create vl s.t. vl[0] = va[imin], vl[1] = va[inext] */
568: nvl = 0;
569: NEXTVL(nvl,tvl); tvl->v = va[imin];
570: NEXTVL(nvl,tvl); tvl->v = va[inext];
571: for ( i = 0; i < n; i++ ) {
572: if ( i == imin || i == inext )
573: continue;
574: NEXTVL(nvl,tvl); tvl->v = va[i];
575: }
576: NEXT(tvl) = 0;
577:
578: reorderp(nvl,vl,f,&g);
579: vx = nvl->v;
580: vy = NEXT(nvl)->v;
1.6 noro 581: MKV(vx,x);
582: MKV(vy,y);
583: /* remaining variables */
584: rvl = NEXT(NEXT(nvl));
585: if ( !rvl ) {
1.5 noro 586: /* bivariate */
587: sfbfctr(g,vx,vy,getdeg(vx,g),&dc1);
588: for ( dc0 = 0; dc1; dc1 = NEXT(dc1) ) {
589: NEXTDC(dc0,dc);
590: DEG(dc) = ONE;
591: reorderp(vl,nvl,COEF(dc1),&COEF(dc));
592: }
593: NEXT(dc) = 0;
594: *dcp = dc0;
595: return;
596: }
1.6 noro 597: /* n >= 3; nvl = (vx,vy,X) */
598: /* find good evaluation pt for X */
599: mev = (int *)CALLOC(n-2,sizeof(int));
600: while ( 1 ) {
1.7 ! noro 601: substvp_sf(nvl,rvl,g,mev,&g0);
1.6 noro 602: pa[0] = g0;
603: diffp(nvl,g0,vx,&pa[1]);
604: if ( pa[1] ) {
605: gcdsf(nvl,pa,2,&gcd);
606: /* XXX maybe we have to accept the case where gcd is a poly of y */
607: if ( NUM(gcd) )
608: break;
609: }
1.7 ! noro 610: /* XXX if generated indices exceed q of GF(q) => error in indextogfs */
! 611: next_evaluation_point(mev,n-2);
1.6 noro 612: }
613: /* g0 = g(x,y,mev) */
614: /* separate content; g0 may have the content wrt x */
615: cont_pp_sfp(nvl,g0,&c0,&pp0);
616:
1.7 ! noro 617: /* factorize pp0; pp0 = pp0(x,y+evy) = prod dc */
! 618: sfbfctr_shift(pp0,vx,vy,getdeg(vx,pp0),&evy,&spp0,&dc); pp0 = spp0;
1.6 noro 619:
620: if ( !NEXT(dc) ) {
621: /* f is irreducible */
622: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0;
623: *dcp = dc;
624: return;
625: }
1.7 ! noro 626: /* ype = y+evy, yme = y-evy */
! 627: addp(nvl,y,(P)evy,&ype); subp(nvl,y,(P)evy,&yme);
! 628:
1.6 noro 629: /* shift c0; c0 <- c0(y+evy) */
1.7 ! noro 630: substp(nvl,c0,vy,ype,&s); c0 = s;
! 631:
! 632: /* shift f; f <- f(y+evy) */
! 633: substp(nvl,f,vy,ype,&s); f = s;
! 634:
! 635: /* now f(x,0,mev) = c0 * prod dc */
1.6 noro 636:
637: /* factorize lc_x(f) */
638: lcf = COEF(DC(f));
1.7 ! noro 639: mfctrsf(nvl,lcf,&dct);
! 640: /* skip the first element (= a number) */
! 641: dct = NEXT(dct);
! 642:
! 643: /* shift lcfdc; c <- c(X+mev) */
! 644: for ( lcfdc = 0; dct; dct = NEXT(dct) ) {
! 645: NEXTDC(lcfdc,dcs);
! 646: DEG(dcs) = DEG(dct);
! 647: shift_sf(nvl,rvl,COEF(dct),mev,1,&COEF(dcs));
! 648: }
! 649: NEXT(dcs) = 0;
1.6 noro 650:
651: /* np = number of bivariate factors */
652: for ( np = 0, dct = dc; dct; dct = NEXT(dct), np++ );
653: fp0 = (P *)ALLOCA((np+1)*sizeof(P));
654: for ( i = 0, dct = dc; i < np; dct = NEXT(dct), i++ )
655: fp0[i] = COEF(dct);
656: fp0[np] = 0;
1.7 ! noro 657: l = tl = (P *)ALLOCA((np+1)*sizeof(P));
1.6 noro 658: win = W_ALLOC(np+1);
1.7 ! noro 659:
! 660: /* f <- f(X+mev) */
! 661: shift_sf(nvl,rvl,f,mev,1,&s); f = s;
! 662:
1.6 noro 663: for ( k = 1, win[0] = 1, --np; ; ) {
664: itogfs(1,&u0);
665: /* u0 = product of selected factors */
666: for ( i = 0; i < k; i++ ) {
667: mulp(nvl,u0,fp0[win[i]],&t); u0 = t;
668: }
669: /* we have to consider the content */
1.7 ! noro 670: /* g0 = c0*u0*v0 */
! 671: mulp(nvl,LC(u0),c0,&c); estimatelc_sf(nvl,rvl,c,lcfdc,&lcu);
1.6 noro 672: divsp(nvl,pp0,u0,&v0);
1.7 ! noro 673: mulp(nvl,LC(v0),c0,&c); estimatelc_sf(nvl,rvl,c,lcfdc,&lcv);
! 674: mfctrsf_hensel(nvl,rvl,f,pp0,u0,v0,lcu,lcv,&u);
! 675: if ( u ) {
! 676: /* save the factor */
! 677: reorderp(vl,nvl,u,&t);
! 678: /* x -> x-mev, y -> y-evy */
! 679: shift_sf(vl,rvl,t,mev,-1,&s); substp(vl,s,vy,yme,tl++);
! 680:
! 681: /* update f,pp0 */
! 682: divsp(nvl,f,u,&t); f = t;
! 683: divsp(nvl,pp0,u0,&t); pp0 = t;
! 684: /* update win, fp0 */
! 685: for ( i = 0; i < k-1; i++ )
! 686: for ( j = win[i]+1; j < win[i+1]; j++ )
! 687: fp0[j-i-1] = fp0[j];
! 688: for ( j = win[k-1]+1; j <= np; j++ )
! 689: fp0[j-k] = fp0[j];
! 690: if ( ( np -= k ) < k ) break;
! 691: if ( np-win[0]+1 < k )
! 692: if ( ++k <= np ) {
! 693: for ( i = 0; i < k; i++ )
! 694: win[i] = i + 1;
! 695: continue;
! 696: } else
! 697: break;
! 698: else
! 699: for ( i = 1; i < k; i++ )
! 700: win[i] = win[0] + i;
! 701: } else {
! 702: if ( ncombi(1,np,k,win) == 0 )
! 703: if ( k == np ) break;
! 704: else
! 705: for ( i = 0, ++k; i < k; i++ )
! 706: win[i] = i + 1;
! 707: }
! 708: reorderp(vl,nvl,f,&t);
! 709: /* x -> x-mev, y -> y-evy */
! 710: shift_sf(vl,rvl,t,mev,-1,&s); substp(vl,s,vy,yme,tl++);
! 711: *tl = 0;
! 712:
! 713: for ( dc0 = 0, i = 0; l[i]; i++ ) {
! 714: NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = l[i];
! 715: }
! 716: NEXT(dc) = 0; *dcp = dc0;
1.6 noro 717: }
718: }
719:
1.7 ! noro 720: void next_evaluation_point(int *e,int n)
! 721: {
! 722: int i,t,j;
! 723:
! 724: for ( i = n-1; i >= 0; i-- )
! 725: if ( e[i] ) break;
! 726: if ( i < 0 ) e[n-1] = 1;
! 727: else if ( i == 0 ) {
! 728: t = e[0]; e[0] = 0; e[n-1] = t+1;
! 729: } else {
! 730: e[i-1]++; t = e[i];
! 731: for ( j = i; j < n-1; j++ )
! 732: e[j] = 0;
! 733: e[n-1] = t-1;
! 734: }
! 735: }
! 736:
! 737: /*
! 738: * dc : f1^E1*...*fk^Ek
! 739: * find e1,...,ek s.t. fi(0)^ei | c
! 740: * and return f1^e1*...*fk^ek
! 741: * vl = (vx,vy,rvl)
! 742: */
! 743:
! 744: void estimatelc_sf(VL vl,VL rvl,P c,DCP dc,P *lcp)
! 745: {
! 746: DCP dct;
! 747: P r,c1,c2,t,s,f;
! 748: int i,d;
! 749: Q q;
! 750:
! 751: for ( dct = dc, r = (P)ONE; dct; dct = NEXT(dct) ) {
! 752: if ( NUM(COEF(dct)) )
! 753: continue;
! 754: /* constant part */
! 755: substvp_sf(vl,rvl,COEF(dct),0,&f);
! 756: d = QTOS(DEG(dct));
! 757: for ( i = 0, c1 = c; i < d; i++ )
! 758: if ( !divtp(vl,c1,f,&c2) )
! 759: break;
! 760: else
! 761: c1 = c2;
! 762: if ( i ) {
! 763: STOQ(i,q);
! 764: pwrp(vl,COEF(dct),q,&s); mulp(vl,r,s,&t); r = t;
! 765: }
! 766: }
! 767: *lcp = r;
! 768: }
! 769:
! 770: void substvp_sf(VL vl,VL rvl,P f,int *mev,P *r)
! 771: {
! 772: int i;
! 773: VL tvl;
! 774: P g,t;
! 775: GFS ev;
! 776:
! 777: for ( g = f, i = 0, tvl = rvl; tvl; tvl = NEXT(tvl), i++ ) {
! 778: if ( !mev )
! 779: ev = 0;
! 780: else
! 781: indextogfs(mev[i],&ev);
! 782: substp(vl,g,tvl->v,(P)ev,&t); g = t;
! 783: }
! 784: *r = g;
! 785: }
! 786:
! 787: /*
! 788: * f <- f(X+sgn*mev)
! 789: */
! 790:
! 791: void shift_sf(VL vl, VL rvl, P f, int *mev, int sgn, P *r)
! 792: {
! 793: VL tvl;
! 794: int i;
! 795: P x,g,t,s;
! 796: GFS ev;
! 797:
! 798: for ( g = f, tvl = rvl, i = 0; tvl; tvl = NEXT(tvl), i++ ) {
! 799: if ( !mev[i] )
! 800: continue;
! 801: indextogfs(mev[i],&ev);
! 802: MKV(tvl->v,x);
! 803: if ( sgn > 0 )
! 804: addp(vl,x,(P)ev,&t);
! 805: else
! 806: subp(vl,x,(P)ev,&t);
! 807: substp(vl,g,tvl->v,t,&s); g = s;
! 808: }
! 809: *r = g;
! 810: }
! 811:
! 812: /*
! 813: * pp(f(0)) = u0*v0
! 814: */
! 815:
! 816: void mfctrsf_hensel(VL vl,VL rvl,P f,P pp0,P u0,P v0,P lcu,P lcv,P *up)
! 817: {
! 818: VL tvl,onevl;
! 819: P t,s,w,u,v,ff,si,wu,wv,fj,cont;
! 820: UM ydy;
! 821: V vx,vy;
! 822: int dy,n,i,dbd,nv,j;
! 823: int *md;
! 824: P *uh,*vh;
! 825: P x,du0,dv0,m,q,r;
! 826: P *cu,*cv;
! 827: GFSN inv;
! 828:
! 829: /* adjust coeffs */
! 830: /* u0 = am x^m+ ... -> lcu*x^m + a(m-1)*(lcu(0)/am)*x^(m-1)+... */
! 831: /* v0 = bm x^l+ ... -> lcv*x^l + b(l-1)*(lcv(0)/bl)*x^(l-1)+... */
! 832: adjust_coef_sf(vl,rvl,lcu,u0,&u);
! 833: adjust_coef_sf(vl,rvl,lcv,v0,&v);
! 834: vx = vl->v; vy = NEXT(vl)->v;
! 835: n = getdeg(vx,f);
! 836: dy = getdeg(vy,f)+1;
! 837: MKV(vx,x);
! 838: cu = (P *)ALLOCA((n+1)*sizeof(P));
! 839: cv = (P *)ALLOCA((n+1)*sizeof(P));
! 840:
! 841: /* ydy = y^dy */
! 842: ydy = C_UMALLOC(dy); COEF(ydy)[dy] = 1;
! 843: setmod_gfsn(ydy);
! 844:
! 845: /* (R[y]/(y^dy))[x,X] */
! 846: poly_to_gfsn_poly(vl,f,vy,&t); ff = t;
! 847: poly_to_gfsn_poly(vl,u,vy,&t); u = t;
! 848: poly_to_gfsn_poly(vl,v,vy,&t); v = t;
! 849: substvp_sf(vl,rvl,u,0,&u0);
! 850: substvp_sf(vl,rvl,v,0,&v0);
! 851:
! 852: /* compute a(x,y), b(x,y) s.t. a*u0+b*v0 = 1 mod y^dy */
! 853: extended_gcd_modyk(u0,v0,&cu[0],&cv[0]);
! 854:
! 855: /* du0 = LC(u0)^(-1)*u0 mod y^dy */
! 856: /* dv0 = LC(v0)^(-1)*v0 mod y^dy */
! 857: invgfsn((GFSN)LC(u0),&inv); mulp(vl,u0,(P)inv,&du0);
! 858: invgfsn((GFSN)LC(v0),&inv); mulp(vl,v0,(P)inv,&dv0);
! 859:
! 860: /* cu[i]*u0+cv[i]*v0 = x^i mod y^dy */
! 861: for ( i = 1; i <= n; i++ ) {
! 862: mulp(vl,x,cu[i-1],&m); divsrp(vl,m,dv0,&q,&cu[i]);
! 863: mulp(vl,x,cv[i-1],&m); divsrp(vl,m,du0,&q,&cv[i]);
! 864: }
! 865: dbd = dbound(vx,f)+1;
! 866:
! 867: /* extract homogeneous parts */
! 868: W_CALLOC(dbd,P,uh); W_CALLOC(dbd,P,vh);
! 869: for ( i = 0; i <= dbd; i++ ) {
! 870: exthpc(vl,vx,u,i,&uh[i]); exthpc(vl,vx,v,i,&vh[i]);
! 871: }
! 872:
! 873: /* register degrees in each variables */
! 874: for ( nv = 0, tvl = rvl; tvl; tvl = NEXT(tvl), nv++ );
! 875: md = (int *)ALLOCA(nv*sizeof(int));
! 876: for ( i = 0, tvl = rvl; i < nv; tvl = NEXT(tvl), i++ )
! 877: md[i] = getdeg(tvl->v,ff);
! 878:
! 879: /* XXX for removing content of factor wrt vx */
! 880: NEWVL(onevl); onevl->v = vx; NEXT(onevl) = 0;
! 881:
! 882: for ( j = 1; j <= dbd; j++ ) {
! 883: for ( i = 0, tvl = rvl; i < nv; tvl = NEXT(tvl), i++ )
! 884: if ( getdeg(tvl->v,u)+getdeg(tvl->v,v) > md[i] ) {
! 885: *up = 0;
! 886: return;
! 887: }
! 888: for ( i = 0, t = 0; i <= j; i++ ) {
! 889: mulp(vl,uh[i],vh[j-i],&s); addp(vl,s,t,&w); t = w;
! 890: }
! 891: /* s = degree j part of (f-uv) */
! 892: exthpc(vl,vx,ff,j,&fj); subp(vl,fj,t,&s);
! 893: for ( i = 0, wu = 0, wv = 0; i <= n; i++ ) {
! 894: if ( s )
! 895: si = 0;
! 896: else if ( VR(s) == vx )
! 897: coefp(s,i,&si);
! 898: else if ( i == 0 )
! 899: si = s;
! 900: else
! 901: si = 0;
! 902: if ( si ) {
! 903: mulp(vl,si,cu[i],&m); addp(vl,wu,m,&t); wu = t;
! 904: mulp(vl,si,cv[i],&m); addp(vl,wv,m,&t); wv = t;
! 905: }
! 906: }
! 907: if ( !wu ) {
! 908: gfsn_poly_to_poly(vl,u,vy,&t); u = t;
! 909: if ( divtp(vl,f,u,&q) ) {
! 910: cont_pp_mv_sf(vl,onevl,u,&cont,up);
! 911: return;
! 912: }
! 913: }
! 914: if ( !wv ) {
! 915: gfsn_poly_to_poly(vl,v,vy,&t); v = t;
! 916: if ( divtp(vl,f,u,&q) ) {
! 917: cont_pp_mv_sf(vl,onevl,q,&cont,up);
! 918: return;
! 919: }
! 920: }
! 921: addp(vl,u,wu,&t); u = t;
! 922: addp(vl,uh[j],wu,&t); uh[j] = t;
! 923: addp(vl,v,wv,&t); v = t;
! 924: addp(vl,vh[j],wv,&t); vh[j] = t;
! 925: }
! 926: }
! 927:
! 928: void adjust_coef_sf(VL vl,VL rvl,P lcu,P u0,P *r)
! 929: {
! 930: P lcu0,cu;
! 931: DCP dc0,dcu,dc;
! 932:
! 933: substvp_sf(vl,rvl,lcu,0,&lcu0);
! 934: divsp(vl,lcu0,LC(u0),&cu);
! 935: for ( dc0 = 0, dcu = DC(u0); dcu; dcu = NEXT(dcu) ) {
! 936: if ( !dc0 ) {
! 937: NEXTDC(dc0,dc);
! 938: COEF(dc) = lcu;
! 939: } else {
! 940: NEXTDC(dc0,dc);
! 941: mulp(vl,cu,COEF(dcu),&COEF(dc));
! 942: }
! 943: DEG(dc) = DEG(dcu);
! 944: }
! 945: NEXT(dc) = 0;
! 946: MKP(VR(u0),dc0,*r);
! 947: }
! 948:
! 949: void extended_gcd_modyk(P u0,P v0,P *cu,P *cv)
1.6 noro 950: {
951: }
952:
1.7 ! noro 953: void poly_to_gfsn_poly(VL vl,P f,V v,P *r)
1.6 noro 954: {
955: }
956:
1.7 ! noro 957: void gfsn_poly_to_poly(VL vl,P f,V v,P *r)
1.6 noro 958: {
1.1 noro 959: }