![[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.3 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/Fgfs.c,v 1.2 2002/09/26 09:07:42 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.1 noro 9:
10: void gcdsf(VL vl,P *pa,int k,P *r)
11: {
1.3 ! noro 12: P *ps,*pl,*pm;
! 13: P **cp;
1.1 noro 14: int *cn;
15: DCP *ml;
16: Obj obj;
17: int i,j,l,m;
18: P mg,mgsf,t;
19: VL avl,nvl,tvl,svl;
20:
21: ps = (P *)ALLOCA(k*sizeof(P));
22: for ( i = 0, m = 0; i < k; i++ ) {
23: simp_ff((Obj)pa[i],&obj);
24: if ( obj )
1.3 ! noro 25: ps[m++] = (P)obj;
1.1 noro 26: }
27: if ( !m ) {
28: *r = 0;
29: return;
30: }
31: if ( m == 1 ) {
1.3 ! noro 32: *r = ps[0];
1.1 noro 33: return;
34: }
35: pl = (P *)ALLOCA(m*sizeof(P));
36: ml = (DCP *)ALLOCA(m*sizeof(DCP));
37: for ( i = 0; i < m; i++ )
38: monomialfctr(vl,ps[i],&pl[i],&ml[i]);
1.3 ! noro 39: gcdmonomial(vl,ml,m,&mg); simp_ff((Obj)mg,&obj); mgsf = (P)obj;
1.1 noro 40: for ( i = 0, nvl = vl, avl = 0; nvl && i < m; i++ ) {
41: clctv(vl,pl[i],&tvl);
42: intersectv(nvl,tvl,&svl); nvl = svl;
43: mergev(vl,avl,tvl,&svl); avl = svl;
44: }
45: if ( !nvl ) {
46: *r = mgsf;
47: return;
48: }
49: if ( !NEXT(avl) ) {
50: ugcdsf(pl,m,&t);
51: mulp(vl,mgsf,t,r);
52: return;
53: }
54: for ( tvl = nvl, i = 0; tvl; tvl = NEXT(tvl), i++ );
55: for ( tvl = avl, j = 0; tvl; tvl = NEXT(tvl), j++ );
56: if ( i == j ) {
57: /* all the pl[i]'s have the same variables */
58: gcdsf_main(avl,pl,m,&t);
59: } else {
60: cp = (P **)ALLOCA(m*sizeof(P *));
61: cn = (int *)ALLOCA(m*sizeof(int));
62: for ( i = 0; i < m; i++ ) {
63: cp[i] = (P *)ALLOCA(lengthp(pl[i])*sizeof(P));
64: cn[i] = pcoef(vl,nvl,pl[i],cp[i]);
65: }
66: for ( i = j = 0; i < m; i++ )
67: j += cn[i];
68: pm = (P *)ALLOCA(j*sizeof(P));
69: for ( i = l = 0; i < m; i++ )
70: for ( j = 0; j < cn[i]; j++ )
71: pm[l++] = cp[i][j];
72: gcdsf(vl,pm,l,&t);
73: }
74: mulp(vl,mgsf,t,r);
75: }
76:
77: /* univariate gcd */
78:
79: void ugcdsf(P *pa,int m,P *r)
80: {
1.3 ! noro 81: P *ps;
1.1 noro 82: int i;
83: UM w1,w2,w3,w;
84: int d;
85: V v;
86:
87: if ( m == 1 ) {
88: *r = pa[0];
89: return;
90: }
1.3 ! noro 91: for ( i = 0; i < m; i++ )
! 92: if ( NUM(pa[i]) ) {
! 93: itogfs(1,r);
! 94: return;
! 95: }
1.1 noro 96: ps = (P *)ALLOCA(m*sizeof(P));
97: sort_by_deg(m,pa,ps);
1.3 ! noro 98: v = VR(ps[m-1]);
! 99: d = getdeg(v,ps[m-1]);
1.1 noro 100: w1 = W_UMALLOC(d);
101: w2 = W_UMALLOC(d);
102: w3 = W_UMALLOC(d);
103: ptosfum(ps[0],w1);
104: for ( i = 1; i < m; i++ ) {
105: ptosfum(ps[i],w2);
106: gcdsfum(w1,w2,w3);
107: w = w1; w1 = w3; w3 = w;
108: if ( !DEG(w1) ) {
1.3 ! noro 109: itogfs(1,r);
1.1 noro 110: return;
111: }
112: }
113: sfumtop(v,w1,r);
114: }
115:
116:
117: /* all the pa[i]'s have the same variables (=vl) */
118:
119: void gcdsf_main(VL vl,P *pa,int m,P *r)
120: {
1.3 ! noro 121: int nv,i,i0,imin,d,d0,d1,d2,dmin,index;
! 122: V v,v0,vmin;
1.2 noro 123: VL tvl,nvl,rvl,nvl0,rvl0;
1.3 ! noro 124: P *pc, *ps, *ph,*lps;
! 125: P x,t,cont,hg,g,hm,mod,s;
! 126: P hge,ge,ce,he,u,cof1e,mode,mod1,adj,cof1,coadj,q;
! 127: GFS sf;
1.2 noro 128:
1.1 noro 129: for ( nv = 0, tvl = vl; tvl; tvl = NEXT(tvl), nv++);
130: if ( nv == 1 ) {
131: ugcdsf(pa,m,r);
132: return;
133: }
134: /* find v s.t. min(deg(pa[i],v)) is minimal */
135: tvl = vl;
136: do {
137: v = tvl->v;
138: i = 0;
139: do {
1.3 ! noro 140: d = getdeg(v,pa[i]);
1.1 noro 141: if ( i == 0 || (d < d0) ) {
142: d0 = d; i0 = i; v0 = v;
143: }
144: } while ( ++i < m );
145: if ( tvl == vl || (d0 < dmin) ) {
146: dmin = d0; imin = i0; vmin = v0;
147: }
148: } while ( tvl = NEXT(tvl) );
149:
150: /* reorder variables so that vmin is the last variable */
151: for ( nvl0 = 0, rvl0 = 0, tvl = vl; tvl; tvl = NEXT(tvl) )
152: if ( tvl->v != vmin ) {
153: NEXTVL(nvl0,nvl); nvl->v = tvl->v;
154: NEXTVL(rvl0,rvl); rvl->v = tvl->v;
155: }
156: /* rvl = remaining variables */
1.3 ! noro 157: NEXT(rvl) = 0; rvl = rvl0;
1.1 noro 158: /* nvl = ...,vmin */
1.3 ! noro 159: NEXTVL(nvl0,nvl); nvl->v = vmin; NEXT(nvl) = 0; nvl = nvl0;
1.2 noro 160: MKV(vmin,x);
1.1 noro 161:
162: /* for content and primitive part */
163: pc = (P *)ALLOCA(m*sizeof(P));
164: ps = (P *)ALLOCA(m*sizeof(P));
165: ph = (P *)ALLOCA(m*sizeof(P));
166: /* separate the contents */
167: for ( i = 0; i < m; i++ ) {
1.3 ! noro 168: reorderp(nvl,vl,pa[i],&t);
1.1 noro 169: cont_pp_mv_sf(nvl,rvl,t,&pc[i],&ps[i]);
170: head_monomial(vmin,ps[i],&ph[i],&t);
171: }
172: ugcdsf(pc,m,&cont);
173: ugcdsf(ph,m,&hg);
174:
175: /* for hg*pp (used in check phase) */
176: lps = (P *)ALLOCA(m*sizeof(P));
177: for ( i = 0; i < m; i++ )
178: mulp(nvl,hg,ps[i],&lps[i]);
179:
180: while ( 1 ) {
181: g = 0;
1.3 ! noro 182: cof1 = 0;
1.1 noro 183: hm = 0;
1.3 ! noro 184: itogfs(1,&mod);
1.1 noro 185: index = 0;
1.3 ! noro 186: for ( index = 0; getdeg(vmin,mod) <= d+1; index++ ) {
1.1 noro 187: /* evaluation pt */
1.3 ! noro 188: indextogfs(index,&s);
1.1 noro 189: substp(nvl,hg,vmin,s,&hge);
190: if ( !hge )
191: continue;
192: for ( i = 0; i < m; i++ )
193: substp(nvl,ps[i],vmin,s,&ph[i]);
194: /* ge = GCD(ps[0]|x=s,...,ps[m-1]|x=s) */
195: gcdsf(nvl,ph,m,&ge);
196: head_monomial(vmin,ge,&ce,&he);
1.3 ! noro 197: if ( NUM(he) ) {
1.1 noro 198: *r = cont;
199: return;
200: }
1.3 ! noro 201: divgfs((GFS)hge,(GFS)ce,&sf); t = (P)sf;
! 202: mulp(nvl,t,ge,&u); ge = u;
1.1 noro 203: divsp(nvl,ph[imin],ge,&t); mulp(nvl,hge,t,&cof1e);
1.2 noro 204: /* hm=0 : reset; he==hm : lucky */
1.3 ! noro 205: if ( !hm || !compp(nvl,he,hm) ) {
1.2 noro 206: substp(nvl,mod,vmin,s,&mode); divsp(nvl,mod,mode,&mod1);
207: /* adj = mod/(mod|x=s)*(ge-g|x=s) */
208: substp(nvl,g,vmin,s,&t);
209: subp(nvl,ge,t,&u); mulp(nvl,mod1,u,&adj);
210: /* coadj = mod/(mod|vmin=s)*(cof1e-cof1e|vmin=s) */
211: substp(nvl,cof1,vmin,s,&t);
1.3 ! noro 212: subp(nvl,cof1e,t,&u); mulp(nvl,mod1,u,&coadj);
1.2 noro 213: if ( !adj ) {
214: /* adj == gcd ? */
215: for ( i = 0; i < m; i++ )
1.3 ! noro 216: if ( !divtp(nvl,lps[i],g,&t) )
1.2 noro 217: break;
218: if ( i == m ) {
1.3 ! noro 219: cont_pp_mv_sf(nvl,rvl,g,&t,&u);
1.2 noro 220: mulp(nvl,cont,u,&t);
1.3 ! noro 221: reorderp(vl,nvl,t,r);
1.2 noro 222: return;
223: }
224: } else if ( !coadj ) {
1.3 ! noro 225: /* ps[imin]/coadj == gcd ? */
! 226: if ( divtp(nvl,lps[imin],cof1,&q) ) {
1.2 noro 227: for ( i = 0; i < m; i++ )
228: if ( !divtp(nvl,lps[i],q,&t) )
229: break;
230: if ( i == m ) {
231: cont_pp_mv_sf(nvl,rvl,q,&t,&u);
232: mulp(nvl,cont,u,&t);
1.3 ! noro 233: reorderp(vl,nvl,t,r);
1.2 noro 234: return;
235: }
236: }
237: }
238: addp(nvl,g,adj,&t); g = t;
239: addp(nvl,cof1,coadj,&t); cof1 = t;
240: subp(nvl,x,s,&t); mulp(nvl,mod,t,&u); mod = u;
241: hm = he;
242: } else {
243: d1 = homdeg(hm); d2 = homdeg(he);
244: if ( d1 < d2 ) /* we use current hm */
245: continue;
246: else if ( d1 > d2 ) {
247: /* use he */
248: g = ge;
249: cof1 = cof1e;
250: hm = he;
251: subp(nvl,x,s,&mod);
252: } else {
253: /* d1==d2, but hm!=he => both are unlucky */
254: g = 0;
255: cof1 = 0;
1.3 ! noro 256: itogfs(1,&mod);
1.2 noro 257: }
1.1 noro 258: }
259: }
260: }
261: }
262:
263: void head_monomial(V v,P p,P *coef,P *term)
264: {
265: P t,s,u;
266: DCP dc;
267: GFS one;
1.3 ! noro 268: VL vl;
1.1 noro 269:
1.3 ! noro 270: itogfs(1,&one);
! 271: t = (P)one;
1.1 noro 272: while ( 1 ) {
273: if ( NUM(p) || VR(p) == v ) {
274: *coef = p;
275: *term = t;
276: return;
277: } else {
1.3 ! noro 278: NEWDC(dc);
! 279: COEF(dc) = (P)one; DEG(dc) = DEG(DC(p));
1.1 noro 280: MKP(VR(p),dc,s);
281: mulp(vl,t,s,&u); t = u;
282: p = COEF(DC(p));
283: }
284: }
285: }
286:
287: void cont_pp_mv_sf(VL vl,VL rvl,P p,P *c,P *pp)
288: {
289: DP dp;
290: MP t;
291: int i,m;
292: P *ps;
293:
294: ptod(vl,rvl,p,&dp);
295: for ( t = BDY(dp), m = 0; t; t = NEXT(t), m++ );
296: ps = (P *)ALLOCA(m*sizeof(P));
1.3 ! noro 297: for ( t = BDY(dp), i = 0; t; t = NEXT(t), i++ )
1.1 noro 298: ps[i] = C(t);
299: ugcdsf(ps,m,c);
1.3 ! noro 300: divsp(vl,p,*c,pp);
1.1 noro 301: }