Annotation of OpenXM_contrib2/asir2018/engine/P.c, Revision 1.2
1.1 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
26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
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.2 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2018/engine/P.c,v 1.1 2018/09/19 05:45:07 noro Exp $
1.1 noro 49: */
50: #ifndef FBASE
51: #define FBASE
52: #endif
53:
54: #include "b.h"
55: #include "ca.h"
56:
57: #include "poly.c"
58:
59: void divcp(P f,Q q,P *rp)
60: {
61: DCP dc,dcr,dcr0;
62:
63: if ( !f )
64: *rp = 0;
65: else if ( NUM(f) )
66: DIVNUM(f,q,rp);
67: else {
68: for ( dc = DC(f), dcr0 = 0; dc; dc = NEXT(dc) ) {
69: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
70: divcp(COEF(dc),q,&COEF(dcr));
71: }
72: NEXT(dcr) = 0; MKP(VR(f),dcr0,*rp);
73: }
74: }
75:
76: void diffp(VL vl,P p,V v,P *r)
77: {
78: P t;
79: DCP dc,dcr,dcr0;
80:
81: if ( !p || NUM(p) )
82: *r = 0;
83: else {
84: if ( v == VR(p) ) {
85: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
86: if ( !DEG(dc) ) continue;
87: MULPQ(COEF(dc),(P)DEG(dc),&t);
88: if ( t ) {
89: NEXTDC(dcr0,dcr); subz(DEG(dc),ONE,&DEG(dcr));
90: COEF(dcr) = t;
91: }
92: }
93: if ( !dcr0 )
94: *r = 0;
95: else {
96: NEXT(dcr) = 0; MKP(v,dcr0,*r);
97: }
98: } else {
99: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
100: diffp(vl,COEF(dc),v,&t);
101: if ( t ) {
102: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
103: }
104: }
105: if ( !dcr0 )
106: *r = 0;
107: else {
108: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
109: }
110: }
111: }
112: }
113:
114: /* Euler derivation */
115: void ediffp(VL vl,P p,V v,P *r)
116: {
117: P t;
118: DCP dc,dcr,dcr0;
119:
120: if ( !p || NUM(p) )
121: *r = 0;
122: else {
123: if ( v == VR(p) ) {
124: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
125: if ( !DEG(dc) ) continue;
126: MULPQ(COEF(dc),(P)DEG(dc),&t);
127: if ( t ) {
128: NEXTDC(dcr0,dcr);
129: DEG(dcr) = DEG(dc);
130: COEF(dcr) = t;
131: }
132: }
133: if ( !dcr0 )
134: *r = 0;
135: else {
136: NEXT(dcr) = 0; MKP(v,dcr0,*r);
137: }
138: } else {
139: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
140: ediffp(vl,COEF(dc),v,&t);
141: if ( t ) {
142: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
143: }
144: }
145: if ( !dcr0 )
146: *r = 0;
147: else {
148: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
149: }
150: }
151: }
152: }
153:
154: void coefp(P p,int d,P *pr)
155: {
156: DCP dc;
157: int sgn;
158: Z dq;
159:
160: if ( NUM(p) )
161: if ( d == 0 )
162: *pr = p;
163: else
164: *pr = 0;
165: else {
1.2 ! noro 166: for ( STOZ(d,dq), dc = DC(p); dc; dc = NEXT(dc) )
1.1 noro 167: if ( (sgn = cmpz(DEG(dc),dq)) > 0 )
168: continue;
169: else if ( sgn == 0 ) {
170: *pr = COEF(dc);
171: return;
172: } else {
173: *pr = 0;
174: break;
175: }
176: *pr = 0;
177: }
178: }
179:
180: int compp(VL vl,P p1,P p2)
181: {
182: DCP dc1,dc2;
183: V v1,v2;
184:
185: if ( !p1 )
186: return p2 ? -1 : 0;
187: else if ( !p2 )
188: return 1;
189: else if ( NUM(p1) )
190: return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;
191: else if ( NUM(p2) )
192: return 1;
193: if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {
194: for ( dc1 = DC(p1), dc2 = DC(p2);
195: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )
196: switch ( cmpz(DEG(dc1),DEG(dc2)) ) {
197: case 1:
198: return 1; break;
199: case -1:
200: return -1; break;
201: default:
202: switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {
203: case 1:
204: return 1; break;
205: case -1:
206: return -1; break;
207: default:
208: break;
209: }
210: break;
211: }
212: return dc1 ? 1 : (dc2 ? -1 : 0);
213: } else {
214: for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );
215: return v1 == VR(vl) ? 1 : -1;
216: }
217: }
218:
219: int equalp(VL vl,P p1,P p2)
220: {
221: DCP dc1,dc2;
222: V v1,v2;
223:
224: if ( !p1 ) {
225: if ( !p2 ) return 1;
226: else return 0;
227: }
228: /* p1 != 0 */
229: if ( !p2 ) return 0;
230:
231: /* p1 != 0, p2 != 0 */
232: if ( NUM(p1) ) {
233: if ( !NUM(p2) ) return 0;
234: /* p1 and p2 are numbers */
235: if ( NID((Num)p1) != NID((Num)p2) ) return 0;
236: if ( !(*cmpnumt[NID(p1),NID(p1)])(p1,p2) ) return 1;
237: return 0;
238: }
239: if ( VR(p1) != VR(p2) ) return 0;
240: for ( dc1 = DC(p1), dc2 = DC(p2);
241: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) ) {
242: if ( cmpz(DEG(dc1),DEG(dc2)) ) return 0;
243: else if ( !equalp(vl,COEF(dc1),COEF(dc2)) ) return 0;
244: }
245: if ( dc1 || dc2 ) return 0;
246: else return 1;
247: }
248:
249: void csump(VL vl,P p,Q *c)
250: {
251: DCP dc;
252: Q s,s1,s2;
253:
254: if ( !p || NUM(p) )
255: *c = (Q)p;
256: else {
257: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
258: csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
259: }
260: *c = s;
261: }
262: }
263:
264: void degp(V v,P p,Z *d)
265: {
266: DCP dc;
267: Z m,m1;
268:
269: if ( !p || NUM(p) )
270: *d = 0;
271: else if ( v == VR(p) )
272: *d = DEG(DC(p));
273: else {
274: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
275: degp(v,COEF(dc),&m1);
276: if ( cmpz(m,m1) < 0 )
277: m = m1;
278: }
279: *d = m;
280: }
281: }
282:
283: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr);
284: void mulpq_trunc(P p,Q q,VN vn,P *pr);
285: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr);
286:
287: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr)
288: {
289: register DCP dc,dct,dcr,dcr0;
290: V v1,v2;
291: P t,s,u;
292: int n1,n2,i,d,d1;
293:
294: if ( !p1 || !p2 ) *pr = 0;
295: else if ( NUM(p1) )
296: mulpq_trunc(p2,(Q)p1,vn,pr);
297: else if ( NUM(p2) )
298: mulpq_trunc(p1,(Q)p2,vn,pr);
299: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
300: for ( ; vn->v && vn->v != v1; vn++ )
301: if ( vn->n ) {
302: /* p1,p2 do not contain vn->v */
303: *pr = 0;
304: return;
305: }
306: if ( !vn->v )
307: error("mulp_trunc : invalid vn");
308: d = vn->n;
309: for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {
310: for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {
1.2 ! noro 311: d1 = ZTOS(DEG(dct))+ZTOS(DEG(dc));
1.1 noro 312: if ( d1 >= d ) {
313: mulp_trunc(vl,COEF(dct),COEF(dc),vn+1,&t);
314: if ( t ) {
315: NEXTDC(dcr0,dcr);
1.2 ! noro 316: STOZ(d1,DEG(dcr));
1.1 noro 317: COEF(dcr) = t;
318: }
319: }
320: }
321: if ( dcr0 ) {
322: NEXT(dcr) = 0; MKP(v1,dcr0,t);
323: addp(vl,s,t,&u); s = u; t = u = 0;
324: }
325: }
326: *pr = s;
327: } else {
328: while ( v1 != VR(vl) && v2 != VR(vl) )
329: vl = NEXT(vl);
330: if ( v1 == VR(vl) )
331: mulpc_trunc(vl,p1,p2,vn,pr);
332: else
333: mulpc_trunc(vl,p2,p1,vn,pr);
334: }
335: }
336:
337: void mulpq_trunc(P p,Q q,VN vn,P *pr)
338: {
339: DCP dc,dcr,dcr0;
340: P t;
341: int i,d;
342: V v;
343:
344: if (!p || !q)
345: *pr = 0;
346: else if ( NUM(p) ) {
347: for ( ; vn->v; vn++ )
348: if ( vn->n ) {
349: *pr = 0;
350: return;
351: }
352: MULNUM(p,q,pr);
353: } else {
354: v = VR(p);
355: for ( ; vn->v && vn->v != v; vn++ ) {
356: if ( vn->n ) {
357: /* p does not contain vn->v */
358: *pr = 0;
359: return;
360: }
361: }
362: if ( !vn->v )
363: error("mulpq_trunc : invalid vn");
364: d = vn->n;
1.2 ! noro 365: for ( dcr0 = 0, dc = DC(p); dc && ZTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
1.1 noro 366: mulpq_trunc(COEF(dc),q,vn+1,&t);
367: if ( t ) {
368: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
369: }
370: }
371: if ( dcr0 ) {
372: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
373: } else
374: *pr = 0;
375: }
376: }
377:
378: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr)
379: {
380: DCP dc,dcr,dcr0;
381: P t;
382: V v;
383: int i,d;
384:
385: if ( NUM(c) )
386: mulpq_trunc(p,(Q)c,vn,pr);
387: else {
388: v = VR(p);
389: for ( ; vn->v && vn->v != v; vn++ )
390: if ( vn->n ) {
391: /* p,c do not contain vn->v */
392: *pr = 0;
393: return;
394: }
395: if ( !vn->v )
396: error("mulpc_trunc : invalid vn");
397: d = vn->n;
1.2 ! noro 398: for ( dcr0 = 0, dc = DC(p); dc && ZTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
1.1 noro 399: mulp_trunc(vl,COEF(dc),c,vn+1,&t);
400: if ( t ) {
401: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
402: }
403: }
404: if ( dcr0 ) {
405: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
406: } else
407: *pr = 0;
408: }
409: }
410:
411: void quop_trunc(VL vl,P p1,P p2,VN vn,P *pr)
412: {
413: DCP dc,dcq0,dcq;
414: P t,s,m,lc2,qt;
415: V v1,v2;
416: Z d2;
417: VN vn1;
418:
419: if ( !p1 )
420: *pr = 0;
421: else if ( NUM(p2) )
422: divsp(vl,p1,p2,pr);
423: else if ( (v1 = VR(p1)) != (v2 = VR(p2)) ) {
424: for ( dcq0 = 0, dc = DC(p1); dc; dc = NEXT(dc) ) {
425: NEXTDC(dcq0,dcq);
426: DEG(dcq) = DEG(dc);
427: quop_trunc(vl,COEF(dc),p2,vn,&COEF(dcq));
428: }
429: NEXT(dcq) = 0;
430: MKP(v1,dcq0,*pr);
431: } else {
432: d2 = DEG(DC(p2));
433: lc2 = COEF(DC(p2));
434: t = p1;
435: dcq0 = 0;
436: /* vn1 = degree list of LC(p2) */
437: for ( vn1 = vn; vn1->v != v1; vn1++ );
438: vn1++;
439: while ( t ) {
440: dc = DC(t);
441: NEXTDC(dcq0,dcq);
442: subz(DEG(dc),d2,&DEG(dcq));
443: quop_trunc(vl,COEF(dc),lc2,vn1,&COEF(dcq));
444: NEXT(dcq) = 0;
445: MKP(v1,dcq,qt);
446: mulp_trunc(vl,p2,qt,vn,&m);
447: subp(vl,t,m,&s); t = s;
448: }
449: NEXT(dcq) = 0;
450: MKP(v1,dcq0,*pr);
451: }
452: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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 / asir2018 / engine