Annotation of OpenXM_contrib2/asir2000/engine/NEZ.c, Revision 1.7
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.7 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/NEZ.c,v 1.6 2010/02/19 08:27:19 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: {
1.7 ! noro 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:
! 68: if ( m == 1 ) {
! 69: *pr = ps[0]; return;
! 70: }
! 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);
1.1 noro 104: #if 0
1.7 ! noro 105: /* search the minimal degree poly */
! 106: for ( i = 0; i < m; i++ ) {
! 107: for ( tvl = nvl; tvl; tvl = NEXT(tvl) ) {
1.1 noro 108: dt = getdeg(tvl->v,pl[i]);
109: if ( tvl == nvl || dt < d1 ) {
110: v1 = tvl->v; d1 = dt;
111: }
1.7 ! noro 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);
1.1 noro 126: #endif
1.7 ! noro 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);
1.1 noro 143: }
144:
145: void sortplistbyhomdeg(p,n)
146: P *p;
147: int n;
148: {
1.7 ! noro 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: }
1.1 noro 162: }
163:
164: void nuezgcdnpzmain(vl,ps,m,r)
165: VL vl;
166: P *ps;
167: int m;
168: P *r;
169: {
1.7 ! noro 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;
1.1 noro 186: }
187:
188: void gcdmonomial(vl,dcl,m,r)
189: VL vl;
190: DCP *dcl;
191: int m;
192: P *r;
193: {
1.7 ! noro 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;
1.1 noro 230: }
231:
232: void nezgcdnpzmain(vl,pl,m,pr)
233: VL vl;
234: P *pl,*pr;
235: int m;
236: {
1.7 ! noro 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);
1.1 noro 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: {
1.7 ! noro 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;
1.1 noro 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: {
1.7 ! noro 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;
! 317: int i,j,k,d0,dg,dg0,dmin,z;
! 318: VN vn1;
! 319: int nv,flag,max;
! 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++ )
! 356: vn1[i].v = tvl->v;
! 357:
! 358: /* main loop */
! 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;
! 364:
! 365: /* find b s.t. LC(p0)(b), LC(tps[i])(b) != 0 */
! 366:
! 367: substvp(nvl,p0,vn1,&p00);
! 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: }
! 400: }
1.1 noro 401: }
402:
403: void intersectv(vl1,vl2,vlp)
404: VL vl1,vl2,*vlp;
405: {
1.7 ! noro 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);
1.1 noro 424: }
425:
426: int pcoef(vl,ivl,p,cp)
427: VL vl,ivl;
428: P p;
429: P *cp;
430: {
1.7 ! noro 431: VL nvl,tvl,svl,mvl,mvl0;
! 432: P t;
1.1 noro 433:
1.7 ! noro 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: }
1.1 noro 458: }
459:
460: int pcoef0(vl,ivl,p,cp)
461: VL vl,ivl;
462: P p;
463: P *cp;
464: {
1.7 ! noro 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: }
1.1 noro 485: }
486:
487: int lengthp(p)
488: P p;
489: {
1.7 ! noro 490: int n;
! 491: DCP dc;
1.1 noro 492:
1.7 ! noro 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;
! 499: }
1.4 noro 500: }
501:
502: int nonzero_const_term(P p)
503: {
1.7 ! noro 504: DCP dc;
1.4 noro 505:
1.7 ! noro 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));
! 517: }
1.1 noro 518: }
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 / asir2000 / engine