Annotation of OpenXM_contrib2/asir2018/engine/EZ.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 ezgcdp(VL vl,P p1,P p2,P *pr)
! 53: {
! 54: P f1,f2;
! 55: Q c;
! 56:
! 57: if ( !p1 )
! 58: if ( !p2 )
! 59: *pr = 0;
! 60: else
! 61: *pr = p2;
! 62: else if ( !p2 )
! 63: *pr = p1;
! 64: else {
! 65: if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
! 66: error("ezgcdp : invalid argument");
! 67: ptozp(p1,1,&c,&f1);
! 68: ptozp(p2,1,&c,&f2);
! 69: ezgcdpz(vl,f1,f2,pr);
! 70: }
! 71: }
! 72:
! 73: /*
! 74: * p1, p2 : integer coefficient
! 75: * returns gcd(pp(p1),pp(p2)) * gcd(cont(p1),cont(p2))
! 76: */
! 77:
! 78: void ezgcdpz(VL vl,P p1,P p2,P *pr)
! 79: {
! 80: P f1,f2,t,s,g,gcd;
! 81: P fc1,fc2,cont;
! 82: VL tvl,svl;
! 83: V v1,v2;
! 84: DCP dc,dc1,dc2;
! 85: Z c1,c2,c;
! 86: P g1,g2;
! 87: P mgcd;
! 88: extern int nez;
! 89: P pm[2];
! 90:
! 91: if ( nez ) {
! 92: pm[0] = p1; pm[1] = p2; nezgcdnpz(vl,pm,2,pr); return;
! 93: }
! 94: monomialfctr(vl,p1,&f1,&dc1); p1 = f1;
! 95: monomialfctr(vl,p2,&f2,&dc2); p2 = f2;
! 96: for ( mgcd = (P)ONE; dc1; dc1 = NEXT(dc1) )
! 97: for ( v1 = VR(dc1->c), dc = dc2; dc; dc = NEXT(dc) )
! 98: if ( v1 == VR(dc->c) ) {
! 99: pwrp(vl,dc->c,cmpz(dc1->d,dc->d)>=0?dc->d:dc1->d,&t);
! 100: mulp(vl,mgcd,t,&s); mgcd = s; break;
! 101: }
! 102:
! 103: if ( NUM(p1) ) {
! 104: if ( NUM(p2) ) {
! 105: gcdz((Z)p1,(Z)p2,&c);
! 106: } else {
! 107: ptozp(p2,1,(Q *)&c2,&t);
! 108: gcdz((Z)p1,c2,&c);
! 109: }
! 110: mulp(vl,(P)c,mgcd,pr);
! 111: return;
! 112: } else if ( NUM(p2) ) {
! 113: ptozp(p1,1,(Q *)&c1,&t);
! 114: gcdz(c1,(Z)p2,&c);
! 115: mulp(vl,(P)c,mgcd,pr);
! 116: return;
! 117: }
! 118:
! 119: ptozp(p1,1,(Q *)&c1,&g1); ptozp(p2,1,(Q *)&c2,&g2);
! 120: gcdz(c1,c2,&c);
! 121:
! 122: if ( ( v1 = VR(f1) ) != ( v2 = VR(f2) ) ) {
! 123: while ( v1 != VR(vl) && v2 != VR(vl) )
! 124: vl = NEXT(vl);
! 125: if ( v1 == VR(vl) ) {
! 126: pcp(vl,f1,&g,&fc1);
! 127: ezgcdpz(vl,fc1,f2,&t);
! 128: } else {
! 129: pcp(vl,f2,&g,&fc2);
! 130: ezgcdpz(vl,fc2,f1,&t);
! 131: }
! 132: mulp(vl,t,(P)c,&s); mulp(vl,s,mgcd,pr);
! 133: return;
! 134: }
! 135:
! 136: pcp(vl,f1,&g1,&fc1); pcp(vl,f2,&g2,&fc2);
! 137: ezgcdpz(vl,fc1,fc2,&cont);
! 138: clctv(vl,g1,&tvl); clctv(vl,g2,&svl);
! 139: if ( !NEXT(tvl) && !NEXT(svl) ) {
! 140: uezgcdpz(vl,g1,g2,&t);
! 141: mulp(vl,t,cont,&s); mulp(vl,s,(P)c,&t); mulp(vl,t,mgcd,pr);
! 142: return;
! 143: }
! 144:
! 145: if ( homdeg(g1) > homdeg(g2) ) {
! 146: t = g1; g1 = g2; g2 = t;
! 147: }
! 148: sqfrp(vl,g1,&dc);
! 149: for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) {
! 150: if ( NUM(COEF(dc)) )
! 151: continue;
! 152: ezgcdpp(vl,dc,g,&t);
! 153: if ( NUM(t) )
! 154: continue;
! 155: mulp(vl,gcd,t,&s); gcd = s;
! 156: divsp(vl,g,t,&s); g = s;
! 157: }
! 158: mulp(vl,gcd,cont,&t); mulp(vl,t,(P)c,&s); mulp(vl,s,mgcd,pr);
! 159: }
! 160:
! 161: void uezgcdpz(VL vl,P p1,P p2,P *pr)
! 162: {
! 163: P g1,g2,gcd,t,s,g;
! 164: DCP dc;
! 165: Z c1,c2,cq;
! 166: P f1,f2;
! 167:
! 168: if ( NUM(p1) ) {
! 169: if ( NUM(p2) ) {
! 170: gcdz((Z)p1,(Z)p2,&cq); *pr = (P)cq;
! 171: return;
! 172: } else {
! 173: ptozp(p2,1,(Q *)&c2,&t);
! 174: gcdz((Z)p1,c2,&cq); *pr = (P)cq;
! 175: return;
! 176: }
! 177: } else if ( NUM(p2) ) {
! 178: ptozp(p1,1,(Q *)&c1,&t);
! 179: gcdz(c1,(Z)p2,&cq); *pr = (P)cq;
! 180: return;
! 181: }
! 182:
! 183: ptozp(p1,1,(Q *)&c1,&f1); ptozp(p2,1,(Q *)&c2,&f2);
! 184: gcdz(c1,c2,&cq);
! 185: if ( UDEG(f1) > UDEG(f2) ) {
! 186: g1 = f2; g2 = f1;
! 187: } else {
! 188: g1 = f1; g2 = f2;
! 189: }
! 190: usqp(g1,&dc);
! 191: for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) {
! 192: if ( NUM(COEF(dc)) )
! 193: continue;
! 194: uezgcdpp(dc,g,&t);
! 195: if ( NUM(t) )
! 196: continue;
! 197: mulp(CO,gcd,t,&s); gcd = s;
! 198: divsp(CO,g,t,&s); g = s;
! 199: }
! 200: mulp(vl,gcd,(P)cq,pr);
! 201: }
! 202:
! 203: /*
! 204: dc : dc pair c^d
! 205: r = gcd(c^d,f)
! 206: */
! 207:
! 208: void ezgcdpp(VL vl,DCP dc,P f,P *r)
! 209: {
! 210: int k;
! 211: P g,h,t,s,gcd;
! 212:
! 213: if ( NUM(f) ) {
! 214: *r = (P)ONE;
! 215: return;
! 216: }
! 217: k = QTOS(DEG(dc)) - 1;
! 218: /* ezgcd1p(vl,COEF(dc),f,&gcd); */
! 219: ezgcdnp(vl,COEF(dc),&f,1,&gcd);
! 220: g = gcd; divsp(vl,f,gcd,&h);
! 221: for ( ; k > 0; k-- ) {
! 222: /* ezgcd1p(vl,g,h,&t); */
! 223: ezgcdnp(vl,g,&h,1,&t);
! 224: if ( NUM(t) )
! 225: break;
! 226: mulp(vl,gcd,t,&s); gcd = s;
! 227: divsp(vl,h,t,&s); h = s;
! 228: }
! 229: *r = gcd;
! 230: }
! 231:
! 232: void ezgcd1p(VL vl,P p0,P p1,P *gcdp)
! 233: {
! 234: VL nvl,tvl,vl0,vl1,vl2;
! 235: P f0,f1,t;
! 236:
! 237: clctv(vl,p0,&vl0); clctv(vl,p1,&vl1); mergev(vl,vl0,vl1,&vl2);
! 238: minchdegp(vl,p0,&tvl,&t); reordvar(vl2,tvl->v,&nvl);
! 239: reorderp(nvl,vl,p0,&f0); reorderp(nvl,vl,p1,&f1);
! 240: ezgcdnp(nvl,f0,&f1,1,&t);
! 241: reorderp(vl,nvl,t,gcdp);
! 242: }
! 243:
! 244: void uezgcdpp(DCP dc,P f,P *r)
! 245: {
! 246: int k;
! 247: P g,h,t,s,gcd;
! 248:
! 249: if ( NUM(f) ) {
! 250: *r = (P)ONE;
! 251: return;
! 252: }
! 253:
! 254: k = QTOS(DEG(dc)) - 1;
! 255: uezgcd1p(COEF(dc),f,&gcd);
! 256: g = gcd; divsp(CO,f,gcd,&h);
! 257: for ( ; k > 0; k-- ) {
! 258: uezgcd1p(g,h,&t);
! 259: if ( NUM(t) )
! 260: break;
! 261: mulp(CO,gcd,t,&s); gcd = s;
! 262: divsp(CO,h,t,&s); h = s;
! 263: }
! 264: *r = gcd;
! 265: }
! 266:
! 267: /*
! 268: * pr = gcd(p0,ps[0],...,ps[m-1])
! 269: *
! 270: * p0 is already square-free and primitive.
! 271: * ps[i] is at least primitive.
! 272: *
! 273: */
! 274:
! 275: void ezgcdnp(VL vl,P p0,P *ps,int m,P *pr)
! 276: {
! 277: /* variables */
! 278: P p00,g,h,g0,h0,a0,b0;
! 279: P lp0,lgp0,lp00,lg,lg0,t;
! 280: Q cbd;
! 281: Z tq;
! 282: P *tps;
! 283: VL nvl1,nvl2,nvl,tvl;
! 284: V v;
! 285: int i,j,k,d0,dg,dg0,dmin;
! 286: VN vn0,vn1,vnt;
! 287: int nv,flag;
! 288:
! 289: /* set-up */
! 290: if ( NUM(p0) ) {
! 291: *pr = (P) ONE;
! 292: return;
! 293: }
! 294: for ( v = VR(p0), i = 0; i < m; i++ )
! 295: if ( NUM(ps[i]) || (v != VR(ps[i])) ) {
! 296: *pr = (P)ONE;
! 297: return;
! 298: }
! 299: tps = (P *) ALLOCA(m*sizeof(P));
! 300: for ( i = 0; i < m; i++ )
! 301: tps[i] = ps[i];
! 302: sortplist(tps,m);
! 303: /* deg(tps[0]) <= deg(tps[1]) <= ... */
! 304:
! 305: if ( !cmpz(DEG(DC(p0)),ONE) ) {
! 306: if ( divcheck(vl,tps,m,(P)ONE,p0) )
! 307: *pr = p0;
! 308: else
! 309: *pr = (P)ONE;
! 310: return;
! 311: }
! 312:
! 313: lp0 = LC(p0); lg = lp0; dmin = d0 = deg(v,p0);
! 314: clctv(vl,p0,&nvl);
! 315: for ( i = 0; i < m; i++ ) {
! 316: clctv(vl,tps[i],&nvl1); mergev(vl,nvl,nvl1,&nvl2);
! 317: nvl = nvl2;
! 318: ezgcdpz(nvl,lg,LC(tps[i]),&t); lg = t;
! 319: }
! 320:
! 321: mulp(nvl,p0,lg,&lgp0);
! 322: k = dbound(v,lgp0) + 1;
! 323: cbound(nvl,lgp0,(Q *)&cbd);
! 324: for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++ );
! 325: W_CALLOC(nv,struct oVN,vn0); W_CALLOC(nv,struct oVN,vnt);
! 326: W_CALLOC(nv,struct oVN,vn1);
! 327: for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
! 328: vn1[i].v = vn0[i].v = tvl->v;
! 329:
! 330: /* main loop */
! 331: for ( dg = deg(v,tps[0]) + 1; ; next(vn0) )
! 332: do {
! 333: for ( i = 0, j = 0; vn0[i].v; i++ )
! 334: if ( vn0[i].n ) vnt[j++].v = (V)((unsigned long)i);
! 335: vnt[j].n = 0;
! 336:
! 337: /* find b s.t. LC(p0)(b), LC(tps[i])(b) != 0 */
! 338: mulsgn(vn0,vnt,j,vn1);
! 339: substvp(nvl,p0,vn1,&p00);
! 340: flag = (!zerovpchk(nvl,lp0,vn1) && sqfrchk(p00));
! 341: for ( i = 0; flag && (i < m); i++ )
! 342: flag &= (!zerovpchk(nvl,LC(tps[i]),vn1));
! 343: if ( !flag )
! 344: continue;
! 345:
! 346: /* substitute y -> b */
! 347: substvp(nvl,lg,vn1,&lg0);
! 348: lp00 = LC(p00);
! 349: /* extended-gcd in 1-variable */
! 350: uexgcdnp(nvl,p00,tps,m,vn1,(Q)cbd,&g0,&h0,&a0,&b0,&tq);
! 351: if ( NUM(g0) ) {
! 352: *pr = (P)ONE;
! 353: return;
! 354: } else if ( dg > ( dg0 = deg(v,g0) ) ) {
! 355: dg = dg0;
! 356: if ( dg == d0 ) {
! 357: if ( divcheck(nvl,tps,m,lp0,p0) ) {
! 358: *pr = p0;
! 359: return;
! 360: }
! 361: } else if ( dg == deg(v,tps[0]) ) {
! 362: if ( divtpz(nvl,p0,tps[0],&t) &&
! 363: divcheck(nvl,tps+1,m-1,LC(tps[0]),tps[0]) ) {
! 364: *pr = tps[0];
! 365: return;
! 366: } else
! 367: break;
! 368: } else {
! 369: henmv(nvl,vn1,lgp0,g0,h0,a0,b0,
! 370: lg,lp0,lg0,lp00,tq,k,&g,&h);
! 371: if ( divtpz(nvl,lgp0,g,&t) &&
! 372: divcheck(nvl,tps,m,lg,g) ) {
! 373: pcp(nvl,g,pr,&t);
! 374: return;
! 375: }
! 376: }
! 377: }
! 378: } while ( !nextbin(vnt,j) );
! 379: }
! 380:
! 381: /*
! 382: f : sqfr
! 383: */
! 384:
! 385: void uezgcd1p(P f,P g,P *r)
! 386: {
! 387: UM wf,wf1,wf2,wg,wg1,wgcd,wcof;
! 388: int d1,d2,dg,m,index,n;
! 389: ML list;
! 390: DCP dc;
! 391: P t;
! 392: Q c;
! 393:
! 394: if ( NUM(f) || NUM(g) ) {
! 395: *r = (P)ONE;
! 396: return;
! 397: }
! 398: ptozp(f,1,&c,&t); f = t; ptozp(g,1,&c,&t); g = t;
! 399: d1 = UDEG(f); d2 = UDEG(g);
! 400: n = MAX(d1,d2)+1;
! 401: wf = W_UMALLOC(n); wf1 = W_UMALLOC(n);
! 402: wf2 = W_UMALLOC(n); wg = W_UMALLOC(n);
! 403: wg1 = W_UMALLOC(n); wgcd = W_UMALLOC(n);
! 404: wcof = W_UMALLOC(n);
! 405:
! 406: for ( index = INDEX, dg = MIN(d1,d2)+1; ; ) {
! 407: m = sprime[index++];
! 408: if ( !remqi((Q)UCOEF(f),m) || !remqi((Q)UCOEF(g),m))
! 409: continue;
! 410: ptoum(m,f,wf); cpyum(wf,wf1);
! 411: diffum(m,wf1,wf2); gcdum(m,wf1,wf2,wgcd);
! 412: if ( DEG(wgcd) > 0 )
! 413: continue;
! 414: ptoum(m,g,wg); cpyum(wg,wg1); cpyum(wf,wf1);
! 415: gcdum(m,wf1,wg1,wgcd);
! 416: if ( dg <= DEG(wgcd) )
! 417: continue;
! 418: else
! 419: dg = DEG(wgcd);
! 420: if ( dg == 0 ) {
! 421: *r = (P)ONE;
! 422: return;
! 423: } else if ( dg == d1 ) {
! 424: if ( divtpz(CO,g,f,&t) ) {
! 425: *r = f;
! 426: return;
! 427: }
! 428: } else if ( dg == d2 ) {
! 429: if ( divtpz(CO,f,g,&t) ) {
! 430: *r = g;
! 431: return;
! 432: }
! 433: } else {
! 434: divum(m,wf,wgcd,wcof);
! 435: ezgcdhensel(f,m,wcof,wgcd,&list);
! 436: dtest(f,list,1,&dc);
! 437: if ( NEXT(dc) ) {
! 438: if ( divtpz(CO,g,COEF(dc),&t) ) {
! 439: *r = COEF(dc);
! 440: return;
! 441: }
! 442: }
! 443: }
! 444: }
! 445: }
! 446:
! 447: void uexgcdnp(VL vl,P p1,P *ps,int n,VN vn,Q cbd,P *g,P *h,P *a,P *b,Z *q)
! 448: {
! 449: UM wg,wh,wg1,wh1,wgcd,wt;
! 450: P t,s,gcd,cof,gk,hk,ak,bk;
! 451: Q c,bb;
! 452: Z tq,tq1,qk;
! 453: int mod,k,i,index,d;
! 454:
! 455: ptozp(p1,1,&c,&gcd);
! 456: for ( i = 0; i < n; i++ ) {
! 457: substvp(vl,ps[i],vn,&t);
! 458: uezgcd1p(gcd,t,&s);
! 459: if ( NUM(s) ) {
! 460: *g = (P)ONE; *h = p1; *a = 0; *b = (P)ONE;
! 461: return;
! 462: } else
! 463: gcd = s;
! 464: }
! 465:
! 466: if ( !dcomp(p1,gcd) ) {
! 467: *g = p1; *h = (P)ONE; *a = 0; *b = (P)ONE;
! 468: return;
! 469: }
! 470:
! 471: divsp(CO,p1,gcd,&cof);
! 472: d = UDEG(p1)+1;
! 473: wg = W_UMALLOC(d); wh = W_UMALLOC(d);
! 474: wg1 = W_UMALLOC(d); wh1 = W_UMALLOC(d);
! 475: wgcd = W_UMALLOC(d); wt = W_UMALLOC(d);
! 476: for ( index = INDEX; ; ) {
! 477: mod = sprime[index++];
! 478: if ( !remqi((Q)LC(p1),mod) )
! 479: continue;
! 480: ptoum(mod,gcd,wg); ptoum(mod,cof,wh);
! 481: cpyum(wg,wg1); cpyum(wh,wh1);
! 482: gcdum(mod,wg,wh,wt); cpyum(wt,wgcd);
! 483: if ( DEG(wgcd) > 0 )
! 484: continue;
! 485: STOQ(mod,tq); addq(cbd,cbd,&bb);
! 486: for ( k = 1; cmpq((Q)tq,bb) < 0; k++ ) {
! 487: mulz(tq,tq,&tq1); tq = tq1;
! 488: }
! 489: henzq(p1,gcd,wg1,cof,wh1,mod,k,&gk,&hk,&bk,&ak,&qk);
! 490: break;
! 491: }
! 492: *g = gk; *h = hk; *a = ak; *b = bk; *q = tq;
! 493: }
! 494:
! 495: void ezgcdhensel(P f,int mod,UM cof,UM gcd,ML *listp)
! 496: {
! 497: register int i,j;
! 498: int q,n,bound;
! 499: int *p;
! 500: int **pp;
! 501: ML blist,clist,bqlist,cqlist,rlist;
! 502: UM *b;
! 503: LUM fl,tl;
! 504: LUM *l;
! 505: LUM LUMALLOC();
! 506:
! 507: blist = W_MLALLOC(2);
! 508: blist->n = 2; blist->mod = mod;
! 509: blist->c[0] = (pointer)cof; blist->c[1] = (pointer)gcd;
! 510: gcdgen(f,blist,&clist);
! 511: henprep(f,blist,clist,&bqlist,&cqlist);
! 512: n = bqlist->n; q = bqlist->mod;
! 513: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 514:
! 515: if ( bound == 1 ) {
! 516: *listp = rlist = MLALLOC(n);
! 517: rlist->n = n;
! 518: rlist->mod = q;
! 519: rlist->bound = bound;
! 520:
! 521: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c;
! 522: i < n; i++ ) {
! 523: tl = LUMALLOC(DEG(b[i]),1);
! 524: l[i] = tl;
! 525: p = COEF(b[i]);
! 526:
! 527: for ( j = 0, pp = COEF(tl);
! 528: j <= DEG(tl); j++ )
! 529: pp[j][0] = p[j];
! 530: }
! 531: } else {
! 532: W_LUMALLOC((int)UDEG(f),bound,fl);
! 533: ptolum(q,bound,f,fl);
! 534: henmain(fl,bqlist,cqlist,listp);
! 535: }
! 536: }
! 537:
! 538: void ezgcdnpz(VL vl,P *ps,int m,P *pr)
! 539: {
! 540: P t,s,gcd;
! 541: P cont;
! 542: VL tvl,svl,nvl;
! 543: int i;
! 544: DCP dc;
! 545: Z cq;
! 546: P *pl,*ppl,*pcl;
! 547: Z *cl;
! 548: V v;
! 549:
! 550: pl = (P *)ALLOCA((m+1)*sizeof(P));
! 551: cl = (Z *)ALLOCA((m+1)*sizeof(Q));
! 552: for ( i = 0; i < m; i++ )
! 553: ptozp(ps[i],1,(Q *)&cl[i],&pl[i]);
! 554: for ( i = 1, cq = cl[0]; i < m; i++ ) {
! 555: gcdz(cl[i],cq,&cq);
! 556: }
! 557: for ( i = 0; i < m; i++ )
! 558: if ( NUM(pl[i]) ) {
! 559: *pr = (P)cq;
! 560: return;
! 561: }
! 562:
! 563: for ( i = 0, nvl = 0; i < m; i++ ) {
! 564: clctv(vl,ps[i],&tvl); mergev(vl,nvl,tvl,&svl); nvl = svl;
! 565: }
! 566:
! 567: ppl = (P *)ALLOCA((m+1)*sizeof(P));
! 568: pcl = (P *)ALLOCA((m+1)*sizeof(P));
! 569: for ( i = 0, v = nvl->v; i < m; i++ )
! 570: if ( v == VR(pl[i]) )
! 571: pcp(nvl,pl[i],&ppl[i],&pcl[i]);
! 572: else {
! 573: ppl[i] = (P)ONE; pcl[i] = pl[i];
! 574: }
! 575: ezgcdnpz(nvl,pcl,m,&cont);
! 576:
! 577: for ( i = 0; i < m; i++ )
! 578: if ( NUM(ppl[i]) ) {
! 579: mulp(nvl,cont,(P)cq,pr);
! 580: return;
! 581: }
! 582: sortplist(ppl,m);
! 583: sqfrp(nvl,ppl[0],&dc);
! 584: for ( dc = NEXT(dc), gcd = (P)ONE; dc; dc = NEXT(dc) ) {
! 585: if ( NUM(COEF(dc)) )
! 586: continue;
! 587: ezgcdnpp(vl,dc,ppl+1,m-1,&t);
! 588: if ( NUM(t) )
! 589: continue;
! 590: mulp(vl,gcd,t,&s); gcd = s;
! 591: }
! 592: mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,pr);
! 593: }
! 594:
! 595: void ezgcdnpp(VL vl,DCP dc,P *pl,int m,P *r)
! 596: {
! 597: int i,k;
! 598: P g,t,s,gcd;
! 599: P *pm;
! 600:
! 601: ezgcdnp(vl,COEF(dc),pl,m,&gcd);
! 602: if ( NUM(gcd) ) {
! 603: *r = (P)ONE;
! 604: return;
! 605: }
! 606: pm = (P *) ALLOCA((m+1)*sizeof(P));
! 607: for ( i = 0; i < m; i++ ) {
! 608: divsp(vl,pl[i],gcd,&pm[i]);
! 609: if ( NUM(pm[i]) ) {
! 610: *r = gcd;
! 611: return;
! 612: }
! 613: }
! 614: for ( g = gcd, k = QTOS(DEG(dc)) - 1; k > 0; k-- ) {
! 615: ezgcdnp(vl,g,pm,m,&t);
! 616: if ( NUM(t) )
! 617: break;
! 618: mulp(vl,gcd,t,&s); gcd = s;
! 619: for ( i = 0; i < m; i++ ) {
! 620: divsp(vl,pm[i],t,&s);
! 621: if ( NUM(s) ) {
! 622: *r = gcd;
! 623: return;
! 624: }
! 625: pm[i] = s;
! 626: }
! 627: }
! 628: *r = gcd;
! 629: }
! 630:
! 631: void pcp(VL vl,P p,P *pp,P *c)
! 632: {
! 633: P f,g,n;
! 634: DCP dc;
! 635: P *pl;
! 636: int i,m;
! 637: extern int nez;
! 638:
! 639: if ( NUM(p) ) {
! 640: *c = p;
! 641: *pp = (P)ONE;
! 642: } else {
! 643: for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
! 644: if ( m == 1 ) {
! 645: *c = COEF(DC(p));
! 646: divsp(vl,p,*c,pp);
! 647: return;
! 648: }
! 649: ptozp(p,1,(Q *)&n,&f);
! 650: pl = (P *)ALLOCA((m+1)*sizeof(P));
! 651: for ( i = 0, dc = DC(f); i < m; dc = NEXT(dc), i++ )
! 652: pl[i] = COEF(dc);
! 653: if ( nez )
! 654: nezgcdnpz(vl,pl,m,&g);
! 655: else
! 656: ezgcdnpz(vl,pl,m,&g);
! 657: mulp(vl,g,n,c); divsp(vl,f,g,pp);
! 658: }
! 659: }
! 660:
! 661: int divcheck(VL vl,P *ps,int m,P l,P p)
! 662: {
! 663: int i;
! 664: P r,t;
! 665:
! 666: for ( i = 0; i < m; i++ ) {
! 667: mulp(vl,ps[i],l,&t);
! 668: if ( !divtpz(vl,t,p,&r) )
! 669: return ( 0 );
! 670: }
! 671: return ( 1 );
! 672: }
! 673:
! 674: void sortplist(P *plist,int n)
! 675: {
! 676: int i,j;
! 677: P t;
! 678:
! 679: for ( i = 0; i < n; i++ )
! 680: for ( j = i + 1; j < n; j++ )
! 681: if ( cmpz(DEG(DC(plist[i])),DEG(DC(plist[j]))) > 0 ) {
! 682: t = plist[i]; plist[i] = plist[j]; plist[j] = t;
! 683: }
! 684: }
! 685:
! 686: int dcomp(P p1,P p2)
! 687: {
! 688: return ( cmpz(DEG(DC(p1)),DEG(DC(p2))) );
! 689: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to EZ.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / engine