![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to PD.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 PD.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.4 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/PD.c,v 1.3 2000/08/22 05:04:04 noro Exp $
1.2 noro 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_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: }
1.4 ! noro 199: }
! 200:
! 201: /* for multivariate polynomials over fields */
! 202:
! 203: int divtp(vl,p1,p2,q)
! 204: VL vl;
! 205: P p1,p2,*q;
! 206: {
! 207: register int i,j,k;
! 208: register DCP dc1,dc2,dc;
! 209: P m,m1,s,dvr,t;
! 210: P *pq,*pr,*pd;
! 211: V v1,v2;
! 212: Q deg1,deg2;
! 213: int d1,d2,sgn;
! 214:
! 215: if ( !p1 ) {
! 216: *q = 0;
! 217: return 1;
! 218: } else if ( NUM(p2) ) {
! 219: divsp(vl,p1,p2,q);
! 220: return 1;
! 221: } else if ( NUM(p1) ) {
! 222: *q = 0;
! 223: return 0;
! 224: } else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
! 225: dc1 = DC(p1); dc2 = DC(p2);
! 226: deg1 = DEG(dc1); deg2 = DEG(dc2);
! 227: sgn = cmpq(deg1,deg2);
! 228: if ( sgn == 0 )
! 229: if ( !divtp(vl,COEF(dc1),COEF(dc2),&m) ) {
! 230: *q = 0;
! 231: return 0;
! 232: } else {
! 233: mulp(vl,p2,m,&m1); subp(vl,p1,m1,&s);
! 234: if ( !s ) {
! 235: *q = m;
! 236: return 1;
! 237: } else {
! 238: *q = 0;
! 239: return 0;
! 240: }
! 241: }
! 242: else if ( sgn < 0 ) {
! 243: *q = 0;
! 244: return 0;
! 245: } else {
! 246: if ( (PL(NM(deg1)) > 1) ) {
! 247: error("divtp : invalid input");
! 248: *q = 0;
! 249: return ( 0 );
! 250: }
! 251: d1 = QTOS(deg1); d2 = QTOS(deg2);
! 252: W_CALLOC(d1-d2,P,pq); W_CALLOC(d1,P,pr); W_CALLOC(d2,P,pd);
! 253: for ( dc = dc1; dc; dc = NEXT(dc) )
! 254: pr[QTOS(DEG(dc))] = COEF(dc);
! 255: for ( dc = dc2; dc; dc = NEXT(dc) )
! 256: pd[QTOS(DEG(dc))] = COEF(dc);
! 257: for ( dvr = COEF(dc2), i = d1 - d2; i >= 0; i-- )
! 258: if ( !pr[i+d2] )
! 259: continue;
! 260: else if ( !divtp(vl,pr[i+d2],dvr,&m) ) {
! 261: *q = 0;
! 262: return 0;
! 263: } else {
! 264: pq[i] = m;
! 265: for ( j = d2; j >= 0; j-- ) {
! 266: mulp(vl,pq[i],pd[j],&m);
! 267: subp(vl,pr[i + j],m,&s); pr[i + j] = s;
! 268: }
! 269: }
! 270: plisttop(pq,v1,d1 - d2,&m); plisttop(pr,v1,d1 - 1,&t);
! 271: if ( t ) {
! 272: *q = 0;
! 273: return 0;
! 274: } else {
! 275: *q = m;
! 276: return 1;
! 277: }
! 278: }
! 279: } else {
! 280: for ( ; (v1 != vl->v) && (v2 != vl->v); vl = NEXT(vl) );
! 281: if ( v2 == vl->v ) {
! 282: *q = 0;
! 283: return 0;
! 284: } else
! 285: return divtdcp(vl,p1,p2,q);
! 286: }
! 287: }
! 288:
! 289: int divtdcp(vl,p1,p2,q)
! 290: VL vl;
! 291: P p1,p2,*q;
! 292: {
! 293:
! 294: P m;
! 295: register DCP dc,dcr,dcr0,dct;
! 296:
! 297: for ( dc = DC(p1), dcr0 = 0; dc; dc = NEXT(dc) )
! 298: if ( !divtp(vl,COEF(dc),p2,&m) ) {
! 299: *q = 0;
! 300: return 0;
! 301: } else {
! 302: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = m; NEXT(dcr) = 0;
! 303: }
! 304: MKP(VR(p1),dcr0,*q);
! 305: return 1;
1.1 noro 306: }
307:
308: int divtpz(vl,p1,p2,q)
309: VL vl;
310: P p1,p2,*q;
311: {
312: register int i,j;
313: register DCP dc1,dc2,dc;
314: P m,m1,s,dvr,t;
315: P *pq,*pr,*pd;
316: V v1,v2;
317: Q deg1,deg2;
318: int d1,d2,sgn;
319:
320: if ( !p1 ) {
321: *q = 0;
322: return ( 1 );
323: } else if ( NUM(p2) )
324: if ( NUM(p1) ) {
325: divq((Q)p1,(Q)p2,(Q *)&s);
326: if ( INT((Q)s) ) {
327: *q = s;
328: return ( 1 );
329: } else {
330: *q = 0;
331: return ( 0 );
332: }
333: } else
334: return ( divtdcpz(vl,p1,p2,q) );
335: else if ( NUM(p1) ) {
336: *q = 0;
337: return ( 0 );
338: } else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
339: Q csum1,csum2;
340:
341: csump(vl,p1,&csum1); csump(vl,p2,&csum2);
342: if ( csum2 && !divtpz(vl,(P)csum1,(P)csum2,&t) ) {
343: *q = 0;
344: return 0;
345: }
346: dc1 = DC(p1); dc2 = DC(p2);
347: deg1 = DEG(dc1); deg2 = DEG(dc2);
348: sgn = cmpq(deg1,deg2);
349: if ( sgn == 0 )
350: if ( !divtpz(vl,COEF(dc1),COEF(dc2),&m) ) {
351: *q = 0;
352: return ( 0 );
353: } else {
354: mulp(vl,p2,m,&m1); subp(vl,p1,m1,&s);
355: if ( !s ) {
356: *q = m;
357: return ( 1 );
358: } else {
359: *q = 0;
360: return ( 0 );
361: }
362: }
363: else if ( sgn < 0 ) {
364: *q = 0;
365: return ( 0 );
366: } else {
367: if ( (PL(NM(deg1)) > 1) ) {
368: error("divtpz : invalid input");
369: *q = 0;
370: return ( 0 );
371: }
372: d1 = QTOS(deg1); d2 = QTOS(deg2);
373: W_CALLOC(d1-d2,P,pq); W_CALLOC(d1,P,pr); W_CALLOC(d2,P,pd);
374: for ( dc = dc1; dc; dc = NEXT(dc) )
375: pr[QTOS(DEG(dc))] = COEF(dc);
376: for ( dc = dc2; dc; dc = NEXT(dc) )
377: pd[QTOS(DEG(dc))] = COEF(dc);
378: for ( dvr = COEF(dc2), i = d1 - d2; i >= 0; i-- )
379: if ( !pr[i+d2] )
380: continue;
381: else if ( !divtpz(vl,pr[i+d2],dvr,&m) ) {
382: *q = 0;
383: return ( 0 );
384: } else {
385: pq[i] = m;
386: for ( j = d2; j >= 0; j-- ) {
387: mulp(vl,pq[i],pd[j],&m);
388: subp(vl,pr[i + j],m,&s); pr[i + j] = s;
389: }
390: }
391: plisttop(pq,v1,d1 - d2,&m); plisttop(pr,v1,d1 - 1,&t);
392: if ( t ) {
393: *q = 0;
394: return ( 0 );
395: } else {
396: *q = m;
397: return ( 1 );
398: }
399: }
400: } else {
401: for ( ; (v1 != vl->v) && (v2 != vl->v); vl = NEXT(vl) );
402: if ( v2 == vl->v ) {
403: *q = 0;
404: return ( 0 );
405: } else
406: return ( divtdcpz(vl,p1,p2,q) ) ;
407: }
408: }
409:
410: int divtdcpz(vl,p1,p2,q)
411: VL vl;
412: P p1,p2,*q;
413: {
414:
415: P m;
416: register DCP dc,dcr,dcr0;
417:
418: for ( dc = DC(p1), dcr0 = 0; dc; dc = NEXT(dc) )
419: if ( !divtpz(vl,COEF(dc),p2,&m) ) {
420: *q = 0;
421: return ( 0 );
422: } else {
423: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = m; NEXT(dcr) = 0;
424: }
425: MKP(VR(p1),dcr0,*q);
426: return ( 1 );
427: }
428:
429: void udivpz(f1,f2,fqp,frp)
430: P f1,f2,*fqp,*frp;
431: {
432: register int n1,n2,i,j;
433: Q *pq,*pr,*pd,d,m,s,q;
434: DCP dc;
435: N qn,rn;
436:
437: if ( !f2 )
438: error("udivpz: division by 0");
439: else if ( !f1 ) {
440: *fqp = *frp = 0;
441: return;
442: } else if ( NUM(f1) )
443: if ( NUM(f2) ) {
444: divn(NM((Q)f1),NM((Q)f2),&qn,&rn);
445: if ( rn ) {
446: *fqp = *frp = 0;
447: } else {
448: NTOQ(qn,SGN((Q)f1)*SGN((Q)f2),q),*fqp = (P)q; *frp = 0;
449: }
450: return;
451: } else {
452: *fqp = 0; *frp = f1;
453: return;
454: }
455: else if ( NUM(f2) ) {
456: n1 = UDEG(f1); W_CALLOC(n1,Q,pq);
457: for ( dc = DC(f1); dc; dc = NEXT(dc) ) {
458: divn(NM((Q)COEF(dc)),NM((Q)f2),&qn,&rn);
459: if ( rn ) {
460: *fqp = *frp = 0;
461: return;
462: } else {
463: NTOQ(qn,SGN((Q)COEF(dc))*SGN((Q)f2),s); pq[QTOS(DEG(dc))] = s;
464: }
465: }
466: plisttop((P *)pq,VR(f1),n1,fqp);
467: return;
468: }
469: n1 = UDEG(f1); n2 = UDEG(f2);
470: if ( n1 < n2 ) {
471: *fqp = NULL; *frp = f1;
472: return;
473: }
474: W_CALLOC(n1-n2,Q,pq); W_CALLOC(n1,Q,pr); W_CALLOC(n2,Q,pd);
475: for ( dc = DC(f1); dc; dc = NEXT(dc) )
476: pr[QTOS(DEG(dc))] = (Q)COEF(dc);
477: for ( dc = DC(f2); dc; dc = NEXT(dc) )
478: pd[QTOS(DEG(dc))] = (Q)COEF(dc);
479: for ( d = (Q)UCOEF(f2), i = n1 - n2; i >= 0; i-- ) {
480: if ( !pr[i+n2] )
481: continue;
482: divn(NM(pr[i+n2]),NM(d),&qn,&rn);
483: if ( rn ) {
484: *fqp = *frp = 0;
485: return;
486: }
487: NTOQ(qn,SGN(pr[i+n2])*SGN(d),pq[i]);
488: for ( j = n2; j >= 0; j-- ) {
489: mulq(pq[i],pd[j],&m); subq(pr[i+j],m,&s); pr[i+j] = s;
490: }
491: }
492: plisttop((P *)pq,VR(f1),n1-n2,fqp); plisttop((P *)pr,VR(f1),n2-1,frp);
493: }
494:
495: void udivpwm(mod,p1,p2,q,r)
496: Q mod;
497: P p1,p2;
498: P *q,*r;
499: {
500: P s,t,u,tq,tr;
501:
502: invl((Q)UCOEF(p2),mod,(Q *)&t); mulpq(p2,t,&s); cmp(mod,s,&u);
503: udivpzwm(mod,p1,u,&tq,&tr);
504: cmp(mod,tr,r); mulpq(tq,t,&s); cmp(mod,s,q);
505: }
506:
507: void udivpzwm(mod,f1,f2,fqp,frp)
508: Q mod;
509: P f1,f2,*fqp,*frp;
510: {
511: register int n1,n2,i,j;
512: Q *pq,*pr,*pd,d,m,s,q;
513: DCP dc;
514: N qn,rn;
515:
516: if ( !f2 )
517: error("udivpz: division by 0");
518: else if ( !f1 ) {
519: *fqp = *frp = 0;
520: return;
521: } else if ( NUM(f1) )
522: if ( NUM(f2) ) {
523: divn(NM((Q)f1),NM((Q)f2),&qn,&rn);
524: if ( rn ) {
525: *fqp = *frp = 0;
526: } else {
527: NTOQ(qn,SGN((Q)f1)*SGN((Q)f2),q),*fqp = (P)q; *frp = 0;
528: }
529: return;
530: } else {
531: *fqp = 0; *frp = f1;
532: return;
533: }
534: else if ( NUM(f2) ) {
535: n1 = UDEG(f1); W_CALLOC(n1,Q,pq);
536: for ( dc = DC(f1); dc; dc = NEXT(dc) ) {
537: divn(NM((Q)COEF(dc)),NM((Q)f2),&qn,&rn);
538: if ( rn ) {
539: *fqp = *frp = 0;
540: return;
541: } else {
542: NTOQ(qn,SGN((Q)COEF(dc))*SGN((Q)f2),s); pq[QTOS(DEG(dc))] = s;
543: }
544: }
545: plisttop((P *)pq,VR(f1),n1,fqp);
546: return;
547: }
548: n1 = UDEG(f1); n2 = UDEG(f2);
549: if ( n1 < n2 ) {
550: *fqp = NULL; *frp = f1;
551: return;
552: }
553: W_CALLOC(n1-n2,Q,pq); W_CALLOC(n1,Q,pr); W_CALLOC(n2,Q,pd);
554: for ( dc = DC(f1); dc; dc = NEXT(dc) )
555: pr[QTOS(DEG(dc))] = (Q)COEF(dc);
556: for ( dc = DC(f2); dc; dc = NEXT(dc) )
557: pd[QTOS(DEG(dc))] = (Q)COEF(dc);
558: for ( d = (Q)UCOEF(f2), i = n1 - n2; i >= 0; i-- ) {
559: if ( !pr[i+n2] )
560: continue;
561: divn(NM(pr[i+n2]),NM(d),&qn,&rn);
562: if ( rn ) {
563: *fqp = *frp = 0;
564: return;
565: }
566: NTOQ(qn,SGN(pr[i+n2])*SGN(d),pq[i]);
567: for ( j = n2; j >= 0; j-- ) {
568: mulq(pq[i],pd[j],&m); remq(m,mod,&s);
569: subq(pr[i+j],s,&m); remq(m,mod,&s); pr[i+j] = s;
570: }
571: }
572: plisttop((P *)pq,VR(f1),n1-n2,fqp); plisttop((P *)pr,VR(f1),n2-1,frp);
573: }
574: #endif