![[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.6 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/PU.c,v 1.5 2001/03/29 09:49:57 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; ) {
1.6 ! noro 503: mod = get_lprime(index++);
1.1 noro 504: ptomp(mod,LC(q1),&tg);
505: if ( !tg )
506: continue;
507: ptomp(mod,LC(q2),&tg);
508: if ( !tg )
509: continue;
510: ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
511: srchmp(nvl,mod,v,tg1,tg2,&resg);
512: chnremp(nvl,mod,c,q,resg,&c1); c = c1;
513: STOQ(mod,t); mulq(q,t,&s); q = s;
514: }
515: *pr = c;
516: }
517:
518: void res_ch_det(vl,v,p1,p2,pr)
519: VL vl;
520: V v;
521: P p1,p2,*pr;
522: {
523: P q1,q2,c,c1;
524: P tg,tg1,tg2,resg;
525: int index,mod;
526: Q a,b,f,q,t,s,u,n,m;
527: VL nvl;
528:
529: reordvar(vl,v,&nvl);
530: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
531:
532: if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
533: *pr = 0;
534: return;
535: }
536: if ( VR(q1) != v ) {
537: pwrp(vl,q1,DEG(DC(q2)),pr);
538: return;
539: }
540: if ( VR(q2) != v ) {
541: pwrp(vl,q2,DEG(DC(q1)),pr);
542: return;
543: }
544: norm1c(q1,&a); norm1c(q2,&b);
545: n = DEG(DC(q1)); m = DEG(DC(q2));
546: pwrq(a,m,&t); pwrq(b,n,&s); mulq(t,s,&u);
547: factorial(QTOS(n)+QTOS(m),&t);
548: mulq(u,t,&s); addq(s,s,&f);
549: for ( index = 0, q = ONE, c = 0; cmpq(f,q) >= 0; ) {
1.6 ! noro 550: mod = get_lprime(index++);
1.1 noro 551: ptomp(mod,LC(q1),&tg);
552: if ( !tg )
553: continue;
554: ptomp(mod,LC(q2),&tg);
555: if ( !tg )
556: continue;
557: ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
558: res_detmp(nvl,mod,v,tg1,tg2,&resg);
559: chnremp(nvl,mod,c,q,resg,&c1); c = c1;
560: STOQ(mod,t); mulq(q,t,&s); q = s;
561: }
562: *pr = c;
563: }
564:
565: void res_detmp(vl,mod,v,p1,p2,dp)
566: VL vl;
567: int mod;
568: V v;
569: P p1,p2;
570: P *dp;
571: {
572: int n1,n2,n,sgn;
573: int i,j,k;
574: P mjj,mij,t,s,u,d;
575: P **mat;
576: P *mi,*mj;
577: DCP dc;
578:
579: n1 = UDEG(p1); n2 = UDEG(p2); n = n1+n2;
580: mat = (P **)almat_pointer(n,n);
581: for ( dc = DC(p1); dc; dc = NEXT(dc) )
582: mat[0][n1-QTOS(DEG(dc))] = COEF(dc);
583: for ( i = 1; i < n2; i++ )
584: for ( j = 0; j <= n1; j++ )
585: mat[i][i+j] = mat[0][j];
586: for ( dc = DC(p2); dc; dc = NEXT(dc) )
587: mat[n2][n2-QTOS(DEG(dc))] = COEF(dc);
588: for ( i = 1; i < n1; i++ )
589: for ( j = 0; j <= n2; j++ )
590: mat[i+n2][i+j] = mat[n2][j];
591: for ( j = 0, d = (P)ONEM, sgn = 1; j < n; j++ ) {
592: for ( i = j; (i < n) && !mat[i][j]; i++ );
593: if ( i == n ) {
594: *dp = 0; return;
595: }
596: for ( k = i; k < n; k++ )
597: if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
598: i = k;
599: if ( j != i ) {
600: mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
601: }
602: for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
603: for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
604: mulmp(vl,mod,mi[k],mjj,&t); mulmp(vl,mod,mj[k],mij,&s);
605: submp(vl,mod,t,s,&u); divsmp(vl,mod,u,d,&mi[k]);
606: }
607: d = mjj;
608: }
609: if ( sgn > 0 )
610: *dp = d;
611: else
612: chsgnmp(mod,d,dp);
613: }
614:
615: #if 0
616: showmat(vl,mat,n)
617: VL vl;
618: P **mat;
619: int n;
620: {
621: int i,j;
622: P t;
623:
624: for ( i = 0; i < n; i++ ) {
625: for ( j = 0; j < n; j++ ) {
626: mptop(mat[i][j],&t); printp(vl,t); fprintf(out," ");
627: }
628: fprintf(out,"\n");
629: }
630: fflush(out);
631: }
632:
633: showmp(vl,p)
634: VL vl;
635: P p;
636: {
637: P t;
638:
639: mptop(p,&t); printp(vl,t); fprintf(out,"\n");
640: }
641: #endif
642:
643: void premp(vl,p1,p2,pr)
644: VL vl;
645: P p1,p2,*pr;
646: {
647: P m,m1,m2;
648: P *pw;
649: DCP dc;
650: V v1,v2;
651: register int i,j;
652: int n1,n2,d;
653:
654: if ( NUM(p1) )
655: if ( NUM(p2) )
656: *pr = 0;
657: else
658: *pr = p1;
659: else if ( NUM(p2) )
660: *pr = 0;
661: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
662: n1 = deg(v1,p1); n2 = deg(v1,p2);
663: pw = (P *)ALLOCA((n1+1)*sizeof(P));
664: bzero((char *)pw,(int)((n1+1)*sizeof(P)));
665:
666: for ( dc = DC(p1); dc; dc = NEXT(dc) )
667: pw[QTOS(DEG(dc))] = COEF(dc);
668:
669: for ( i = n1; i >= n2; i-- ) {
670: if ( pw[i] ) {
671: m = pw[i];
672: for ( j = i; j >= 0; j-- ) {
673: mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
674: }
675:
676: for ( dc = DC(p2), d = i - n2; dc; dc = NEXT(dc) ) {
677: mulp(vl,COEF(dc),m,&m1);
678: subp(vl,pw[QTOS(DEG(dc))+d],m1,&m2);
679: pw[QTOS(DEG(dc))+d] = m2;
680: }
681: } else
682: for ( j = i; j >= 0; j-- )
683: if ( pw[j] ) {
684: mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
685: }
686: }
687: plisttop(pw,v1,n2-1,pr);
688: } else {
689: while ( v1 != vl->v && v2 != vl->v )
690: vl = NEXT(vl);
691: if ( v1 == vl->v )
692: *pr = 0;
693: else
694: *pr = p1;
695: }
696: }
697:
698: void ptozp0(p,pr)
699: P p;
700: P *pr;
701: {
702: Q c;
703:
1.5 noro 704: if ( qpcheck((Obj)p) )
705: ptozp(p,1,&c,pr);
706: else
707: *pr = p;
1.1 noro 708: }
709:
710: void mindegp(vl,p,mvlp,pr)
711: VL vl,*mvlp;
712: P p,*pr;
713: {
714: P t;
715: VL nvl,tvl,avl;
716: V v;
717: int n,d;
718:
719: clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
720: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
721: n = getdeg(avl->v,p);
722: if ( n < d ) {
723: v = avl->v; d = n;
724: }
725: }
726: if ( v != nvl->v ) {
727: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
728: *pr = t; *mvlp = tvl;
729: } else {
730: *pr = p; *mvlp = nvl;
731: }
732: }
733:
734: void maxdegp(vl,p,mvlp,pr)
735: VL vl,*mvlp;
736: P p,*pr;
737: {
738: P t;
739: VL nvl,tvl,avl;
740: V v;
741: int n,d;
742:
743: clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
744: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
745: n = getdeg(avl->v,p);
746: if ( n > d ) {
747: v = avl->v; d = n;
748: }
749: }
750: if ( v != nvl->v ) {
751: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
752: *pr = t; *mvlp = tvl;
753: } else {
754: *pr = p; *mvlp = nvl;
755: }
756: }
757:
758: void min_common_vars_in_coefp(vl,p,mvlp,pr)
759: VL vl,*mvlp;
760: P p,*pr;
761: {
762: P u,p0;
763: VL tvl,cvl,svl,uvl,avl,vl0;
764: int n,n0,d,d0;
765: DCP dc;
766:
767: clctv(vl,p,&tvl); vl = tvl;
768: for ( tvl = vl, n0 = 0; tvl; tvl = NEXT(tvl), n0++ );
769: for ( avl = vl; avl; avl = NEXT(avl) ) {
770: if ( VR(p) != avl->v ) {
771: reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
772: } else {
773: uvl = vl; u = p;
774: }
775: for ( cvl = NEXT(uvl), dc = DC(u); dc && cvl; dc = NEXT(dc) ) {
776: clctv(uvl,COEF(dc),&tvl); intersectv(cvl,tvl,&svl); cvl = svl;
777: }
778: if ( !cvl ) {
779: *mvlp = uvl; *pr = u; return;
780: } else {
781: for ( tvl = cvl, n = 0; tvl; tvl = NEXT(tvl), n++ );
782: if ( n < n0 ) {
783: vl0 = uvl; p0 = u; n0 = n; d0 = getdeg(uvl->v,u);
784: } else if ( (n == n0) && ((d = getdeg(uvl->v,u)) < d0) ) {
785: vl0 = uvl; p0 = u; n0 = n; d0 = d;
786: }
787: }
788: }
789: *pr = p0; *mvlp = vl0;
790: }
791:
792: void minlcdegp(vl,p,mvlp,pr)
793: VL vl,*mvlp;
794: P p,*pr;
795: {
796: P u,p0;
797: VL tvl,uvl,avl,vl0;
798: int d,d0;
799:
800: clctv(vl,p,&tvl); vl = tvl;
801: vl0 = vl; p0 = p; d0 = homdeg(COEF(DC(p)));
802: for ( avl = NEXT(vl); avl; avl = NEXT(avl) ) {
803: reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
804: d = homdeg(COEF(DC(u)));
805: if ( d < d0 ) {
806: vl0 = uvl; p0 = u; d0 = d;
807: }
808: }
809: *pr = p0; *mvlp = vl0;
810: }
811:
812: void sort_by_deg(n,p,pr)
813: int n;
814: P *p,*pr;
815: {
816: int j,k,d,k0;
817: V v;
818:
819: for ( j = 0; j < n; j++ ) {
820: for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
821: k < n; k++ )
822: if ( deg(v,p[k]) < d ) {
823: k0 = k;
824: d = deg(v,p[k]);
825: }
826: pr[j] = p[k0];
827: if ( j != k0 )
828: p[k0] = p[j];
829: }
830: }
831:
832: void sort_by_deg_rev(n,p,pr)
833: int n;
834: P *p,*pr;
835: {
836: int j,k,d,k0;
837: V v;
838:
839: for ( j = 0; j < n; j++ ) {
840: for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
841: k < n; k++ )
842: if ( deg(v,p[k]) > d ) {
843: k0 = k;
844: d = deg(v,p[k]);
845: }
846: pr[j] = p[k0];
847: if ( j != k0 )
848: p[k0] = p[j];
849: }
850: }
851:
852:
853: void getmindeg(v,p,dp)
854: V v;
855: P p;
856: Q *dp;
857: {
858: Q dt,d;
859: DCP dc;
860:
861: if ( !p || NUM(p) )
862: *dp = 0;
863: else if ( v == VR(p) ) {
864: for ( dc = DC(p); NEXT(dc); dc = NEXT(dc) );
865: *dp = DEG(dc);
866: } else {
867: dc = DC(p);
868: getmindeg(v,COEF(dc),&d);
869: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
870: getmindeg(v,COEF(dc),&dt);
871: if ( cmpq(dt,d) < 0 )
872: d = dt;
873: }
874: *dp = d;
875: }
876: }
877:
878: void minchdegp(vl,p,mvlp,pr)
879: VL vl,*mvlp;
880: P p,*pr;
881: {
882: P t;
883: VL tvl,nvl,avl;
884: int n,d,z;
885: V v;
886:
887: if ( NUM(p) ) {
888: *mvlp = vl;
889: *pr = p;
890: return;
891: }
892: clctv(vl,p,&nvl);
893: v = nvl->v;
894: d = getchomdeg(v,p) + getlchomdeg(v,p,&z);
895: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
896: n = getchomdeg(avl->v,p) + getlchomdeg(avl->v,p,&z);
897: if ( n < d ) {
898: v = avl->v; d = n;
899: }
900: }
901: if ( v != nvl->v ) {
902: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
903: *pr = t; *mvlp = tvl;
904: } else {
905: *pr = p; *mvlp = nvl;
906: }
907: }
908:
909: int getchomdeg(v,p)
910: V v;
911: P p;
912: {
913: int m,m1;
914: DCP dc;
915:
916: if ( !p || NUM(p) )
917: return ( 0 );
918: else if ( VR(p) == v )
919: return ( dbound(v,p) );
920: else {
921: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
922: m1 = getchomdeg(v,COEF(dc))+QTOS(DEG(dc));
923: m = MAX(m,m1);
924: }
925: return ( m );
926: }
927: }
928:
929: int getlchomdeg(v,p,d)
930: V v;
931: P p;
932: int *d;
933: {
934: int m0,m1,d0,d1;
935: DCP dc;
936:
937: if ( !p || NUM(p) ) {
938: *d = 0;
939: return ( 0 );
940: } else if ( VR(p) == v ) {
941: *d = QTOS(DEG(DC(p)));
942: return ( homdeg(LC(p)) );
943: } else {
944: for ( dc = DC(p), m0 = 0, d0 = 0; dc; dc = NEXT(dc) ) {
945: m1 = getlchomdeg(v,COEF(dc),&d1)+QTOS(DEG(dc));
946: if ( d1 > d0 ) {
947: m0 = m1;
948: d0 = d1;
949: } else if ( d1 == d0 )
950: m0 = MAX(m0,m1);
951: }
952: *d = d0;
953: return ( m0 );
954: }
955: }
956:
957: int nmonop(p)
958: P p;
959: {
960: int s;
961: DCP dc;
962:
963: if ( !p )
964: return 0;
965: else
966: switch ( OID(p) ) {
967: case O_N:
968: return 1; break;
969: case O_P:
970: for ( dc = DC((P)p), s = 0; dc; dc = NEXT(dc) )
971: s += nmonop(COEF(dc));
972: return s; break;
973: default:
974: return 0;
975: }
976: }
977:
978: int qpcheck(p)
979: Obj p;
980: {
981: DCP dc;
982:
983: if ( !p )
984: return 1;
985: else
986: switch ( OID(p) ) {
987: case O_N:
988: return RATN((Num)p)?1:0;
989: case O_P:
990: for ( dc = DC((P)p); dc; dc = NEXT(dc) )
991: if ( !qpcheck((Obj)COEF(dc)) )
992: return 0;
993: return 1;
994: default:
995: return 0;
996: }
997: }
998:
999: /* check if p is univariate and all coeffs are INT or LM */
1000:
1001: int uzpcheck(p)
1002: Obj p;
1003: {
1004: DCP dc;
1005: P c;
1006:
1007: if ( !p )
1008: return 1;
1009: else
1010: switch ( OID(p) ) {
1011: case O_N:
1012: return (RATN((Num)p)&&INT(p))||(NID((Num)p)==N_LM);
1013: case O_P:
1014: for ( dc = DC((P)p); dc; dc = NEXT(dc) ) {
1015: c = COEF(dc);
1016: if ( !NUM(c) || !uzpcheck(c) )
1017: return 0;
1018: }
1019: return 1;
1020: default:
1021: return 0;
1022: }
1023: }
1024:
1025: int p_mag(p)
1026: P p;
1027: {
1028: int s;
1029: DCP dc;
1030:
1031: if ( !p )
1032: return 0;
1033: else if ( OID(p) == O_N )
1034: return n_bits(NM((Q)p))+(INT((Q)p)?0:n_bits(DN((Q)p)));
1035: else {
1036: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) )
1037: s += p_mag(COEF(dc));
1038: return s;
1039: }
1040: }
1.2 noro 1041:
1042: int maxblenp(p)
1043: P p;
1044: {
1045: int s,t;
1046: DCP dc;
1047:
1048: if ( !p )
1049: return 0;
1050: else if ( OID(p) == O_N )
1051: return n_bits(NM((Q)p))+(INT((Q)p)?0:n_bits(DN((Q)p)));
1052: else {
1053: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
1054: t = p_mag(COEF(dc));
1055: s = MAX(t,s);
1056: }
1057: return s;
1058: }
1059: }