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