![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to PDM.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 PDM.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/PDM.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_DIVSRP(vl,p1,p2,q,r)
10: VL vl;
11: P p1,p2;
12: P *q,*r;
13: {
14: register int i,j;
15: register DCP dc1,dc2,dc;
16: P m,s,dvr;
17: P *pq,*pr,*pd;
18: V v1,v2;
19: Q deg1,deg2;
20: int d1,d2,sgn;
21:
22: if ( !p1 ) {
23: *q = 0; *r = 0;
24: } else if ( NUM(p2) )
25: if ( NUM(p1) ) {
26: DIVNUM(p1,p2,q); *r = 0;
27: } else {
28: DIVSDCP(vl,p1,p2,q); *r = 0;
29: }
30: else if ( NUM(p1) ) {
31: *q = 0; *r = p1;
32: } else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
33: dc1 = DC(p1); dc2 = DC(p2);
34: deg1 = DEG(dc1); deg2 = DEG(dc2);
35: sgn = cmpq(deg1,deg2);
36: if ( sgn == 0 ) {
37: DIVSP(vl,COEF(dc1),COEF(dc2),q);
38: MULP(vl,p2,*q,&m); SUBP(vl,p1,m,r);
39: } else if ( sgn < 0 ) {
40: *q = 0; *r = p1;
41: } else {
42: if ( (PL(NM(deg1)) > 1) )
43: error("divsrp : invalid input");
44: d1 = QTOS(deg1); d2 = QTOS(deg2);
45: W_CALLOC(d1-d2,P,pq); W_CALLOC(d1,P,pr); W_CALLOC(d2,P,pd);
46: for ( dc = dc1; dc; dc = NEXT(dc) )
47: pr[QTOS(DEG(dc))] = COEF(dc);
48: for ( dc = dc2; dc; dc = NEXT(dc) )
49: pd[QTOS(DEG(dc))] = COEF(dc);
50: for ( dvr = COEF(dc2), i = d1 - d2; i >= 0; i-- ) {
51: if ( !pr[i+d2] )
52: continue;
53: DIVSP(vl,pr[i+d2],dvr,&pq[i]);
54: for ( j = d2; j >= 0; j-- ) {
55: MULP(vl,pq[i],pd[j],&m);
56: SUBP(vl,pr[i + j],m,&s); pr[i + j] = s;
57: }
58: }
59: plisttop(pq,v1,d1 - d2,q); plisttop(pr,v1,d1 - 1,r);
60: }
61: } else {
62: for ( ; (v1 != vl->v) && (v2 != vl->v); vl = NEXT(vl) );
63: if ( v2 == vl->v ) {
64: *q = 0; *r = p1;
65: } else
66: DIVSRDCP(vl,p1,p2,q,r);
67: }
68: }
69:
70: void D_DIVSRDCP(vl,p1,p2,q,r)
71: VL vl;
72: P p1,p2,*q,*r;
73: {
74:
75: P qc,rc;
76: DCP dc,dcq,dcq0,dcr,dcr0;
77:
78: for ( dc = DC(p1), dcq0 = 0, dcr0 = 0; dc; dc = NEXT(dc) ) {
79: DIVSRP(vl,COEF(dc),p2,&qc,&rc);
80: if ( qc ) {
81: NEXTDC(dcq0,dcq); DEG(dcq) = DEG(dc); COEF(dcq) = qc;
82: }
83: if ( rc ) {
84: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = rc;
85: }
86: }
87: if ( dcq0 ) {
88: NEXT(dcq) = 0; MKP(VR(p1),dcq0,*q);
89: } else
90: *q = 0;
91: if ( dcr0 ) {
92: NEXT(dcr) = 0; MKP(VR(p1),dcr0,*r);
93: } else
94: *r = 0;
95: }
96:
97: void D_DIVSP(vl,p1,p2,q)
98: VL vl;
99: P p1,p2,*q;
100: {
101: P t;
102:
103: DIVSRP(vl,p1,p2,q,&t);
104: if ( t )
105: error("divsp: cannot happen");
106: }
107:
108: void D_DIVSDCP(vl,p1,p2,q)
109: VL vl;
110: P p1,p2,*q;
111: {
112:
113: P m;
114: register DCP dc,dcr,dcr0;
115:
116: for ( dc = DC(p1), dcr0 = 0; dc; dc = NEXT(dc) ) {
117: DIVSP(vl,COEF(dc),p2,&m);
118: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = m; NEXT(dcr) = 0;
119: }
120: MKP(VR(p1),dcr0,*q);
121: }
122:
123: #ifdef FBASE
124: void plisttop(f,v,n,gp)
125: P *f;
126: V v;
127: int n;
128: P *gp;
129: {
130: int i;
131: DCP dc,dc0;
132:
133: for ( i = n; (i >= 0) && !f[i]; i-- );
134: if ( i < 0 )
135: *gp = 0;
136: else if ( i == 0 )
137: *gp = f[0];
138: else {
139: for ( dc0 = 0; i >= 0; i-- ) {
140: if ( !f[i] )
141: continue;
142: NEXTDC(dc0,dc);
143: if ( i )
144: STOQ(i,DEG(dc));
145: else
146: DEG(dc) = 0;
147: COEF(dc) = f[i];
148: }
149: NEXT(dc) = 0; MKP(v,dc0,*gp);
150: }
151: }
152:
153: int divtpz(vl,p1,p2,q)
154: VL vl;
155: P p1,p2,*q;
156: {
157: register int i,j,k;
158: register DCP dc1,dc2,dc;
159: P m,m1,s,dvr,t;
160: P *pq,*pr,*pd;
161: V v1,v2;
162: Q deg1,deg2;
163: int d1,d2,sgn;
164:
165: if ( !p1 ) {
166: *q = 0;
167: return ( 1 );
168: } else if ( NUM(p2) )
169: if ( NUM(p1) ) {
170: divq((Q)p1,(Q)p2,(Q *)&s);
171: if ( INT((Q)s) ) {
172: *q = s;
173: return ( 1 );
174: } else {
175: *q = 0;
176: return ( 0 );
177: }
178: } else
179: return ( divtdcpz(vl,p1,p2,q) );
180: else if ( NUM(p1) ) {
181: *q = 0;
182: return ( 0 );
183: } else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
184: Q csum1,csum2;
185:
186: csump(vl,p1,&csum1); csump(vl,p2,&csum2);
187: if ( csum2 && !divtpz(vl,(P)csum1,(P)csum2,&t) ) {
188: *q = 0;
189: return 0;
190: }
191: dc1 = DC(p1); dc2 = DC(p2);
192: deg1 = DEG(dc1); deg2 = DEG(dc2);
193: sgn = cmpq(deg1,deg2);
194: if ( sgn == 0 )
195: if ( !divtpz(vl,COEF(dc1),COEF(dc2),&m) ) {
196: *q = 0;
197: return ( 0 );
198: } else {
199: mulp(vl,p2,m,&m1); subp(vl,p1,m1,&s);
200: if ( !s ) {
201: *q = m;
202: return ( 1 );
203: } else {
204: *q = 0;
205: return ( 0 );
206: }
207: }
208: else if ( sgn < 0 ) {
209: *q = 0;
210: return ( 0 );
211: } else {
212: if ( (PL(NM(deg1)) > 1) ) {
213: error("divtpz : invalid input");
214: *q = 0;
215: return ( 0 );
216: }
217: d1 = QTOS(deg1); d2 = QTOS(deg2);
218: W_CALLOC(d1-d2,P,pq); W_CALLOC(d1,P,pr); W_CALLOC(d2,P,pd);
219: for ( dc = dc1; dc; dc = NEXT(dc) )
220: pr[QTOS(DEG(dc))] = COEF(dc);
221: for ( dc = dc2; dc; dc = NEXT(dc) )
222: pd[QTOS(DEG(dc))] = COEF(dc);
223: for ( dvr = COEF(dc2), i = d1 - d2; i >= 0; i-- )
224: if ( !pr[i+d2] )
225: continue;
226: else if ( !divtpz(vl,pr[i+d2],dvr,&m) ) {
227: *q = 0;
228: return ( 0 );
229: } else {
230: pq[i] = m;
231: for ( j = d2; j >= 0; j-- ) {
232: mulp(vl,pq[i],pd[j],&m);
233: subp(vl,pr[i + j],m,&s); pr[i + j] = s;
234: }
235: }
236: plisttop(pq,v1,d1 - d2,&m); plisttop(pr,v1,d1 - 1,&t);
237: if ( t ) {
238: *q = 0;
239: return ( 0 );
240: } else {
241: *q = m;
242: return ( 1 );
243: }
244: }
245: } else {
246: for ( ; (v1 != vl->v) && (v2 != vl->v); vl = NEXT(vl) );
247: if ( v2 == vl->v ) {
248: *q = 0;
249: return ( 0 );
250: } else
251: return ( divtdcpz(vl,p1,p2,q) ) ;
252: }
253: }
254:
255: int divtdcpz(vl,p1,p2,q)
256: VL vl;
257: P p1,p2,*q;
258: {
259:
260: P m;
261: register DCP dc,dcr,dcr0,dct;
262:
263: for ( dc = DC(p1), dcr0 = 0; dc; dc = NEXT(dc) )
264: if ( !divtpz(vl,COEF(dc),p2,&m) ) {
265: *q = 0;
266: return ( 0 );
267: } else {
268: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = m; NEXT(dcr) = 0;
269: }
270: MKP(VR(p1),dcr0,*q);
271: return ( 1 );
272: }
273:
274: void udivpz(f1,f2,fqp,frp)
275: P f1,f2,*fqp,*frp;
276: {
277: register int n1,n2,i,j;
278: Q *pq,*pr,*pd,d,m,s,q,r;
279: DCP dc;
280: N qn,rn;
281:
282: if ( !f2 )
283: error("udivpz: division by 0");
284: else if ( !f1 ) {
285: *fqp = *frp = 0;
286: return;
287: } else if ( NUM(f1) )
288: if ( NUM(f2) ) {
289: divn(NM((Q)f1),NM((Q)f2),&qn,&rn);
290: if ( rn ) {
291: *fqp = *frp = 0;
292: } else {
293: NTOQ(qn,SGN((Q)f1)*SGN((Q)f2),q),*fqp = (P)q; *frp = 0;
294: }
295: return;
296: } else {
297: *fqp = 0; *frp = f1;
298: return;
299: }
300: else if ( NUM(f2) ) {
301: n1 = UDEG(f1); W_CALLOC(n1,Q,pq);
302: for ( dc = DC(f1); dc; dc = NEXT(dc) ) {
303: divn(NM((Q)COEF(dc)),NM((Q)f2),&qn,&rn);
304: if ( rn ) {
305: *fqp = *frp = 0;
306: return;
307: } else {
308: NTOQ(qn,SGN((Q)COEF(dc))*SGN((Q)f2),s); pq[QTOS(DEG(dc))] = s;
309: }
310: }
311: plisttop((P *)pq,VR(f1),n1,fqp);
312: return;
313: }
314: n1 = UDEG(f1); n2 = UDEG(f2);
315: if ( n1 < n2 ) {
316: *fqp = NULL; *frp = f1;
317: return;
318: }
319: W_CALLOC(n1-n2,Q,pq); W_CALLOC(n1,Q,pr); W_CALLOC(n2,Q,pd);
320: for ( dc = DC(f1); dc; dc = NEXT(dc) )
321: pr[QTOS(DEG(dc))] = (Q)COEF(dc);
322: for ( dc = DC(f2); dc; dc = NEXT(dc) )
323: pd[QTOS(DEG(dc))] = (Q)COEF(dc);
324: for ( d = (Q)UCOEF(f2), i = n1 - n2; i >= 0; i-- ) {
325: if ( !pr[i+n2] )
326: continue;
327: divn(NM(pr[i+n2]),NM(d),&qn,&rn);
328: if ( rn ) {
329: *fqp = *frp = 0;
330: return;
331: }
332: NTOQ(qn,SGN(pr[i+n2])*SGN(d),pq[i]);
333: for ( j = n2; j >= 0; j-- ) {
334: mulq(pq[i],pd[j],&m); subq(pr[i+j],m,&s); pr[i+j] = s;
335: }
336: }
337: plisttop((P *)pq,VR(f1),n1-n2,fqp); plisttop((P *)pr,VR(f1),n2-1,frp);
338: }
339:
340: void udivpwm(mod,p1,p2,q,r)
341: Q mod;
342: P p1,p2;
343: P *q,*r;
344: {
345: P s,t,u,tq,tr;
346:
347: invl((Q)UCOEF(p2),mod,(Q *)&t); mulpq(p2,t,&s); cmp(mod,s,&u);
348: udivpzwm(mod,p1,u,&tq,&tr);
349: cmp(mod,tr,r); mulpq(tq,t,&s); cmp(mod,s,q);
350: }
351:
352: void udivpzwm(mod,f1,f2,fqp,frp)
353: Q mod;
354: P f1,f2,*fqp,*frp;
355: {
356: register int n1,n2,i,j;
357: Q *pq,*pr,*pd,d,m,s,q,r;
358: DCP dc;
359: N qn,rn;
360:
361: if ( !f2 )
362: error("udivpz: division by 0");
363: else if ( !f1 ) {
364: *fqp = *frp = 0;
365: return;
366: } else if ( NUM(f1) )
367: if ( NUM(f2) ) {
368: divn(NM((Q)f1),NM((Q)f2),&qn,&rn);
369: if ( rn ) {
370: *fqp = *frp = 0;
371: } else {
372: NTOQ(qn,SGN((Q)f1)*SGN((Q)f2),q),*fqp = (P)q; *frp = 0;
373: }
374: return;
375: } else {
376: *fqp = 0; *frp = f1;
377: return;
378: }
379: else if ( NUM(f2) ) {
380: n1 = UDEG(f1); W_CALLOC(n1,Q,pq);
381: for ( dc = DC(f1); dc; dc = NEXT(dc) ) {
382: divn(NM((Q)COEF(dc)),NM((Q)f2),&qn,&rn);
383: if ( rn ) {
384: *fqp = *frp = 0;
385: return;
386: } else {
387: NTOQ(qn,SGN((Q)COEF(dc))*SGN((Q)f2),s); pq[QTOS(DEG(dc))] = s;
388: }
389: }
390: plisttop((P *)pq,VR(f1),n1,fqp);
391: return;
392: }
393: n1 = UDEG(f1); n2 = UDEG(f2);
394: if ( n1 < n2 ) {
395: *fqp = NULL; *frp = f1;
396: return;
397: }
398: W_CALLOC(n1-n2,Q,pq); W_CALLOC(n1,Q,pr); W_CALLOC(n2,Q,pd);
399: for ( dc = DC(f1); dc; dc = NEXT(dc) )
400: pr[QTOS(DEG(dc))] = (Q)COEF(dc);
401: for ( dc = DC(f2); dc; dc = NEXT(dc) )
402: pd[QTOS(DEG(dc))] = (Q)COEF(dc);
403: for ( d = (Q)UCOEF(f2), i = n1 - n2; i >= 0; i-- ) {
404: if ( !pr[i+n2] )
405: continue;
406: divn(NM(pr[i+n2]),NM(d),&qn,&rn);
407: if ( rn ) {
408: *fqp = *frp = 0;
409: return;
410: }
411: NTOQ(qn,SGN(pr[i+n2])*SGN(d),pq[i]);
412: for ( j = n2; j >= 0; j-- ) {
413: mulq(pq[i],pd[j],&m); remq(m,mod,&s);
414: subq(pr[i+j],s,&m); remq(m,mod,&s); pr[i+j] = s;
415: }
416: }
417: plisttop((P *)pq,VR(f1),n1-n2,fqp); plisttop((P *)pr,VR(f1),n2-1,frp);
418: }
419: #endif