Annotation of OpenXM_contrib2/asir2018/engine/NEZ.c, Revision 1.2
1.1 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
26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
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.2 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2018/engine/NEZ.c,v 1.1 2018/09/19 05:45:07 noro Exp $
1.1 noro 49: */
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: Z cq,cq1;
61: P *pl,*pm;
62: DCP *ml;
63: Z *cl;
64: P **cp;
65: int *cn;
66:
67: if ( m == 1 ) {
68: *pr = ps[0]; return;
69: }
70: pl = (P *)ALLOCA(m*sizeof(P));
71: ml = (DCP *)ALLOCA(m*sizeof(DCP));
72: cl = (Z *)ALLOCA(m*sizeof(Z));
73: for ( i = 0; i < m; i++ )
74: monomialfctr(vl,ps[i],&pl[i],&ml[i]);
75: gcdmonomial(vl,ml,m,&mg);
76: for ( i = 0; i < m; i++ ) {
77: ptozp(pl[i],1,(Q *)&cl[i],&t); pl[i] = t;
78: }
79: for ( i = 1, cq = cl[0]; i < m; i++ ) {
80: gcdz(cl[i],cq,&cq1); cq = cq1;
81: }
82: for ( i = 0; i < m; i++ )
83: if ( NUM(pl[i]) ) {
84: mulp(vl,(P)cq,mg,pr); return;
85: }
86: for ( i = 0, nvl = vl, avl = 0; nvl && i < m; i++ ) {
87: clctv(vl,pl[i],&tvl);
88: intersectv(nvl,tvl,&svl); nvl = svl;
89: mergev(vl,avl,tvl,&svl); avl = svl;
90: }
91: if ( !nvl ) {
92: mulp(vl,(P)cq,mg,pr); return;
93: }
94: if ( !NEXT(avl) ) {
95: nuezgcdnpzmain(vl,pl,m,&t); mulp(vl,t,(P)cq,&s); mulp(vl,s,mg,pr);
96: return;
97: }
98: for ( tvl = nvl, i = 0; tvl; tvl = NEXT(tvl), i++ );
99: for ( tvl = avl, j = 0; tvl; tvl = NEXT(tvl), j++ );
100: if ( i == j ) {
101: /* all the pl[i]'s have the same variables */
102: sortplistbyhomdeg(pl,m); nezgcdnpzmain(nvl,pl,m,&t);
103: #if 0
104: /* search the minimal degree poly */
105: for ( i = 0; i < m; i++ ) {
106: for ( tvl = nvl; tvl; tvl = NEXT(tvl) ) {
107: dt = getdeg(tvl->v,pl[i]);
108: if ( tvl == nvl || dt < d1 ) {
109: v1 = tvl->v; d1 = dt;
110: }
111: }
112: if ( i == 0 || d1 < d ) {
113: v = v1; d = d1; j = i;
114: }
115: }
116: t = pl[0]; pl[0] = pl[j]; pl[j] = t;
117: if ( v != nvl->v ) {
118: reordvar(nvl,v,&mvl);
119: for ( i = 0; i < m; i++ ) {
120: reorderp(mvl,nvl,pl[i],&t); pl[i] = t;
121: }
122: nezgcdnpzmain(mvl,pl,m,&s); reorderp(nvl,mvl,s,&t);
123: } else
124: nezgcdnpzmain(nvl,pl,m,&t);
125: #endif
126: } else {
127: cp = (P **)ALLOCA(m*sizeof(P *));
128: cn = (int *)ALLOCA(m*sizeof(int));
129: for ( i = 0; i < m; i++ ) {
130: cp[i] = (P *)ALLOCA(lengthp(pl[i])*sizeof(P));
131: cn[i] = pcoef(vl,nvl,pl[i],cp[i]);
132: }
133: for ( i = j = 0; i < m; i++ )
134: j += cn[i];
135: pm = (P *)ALLOCA(j*sizeof(P));
136: for ( i = k = 0; i < m; i++ )
137: for ( j = 0; j < cn[i]; j++ )
138: pm[k++] = cp[i][j];
139: nezgcdnpz(vl,pm,k,&t);
140: }
141: mulp(vl,t,(P)cq,&s); mulp(vl,s,mg,pr);
142: }
143:
144: void sortplistbyhomdeg(p,n)
145: P *p;
146: int n;
147: {
148: int i,j,k;
149: int *l;
150: P t;
151:
152: l = (int *)ALLOCA(n*sizeof(int));
153: for ( i = 0; i < n; i++ )
154: l[i] = homdeg(p[i]);
155: for ( i = 0; i < n; i++ )
156: for ( j = i + 1; j < n; j++ )
157: if ( l[j] < l[i] ) {
158: t = p[i]; p[i] = p[j]; p[j] = t;
159: k = l[i]; l[i] = l[j]; l[j] = k;
160: }
161: }
162:
163: void nuezgcdnpzmain(vl,ps,m,r)
164: VL vl;
165: P *ps;
166: int m;
167: P *r;
168: {
169: P *tps;
170: P f,t;
171: int i;
172:
173: for ( i = 0; i < m; i++ )
174: if ( NUM(ps[i]) ) {
175: *r = (P)ONE; return;
176: }
177: tps = (P *) ALLOCA(m*sizeof(P));
178: for ( i = 0; i < m; i++ )
179: tps[i] = ps[i];
180: sortplist(tps,m);
181: for ( i = 1, f = tps[0]; i < m && !NUM(f); i++ ) {
182: uezgcdpz(vl,f,tps[i],&t); f = t;
183: }
184: *r = f;
185: }
186:
187: void gcdmonomial(vl,dcl,m,r)
188: VL vl;
189: DCP *dcl;
190: int m;
191: P *r;
192: {
193: int i,j,n;
194: P g,x,s,t;
195: Z d;
196: DCP dc;
197: VN vn;
198:
199: for ( i = 0; i < m; i++ )
200: if ( !dcl[i] ) {
201: *r = (P)ONE; return;
202: }
203: for ( n = 0, dc = dcl[0]; dc; dc = NEXT(dc), n++ );
204: vn = (VN)ALLOCA(n*sizeof(struct oVN));
205: for ( i = 0, dc = dcl[0]; i < n; dc = NEXT(dc), i++ ) {
1.2 ! noro 206: vn[i].v = VR(COEF(dc)); vn[i].n = ZTOS(DEG(dc));
1.1 noro 207: }
208: for ( i = 1; i < m; i++ ) {
209: for ( j = 0; j < n; j++ ) {
210: for ( dc = dcl[i]; dc; dc = NEXT(dc) )
211: if ( VR(COEF(dc)) == vn[j].v ) {
1.2 ! noro 212: vn[j].n = MIN(vn[j].n,ZTOS(DEG(dc))); break;
1.1 noro 213: }
214: if ( !dc )
215: vn[j].n = 0;
216: }
217: for ( j = n-1; j >= 0 && !vn[j].n; j-- );
218: if ( j < 0 ) {
219: *r = (P)ONE; return;
220: } else
221: n = j+1;
222: }
223: for ( j = 0, g = (P)ONE; j < n; j++ )
224: if ( vn[j].n ) {
1.2 ! noro 225: MKV(vn[j].v,x); STOZ(vn[j].n,d);
1.1 noro 226: pwrp(vl,x,d,&t); mulp(vl,t,g,&s); g = s;
227: }
228: *r = g;
229: }
230:
231: void nezgcdnpzmain(vl,pl,m,pr)
232: VL vl;
233: P *pl,*pr;
234: int m;
235: {
236: P *ppl,*pcl;
237: int i;
238: P cont,gcd,t,s;
239: DCP dc;
240:
241: ppl = (P *)ALLOCA(m*sizeof(P));
242: pcl = (P *)ALLOCA(m*sizeof(P));
243: for ( i = 0; i < m; i++ )
244: pcp(vl,pl[i],&ppl[i],&pcl[i]);
245: nezgcdnpz(vl,pcl,m,&cont);
246: sqfrp(vl,ppl[0],&dc);
247: for ( dc = NEXT(dc), gcd = (P)ONE; dc; dc = NEXT(dc) ) {
248: if ( NUM(COEF(dc)) )
249: continue;
250: nezgcdnpp(vl,dc,ppl+1,m-1,&t);
251: if ( NUM(t) )
252: continue;
253: mulp(vl,gcd,t,&s); gcd = s;
254: }
255: mulp(vl,gcd,cont,pr);
256: }
257:
258: void nezgcdnpp(vl,dc,pl,m,r)
259: VL vl;
260: DCP dc;
261: P *pl;
262: int m;
263: P *r;
264: {
265: int i,k;
266: P g,t,s,gcd;
267: P *pm;
268:
269: nezgcdnp_sqfr_primitive(vl,COEF(dc),pl,m,&gcd);
270: if ( NUM(gcd) ) {
271: *r = (P)ONE; return;
272: }
273: pm = (P *) ALLOCA(m*sizeof(P));
274: for ( i = 0; i < m; i++ ) {
275: divsp(vl,pl[i],gcd,&pm[i]);
276: if ( NUM(pm[i]) ) {
277: *r = gcd; return;
278: }
279: }
1.2 ! noro 280: for ( g = gcd, k = ZTOS(DEG(dc))-1; k > 0; k-- ) {
1.1 noro 281: nezgcdnp_sqfr_primitive(vl,g,pm,m,&t);
282: if ( NUM(t) )
283: break;
284: mulp(vl,gcd,t,&s); gcd = s;
285: for ( i = 0; i < m; i++ ) {
286: divsp(vl,pm[i],t,&s);
287: if ( NUM(s) ) {
288: *r = gcd; return;
289: }
290: pm[i] = s;
291: }
292: }
293: *r = gcd;
294: }
295:
296: /*
297: * pr = gcd(p0,ps[0],...,ps[m-1])
298: *
299: * p0 is already square-free and primitive.
300: * ps[i] is at least primitive.
301: *
302: */
303:
304: void nezgcdnp_sqfr_primitive(vl,p0,ps,m,pr)
305: VL vl;
306: int m;
307: P p0,*ps,*pr;
308: {
309: /* variables */
310: P p00,g,h,g0,h0,a0,b0;
311: P lp0,lgp0,lp00,lg,lg0,cbd,tq,t;
312: P *lc;
313: P *tps;
314: VL nvl1,nvl2,nvl,tvl;
315: V v;
316: int i,j,k,d0,dg,dg0,dmin,z;
317: VN vn1;
318: int nv,flag,max;
319:
320: /* set-up */
321: if ( NUM(p0) ) {
322: *pr = (P) ONE; return;
323: }
324: for ( v = VR(p0), i = 0; i < m; i++ )
325: if ( NUM(ps[i]) || (v != VR(ps[i])) ) {
326: *pr = (P)ONE; return;
327: }
328: tps = (P *) ALLOCA(m*sizeof(P));
329: for ( i = 0; i < m; i++ )
330: tps[i] = ps[i];
331: sortplist(tps,m);
332: /* deg(tps[0]) <= deg(tps[1]) <= ... */
333:
334: if ( !cmpz(DEG(DC(p0)),ONE) ) {
335: if ( divcheck(vl,tps,m,(P)ONE,p0) )
336: *pr = p0;
337: else
338: *pr = (P)ONE;
339: return;
340: }
341:
342: lp0 = LC(p0); dmin = d0 = deg(v,p0); lc = (P *)ALLOCA((m+1)*sizeof(P));
343: for ( lc[0] = lp0, i = 0; i < m; i++ )
344: lc[i+1] = LC(tps[i]);
345: clctv(vl,p0,&nvl);
346: for ( i = 0; i < m; i++ ) {
347: clctv(vl,tps[i],&nvl1); mergev(vl,nvl,nvl1,&nvl2); nvl = nvl2;
348: }
349: nezgcdnpz(nvl,lc,m+1,&lg);
350:
351: mulp(nvl,p0,lg,&lgp0); k = dbound(v,lgp0)+1; cbound(nvl,lgp0,(Q *)&cbd);
352: for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++ );
353: W_CALLOC(nv,struct oVN,vn1);
354: for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
355: vn1[i].v = tvl->v;
356:
357: /* main loop */
358: /* finally, 'max' random evaluations will be generated */
359: for ( max = 1, dg = deg(v,tps[0]) + 1; ; max = 2*max )
360: for ( z = 0; z < max; z++ ) {
361: for ( i = 0; vn1[i].v; i++ )
362: vn1[i].n = mt_genrand()%max;
363:
364: /* find b s.t. LC(p0)(b), LC(tps[i])(b) != 0 */
365:
366: substvp(nvl,p0,vn1,&p00);
367: flag = (!zerovpchk(nvl,lp0,vn1) && sqfrchk(p00));
368: for ( i = 0; flag && (i < m); i++ )
369: flag &= (!zerovpchk(nvl,LC(tps[i]),vn1));
370: if ( !flag )
371: continue;
372:
373: /* substitute y -> b */
374: substvp(nvl,lg,vn1,&lg0); lp00 = LC(p00);
375: /* extended-gcd in 1-variable */
376: uexgcdnp(nvl,p00,tps,m,vn1,(Q)cbd,&g0,&h0,&a0,&b0,(Z *)&tq);
377: if ( NUM(g0) ) {
378: *pr = (P)ONE; return;
379: } else if ( dg > ( dg0 = deg(v,g0) ) ) {
380: dg = dg0;
381: if ( dg == d0 ) {
382: if ( divcheck(nvl,tps,m,lp0,p0) ) {
383: *pr = p0; return;
384: }
385: } else if ( dg == deg(v,tps[0]) ) {
386: if ( divtpz(nvl,p0,tps[0],&t) &&
387: divcheck(nvl,tps+1,m-1,LC(tps[0]),tps[0]) ) {
388: *pr = tps[0]; return;
389: } else
390: break;
391: } else {
392: henmv(nvl,vn1,lgp0,g0,h0,a0,b0,lg,lp0,lg0,lp00,(Z)tq,k,&g,&h);
393: if ( divtpz(nvl,lgp0,g,&t) &&
394: divcheck(nvl,tps,m,lg,g) ) {
395: pcp(nvl,g,pr,&t); return;
396: }
397: }
398: }
399: }
400: }
401:
402: void intersectv(vl1,vl2,vlp)
403: VL vl1,vl2,*vlp;
404: {
405: int i,n;
406: VL tvl;
407: VN tvn;
408:
409: if ( !vl1 || !vl2 ) {
410: *vlp = 0; return;
411: }
412: for ( n = 0, tvl = vl1; tvl; tvl = NEXT(tvl), n++ );
413: tvn = (VN) ALLOCA(n*sizeof(struct oVN));
414: for ( i = 0, tvl = vl1; i < n; tvl = NEXT(tvl), i++ ) {
415: tvn[i].v = tvl->v; tvn[i].n = 0;
416: }
417: for ( tvl = vl2; tvl; tvl = NEXT(tvl) )
418: for ( i = 0; i < n; i++ )
419: if ( tvn[i].v == tvl->v ) {
420: tvn[i].n = 1; break;
421: }
422: vntovl(tvn,n,vlp);
423: }
424:
425: int pcoef(vl,ivl,p,cp)
426: VL vl,ivl;
427: P p;
428: P *cp;
429: {
430: VL nvl,tvl,svl,mvl,mvl0;
431: P t;
432:
433: if ( NUM(p) ) {
434: cp[0] = p; return 1;
435: } else {
436: clctv(vl,p,&nvl);
437: for ( tvl = nvl, mvl0 = 0; tvl; tvl = NEXT(tvl) ) {
438: for ( svl = ivl; svl; svl = NEXT(svl) )
439: if ( tvl->v == svl->v )
440: break;
441: if ( !svl ) {
442: if ( !mvl0 ) {
443: NEWVL(mvl0); mvl = mvl0;
444: } else {
445: NEWVL(NEXT(mvl)); mvl = NEXT(mvl);
446: }
447: mvl->v = tvl->v;
448: }
449: }
450: if ( mvl0 )
451: NEXT(mvl) = ivl;
452: else
453: mvl0 = ivl;
454: reorderp(mvl0,nvl,p,&t);
455: return pcoef0(mvl0,ivl,t,cp);
456: }
457: }
458:
459: int pcoef0(vl,ivl,p,cp)
460: VL vl,ivl;
461: P p;
462: P *cp;
463: {
464: int cn,n;
465: DCP dc;
466: V v;
467: VL tvl;
468:
469: if ( NUM(p) ) {
470: cp[0] = p; return 1;
471: } else {
472: for ( v = VR(p), tvl = ivl; tvl; tvl = NEXT(tvl) )
473: if ( v == tvl->v )
474: break;
475: if ( tvl ) {
476: cp[0] = p; return 1;
477: } else {
478: for ( dc = DC(p), n = 0; dc; dc = NEXT(dc) ) {
479: cn = pcoef0(vl,ivl,COEF(dc),cp); cp += cn; n += cn;
480: }
481: return n;
482: }
483: }
484: }
485:
486: int lengthp(p)
487: P p;
488: {
489: int n;
490: DCP dc;
491:
492: if ( NUM(p) )
493: return 1;
494: else {
495: for ( dc = DC(p), n = 0; dc; dc = NEXT(dc) )
496: n += lengthp(COEF(dc));
497: return n;
498: }
499: }
500:
501: int nonzero_const_term(P p)
502: {
503: DCP dc;
504:
505: if ( !p )
506: return 0;
507: else if ( NUM(p) )
508: return 1;
509: else {
510: dc = DC(p);
511: for ( ; NEXT(dc); dc = NEXT(dc) );
512: if ( DEG(dc) )
513: return 0;
514: else
515: return nonzero_const_term(COEF(dc));
516: }
517: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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 / asir2018 / engine