Annotation of OpenXM_contrib2/asir2000/builtin/pdiv.c, Revision 1.8
1.3 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.4 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.3 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.8 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/pdiv.c,v 1.7 2002/02/08 08:28:00 noro Exp $
1.3 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "parse.h"
52:
1.8 ! noro 53: void Psdiv(), Psrem(), Ptdiv(), Psqr(), Pinva_mod(), Pprem();
1.7 noro 54: void Psdiv_gf2n(), Psrem_gf2n(), Pgcd_gf2n();
1.1 noro 55: void Psdivm(), Psremm(), Psqrm();
56: void Psrem_mod();
57: void Pugcd();
58: void Purem();
59: void Pudiv();
60:
61: struct ftab pdiv_tab[] = {
62: {"sdiv",Psdiv,-3},
63: {"srem",Psrem,-3},
1.8 ! noro 64: {"prem",Pprem,-3},
1.1 noro 65: {"sdiv_gf2n",Psdiv_gf2n,2},
66: {"srem_gf2n",Psrem_gf2n,2},
1.7 noro 67: {"gcd_gf2n",Pgcd_gf2n,2},
1.1 noro 68: {"sqr",Psqr,-3},
1.5 noro 69: {"tdiv",Ptdiv,-3},
1.1 noro 70: {"udiv",Pudiv,2},
71: {"sdivm",Psdivm,-4},
72: {"sremm",Psremm,-4},
73: {"sqrm",Psqrm,-4},
74: {"inva_mod",Pinva_mod,3},
75: {"srem_mod",Psrem_mod,3},
76: {"ugcd",Pugcd,2},
77: {"urem",Purem,2},
78: {0,0,0},
79: };
80:
81: void Psdiv(arg,rp)
82: NODE arg;
83: Obj *rp;
84: {
85: P q,r,dnd,dnd1,dvr,dvr1;
86: V v;
87: VL vl;
88:
89: asir_assert(ARG0(arg),O_P,"sdiv");
90: asir_assert(ARG1(arg),O_P,"sdiv");
91: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg);
92: if ( argc(arg) == 3 ) {
93: v = VR((P)ARG2(arg));
94: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
95: reordvar(CO,v,&vl);
96: divsrp(vl,dnd1,dvr1,&q,&r);
97: restore_mvar(CO,q,v,(P *)rp);
98: } else
99: divsrp(CO,dnd,dvr,(P *)rp,&r);
100: }
101:
102: void Psrem(arg,rp)
103: NODE arg;
104: Obj *rp;
105: {
106: P q,r,dnd,dnd1,dvr,dvr1;
107: V v;
108: VL vl;
109:
110: asir_assert(ARG0(arg),O_P,"srem");
111: asir_assert(ARG1(arg),O_P,"srem");
112: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg);
113: if ( argc(arg) == 3 ) {
114: v = VR((P)ARG2(arg));
115: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
116: reordvar(CO,v,&vl);
117: divsrp(vl,dnd1,dvr1,&q,&r);
118: restore_mvar(CO,r,v,(P *)rp);
119: } else
120: divsrp(CO,dnd,dvr,&q,(P *)rp);
1.8 ! noro 121: }
! 122:
! 123: void Pprem(arg,rp)
! 124: NODE arg;
! 125: P *rp;
! 126: {
! 127: P q,r,dnd,dnd1,dvr,dvr1;
! 128: V v;
! 129: VL vl;
! 130:
! 131: asir_assert(ARG0(arg),O_P,"prem");
! 132: asir_assert(ARG1(arg),O_P,"prem");
! 133: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg);
! 134: if ( argc(arg) == 3 ) {
! 135: v = VR((P)ARG2(arg));
! 136: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
! 137: reordvar(CO,v,&vl);
! 138: premp(vl,dnd1,dvr1,&r);
! 139: restore_mvar(CO,r,v,rp);
! 140: } else
! 141: premp(CO,dnd,dvr,rp);
1.1 noro 142: }
143:
144: void Psqr(arg,rp)
145: NODE arg;
146: LIST *rp;
147: {
148: P q,q1,r,r1,dnd,dnd1,dvr,dvr1;
149: NODE n,tn;
150: V v;
151: VL vl;
152:
153: asir_assert(ARG0(arg),O_P,"sqr");
154: asir_assert(ARG1(arg),O_P,"sqr");
155: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg);
156: if ( argc(arg) == 3 ) {
157: v = VR((P)ARG2(arg));
158: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
159: reordvar(CO,v,&vl);
160: divsrp(vl,dnd1,dvr1,&q1,&r1);
161: restore_mvar(CO,q1,v,&q); restore_mvar(CO,r1,v,&r);
162: } else
163: divsrp(CO,dnd,dvr,&q,&r);
164: MKNODE(tn,r,0); MKNODE(n,q,tn); MKLIST(*rp,n);
165: }
166:
167: void Psdiv_gf2n(arg,rp)
168: NODE arg;
169: GF2N *rp;
170: {
171: GF2N dnd,dvr;
172: UP2 q,r;
173:
174: dnd = (GF2N)ARG0(arg); dvr = (GF2N)ARG1(arg);
175: if ( !dvr )
176: error("sdiv_gf2n : division by 0");
177: else if ( !dnd )
178: *rp = 0;
179: else {
180: qrup2(dnd->body,dvr->body,&q,&r);
181: MKGF2N(q,*rp);
182: }
183: }
184:
185: void Psrem_gf2n(arg,rp)
186: NODE arg;
187: GF2N *rp;
188: {
189: GF2N dnd,dvr;
190: UP2 q,r;
191:
192: dnd = (GF2N)ARG0(arg); dvr = (GF2N)ARG1(arg);
193: if ( !dvr )
194: error("srem_gf2n : division by 0");
195: else if ( !dnd )
196: *rp = 0;
197: else {
198: qrup2(dnd->body,dvr->body,&q,&r);
199: MKGF2N(r,*rp);
1.7 noro 200: }
201: }
202:
203: void Pgcd_gf2n(arg,rp)
204: NODE arg;
205: GF2N *rp;
206: {
207: GF2N p1,p2;
208: UP2 gcd;
209:
210: p1 = (GF2N)ARG0(arg); p2 = (GF2N)ARG1(arg);
211: if ( !p1 )
212: *rp = p2;
213: else if ( !p2 )
214: *rp = p1;
215: else {
216: gcdup2(p1->body,p2->body,&gcd);
217: MKGF2N(gcd,*rp);
1.1 noro 218: }
219: }
220:
221: void Ptdiv(arg,rp)
222: NODE arg;
223: P *rp;
224: {
225: P p1,p2,q1,q2,q,c1,c2,c;
1.5 noro 226: int m;
1.1 noro 227:
228: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
229: asir_assert(p1,O_P,"tdiv");
230: asir_assert(p2,O_P,"tdiv");
231: if ( !p1 || !p2 )
232: *rp = 0;
233: else if ( (OID(p1) > O_P) || (OID(p2) > O_P ) )
234: *rp = 0;
1.5 noro 235: else if ( argc(arg) == 3 ) {
236: m = QTOS((Q)ARG2(arg));
237: ptomp(m,p1,&q1); ptomp(m,p2,&q2);
238: if ( divtmp(CO,m,q1,q2,&q) )
239: mptop(q,rp);
240: else
241: *rp = 0;
1.6 noro 242: } else if ( qpcheck((Obj)p1) && qpcheck((Obj)p2) ) {
1.1 noro 243: ptozp(p1,1,(Q *)&c1,&q1); ptozp(p2,1,(Q *)&c2,&q2);
244: if ( divtpz(CO,q1,q2,&q) ) {
245: divq((Q)c1,(Q)c2,(Q *)&c); mulp(CO,q,c,rp);
246: } else
1.6 noro 247: *rp = 0;
248: } else {
249: if ( !divtp(CO,p1,p2,rp) )
1.1 noro 250: *rp = 0;
251: }
252: }
253:
254: void Pudiv(arg,rp)
255: NODE arg;
256: LIST *rp;
257: {
258: P q,r,dnd,dvr;
259: NODE n,tn;
260:
261: asir_assert(ARG0(arg),O_P,"udiv");
262: asir_assert(ARG1(arg),O_P,"udiv");
263: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg);
264: udivpz(dnd,dvr,&q,&r);
265: MKNODE(tn,r,0); MKNODE(n,q,tn); MKLIST(*rp,n);
266: }
267:
268: void Psdivm(arg,rp)
269: NODE arg;
270: Obj *rp;
271: {
272: P q,r,dnd,dnd1,dndm,dvr,dvr1,dvrm,t;
273: V v;
274: VL vl;
275: int m;
276:
277: asir_assert(ARG0(arg),O_P,"sdivm");
278: asir_assert(ARG1(arg),O_P,"sdivm");
279: asir_assert(ARG2(arg),O_N,"sdivm");
280: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg); m = QTOS((Q)ARG2(arg));
281: if ( argc(arg) == 4 ) {
282: v = VR((P)ARG3(arg));
283: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
284: reordvar(CO,v,&vl);
285: ptomp(m,dnd1,&dndm); ptomp(m,dvr1,&dvrm);
286: divsrmp(vl,m,dndm,dvrm,&t,&r); mptop(t,&q);
287: restore_mvar(CO,q,v,(P *)rp);
288: } else {
289: ptomp(m,dnd,&dndm); ptomp(m,dvr,&dvrm);
290: divsrmp(CO,m,dndm,dvrm,&t,&r); mptop(t,(P *)rp);
291: }
292: }
293:
294: void Psremm(arg,rp)
295: NODE arg;
296: Obj *rp;
297: {
298: P q,r,dnd,dnd1,dndm,dvr,dvr1,dvrm,t;
299: V v;
300: VL vl;
301: int m;
302:
303: asir_assert(ARG0(arg),O_P,"sremm");
304: asir_assert(ARG1(arg),O_P,"sremm");
305: asir_assert(ARG2(arg),O_N,"sremm");
306: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg); m = QTOS((Q)ARG2(arg));
307: if ( argc(arg) == 4 ) {
308: v = VR((P)ARG3(arg));
309: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
310: reordvar(CO,v,&vl);
311: ptomp(m,dnd1,&dndm); ptomp(m,dvr1,&dvrm);
312: divsrmp(vl,m,dndm,dvrm,&q,&t); mptop(t,&r);
313: restore_mvar(CO,r,v,(P *)rp);
314: } else {
315: ptomp(m,dnd,&dndm); ptomp(m,dvr,&dvrm);
316: divsrmp(CO,m,dndm,dvrm,&q,&t); mptop(t,(P *)rp);
317: }
318: }
319:
320: void Psqrm(arg,rp)
321: NODE arg;
322: LIST *rp;
323: {
324: P q,q1,r,r1,dnd,dnd1,dndm,dvr,dvr1,dvrm;
325: NODE n,tn;
326: V v;
327: VL vl;
328: int m;
329:
330: asir_assert(ARG0(arg),O_P,"sqrm");
331: asir_assert(ARG1(arg),O_P,"sqrm");
332: asir_assert(ARG2(arg),O_N,"sqrm");
333: dnd = (P)ARG0(arg); dvr = (P)ARG1(arg); m = QTOS((Q)ARG2(arg));
334: if ( argc(arg) == 4 ) {
335: v = VR((P)ARG3(arg));
336: change_mvar(CO,dnd,v,&dnd1); change_mvar(CO,dvr,v,&dvr1);
337: reordvar(CO,v,&vl);
338: ptomp(m,dnd1,&dndm); ptomp(m,dvr1,&dvrm);
1.2 noro 339: divsrmp(vl,m,dndm,dvrm,&q,&r); mptop(q,&q1); mptop(r,&r1);
1.1 noro 340: restore_mvar(CO,q1,v,&q); restore_mvar(CO,r1,v,&r);
341: } else {
342: ptomp(m,dnd,&dndm); ptomp(m,dvr,&dvrm);
1.2 noro 343: divsrmp(CO,m,dndm,dvrm,&q1,&r1); mptop(q1,&q); mptop(r1,&r);
1.1 noro 344: }
345: MKNODE(tn,r,0); MKNODE(n,q,tn); MKLIST(*rp,n);
346: }
347:
348: void Pinva_mod(arg,rp)
349: NODE arg;
350: P *rp;
351: {
352: P dp,f;
353: Q q;
354: int n,i;
355: int mod;
356: V v;
357: UM wf,wdp,winv;
358:
359: asir_assert(ARG0(arg),O_P,"gcda_mod");
360: asir_assert(ARG1(arg),O_N,"gcda_mod");
361: asir_assert(ARG2(arg),O_P,"gcda_mod");
362: dp = (P)ARG0(arg);
363: mod = QTOS((Q)ARG1(arg));
364: f = (P)ARG2(arg);
365: if ( NUM(f) ) {
366: i = invm(rem(NM((Q)f),mod),mod);
367: STOQ(i,q); *rp = (P)q;
368: } else {
369: v = VR(dp);
370: n = MAX(UDEG(dp),UDEG(f));
371: wf = W_UMALLOC(n); wdp = W_UMALLOC(n);
372: winv = W_UMALLOC(n);
373: ptoum(mod,f,wf); ptoum(mod,dp,wdp);
374: invum(mod,wdp,wf,winv);
375: if ( DEG(winv) < 0 )
376: *rp = 0;
377: else {
378: umtop(v,winv,rp);
379: }
380: }
381: }
382:
383: void Psrem_mod(arg,rp)
384: NODE arg;
385: P *rp;
386: {
387: P p1,p2;
388: int n,dr;
389: int mod;
390: V v;
391: UM wp1,wp2,q;
392:
393: asir_assert(ARG0(arg),O_P,"srem_mod");
394: asir_assert(ARG1(arg),O_P,"srem_mod");
395: asir_assert(ARG2(arg),O_N,"srem_mod");
396: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg); mod = QTOS((Q)ARG2(arg));
397: if ( !p1 || NUM(p1) )
398: *rp = p1;
399: else {
400: v = VR(p1);
401: n = MAX(UDEG(p1),UDEG(p2));
402: wp1 = W_UMALLOC(n); wp2 = W_UMALLOC(n); q = W_UMALLOC(n);
403: ptoum(mod,p1,wp1); ptoum(mod,p2,wp2);
404: dr = divum(mod,wp1,wp2,q);
405: if ( ( DEG(wp1) = dr ) == -1 )
406: *rp = 0;
407: else
408: umtop(v,wp1,rp);
409: }
410: }
411:
412: void Purem(arg,rp)
413: NODE arg;
414: P *rp;
415: {
416: asir_assert(ARG0(arg),O_P,"urem");
417: asir_assert(ARG1(arg),O_P,"urem");
418: uremp((P)ARG0(arg),(P)ARG1(arg),rp);
419: }
420:
421: void Pugcd(arg,rp)
422: NODE arg;
423: P *rp;
424: {
425: asir_assert(ARG0(arg),O_P,"ugcd");
426: asir_assert(ARG1(arg),O_P,"ugcd");
427: ugcdp((P)ARG0(arg),(P)ARG1(arg),rp);
428: }
429:
430: void invum(mod,dp,f,inv)
431: int mod;
432: UM dp,f,inv;
433: {
434: UM g1,g2,a1,a2,a3,wm,q,tum;
435: int d,dr;
436:
437: d = DEG(dp)+DEG(f)+10;
438: g1 = W_UMALLOC(d); g2 = W_UMALLOC(d); a1 = W_UMALLOC(d);
439: a2 = W_UMALLOC(d); a3 = W_UMALLOC(d); wm = W_UMALLOC(d);
440: q = W_UMALLOC(d);
441: DEG(a1) = 0; COEF(a1)[0] = 1; DEG(a2) = -1;
442: cpyum(f,g1); cpyum(dp,g2);
443: while ( 1 ) {
444: dr = divum(mod,g1,g2,q); tum = g1; g1 = g2; g2 = tum;
445: if ( ( DEG(g2) = dr ) == -1 )
446: break;
447: mulum(mod,a2,q,wm); subum(mod,a1,wm,a3); dr = divum(mod,a3,dp,q);
448: tum = a1; a1 = a2; a2 = a3; a3 = tum; DEG(a3) = dr;
449: }
450: if ( DEG(g1) != 0 )
451: DEG(inv) = -1;
452: else if ( COEF(g1)[0] != 1 )
453: mulsum(mod,a2,invm(COEF(g1)[0],mod),inv);
454: else
455: cpyum(a2,inv);
456: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>