![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to PU.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 PU.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/PU.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3:
4: void reorderp(nvl,ovl,p,pr)
5: VL nvl,ovl;
6: P p;
7: P *pr;
8: {
9: DCP dc;
10: P x,m,s,t,c;
11:
12: if ( !p )
13: *pr = 0;
14: else if ( NUM(p) )
15: *pr = p;
16: else {
17: MKV(VR(p),x);
18: for ( s = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
19: reorderp(nvl,ovl,COEF(dc),&c);
20: if ( DEG(dc) ) {
21: pwrp(nvl,x,DEG(dc),&t); mulp(nvl,c,t,&m);
22: addp(nvl,m,s,&t);
23: } else
24: addp(nvl,s,c,&t);
25: s = t;
26: }
27: *pr = s;
28: }
29: }
30:
31: void substp(vl,p,v0,p0,pr)
32: VL vl;
33: V v0;
34: P p,p0;
35: P *pr;
36: {
37: P x,t,m,c,s,a;
38: DCP dc;
39: Q d;
40:
41: if ( !p )
42: *pr = 0;
43: else if ( NUM(p) )
44: *pr = p;
45: else if ( VR(p) != v0 ) {
46: MKV(VR(p),x);
47: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
48: substp(vl,COEF(dc),v0,p0,&t);
49: if ( DEG(dc) ) {
50: pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
51: addp(vl,m,c,&a);
52: c = a;
53: } else {
54: addp(vl,t,c,&a);
55: c = a;
56: }
57: }
58: *pr = c;
59: } else {
60: dc = DC(p);
61: c = COEF(dc);
62: for ( d = DEG(dc), dc = NEXT(dc);
63: dc; d = DEG(dc), dc = NEXT(dc) ) {
64: subq(d,DEG(dc),(Q *)&t); pwrp(vl,p0,(Q)t,&s);
65: mulp(vl,s,c,&m);
66: addp(vl,m,COEF(dc),&c);
67: }
68: if ( d ) {
69: pwrp(vl,p0,d,&t); mulp(vl,t,c,&m);
70: c = m;
71: }
72: *pr = c;
73: }
74: }
75:
76: void detp(vl,rmat,n,dp)
77: VL vl;
78: P **rmat;
79: int n;
80: P *dp;
81: {
82: int i,j,k,sgn;
83: P mjj,mij,t,s,u,d;
84: P **mat;
85: P *mi,*mj;
86:
87: mat = (P **)almat_pointer(n,n);
88: for ( i = 0; i < n; i++ )
89: for ( j = 0; j < n; j++ )
90: mat[i][j] = rmat[i][j];
91: for ( j = 0, d = (P)ONE, sgn = 1; j < n; j++ ) {
92: for ( i = j; (i < n) && !mat[i][j]; i++ );
93: if ( i == n ) {
94: *dp = 0; return;
95: }
96: for ( k = i; k < n; k++ )
97: if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
98: i = k;
99: if ( j != i ) {
100: mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
101: }
102: for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
103: for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
104: mulp(vl,mi[k],mjj,&t); mulp(vl,mj[k],mij,&s);
105: subp(vl,t,s,&u); divsp(vl,u,d,&mi[k]);
106: }
107: d = mjj;
108: }
109: if ( sgn < 0 )
110: chsgnp(d,dp);
111: else
112: *dp = d;
113: }
114:
115: void reordvar(vl,v,nvlp)
116: VL vl;
117: V v;
118: VL *nvlp;
119: {
120: VL nvl,nvl0;
121:
122: for ( NEWVL(nvl0), nvl0->v = v, nvl = nvl0;
123: vl; vl = NEXT(vl) )
124: if ( vl->v == v )
125: continue;
126: else {
127: NEWVL(NEXT(nvl));
128: nvl = NEXT(nvl);
129: nvl->v = vl->v;
130: }
131: NEXT(nvl) = 0;
132: *nvlp = nvl0;
133: }
134:
135: void gcdprsp(vl,p1,p2,pr)
136: VL vl;
137: P p1,p2,*pr;
138: {
139: P g1,g2,gc1,gc2,gp1,gp2,g,gc,gp,gcr;
140: V v1,v2;
141:
142: if ( !p1 )
143: ptozp0(p2,pr);
144: else if ( !p2 )
145: ptozp0(p1,pr);
146: else if ( NUM(p1) || NUM(p2) )
147: *pr = (P)ONE;
148: else {
149: ptozp0(p1,&g1); ptozp0(p2,&g2);
150: if ( ( v1 = VR(g1) ) == ( v2 = VR(g2) ) ) {
151: gcdcp(vl,g1,&gc1); divsp(vl,g1,gc1,&gp1);
152: gcdcp(vl,g2,&gc2); divsp(vl,g2,gc2,&gp2);
153: gcdprsp(vl,gc1,gc2,&gcr);
154: sprs(vl,v1,gp1,gp2,&g);
155:
156: if ( VR(g) == v1 ) {
157: ptozp0(g,&gp);
158: gcdcp(vl,gp,&gc); divsp(vl,gp,gc,&gp1);
159: mulp(vl,gp1,gcr,pr);
160:
161: } else
162: *pr = gcr;
163: } else {
164: while ( v1 != vl->v && v2 != vl->v )
165: vl = NEXT(vl);
166: if ( v1 == vl->v ) {
167: gcdcp(vl,g1,&gc1); gcdprsp(vl,gc1,g2,pr);
168: } else {
169: gcdcp(vl,g2,&gc2); gcdprsp(vl,gc2,g1,pr);
170: }
171: }
172: }
173: }
174:
175: void gcdcp(vl,p,pr)
176: VL vl;
177: P p,*pr;
178: {
179: P g,g1;
180: DCP dc;
181:
182: if ( NUM(p) )
183: *pr = (P)ONE;
184: else {
185: dc = DC(p);
186: g = COEF(dc);
187: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
188: gcdprsp(vl,g,COEF(dc),&g1);
189: g = g1;
190: }
191: *pr = g;
192: }
193: }
194:
195: void sprs(vl,v,p1,p2,pr)
196: VL vl;
197: V v;
198: P p1,p2,*pr;
199: {
200: P q1,q2,m,m1,m2,x,h,r,g1,g2;
201: int d;
202: Q dq;
203: VL nvl;
204:
205: reordvar(vl,v,&nvl);
206: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
207:
208: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
209: *pr = 0;
210: return;
211: }
212:
213: if ( deg(v,q1) >= deg(v,q2) ) {
214: g1 = q1; g2 = q2;
215: } else {
216: g2 = q1; g1 = q2;
217: }
218:
219: for ( h = (P)ONE, x = (P)ONE; ; ) {
220: if ( !deg(v,g2) )
221: break;
222:
223: premp(nvl,g1,g2,&r);
224: if ( !r )
225: break;
226:
227: d = deg(v,g1) - deg(v,g2); STOQ(d,dq);
228: pwrp(nvl,h,dq,&m); mulp(nvl,m,x,&m1); g1 = g2;
229: divsp(nvl,r,m1,&g2); x = LC(g1); /* g1 is not const w.r.t v */
230: pwrp(nvl,x,dq,&m1); mulp(nvl,m1,h,&m2);
231: divsp(nvl,m2,m,&h);
232: }
233: *pr = g2;
234: }
235:
236: void resultp(vl,v,p1,p2,pr)
237: VL vl;
238: V v;
239: P p1,p2,*pr;
240: {
241: P q1,q2,m,m1,m2,lc,q,r,t,g1,g2,adj;
242: int d,d1,d2,j,k;
243: VL nvl;
244: Q dq;
245:
246: if ( !p1 || !p2 ) {
247: *pr = 0; return;
248: }
249: reordvar(vl,v,&nvl);
250: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
251:
252: if ( VR(q1) != v )
253: if ( VR(q2) != v ) {
254: *pr = 0;
255: return;
256: } else {
257: d = deg(v,q2); STOQ(d,dq);
258: pwrp(vl,q1,dq,pr);
259: return;
260: }
261: else if ( VR(q2) != v ) {
262: d = deg(v,q1); STOQ(d,dq);
263: pwrp(vl,q2,dq,pr);
264: return;
265: }
266:
267: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
268: *pr = 0;
269: return;
270: }
271:
272: d1 = deg(v,q1); d2 = deg(v,q2);
273: if ( d1 > d2 ) {
274: g1 = q1; g2 = q2; adj = (P)ONE;
275: } else if ( d1 < d2 ) {
276: g2 = q1; g1 = q2;
277: if ( (d1 % 2) && (d2 % 2) ) {
278: STOQ(-1,dq); adj = (P)dq;
279: } else
280: adj = (P)ONE;
281: } else {
282: premp(nvl,q1,q2,&t);
283: d = deg(v,t); STOQ(d,dq); pwrp(nvl,LC(q2),dq,&adj);
284: g1 = q2; g2 = t;
285: if ( d1 % 2 ) {
286: chsgnp(adj,&t); adj = t;
287: }
288: }
289: d1 = deg(v,g1); j = d1 - 1;
290:
291: for ( lc = (P)ONE; ; ) {
292: if ( ( k = deg(v,g2) ) < 0 ) {
293: *pr = 0;
294: return;
295: }
296:
297: if ( k == j )
298: if ( !k ) {
299: divsp(nvl,g2,adj,pr);
300: return;
301: } else {
302: premp(nvl,g1,g2,&r); mulp(nvl,lc,lc,&m);
303: divsp(nvl,r,m,&q);
304: g1 = g2; g2 = q;
305: lc = LC(g1); /* g1 is not const */
306: j = k - 1;
307: }
308: else {
309: d = j - k; STOQ(d,dq);
310: pwrp(nvl,(VR(g2)==v?LC(g2):g2),dq,&m);
311: mulp(nvl,g2,m,&m1);
312: pwrp(nvl,lc,dq,&m); divsp(nvl,m1,m,&t);
313: if ( k == 0 ) {
314: divsp(nvl,t,adj,pr);
315: return;
316: } else {
317: premp(nvl,g1,g2,&r);
318: mulp(nvl,lc,lc,&m1); mulp(nvl,m,m1,&m2);
319: divsp(nvl,r,m2,&q);
320: g1 = t; g2 = q;
321: if ( d % 2 ) {
322: chsgnp(g2,&t); g2 = t;
323: }
324: lc = LC(g1); /* g1 is not const */
325: j = k - 1;
326: }
327: }
328: }
329: }
330:
331: void srch2(vl,v,p1,p2,pr)
332: VL vl;
333: V v;
334: P p1,p2,*pr;
335: {
336: P q1,q2,m,m1,m2,lc,q,r,t,s,g1,g2,adj;
337: int d,d1,d2,j,k;
338: VL nvl;
339: Q dq;
340:
341: reordvar(vl,v,&nvl);
342: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
343:
344: if ( VR(q1) != v )
345: if ( VR(q2) != v ) {
346: *pr = 0;
347: return;
348: } else {
349: d = deg(v,q2); STOQ(d,dq);
350: pwrp(vl,q1,dq,pr);
351: return;
352: }
353: else if ( VR(q2) != v ) {
354: d = deg(v,q1); STOQ(d,dq);
355: pwrp(vl,q2,dq,pr);
356: return;
357: }
358:
359: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
360: *pr = 0;
361: return;
362: }
363:
364: if ( deg(v,q1) >= deg(v,q2) ) {
365: g1 = q1; g2 = q2;
366: } else {
367: g2 = q1; g1 = q2;
368: }
369:
370:
371: if ( ( d1 = deg(v,g1) ) > ( d2 = deg(v,g2) ) ) {
372: j = d1 - 1;
373: adj = (P)ONE;
374: } else {
375: premp(nvl,g1,g2,&t);
376: d = deg(v,t); STOQ(d,dq);
377: pwrp(nvl,LC(g2),dq,&adj);
378: g1 = g2; g2 = t;
379: j = deg(v,g1) - 1;
380: }
381:
382: for ( lc = (P)ONE; ; ) {
383: if ( ( k = deg(v,g2) ) < 0 ) {
384: *pr = 0;
385: return;
386: }
387:
388: ptozp(g1,1,(Q *)&t,&s); g1 = s; ptozp(g2,1,(Q *)&t,&s); g2 = s;
389: if ( k == j )
390: if ( !k ) {
391: divsp(nvl,g2,adj,pr);
392: return;
393: } else {
394: premp(nvl,g1,g2,&r); mulp(nvl,lc,lc,&m);
395: divsp(nvl,r,m,&q);
396: g1 = g2; g2 = q;
397: lc = LC(g1); /* g1 is not const */
398: j = k - 1;
399: }
400: else {
401: d = j - k; STOQ(d,dq);
402: pwrp(nvl,(VR(g2)==v?LC(g2):g2),dq,&m);
403: mulp(nvl,g2,m,&m1);
404: pwrp(nvl,lc,dq,&m); divsp(nvl,m1,m,&t);
405: if ( k == 0 ) {
406: divsp(nvl,t,adj,pr);
407: return;
408: } else {
409: premp(nvl,g1,g2,&r);
410: mulp(nvl,lc,lc,&m1); mulp(nvl,m,m1,&m2);
411: divsp(nvl,r,m2,&q);
412: g1 = t; g2 = q;
413: if ( d % 2 ) {
414: chsgnp(g2,&t); g2 = t;
415: }
416: lc = LC(g1); /* g1 is not const */
417: j = k - 1;
418: }
419: }
420: }
421: }
422:
423: void srcr(vl,v,p1,p2,pr)
424: VL vl;
425: V v;
426: P p1,p2,*pr;
427: {
428: P q1,q2,c,c1;
429: P tg,tg1,tg2,resg;
430: int index,mod;
431: Q a,b,f,q,t,s,u,n,m;
432: VL nvl;
433:
434: reordvar(vl,v,&nvl);
435: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
436:
437: if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
438: *pr = 0;
439: return;
440: }
441: if ( VR(q1) != v ) {
442: pwrp(vl,q1,DEG(DC(q2)),pr);
443: return;
444: }
445: if ( VR(q2) != v ) {
446: pwrp(vl,q2,DEG(DC(q1)),pr);
447: return;
448: }
449: norm1c(q1,&a); norm1c(q2,&b);
450: n = DEG(DC(q1)); m = DEG(DC(q2));
451: pwrq(a,m,&t); pwrq(b,n,&s); mulq(t,s,&u);
452: factorial(QTOS(n)+QTOS(m),&t);
453: mulq(u,t,&s); addq(s,s,&f);
454: for ( index = 0, q = ONE, c = 0; cmpq(f,q) >= 0; ) {
455: mod = lprime[index++];
456: if ( !mod )
457: error("sqfrum : lprime[] exhausted.");
458: ptomp(mod,LC(q1),&tg);
459: if ( !tg )
460: continue;
461: ptomp(mod,LC(q2),&tg);
462: if ( !tg )
463: continue;
464: ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
465: srchmp(nvl,mod,v,tg1,tg2,&resg);
466: chnremp(nvl,mod,c,q,resg,&c1); c = c1;
467: STOQ(mod,t); mulq(q,t,&s); q = s;
468: }
469: *pr = c;
470: }
471:
472: void res_ch_det(vl,v,p1,p2,pr)
473: VL vl;
474: V v;
475: P p1,p2,*pr;
476: {
477: P q1,q2,c,c1;
478: P tg,tg1,tg2,resg;
479: int index,mod;
480: Q a,b,f,q,t,s,u,n,m;
481: VL nvl;
482:
483: reordvar(vl,v,&nvl);
484: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
485:
486: if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
487: *pr = 0;
488: return;
489: }
490: if ( VR(q1) != v ) {
491: pwrp(vl,q1,DEG(DC(q2)),pr);
492: return;
493: }
494: if ( VR(q2) != v ) {
495: pwrp(vl,q2,DEG(DC(q1)),pr);
496: return;
497: }
498: norm1c(q1,&a); norm1c(q2,&b);
499: n = DEG(DC(q1)); m = DEG(DC(q2));
500: pwrq(a,m,&t); pwrq(b,n,&s); mulq(t,s,&u);
501: factorial(QTOS(n)+QTOS(m),&t);
502: mulq(u,t,&s); addq(s,s,&f);
503: for ( index = 0, q = ONE, c = 0; cmpq(f,q) >= 0; ) {
504: mod = lprime[index++];
505: if ( !mod )
506: error("sqfrum : lprime[] exhausted.");
507: ptomp(mod,LC(q1),&tg);
508: if ( !tg )
509: continue;
510: ptomp(mod,LC(q2),&tg);
511: if ( !tg )
512: continue;
513: ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
514: res_detmp(nvl,mod,v,tg1,tg2,&resg);
515: chnremp(nvl,mod,c,q,resg,&c1); c = c1;
516: STOQ(mod,t); mulq(q,t,&s); q = s;
517: }
518: *pr = c;
519: }
520:
521: void res_detmp(vl,mod,v,p1,p2,dp)
522: VL vl;
523: int mod;
524: V v;
525: P p1,p2;
526: P *dp;
527: {
528: int n1,n2,n,sgn;
529: int i,j,k;
530: P mjj,mij,t,s,u,d;
531: P **mat;
532: P *mi,*mj;
533: DCP dc;
534:
535: n1 = UDEG(p1); n2 = UDEG(p2); n = n1+n2;
536: mat = (P **)almat_pointer(n,n);
537: for ( dc = DC(p1); dc; dc = NEXT(dc) )
538: mat[0][n1-QTOS(DEG(dc))] = COEF(dc);
539: for ( i = 1; i < n2; i++ )
540: for ( j = 0; j <= n1; j++ )
541: mat[i][i+j] = mat[0][j];
542: for ( dc = DC(p2); dc; dc = NEXT(dc) )
543: mat[n2][n2-QTOS(DEG(dc))] = COEF(dc);
544: for ( i = 1; i < n1; i++ )
545: for ( j = 0; j <= n2; j++ )
546: mat[i+n2][i+j] = mat[n2][j];
547: for ( j = 0, d = (P)ONEM, sgn = 1; j < n; j++ ) {
548: for ( i = j; (i < n) && !mat[i][j]; i++ );
549: if ( i == n ) {
550: *dp = 0; return;
551: }
552: for ( k = i; k < n; k++ )
553: if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
554: i = k;
555: if ( j != i ) {
556: mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
557: }
558: for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
559: for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
560: mulmp(vl,mod,mi[k],mjj,&t); mulmp(vl,mod,mj[k],mij,&s);
561: submp(vl,mod,t,s,&u); divsmp(vl,mod,u,d,&mi[k]);
562: }
563: d = mjj;
564: }
565: if ( sgn > 0 )
566: *dp = d;
567: else
568: chsgnmp(mod,d,dp);
569: }
570:
571: #if 0
572: showmat(vl,mat,n)
573: VL vl;
574: P **mat;
575: int n;
576: {
577: int i,j;
578: P t;
579:
580: for ( i = 0; i < n; i++ ) {
581: for ( j = 0; j < n; j++ ) {
582: mptop(mat[i][j],&t); printp(vl,t); fprintf(out," ");
583: }
584: fprintf(out,"\n");
585: }
586: fflush(out);
587: }
588:
589: showmp(vl,p)
590: VL vl;
591: P p;
592: {
593: P t;
594:
595: mptop(p,&t); printp(vl,t); fprintf(out,"\n");
596: }
597: #endif
598:
599: void premp(vl,p1,p2,pr)
600: VL vl;
601: P p1,p2,*pr;
602: {
603: P m,m1,m2;
604: P *pw;
605: DCP dc;
606: V v1,v2;
607: register int i,j;
608: int n1,n2,d;
609:
610: if ( NUM(p1) )
611: if ( NUM(p2) )
612: *pr = 0;
613: else
614: *pr = p1;
615: else if ( NUM(p2) )
616: *pr = 0;
617: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
618: n1 = deg(v1,p1); n2 = deg(v1,p2);
619: pw = (P *)ALLOCA((n1+1)*sizeof(P));
620: bzero((char *)pw,(int)((n1+1)*sizeof(P)));
621:
622: for ( dc = DC(p1); dc; dc = NEXT(dc) )
623: pw[QTOS(DEG(dc))] = COEF(dc);
624:
625: for ( i = n1; i >= n2; i-- ) {
626: if ( pw[i] ) {
627: m = pw[i];
628: for ( j = i; j >= 0; j-- ) {
629: mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
630: }
631:
632: for ( dc = DC(p2), d = i - n2; dc; dc = NEXT(dc) ) {
633: mulp(vl,COEF(dc),m,&m1);
634: subp(vl,pw[QTOS(DEG(dc))+d],m1,&m2);
635: pw[QTOS(DEG(dc))+d] = m2;
636: }
637: } else
638: for ( j = i; j >= 0; j-- )
639: if ( pw[j] ) {
640: mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
641: }
642: }
643: plisttop(pw,v1,n2-1,pr);
644: } else {
645: while ( v1 != vl->v && v2 != vl->v )
646: vl = NEXT(vl);
647: if ( v1 == vl->v )
648: *pr = 0;
649: else
650: *pr = p1;
651: }
652: }
653:
654: void ptozp0(p,pr)
655: P p;
656: P *pr;
657: {
658: Q c;
659:
660: ptozp(p,1,&c,pr);
661: }
662:
663: void mindegp(vl,p,mvlp,pr)
664: VL vl,*mvlp;
665: P p,*pr;
666: {
667: P t;
668: VL nvl,tvl,avl;
669: V v;
670: int n,d;
671:
672: clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
673: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
674: n = getdeg(avl->v,p);
675: if ( n < d ) {
676: v = avl->v; d = n;
677: }
678: }
679: if ( v != nvl->v ) {
680: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
681: *pr = t; *mvlp = tvl;
682: } else {
683: *pr = p; *mvlp = nvl;
684: }
685: }
686:
687: void maxdegp(vl,p,mvlp,pr)
688: VL vl,*mvlp;
689: P p,*pr;
690: {
691: P t;
692: VL nvl,tvl,avl;
693: V v;
694: int n,d;
695:
696: clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
697: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
698: n = getdeg(avl->v,p);
699: if ( n > d ) {
700: v = avl->v; d = n;
701: }
702: }
703: if ( v != nvl->v ) {
704: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
705: *pr = t; *mvlp = tvl;
706: } else {
707: *pr = p; *mvlp = nvl;
708: }
709: }
710:
711: void min_common_vars_in_coefp(vl,p,mvlp,pr)
712: VL vl,*mvlp;
713: P p,*pr;
714: {
715: P u,p0;
716: VL tvl,cvl,svl,uvl,avl,vl0;
717: int n,n0,d,d0;
718: DCP dc;
719:
720: clctv(vl,p,&tvl); vl = tvl;
721: for ( tvl = vl, n0 = 0; tvl; tvl = NEXT(tvl), n0++ );
722: for ( avl = vl; avl; avl = NEXT(avl) ) {
723: if ( VR(p) != avl->v ) {
724: reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
725: } else {
726: uvl = vl; u = p;
727: }
728: for ( cvl = NEXT(uvl), dc = DC(u); dc && cvl; dc = NEXT(dc) ) {
729: clctv(uvl,COEF(dc),&tvl); intersectv(cvl,tvl,&svl); cvl = svl;
730: }
731: if ( !cvl ) {
732: *mvlp = uvl; *pr = u; return;
733: } else {
734: for ( tvl = cvl, n = 0; tvl; tvl = NEXT(tvl), n++ );
735: if ( n < n0 ) {
736: vl0 = uvl; p0 = u; n0 = n; d0 = getdeg(uvl->v,u);
737: } else if ( (n == n0) && ((d = getdeg(uvl->v,u)) < d0) ) {
738: vl0 = uvl; p0 = u; n0 = n; d0 = d;
739: }
740: }
741: }
742: *pr = p0; *mvlp = vl0;
743: }
744:
745: void minlcdegp(vl,p,mvlp,pr)
746: VL vl,*mvlp;
747: P p,*pr;
748: {
749: P u,p0;
750: VL tvl,uvl,avl,vl0;
751: int d,d0;
752:
753: clctv(vl,p,&tvl); vl = tvl;
754: vl0 = vl; p0 = p; d0 = homdeg(COEF(DC(p)));
755: for ( avl = NEXT(vl); avl; avl = NEXT(avl) ) {
756: reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
757: d = homdeg(COEF(DC(u)));
758: if ( d < d0 ) {
759: vl0 = uvl; p0 = u; d0 = d;
760: }
761: }
762: *pr = p0; *mvlp = vl0;
763: }
764:
765: void sort_by_deg(n,p,pr)
766: int n;
767: P *p,*pr;
768: {
769: int j,k,d,k0;
770: V v;
771:
772: for ( j = 0; j < n; j++ ) {
773: for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
774: k < n; k++ )
775: if ( deg(v,p[k]) < d ) {
776: k0 = k;
777: d = deg(v,p[k]);
778: }
779: pr[j] = p[k0];
780: if ( j != k0 )
781: p[k0] = p[j];
782: }
783: }
784:
785: void sort_by_deg_rev(n,p,pr)
786: int n;
787: P *p,*pr;
788: {
789: int j,k,d,k0;
790: V v;
791:
792: for ( j = 0; j < n; j++ ) {
793: for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
794: k < n; k++ )
795: if ( deg(v,p[k]) > d ) {
796: k0 = k;
797: d = deg(v,p[k]);
798: }
799: pr[j] = p[k0];
800: if ( j != k0 )
801: p[k0] = p[j];
802: }
803: }
804:
805:
806: void getmindeg(v,p,dp)
807: V v;
808: P p;
809: Q *dp;
810: {
811: Q dt,d;
812: DCP dc;
813:
814: if ( !p || NUM(p) )
815: *dp = 0;
816: else if ( v == VR(p) ) {
817: for ( dc = DC(p); NEXT(dc); dc = NEXT(dc) );
818: *dp = DEG(dc);
819: } else {
820: dc = DC(p);
821: getmindeg(v,COEF(dc),&d);
822: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
823: getmindeg(v,COEF(dc),&dt);
824: if ( cmpq(dt,d) < 0 )
825: d = dt;
826: }
827: *dp = d;
828: }
829: }
830:
831: void minchdegp(vl,p,mvlp,pr)
832: VL vl,*mvlp;
833: P p,*pr;
834: {
835: P t;
836: VL tvl,nvl,avl;
837: int n,d,z;
838: V v;
839:
840: if ( NUM(p) ) {
841: *mvlp = vl;
842: *pr = p;
843: return;
844: }
845: clctv(vl,p,&nvl);
846: v = nvl->v;
847: d = getchomdeg(v,p) + getlchomdeg(v,p,&z);
848: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
849: n = getchomdeg(avl->v,p) + getlchomdeg(avl->v,p,&z);
850: if ( n < d ) {
851: v = avl->v; d = n;
852: }
853: }
854: if ( v != nvl->v ) {
855: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
856: *pr = t; *mvlp = tvl;
857: } else {
858: *pr = p; *mvlp = nvl;
859: }
860: }
861:
862: int getchomdeg(v,p)
863: V v;
864: P p;
865: {
866: int m,m1;
867: DCP dc;
868:
869: if ( !p || NUM(p) )
870: return ( 0 );
871: else if ( VR(p) == v )
872: return ( dbound(v,p) );
873: else {
874: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
875: m1 = getchomdeg(v,COEF(dc))+QTOS(DEG(dc));
876: m = MAX(m,m1);
877: }
878: return ( m );
879: }
880: }
881:
882: int getlchomdeg(v,p,d)
883: V v;
884: P p;
885: int *d;
886: {
887: int m0,m1,d0,d1;
888: DCP dc;
889:
890: if ( !p || NUM(p) ) {
891: *d = 0;
892: return ( 0 );
893: } else if ( VR(p) == v ) {
894: *d = QTOS(DEG(DC(p)));
895: return ( homdeg(LC(p)) );
896: } else {
897: for ( dc = DC(p), m0 = 0, d0 = 0; dc; dc = NEXT(dc) ) {
898: m1 = getlchomdeg(v,COEF(dc),&d1)+QTOS(DEG(dc));
899: if ( d1 > d0 ) {
900: m0 = m1;
901: d0 = d1;
902: } else if ( d1 == d0 )
903: m0 = MAX(m0,m1);
904: }
905: *d = d0;
906: return ( m0 );
907: }
908: }
909:
910: int nmonop(p)
911: P p;
912: {
913: int s;
914: DCP dc;
915:
916: if ( !p )
917: return 0;
918: else
919: switch ( OID(p) ) {
920: case O_N:
921: return 1; break;
922: case O_P:
923: for ( dc = DC((P)p), s = 0; dc; dc = NEXT(dc) )
924: s += nmonop(COEF(dc));
925: return s; break;
926: default:
927: return 0;
928: }
929: }
930:
931: int qpcheck(p)
932: Obj p;
933: {
934: DCP dc;
935:
936: if ( !p )
937: return 1;
938: else
939: switch ( OID(p) ) {
940: case O_N:
941: return RATN((Num)p)?1:0;
942: case O_P:
943: for ( dc = DC((P)p); dc; dc = NEXT(dc) )
944: if ( !qpcheck((Obj)COEF(dc)) )
945: return 0;
946: return 1;
947: default:
948: return 0;
949: }
950: }
951:
952: /* check if p is univariate and all coeffs are INT or LM */
953:
954: int uzpcheck(p)
955: Obj p;
956: {
957: DCP dc;
958: P c;
959:
960: if ( !p )
961: return 1;
962: else
963: switch ( OID(p) ) {
964: case O_N:
965: return (RATN((Num)p)&&INT(p))||(NID((Num)p)==N_LM);
966: case O_P:
967: for ( dc = DC((P)p); dc; dc = NEXT(dc) ) {
968: c = COEF(dc);
969: if ( !NUM(c) || !uzpcheck(c) )
970: return 0;
971: }
972: return 1;
973: default:
974: return 0;
975: }
976: }
977:
978: int p_mag(p)
979: P p;
980: {
981: int s;
982: DCP dc;
983:
984: if ( !p )
985: return 0;
986: else if ( OID(p) == O_N )
987: return n_bits(NM((Q)p))+(INT((Q)p)?0:n_bits(DN((Q)p)));
988: else {
989: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) )
990: s += p_mag(COEF(dc));
991: return s;
992: }
993: }