Annotation of OpenXM_contrib2/asir2000/engine/gfspn.c, Revision 1.3
1.3 ! noro 1: /* $OpenXM$ */
1.2 noro 2:
1.1 noro 3: #include "ca.h"
4: #include "base.h"
5:
6: UM current_mod_gfspn;
7:
8: void setmod_gfspn(p)
9: UM p;
10: {
11: current_mod_gfspn = p;
12: }
13:
14: void getmod_gfspn(up)
15: UM *up;
16: {
17: *up = current_mod_gfspn;
18: }
19:
20: void simpgfspn(n,r)
21: GFSPN n;
22: GFSPN *r;
23: {
24: UM t,q;
25:
26: if ( !n )
27: *r = 0;
28: else if ( NID(n) != N_GFSPN )
29: *r = n;
30: else {
31: t = UMALLOC(DEG(BDY(n)));
32: q = W_UMALLOC(DEG(BDY(n)));
33: cpyum(BDY(n),t);
34: divsfum(t,current_mod_gfspn,q);
35: MKGFSPN(t,*r);
36: }
37: }
38:
39: #define NZGFSPN(a) ((a)&&(OID(a)==O_N)&&(NID(a)==N_GFSPN))
40:
41: void ntogfspn(a,b)
42: Obj a;
43: GFSPN *b;
44: {
45: UM t;
46: GFS c;
47:
48: if ( !a || (OID(a)==O_N && NID(a) == N_GFSPN) )
49: *b = (GFSPN)a;
50: else if ( OID(a) == O_N && (NID(a) == N_GFS || NID(a) == N_Q) ) {
51: ntogfs((Num)a,&c);
52: if ( !b )
53: *b = 0;
54: else {
55: t = UMALLOC(0); ptosfum((P)c,t);
56: MKGFSPN(t,*b);
57: }
58: } else
59: error("ntogfspn : invalid argument");
60: }
61:
62: void addgfspn(a,b,c)
63: GFSPN a,b;
64: GFSPN *c;
65: {
66: UM t,q;
67: GFSPN z;
68: int d;
69:
70: ntogfspn((Obj)a,&z); a = z; ntogfspn((Obj)b,&z); b = z;
71: if ( !a )
72: *c = b;
73: else if ( !b )
74: *c = a;
75: else {
76: d = MAX(DEG(BDY(a)),DEG(BDY(b)));
77: t = UMALLOC(d);
78: q = W_UMALLOC(d);
79: addsfum(BDY(a),BDY(b),t);
80: divsfum(t,current_mod_gfspn,q);
81: MKGFSPN(t,*c);
82: }
83: }
84:
85: void subgfspn(a,b,c)
86: GFSPN a,b;
87: GFSPN *c;
88: {
89: UM t,q;
90: GFSPN z;
91: int d;
92:
93: ntogfspn((Obj)a,&z); a = z; ntogfspn((Obj)b,&z); b = z;
94: if ( !a )
95: chsgngfspn(b,c);
96: else if ( !b )
97: *c = a;
98: else {
99: d = MAX(DEG(BDY(a)),DEG(BDY(b)));
100: t = UMALLOC(d);
101: q = W_UMALLOC(d);
102: subsfum(BDY(a),BDY(b),t);
103: divsfum(t,current_mod_gfspn,q);
104: MKGFSPN(t,*c);
105: }
106: }
107:
108: extern int up_lazy;
109:
110: void mulgfspn(a,b,c)
111: GFSPN a,b;
112: GFSPN *c;
113: {
114: UM t,q;
115: GFSPN z;
116: int d;
117:
118: ntogfspn((Obj)a,&z); a = z; ntogfspn((Obj)b,&z); b = z;
119: if ( !a || !b )
120: *c = 0;
121: else {
122: d = DEG(BDY(a))+DEG(BDY(b));
123: t = UMALLOC(d);
124: q = W_UMALLOC(d);
125: mulsfum(BDY(a),BDY(b),t);
126: divsfum(t,current_mod_gfspn,q);
127: MKGFSPN(t,*c);
128: }
129: }
130:
131: void divgfspn(a,b,c)
132: GFSPN a,b;
133: GFSPN *c;
134: {
135: GFSPN z;
136: int d;
137: UM wb,wc,wd,we,t,q;
138:
139: ntogfspn((Obj)a,&z); a = z; ntogfspn((Obj)b,&z); b = z;
140: if ( !b )
141: error("divgfspn: division by 0");
142: else if ( !a )
143: *c = 0;
144: else {
145: wb = W_UMALLOC(DEG(BDY(b))); cpyum(BDY(b),wb);
146: d = DEG(current_mod_gfspn);
147: wc = W_UMALLOC(d); cpyum(current_mod_gfspn,wc);
148: wd = W_UMALLOC(2*d); we = W_UMALLOC(2*d);
149: /* wd*wb+we*wc=1 */
150: eucsfum(wb,wc,wd,we);
151: d = DEG(BDY(a))+DEG(wd);
152: t = UMALLOC(d);
153: q = W_UMALLOC(d);
154: mulsfum(a,wd,t);
155: divsfum(t,current_mod_gfspn,q);
156: MKGFSPN(t,*c);
157: }
158: }
159:
160: void invgfspn(b,c)
161: GFSPN b;
162: GFSPN *c;
163: {
164: GFSPN z;
165: int d;
166: UM wb,wc,wd,we,t;
167:
168: ntogfspn((Obj)b,&z); b = z;
169: if ( !b )
170: error("divgfspn: division by 0");
171: else {
172: wb = W_UMALLOC(DEG(BDY(b))); cpyum(BDY(b),wb);
173: d = DEG(current_mod_gfspn);
174: wc = W_UMALLOC(d); cpyum(current_mod_gfspn,wc);
175: wd = W_UMALLOC(2*d); we = W_UMALLOC(2*d);
176: /* wd*wb+we*wc=1 */
177: eucsfum(wb,wc,wd,we);
178: d = DEG(wd);
179: t = UMALLOC(d);
180: cpyum(wd,t);
181: MKGFSPN(t,*c);
182: }
183: }
184:
185: void chsgngfspn(a,c)
186: GFSPN a,*c;
187: {
188: GFSPN z;
189: int d;
190: struct oUM zero;
191: UM t;
192:
193: ntogfspn((Obj)a,&z); a = z;
194: if ( !a )
195: *c = 0;
196: else {
197: d = DEG(BDY(a));
198: t = UMALLOC(d);
199: DEG(&zero) = -1;
200: subsfum(&zero,BDY(a),t);
201: MKGFSPN(t,*c);
202: }
203: }
204:
205: void pwrgfspn(a,b,c)
206: GFSPN a;
207: Q b;
208: GFSPN *c;
209: {
210: GFSPN z;
211: UM t,x,y,q;
212: int d,k;
213: N e;
214:
215: ntogfspn((Obj)a,&z); a = z;
216: if ( !b ) {
217: t = UMALLOC(0); DEG(t) = 0; COEF(t)[0] = _onesf();
218: MKGFSPN(t,*c);
219: } else if ( !a )
220: *c = 0;
221: else {
222: d = DEG(current_mod_gfspn);
223:
224: /* y = 1 */
225: y = UMALLOC(d); DEG(y) = 0; COEF(y)[0] = _onesf();
226:
227: t = W_UMALLOC(2*d); q = W_UMALLOC(2*d);
228:
229: /* x = simplify(a) */
230: x = W_UMALLOC(DEG(BDY(a))); cpyum(BDY(a),x);
231: divsfum(x,current_mod_gfspn,q);
232: if ( DEG(x) < 0 ) {
233: *c = 0;
234: } else {
235: e = NM(b);
236: for ( k = n_bits(e)-1; k >= 0; k-- ) {
237: mulsfum(y,y,t);
238: divsfum(t,current_mod_gfspn,q);
239: cpyum(t,y);
240: if ( e->b[k/32] & (1<<(k%32)) ) {
241: mulsfum(y,x,t);
242: divsfum(t,current_mod_gfspn,q);
243: cpyum(t,y);
244: }
245: }
246: MKGFSPN(y,*c);
247: }
248: }
249: }
250:
251: int cmpgfspn(a,b)
252: GFSPN a,b;
253: {
254: GFSPN z;
255:
256: ntogfspn((Obj)a,&z); a = z; ntogfspn((Obj)b,&z); b = z;
257: if ( !a )
258: if ( !b )
259: return 0;
260: else
261: return -1;
262: else if ( !b )
263: return 1;
264: else
265: return compsfum(BDY(a),BDY(b));
266: }
267:
268: void randomgfspn(r)
269: GFSPN *r;
270: {
271: int i,d;
272: UM t;
273:
274: if ( !current_mod_gfspn )
275: error("randomgfspn : current_mod_gfspn is not set");
276: d = DEG(current_mod_gfspn);
277: t = UMALLOC(d-1);
278: randsfum(d,t);
279: degum(t,d-1);
280: if ( DEG(t) < 0 )
281: *r = 0;
282: else {
283: MKGFSPN(t,*r);
284: }
285: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>