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