![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to NEZ.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 NEZ.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/NEZ.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3:
4: void nezgcdnpz(vl,ps,m,pr)
5: VL vl;
6: P *ps,*pr;
7: int m;
8: {
9: P t,s,mg;
10: VL tvl,svl,avl,nvl;
11: int i,j,k;
12: N c;
13: Q cq;
14: P *pl,*pm;
15: DCP *ml;
16: Q *cl;
17: P **cp;
18: int *cn;
19:
20: pl = (P *)ALLOCA(m*sizeof(P));
21: ml = (DCP *)ALLOCA(m*sizeof(DCP));
22: cl = (Q *)ALLOCA(m*sizeof(Q));
23: for ( i = 0; i < m; i++ )
24: monomialfctr(vl,ps[i],&pl[i],&ml[i]);
25: gcdmonomial(vl,ml,m,&mg);
26: for ( i = 0; i < m; i++ ) {
27: ptozp(pl[i],1,&cl[i],&t); pl[i] = t;
28: }
29: for ( i = 1, cq = cl[0]; i < m; i++ ) {
30: gcdn(NM(cl[i]),NM(cq),&c); NTOQ(c,1,cq);
31: }
32: for ( i = 0; i < m; i++ )
33: if ( NUM(pl[i]) ) {
34: mulp(vl,(P)cq,mg,pr); return;
35: }
36: for ( i = 0, nvl = vl, avl = 0; nvl && i < m; i++ ) {
37: clctv(vl,pl[i],&tvl);
38: intersectv(nvl,tvl,&svl); nvl = svl;
39: mergev(vl,avl,tvl,&svl); avl = svl;
40: }
41: if ( !nvl ) {
42: mulp(vl,(P)cq,mg,pr); return;
43: }
44: if ( !NEXT(avl) ) {
45: nuezgcdnpzmain(vl,pl,m,&t); mulp(vl,t,(P)cq,&s); mulp(vl,s,mg,pr);
46: return;
47: }
48: for ( tvl = nvl, i = 0; tvl; tvl = NEXT(tvl), i++ );
49: for ( tvl = avl, j = 0; tvl; tvl = NEXT(tvl), j++ );
50: if ( i == j ) {
51: /* all the pl[i]'s have the same variables */
52: sortplistbyhomdeg(pl,m); nezgcdnpzmain(nvl,pl,m,&t);
53: #if 0
54: /* search the minimal degree poly */
55: for ( i = 0; i < m; i++ ) {
56: for ( tvl = nvl; tvl; tvl = NEXT(tvl) ) {
57: dt = getdeg(tvl->v,pl[i]);
58: if ( tvl == nvl || dt < d1 ) {
59: v1 = tvl->v; d1 = dt;
60: }
61: }
62: if ( i == 0 || d1 < d ) {
63: v = v1; d = d1; j = i;
64: }
65: }
66: t = pl[0]; pl[0] = pl[j]; pl[j] = t;
67: if ( v != nvl->v ) {
68: reordvar(nvl,v,&mvl);
69: for ( i = 0; i < m; i++ ) {
70: reorderp(mvl,nvl,pl[i],&t); pl[i] = t;
71: }
72: nezgcdnpzmain(mvl,pl,m,&s); reorderp(nvl,mvl,s,&t);
73: } else
74: nezgcdnpzmain(nvl,pl,m,&t);
75: #endif
76: } else {
77: cp = (P **)ALLOCA(m*sizeof(P *));
78: cn = (int *)ALLOCA(m*sizeof(int));
79: for ( i = 0; i < m; i++ ) {
80: cp[i] = (P *)ALLOCA(lengthp(pl[i])*sizeof(P));
81: cn[i] = pcoef(vl,nvl,pl[i],cp[i]);
82: }
83: for ( i = j = 0; i < m; i++ )
84: j += cn[i];
85: pm = (P *)ALLOCA(j*sizeof(P));
86: for ( i = k = 0; i < m; i++ )
87: for ( j = 0; j < cn[i]; j++ )
88: pm[k++] = cp[i][j];
89: nezgcdnpz(vl,pm,k,&t);
90: }
91: mulp(vl,t,(P)cq,&s); mulp(vl,s,mg,pr);
92: }
93:
94: void sortplistbyhomdeg(p,n)
95: P *p;
96: int n;
97: {
98: int i,j,k;
99: int *l;
100: P t;
101:
102: l = (int *)ALLOCA(n*sizeof(int));
103: for ( i = 0; i < n; i++ )
104: l[i] = homdeg(p[i]);
105: for ( i = 0; i < n; i++ )
106: for ( j = i + 1; j < n; j++ )
107: if ( l[j] < l[i] ) {
108: t = p[i]; p[i] = p[j]; p[j] = t;
109: k = l[i]; l[i] = l[j]; l[j] = k;
110: }
111: }
112:
113: void nuezgcdnpzmain(vl,ps,m,r)
114: VL vl;
115: P *ps;
116: int m;
117: P *r;
118: {
119: P *tps;
120: P f,t;
121: int i;
122:
123: for ( i = 0; i < m; i++ )
124: if ( NUM(ps[i]) ) {
125: *r = (P)ONE; return;
126: }
127: tps = (P *) ALLOCA(m*sizeof(P));
128: for ( i = 0; i < m; i++ )
129: tps[i] = ps[i];
130: sortplist(tps,m);
131: for ( i = 1, f = tps[0]; i < m && !NUM(f); i++ ) {
132: uezgcdpz(vl,f,tps[i],&t); f = t;
133: }
134: *r = f;
135: }
136:
137: void gcdmonomial(vl,dcl,m,r)
138: VL vl;
139: DCP *dcl;
140: int m;
141: P *r;
142: {
143: int i,j,n;
144: P g,x,s,t;
145: Q d;
146: DCP dc;
147: VN vn;
148:
149: for ( i = 0; i < m; i++ )
150: if ( !dcl[i] ) {
151: *r = (P)ONE; return;
152: }
153: for ( n = 0, dc = dcl[0]; dc; dc = NEXT(dc), n++ );
154: vn = (VN)ALLOCA(n*sizeof(struct oVN));
155: for ( i = 0, dc = dcl[0]; i < n; dc = NEXT(dc), i++ ) {
156: vn[i].v = VR(COEF(dc)); vn[i].n = QTOS(DEG(dc));
157: }
158: for ( i = 1; i < m; i++ ) {
159: for ( j = 0; j < n; j++ ) {
160: for ( dc = dcl[i]; dc; dc = NEXT(dc) )
161: if ( VR(COEF(dc)) == vn[j].v ) {
162: vn[j].n = MIN(vn[j].n,QTOS(DEG(dc))); break;
163: }
164: if ( !dc )
165: vn[j].n = 0;
166: }
167: for ( j = n-1; j >= 0 && !vn[j].n; j-- );
168: if ( j < 0 ) {
169: *r = (P)ONE; return;
170: } else
171: n = j+1;
172: }
173: for ( j = 0, g = (P)ONE; j < n; j++ )
174: if ( vn[j].n ) {
175: MKV(vn[j].v,x); STOQ(vn[j].n,d);
176: pwrp(vl,x,d,&t); mulp(vl,t,g,&s); g = s;
177: }
178: *r = g;
179: }
180:
181: void nezgcdnpzmain(vl,pl,m,pr)
182: VL vl;
183: P *pl,*pr;
184: int m;
185: {
186: P *ppl,*pcl;
187: int i;
188: P cont,gcd,t,s;
189: DCP dc;
190:
191: ppl = (P *)ALLOCA(m*sizeof(P));
192: pcl = (P *)ALLOCA(m*sizeof(P));
193: for ( i = 0; i < m; i++ )
194: pcp(vl,pl[i],&ppl[i],&pcl[i]);
195: nezgcdnpz(vl,pcl,m,&cont);
196: sqfrp(vl,ppl[0],&dc);
197: for ( dc = NEXT(dc), gcd = (P)ONE; dc; dc = NEXT(dc) ) {
198: if ( NUM(COEF(dc)) )
199: continue;
200: nezgcdnpp(vl,dc,ppl+1,m-1,&t);
201: if ( NUM(t) )
202: continue;
203: mulp(vl,gcd,t,&s); gcd = s;
204: }
205: mulp(vl,gcd,cont,pr);
206: }
207:
208: void nezgcdnpp(vl,dc,pl,m,r)
209: VL vl;
210: DCP dc;
211: P *pl;
212: int m;
213: P *r;
214: {
215: int i,k;
216: P g,t,s,gcd;
217: P *pm;
218:
219: nezgcdnp_sqfr_primitive(vl,COEF(dc),pl,m,&gcd);
220: if ( NUM(gcd) ) {
221: *r = (P)ONE; return;
222: }
223: pm = (P *) ALLOCA(m*sizeof(P));
224: for ( i = 0; i < m; i++ ) {
225: divsp(vl,pl[i],gcd,&pm[i]);
226: if ( NUM(pm[i]) ) {
227: *r = gcd; return;
228: }
229: }
230: for ( g = gcd, k = QTOS(DEG(dc))-1; k > 0; k-- ) {
231: nezgcdnp_sqfr_primitive(vl,g,pm,m,&t);
232: if ( NUM(t) )
233: break;
234: mulp(vl,gcd,t,&s); gcd = s;
235: for ( i = 0; i < m; i++ ) {
236: divsp(vl,pm[i],t,&s);
237: if ( NUM(s) ) {
238: *r = gcd; return;
239: }
240: pm[i] = s;
241: }
242: }
243: *r = gcd;
244: }
245:
246: /*
247: * pr = gcd(p0,ps[0],...,ps[m-1])
248: *
249: * p0 is already square-free and primitive.
250: * ps[i] is at least primitive.
251: *
252: */
253:
254: void nezgcdnp_sqfr_primitive(vl,p0,ps,m,pr)
255: VL vl;
256: int m;
257: P p0,*ps,*pr;
258: {
259: /* variables */
260: P p00,g,h,g0,h0,a0,b0;
261: P lp0,lgp0,lp00,lg,lg0,cbd,tq,t;
262: P *lc;
263: P *tps;
264: VL nvl1,nvl2,nvl,tvl;
265: V v;
266: int i,j,k,d0,dg,dg0,dmin;
267: VN vn0,vn1,vnt;
268: int nv,flag;
269:
270: /* set-up */
271: if ( NUM(p0) ) {
272: *pr = (P) ONE; return;
273: }
274: for ( v = VR(p0), i = 0; i < m; i++ )
275: if ( NUM(ps[i]) || (v != VR(ps[i])) ) {
276: *pr = (P)ONE; return;
277: }
278: tps = (P *) ALLOCA(m*sizeof(P));
279: for ( i = 0; i < m; i++ )
280: tps[i] = ps[i];
281: sortplist(tps,m);
282: /* deg(tps[0]) <= deg(tps[1]) <= ... */
283:
284: if ( !cmpq(DEG(DC(p0)),ONE) ) {
285: if ( divcheck(vl,tps,m,(P)ONE,p0) )
286: *pr = p0;
287: else
288: *pr = (P)ONE;
289: return;
290: }
291:
292: lp0 = LC(p0); dmin = d0 = deg(v,p0); lc = (P *)ALLOCA((m+1)*sizeof(P));
293: for ( lc[0] = lp0, i = 0; i < m; i++ )
294: lc[i+1] = LC(tps[i]);
295: clctv(vl,p0,&nvl);
296: for ( i = 0; i < m; i++ ) {
297: clctv(vl,tps[i],&nvl1); mergev(vl,nvl,nvl1,&nvl2); nvl = nvl2;
298: }
299: nezgcdnpz(nvl,lc,m+1,&lg);
300:
301: mulp(nvl,p0,lg,&lgp0); k = dbound(v,lgp0)+1; cbound(nvl,lgp0,(Q *)&cbd);
302: for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++ );
303: W_CALLOC(nv,struct oVN,vn0); W_CALLOC(nv,struct oVN,vnt);
304: W_CALLOC(nv,struct oVN,vn1);
305: for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
306: vn1[i].v = vn0[i].v = tvl->v;
307:
308: /* main loop */
309: for ( dg = deg(v,tps[0]) + 1; ; next(vn0) )
310: do {
311: for ( i = 0, j = 0; vn0[i].v; i++ )
312: if ( vn0[i].n ) vnt[j++].v = (V)i;
313: vnt[j].n = 0;
314:
315: /* find b s.t. LC(p0)(b), LC(tps[i])(b) != 0 */
316: mulsgn(vn0,vnt,j,vn1); substvp(nvl,p0,vn1,&p00);
317: flag = (!zerovpchk(nvl,lp0,vn1) && sqfrchk(p00));
318: for ( i = 0; flag && (i < m); i++ )
319: flag &= (!zerovpchk(nvl,LC(tps[i]),vn1));
320: if ( !flag )
321: continue;
322:
323: /* substitute y -> b */
324: substvp(nvl,lg,vn1,&lg0); lp00 = LC(p00);
325: /* extended-gcd in 1-variable */
326: uexgcdnp(nvl,p00,tps,m,vn1,(Q)cbd,&g0,&h0,&a0,&b0,(Q *)&tq);
327: if ( NUM(g0) ) {
328: *pr = (P)ONE; return;
329: } else if ( dg > ( dg0 = deg(v,g0) ) ) {
330: dg = dg0;
331: if ( dg == d0 ) {
332: if ( divcheck(nvl,tps,m,lp0,p0) ) {
333: *pr = p0; return;
334: }
335: } else if ( dg == deg(v,tps[0]) ) {
336: if ( divtpz(nvl,p0,tps[0],&t) &&
337: divcheck(nvl,tps+1,m-1,LC(tps[0]),tps[0]) ) {
338: *pr = tps[0]; return;
339: } else
340: break;
341: } else {
342: henmv(nvl,vn1,lgp0,g0,h0,a0,b0,lg,lp0,lg0,lp00,(Q)tq,k,&g,&h);
343: if ( divtpz(nvl,lgp0,g,&t) &&
344: divcheck(nvl,tps,m,lg,g) ) {
345: pcp(nvl,g,pr,&t); return;
346: }
347: }
348: }
349: } while ( !nextbin(vnt,j) );
350: }
351:
352: void intersectv(vl1,vl2,vlp)
353: VL vl1,vl2,*vlp;
354: {
355: int i,n;
356: VL tvl;
357: VN tvn;
358:
359: if ( !vl1 || !vl2 ) {
360: *vlp = 0; return;
361: }
362: for ( n = 0, tvl = vl1; tvl; tvl = NEXT(tvl), n++ );
363: tvn = (VN) ALLOCA(n*sizeof(struct oVN));
364: for ( i = 0, tvl = vl1; i < n; tvl = NEXT(tvl), i++ ) {
365: tvn[i].v = tvl->v; tvn[i].n = 0;
366: }
367: for ( tvl = vl2; tvl; tvl = NEXT(tvl) )
368: for ( i = 0; i < n; i++ )
369: if ( tvn[i].v == tvl->v ) {
370: tvn[i].n = 1; break;
371: }
372: vntovl(tvn,n,vlp);
373: }
374:
375: int pcoef(vl,ivl,p,cp)
376: VL vl,ivl;
377: P p;
378: P *cp;
379: {
380: VL nvl,tvl,svl,mvl,mvl0;
381: P t;
382:
383: if ( NUM(p) ) {
384: cp[0] = p; return 1;
385: } else {
386: clctv(vl,p,&nvl);
387: for ( tvl = nvl, mvl0 = 0; tvl; tvl = NEXT(tvl) ) {
388: for ( svl = ivl; svl; svl = NEXT(svl) )
389: if ( tvl->v == svl->v )
390: break;
391: if ( !svl ) {
392: if ( !mvl0 ) {
393: NEWVL(mvl0); mvl = mvl0;
394: } else {
395: NEWVL(NEXT(mvl)); mvl = NEXT(mvl);
396: }
397: mvl->v = tvl->v;
398: }
399: }
400: if ( mvl0 )
401: NEXT(mvl) = ivl;
402: else
403: mvl0 = ivl;
404: reorderp(mvl0,nvl,p,&t);
405: return pcoef0(mvl0,ivl,t,cp);
406: }
407: }
408:
409: int pcoef0(vl,ivl,p,cp)
410: VL vl,ivl;
411: P p;
412: P *cp;
413: {
414: int cn,n;
415: DCP dc;
416: V v;
417: VL tvl;
418:
419: if ( NUM(p) ) {
420: cp[0] = p; return 1;
421: } else {
422: for ( v = VR(p), tvl = ivl; tvl; tvl = NEXT(tvl) )
423: if ( v == tvl->v )
424: break;
425: if ( tvl ) {
426: cp[0] = p; return 1;
427: } else {
428: for ( dc = DC(p), n = 0; dc; dc = NEXT(dc) ) {
429: cn = pcoef0(vl,ivl,COEF(dc),cp); cp += cn; n += cn;
430: }
431: return n;
432: }
433: }
434: }
435:
436: int lengthp(p)
437: P p;
438: {
439: int n;
440: DCP dc;
441:
442: if ( NUM(p) )
443: return 1;
444: else {
445: for ( dc = DC(p), n = 0; dc; dc = NEXT(dc) )
446: n += lengthp(COEF(dc));
447: return n;
448: }
449: }