![[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.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
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 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: *
1.5 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/PDM.c,v 1.4 2001/02/21 07:10:18 noro Exp $
1.2 noro 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_DIVSRP(vl,p1,p2,q,r)
58: VL vl;
59: P p1,p2;
60: P *q,*r;
61: {
62: register int i,j;
63: register DCP dc1,dc2,dc;
64: P m,s,dvr;
65: P *pq,*pr,*pd;
66: V v1,v2;
67: Q deg1,deg2;
68: int d1,d2,sgn;
69:
70: if ( !p1 ) {
71: *q = 0; *r = 0;
72: } else if ( NUM(p2) )
73: if ( NUM(p1) ) {
74: DIVNUM(p1,p2,q); *r = 0;
75: } else {
76: DIVSDCP(vl,p1,p2,q); *r = 0;
77: }
78: else if ( NUM(p1) ) {
79: *q = 0; *r = p1;
80: } else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
81: dc1 = DC(p1); dc2 = DC(p2);
82: deg1 = DEG(dc1); deg2 = DEG(dc2);
83: sgn = cmpq(deg1,deg2);
84: if ( sgn == 0 ) {
85: DIVSP(vl,COEF(dc1),COEF(dc2),q);
86: MULP(vl,p2,*q,&m); SUBP(vl,p1,m,r);
87: } else if ( sgn < 0 ) {
88: *q = 0; *r = p1;
89: } else {
90: if ( (PL(NM(deg1)) > 1) )
91: error("divsrp : invalid input");
92: d1 = QTOS(deg1); d2 = QTOS(deg2);
93: W_CALLOC(d1-d2,P,pq); W_CALLOC(d1,P,pr); W_CALLOC(d2,P,pd);
94: for ( dc = dc1; dc; dc = NEXT(dc) )
95: pr[QTOS(DEG(dc))] = COEF(dc);
96: for ( dc = dc2; dc; dc = NEXT(dc) )
97: pd[QTOS(DEG(dc))] = COEF(dc);
98: for ( dvr = COEF(dc2), i = d1 - d2; i >= 0; i-- ) {
99: if ( !pr[i+d2] )
100: continue;
101: DIVSP(vl,pr[i+d2],dvr,&pq[i]);
102: for ( j = d2; j >= 0; j-- ) {
103: MULP(vl,pq[i],pd[j],&m);
104: SUBP(vl,pr[i + j],m,&s); pr[i + j] = s;
105: }
106: }
107: plisttop(pq,v1,d1 - d2,q); plisttop(pr,v1,d1 - 1,r);
108: }
109: } else {
110: for ( ; (v1 != vl->v) && (v2 != vl->v); vl = NEXT(vl) );
111: if ( v2 == vl->v ) {
112: *q = 0; *r = p1;
113: } else
114: DIVSRDCP(vl,p1,p2,q,r);
115: }
116: }
117:
118: void D_DIVSRDCP(vl,p1,p2,q,r)
119: VL vl;
120: P p1,p2,*q,*r;
121: {
122:
123: P qc,rc;
124: DCP dc,dcq,dcq0,dcr,dcr0;
125:
126: for ( dc = DC(p1), dcq0 = 0, dcr0 = 0; dc; dc = NEXT(dc) ) {
127: DIVSRP(vl,COEF(dc),p2,&qc,&rc);
128: if ( qc ) {
129: NEXTDC(dcq0,dcq); DEG(dcq) = DEG(dc); COEF(dcq) = qc;
130: }
131: if ( rc ) {
132: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = rc;
133: }
134: }
135: if ( dcq0 ) {
136: NEXT(dcq) = 0; MKP(VR(p1),dcq0,*q);
137: } else
138: *q = 0;
139: if ( dcr0 ) {
140: NEXT(dcr) = 0; MKP(VR(p1),dcr0,*r);
141: } else
142: *r = 0;
143: }
144:
145: void D_DIVSP(vl,p1,p2,q)
146: VL vl;
147: P p1,p2,*q;
148: {
149: P t;
150:
151: DIVSRP(vl,p1,p2,q,&t);
152: if ( t )
153: error("divsp: cannot happen");
154: }
155:
156: void D_DIVSDCP(vl,p1,p2,q)
157: VL vl;
158: P p1,p2,*q;
159: {
160:
161: P m;
162: register DCP dc,dcr,dcr0;
163:
164: for ( dc = DC(p1), dcr0 = 0; dc; dc = NEXT(dc) ) {
165: DIVSP(vl,COEF(dc),p2,&m);
166: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = m; NEXT(dcr) = 0;
167: }
168: MKP(VR(p1),dcr0,*q);
169: }
170:
171: #ifdef FBASE
172: void plisttop(f,v,n,gp)
173: P *f;
174: V v;
175: int n;
176: P *gp;
177: {
178: int i;
179: DCP dc,dc0;
180:
181: for ( i = n; (i >= 0) && !f[i]; i-- );
182: if ( i < 0 )
183: *gp = 0;
184: else if ( i == 0 )
185: *gp = f[0];
186: else {
187: for ( dc0 = 0; i >= 0; i-- ) {
188: if ( !f[i] )
189: continue;
190: NEXTDC(dc0,dc);
191: if ( i )
192: STOQ(i,DEG(dc));
193: else
194: DEG(dc) = 0;
195: COEF(dc) = f[i];
196: }
197: NEXT(dc) = 0; MKP(v,dc0,*gp);
198: }
199: }
200:
201: int divtpz(vl,p1,p2,q)
202: VL vl;
203: P p1,p2,*q;
204: {
205: register int i,j,k;
206: register DCP dc1,dc2,dc;
207: P m,m1,s,dvr,t;
208: P *pq,*pr,*pd;
209: V v1,v2;
210: Q deg1,deg2;
211: int d1,d2,sgn;
212:
213: if ( !p1 ) {
214: *q = 0;
215: return ( 1 );
216: } else if ( NUM(p2) )
217: if ( NUM(p1) ) {
218: divq((Q)p1,(Q)p2,(Q *)&s);
219: if ( INT((Q)s) ) {
220: *q = s;
221: return ( 1 );
222: } else {
223: *q = 0;
224: return ( 0 );
225: }
226: } else
227: return ( divtdcpz(vl,p1,p2,q) );
228: else if ( NUM(p1) ) {
229: *q = 0;
230: return ( 0 );
231: } else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
232: Q csum1,csum2;
233:
234: csump(vl,p1,&csum1); csump(vl,p2,&csum2);
235: if ( csum2 && !divtpz(vl,(P)csum1,(P)csum2,&t) ) {
236: *q = 0;
237: return 0;
238: }
239: dc1 = DC(p1); dc2 = DC(p2);
240: deg1 = DEG(dc1); deg2 = DEG(dc2);
241: sgn = cmpq(deg1,deg2);
242: if ( sgn == 0 )
243: if ( !divtpz(vl,COEF(dc1),COEF(dc2),&m) ) {
244: *q = 0;
245: return ( 0 );
246: } else {
247: mulp(vl,p2,m,&m1); subp(vl,p1,m1,&s);
248: if ( !s ) {
249: *q = m;
250: return ( 1 );
251: } else {
252: *q = 0;
253: return ( 0 );
254: }
255: }
256: else if ( sgn < 0 ) {
257: *q = 0;
258: return ( 0 );
259: } else {
260: if ( (PL(NM(deg1)) > 1) ) {
261: error("divtpz : invalid input");
262: *q = 0;
263: return ( 0 );
264: }
265: d1 = QTOS(deg1); d2 = QTOS(deg2);
266: W_CALLOC(d1-d2,P,pq); W_CALLOC(d1,P,pr); W_CALLOC(d2,P,pd);
267: for ( dc = dc1; dc; dc = NEXT(dc) )
268: pr[QTOS(DEG(dc))] = COEF(dc);
269: for ( dc = dc2; dc; dc = NEXT(dc) )
270: pd[QTOS(DEG(dc))] = COEF(dc);
271: for ( dvr = COEF(dc2), i = d1 - d2; i >= 0; i-- )
272: if ( !pr[i+d2] )
273: continue;
274: else if ( !divtpz(vl,pr[i+d2],dvr,&m) ) {
275: *q = 0;
276: return ( 0 );
277: } else {
278: pq[i] = m;
279: for ( j = d2; j >= 0; j-- ) {
280: mulp(vl,pq[i],pd[j],&m);
281: subp(vl,pr[i + j],m,&s); pr[i + j] = s;
282: }
283: }
284: plisttop(pq,v1,d1 - d2,&m); plisttop(pr,v1,d1 - 1,&t);
285: if ( t ) {
286: *q = 0;
287: return ( 0 );
288: } else {
289: *q = m;
290: return ( 1 );
291: }
292: }
293: } else {
294: for ( ; (v1 != vl->v) && (v2 != vl->v); vl = NEXT(vl) );
295: if ( v2 == vl->v ) {
296: *q = 0;
297: return ( 0 );
298: } else
299: return ( divtdcpz(vl,p1,p2,q) ) ;
300: }
301: }
302:
303: int divtdcpz(vl,p1,p2,q)
304: VL vl;
305: P p1,p2,*q;
306: {
307:
308: P m;
309: register DCP dc,dcr,dcr0,dct;
310:
311: for ( dc = DC(p1), dcr0 = 0; dc; dc = NEXT(dc) )
312: if ( !divtpz(vl,COEF(dc),p2,&m) ) {
313: *q = 0;
314: return ( 0 );
315: } else {
316: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = m; NEXT(dcr) = 0;
317: }
318: MKP(VR(p1),dcr0,*q);
319: return ( 1 );
320: }
321:
322: void udivpz(f1,f2,fqp,frp)
323: P f1,f2,*fqp,*frp;
324: {
325: register int n1,n2,i,j;
326: Q *pq,*pr,*pd,d,m,s,q,r;
327: DCP dc;
328: N qn,rn;
329:
330: if ( !f2 )
331: error("udivpz: division by 0");
332: else if ( !f1 ) {
333: *fqp = *frp = 0;
334: return;
335: } else if ( NUM(f1) )
336: if ( NUM(f2) ) {
337: divn(NM((Q)f1),NM((Q)f2),&qn,&rn);
338: if ( rn ) {
339: *fqp = *frp = 0;
340: } else {
341: NTOQ(qn,SGN((Q)f1)*SGN((Q)f2),q),*fqp = (P)q; *frp = 0;
342: }
343: return;
344: } else {
345: *fqp = 0; *frp = f1;
346: return;
347: }
348: else if ( NUM(f2) ) {
349: n1 = UDEG(f1); W_CALLOC(n1,Q,pq);
350: for ( dc = DC(f1); dc; dc = NEXT(dc) ) {
351: divn(NM((Q)COEF(dc)),NM((Q)f2),&qn,&rn);
352: if ( rn ) {
353: *fqp = *frp = 0;
354: return;
355: } else {
356: NTOQ(qn,SGN((Q)COEF(dc))*SGN((Q)f2),s); pq[QTOS(DEG(dc))] = s;
357: }
358: }
359: plisttop((P *)pq,VR(f1),n1,fqp);
360: return;
361: }
362: n1 = UDEG(f1); n2 = UDEG(f2);
363: if ( n1 < n2 ) {
364: *fqp = NULL; *frp = f1;
365: return;
366: }
367: W_CALLOC(n1-n2,Q,pq); W_CALLOC(n1,Q,pr); W_CALLOC(n2,Q,pd);
368: for ( dc = DC(f1); dc; dc = NEXT(dc) )
369: pr[QTOS(DEG(dc))] = (Q)COEF(dc);
370: for ( dc = DC(f2); dc; dc = NEXT(dc) )
371: pd[QTOS(DEG(dc))] = (Q)COEF(dc);
372: for ( d = (Q)UCOEF(f2), i = n1 - n2; i >= 0; i-- ) {
373: if ( !pr[i+n2] )
374: continue;
375: divn(NM(pr[i+n2]),NM(d),&qn,&rn);
376: if ( rn ) {
377: *fqp = *frp = 0;
378: return;
379: }
380: NTOQ(qn,SGN(pr[i+n2])*SGN(d),pq[i]);
381: for ( j = n2; j >= 0; j-- ) {
382: mulq(pq[i],pd[j],&m); subq(pr[i+j],m,&s); pr[i+j] = s;
383: }
384: }
385: plisttop((P *)pq,VR(f1),n1-n2,fqp); plisttop((P *)pr,VR(f1),n2-1,frp);
386: }
387:
388: void udivpwm(mod,p1,p2,q,r)
389: Q mod;
390: P p1,p2;
391: P *q,*r;
392: {
393: P s,t,u,tq,tr;
394:
395: invl((Q)UCOEF(p2),mod,(Q *)&t); mulpq(p2,t,&s); cmp(mod,s,&u);
396: udivpzwm(mod,p1,u,&tq,&tr);
397: cmp(mod,tr,r); mulpq(tq,t,&s); cmp(mod,s,q);
398: }
399:
400: void udivpzwm(mod,f1,f2,fqp,frp)
401: Q mod;
402: P f1,f2,*fqp,*frp;
403: {
404: register int n1,n2,i,j;
405: Q *pq,*pr,*pd,d,m,s,q,r;
406: DCP dc;
407: N qn,rn;
408:
409: if ( !f2 )
410: error("udivpz: division by 0");
411: else if ( !f1 ) {
412: *fqp = *frp = 0;
413: return;
414: } else if ( NUM(f1) )
415: if ( NUM(f2) ) {
416: divn(NM((Q)f1),NM((Q)f2),&qn,&rn);
417: if ( rn ) {
418: *fqp = *frp = 0;
419: } else {
420: NTOQ(qn,SGN((Q)f1)*SGN((Q)f2),q),*fqp = (P)q; *frp = 0;
421: }
422: return;
423: } else {
424: *fqp = 0; *frp = f1;
425: return;
426: }
427: else if ( NUM(f2) ) {
428: n1 = UDEG(f1); W_CALLOC(n1,Q,pq);
429: for ( dc = DC(f1); dc; dc = NEXT(dc) ) {
430: divn(NM((Q)COEF(dc)),NM((Q)f2),&qn,&rn);
431: if ( rn ) {
432: *fqp = *frp = 0;
433: return;
434: } else {
435: NTOQ(qn,SGN((Q)COEF(dc))*SGN((Q)f2),s); pq[QTOS(DEG(dc))] = s;
436: }
437: }
438: plisttop((P *)pq,VR(f1),n1,fqp);
439: return;
440: }
441: n1 = UDEG(f1); n2 = UDEG(f2);
442: if ( n1 < n2 ) {
443: *fqp = NULL; *frp = f1;
444: return;
445: }
446: W_CALLOC(n1-n2,Q,pq); W_CALLOC(n1,Q,pr); W_CALLOC(n2,Q,pd);
447: for ( dc = DC(f1); dc; dc = NEXT(dc) )
448: pr[QTOS(DEG(dc))] = (Q)COEF(dc);
449: for ( dc = DC(f2); dc; dc = NEXT(dc) )
450: pd[QTOS(DEG(dc))] = (Q)COEF(dc);
451: for ( d = (Q)UCOEF(f2), i = n1 - n2; i >= 0; i-- ) {
452: if ( !pr[i+n2] )
453: continue;
454: divn(NM(pr[i+n2]),NM(d),&qn,&rn);
455: if ( rn ) {
456: *fqp = *frp = 0;
457: return;
458: }
459: NTOQ(qn,SGN(pr[i+n2])*SGN(d),pq[i]);
460: for ( j = n2; j >= 0; j-- ) {
461: mulq(pq[i],pd[j],&m); remq(m,mod,&s);
462: subq(pr[i+j],s,&m); remq(m,mod,&s); pr[i+j] = s;
463: }
464: }
465: plisttop((P *)pq,VR(f1),n1-n2,fqp); plisttop((P *)pr,VR(f1),n2-1,frp);
466: }
467: #endif
1.4 noro 468:
469: #ifdef MODULAR
470:
1.5 ! noro 471: int divtmp(VL vl,int mod,P p1,P p2,P *q)
1.4 noro 472: {
1.5 ! noro 473: register int i,j;
1.4 noro 474: register DCP dc1,dc2,dc;
475: P m,m1,s,dvr,t;
476: P *pq,*pr,*pd;
477: V v1,v2;
478: Q deg1,deg2;
479: int d1,d2,sgn;
480:
481: if ( !p1 ) {
482: *q = 0;
483: return 1;
484: } else if ( NUM(p2) ) {
485: divsmp(vl,mod,p1,p2,q);
486: return 1;
487: } else if ( NUM(p1) ) {
488: *q = 0;
489: return 0;
490: } else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
491: dc1 = DC(p1); dc2 = DC(p2);
492: deg1 = DEG(dc1); deg2 = DEG(dc2);
493: sgn = cmpq(deg1,deg2);
494: if ( sgn == 0 )
495: if ( !divtmp(vl,mod,COEF(dc1),COEF(dc2),&m) ) {
496: *q = 0;
497: return 0;
498: } else {
499: mulmp(vl,mod,p2,m,&m1); submp(vl,mod,p1,m1,&s);
500: if ( !s ) {
501: *q = m;
502: return 1;
503: } else {
504: *q = 0;
505: return 0;
506: }
507: }
508: else if ( sgn < 0 ) {
509: *q = 0;
510: return 0;
511: } else {
512: if ( (PL(NM(deg1)) > 1) ) {
513: error("divtmp : invalid input");
514: *q = 0;
515: return ( 0 );
516: }
517: d1 = QTOS(deg1); d2 = QTOS(deg2);
518: W_CALLOC(d1-d2,P,pq); W_CALLOC(d1,P,pr); W_CALLOC(d2,P,pd);
519: for ( dc = dc1; dc; dc = NEXT(dc) )
520: pr[QTOS(DEG(dc))] = COEF(dc);
521: for ( dc = dc2; dc; dc = NEXT(dc) )
522: pd[QTOS(DEG(dc))] = COEF(dc);
523: for ( dvr = COEF(dc2), i = d1 - d2; i >= 0; i-- )
524: if ( !pr[i+d2] )
525: continue;
526: else if ( !divtmp(vl,mod,pr[i+d2],dvr,&m) ) {
527: *q = 0;
528: return 0;
529: } else {
530: pq[i] = m;
531: for ( j = d2; j >= 0; j-- ) {
532: mulmp(vl,mod,pq[i],pd[j],&m);
533: submp(vl,mod,pr[i + j],m,&s); pr[i + j] = s;
534: }
535: }
536: plisttop(pq,v1,d1 - d2,&m); plisttop(pr,v1,d1 - 1,&t);
537: if ( t ) {
538: *q = 0;
539: return 0;
540: } else {
541: *q = m;
542: return 1;
543: }
544: }
545: } else {
546: for ( ; (v1 != vl->v) && (v2 != vl->v); vl = NEXT(vl) );
547: if ( v2 == vl->v ) {
548: *q = 0;
549: return 0;
550: } else
551: return divtdcmp(vl,mod,p1,p2,q);
552: }
553: }
554:
1.5 ! noro 555: int divtdcmp(VL vl,int mod,P p1,P p2,P *q)
1.4 noro 556: {
557:
558: P m;
1.5 ! noro 559: register DCP dc,dcr,dcr0;
1.4 noro 560:
561: for ( dc = DC(p1), dcr0 = 0; dc; dc = NEXT(dc) )
562: if ( !divtmp(vl,mod,COEF(dc),p2,&m) ) {
563: *q = 0;
564: return 0;
565: } else {
566: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = m; NEXT(dcr) = 0;
567: }
568: MKP(VR(p1),dcr0,*q);
569: return 1;
570: }
571: #endif /* MODULAR */