![[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|
Return to EZ.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
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.3 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/EZ.c,v 1.2 2000/08/21 08:31:24 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: {
57: P f1,f2;
58: Q c;
59:
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: }
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: {
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);
166: }
167:
168: void uezgcdpz(vl,p1,p2,pr)
169: VL vl;
170: P p1,p2,*pr;
171: {
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);
211: }
212:
213: /*
214: dc : dc pair c^d
215: r = gcd(c^d,f)
216: */
217:
218: void ezgcdpp(vl,dc,f,r)
219: VL vl;
220: DCP dc;
221: P f;
222: P *r;
223: {
224: int k;
225: P g,h,t,s,gcd;
226:
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;
244: }
245:
246: void ezgcd1p(vl,p0,p1,gcdp)
247: VL vl;
248: P p0,p1;
249: P *gcdp;
250: {
251: VL nvl,tvl,vl0,vl1,vl2;
252: P f0,f1,t;
253:
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);
259: }
260:
261: void uezgcdpp(dc,f,r)
262: DCP dc;
263: P f;
264: P *r;
265: {
266: int k;
267: P g,h,t,s,gcd;
268:
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;
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: {
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) );
400: }
401:
402: /*
403: f : sqfr
404: */
405:
406: void uezgcd1p(f,g,r)
407: P f,g;
408: P *r;
409: {
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: }
468: }
469:
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: {
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;
523: }
524:
525: void ezgcdhensel(f,mod,cof,gcd,listp)
526: P f;
527: int mod;
528: UM gcd,cof;
529: ML *listp;
530: {
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: }
570: }
571:
572: void ezgcdnpz(vl,ps,m,pr)
573: VL vl;
574: P *ps,*pr;
575: int m;
576: {
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);
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: {
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;
672: }
673:
674: void pcp(vl,p,pp,c)
675: VL vl;
676: P p;
677: P *pp,*c;
678: {
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: }
705: }
706:
707: #if 0
708: pcp(vl,p,pp,c)
709: VL vl;
710: P p;
711: P *pp,*c;
712: {
713: P f,g,n,t;
714: DCP dc;
715:
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: }
727: }
728: #endif
729:
730: int divcheck(vl,ps,m,l,p)
731: VL vl;
732: int m;
733: P *ps,l,p;
734: {
735: int i;
736: P r,t;
737:
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 );
744: }
745:
746: void sortplist(plist,n)
747: P *plist;
748: int n;
749: {
750: int i,j;
751: P t;
752:
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: }
758: }
759:
760: int dcomp(p1,p2)
761: P p1,p2;
762: {
763: return ( cmpq(DEG(DC(p1)),DEG(DC(p2))) );
764: }