![[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.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/P.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #ifndef FBASE
3: #define FBASE
4: #endif
5:
6: #include "b.h"
7: #include "ca.h"
8:
9: void D_ADDP(vl,p1,p2,pr)
10: VL vl;
11: P p1,p2,*pr;
12: {
13: register DCP dc1,dc2,dcr0,dcr;
14: V v1,v2;
15: P t;
16:
17: if ( !p1 )
18: *pr = p2;
19: else if ( !p2 )
20: *pr = p1;
21: else if ( NUM(p1) )
22: if ( NUM(p2) )
23: ADDNUM(p1,p2,pr);
24: else
25: ADDPQ(p2,p1,pr);
26: else if ( NUM(p2) )
27: ADDPQ(p1,p2,pr);
28: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
29: for ( dc1 = DC(p1), dc2 = DC(p2), dcr0 = 0; dc1 && dc2; )
30: switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
31: case 0:
32: ADDP(vl,COEF(dc1),COEF(dc2),&t);
33: if ( t ) {
34: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc1); COEF(dcr) = t;
35: }
36: dc1 = NEXT(dc1); dc2 = NEXT(dc2); break;
37: case 1:
38: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc1); COEF(dcr) = COEF(dc1);
39: dc1 = NEXT(dc1); break;
40: case -1:
41: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc2); COEF(dcr) = COEF(dc2);
42: dc2 = NEXT(dc2); break;
43: }
44: if ( !dcr0 )
45: if ( dc1 )
46: dcr0 = dc1;
47: else if ( dc2 )
48: dcr0 = dc2;
49: else {
50: *pr = 0;
51: return;
52: }
53: else
54: if ( dc1 )
55: NEXT(dcr) = dc1;
56: else if ( dc2 )
57: NEXT(dcr) = dc2;
58: else
59: NEXT(dcr) = 0;
60: MKP(v1,dcr0,*pr);
61: } else {
62: while ( v1 != VR(vl) && v2 != VR(vl) )
63: vl = NEXT(vl);
64: if ( v1 == VR(vl) )
65: ADDPTOC(vl,p1,p2,pr);
66: else
67: ADDPTOC(vl,p2,p1,pr);
68: }
69: }
70:
71: void D_ADDPQ(p,q,pr)
72: P p,q,*pr;
73: {
74: DCP dc,dcr,dcr0;
75: P t;
76:
77: if ( NUM(p) )
78: ADDNUM(p,q,pr);
79: else {
80: for ( dcr0 = 0, dc = DC(p); dc && DEG(dc); dc = NEXT(dc) ) {
81: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = COEF(dc);
82: }
83: if ( !dc ) {
84: NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = q;
85: } else {
86: ADDPQ(COEF(dc),q,&t);
87: if ( t ) {
88: NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = t;
89: }
90: }
91: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
92: }
93: }
94:
95: void D_ADDPTOC(vl,p,c,pr)
96: VL vl;
97: P p,c,*pr;
98: {
99: DCP dc,dcr,dcr0;
100: P t;
101:
102: for ( dcr0 = 0, dc = DC(p); dc && DEG(dc); dc = NEXT(dc) ) {
103: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = COEF(dc);
104: }
105: if ( !dc ) {
106: NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = c;
107: } else {
108: ADDP(vl,COEF(dc),c,&t);
109: if ( t ) {
110: NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = t;
111: }
112: }
113: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
114: }
115:
116: void D_SUBP(vl,p1,p2,pr)
117: VL vl;
118: P p1,p2,*pr;
119: {
120: P t;
121:
122: if ( !p2 )
123: *pr = p1;
124: else {
125: CHSGNP(p2,&t); ADDP(vl,p1,t,pr);
126: }
127: }
128:
129: void D_MULP(vl,p1,p2,pr)
130: VL vl;
131: P p1,p2,*pr;
132: {
133: register DCP dc,dct,dcr,dcr0;
134: V v1,v2;
135: P t,s,u;
136: int n1,n2;
137:
138: if ( !p1 || !p2 ) *pr = 0;
139: else if ( NUM(p1) )
140: MULPQ(p2,p1,pr);
141: else if ( NUM(p2) )
142: MULPQ(p1,p2,pr);
143: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
144: for ( dc = DC(p1), n1 = 0; dc; dc = NEXT(dc), n1++ );
145: for ( dc = DC(p2), n2 = 0; dc; dc = NEXT(dc), n2++ );
146: if ( n1 > n2 )
147: for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {
148: for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {
149: NEXTDC(dcr0,dcr); MULP(vl,COEF(dct),COEF(dc),&COEF(dcr));
150: addq(DEG(dct),DEG(dc),&DEG(dcr));
151: }
152: NEXT(dcr) = 0; MKP(v1,dcr0,t);
153: ADDP(vl,s,t,&u); s = u; t = u = 0;
154: }
155: else
156: for ( dc = DC(p1), s = 0; dc; dc = NEXT(dc) ) {
157: for ( dcr0 = 0, dct = DC(p2); dct; dct = NEXT(dct) ) {
158: NEXTDC(dcr0,dcr); MULP(vl,COEF(dct),COEF(dc),&COEF(dcr));
159: addq(DEG(dct),DEG(dc),&DEG(dcr));
160: }
161: NEXT(dcr) = 0; MKP(v1,dcr0,t);
162: ADDP(vl,s,t,&u); s = u; t = u = 0;
163: }
164: *pr = s;
165: } else {
166: while ( v1 != VR(vl) && v2 != VR(vl) )
167: vl = NEXT(vl);
168: if ( v1 == VR(vl) )
169: MULPC(vl,p1,p2,pr);
170: else
171: MULPC(vl,p2,p1,pr);
172: }
173: }
174:
175: void D_MULPQ(p,q,pr)
176: P p,q,*pr;
177: {
178: DCP dc,dcr,dcr0;
179: P t;
180:
181: if (!p || !q)
182: *pr = 0;
183: else if ( Uniq(q) )
184: *pr = p;
185: else if ( NUM(p) )
186: MULNUM(p,q,pr);
187: else {
188: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
189: MULPQ(COEF(dc),q,&t);
190: if ( t ) {
191: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
192: }
193: }
194: if ( dcr0 ) {
195: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
196: } else
197: *pr = 0;
198: }
199: }
200:
201: void D_MULPC(vl,p,c,pr)
202: VL vl;
203: P p,c,*pr;
204: {
205: DCP dc,dcr,dcr0;
206: P t;
207:
208: if ( NUM(c) )
209: MULPQ(p,c,pr);
210: else {
211: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
212: MULP(vl,COEF(dc),c,&t);
213: if ( t ) {
214: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
215: }
216: }
217: if ( dcr0 ) {
218: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
219: } else
220: *pr = 0;
221: }
222: }
223:
224: void D_PWRP(vl,p,q,pr)
225: VL vl;
226: P p,*pr;
227: Q q;
228: {
229: DCP dc,dcr;
230: int n,i;
231: P *x,*y;
232: P t,s,u;
233: DCP dct;
234: P *pt;
235:
236: if ( !q ) {
237: *pr = (P)One;
238: } else if ( !p )
239: *pr = 0;
240: else if ( UNIQ(q) )
241: *pr = p;
242: else if ( NUM(p) )
243: PWRNUM(p,q,pr);
244: else {
245: dc = DC(p);
246: if ( !NEXT(dc) ) {
247: NEWDC(dcr);
248: PWRP(vl,COEF(dc),q,&COEF(dcr)); mulq(DEG(dc),q,&DEG(dcr));
249: NEXT(dcr) = 0; MKP(VR(p),dcr,*pr);
250: } else if ( !INT(q) ) {
251: error("pwrp: can't calculate fractional power."); *pr = 0;
252: } else if ( PL(NM(q)) == 1 ) {
253: n = QTOS(q); x = (P *)ALLOCA((n+1)*sizeof(pointer));
254: NEWDC(dct); DEG(dct) = DEG(dc); COEF(dct) = COEF(dc);
255: NEXT(dct) = 0; MKP(VR(p),dct,t);
256: for ( i = 0, u = (P)One; i < n; i++ ) {
257: x[i] = u; MULP(vl,u,t,&s); u = s;
258: }
259: x[n] = u; y = (P *)ALLOCA((n+1)*sizeof(pointer));
260: if ( DEG(NEXT(dc)) ) {
261: dct = NEXT(dc); MKP(VR(p),dct,t);
262: } else
263: t = COEF(NEXT(dc));
264: for ( i = 0, u = (P)One; i < n; i++ ) {
265: y[i] = u; MULP(vl,u,t,&s); u = s;
266: }
267: y[n] = u;
268: pt = (P *)ALLOCA((n+1)*sizeof(pointer)); MKBC(n,pt);
269: for ( i = 0, u = 0; i <= n; i++ ) {
270: MULP(vl,x[i],y[n-i],&t); MULP(vl,t,pt[i],&s);
271: ADDP(vl,u,s,&t); u = t;
272: }
273: *pr = u;
274: } else {
275: error("exponent too big");
276: *pr = 0;
277: }
278: }
279: }
280:
281: void D_CHSGNP(p,pr)
282: P p,*pr;
283: {
284: register DCP dc,dcr,dcr0;
285: P t;
286:
287: if ( !p )
288: *pr = NULL;
289: else if ( NUM(p) ) {
290: #if defined(THINK_C) || defined(_PA_RISC1_1) || defined(__alpha) || defined(mips)
291: #ifdef FBASE
292: chsgnnum((Num)p,(Num *)pr);
293: #else
294: MQ mq;
295:
296: NEWMQ(mq); CONT(mq)=mod-CONT((MQ)p); *pr = (P)mq;
297: #endif
298: #else
299: CHSGNNUM(p,t); *pr = t;
300: #endif
301: } else {
302: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
303: NEXTDC(dcr0,dcr); CHSGNP(COEF(dc),&COEF(dcr)); DEG(dcr) = DEG(dc);
304: }
305: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
306: }
307: }
308:
309:
310: #ifdef FBASE
311: void ptozp(p,sgn,c,pr)
312: P p;
313: int sgn;
314: Q *c;
315: P *pr;
316: {
317: N nm,dn;
318:
319: if ( !p ) {
320: *c = 0; *pr = 0;
321: } else {
322: lgp(p,&nm,&dn);
323: if ( UNIN(dn) )
324: NTOQ(nm,sgn,*c);
325: else
326: NDTOQ(nm,dn,sgn,*c);
327: divcp(p,*c,pr);
328: }
329: }
330:
331: void lgp(p,g,l)
332: P p;
333: N *g,*l;
334: {
335: N g1,g2,l1,l2,l3,l4,tmp;
336: DCP dc;
337:
338: if ( NUM(p) ) {
339: *g = NM((Q)p);
340: if ( INT((Q)p) )
341: *l = ONEN;
342: else
343: *l = DN((Q)p);
344: } else {
345: dc = DC(p); lgp(COEF(dc),g,l);
346: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
347: lgp(COEF(dc),&g1,&l1); gcdn(*g,g1,&g2); *g = g2;
348: gcdn(*l,l1,&l2); kmuln(*l,l1,&l3); divn(l3,l2,&l4,&tmp); *l = l4;
349: }
350: }
351: }
352:
353: void divcp(f,q,rp)
354: P f;
355: Q q;
356: P *rp;
357: {
358: DCP dc,dcr,dcr0;
359:
360: if ( !f )
361: *rp = 0;
362: else if ( NUM(f) )
363: DIVNUM(f,q,rp);
364: else {
365: for ( dc = DC(f), dcr0 = 0; dc; dc = NEXT(dc) ) {
366: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
367: divcp(COEF(dc),q,&COEF(dcr));
368: }
369: NEXT(dcr) = 0; MKP(VR(f),dcr0,*rp);
370: }
371: }
372:
373: void diffp(vl,p,v,r)
374: VL vl;
375: P p;
376: V v;
377: P *r;
378: {
379: P t;
380: DCP dc,dcr,dcr0;
381:
382: if ( !p || NUM(p) )
383: *r = 0;
384: else {
385: if ( v == VR(p) ) {
386: for ( dc = DC(p), dcr0 = 0;
387: dc && DEG(dc); dc = NEXT(dc) ) {
388: MULPQ(COEF(dc),(P)DEG(dc),&t);
389: if ( t ) {
390: NEXTDC(dcr0,dcr); SUBQ(DEG(dc),ONE,&DEG(dcr));
391: COEF(dcr) = t;
392: }
393: }
394: if ( !dcr0 )
395: *r = 0;
396: else {
397: NEXT(dcr) = 0; MKP(v,dcr0,*r);
398: }
399: } else {
400: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
401: diffp(vl,COEF(dc),v,&t);
402: if ( t ) {
403: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
404: }
405: }
406: if ( !dcr0 )
407: *r = 0;
408: else {
409: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
410: }
411: }
412: }
413: }
414:
415: void coefp(p,d,pr)
416: P p;
417: int d;
418: P *pr;
419: {
420: DCP dc;
421: int sgn;
422: Q dq;
423:
424: if ( NUM(p) )
425: if ( d == 0 )
426: *pr = p;
427: else
428: *pr = 0;
429: else {
430: for ( STOQ(d,dq), dc = DC(p); dc; dc = NEXT(dc) )
431: if ( (sgn = cmpq(DEG(dc),dq)) > 0 )
432: continue;
433: else if ( sgn == 0 ) {
434: *pr = COEF(dc);
435: return;
436: } else {
437: *pr = 0;
438: break;
439: }
440: *pr = 0;
441: }
442: }
443:
444: int compp(vl,p1,p2)
445: VL vl;
446: P p1,p2;
447: {
448: DCP dc1,dc2;
449: V v1,v2;
450:
451: if ( !p1 )
452: return p2 ? -1 : 0;
453: else if ( !p2 )
454: return 1;
455: else if ( NUM(p1) )
456: return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;
457: else if ( NUM(p2) )
458: return 1;
459: if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {
460: for ( dc1 = DC(p1), dc2 = DC(p2);
461: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )
462: switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
463: case 1:
464: return 1; break;
465: case -1:
466: return -1; break;
467: default:
468: switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {
469: case 1:
470: return 1; break;
471: case -1:
472: return -1; break;
473: default:
474: break;
475: }
476: break;
477: }
478: return dc1 ? 1 : (dc2 ? -1 : 0);
479: } else {
480: for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );
481: return v1 == VR(vl) ? 1 : -1;
482: }
483: }
484:
485: void csump(vl,p,c)
486: VL vl;
487: P p;
488: Q *c;
489: {
490: DCP dc;
491: Q s,s1,s2;
492:
493: if ( !p || NUM(p) )
494: *c = (Q)p;
495: else {
496: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
497: csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
498: }
499: *c = s;
500: }
501: }
502:
503: void degp(v,p,d)
504: V v;
505: P p;
506: Q *d;
507: {
508: DCP dc;
509: Q m,m1;
510:
511: if ( !p || NUM(p) )
512: *d = 0;
513: else if ( v == VR(p) )
514: *d = DEG(DC(p));
515: else {
516: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
517: degp(v,COEF(dc),&m1);
518: if ( cmpq(m,m1) < 0 )
519: m = m1;
520: }
521: *d = m;
522: }
523: }
524: #endif