![[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.2 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.6 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/NEZ.c,v 1.5 2006/12/03 01:04:10 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51:
52: void nezgcdnpz(vl,ps,m,pr)
53: VL vl;
54: P *ps,*pr;
55: int m;
56: {
57: P t,s,mg;
58: VL tvl,svl,avl,nvl;
59: int i,j,k;
60: N c;
61: Q cq;
62: P *pl,*pm;
63: DCP *ml;
64: Q *cl;
65: P **cp;
66: int *cn;
67:
1.6 ! noro 68: if ( m == 1 ) {
! 69: *pr = ps[0]; return;
! 70: }
1.1 noro 71: pl = (P *)ALLOCA(m*sizeof(P));
72: ml = (DCP *)ALLOCA(m*sizeof(DCP));
73: cl = (Q *)ALLOCA(m*sizeof(Q));
74: for ( i = 0; i < m; i++ )
75: monomialfctr(vl,ps[i],&pl[i],&ml[i]);
76: gcdmonomial(vl,ml,m,&mg);
77: for ( i = 0; i < m; i++ ) {
78: ptozp(pl[i],1,&cl[i],&t); pl[i] = t;
79: }
80: for ( i = 1, cq = cl[0]; i < m; i++ ) {
81: gcdn(NM(cl[i]),NM(cq),&c); NTOQ(c,1,cq);
82: }
83: for ( i = 0; i < m; i++ )
84: if ( NUM(pl[i]) ) {
85: mulp(vl,(P)cq,mg,pr); return;
86: }
87: for ( i = 0, nvl = vl, avl = 0; nvl && i < m; i++ ) {
88: clctv(vl,pl[i],&tvl);
89: intersectv(nvl,tvl,&svl); nvl = svl;
90: mergev(vl,avl,tvl,&svl); avl = svl;
91: }
92: if ( !nvl ) {
93: mulp(vl,(P)cq,mg,pr); return;
94: }
95: if ( !NEXT(avl) ) {
96: nuezgcdnpzmain(vl,pl,m,&t); mulp(vl,t,(P)cq,&s); mulp(vl,s,mg,pr);
97: return;
98: }
99: for ( tvl = nvl, i = 0; tvl; tvl = NEXT(tvl), i++ );
100: for ( tvl = avl, j = 0; tvl; tvl = NEXT(tvl), j++ );
101: if ( i == j ) {
102: /* all the pl[i]'s have the same variables */
103: sortplistbyhomdeg(pl,m); nezgcdnpzmain(nvl,pl,m,&t);
104: #if 0
105: /* search the minimal degree poly */
106: for ( i = 0; i < m; i++ ) {
107: for ( tvl = nvl; tvl; tvl = NEXT(tvl) ) {
108: dt = getdeg(tvl->v,pl[i]);
109: if ( tvl == nvl || dt < d1 ) {
110: v1 = tvl->v; d1 = dt;
111: }
112: }
113: if ( i == 0 || d1 < d ) {
114: v = v1; d = d1; j = i;
115: }
116: }
117: t = pl[0]; pl[0] = pl[j]; pl[j] = t;
118: if ( v != nvl->v ) {
119: reordvar(nvl,v,&mvl);
120: for ( i = 0; i < m; i++ ) {
121: reorderp(mvl,nvl,pl[i],&t); pl[i] = t;
122: }
123: nezgcdnpzmain(mvl,pl,m,&s); reorderp(nvl,mvl,s,&t);
124: } else
125: nezgcdnpzmain(nvl,pl,m,&t);
126: #endif
127: } else {
128: cp = (P **)ALLOCA(m*sizeof(P *));
129: cn = (int *)ALLOCA(m*sizeof(int));
130: for ( i = 0; i < m; i++ ) {
131: cp[i] = (P *)ALLOCA(lengthp(pl[i])*sizeof(P));
132: cn[i] = pcoef(vl,nvl,pl[i],cp[i]);
133: }
134: for ( i = j = 0; i < m; i++ )
135: j += cn[i];
136: pm = (P *)ALLOCA(j*sizeof(P));
137: for ( i = k = 0; i < m; i++ )
138: for ( j = 0; j < cn[i]; j++ )
139: pm[k++] = cp[i][j];
140: nezgcdnpz(vl,pm,k,&t);
141: }
142: mulp(vl,t,(P)cq,&s); mulp(vl,s,mg,pr);
143: }
144:
145: void sortplistbyhomdeg(p,n)
146: P *p;
147: int n;
148: {
149: int i,j,k;
150: int *l;
151: P t;
152:
153: l = (int *)ALLOCA(n*sizeof(int));
154: for ( i = 0; i < n; i++ )
155: l[i] = homdeg(p[i]);
156: for ( i = 0; i < n; i++ )
157: for ( j = i + 1; j < n; j++ )
158: if ( l[j] < l[i] ) {
159: t = p[i]; p[i] = p[j]; p[j] = t;
160: k = l[i]; l[i] = l[j]; l[j] = k;
161: }
162: }
163:
164: void nuezgcdnpzmain(vl,ps,m,r)
165: VL vl;
166: P *ps;
167: int m;
168: P *r;
169: {
170: P *tps;
171: P f,t;
172: int i;
173:
174: for ( i = 0; i < m; i++ )
175: if ( NUM(ps[i]) ) {
176: *r = (P)ONE; return;
177: }
178: tps = (P *) ALLOCA(m*sizeof(P));
179: for ( i = 0; i < m; i++ )
180: tps[i] = ps[i];
181: sortplist(tps,m);
182: for ( i = 1, f = tps[0]; i < m && !NUM(f); i++ ) {
183: uezgcdpz(vl,f,tps[i],&t); f = t;
184: }
185: *r = f;
186: }
187:
188: void gcdmonomial(vl,dcl,m,r)
189: VL vl;
190: DCP *dcl;
191: int m;
192: P *r;
193: {
194: int i,j,n;
195: P g,x,s,t;
196: Q d;
197: DCP dc;
198: VN vn;
199:
200: for ( i = 0; i < m; i++ )
201: if ( !dcl[i] ) {
202: *r = (P)ONE; return;
203: }
204: for ( n = 0, dc = dcl[0]; dc; dc = NEXT(dc), n++ );
205: vn = (VN)ALLOCA(n*sizeof(struct oVN));
206: for ( i = 0, dc = dcl[0]; i < n; dc = NEXT(dc), i++ ) {
207: vn[i].v = VR(COEF(dc)); vn[i].n = QTOS(DEG(dc));
208: }
209: for ( i = 1; i < m; i++ ) {
210: for ( j = 0; j < n; j++ ) {
211: for ( dc = dcl[i]; dc; dc = NEXT(dc) )
212: if ( VR(COEF(dc)) == vn[j].v ) {
213: vn[j].n = MIN(vn[j].n,QTOS(DEG(dc))); break;
214: }
215: if ( !dc )
216: vn[j].n = 0;
217: }
218: for ( j = n-1; j >= 0 && !vn[j].n; j-- );
219: if ( j < 0 ) {
220: *r = (P)ONE; return;
221: } else
222: n = j+1;
223: }
224: for ( j = 0, g = (P)ONE; j < n; j++ )
225: if ( vn[j].n ) {
226: MKV(vn[j].v,x); STOQ(vn[j].n,d);
227: pwrp(vl,x,d,&t); mulp(vl,t,g,&s); g = s;
228: }
229: *r = g;
230: }
231:
232: void nezgcdnpzmain(vl,pl,m,pr)
233: VL vl;
234: P *pl,*pr;
235: int m;
236: {
237: P *ppl,*pcl;
238: int i;
239: P cont,gcd,t,s;
240: DCP dc;
241:
242: ppl = (P *)ALLOCA(m*sizeof(P));
243: pcl = (P *)ALLOCA(m*sizeof(P));
244: for ( i = 0; i < m; i++ )
245: pcp(vl,pl[i],&ppl[i],&pcl[i]);
246: nezgcdnpz(vl,pcl,m,&cont);
247: sqfrp(vl,ppl[0],&dc);
248: for ( dc = NEXT(dc), gcd = (P)ONE; dc; dc = NEXT(dc) ) {
249: if ( NUM(COEF(dc)) )
250: continue;
251: nezgcdnpp(vl,dc,ppl+1,m-1,&t);
252: if ( NUM(t) )
253: continue;
254: mulp(vl,gcd,t,&s); gcd = s;
255: }
256: mulp(vl,gcd,cont,pr);
257: }
258:
259: void nezgcdnpp(vl,dc,pl,m,r)
260: VL vl;
261: DCP dc;
262: P *pl;
263: int m;
264: P *r;
265: {
266: int i,k;
267: P g,t,s,gcd;
268: P *pm;
269:
270: nezgcdnp_sqfr_primitive(vl,COEF(dc),pl,m,&gcd);
271: if ( NUM(gcd) ) {
272: *r = (P)ONE; return;
273: }
274: pm = (P *) ALLOCA(m*sizeof(P));
275: for ( i = 0; i < m; i++ ) {
276: divsp(vl,pl[i],gcd,&pm[i]);
277: if ( NUM(pm[i]) ) {
278: *r = gcd; return;
279: }
280: }
281: for ( g = gcd, k = QTOS(DEG(dc))-1; k > 0; k-- ) {
282: nezgcdnp_sqfr_primitive(vl,g,pm,m,&t);
283: if ( NUM(t) )
284: break;
285: mulp(vl,gcd,t,&s); gcd = s;
286: for ( i = 0; i < m; i++ ) {
287: divsp(vl,pm[i],t,&s);
288: if ( NUM(s) ) {
289: *r = gcd; return;
290: }
291: pm[i] = s;
292: }
293: }
294: *r = gcd;
295: }
296:
297: /*
298: * pr = gcd(p0,ps[0],...,ps[m-1])
299: *
300: * p0 is already square-free and primitive.
301: * ps[i] is at least primitive.
302: *
303: */
304:
305: void nezgcdnp_sqfr_primitive(vl,p0,ps,m,pr)
306: VL vl;
307: int m;
308: P p0,*ps,*pr;
309: {
310: /* variables */
311: P p00,g,h,g0,h0,a0,b0;
312: P lp0,lgp0,lp00,lg,lg0,cbd,tq,t;
313: P *lc;
314: P *tps;
315: VL nvl1,nvl2,nvl,tvl;
316: V v;
1.5 noro 317: int i,j,k,d0,dg,dg0,dmin,z;
318: VN vn1;
319: int nv,flag,max;
1.1 noro 320:
321: /* set-up */
322: if ( NUM(p0) ) {
323: *pr = (P) ONE; return;
324: }
325: for ( v = VR(p0), i = 0; i < m; i++ )
326: if ( NUM(ps[i]) || (v != VR(ps[i])) ) {
327: *pr = (P)ONE; return;
328: }
329: tps = (P *) ALLOCA(m*sizeof(P));
330: for ( i = 0; i < m; i++ )
331: tps[i] = ps[i];
332: sortplist(tps,m);
333: /* deg(tps[0]) <= deg(tps[1]) <= ... */
334:
335: if ( !cmpq(DEG(DC(p0)),ONE) ) {
336: if ( divcheck(vl,tps,m,(P)ONE,p0) )
337: *pr = p0;
338: else
339: *pr = (P)ONE;
340: return;
341: }
342:
343: lp0 = LC(p0); dmin = d0 = deg(v,p0); lc = (P *)ALLOCA((m+1)*sizeof(P));
344: for ( lc[0] = lp0, i = 0; i < m; i++ )
345: lc[i+1] = LC(tps[i]);
346: clctv(vl,p0,&nvl);
347: for ( i = 0; i < m; i++ ) {
348: clctv(vl,tps[i],&nvl1); mergev(vl,nvl,nvl1,&nvl2); nvl = nvl2;
349: }
350: nezgcdnpz(nvl,lc,m+1,&lg);
351:
352: mulp(nvl,p0,lg,&lgp0); k = dbound(v,lgp0)+1; cbound(nvl,lgp0,(Q *)&cbd);
353: for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++ );
354: W_CALLOC(nv,struct oVN,vn1);
355: for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
1.5 noro 356: vn1[i].v = tvl->v;
1.1 noro 357:
358: /* main loop */
1.5 noro 359: /* finally, 'max' random evaluations will be generated */
360: for ( max = 1, dg = deg(v,tps[0]) + 1; ; max = 2*max )
361: for ( z = 0; z < max; z++ ) {
362: for ( i = 0; vn1[i].v; i++ )
363: vn1[i].n = mt_genrand()%max;
1.1 noro 364:
365: /* find b s.t. LC(p0)(b), LC(tps[i])(b) != 0 */
1.5 noro 366:
367: substvp(nvl,p0,vn1,&p00);
1.1 noro 368: flag = (!zerovpchk(nvl,lp0,vn1) && sqfrchk(p00));
369: for ( i = 0; flag && (i < m); i++ )
370: flag &= (!zerovpchk(nvl,LC(tps[i]),vn1));
371: if ( !flag )
372: continue;
373:
374: /* substitute y -> b */
375: substvp(nvl,lg,vn1,&lg0); lp00 = LC(p00);
376: /* extended-gcd in 1-variable */
377: uexgcdnp(nvl,p00,tps,m,vn1,(Q)cbd,&g0,&h0,&a0,&b0,(Q *)&tq);
378: if ( NUM(g0) ) {
379: *pr = (P)ONE; return;
380: } else if ( dg > ( dg0 = deg(v,g0) ) ) {
381: dg = dg0;
382: if ( dg == d0 ) {
383: if ( divcheck(nvl,tps,m,lp0,p0) ) {
384: *pr = p0; return;
385: }
386: } else if ( dg == deg(v,tps[0]) ) {
387: if ( divtpz(nvl,p0,tps[0],&t) &&
388: divcheck(nvl,tps+1,m-1,LC(tps[0]),tps[0]) ) {
389: *pr = tps[0]; return;
390: } else
391: break;
392: } else {
393: henmv(nvl,vn1,lgp0,g0,h0,a0,b0,lg,lp0,lg0,lp00,(Q)tq,k,&g,&h);
394: if ( divtpz(nvl,lgp0,g,&t) &&
395: divcheck(nvl,tps,m,lg,g) ) {
396: pcp(nvl,g,pr,&t); return;
397: }
398: }
399: }
1.5 noro 400: }
1.1 noro 401: }
402:
403: void intersectv(vl1,vl2,vlp)
404: VL vl1,vl2,*vlp;
405: {
406: int i,n;
407: VL tvl;
408: VN tvn;
409:
410: if ( !vl1 || !vl2 ) {
411: *vlp = 0; return;
412: }
413: for ( n = 0, tvl = vl1; tvl; tvl = NEXT(tvl), n++ );
414: tvn = (VN) ALLOCA(n*sizeof(struct oVN));
415: for ( i = 0, tvl = vl1; i < n; tvl = NEXT(tvl), i++ ) {
416: tvn[i].v = tvl->v; tvn[i].n = 0;
417: }
418: for ( tvl = vl2; tvl; tvl = NEXT(tvl) )
419: for ( i = 0; i < n; i++ )
420: if ( tvn[i].v == tvl->v ) {
421: tvn[i].n = 1; break;
422: }
423: vntovl(tvn,n,vlp);
424: }
425:
426: int pcoef(vl,ivl,p,cp)
427: VL vl,ivl;
428: P p;
429: P *cp;
430: {
431: VL nvl,tvl,svl,mvl,mvl0;
432: P t;
433:
434: if ( NUM(p) ) {
435: cp[0] = p; return 1;
436: } else {
437: clctv(vl,p,&nvl);
438: for ( tvl = nvl, mvl0 = 0; tvl; tvl = NEXT(tvl) ) {
439: for ( svl = ivl; svl; svl = NEXT(svl) )
440: if ( tvl->v == svl->v )
441: break;
442: if ( !svl ) {
443: if ( !mvl0 ) {
444: NEWVL(mvl0); mvl = mvl0;
445: } else {
446: NEWVL(NEXT(mvl)); mvl = NEXT(mvl);
447: }
448: mvl->v = tvl->v;
449: }
450: }
451: if ( mvl0 )
452: NEXT(mvl) = ivl;
453: else
454: mvl0 = ivl;
455: reorderp(mvl0,nvl,p,&t);
456: return pcoef0(mvl0,ivl,t,cp);
457: }
458: }
459:
460: int pcoef0(vl,ivl,p,cp)
461: VL vl,ivl;
462: P p;
463: P *cp;
464: {
465: int cn,n;
466: DCP dc;
467: V v;
468: VL tvl;
469:
470: if ( NUM(p) ) {
471: cp[0] = p; return 1;
472: } else {
473: for ( v = VR(p), tvl = ivl; tvl; tvl = NEXT(tvl) )
474: if ( v == tvl->v )
475: break;
476: if ( tvl ) {
477: cp[0] = p; return 1;
478: } else {
479: for ( dc = DC(p), n = 0; dc; dc = NEXT(dc) ) {
480: cn = pcoef0(vl,ivl,COEF(dc),cp); cp += cn; n += cn;
481: }
482: return n;
483: }
484: }
485: }
486:
487: int lengthp(p)
488: P p;
489: {
490: int n;
491: DCP dc;
492:
493: if ( NUM(p) )
494: return 1;
495: else {
496: for ( dc = DC(p), n = 0; dc; dc = NEXT(dc) )
497: n += lengthp(COEF(dc));
498: return n;
1.4 noro 499: }
500: }
501:
502: int nonzero_const_term(P p)
503: {
504: DCP dc;
505:
506: if ( !p )
507: return 0;
508: else if ( NUM(p) )
509: return 1;
510: else {
511: dc = DC(p);
512: for ( ; NEXT(dc); dc = NEXT(dc) );
513: if ( DEG(dc) )
514: return 0;
515: else
516: return nonzero_const_term(COEF(dc));
1.1 noro 517: }
518: }