![[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: }