![[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.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.11 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/P.c,v 1.10 2004/08/18 00:17:02 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:
1.10 noro 57: #include "poly.c"
1.1 noro 58:
59: void ptozp(p,sgn,c,pr)
60: P p;
61: int sgn;
62: Q *c;
63: P *pr;
64: {
65: N nm,dn;
66:
67: if ( !p ) {
68: *c = 0; *pr = 0;
69: } else {
70: lgp(p,&nm,&dn);
71: if ( UNIN(dn) )
72: NTOQ(nm,sgn,*c);
73: else
74: NDTOQ(nm,dn,sgn,*c);
75: divcp(p,*c,pr);
76: }
77: }
78:
79: void lgp(p,g,l)
80: P p;
81: N *g,*l;
82: {
83: N g1,g2,l1,l2,l3,l4,tmp;
84: DCP dc;
85:
86: if ( NUM(p) ) {
87: *g = NM((Q)p);
88: if ( INT((Q)p) )
89: *l = ONEN;
90: else
91: *l = DN((Q)p);
92: } else {
93: dc = DC(p); lgp(COEF(dc),g,l);
94: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
95: lgp(COEF(dc),&g1,&l1); gcdn(*g,g1,&g2); *g = g2;
96: gcdn(*l,l1,&l2); kmuln(*l,l1,&l3); divn(l3,l2,&l4,&tmp); *l = l4;
97: }
98: }
99: }
100:
101: void divcp(f,q,rp)
102: P f;
103: Q q;
104: P *rp;
105: {
106: DCP dc,dcr,dcr0;
107:
108: if ( !f )
109: *rp = 0;
110: else if ( NUM(f) )
111: DIVNUM(f,q,rp);
112: else {
113: for ( dc = DC(f), dcr0 = 0; dc; dc = NEXT(dc) ) {
114: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
115: divcp(COEF(dc),q,&COEF(dcr));
116: }
117: NEXT(dcr) = 0; MKP(VR(f),dcr0,*rp);
118: }
119: }
120:
121: void diffp(vl,p,v,r)
122: VL vl;
123: P p;
124: V v;
125: P *r;
126: {
127: P t;
128: DCP dc,dcr,dcr0;
129:
130: if ( !p || NUM(p) )
131: *r = 0;
132: else {
133: if ( v == VR(p) ) {
1.11 ! noro 134: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 135: if ( !DEG(dc) ) continue;
1.1 noro 136: MULPQ(COEF(dc),(P)DEG(dc),&t);
137: if ( t ) {
138: NEXTDC(dcr0,dcr); SUBQ(DEG(dc),ONE,&DEG(dcr));
139: COEF(dcr) = t;
140: }
141: }
142: if ( !dcr0 )
143: *r = 0;
144: else {
145: NEXT(dcr) = 0; MKP(v,dcr0,*r);
146: }
147: } else {
148: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
149: diffp(vl,COEF(dc),v,&t);
1.8 ohara 150: if ( t ) {
151: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
152: }
153: }
154: if ( !dcr0 )
155: *r = 0;
156: else {
157: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
158: }
159: }
160: }
161: }
162:
163: /* Euler derivation */
164: void ediffp(vl,p,v,r)
165: VL vl;
166: P p;
167: V v;
168: P *r;
169: {
170: P t;
171: DCP dc,dcr,dcr0;
172:
173: if ( !p || NUM(p) )
174: *r = 0;
175: else {
176: if ( v == VR(p) ) {
1.11 ! noro 177: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 178: if ( !DEG(dc) ) continue;
1.8 ohara 179: MULPQ(COEF(dc),(P)DEG(dc),&t);
180: if ( t ) {
181: NEXTDC(dcr0,dcr);
182: DEG(dcr) = DEG(dc);
183: COEF(dcr) = t;
184: }
185: }
186: if ( !dcr0 )
187: *r = 0;
188: else {
189: NEXT(dcr) = 0; MKP(v,dcr0,*r);
190: }
191: } else {
192: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
193: ediffp(vl,COEF(dc),v,&t);
1.1 noro 194: if ( t ) {
195: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
196: }
197: }
198: if ( !dcr0 )
199: *r = 0;
200: else {
201: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
202: }
203: }
204: }
205: }
206:
207: void coefp(p,d,pr)
208: P p;
209: int d;
210: P *pr;
211: {
212: DCP dc;
213: int sgn;
214: Q dq;
215:
216: if ( NUM(p) )
217: if ( d == 0 )
218: *pr = p;
219: else
220: *pr = 0;
221: else {
222: for ( STOQ(d,dq), dc = DC(p); dc; dc = NEXT(dc) )
223: if ( (sgn = cmpq(DEG(dc),dq)) > 0 )
224: continue;
225: else if ( sgn == 0 ) {
226: *pr = COEF(dc);
227: return;
228: } else {
229: *pr = 0;
230: break;
231: }
232: *pr = 0;
233: }
234: }
235:
236: int compp(vl,p1,p2)
237: VL vl;
238: P p1,p2;
239: {
240: DCP dc1,dc2;
241: V v1,v2;
242:
243: if ( !p1 )
244: return p2 ? -1 : 0;
245: else if ( !p2 )
246: return 1;
247: else if ( NUM(p1) )
248: return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;
249: else if ( NUM(p2) )
250: return 1;
251: if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {
252: for ( dc1 = DC(p1), dc2 = DC(p2);
253: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )
254: switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
255: case 1:
256: return 1; break;
257: case -1:
258: return -1; break;
259: default:
260: switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {
261: case 1:
262: return 1; break;
263: case -1:
264: return -1; break;
265: default:
266: break;
267: }
268: break;
269: }
270: return dc1 ? 1 : (dc2 ? -1 : 0);
271: } else {
272: for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );
273: return v1 == VR(vl) ? 1 : -1;
274: }
1.9 noro 275: }
276:
277: int equalp(vl,p1,p2)
278: VL vl;
279: P p1,p2;
280: {
281: DCP dc1,dc2;
282: V v1,v2;
283:
284: if ( !p1 ) {
285: if ( !p2 ) return 1;
286: else return 0;
287: }
288: /* p1 != 0 */
289: if ( !p2 ) return 0;
290:
291: /* p1 != 0, p2 != 0 */
292: if ( NUM(p1) ) {
293: if ( !NUM(p2) ) return 0;
294: /* p1 and p2 are numbers */
295: if ( NID((Num)p1) != NID((Num)p2) ) return 0;
296: if ( !(*cmpnumt[NID(p1),NID(p1)])(p1,p2) ) return 1;
297: return 0;
298: }
299: if ( VR(p1) != VR(p2) ) return 0;
300: for ( dc1 = DC(p1), dc2 = DC(p2);
301: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) ) {
302: if ( cmpq(DEG(dc1),DEG(dc2)) ) return 0;
303: else if ( !equalp(vl,COEF(dc1),COEF(dc2)) ) return 0;
304: }
305: if ( dc1 || dc2 ) return 0;
306: else return 1;
1.1 noro 307: }
308:
309: void csump(vl,p,c)
310: VL vl;
311: P p;
312: Q *c;
313: {
314: DCP dc;
315: Q s,s1,s2;
316:
317: if ( !p || NUM(p) )
318: *c = (Q)p;
319: else {
320: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
321: csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
322: }
323: *c = s;
324: }
325: }
326:
327: void degp(v,p,d)
328: V v;
329: P p;
330: Q *d;
331: {
332: DCP dc;
333: Q m,m1;
334:
335: if ( !p || NUM(p) )
336: *d = 0;
337: else if ( v == VR(p) )
338: *d = DEG(DC(p));
339: else {
340: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
341: degp(v,COEF(dc),&m1);
342: if ( cmpq(m,m1) < 0 )
343: m = m1;
344: }
345: *d = m;
1.6 noro 346: }
347: }
348:
349: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr);
350: void mulpq_trunc(P p,Q q,VN vn,P *pr);
351: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr);
352:
353: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr)
354: {
355: register DCP dc,dct,dcr,dcr0;
356: V v1,v2;
357: P t,s,u;
358: int n1,n2,i,d,d1;
359:
360: if ( !p1 || !p2 ) *pr = 0;
361: else if ( NUM(p1) )
362: mulpq_trunc(p2,(Q)p1,vn,pr);
363: else if ( NUM(p2) )
364: mulpq_trunc(p1,(Q)p2,vn,pr);
365: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
366: for ( ; vn->v && vn->v != v1; vn++ )
367: if ( vn->n ) {
368: /* p1,p2 do not contain vn->v */
369: *pr = 0;
370: return;
371: }
372: if ( !vn->v )
373: error("mulp_trunc : invalid vn");
374: d = vn->n;
375: for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {
376: for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {
377: d1 = QTOS(DEG(dct))+QTOS(DEG(dc));
378: if ( d1 >= d ) {
379: mulp_trunc(vl,COEF(dct),COEF(dc),vn+1,&t);
380: if ( t ) {
381: NEXTDC(dcr0,dcr);
382: STOQ(d1,DEG(dcr));
383: COEF(dcr) = t;
384: }
385: }
386: }
387: if ( dcr0 ) {
388: NEXT(dcr) = 0; MKP(v1,dcr0,t);
389: addp(vl,s,t,&u); s = u; t = u = 0;
390: }
391: }
392: *pr = s;
393: } else {
394: while ( v1 != VR(vl) && v2 != VR(vl) )
395: vl = NEXT(vl);
396: if ( v1 == VR(vl) )
397: mulpc_trunc(vl,p1,p2,vn,pr);
398: else
399: mulpc_trunc(vl,p2,p1,vn,pr);
400: }
401: }
402:
403: void mulpq_trunc(P p,Q q,VN vn,P *pr)
404: {
405: DCP dc,dcr,dcr0;
406: P t;
407: int i,d;
408: V v;
409:
410: if (!p || !q)
411: *pr = 0;
412: else if ( NUM(p) ) {
413: for ( ; vn->v; vn++ )
414: if ( vn->n ) {
415: *pr = 0;
416: return;
417: }
418: MULNUM(p,q,pr);
419: } else {
420: v = VR(p);
421: for ( ; vn->v && vn->v != v; vn++ ) {
422: if ( vn->n ) {
423: /* p does not contain vn->v */
424: *pr = 0;
425: return;
426: }
427: }
428: if ( !vn->v )
429: error("mulpq_trunc : invalid vn");
430: d = vn->n;
431: for ( dcr0 = 0, dc = DC(p); dc && QTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
432: mulpq_trunc(COEF(dc),q,vn+1,&t);
433: if ( t ) {
434: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
435: }
436: }
437: if ( dcr0 ) {
438: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
439: } else
440: *pr = 0;
441: }
442: }
443:
444: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr)
445: {
446: DCP dc,dcr,dcr0;
447: P t;
448: V v;
449: int i,d;
450:
451: if ( NUM(c) )
452: mulpq_trunc(p,(Q)c,vn,pr);
453: else {
454: v = VR(p);
455: for ( ; vn->v && vn->v != v; vn++ )
456: if ( vn->n ) {
457: /* p,c do not contain vn->v */
458: *pr = 0;
459: return;
460: }
461: if ( !vn->v )
462: error("mulpc_trunc : invalid vn");
463: d = vn->n;
464: for ( dcr0 = 0, dc = DC(p); dc && QTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
465: mulp_trunc(vl,COEF(dc),c,vn+1,&t);
466: if ( t ) {
467: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
468: }
469: }
470: if ( dcr0 ) {
471: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
472: } else
473: *pr = 0;
1.7 noro 474: }
475: }
476:
477: void quop_trunc(VL vl,P p1,P p2,VN vn,P *pr)
478: {
479: DCP dc,dcq0,dcq;
480: P t,s,m,lc2,qt;
481: V v1,v2;
482: Q d2;
483: VN vn1;
484:
485: if ( !p1 )
486: *pr = 0;
487: else if ( NUM(p2) )
488: divsp(vl,p1,p2,pr);
489: else if ( (v1 = VR(p1)) != (v2 = VR(p2)) ) {
490: for ( dcq0 = 0, dc = DC(p1); dc; dc = NEXT(dc) ) {
491: NEXTDC(dcq0,dcq);
492: DEG(dcq) = DEG(dc);
493: quop_trunc(vl,COEF(dc),p2,vn,&COEF(dcq));
494: }
495: NEXT(dcq) = 0;
496: MKP(v1,dcq0,*pr);
497: } else {
498: d2 = DEG(DC(p2));
499: lc2 = COEF(DC(p2));
500: t = p1;
501: dcq0 = 0;
502: /* vn1 = degree list of LC(p2) */
503: for ( vn1 = vn; vn1->v != v1; vn1++ );
504: vn1++;
505: while ( t ) {
506: dc = DC(t);
507: NEXTDC(dcq0,dcq);
508: subq(DEG(dc),d2,&DEG(dcq));
509: quop_trunc(vl,COEF(dc),lc2,vn1,&COEF(dcq));
510: NEXT(dcq) = 0;
511: MKP(v1,dcq,qt);
512: mulp_trunc(vl,p2,qt,vn,&m);
513: subp(vl,t,m,&s); t = s;
514: }
515: NEXT(dcq) = 0;
516: MKP(v1,dcq0,*pr);
1.1 noro 517: }
518: }