![[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.9 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/P.c,v 1.8 2003/12/23 10:39:57 ohara 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: }
1.9 ! noro 575: }
! 576:
! 577: int equalp(vl,p1,p2)
! 578: VL vl;
! 579: P p1,p2;
! 580: {
! 581: DCP dc1,dc2;
! 582: V v1,v2;
! 583:
! 584: if ( !p1 ) {
! 585: if ( !p2 ) return 1;
! 586: else return 0;
! 587: }
! 588: /* p1 != 0 */
! 589: if ( !p2 ) return 0;
! 590:
! 591: /* p1 != 0, p2 != 0 */
! 592: if ( NUM(p1) ) {
! 593: if ( !NUM(p2) ) return 0;
! 594: /* p1 and p2 are numbers */
! 595: if ( NID((Num)p1) != NID((Num)p2) ) return 0;
! 596: if ( !(*cmpnumt[NID(p1),NID(p1)])(p1,p2) ) return 1;
! 597: return 0;
! 598: }
! 599: if ( VR(p1) != VR(p2) ) return 0;
! 600: for ( dc1 = DC(p1), dc2 = DC(p2);
! 601: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) ) {
! 602: if ( cmpq(DEG(dc1),DEG(dc2)) ) return 0;
! 603: else if ( !equalp(vl,COEF(dc1),COEF(dc2)) ) return 0;
! 604: }
! 605: if ( dc1 || dc2 ) return 0;
! 606: else return 1;
1.1 noro 607: }
608:
609: void csump(vl,p,c)
610: VL vl;
611: P p;
612: Q *c;
613: {
614: DCP dc;
615: Q s,s1,s2;
616:
617: if ( !p || NUM(p) )
618: *c = (Q)p;
619: else {
620: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
621: csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
622: }
623: *c = s;
624: }
625: }
626:
627: void degp(v,p,d)
628: V v;
629: P p;
630: Q *d;
631: {
632: DCP dc;
633: Q m,m1;
634:
635: if ( !p || NUM(p) )
636: *d = 0;
637: else if ( v == VR(p) )
638: *d = DEG(DC(p));
639: else {
640: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
641: degp(v,COEF(dc),&m1);
642: if ( cmpq(m,m1) < 0 )
643: m = m1;
644: }
645: *d = m;
1.6 noro 646: }
647: }
648:
649: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr);
650: void mulpq_trunc(P p,Q q,VN vn,P *pr);
651: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr);
652:
653: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr)
654: {
655: register DCP dc,dct,dcr,dcr0;
656: V v1,v2;
657: P t,s,u;
658: int n1,n2,i,d,d1;
659:
660: if ( !p1 || !p2 ) *pr = 0;
661: else if ( NUM(p1) )
662: mulpq_trunc(p2,(Q)p1,vn,pr);
663: else if ( NUM(p2) )
664: mulpq_trunc(p1,(Q)p2,vn,pr);
665: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
666: for ( ; vn->v && vn->v != v1; vn++ )
667: if ( vn->n ) {
668: /* p1,p2 do not contain vn->v */
669: *pr = 0;
670: return;
671: }
672: if ( !vn->v )
673: error("mulp_trunc : invalid vn");
674: d = vn->n;
675: for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {
676: for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {
677: d1 = QTOS(DEG(dct))+QTOS(DEG(dc));
678: if ( d1 >= d ) {
679: mulp_trunc(vl,COEF(dct),COEF(dc),vn+1,&t);
680: if ( t ) {
681: NEXTDC(dcr0,dcr);
682: STOQ(d1,DEG(dcr));
683: COEF(dcr) = t;
684: }
685: }
686: }
687: if ( dcr0 ) {
688: NEXT(dcr) = 0; MKP(v1,dcr0,t);
689: addp(vl,s,t,&u); s = u; t = u = 0;
690: }
691: }
692: *pr = s;
693: } else {
694: while ( v1 != VR(vl) && v2 != VR(vl) )
695: vl = NEXT(vl);
696: if ( v1 == VR(vl) )
697: mulpc_trunc(vl,p1,p2,vn,pr);
698: else
699: mulpc_trunc(vl,p2,p1,vn,pr);
700: }
701: }
702:
703: void mulpq_trunc(P p,Q q,VN vn,P *pr)
704: {
705: DCP dc,dcr,dcr0;
706: P t;
707: int i,d;
708: V v;
709:
710: if (!p || !q)
711: *pr = 0;
712: else if ( NUM(p) ) {
713: for ( ; vn->v; vn++ )
714: if ( vn->n ) {
715: *pr = 0;
716: return;
717: }
718: MULNUM(p,q,pr);
719: } else {
720: v = VR(p);
721: for ( ; vn->v && vn->v != v; vn++ ) {
722: if ( vn->n ) {
723: /* p does not contain vn->v */
724: *pr = 0;
725: return;
726: }
727: }
728: if ( !vn->v )
729: error("mulpq_trunc : invalid vn");
730: d = vn->n;
731: for ( dcr0 = 0, dc = DC(p); dc && QTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
732: mulpq_trunc(COEF(dc),q,vn+1,&t);
733: if ( t ) {
734: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
735: }
736: }
737: if ( dcr0 ) {
738: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
739: } else
740: *pr = 0;
741: }
742: }
743:
744: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr)
745: {
746: DCP dc,dcr,dcr0;
747: P t;
748: V v;
749: int i,d;
750:
751: if ( NUM(c) )
752: mulpq_trunc(p,(Q)c,vn,pr);
753: else {
754: v = VR(p);
755: for ( ; vn->v && vn->v != v; vn++ )
756: if ( vn->n ) {
757: /* p,c do not contain vn->v */
758: *pr = 0;
759: return;
760: }
761: if ( !vn->v )
762: error("mulpc_trunc : invalid vn");
763: d = vn->n;
764: for ( dcr0 = 0, dc = DC(p); dc && QTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
765: mulp_trunc(vl,COEF(dc),c,vn+1,&t);
766: if ( t ) {
767: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
768: }
769: }
770: if ( dcr0 ) {
771: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
772: } else
773: *pr = 0;
1.7 noro 774: }
775: }
776:
777: void quop_trunc(VL vl,P p1,P p2,VN vn,P *pr)
778: {
779: DCP dc,dcq0,dcq;
780: P t,s,m,lc2,qt;
781: V v1,v2;
782: Q d2;
783: VN vn1;
784:
785: if ( !p1 )
786: *pr = 0;
787: else if ( NUM(p2) )
788: divsp(vl,p1,p2,pr);
789: else if ( (v1 = VR(p1)) != (v2 = VR(p2)) ) {
790: for ( dcq0 = 0, dc = DC(p1); dc; dc = NEXT(dc) ) {
791: NEXTDC(dcq0,dcq);
792: DEG(dcq) = DEG(dc);
793: quop_trunc(vl,COEF(dc),p2,vn,&COEF(dcq));
794: }
795: NEXT(dcq) = 0;
796: MKP(v1,dcq0,*pr);
797: } else {
798: d2 = DEG(DC(p2));
799: lc2 = COEF(DC(p2));
800: t = p1;
801: dcq0 = 0;
802: /* vn1 = degree list of LC(p2) */
803: for ( vn1 = vn; vn1->v != v1; vn1++ );
804: vn1++;
805: while ( t ) {
806: dc = DC(t);
807: NEXTDC(dcq0,dcq);
808: subq(DEG(dc),d2,&DEG(dcq));
809: quop_trunc(vl,COEF(dc),lc2,vn1,&COEF(dcq));
810: NEXT(dcq) = 0;
811: MKP(v1,dcq,qt);
812: mulp_trunc(vl,p2,qt,vn,&m);
813: subp(vl,t,m,&s); t = s;
814: }
815: NEXT(dcq) = 0;
816: MKP(v1,dcq0,*pr);
1.1 noro 817: }
818: }
819: #endif