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