![[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.4 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/PU.c,v 1.3 2000/08/21 08:31:26 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:
708: ptozp(p,1,&c,pr);
709: }
710:
711: void mindegp(vl,p,mvlp,pr)
712: VL vl,*mvlp;
713: P p,*pr;
714: {
715: P t;
716: VL nvl,tvl,avl;
717: V v;
718: int n,d;
719:
720: clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
721: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
722: n = getdeg(avl->v,p);
723: if ( n < d ) {
724: v = avl->v; d = n;
725: }
726: }
727: if ( v != nvl->v ) {
728: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
729: *pr = t; *mvlp = tvl;
730: } else {
731: *pr = p; *mvlp = nvl;
732: }
733: }
734:
735: void maxdegp(vl,p,mvlp,pr)
736: VL vl,*mvlp;
737: P p,*pr;
738: {
739: P t;
740: VL nvl,tvl,avl;
741: V v;
742: int n,d;
743:
744: clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
745: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
746: n = getdeg(avl->v,p);
747: if ( n > d ) {
748: v = avl->v; d = n;
749: }
750: }
751: if ( v != nvl->v ) {
752: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
753: *pr = t; *mvlp = tvl;
754: } else {
755: *pr = p; *mvlp = nvl;
756: }
757: }
758:
759: void min_common_vars_in_coefp(vl,p,mvlp,pr)
760: VL vl,*mvlp;
761: P p,*pr;
762: {
763: P u,p0;
764: VL tvl,cvl,svl,uvl,avl,vl0;
765: int n,n0,d,d0;
766: DCP dc;
767:
768: clctv(vl,p,&tvl); vl = tvl;
769: for ( tvl = vl, n0 = 0; tvl; tvl = NEXT(tvl), n0++ );
770: for ( avl = vl; avl; avl = NEXT(avl) ) {
771: if ( VR(p) != avl->v ) {
772: reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
773: } else {
774: uvl = vl; u = p;
775: }
776: for ( cvl = NEXT(uvl), dc = DC(u); dc && cvl; dc = NEXT(dc) ) {
777: clctv(uvl,COEF(dc),&tvl); intersectv(cvl,tvl,&svl); cvl = svl;
778: }
779: if ( !cvl ) {
780: *mvlp = uvl; *pr = u; return;
781: } else {
782: for ( tvl = cvl, n = 0; tvl; tvl = NEXT(tvl), n++ );
783: if ( n < n0 ) {
784: vl0 = uvl; p0 = u; n0 = n; d0 = getdeg(uvl->v,u);
785: } else if ( (n == n0) && ((d = getdeg(uvl->v,u)) < d0) ) {
786: vl0 = uvl; p0 = u; n0 = n; d0 = d;
787: }
788: }
789: }
790: *pr = p0; *mvlp = vl0;
791: }
792:
793: void minlcdegp(vl,p,mvlp,pr)
794: VL vl,*mvlp;
795: P p,*pr;
796: {
797: P u,p0;
798: VL tvl,uvl,avl,vl0;
799: int d,d0;
800:
801: clctv(vl,p,&tvl); vl = tvl;
802: vl0 = vl; p0 = p; d0 = homdeg(COEF(DC(p)));
803: for ( avl = NEXT(vl); avl; avl = NEXT(avl) ) {
804: reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
805: d = homdeg(COEF(DC(u)));
806: if ( d < d0 ) {
807: vl0 = uvl; p0 = u; d0 = d;
808: }
809: }
810: *pr = p0; *mvlp = vl0;
811: }
812:
813: void sort_by_deg(n,p,pr)
814: int n;
815: P *p,*pr;
816: {
817: int j,k,d,k0;
818: V v;
819:
820: for ( j = 0; j < n; j++ ) {
821: for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
822: k < n; k++ )
823: if ( deg(v,p[k]) < d ) {
824: k0 = k;
825: d = deg(v,p[k]);
826: }
827: pr[j] = p[k0];
828: if ( j != k0 )
829: p[k0] = p[j];
830: }
831: }
832:
833: void sort_by_deg_rev(n,p,pr)
834: int n;
835: P *p,*pr;
836: {
837: int j,k,d,k0;
838: V v;
839:
840: for ( j = 0; j < n; j++ ) {
841: for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
842: k < n; k++ )
843: if ( deg(v,p[k]) > d ) {
844: k0 = k;
845: d = deg(v,p[k]);
846: }
847: pr[j] = p[k0];
848: if ( j != k0 )
849: p[k0] = p[j];
850: }
851: }
852:
853:
854: void getmindeg(v,p,dp)
855: V v;
856: P p;
857: Q *dp;
858: {
859: Q dt,d;
860: DCP dc;
861:
862: if ( !p || NUM(p) )
863: *dp = 0;
864: else if ( v == VR(p) ) {
865: for ( dc = DC(p); NEXT(dc); dc = NEXT(dc) );
866: *dp = DEG(dc);
867: } else {
868: dc = DC(p);
869: getmindeg(v,COEF(dc),&d);
870: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
871: getmindeg(v,COEF(dc),&dt);
872: if ( cmpq(dt,d) < 0 )
873: d = dt;
874: }
875: *dp = d;
876: }
877: }
878:
879: void minchdegp(vl,p,mvlp,pr)
880: VL vl,*mvlp;
881: P p,*pr;
882: {
883: P t;
884: VL tvl,nvl,avl;
885: int n,d,z;
886: V v;
887:
888: if ( NUM(p) ) {
889: *mvlp = vl;
890: *pr = p;
891: return;
892: }
893: clctv(vl,p,&nvl);
894: v = nvl->v;
895: d = getchomdeg(v,p) + getlchomdeg(v,p,&z);
896: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
897: n = getchomdeg(avl->v,p) + getlchomdeg(avl->v,p,&z);
898: if ( n < d ) {
899: v = avl->v; d = n;
900: }
901: }
902: if ( v != nvl->v ) {
903: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
904: *pr = t; *mvlp = tvl;
905: } else {
906: *pr = p; *mvlp = nvl;
907: }
908: }
909:
910: int getchomdeg(v,p)
911: V v;
912: P p;
913: {
914: int m,m1;
915: DCP dc;
916:
917: if ( !p || NUM(p) )
918: return ( 0 );
919: else if ( VR(p) == v )
920: return ( dbound(v,p) );
921: else {
922: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
923: m1 = getchomdeg(v,COEF(dc))+QTOS(DEG(dc));
924: m = MAX(m,m1);
925: }
926: return ( m );
927: }
928: }
929:
930: int getlchomdeg(v,p,d)
931: V v;
932: P p;
933: int *d;
934: {
935: int m0,m1,d0,d1;
936: DCP dc;
937:
938: if ( !p || NUM(p) ) {
939: *d = 0;
940: return ( 0 );
941: } else if ( VR(p) == v ) {
942: *d = QTOS(DEG(DC(p)));
943: return ( homdeg(LC(p)) );
944: } else {
945: for ( dc = DC(p), m0 = 0, d0 = 0; dc; dc = NEXT(dc) ) {
946: m1 = getlchomdeg(v,COEF(dc),&d1)+QTOS(DEG(dc));
947: if ( d1 > d0 ) {
948: m0 = m1;
949: d0 = d1;
950: } else if ( d1 == d0 )
951: m0 = MAX(m0,m1);
952: }
953: *d = d0;
954: return ( m0 );
955: }
956: }
957:
958: int nmonop(p)
959: P p;
960: {
961: int s;
962: DCP dc;
963:
964: if ( !p )
965: return 0;
966: else
967: switch ( OID(p) ) {
968: case O_N:
969: return 1; break;
970: case O_P:
971: for ( dc = DC((P)p), s = 0; dc; dc = NEXT(dc) )
972: s += nmonop(COEF(dc));
973: return s; break;
974: default:
975: return 0;
976: }
977: }
978:
979: int qpcheck(p)
980: Obj p;
981: {
982: DCP dc;
983:
984: if ( !p )
985: return 1;
986: else
987: switch ( OID(p) ) {
988: case O_N:
989: return RATN((Num)p)?1:0;
990: case O_P:
991: for ( dc = DC((P)p); dc; dc = NEXT(dc) )
992: if ( !qpcheck((Obj)COEF(dc)) )
993: return 0;
994: return 1;
995: default:
996: return 0;
997: }
998: }
999:
1000: /* check if p is univariate and all coeffs are INT or LM */
1001:
1002: int uzpcheck(p)
1003: Obj p;
1004: {
1005: DCP dc;
1006: P c;
1007:
1008: if ( !p )
1009: return 1;
1010: else
1011: switch ( OID(p) ) {
1012: case O_N:
1013: return (RATN((Num)p)&&INT(p))||(NID((Num)p)==N_LM);
1014: case O_P:
1015: for ( dc = DC((P)p); dc; dc = NEXT(dc) ) {
1016: c = COEF(dc);
1017: if ( !NUM(c) || !uzpcheck(c) )
1018: return 0;
1019: }
1020: return 1;
1021: default:
1022: return 0;
1023: }
1024: }
1025:
1026: int p_mag(p)
1027: P p;
1028: {
1029: int s;
1030: DCP dc;
1031:
1032: if ( !p )
1033: return 0;
1034: else if ( OID(p) == O_N )
1035: return n_bits(NM((Q)p))+(INT((Q)p)?0:n_bits(DN((Q)p)));
1036: else {
1037: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) )
1038: s += p_mag(COEF(dc));
1039: return s;
1040: }
1041: }
1.2 noro 1042:
1043: int maxblenp(p)
1044: P p;
1045: {
1046: int s,t;
1047: DCP dc;
1048:
1049: if ( !p )
1050: return 0;
1051: else if ( OID(p) == O_N )
1052: return n_bits(NM((Q)p))+(INT((Q)p)?0:n_bits(DN((Q)p)));
1053: else {
1054: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
1055: t = p_mag(COEF(dc));
1056: s = MAX(t,s);
1057: }
1058: return s;
1059: }
1060: }