![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to pdiv.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin|
Return to pdiv.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/builtin/pdiv.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3: #include "parse.h"
4:
5: void Psdiv(), Psrem(), Ptdiv(), Psqr(), Pinva_mod();
6: void Psdiv_gf2n(), Psrem_gf2n();
7: void Psdivm(), Psremm(), Psqrm();
8: void Psrem_mod();
9: void Pugcd();
10: void Purem();
11: void Pudiv();
12:
13: struct ftab pdiv_tab[] = {
14: {"sdiv",Psdiv,-3},
15: {"srem",Psrem,-3},
16: {"sdiv_gf2n",Psdiv_gf2n,2},
17: {"srem_gf2n",Psrem_gf2n,2},
18: {"sqr",Psqr,-3},
19: {"tdiv",Ptdiv,2},
20: {"udiv",Pudiv,2},
21: {"sdivm",Psdivm,-4},
22: {"sremm",Psremm,-4},
23: {"sqrm",Psqrm,-4},
24: {"inva_mod",Pinva_mod,3},
25: {"srem_mod",Psrem_mod,3},
26: {"ugcd",Pugcd,2},
27: {"urem",Purem,2},
28: {0,0,0},
29: };
30:
31: void Psdiv(arg,rp)
32: NODE arg;
33: Obj *rp;
34: {
35: P q,r,dnd,dnd1,dvr,dvr1;
36: V v;
37: VL vl;
38:
39: asir_assert(ARG0(arg),O_P,"sdiv");
40: asir_assert(ARG1(arg),O_P,"sdiv");
41: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg);
42: if ( argc(arg) == 3 ) {
43: v = VR((P)ARG2(arg));
44: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
45: reordvar(CO,v,&vl);
46: divsrp(vl,dnd1,dvr1,&q,&r);
47: restore_mvar(CO,q,v,(P *)rp);
48: } else
49: divsrp(CO,dnd,dvr,(P *)rp,&r);
50: }
51:
52: void Psrem(arg,rp)
53: NODE arg;
54: Obj *rp;
55: {
56: P q,r,dnd,dnd1,dvr,dvr1;
57: V v;
58: VL vl;
59:
60: asir_assert(ARG0(arg),O_P,"srem");
61: asir_assert(ARG1(arg),O_P,"srem");
62: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg);
63: if ( argc(arg) == 3 ) {
64: v = VR((P)ARG2(arg));
65: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
66: reordvar(CO,v,&vl);
67: divsrp(vl,dnd1,dvr1,&q,&r);
68: restore_mvar(CO,r,v,(P *)rp);
69: } else
70: divsrp(CO,dnd,dvr,&q,(P *)rp);
71: }
72:
73: void Psqr(arg,rp)
74: NODE arg;
75: LIST *rp;
76: {
77: P q,q1,r,r1,dnd,dnd1,dvr,dvr1;
78: NODE n,tn;
79: V v;
80: VL vl;
81:
82: asir_assert(ARG0(arg),O_P,"sqr");
83: asir_assert(ARG1(arg),O_P,"sqr");
84: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg);
85: if ( argc(arg) == 3 ) {
86: v = VR((P)ARG2(arg));
87: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
88: reordvar(CO,v,&vl);
89: divsrp(vl,dnd1,dvr1,&q1,&r1);
90: restore_mvar(CO,q1,v,&q); restore_mvar(CO,r1,v,&r);
91: } else
92: divsrp(CO,dnd,dvr,&q,&r);
93: MKNODE(tn,r,0); MKNODE(n,q,tn); MKLIST(*rp,n);
94: }
95:
96: void Psdiv_gf2n(arg,rp)
97: NODE arg;
98: GF2N *rp;
99: {
100: GF2N dnd,dvr;
101: UP2 q,r;
102:
103: dnd = (GF2N)ARG0(arg); dvr = (GF2N)ARG1(arg);
104: if ( !dvr )
105: error("sdiv_gf2n : division by 0");
106: else if ( !dnd )
107: *rp = 0;
108: else {
109: qrup2(dnd->body,dvr->body,&q,&r);
110: MKGF2N(q,*rp);
111: }
112: }
113:
114: void Psrem_gf2n(arg,rp)
115: NODE arg;
116: GF2N *rp;
117: {
118: GF2N dnd,dvr;
119: UP2 q,r;
120:
121: dnd = (GF2N)ARG0(arg); dvr = (GF2N)ARG1(arg);
122: if ( !dvr )
123: error("srem_gf2n : division by 0");
124: else if ( !dnd )
125: *rp = 0;
126: else {
127: qrup2(dnd->body,dvr->body,&q,&r);
128: MKGF2N(r,*rp);
129: }
130: }
131:
132: void Ptdiv(arg,rp)
133: NODE arg;
134: P *rp;
135: {
136: P p1,p2,q1,q2,q,c1,c2,c;
137:
138: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
139: asir_assert(p1,O_P,"tdiv");
140: asir_assert(p2,O_P,"tdiv");
141: if ( !p1 || !p2 )
142: *rp = 0;
143: else if ( (OID(p1) > O_P) || (OID(p2) > O_P ) )
144: *rp = 0;
145: else {
146: ptozp(p1,1,(Q *)&c1,&q1); ptozp(p2,1,(Q *)&c2,&q2);
147: if ( divtpz(CO,q1,q2,&q) ) {
148: divq((Q)c1,(Q)c2,(Q *)&c); mulp(CO,q,c,rp);
149: } else
150: *rp = 0;
151: }
152: }
153:
154: void Pudiv(arg,rp)
155: NODE arg;
156: LIST *rp;
157: {
158: P q,r,dnd,dvr;
159: NODE n,tn;
160:
161: asir_assert(ARG0(arg),O_P,"udiv");
162: asir_assert(ARG1(arg),O_P,"udiv");
163: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg);
164: udivpz(dnd,dvr,&q,&r);
165: MKNODE(tn,r,0); MKNODE(n,q,tn); MKLIST(*rp,n);
166: }
167:
168: void Psdivm(arg,rp)
169: NODE arg;
170: Obj *rp;
171: {
172: P q,r,dnd,dnd1,dndm,dvr,dvr1,dvrm,t;
173: V v;
174: VL vl;
175: int m;
176:
177: asir_assert(ARG0(arg),O_P,"sdivm");
178: asir_assert(ARG1(arg),O_P,"sdivm");
179: asir_assert(ARG2(arg),O_N,"sdivm");
180: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg); m = QTOS((Q)ARG2(arg));
181: if ( argc(arg) == 4 ) {
182: v = VR((P)ARG3(arg));
183: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
184: reordvar(CO,v,&vl);
185: ptomp(m,dnd1,&dndm); ptomp(m,dvr1,&dvrm);
186: divsrmp(vl,m,dndm,dvrm,&t,&r); mptop(t,&q);
187: restore_mvar(CO,q,v,(P *)rp);
188: } else {
189: ptomp(m,dnd,&dndm); ptomp(m,dvr,&dvrm);
190: divsrmp(CO,m,dndm,dvrm,&t,&r); mptop(t,(P *)rp);
191: }
192: }
193:
194: void Psremm(arg,rp)
195: NODE arg;
196: Obj *rp;
197: {
198: P q,r,dnd,dnd1,dndm,dvr,dvr1,dvrm,t;
199: V v;
200: VL vl;
201: int m;
202:
203: asir_assert(ARG0(arg),O_P,"sremm");
204: asir_assert(ARG1(arg),O_P,"sremm");
205: asir_assert(ARG2(arg),O_N,"sremm");
206: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg); m = QTOS((Q)ARG2(arg));
207: if ( argc(arg) == 4 ) {
208: v = VR((P)ARG3(arg));
209: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
210: reordvar(CO,v,&vl);
211: ptomp(m,dnd1,&dndm); ptomp(m,dvr1,&dvrm);
212: divsrmp(vl,m,dndm,dvrm,&q,&t); mptop(t,&r);
213: restore_mvar(CO,r,v,(P *)rp);
214: } else {
215: ptomp(m,dnd,&dndm); ptomp(m,dvr,&dvrm);
216: divsrmp(CO,m,dndm,dvrm,&q,&t); mptop(t,(P *)rp);
217: }
218: }
219:
220: void Psqrm(arg,rp)
221: NODE arg;
222: LIST *rp;
223: {
224: P q,q1,r,r1,dnd,dnd1,dndm,dvr,dvr1,dvrm;
225: NODE n,tn;
226: V v;
227: VL vl;
228: int m;
229:
230: asir_assert(ARG0(arg),O_P,"sqrm");
231: asir_assert(ARG1(arg),O_P,"sqrm");
232: asir_assert(ARG2(arg),O_N,"sqrm");
233: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg); m = QTOS((Q)ARG2(arg));
234: if ( argc(arg) == 4 ) {
235: v = VR((P)ARG3(arg));
236: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
237: reordvar(CO,v,&vl);
238: ptomp(m,dnd1,&dndm); ptomp(m,dvr1,&dvrm);
239: divsrp(vl,dndm,dvrm,&q,&r); mptop(q,&q1); mptop(r,&r1);
240: restore_mvar(CO,q1,v,&q); restore_mvar(CO,r1,v,&r);
241: } else {
242: ptomp(m,dnd,&dndm); ptomp(m,dvr,&dvrm);
243: divsrp(CO,dnd,dvr,&q1,&r1); mptop(q1,&q); mptop(r1,&r);
244: }
245: MKNODE(tn,r,0); MKNODE(n,q,tn); MKLIST(*rp,n);
246: }
247:
248: void Pinva_mod(arg,rp)
249: NODE arg;
250: P *rp;
251: {
252: P dp,f;
253: Q q;
254: int n,i;
255: int mod;
256: V v;
257: UM wf,wdp,winv;
258:
259: asir_assert(ARG0(arg),O_P,"gcda_mod");
260: asir_assert(ARG1(arg),O_N,"gcda_mod");
261: asir_assert(ARG2(arg),O_P,"gcda_mod");
262: dp = (P)ARG0(arg);
263: mod = QTOS((Q)ARG1(arg));
264: f = (P)ARG2(arg);
265: if ( NUM(f) ) {
266: i = invm(rem(NM((Q)f),mod),mod);
267: STOQ(i,q); *rp = (P)q;
268: } else {
269: v = VR(dp);
270: n = MAX(UDEG(dp),UDEG(f));
271: wf = W_UMALLOC(n); wdp = W_UMALLOC(n);
272: winv = W_UMALLOC(n);
273: ptoum(mod,f,wf); ptoum(mod,dp,wdp);
274: invum(mod,wdp,wf,winv);
275: if ( DEG(winv) < 0 )
276: *rp = 0;
277: else {
278: umtop(v,winv,rp);
279: }
280: }
281: }
282:
283: void Psrem_mod(arg,rp)
284: NODE arg;
285: P *rp;
286: {
287: P p1,p2;
288: int n,dr;
289: int mod;
290: V v;
291: UM wp1,wp2,q;
292:
293: asir_assert(ARG0(arg),O_P,"srem_mod");
294: asir_assert(ARG1(arg),O_P,"srem_mod");
295: asir_assert(ARG2(arg),O_N,"srem_mod");
296: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg); mod = QTOS((Q)ARG2(arg));
297: if ( !p1 || NUM(p1) )
298: *rp = p1;
299: else {
300: v = VR(p1);
301: n = MAX(UDEG(p1),UDEG(p2));
302: wp1 = W_UMALLOC(n); wp2 = W_UMALLOC(n); q = W_UMALLOC(n);
303: ptoum(mod,p1,wp1); ptoum(mod,p2,wp2);
304: dr = divum(mod,wp1,wp2,q);
305: if ( ( DEG(wp1) = dr ) == -1 )
306: *rp = 0;
307: else
308: umtop(v,wp1,rp);
309: }
310: }
311:
312: void Purem(arg,rp)
313: NODE arg;
314: P *rp;
315: {
316: asir_assert(ARG0(arg),O_P,"urem");
317: asir_assert(ARG1(arg),O_P,"urem");
318: uremp((P)ARG0(arg),(P)ARG1(arg),rp);
319: }
320:
321: void Pugcd(arg,rp)
322: NODE arg;
323: P *rp;
324: {
325: asir_assert(ARG0(arg),O_P,"ugcd");
326: asir_assert(ARG1(arg),O_P,"ugcd");
327: ugcdp((P)ARG0(arg),(P)ARG1(arg),rp);
328: }
329:
330: void invum(mod,dp,f,inv)
331: int mod;
332: UM dp,f,inv;
333: {
334: UM g1,g2,a1,a2,a3,wm,q,tum;
335: int d,dr;
336:
337: d = DEG(dp)+DEG(f)+10;
338: g1 = W_UMALLOC(d); g2 = W_UMALLOC(d); a1 = W_UMALLOC(d);
339: a2 = W_UMALLOC(d); a3 = W_UMALLOC(d); wm = W_UMALLOC(d);
340: q = W_UMALLOC(d);
341: DEG(a1) = 0; COEF(a1)[0] = 1; DEG(a2) = -1;
342: cpyum(f,g1); cpyum(dp,g2);
343: while ( 1 ) {
344: dr = divum(mod,g1,g2,q); tum = g1; g1 = g2; g2 = tum;
345: if ( ( DEG(g2) = dr ) == -1 )
346: break;
347: mulum(mod,a2,q,wm); subum(mod,a1,wm,a3); dr = divum(mod,a3,dp,q);
348: tum = a1; a1 = a2; a2 = a3; a3 = tum; DEG(a3) = dr;
349: }
350: if ( DEG(g1) != 0 )
351: DEG(inv) = -1;
352: else if ( COEF(g1)[0] != 1 )
353: mulsum(mod,a2,invm(COEF(g1)[0],mod),inv);
354: else
355: cpyum(a2,inv);
356: }