Annotation of OpenXM_contrib2/asir2018/engine/poly.c, Revision 1.2
1.1 noro 1: #ifdef MODULAR
2: void addmp(VL vl,int mod,P p1,P p2,P *pr)
3: #else
4: void addp(VL vl,P p1,P p2,P *pr)
5: #endif
6: {
7: DCP dc1,dc2,dcr0,dcr;
8: V v1,v2;
9: P t;
10:
11: if ( !p1 )
12: *pr = p2;
13: else if ( !p2 )
14: *pr = p1;
15: else if ( NUM(p1) )
16: if ( NUM(p2) )
17: ADDNUM(p1,p2,pr);
18: else
19: ADDPQ(p2,p1,pr);
20: else if ( NUM(p2) )
21: ADDPQ(p1,p2,pr);
22: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
23: for ( dc1 = DC(p1), dc2 = DC(p2), dcr0 = 0; dc1 && dc2; )
24: switch ( cmpz(DEG(dc1),DEG(dc2)) ) {
25: case 0:
26: ADDP(vl,COEF(dc1),COEF(dc2),&t);
27: if ( t ) {
28: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc1); COEF(dcr) = t;
29: }
30: dc1 = NEXT(dc1); dc2 = NEXT(dc2); break;
31: case 1:
32: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc1); COEF(dcr) = COEF(dc1);
33: dc1 = NEXT(dc1); break;
34: case -1:
35: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc2); COEF(dcr) = COEF(dc2);
36: dc2 = NEXT(dc2); break;
37: }
38: if ( !dcr0 )
39: if ( dc1 )
40: dcr0 = dc1;
41: else if ( dc2 )
42: dcr0 = dc2;
43: else {
44: *pr = 0;
45: return;
46: }
47: else
48: if ( dc1 )
49: NEXT(dcr) = dc1;
50: else if ( dc2 )
51: NEXT(dcr) = dc2;
52: else
53: NEXT(dcr) = 0;
54: MKP(v1,dcr0,*pr);
55: } else {
56: while ( v1 != VR(vl) && v2 != VR(vl) )
57: vl = NEXT(vl);
58: if ( v1 == VR(vl) )
59: ADDPTOC(vl,p1,p2,pr);
60: else
61: ADDPTOC(vl,p2,p1,pr);
62: }
63: }
64:
65: #ifdef MODULAR
66: void addmpq(int mod,P p,P q,P *pr)
67: #else
68: void addpq(P p,P q,P *pr)
69: #endif
70: {
71: DCP dc,dcr,dcr0;
72: P t;
73:
74: if ( NUM(p) )
75: ADDNUM(p,q,pr);
76: else {
77: dc = DC(p);
78: for ( dcr0 = 0; dc && cmpz(DEG(dc),0) > 0; dc = NEXT(dc) ) {
79: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = COEF(dc);
80: }
81: if ( !dc ) {
82: NEXTDC(dcr0,dcr);
83: DEG(dcr) = 0; COEF(dcr) = q; NEXT(dcr) = 0;
84: } else if ( !DEG(dc) ) {
85: ADDPQ(COEF(dc),q,&t);
86: if ( t ) {
87: NEXTDC(dcr0,dcr);
88: DEG(dcr) = 0; COEF(dcr) = t;
89: }
90: NEXT(dcr) = NEXT(dc);
91: } else {
92: NEXTDC(dcr0,dcr);
93: DEG(dcr) = 0; COEF(dcr) = q; NEXT(dcr) = dc;
94: }
95: MKP(VR(p),dcr0,*pr);
96: }
97: }
98:
99: #ifdef MODULAR
100: void addmptoc(VL vl,int mod,P p,P c,P *pr)
101: #else
102: void addptoc(VL vl,P p,P c,P *pr)
103: #endif
104: {
105: DCP dc,dcr,dcr0;
106: P t;
107:
108: for ( dcr0 = 0, dc = DC(p); dc && cmpz(DEG(dc),0) > 0; dc = NEXT(dc) ) {
109: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = COEF(dc);
110: }
111: if ( !dc ) {
112: NEXTDC(dcr0,dcr);
113: DEG(dcr) = 0; COEF(dcr) = c; NEXT(dcr) = 0;
114: } else if ( !DEG(dc) ) {
115: ADDP(vl,COEF(dc),c,&t);
116: if ( t ) {
117: NEXTDC(dcr0,dcr);
118: DEG(dcr) = 0; COEF(dcr) = t;
119: }
120: NEXT(dcr) = NEXT(dc);
121: } else {
122: NEXTDC(dcr0,dcr);
123: DEG(dcr) = 0; COEF(dcr) = c; NEXT(dcr) = dc;
124: }
125: MKP(VR(p),dcr0,*pr);
126: }
127:
128: #ifdef MODULAR
129: void submp(VL vl,int mod,P p1,P p2,P *pr)
130: #else
131: void subp(VL vl,P p1,P p2,P *pr)
132: #endif
133: {
134: P t;
135:
136: if ( !p2 )
137: *pr = p1;
138: else {
139: CHSGNP(p2,&t); ADDP(vl,p1,t,pr);
140: }
141: }
142:
143: #ifdef MODULAR
144: void mulmp(VL vl,int mod,P p1,P p2,P *pr)
145: #else
146: void mulp(VL vl,P p1,P p2,P *pr)
147: #endif
148: {
149: DCP dc,dct,dcr,dcr0;
150: V v1,v2;
151: P t,s,u;
152: int n1,n2;
153:
154: if ( !p1 || !p2 ) *pr = 0;
155: else if ( NUM(p1) )
156: MULPQ(p2,p1,pr);
157: else if ( NUM(p2) )
158: MULPQ(p1,p2,pr);
159: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
160: for ( dc = DC(p1), n1 = 0; dc; dc = NEXT(dc), n1++ );
161: for ( dc = DC(p2), n2 = 0; dc; dc = NEXT(dc), n2++ );
162: if ( n1 > n2 )
163: for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {
164: for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {
165: NEXTDC(dcr0,dcr); MULP(vl,COEF(dct),COEF(dc),&COEF(dcr));
166: addz(DEG(dct),DEG(dc),&DEG(dcr));
167: }
168: NEXT(dcr) = 0; MKP(v1,dcr0,t);
169: ADDP(vl,s,t,&u); s = u; t = u = 0;
170: }
171: else
172: for ( dc = DC(p1), s = 0; dc; dc = NEXT(dc) ) {
173: for ( dcr0 = 0, dct = DC(p2); dct; dct = NEXT(dct) ) {
174: NEXTDC(dcr0,dcr); MULP(vl,COEF(dct),COEF(dc),&COEF(dcr));
175: addz(DEG(dct),DEG(dc),&DEG(dcr));
176: }
177: NEXT(dcr) = 0; MKP(v1,dcr0,t);
178: ADDP(vl,s,t,&u); s = u; t = u = 0;
179: }
180: *pr = s;
181: } else {
182: while ( v1 != VR(vl) && v2 != VR(vl) )
183: vl = NEXT(vl);
184: if ( v1 == VR(vl) )
185: MULPC(vl,p1,p2,pr);
186: else
187: MULPC(vl,p2,p1,pr);
188: }
189: }
190:
191: #ifdef MODULAR
192: void mulmpq(int mod,P p,P q,P *pr)
193: #else
194: void mulpq(P p,P q,P *pr)
195: #endif
196: {
197: DCP dc,dcr,dcr0;
198: P t;
199:
200: if (!p || !q)
201: *pr = 0;
202: else if ( Uniq(q) )
203: *pr = p;
204: else if ( NUM(p) )
205: MULNUM(p,q,pr);
206: else {
207: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
208: MULPQ(COEF(dc),q,&t);
209: if ( t ) {
210: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
211: }
212: }
213: if ( dcr0 ) {
214: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
215: } else
216: *pr = 0;
217: }
218: }
219:
220: #ifdef MODULAR
221: void mulmpc(VL vl,int mod,P p,P c,P *pr)
222: #else
223: void mulpc(VL vl,P p,P c,P *pr)
224: #endif
225: {
226: DCP dc,dcr,dcr0;
227: P t;
228:
229: if ( NUM(c) )
230: MULPQ(p,c,pr);
231: else {
232: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
233: MULP(vl,COEF(dc),c,&t);
234: if ( t ) {
235: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
236: }
237: }
238: if ( dcr0 ) {
239: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
240: } else
241: *pr = 0;
242: }
243: }
244:
245: #ifdef MODULAR
246: void pwrmp(VL vl,int mod,P p,Z q,P *pr)
247: #else
248: void pwrp(VL vl,P p,Z q,P *pr)
249: #endif
250: {
251: DCP dc,dcr;
252: int n,i;
253: P *x,*y;
254: P t,s,u;
255: DCP dct;
256: P *pt;
257:
258: if ( !q ) {
259: *pr = (P)One;
260: } else if ( !p )
261: *pr = 0;
262: else if ( UNIQ(q) )
263: *pr = p;
264: else if ( NUM(p) )
265: PWRNUM(p,(Q)q,pr);
266: else {
267: dc = DC(p);
268: if ( !NEXT(dc) ) {
269: NEWDC(dcr);
270: PWRP(vl,COEF(dc),q,&COEF(dcr)); mulz(DEG(dc),(Z)q,&DEG(dcr));
271: NEXT(dcr) = 0; MKP(VR(p),dcr,*pr);
272: } else if ( !INT(q) ) {
273: error("pwrp: can't calculate fractional power."); *pr = 0;
274: } else if ( smallz(q) ) {
1.2 ! noro 275: n = ZTOS(q); x = (P *)ALLOCA((n+1)*sizeof(pointer));
1.1 noro 276: NEWDC(dct); DEG(dct) = DEG(dc); COEF(dct) = COEF(dc);
277: NEXT(dct) = 0; MKP(VR(p),dct,t);
278: for ( i = 0, u = (P)One; i < n; i++ ) {
279: x[i] = u; MULP(vl,u,t,&s); u = s;
280: }
281: x[n] = u; y = (P *)ALLOCA((n+1)*sizeof(pointer));
282:
283: MKP(VR(p),NEXT(dc),t);
284: for ( i = 0, u = (P)One; i < n; i++ ) {
285: y[i] = u; MULP(vl,u,t,&s); u = s;
286: }
287: y[n] = u;
288: pt = (P *)ALLOCA((n+1)*sizeof(pointer)); MKBC(n,pt);
289: for ( i = 0, u = 0; i <= n; i++ ) {
290: MULP(vl,x[i],y[n-i],&t); MULP(vl,t,pt[i],&s);
291: ADDP(vl,u,s,&t); u = t;
292: }
293: *pr = u;
294: } else {
295: error("exponent too big");
296: *pr = 0;
297: }
298: }
299: }
300:
301: #ifdef MODULAR
302: void chsgnmp(int mod,P p,P *pr)
303: #else
304: void chsgnp(P p,P *pr)
305: #endif
306: {
307: register DCP dc,dcr,dcr0;
308:
309: if ( !p )
310: *pr = NULL;
311: else if ( NUM(p) ) {
312: #if defined(_PA_RISC1_1) || defined(__alpha) || defined(mips) || defined(_IBMR2)
313: #ifdef FBASE
314: chsgnnum((Num)p,(Num *)pr);
315: #else
316: MQ mq;
317:
318: NEWMQ(mq); CONT(mq)=mod-CONT((MQ)p); *pr = (P)mq;
319: #endif
320: #else
321: CHSGNNUM(p,*pr);
322: #endif
323: } else {
324: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
325: NEXTDC(dcr0,dcr); CHSGNP(COEF(dc),&COEF(dcr)); DEG(dcr) = DEG(dc);
326: }
327: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
328: }
329: }
330:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>