![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to PUM.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 PUM.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/PUM.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3:
4: void gcdcmp(), sprsm();
5:
6: void detmp(vl,mod,rmat,n,dp)
7: VL vl;
8: int mod;
9: P **rmat;
10: int n;
11: P *dp;
12: {
13: int i,j,k,sgn;
14: P mjj,mij,t,s,u,d;
15: P **mat;
16: P *mi,*mj;
17:
18: mat = (P **)almat_pointer(n,n);
19: for ( i = 0; i < n; i++ )
20: for ( j = 0; j < n; j++ )
21: mat[i][j] = rmat[i][j];
22: for ( j = 0, d = (P)ONEM, sgn = 1; j < n; j++ ) {
23: for ( i = j; (i < n) && !mat[i][j]; i++ );
24: if ( i == n ) {
25: *dp = 0; return;
26: }
27: for ( k = i; k < n; k++ )
28: if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
29: i = k;
30: if ( j != i ) {
31: mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
32: }
33: for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
34: for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
35: mulmp(vl,mod,mi[k],mjj,&t); mulmp(vl,mod,mj[k],mij,&s);
36: submp(vl,mod,t,s,&u); divsmp(vl,mod,u,d,&mi[k]);
37: }
38: d = mjj;
39: }
40: if ( sgn < 0 )
41: chsgnmp(mod,d,dp);
42: else
43: *dp = d;
44: }
45:
46: void resultmp(vl,mod,v,p1,p2,pr)
47: VL vl;
48: int mod;
49: V v;
50: P p1,p2,*pr;
51: {
52: P q1,q2,m,m1,m2,lc,q,r,t,g1,g2,adj;
53: int d,d1,d2,j,k;
54: VL nvl;
55: Q dq;
56: MQ mq;
57:
58: if ( !p1 || !p2 ) {
59: *pr = 0; return;
60: }
61: reordvar(vl,v,&nvl);
62: reordermp(nvl,mod,vl,p1,&q1); reordermp(nvl,mod,vl,p2,&q2);
63:
64: if ( VR(q1) != v )
65: if ( VR(q2) != v ) {
66: *pr = 0;
67: return;
68: } else {
69: d = deg(v,q2); STOQ(d,dq);
70: pwrmp(vl,mod,q1,dq,pr);
71: return;
72: }
73: else if ( VR(q2) != v ) {
74: d = deg(v,q1); STOQ(d,dq);
75: pwrmp(vl,mod,q2,dq,pr);
76: return;
77: }
78:
79: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
80: *pr = 0;
81: return;
82: }
83:
84: d1 = deg(v,q1); d2 = deg(v,q2);
85: if ( d1 > d2 ) {
86: g1 = q1; g2 = q2; adj = (P)ONEM;
87: } else if ( d1 < d2 ) {
88: g2 = q1; g1 = q2;
89: if ( (d1 % 2) && (d2 % 2) ) {
90: STOMQ(-1,mq); adj = (P)mq;
91: } else
92: adj = (P)ONEM;
93: } else {
94: premmp(nvl,mod,q1,q2,&t);
95: d = deg(v,t); STOQ(d,dq); pwrmp(nvl,mod,LC(q2),dq,&adj);
96: g1 = q2; g2 = t;
97: if ( d1 % 2 ) {
98: chsgnmp(mod,adj,&t); adj = t;
99: }
100: }
101: d1 = deg(v,g1); j = d1 - 1;
102:
103: for ( lc = (P)ONEM; ; ) {
104: if ( ( k = deg(v,g2) ) < 0 ) {
105: *pr = 0;
106: return;
107: }
108:
109: if ( k == j )
110: if ( !k ) {
111: divsmp(nvl,mod,g2,adj,pr);
112: return;
113: } else {
114: premmp(nvl,mod,g1,g2,&r); mulmp(nvl,mod,lc,lc,&m);
115: divsmp(nvl,mod,r,m,&q);
116: g1 = g2; g2 = q;
117: lc = LC(g1); /* g1 is not const */
118: j = k - 1;
119: }
120: else {
121: d = j - k; STOQ(d,dq);
122: pwrmp(nvl,mod,(VR(g2)==v?LC(g2):g2),dq,&m);
123: mulmp(nvl,mod,g2,m,&m1);
124: pwrmp(nvl,mod,lc,dq,&m); divsmp(nvl,mod,m1,m,&t);
125: if ( k == 0 ) {
126: divsmp(nvl,mod,t,adj,pr);
127: return;
128: } else {
129: premmp(nvl,mod,g1,g2,&r);
130: mulmp(nvl,mod,lc,lc,&m1); mulmp(nvl,mod,m,m1,&m2);
131: divsmp(nvl,mod,r,m2,&q);
132: g1 = t; g2 = q;
133: if ( d % 2 ) {
134: chsgnmp(mod,g2,&t); g2 = t;
135: }
136: lc = LC(g1); /* g1 is not const */
137: j = k - 1;
138: }
139: }
140: }
141: }
142:
143: void premmp(vl,mod,p1,p2,pr)
144: VL vl;
145: int mod;
146: P p1,p2,*pr;
147: {
148: P m,m1,m2;
149: P *pw;
150: DCP dc;
151: V v1,v2;
152: register int i,j;
153: int n1,n2,d;
154:
155: if ( NUM(p1) )
156: if ( NUM(p2) )
157: *pr = 0;
158: else
159: *pr = p1;
160: else if ( NUM(p2) )
161: *pr = 0;
162: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
163: n1 = deg(v1,p1); n2 = deg(v1,p2);
164: pw = (P *)ALLOCA((n1+1)*sizeof(P));
165: bzero((char *)pw,(int)((n1+1)*sizeof(P)));
166:
167: for ( dc = DC(p1); dc; dc = NEXT(dc) )
168: pw[QTOS(DEG(dc))] = COEF(dc);
169:
170: for ( i = n1; i >= n2; i-- ) {
171: if ( pw[i] ) {
172: m = pw[i];
173: for ( j = i; j >= 0; j-- ) {
174: mulmp(vl,mod,pw[j],LC(p2),&m1); pw[j] = m1;
175: }
176:
177: for ( dc = DC(p2), d = i - n2; dc; dc = NEXT(dc) ) {
178: mulmp(vl,mod,COEF(dc),m,&m1);
179: submp(vl,mod,pw[QTOS(DEG(dc))+d],m1,&m2);
180: pw[QTOS(DEG(dc))+d] = m2;
181: }
182: } else
183: for ( j = i; j >= 0; j-- )
184: if ( pw[j] ) {
185: mulmp(vl,mod,pw[j],LC(p2),&m1); pw[j] = m1;
186: }
187: }
188: plisttop(pw,v1,n2-1,pr);
189: } else {
190: while ( v1 != vl->v && v2 != vl->v )
191: vl = NEXT(vl);
192: if ( v1 == vl->v )
193: *pr = 0;
194: else
195: *pr = p1;
196: }
197: }
198:
199: void srchmp(vl,mod,v,p1,p2,pr)
200: VL vl;
201: int mod;
202: V v;
203: P p1,p2,*pr;
204: {
205: P a,b,q1,q2,x,t,s,d,bg,c,c0,db;
206: int i,m,k;
207: V v0;
208: VL nvl,tvl,nvl1,nvl2;
209: Q dq;
210: MQ q;
211:
212: if ( vl->v != v ) {
213: reordvar(vl,v,&tvl);
214: reordermp(tvl,mod,vl,p1,&q1); reordermp(tvl,mod,vl,p2,&q2);
215: } else {
216: q1 = p1; q2 = p2; tvl = vl;
217: }
218: clctv(tvl,q1,&nvl1); clctv(tvl,q2,&nvl2); mergev(tvl,nvl1,nvl2,&nvl);
219: if ( VR(q1) != v )
220: if ( VR(q2) != v )
221: *pr = 0;
222: else {
223: m = getdeg(v,q2); STOQ(m,dq); pwrmp(vl,mod,q1,dq,pr);
224: }
225: else if ( VR(q2) != v ) {
226: m = getdeg(v,q1); STOQ(m,dq); pwrmp(vl,mod,q2,dq,pr);
227: } else if ( !NEXT(nvl) )
228: srchump(mod,p1,p2,pr);
229: else {
230: v0 = NEXT(nvl)->v;
231: k = getdeg(v,q1)*getdeg(v0,q2)+getdeg(v,q2)*getdeg(v0,q1)+1;
232: for ( i = 0, c = 0, d = (P)ONEM, MKMV(v0,x);
233: ( i < mod ) && (getdeg(v0,d) < k) ; i++ ) {
234: STOMQ(i,q),bg = (P)q; substmp(nvl,mod,LC(q1),v0,bg,&t);
235: if ( !t )
236: continue;
237: substmp(nvl,mod,LC(q2),v0,bg,&t);
238: if ( !t )
239: continue;
240: substmp(nvl,mod,q1,v0,bg,&a); substmp(nvl,mod,q2,v0,bg,&b);
241: srchmp(nvl,mod,v,a,b,&c0); substmp(nvl,mod,c,v0,bg,&t);
242: submp(nvl,mod,c0,t,&s); mulmp(nvl,mod,s,d,&t);
243: substmp(nvl,mod,d,v0,bg,&db);
244: divsmp(nvl,mod,t,db,&s); addmp(nvl,mod,s,c,&t); c = t;
245: submp(nvl,mod,x,bg,&t); mulmp(nvl,mod,d,t,&s); d = s;
246: }
247: if ( i == mod )
248: error("srchmp : ???");
249: *pr = c;
250: }
251: }
252:
253: int ucmpp(p,q)
254: P p,q;
255: {
256: DCP dcp,dcq;
257:
258: if ( !p )
259: if ( !q )
260: return ( 0 );
261: else
262: return ( 1 );
263: else if ( !q )
264: return ( 1 );
265: else if ( NUM(p) )
266: if ( !NUM(q) )
267: return ( 1 );
268: else
269: return ( cmpq((Q)p,(Q)q) );
270: else if ( NUM(q) )
271: return ( 1 );
272: else {
273: for ( dcp = DC(p), dcq = DC(q); dcp && dcq;
274: dcp = NEXT(dcp), dcq = NEXT(dcq) )
275: if ( cmpq(DEG(dcp),DEG(dcq) ) )
276: return ( 1 );
277: else if ( cmpq((Q)COEF(dcp),(Q)COEF(dcq) ) )
278: return ( 1 );
279: if ( dcp || dcq )
280: return ( 1 );
281: else
282: return ( 0 );
283: }
284: }
285:
286: #if 0
287: srchump(mod,p1,p2,pr)
288: int mod;
289: P p1,p2,*pr;
290: {
291: int r;
292: MQ mq;
293:
294: r = eucresum(mod,p1,p2);
295: STOMQ(r,mq); *pr = (P)mq;
296: }
297: #endif
298:
299: void srchump(mod,p1,p2,pr)
300: int mod;
301: P p1,p2,*pr;
302: {
303: UM m,m1,q,r,t,g1,g2;
304: int lc,d,d1,d2,i,j,k,l,l1,l2,tmp,adj;
305: V v;
306:
307: v = VR(p1); d = MAX(UDEG(p1),UDEG(p2));
308: g1 = W_UMALLOC(d); g2 = W_UMALLOC(d);
309: bzero((char *)g1,(int)((d+2)*sizeof(int))); bzero((char *)g2,(int)((d+2)*sizeof(int)));
310: if ( d == (int)UDEG(p1) ) {
311: mptoum(p1,g1); mptoum(p2,g2);
312: } else {
313: mptoum(p2,g1); mptoum(p1,g2);
314: }
315: if ( ( d1 = DEG(g1) ) > ( d2 = DEG(g2) ) ) {
316: j = d1 - 1; adj = 1;
317: } else
318: j = d2;
319: lc = 1;
320: r = W_UMALLOC(d1+d2); q = W_UMALLOC(d1+d2);
321: m1 = W_UMALLOC(d1+d2); t = W_UMALLOC(d1+d2);
322: bzero((char *)r,(int)((d1+d2+2)*sizeof(int))); bzero((char *)q,(int)((d1+d2+2)*sizeof(int)));
323: bzero((char *)m1,(int)((d1+d2+2)*sizeof(int))); bzero((char *)t,(int)((d1+d2+2)*sizeof(int)));
324: m = W_UMALLOC(0); bzero((char *)m,(int)(2*sizeof(int)));
325: adj = pwrm(mod,COEF(g2)[DEG(g2)],DEG(g1));
326: DEG(m) = 0; COEF(m)[0] = invm(COEF(g2)[DEG(g2)],mod);
327: mulum(mod,g2,m,r); cpyum(r,g2);
328: while ( 1 ) {
329: if ( ( k = DEG(g2) ) < 0 ) {
330: *pr = 0;
331: return;
332: }
333: if ( k == j ) {
334: if ( k == 0 ) {
335: DEG(m) = 0; COEF(m)[0] = adj;
336: mulum(mod,g2,m,r); umtomp(v,r,pr);
337: return;
338: } else {
339: DEG(m) = 0;
340: COEF(m)[0] = pwrm(mod,COEF(g2)[k],DEG(g1)-k+1);
341: mulum(mod,g1,m,r); DEG(r) = divum(mod,r,g2,t);
342: DEG(m) = 0; COEF(m)[0] = dmb(mod,lc,lc,&tmp);
343: divum(mod,r,m,q); cpyum(g2,g1); cpyum(q,g2);
344: lc = COEF(g1)[DEG(g1)]; j = k - 1;
345: }
346: } else {
347: d = j - k;
348: DEG(m) = 0; COEF(m)[0] = pwrm(mod,COEF(g2)[DEG(g2)],d);
349: mulum(mod,g2,m,m1); l = pwrm(mod,lc,d);
350: DEG(m) = 0; COEF(m)[0] = l; divum(mod,m1,m,t);
351: if ( k == 0 ) {
352: DEG(m) = 0; COEF(m)[0] = adj;
353: mulum(mod,t,m,r); umtomp(v,r,pr);
354: return;
355: } else {
356: DEG(m) = 0;
357: COEF(m)[0] = pwrm(mod,COEF(g2)[k],DEG(g1)-k+1);
358: mulum(mod,g1,m,r); DEG(r) = divum(mod,r,g2,q);
359: l1 = dmb(mod,lc,lc,&tmp); l2 = dmb(mod,l,l1,&tmp);
360: DEG(m) = 0; COEF(m)[0] = l2;
361: divum(mod,r,m,q); cpyum(t,g1); cpyum(q,g2);
362: if ( d % 2 )
363: for ( i = DEG(g2); i >= 0; i-- )
364: COEF(g2)[i] = ( mod - COEF(g2)[i] ) % mod;
365: lc = COEF(g1)[DEG(g1)]; j = k - 1;
366: }
367: }
368: }
369: }
370:
371: void substmp(vl,mod,p,v0,p0,pr)
372: VL vl;
373: int mod;
374: V v0;
375: P p,p0;
376: P *pr;
377: {
378: P x,t,m,c,s,a;
379: DCP dc;
380: Q d;
381:
382: if ( !p )
383: *pr = 0;
384: else if ( NUM(p) )
385: *pr = p;
386: else if ( VR(p) != v0 ) {
387: MKMV(VR(p),x);
388: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
389: substmp(vl,mod,COEF(dc),v0,p0,&t);
390: if ( DEG(dc) ) {
391: pwrmp(vl,mod,x,DEG(dc),&s); mulmp(vl,mod,s,t,&m);
392: addmp(vl,mod,m,c,&a);
393: c = a;
394: } else {
395: addmp(vl,mod,t,c,&a);
396: c = a;
397: }
398: }
399: *pr = c;
400: } else {
401: dc = DC(p);
402: c = COEF(dc);
403: for ( d = DEG(dc), dc = NEXT(dc);
404: dc; d = DEG(dc), dc = NEXT(dc) ) {
405: subq(d,DEG(dc),(Q *)&t); pwrmp(vl,mod,p0,(Q)t,&s);
406: mulmp(vl,mod,s,c,&m);
407: addmp(vl,mod,m,COEF(dc),&c);
408: }
409: if ( d ) {
410: pwrmp(vl,mod,p0,d,&t); mulmp(vl,mod,t,c,&m);
411: c = m;
412: }
413: *pr = c;
414: }
415: }
416:
417: void reordermp(nvl,mod,ovl,p,pr)
418: VL nvl,ovl;
419: int mod;
420: P p;
421: P *pr;
422: {
423: DCP dc;
424: P x,m,s,t,c;
425:
426: if ( !p )
427: *pr = 0;
428: else if ( NUM(p) )
429: *pr = p;
430: else {
431: MKMV(VR(p),x);
432: for ( s = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
433: reordermp(nvl,mod,ovl,COEF(dc),&c);
434: if ( DEG(dc) ) {
435: pwrmp(nvl,mod,x,DEG(dc),&t); mulmp(nvl,mod,c,t,&m);
436: addmp(nvl,mod,m,s,&t);
437: } else
438: addmp(nvl,mod,s,c,&t);
439: s = t;
440: }
441: *pr = s;
442: }
443: }
444:
445: void chnremp(vl,mod,p,q,c,r)
446: VL vl;
447: int mod;
448: P p;
449: Q q;
450: P c;
451: P *r;
452: {
453: P tg,sg,ug;
454: P t,u;
455: MQ mq;
456:
457: ptomp(mod,p,&tg); submp(vl,mod,c,tg,&sg);
458: UTOMQ(rem(NM(q),mod),mq),tg = (P)mq; divsmp(vl,mod,sg,tg,&ug);
459: normalizemp(mod,ug);
460: mptop(ug,&u); mulp(vl,u,(P)q,&t); addp(vl,t,p,r);
461: }
462:
463: /* XXX strange behavior of invm() on SPARC */
464:
465: void chnrem(mod,v,c,q,t,cr,qr)
466: int mod;
467: V v;
468: UM t;
469: Q q,*qr;
470: P c,*cr;
471: {
472: int n,m,i,d,a,sd,tmp;
473: Q b,s,z;
474: Q *pc,*pcr;
475: DCP dc;
476:
477: if ( !c || NUM(c) )
478: n = 0;
479: else
480: n = UDEG(c);
481: m = DEG(t); d = MAX(n,m); W_CALLOC(n,Q,pc); W_CALLOC(d,Q,pcr);
482: if ( !c )
483: pc[0] = 0;
484: else if ( NUM(c) )
485: pc[0] = (Q)c;
486: else
487: for ( dc = DC(c); dc; dc = NEXT(dc) )
488: pc[QTOS(DEG(dc))] = (Q)COEF(dc);
489: for ( i = 0; i <= d; i++ ) {
490: b = (i>n?0:pc[i]); a = (i>m?0:COEF(t)[i]);
491: if ( b )
492: a = (a-SGN(pc[i])*((int)rem(NM(pc[i]),mod)))%mod;
493: sd = dmb(mod,(a>=0?a:a+mod),invm(rem(NM(q),mod),mod),&tmp);
494: if ( ( 2 * sd ) > mod )
495: sd -= mod;
496: STOQ(sd,z); mulq(z,q,&s); addq(s,b,&pcr[i]);
497: }
498: STOQ(mod,z); mulq(q,z,qr); plisttop((P *)pcr,v,d,cr);
499: }
500:
501: void normalizemp(mod,g)
502: int mod;
503: P g;
504: {
505: DCP dc;
506:
507: if ( !g )
508: return;
509: else if ( NUM(g) ) {
510: if ( 2 * CONT((MQ)g) > mod )
511: CONT((MQ)g) -= mod;
512: return;
513: } else
514: for ( dc = DC(g); dc; dc = NEXT(dc) )
515: normalizemp(mod,COEF(dc));
516: }
517:
518: void norm(p,r)
519: P p;
520: Q *r;
521: {
522: N t;
523: DCP dc;
524:
525: for ( dc = DC(p), t = ONEN; dc; dc = NEXT(dc) )
526: if ( cmpn(NM((Q)COEF(dc)),t) > 0 ) {
527: t = NM((Q)COEF(dc));
528: }
529: NTOQ(t,1,*r);
530: }
531:
532: void norm1(p,r)
533: P p,*r;
534: {
535: DCP dc;
536: P t,s,u;
537: Q q;
538:
539: if ( NUM(p) )
540: NTOQ(NM((Q)p),1,q),*r = (P)q;
541: else {
542: for ( t = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
543: norm1(COEF(dc),&s); addq((Q)t,(Q)s,(Q *)&u); t = u;
544: }
545: *r = t;
546: }
547: }
548:
549: void norm1c(p,r)
550: P p;
551: Q *r;
552: {
553: N t;
554: Q s;
555: DCP dc;
556:
557: if ( NUM(p) )
558: norm1(p,(P *)r);
559: else {
560: for ( dc = DC(p), t = ONEN; dc; dc = NEXT(dc) ) {
561: norm1(COEF(dc),(P *)&s);
562: if ( cmpn(NM(s),t) > 0 )
563: t = NM(s);
564: }
565: NTOQ(t,1,*r);
566: }
567: }
568:
569: void gcdprsmp(vl,mod,p1,p2,pr)
570: VL vl;
571: int mod;
572: P p1,p2,*pr;
573: {
574: P g1,g2,gc1,gc2,gp1,gp2,g,gc,gp,gcr;
575: V v1,v2;
576:
577: if ( !p1 )
578: *pr = p2;
579: else if ( !p2 )
580: *pr = p1;
581: else if ( NUM(p1) || NUM(p2) )
582: *pr = (P)ONEM;
583: else {
584: g1 = p1; g2 = p2;
585: if ( ( v1 = VR(g1) ) == ( v2 = VR(g2) ) ) {
586: gcdcmp(vl,mod,g1,&gc1); divsmp(vl,mod,g1,gc1,&gp1);
587: gcdcmp(vl,mod,g2,&gc2); divsmp(vl,mod,g2,gc2,&gp2);
588: gcdprsmp(vl,mod,gc1,gc2,&gcr);
589: sprsm(vl,mod,v1,gp1,gp2,&g);
590:
591: if ( VR(g) == v1 ) {
592: gp = g;
593: gcdcmp(vl,mod,gp,&gc); divsmp(vl,mod,gp,gc,&gp1);
594: mulmp(vl,mod,gp1,gcr,pr);
595: } else
596: *pr = gcr;
597: } else {
598: while ( v1 != vl->v && v2 != vl->v )
599: vl = NEXT(vl);
600: if ( v1 == vl->v ) {
601: gcdcmp(vl,mod,g1,&gc1); gcdprsmp(vl,mod,gc1,g2,pr);
602: } else {
603: gcdcmp(vl,mod,g2,&gc2); gcdprsmp(vl,mod,gc2,g1,pr);
604: }
605: }
606: }
607: }
608:
609: void gcdcmp(vl,mod,p,pr)
610: VL vl;
611: int mod;
612: P p,*pr;
613: {
614: P g,g1;
615: DCP dc;
616:
617: if ( NUM(p) )
618: *pr = (P)ONEM;
619: else {
620: dc = DC(p);
621: g = COEF(dc);
622: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
623: gcdprsmp(vl,mod,g,COEF(dc),&g1);
624: g = g1;
625: }
626: *pr = g;
627: }
628: }
629:
630: void sprsm(vl,mod,v,p1,p2,pr)
631: VL vl;
632: int mod;
633: V v;
634: P p1,p2,*pr;
635: {
636: P q1,q2,m,m1,m2,x,h,r,g1,g2;
637: int d;
638: Q dq;
639: VL nvl;
640:
641: reordvar(vl,v,&nvl);
642: reordermp(nvl,mod,vl,p1,&q1); reordermp(nvl,mod,vl,p2,&q2);
643:
644: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
645: *pr = 0;
646: return;
647: }
648:
649: if ( deg(v,q1) >= deg(v,q2) ) {
650: g1 = q1; g2 = q2;
651: } else {
652: g2 = q1; g1 = q2;
653: }
654:
655: for ( h = (P)ONEM, x = (P)ONEM; ; ) {
656: if ( !deg(v,g2) )
657: break;
658:
659: premmp(nvl,mod,g1,g2,&r);
660: if ( !r )
661: break;
662:
663: d = deg(v,g1) - deg(v,g2); STOQ(d,dq);
664: pwrmp(nvl,mod,h,dq,&m); mulmp(nvl,mod,m,x,&m1); g1 = g2;
665: divsmp(nvl,mod,r,m1,&g2); x = LC(g1); /* g1 is not const w.r.t v */
666: pwrmp(nvl,mod,x,dq,&m1); mulmp(nvl,mod,m1,h,&m2);
667: divsmp(nvl,mod,m2,m,&h);
668: }
669: *pr = g2;
670: }