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