![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to P.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 P.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.2 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.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 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.8 ! ohara 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/P.c,v 1.7 2003/06/24 09:49:36 noro Exp $
1.2 noro 49: */
1.1 noro 50: #ifndef FBASE
51: #define FBASE
52: #endif
53:
54: #include "b.h"
55: #include "ca.h"
56:
57: void D_ADDP(vl,p1,p2,pr)
58: VL vl;
59: P p1,p2,*pr;
60: {
61: register DCP dc1,dc2,dcr0,dcr;
62: V v1,v2;
63: P t;
64:
65: if ( !p1 )
66: *pr = p2;
67: else if ( !p2 )
68: *pr = p1;
69: else if ( NUM(p1) )
70: if ( NUM(p2) )
71: ADDNUM(p1,p2,pr);
72: else
73: ADDPQ(p2,p1,pr);
74: else if ( NUM(p2) )
75: ADDPQ(p1,p2,pr);
76: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
77: for ( dc1 = DC(p1), dc2 = DC(p2), dcr0 = 0; dc1 && dc2; )
78: switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
79: case 0:
80: ADDP(vl,COEF(dc1),COEF(dc2),&t);
81: if ( t ) {
82: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc1); COEF(dcr) = t;
83: }
84: dc1 = NEXT(dc1); dc2 = NEXT(dc2); break;
85: case 1:
86: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc1); COEF(dcr) = COEF(dc1);
87: dc1 = NEXT(dc1); break;
88: case -1:
89: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc2); COEF(dcr) = COEF(dc2);
90: dc2 = NEXT(dc2); break;
91: }
92: if ( !dcr0 )
93: if ( dc1 )
94: dcr0 = dc1;
95: else if ( dc2 )
96: dcr0 = dc2;
97: else {
98: *pr = 0;
99: return;
100: }
101: else
102: if ( dc1 )
103: NEXT(dcr) = dc1;
104: else if ( dc2 )
105: NEXT(dcr) = dc2;
106: else
107: NEXT(dcr) = 0;
108: MKP(v1,dcr0,*pr);
109: } else {
110: while ( v1 != VR(vl) && v2 != VR(vl) )
111: vl = NEXT(vl);
112: if ( v1 == VR(vl) )
113: ADDPTOC(vl,p1,p2,pr);
114: else
115: ADDPTOC(vl,p2,p1,pr);
116: }
117: }
118:
119: void D_ADDPQ(p,q,pr)
120: P p,q,*pr;
121: {
122: DCP dc,dcr,dcr0;
123: P t;
124:
125: if ( NUM(p) )
126: ADDNUM(p,q,pr);
127: else {
128: for ( dcr0 = 0, dc = DC(p); dc && DEG(dc); dc = NEXT(dc) ) {
129: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = COEF(dc);
130: }
131: if ( !dc ) {
132: NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = q;
133: } else {
134: ADDPQ(COEF(dc),q,&t);
135: if ( t ) {
136: NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = t;
137: }
138: }
139: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
140: }
141: }
142:
143: void D_ADDPTOC(vl,p,c,pr)
144: VL vl;
145: P p,c,*pr;
146: {
147: DCP dc,dcr,dcr0;
148: P t;
149:
150: for ( dcr0 = 0, dc = DC(p); dc && DEG(dc); dc = NEXT(dc) ) {
151: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = COEF(dc);
152: }
153: if ( !dc ) {
154: NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = c;
155: } else {
156: ADDP(vl,COEF(dc),c,&t);
157: if ( t ) {
158: NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = t;
159: }
160: }
161: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
162: }
163:
164: void D_SUBP(vl,p1,p2,pr)
165: VL vl;
166: P p1,p2,*pr;
167: {
168: P t;
169:
170: if ( !p2 )
171: *pr = p1;
172: else {
173: CHSGNP(p2,&t); ADDP(vl,p1,t,pr);
174: }
175: }
176:
177: void D_MULP(vl,p1,p2,pr)
178: VL vl;
179: P p1,p2,*pr;
180: {
181: register DCP dc,dct,dcr,dcr0;
182: V v1,v2;
183: P t,s,u;
184: int n1,n2;
185:
186: if ( !p1 || !p2 ) *pr = 0;
187: else if ( NUM(p1) )
188: MULPQ(p2,p1,pr);
189: else if ( NUM(p2) )
190: MULPQ(p1,p2,pr);
191: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
192: for ( dc = DC(p1), n1 = 0; dc; dc = NEXT(dc), n1++ );
193: for ( dc = DC(p2), n2 = 0; dc; dc = NEXT(dc), n2++ );
194: if ( n1 > n2 )
195: for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {
196: for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {
197: NEXTDC(dcr0,dcr); MULP(vl,COEF(dct),COEF(dc),&COEF(dcr));
198: addq(DEG(dct),DEG(dc),&DEG(dcr));
199: }
200: NEXT(dcr) = 0; MKP(v1,dcr0,t);
201: ADDP(vl,s,t,&u); s = u; t = u = 0;
202: }
203: else
204: for ( dc = DC(p1), s = 0; dc; dc = NEXT(dc) ) {
205: for ( dcr0 = 0, dct = DC(p2); dct; dct = NEXT(dct) ) {
206: NEXTDC(dcr0,dcr); MULP(vl,COEF(dct),COEF(dc),&COEF(dcr));
207: addq(DEG(dct),DEG(dc),&DEG(dcr));
208: }
209: NEXT(dcr) = 0; MKP(v1,dcr0,t);
210: ADDP(vl,s,t,&u); s = u; t = u = 0;
211: }
212: *pr = s;
213: } else {
214: while ( v1 != VR(vl) && v2 != VR(vl) )
215: vl = NEXT(vl);
216: if ( v1 == VR(vl) )
217: MULPC(vl,p1,p2,pr);
218: else
219: MULPC(vl,p2,p1,pr);
220: }
221: }
222:
223: void D_MULPQ(p,q,pr)
224: P p,q,*pr;
225: {
226: DCP dc,dcr,dcr0;
227: P t;
228:
229: if (!p || !q)
230: *pr = 0;
231: else if ( Uniq(q) )
232: *pr = p;
233: else if ( NUM(p) )
234: MULNUM(p,q,pr);
235: else {
236: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
237: MULPQ(COEF(dc),q,&t);
238: if ( t ) {
239: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
240: }
241: }
242: if ( dcr0 ) {
243: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
244: } else
245: *pr = 0;
246: }
247: }
248:
249: void D_MULPC(vl,p,c,pr)
250: VL vl;
251: P p,c,*pr;
252: {
253: DCP dc,dcr,dcr0;
254: P t;
255:
256: if ( NUM(c) )
257: MULPQ(p,c,pr);
258: else {
259: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
260: MULP(vl,COEF(dc),c,&t);
261: if ( t ) {
262: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
263: }
264: }
265: if ( dcr0 ) {
266: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
267: } else
268: *pr = 0;
269: }
270: }
271:
272: void D_PWRP(vl,p,q,pr)
273: VL vl;
274: P p,*pr;
275: Q q;
276: {
277: DCP dc,dcr;
278: int n,i;
279: P *x,*y;
280: P t,s,u;
281: DCP dct;
282: P *pt;
283:
284: if ( !q ) {
285: *pr = (P)One;
286: } else if ( !p )
287: *pr = 0;
288: else if ( UNIQ(q) )
289: *pr = p;
290: else if ( NUM(p) )
291: PWRNUM(p,q,pr);
292: else {
293: dc = DC(p);
294: if ( !NEXT(dc) ) {
295: NEWDC(dcr);
296: PWRP(vl,COEF(dc),q,&COEF(dcr)); mulq(DEG(dc),q,&DEG(dcr));
297: NEXT(dcr) = 0; MKP(VR(p),dcr,*pr);
298: } else if ( !INT(q) ) {
299: error("pwrp: can't calculate fractional power."); *pr = 0;
300: } else if ( PL(NM(q)) == 1 ) {
301: n = QTOS(q); x = (P *)ALLOCA((n+1)*sizeof(pointer));
302: NEWDC(dct); DEG(dct) = DEG(dc); COEF(dct) = COEF(dc);
303: NEXT(dct) = 0; MKP(VR(p),dct,t);
304: for ( i = 0, u = (P)One; i < n; i++ ) {
305: x[i] = u; MULP(vl,u,t,&s); u = s;
306: }
307: x[n] = u; y = (P *)ALLOCA((n+1)*sizeof(pointer));
308: if ( DEG(NEXT(dc)) ) {
309: dct = NEXT(dc); MKP(VR(p),dct,t);
310: } else
311: t = COEF(NEXT(dc));
312: for ( i = 0, u = (P)One; i < n; i++ ) {
313: y[i] = u; MULP(vl,u,t,&s); u = s;
314: }
315: y[n] = u;
316: pt = (P *)ALLOCA((n+1)*sizeof(pointer)); MKBC(n,pt);
317: for ( i = 0, u = 0; i <= n; i++ ) {
318: MULP(vl,x[i],y[n-i],&t); MULP(vl,t,pt[i],&s);
319: ADDP(vl,u,s,&t); u = t;
320: }
321: *pr = u;
322: } else {
323: error("exponent too big");
324: *pr = 0;
325: }
326: }
327: }
328:
329: void D_CHSGNP(p,pr)
330: P p,*pr;
331: {
332: register DCP dc,dcr,dcr0;
333: P t;
334:
335: if ( !p )
336: *pr = NULL;
337: else if ( NUM(p) ) {
1.5 noro 338: #if defined(_PA_RISC1_1) || defined(__alpha) || defined(mips) || defined(_IBMR2)
1.1 noro 339: #ifdef FBASE
340: chsgnnum((Num)p,(Num *)pr);
341: #else
342: MQ mq;
343:
344: NEWMQ(mq); CONT(mq)=mod-CONT((MQ)p); *pr = (P)mq;
345: #endif
346: #else
347: CHSGNNUM(p,t); *pr = t;
348: #endif
349: } else {
350: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
351: NEXTDC(dcr0,dcr); CHSGNP(COEF(dc),&COEF(dcr)); DEG(dcr) = DEG(dc);
352: }
353: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
354: }
355: }
356:
357:
358: #ifdef FBASE
359: void ptozp(p,sgn,c,pr)
360: P p;
361: int sgn;
362: Q *c;
363: P *pr;
364: {
365: N nm,dn;
366:
367: if ( !p ) {
368: *c = 0; *pr = 0;
369: } else {
370: lgp(p,&nm,&dn);
371: if ( UNIN(dn) )
372: NTOQ(nm,sgn,*c);
373: else
374: NDTOQ(nm,dn,sgn,*c);
375: divcp(p,*c,pr);
376: }
377: }
378:
379: void lgp(p,g,l)
380: P p;
381: N *g,*l;
382: {
383: N g1,g2,l1,l2,l3,l4,tmp;
384: DCP dc;
385:
386: if ( NUM(p) ) {
387: *g = NM((Q)p);
388: if ( INT((Q)p) )
389: *l = ONEN;
390: else
391: *l = DN((Q)p);
392: } else {
393: dc = DC(p); lgp(COEF(dc),g,l);
394: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
395: lgp(COEF(dc),&g1,&l1); gcdn(*g,g1,&g2); *g = g2;
396: gcdn(*l,l1,&l2); kmuln(*l,l1,&l3); divn(l3,l2,&l4,&tmp); *l = l4;
397: }
398: }
399: }
400:
401: void divcp(f,q,rp)
402: P f;
403: Q q;
404: P *rp;
405: {
406: DCP dc,dcr,dcr0;
407:
408: if ( !f )
409: *rp = 0;
410: else if ( NUM(f) )
411: DIVNUM(f,q,rp);
412: else {
413: for ( dc = DC(f), dcr0 = 0; dc; dc = NEXT(dc) ) {
414: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
415: divcp(COEF(dc),q,&COEF(dcr));
416: }
417: NEXT(dcr) = 0; MKP(VR(f),dcr0,*rp);
418: }
419: }
420:
421: void diffp(vl,p,v,r)
422: VL vl;
423: P p;
424: V v;
425: P *r;
426: {
427: P t;
428: DCP dc,dcr,dcr0;
429:
430: if ( !p || NUM(p) )
431: *r = 0;
432: else {
433: if ( v == VR(p) ) {
434: for ( dc = DC(p), dcr0 = 0;
435: dc && DEG(dc); dc = NEXT(dc) ) {
436: MULPQ(COEF(dc),(P)DEG(dc),&t);
437: if ( t ) {
438: NEXTDC(dcr0,dcr); SUBQ(DEG(dc),ONE,&DEG(dcr));
439: COEF(dcr) = t;
440: }
441: }
442: if ( !dcr0 )
443: *r = 0;
444: else {
445: NEXT(dcr) = 0; MKP(v,dcr0,*r);
446: }
447: } else {
448: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
449: diffp(vl,COEF(dc),v,&t);
1.8 ! ohara 450: if ( t ) {
! 451: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
! 452: }
! 453: }
! 454: if ( !dcr0 )
! 455: *r = 0;
! 456: else {
! 457: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
! 458: }
! 459: }
! 460: }
! 461: }
! 462:
! 463: /* Euler derivation */
! 464: void ediffp(vl,p,v,r)
! 465: VL vl;
! 466: P p;
! 467: V v;
! 468: P *r;
! 469: {
! 470: P t;
! 471: DCP dc,dcr,dcr0;
! 472:
! 473: if ( !p || NUM(p) )
! 474: *r = 0;
! 475: else {
! 476: if ( v == VR(p) ) {
! 477: for ( dc = DC(p), dcr0 = 0;
! 478: dc && DEG(dc); dc = NEXT(dc) ) {
! 479: MULPQ(COEF(dc),(P)DEG(dc),&t);
! 480: if ( t ) {
! 481: NEXTDC(dcr0,dcr);
! 482: DEG(dcr) = DEG(dc);
! 483: COEF(dcr) = t;
! 484: }
! 485: }
! 486: if ( !dcr0 )
! 487: *r = 0;
! 488: else {
! 489: NEXT(dcr) = 0; MKP(v,dcr0,*r);
! 490: }
! 491: } else {
! 492: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 493: ediffp(vl,COEF(dc),v,&t);
1.1 noro 494: if ( t ) {
495: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
496: }
497: }
498: if ( !dcr0 )
499: *r = 0;
500: else {
501: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
502: }
503: }
504: }
505: }
506:
507: void coefp(p,d,pr)
508: P p;
509: int d;
510: P *pr;
511: {
512: DCP dc;
513: int sgn;
514: Q dq;
515:
516: if ( NUM(p) )
517: if ( d == 0 )
518: *pr = p;
519: else
520: *pr = 0;
521: else {
522: for ( STOQ(d,dq), dc = DC(p); dc; dc = NEXT(dc) )
523: if ( (sgn = cmpq(DEG(dc),dq)) > 0 )
524: continue;
525: else if ( sgn == 0 ) {
526: *pr = COEF(dc);
527: return;
528: } else {
529: *pr = 0;
530: break;
531: }
532: *pr = 0;
533: }
534: }
535:
536: int compp(vl,p1,p2)
537: VL vl;
538: P p1,p2;
539: {
540: DCP dc1,dc2;
541: V v1,v2;
542:
543: if ( !p1 )
544: return p2 ? -1 : 0;
545: else if ( !p2 )
546: return 1;
547: else if ( NUM(p1) )
548: return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;
549: else if ( NUM(p2) )
550: return 1;
551: if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {
552: for ( dc1 = DC(p1), dc2 = DC(p2);
553: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )
554: switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
555: case 1:
556: return 1; break;
557: case -1:
558: return -1; break;
559: default:
560: switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {
561: case 1:
562: return 1; break;
563: case -1:
564: return -1; break;
565: default:
566: break;
567: }
568: break;
569: }
570: return dc1 ? 1 : (dc2 ? -1 : 0);
571: } else {
572: for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );
573: return v1 == VR(vl) ? 1 : -1;
574: }
575: }
576:
577: void csump(vl,p,c)
578: VL vl;
579: P p;
580: Q *c;
581: {
582: DCP dc;
583: Q s,s1,s2;
584:
585: if ( !p || NUM(p) )
586: *c = (Q)p;
587: else {
588: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
589: csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
590: }
591: *c = s;
592: }
593: }
594:
595: void degp(v,p,d)
596: V v;
597: P p;
598: Q *d;
599: {
600: DCP dc;
601: Q m,m1;
602:
603: if ( !p || NUM(p) )
604: *d = 0;
605: else if ( v == VR(p) )
606: *d = DEG(DC(p));
607: else {
608: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
609: degp(v,COEF(dc),&m1);
610: if ( cmpq(m,m1) < 0 )
611: m = m1;
612: }
613: *d = m;
1.6 noro 614: }
615: }
616:
617: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr);
618: void mulpq_trunc(P p,Q q,VN vn,P *pr);
619: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr);
620:
621: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr)
622: {
623: register DCP dc,dct,dcr,dcr0;
624: V v1,v2;
625: P t,s,u;
626: int n1,n2,i,d,d1;
627:
628: if ( !p1 || !p2 ) *pr = 0;
629: else if ( NUM(p1) )
630: mulpq_trunc(p2,(Q)p1,vn,pr);
631: else if ( NUM(p2) )
632: mulpq_trunc(p1,(Q)p2,vn,pr);
633: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
634: for ( ; vn->v && vn->v != v1; vn++ )
635: if ( vn->n ) {
636: /* p1,p2 do not contain vn->v */
637: *pr = 0;
638: return;
639: }
640: if ( !vn->v )
641: error("mulp_trunc : invalid vn");
642: d = vn->n;
643: for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {
644: for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {
645: d1 = QTOS(DEG(dct))+QTOS(DEG(dc));
646: if ( d1 >= d ) {
647: mulp_trunc(vl,COEF(dct),COEF(dc),vn+1,&t);
648: if ( t ) {
649: NEXTDC(dcr0,dcr);
650: STOQ(d1,DEG(dcr));
651: COEF(dcr) = t;
652: }
653: }
654: }
655: if ( dcr0 ) {
656: NEXT(dcr) = 0; MKP(v1,dcr0,t);
657: addp(vl,s,t,&u); s = u; t = u = 0;
658: }
659: }
660: *pr = s;
661: } else {
662: while ( v1 != VR(vl) && v2 != VR(vl) )
663: vl = NEXT(vl);
664: if ( v1 == VR(vl) )
665: mulpc_trunc(vl,p1,p2,vn,pr);
666: else
667: mulpc_trunc(vl,p2,p1,vn,pr);
668: }
669: }
670:
671: void mulpq_trunc(P p,Q q,VN vn,P *pr)
672: {
673: DCP dc,dcr,dcr0;
674: P t;
675: int i,d;
676: V v;
677:
678: if (!p || !q)
679: *pr = 0;
680: else if ( NUM(p) ) {
681: for ( ; vn->v; vn++ )
682: if ( vn->n ) {
683: *pr = 0;
684: return;
685: }
686: MULNUM(p,q,pr);
687: } else {
688: v = VR(p);
689: for ( ; vn->v && vn->v != v; vn++ ) {
690: if ( vn->n ) {
691: /* p does not contain vn->v */
692: *pr = 0;
693: return;
694: }
695: }
696: if ( !vn->v )
697: error("mulpq_trunc : invalid vn");
698: d = vn->n;
699: for ( dcr0 = 0, dc = DC(p); dc && QTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
700: mulpq_trunc(COEF(dc),q,vn+1,&t);
701: if ( t ) {
702: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
703: }
704: }
705: if ( dcr0 ) {
706: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
707: } else
708: *pr = 0;
709: }
710: }
711:
712: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr)
713: {
714: DCP dc,dcr,dcr0;
715: P t;
716: V v;
717: int i,d;
718:
719: if ( NUM(c) )
720: mulpq_trunc(p,(Q)c,vn,pr);
721: else {
722: v = VR(p);
723: for ( ; vn->v && vn->v != v; vn++ )
724: if ( vn->n ) {
725: /* p,c do not contain vn->v */
726: *pr = 0;
727: return;
728: }
729: if ( !vn->v )
730: error("mulpc_trunc : invalid vn");
731: d = vn->n;
732: for ( dcr0 = 0, dc = DC(p); dc && QTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
733: mulp_trunc(vl,COEF(dc),c,vn+1,&t);
734: if ( t ) {
735: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
736: }
737: }
738: if ( dcr0 ) {
739: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
740: } else
741: *pr = 0;
1.7 noro 742: }
743: }
744:
745: void quop_trunc(VL vl,P p1,P p2,VN vn,P *pr)
746: {
747: DCP dc,dcq0,dcq;
748: P t,s,m,lc2,qt;
749: V v1,v2;
750: Q d2;
751: VN vn1;
752:
753: if ( !p1 )
754: *pr = 0;
755: else if ( NUM(p2) )
756: divsp(vl,p1,p2,pr);
757: else if ( (v1 = VR(p1)) != (v2 = VR(p2)) ) {
758: for ( dcq0 = 0, dc = DC(p1); dc; dc = NEXT(dc) ) {
759: NEXTDC(dcq0,dcq);
760: DEG(dcq) = DEG(dc);
761: quop_trunc(vl,COEF(dc),p2,vn,&COEF(dcq));
762: }
763: NEXT(dcq) = 0;
764: MKP(v1,dcq0,*pr);
765: } else {
766: d2 = DEG(DC(p2));
767: lc2 = COEF(DC(p2));
768: t = p1;
769: dcq0 = 0;
770: /* vn1 = degree list of LC(p2) */
771: for ( vn1 = vn; vn1->v != v1; vn1++ );
772: vn1++;
773: while ( t ) {
774: dc = DC(t);
775: NEXTDC(dcq0,dcq);
776: subq(DEG(dc),d2,&DEG(dcq));
777: quop_trunc(vl,COEF(dc),lc2,vn1,&COEF(dcq));
778: NEXT(dcq) = 0;
779: MKP(v1,dcq,qt);
780: mulp_trunc(vl,p2,qt,vn,&m);
781: subp(vl,t,m,&s); t = s;
782: }
783: NEXT(dcq) = 0;
784: MKP(v1,dcq0,*pr);
1.1 noro 785: }
786: }
787: #endif