Annotation of OpenXM_contrib2/asir2018/engine/EZ.c, Revision 1.2
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: *
1.2 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2018/engine/EZ.c,v 1.1 2018/09/19 05:45:06 noro Exp $
1.1 noro 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: }
1.2 ! noro 217: k = ZTOS(DEG(dc)) - 1;
1.1 noro 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:
1.2 ! noro 254: k = ZTOS(DEG(dc)) - 1;
1.1 noro 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;
1.2 ! noro 485: STOZ(mod,tq); addq(cbd,cbd,&bb);
1.1 noro 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: }
1.2 ! noro 614: for ( g = gcd, k = ZTOS(DEG(dc)) - 1; k > 0; k-- ) {
1.1 noro 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