![[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
! 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/P.c,v 1.1.1.1 1999/12/03 07:39:08 noro Exp $
! 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) ) {
338: #if defined(THINK_C) || defined(_PA_RISC1_1) || defined(__alpha) || defined(mips)
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);
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: void coefp(p,d,pr)
464: P p;
465: int d;
466: P *pr;
467: {
468: DCP dc;
469: int sgn;
470: Q dq;
471:
472: if ( NUM(p) )
473: if ( d == 0 )
474: *pr = p;
475: else
476: *pr = 0;
477: else {
478: for ( STOQ(d,dq), dc = DC(p); dc; dc = NEXT(dc) )
479: if ( (sgn = cmpq(DEG(dc),dq)) > 0 )
480: continue;
481: else if ( sgn == 0 ) {
482: *pr = COEF(dc);
483: return;
484: } else {
485: *pr = 0;
486: break;
487: }
488: *pr = 0;
489: }
490: }
491:
492: int compp(vl,p1,p2)
493: VL vl;
494: P p1,p2;
495: {
496: DCP dc1,dc2;
497: V v1,v2;
498:
499: if ( !p1 )
500: return p2 ? -1 : 0;
501: else if ( !p2 )
502: return 1;
503: else if ( NUM(p1) )
504: return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;
505: else if ( NUM(p2) )
506: return 1;
507: if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {
508: for ( dc1 = DC(p1), dc2 = DC(p2);
509: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )
510: switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
511: case 1:
512: return 1; break;
513: case -1:
514: return -1; break;
515: default:
516: switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {
517: case 1:
518: return 1; break;
519: case -1:
520: return -1; break;
521: default:
522: break;
523: }
524: break;
525: }
526: return dc1 ? 1 : (dc2 ? -1 : 0);
527: } else {
528: for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );
529: return v1 == VR(vl) ? 1 : -1;
530: }
531: }
532:
533: void csump(vl,p,c)
534: VL vl;
535: P p;
536: Q *c;
537: {
538: DCP dc;
539: Q s,s1,s2;
540:
541: if ( !p || NUM(p) )
542: *c = (Q)p;
543: else {
544: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
545: csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
546: }
547: *c = s;
548: }
549: }
550:
551: void degp(v,p,d)
552: V v;
553: P p;
554: Q *d;
555: {
556: DCP dc;
557: Q m,m1;
558:
559: if ( !p || NUM(p) )
560: *d = 0;
561: else if ( v == VR(p) )
562: *d = DEG(DC(p));
563: else {
564: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
565: degp(v,COEF(dc),&m1);
566: if ( cmpq(m,m1) < 0 )
567: m = m1;
568: }
569: *d = m;
570: }
571: }
572: #endif