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