![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to PM.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 PM.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
! 26: * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification
! 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: *
! 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/PM.c,v 1.1.1.1 1999/12/03 07:39:08 noro Exp $
! 49: */
1.1 noro 50: #ifndef MODULAR
51: #define MODULAR
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:
228: if (!p || !q)
229: *pr = 0;
230: else if ( Uniq(q) )
231: *pr = p;
232: else if ( NUM(p) )
233: MULNUM(p,q,pr);
234: else {
235: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
236: NEXTDC(dcr0,dcr); MULPQ(COEF(dc),q,&COEF(dcr)); DEG(dcr) = DEG(dc);
237: }
238: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
239: }
240: }
241:
242: void D_MULPC(vl,p,c,pr)
243: VL vl;
244: P p,c,*pr;
245: {
246: DCP dc,dcr,dcr0;
247:
248: if ( NUM(c) )
249: MULPQ(p,c,pr);
250: else {
251: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
252: NEXTDC(dcr0,dcr); MULP(vl,COEF(dc),c,&COEF(dcr)); DEG(dcr) = DEG(dc);
253: }
254: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
255: }
256: }
257:
258: void D_PWRP(vl,p,q,pr)
259: VL vl;
260: P p,*pr;
261: Q q;
262: {
263: register DCP dc,dcr;
264: int n,i;
265: P *x,*y;
266: P t,s,u;
267: DCP dct;
268: P *pt;
269:
270: if ( !q ) {
271: *pr = (P)One;
272: } else if ( !p )
273: *pr = 0;
274: else if ( UNIQ(q) )
275: *pr = p;
276: else if ( NUM(p) )
277: PWRNUM(p,q,pr);
278: else {
279: dc = DC(p);
280: if ( !NEXT(dc) ) {
281: NEWDC(dcr);
282: PWRP(vl,COEF(dc),q,&COEF(dcr)); mulq(DEG(dc),q,&DEG(dcr));
283: NEXT(dcr) = 0; MKP(VR(p),dcr,*pr);
284: } else if ( !INT(q) ) {
285: error("pwrp: can't calculate fractional power."); *pr = 0;
286: } else if ( PL(NM(q)) == 1 ) {
287: n = QTOS(q); x = (P *)ALLOCA((n+1)*sizeof(pointer));
288: NEWDC(dct); DEG(dct) = DEG(dc); COEF(dct) = COEF(dc);
289: NEXT(dct) = 0; MKP(VR(p),dct,t);
290: for ( i = 0, u = (P)One; i < n; i++ ) {
291: x[i] = u; MULP(vl,u,t,&s); u = s;
292: }
293: x[n] = u; y = (P *)ALLOCA((n+1)*sizeof(pointer));
294: if ( DEG(NEXT(dc)) ) {
295: dct = NEXT(dc); MKP(VR(p),dct,t);
296: } else
297: t = COEF(NEXT(dc));
298: for ( i = 0, u = (P)One; i < n; i++ ) {
299: y[i] = u; MULP(vl,u,t,&s); u = s;
300: }
301: y[n] = u;
302: pt = (P *)ALLOCA((n+1)*sizeof(pointer)); MKBC(n,pt);
303: for ( i = 0, u = 0; i <= n; i++ ) {
304: MULP(vl,x[i],y[n-i],&t); MULP(vl,t,pt[i],&s);
305: ADDP(vl,u,s,&t); u = t;
306: }
307: *pr = u;
308: } else {
309: error("exponent too big");
310: *pr = 0;
311: }
312: }
313: }
314:
315: void D_CHSGNP(p,pr)
316: P p,*pr;
317: {
318: register DCP dc,dcr,dcr0;
319: P t;
320:
321: if ( !p )
322: *pr = NULL;
323: else if ( NUM(p) ) {
324: #if defined(THINK_C) || defined(_PA_RISC1_1) || defined(__alpha) || defined(mips)
325: #ifdef FBASE
326: chsgnnum((Num)p,(Num *)pr);
327: #else
328: MQ mq;
329:
330: NEWMQ(mq); CONT(mq)=mod-CONT((MQ)p); *pr = (P)mq;
331: #endif
332: #else
333: CHSGNNUM(p,t); *pr = t;
334: #endif
335: } else {
336: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
337: NEXTDC(dcr0,dcr); CHSGNP(COEF(dc),&COEF(dcr)); DEG(dcr) = DEG(dc);
338: }
339: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
340: }
341: }
342:
343:
344: #ifdef FBASE
345: void ptozp(p,sgn,c,pr)
346: P p;
347: int sgn;
348: Q *c;
349: P *pr;
350: {
351: N nm,dn;
352:
353: if ( !p ) {
354: *c = 0; *pr = 0;
355: } else {
356: lgp(p,&nm,&dn);
357: if ( UNIN(dn) )
358: NTOQ(nm,sgn,*c);
359: else
360: NDTOQ(nm,dn,sgn,*c);
361: divcp(p,*c,pr);
362: }
363: }
364:
365: void lgp(p,g,l)
366: P p;
367: N *g,*l;
368: {
369: N g1,g2,l1,l2,l3,l4,tmp;
370: DCP dc;
371:
372: if ( NUM(p) ) {
373: *g = NM((Q)p);
374: if ( INT((Q)p) )
375: *l = ONEN;
376: else
377: *l = DN((Q)p);
378: } else {
379: dc = DC(p); lgp(COEF(dc),g,l);
380: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
381: lgp(COEF(dc),&g1,&l1); gcdn(*g,g1,&g2); *g = g2;
382: gcdn(*l,l1,&l2); kmuln(*l,l1,&l3); divn(l3,l2,&l4,&tmp); *l = l4;
383: }
384: }
385: }
386:
387: void divcp(f,q,rp)
388: P f;
389: Q q;
390: P *rp;
391: {
392: DCP dc,dcr,dcr0;
393:
394: if ( !f )
395: *rp = 0;
396: else if ( NUM(f) )
397: DIVNUM(f,q,rp);
398: else {
399: for ( dc = DC(f), dcr0 = 0; dc; dc = NEXT(dc) ) {
400: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
401: divcp(COEF(dc),q,&COEF(dcr));
402: }
403: NEXT(dcr) = 0; MKP(VR(f),dcr0,*rp);
404: }
405: }
406:
407: void diffp(vl,p,v,r)
408: VL vl;
409: P p;
410: V v;
411: P *r;
412: {
413: P t;
414: DCP dc,dcr,dcr0;
415:
416: if ( !p || NUM(p) )
417: *r = 0;
418: else {
419: if ( v == VR(p) ) {
420: for ( dc = DC(p), dcr0 = 0;
421: dc && DEG(dc); dc = NEXT(dc) ) {
422: NEXTDC(dcr0,dcr); SUBQ(DEG(dc),ONE,&DEG(dcr));
423: MULPQ(COEF(dc),(P)DEG(dc),&COEF(dcr));
424: }
425: NEXT(dcr) = 0; MKP(v,dcr0,*r);
426: } else {
427: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
428: diffp(vl,COEF(dc),v,&t);
429: if ( t ) {
430: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
431: }
432: }
433: if ( !dcr0 )
434: *r = 0;
435: else {
436: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
437: }
438: }
439: }
440: }
441:
442: void coefp(p,d,pr)
443: P p;
444: int d;
445: P *pr;
446: {
447: DCP dc;
448: int sgn;
449: Q dq;
450:
451: if ( NUM(p) )
452: if ( d == 0 )
453: *pr = p;
454: else
455: *pr = 0;
456: else {
457: for ( STOQ(d,dq), dc = DC(p); dc; dc = NEXT(dc) )
458: if ( (sgn = cmpq(DEG(dc),dq)) > 0 )
459: continue;
460: else if ( sgn == 0 ) {
461: *pr = COEF(dc);
462: return;
463: } else {
464: *pr = 0;
465: break;
466: }
467: *pr = 0;
468: }
469: }
470:
471: int compp(vl,p1,p2)
472: VL vl;
473: P p1,p2;
474: {
475: DCP dc1,dc2;
476: V v1,v2;
477:
478: if ( !p1 )
479: return p2 ? -1 : 0;
480: else if ( !p2 )
481: return 1;
482: else if ( NUM(p1) )
483: return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;
484: else if ( NUM(p2) )
485: return 1;
486: if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {
487: for ( dc1 = DC(p1), dc2 = DC(p2);
488: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )
489: switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
490: case 1:
491: return 1; break;
492: case -1:
493: return -1; break;
494: default:
495: switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {
496: case 1:
497: return 1; break;
498: case -1:
499: return -1; break;
500: default:
501: break;
502: }
503: break;
504: }
505: return dc1 ? 1 : (dc2 ? -1 : 0);
506: } else {
507: for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );
508: return v1 == VR(vl) ? 1 : -1;
509: }
510: }
511:
512: void csump(vl,p,c)
513: VL vl;
514: P p;
515: Q *c;
516: {
517: DCP dc;
518: Q s,s1,s2;
519:
520: if ( !p || NUM(p) )
521: *c = (Q)p;
522: else {
523: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
524: csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
525: }
526: *c = s;
527: }
528: }
529:
530: void degp(v,p,d)
531: V v;
532: P p;
533: Q *d;
534: {
535: DCP dc;
536: Q m,m1;
537:
538: if ( !p || NUM(p) )
539: *d = 0;
540: else if ( v == VR(p) )
541: *d = DEG(DC(p));
542: else {
543: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
544: degp(v,COEF(dc),&m1);
545: if ( cmpq(m,m1) < 0 )
546: m = m1;
547: }
548: *d = m;
549: }
550: }
551: #endif