![[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.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/PM.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #ifndef MODULAR
3: #define MODULAR
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:
180: if (!p || !q)
181: *pr = 0;
182: else if ( Uniq(q) )
183: *pr = p;
184: else if ( NUM(p) )
185: MULNUM(p,q,pr);
186: else {
187: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
188: NEXTDC(dcr0,dcr); MULPQ(COEF(dc),q,&COEF(dcr)); DEG(dcr) = DEG(dc);
189: }
190: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
191: }
192: }
193:
194: void D_MULPC(vl,p,c,pr)
195: VL vl;
196: P p,c,*pr;
197: {
198: DCP dc,dcr,dcr0;
199:
200: if ( NUM(c) )
201: MULPQ(p,c,pr);
202: else {
203: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
204: NEXTDC(dcr0,dcr); MULP(vl,COEF(dc),c,&COEF(dcr)); DEG(dcr) = DEG(dc);
205: }
206: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
207: }
208: }
209:
210: void D_PWRP(vl,p,q,pr)
211: VL vl;
212: P p,*pr;
213: Q q;
214: {
215: register DCP dc,dcr;
216: int n,i;
217: P *x,*y;
218: P t,s,u;
219: DCP dct;
220: P *pt;
221:
222: if ( !q ) {
223: *pr = (P)One;
224: } else if ( !p )
225: *pr = 0;
226: else if ( UNIQ(q) )
227: *pr = p;
228: else if ( NUM(p) )
229: PWRNUM(p,q,pr);
230: else {
231: dc = DC(p);
232: if ( !NEXT(dc) ) {
233: NEWDC(dcr);
234: PWRP(vl,COEF(dc),q,&COEF(dcr)); mulq(DEG(dc),q,&DEG(dcr));
235: NEXT(dcr) = 0; MKP(VR(p),dcr,*pr);
236: } else if ( !INT(q) ) {
237: error("pwrp: can't calculate fractional power."); *pr = 0;
238: } else if ( PL(NM(q)) == 1 ) {
239: n = QTOS(q); x = (P *)ALLOCA((n+1)*sizeof(pointer));
240: NEWDC(dct); DEG(dct) = DEG(dc); COEF(dct) = COEF(dc);
241: NEXT(dct) = 0; MKP(VR(p),dct,t);
242: for ( i = 0, u = (P)One; i < n; i++ ) {
243: x[i] = u; MULP(vl,u,t,&s); u = s;
244: }
245: x[n] = u; y = (P *)ALLOCA((n+1)*sizeof(pointer));
246: if ( DEG(NEXT(dc)) ) {
247: dct = NEXT(dc); MKP(VR(p),dct,t);
248: } else
249: t = COEF(NEXT(dc));
250: for ( i = 0, u = (P)One; i < n; i++ ) {
251: y[i] = u; MULP(vl,u,t,&s); u = s;
252: }
253: y[n] = u;
254: pt = (P *)ALLOCA((n+1)*sizeof(pointer)); MKBC(n,pt);
255: for ( i = 0, u = 0; i <= n; i++ ) {
256: MULP(vl,x[i],y[n-i],&t); MULP(vl,t,pt[i],&s);
257: ADDP(vl,u,s,&t); u = t;
258: }
259: *pr = u;
260: } else {
261: error("exponent too big");
262: *pr = 0;
263: }
264: }
265: }
266:
267: void D_CHSGNP(p,pr)
268: P p,*pr;
269: {
270: register DCP dc,dcr,dcr0;
271: P t;
272:
273: if ( !p )
274: *pr = NULL;
275: else if ( NUM(p) ) {
276: #if defined(THINK_C) || defined(_PA_RISC1_1) || defined(__alpha) || defined(mips)
277: #ifdef FBASE
278: chsgnnum((Num)p,(Num *)pr);
279: #else
280: MQ mq;
281:
282: NEWMQ(mq); CONT(mq)=mod-CONT((MQ)p); *pr = (P)mq;
283: #endif
284: #else
285: CHSGNNUM(p,t); *pr = t;
286: #endif
287: } else {
288: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
289: NEXTDC(dcr0,dcr); CHSGNP(COEF(dc),&COEF(dcr)); DEG(dcr) = DEG(dc);
290: }
291: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
292: }
293: }
294:
295:
296: #ifdef FBASE
297: void ptozp(p,sgn,c,pr)
298: P p;
299: int sgn;
300: Q *c;
301: P *pr;
302: {
303: N nm,dn;
304:
305: if ( !p ) {
306: *c = 0; *pr = 0;
307: } else {
308: lgp(p,&nm,&dn);
309: if ( UNIN(dn) )
310: NTOQ(nm,sgn,*c);
311: else
312: NDTOQ(nm,dn,sgn,*c);
313: divcp(p,*c,pr);
314: }
315: }
316:
317: void lgp(p,g,l)
318: P p;
319: N *g,*l;
320: {
321: N g1,g2,l1,l2,l3,l4,tmp;
322: DCP dc;
323:
324: if ( NUM(p) ) {
325: *g = NM((Q)p);
326: if ( INT((Q)p) )
327: *l = ONEN;
328: else
329: *l = DN((Q)p);
330: } else {
331: dc = DC(p); lgp(COEF(dc),g,l);
332: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
333: lgp(COEF(dc),&g1,&l1); gcdn(*g,g1,&g2); *g = g2;
334: gcdn(*l,l1,&l2); kmuln(*l,l1,&l3); divn(l3,l2,&l4,&tmp); *l = l4;
335: }
336: }
337: }
338:
339: void divcp(f,q,rp)
340: P f;
341: Q q;
342: P *rp;
343: {
344: DCP dc,dcr,dcr0;
345:
346: if ( !f )
347: *rp = 0;
348: else if ( NUM(f) )
349: DIVNUM(f,q,rp);
350: else {
351: for ( dc = DC(f), dcr0 = 0; dc; dc = NEXT(dc) ) {
352: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
353: divcp(COEF(dc),q,&COEF(dcr));
354: }
355: NEXT(dcr) = 0; MKP(VR(f),dcr0,*rp);
356: }
357: }
358:
359: void diffp(vl,p,v,r)
360: VL vl;
361: P p;
362: V v;
363: P *r;
364: {
365: P t;
366: DCP dc,dcr,dcr0;
367:
368: if ( !p || NUM(p) )
369: *r = 0;
370: else {
371: if ( v == VR(p) ) {
372: for ( dc = DC(p), dcr0 = 0;
373: dc && DEG(dc); dc = NEXT(dc) ) {
374: NEXTDC(dcr0,dcr); SUBQ(DEG(dc),ONE,&DEG(dcr));
375: MULPQ(COEF(dc),(P)DEG(dc),&COEF(dcr));
376: }
377: NEXT(dcr) = 0; MKP(v,dcr0,*r);
378: } else {
379: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
380: diffp(vl,COEF(dc),v,&t);
381: if ( t ) {
382: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
383: }
384: }
385: if ( !dcr0 )
386: *r = 0;
387: else {
388: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
389: }
390: }
391: }
392: }
393:
394: void coefp(p,d,pr)
395: P p;
396: int d;
397: P *pr;
398: {
399: DCP dc;
400: int sgn;
401: Q dq;
402:
403: if ( NUM(p) )
404: if ( d == 0 )
405: *pr = p;
406: else
407: *pr = 0;
408: else {
409: for ( STOQ(d,dq), dc = DC(p); dc; dc = NEXT(dc) )
410: if ( (sgn = cmpq(DEG(dc),dq)) > 0 )
411: continue;
412: else if ( sgn == 0 ) {
413: *pr = COEF(dc);
414: return;
415: } else {
416: *pr = 0;
417: break;
418: }
419: *pr = 0;
420: }
421: }
422:
423: int compp(vl,p1,p2)
424: VL vl;
425: P p1,p2;
426: {
427: DCP dc1,dc2;
428: V v1,v2;
429:
430: if ( !p1 )
431: return p2 ? -1 : 0;
432: else if ( !p2 )
433: return 1;
434: else if ( NUM(p1) )
435: return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;
436: else if ( NUM(p2) )
437: return 1;
438: if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {
439: for ( dc1 = DC(p1), dc2 = DC(p2);
440: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )
441: switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
442: case 1:
443: return 1; break;
444: case -1:
445: return -1; break;
446: default:
447: switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {
448: case 1:
449: return 1; break;
450: case -1:
451: return -1; break;
452: default:
453: break;
454: }
455: break;
456: }
457: return dc1 ? 1 : (dc2 ? -1 : 0);
458: } else {
459: for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );
460: return v1 == VR(vl) ? 1 : -1;
461: }
462: }
463:
464: void csump(vl,p,c)
465: VL vl;
466: P p;
467: Q *c;
468: {
469: DCP dc;
470: Q s,s1,s2;
471:
472: if ( !p || NUM(p) )
473: *c = (Q)p;
474: else {
475: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
476: csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
477: }
478: *c = s;
479: }
480: }
481:
482: void degp(v,p,d)
483: V v;
484: P p;
485: Q *d;
486: {
487: DCP dc;
488: Q m,m1;
489:
490: if ( !p || NUM(p) )
491: *d = 0;
492: else if ( v == VR(p) )
493: *d = DEG(DC(p));
494: else {
495: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
496: degp(v,COEF(dc),&m1);
497: if ( cmpq(m,m1) < 0 )
498: m = m1;
499: }
500: *d = m;
501: }
502: }
503: #endif