Annotation of OpenXM_contrib2/asir2018/engine/NEZ.c, Revision 1.1
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: *
! 48: * $OpenXM$
! 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++ ) {
! 206: vn[i].v = VR(COEF(dc)); vn[i].n = QTOS(DEG(dc));
! 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 ) {
! 212: vn[j].n = MIN(vn[j].n,QTOS(DEG(dc))); break;
! 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 ) {
! 225: MKV(vn[j].v,x); STOQ(vn[j].n,d);
! 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: }
! 280: for ( g = gcd, k = QTOS(DEG(dc))-1; k > 0; k-- ) {
! 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