Annotation of OpenXM_contrib2/asir2000/engine/Q.c, Revision 1.1.1.1
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/Q.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3: #include "base.h"
4: #include "inline.h"
5:
6: void addq(n1,n2,nr)
7: Q n1,n2,*nr;
8: {
9: N nm,dn,nm1,nm2,nm3,dn0,dn1,dn2,g,g1,g0,m;
10: int sgn;
11: unsigned t,t1;
12:
13: if ( !n1 )
14: *nr = n2;
15: else if ( !n2 )
16: *nr = n1;
17: else if ( INT(n1) )
18: if ( INT(n2) ) {
19: nm1 = NM(n1); nm2 = NM(n2);
20: if ( SGN(n1) == SGN(n2) ) {
21: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
22: t1 = BD(nm1)[0]; t = t1+BD(nm2)[0];
23: if ( t < t1 ) {
24: m = NALLOC(2); PL(m) = 2; BD(m)[0] = t; BD(m)[1] = 1;
25: } else
26: UTON(t,m);
27: } else
28: addn(NM(n1),NM(n2),&m);
29: NTOQ(m,SGN(n1),*nr);
30: } else {
31: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
32: t1 = BD(nm1)[0]; t = t1-BD(nm2)[0];
33: if ( !t )
34: sgn = 0;
35: else if ( t > t1 ) {
36: sgn = -1; t = -((int)t); UTON(t,m);
37: } else {
38: sgn = 1; UTON(t,m);
39: }
40: } else
41: sgn = subn(NM(n1),NM(n2),&m);
42: if ( sgn ) {
43: if ( SGN(n1) == sgn )
44: NTOQ(m,1,*nr);
45: else
46: NTOQ(m,-1,*nr);
47: } else
48: *nr = 0;
49: }
50: } else {
51: kmuln(NM(n1),DN(n2),&m);
52: if ( SGN(n1) == SGN(n2) ) {
53: sgn = SGN(n1); addn(m,NM(n2),&nm);
54: } else
55: sgn = SGN(n1)*subn(m,NM(n2),&nm);
56: if ( sgn ) {
57: dn = DN(n2); NDTOQ(nm,dn,sgn,*nr);
58: } else
59: *nr = 0;
60: }
61: else if ( INT(n2) ) {
62: kmuln(NM(n2),DN(n1),&m);
63: if ( SGN(n1) == SGN(n2) ) {
64: sgn = SGN(n1); addn(m,NM(n1),&nm);
65: } else
66: sgn = SGN(n1)*subn(NM(n1),m,&nm);
67: if ( sgn ) {
68: dn = DN(n1); NDTOQ(nm,dn,sgn,*nr);
69: } else
70: *nr = 0;
71: } else {
72: gcdn(DN(n1),DN(n2),&g); divsn(DN(n1),g,&dn1); divsn(DN(n2),g,&dn2);
73: kmuln(NM(n1),dn2,&nm1); kmuln(NM(n2),dn1,&nm2);
74: if ( SGN(n1) == SGN(n2) ) {
75: sgn = SGN(n1); addn(nm1,nm2,&nm3);
76: } else
77: sgn = SGN(n1)*subn(nm1,nm2,&nm3);
78: if ( sgn ) {
79: gcdn(nm3,g,&g1); divsn(nm3,g1,&nm); divsn(g,g1,&g0);
80: kmuln(dn1,dn2,&dn0); kmuln(g0,dn0,&dn);
81: if ( UNIN(dn) )
82: NTOQ(nm,sgn,*nr);
83: else
84: NDTOQ(nm,dn,sgn,*nr);
85: } else
86: *nr = 0;
87: }
88: }
89:
90: void subq(n1,n2,nr)
91: Q n1,n2,*nr;
92: {
93: Q m;
94:
95: if ( !n1 )
96: if ( !n2 )
97: *nr = 0;
98: else {
99: DUPQ(n2,*nr); SGN(*nr) = -SGN(n2);
100: }
101: else if ( !n2 )
102: *nr = n1;
103: else if ( n1 == n2 )
104: *nr = 0;
105: else {
106: DUPQ(n2,m); SGN(m) = -SGN(n2); addq(n1,m,nr);
107: }
108: }
109:
110: void mulq(n1,n2,nr)
111: Q n1,n2,*nr;
112: {
113: N nm,nm1,nm2,dn,dn1,dn2,g;
114: int sgn;
115: unsigned u,l;
116:
117: if ( !n1 || !n2 ) *nr = 0;
118: else if ( INT(n1) )
119: if ( INT(n2) ) {
120: nm1 = NM(n1); nm2 = NM(n2);
121: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
122: DM(BD(NM(n1))[0],BD(NM(n2))[0],u,l)
123: if ( u ) {
124: nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = l; BD(nm)[1] = u;
125: } else {
126: nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = l;
127: }
128: } else
129: kmuln(nm1,nm2,&nm);
130: sgn = (SGN(n1)==SGN(n2)?1:-1); NTOQ(nm,sgn,*nr);
131: } else {
132: gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn);
133: kmuln(nm1,NM(n2),&nm); sgn = SGN(n1)*SGN(n2);
134: if ( UNIN(dn) )
135: NTOQ(nm,sgn,*nr);
136: else
137: NDTOQ(nm,dn,sgn,*nr);
138: }
139: else if ( INT(n2) ) {
140: gcdn(NM(n2),DN(n1),&g); divsn(NM(n2),g,&nm2); divsn(DN(n1),g,&dn);
141: kmuln(nm2,NM(n1),&nm); sgn = SGN(n1)*SGN(n2);
142: if ( UNIN(dn) )
143: NTOQ(nm,sgn,*nr);
144: else
145: NDTOQ(nm,dn,sgn,*nr);
146: } else {
147: gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn2);
148: gcdn(DN(n1),NM(n2),&g); divsn(DN(n1),g,&dn1); divsn(NM(n2),g,&nm2);
149: kmuln(nm1,nm2,&nm); kmuln(dn1,dn2,&dn); sgn = SGN(n1) * SGN(n2);
150: if ( UNIN(dn) )
151: NTOQ(nm,sgn,*nr);
152: else
153: NDTOQ(nm,dn,sgn,*nr);
154: }
155: }
156:
157: void divq(n1,n2,nq)
158: Q n1,n2,*nq;
159: {
160: Q m;
161:
162: if ( !n2 ) {
163: error("division by 0");
164: *nq = 0;
165: return;
166: } else if ( !n1 )
167: *nq = 0;
168: else if ( n1 == n2 )
169: *nq = ONE;
170: else {
171: invq(n2,&m); mulq(n1,m,nq);
172: }
173: }
174:
175: void invq(n,nr)
176: Q n,*nr;
177: {
178: N nm,dn;
179:
180: if ( INT(n) )
181: if ( UNIN(NM(n)) )
182: if ( SGN(n) > 0 )
183: *nr = ONE;
184: else {
185: nm = ONEN; NTOQ(nm,SGN(n),*nr);
186: }
187: else {
188: nm = ONEN; dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
189: }
190: else if ( UNIN(NM(n)) ) {
191: nm = DN(n); NTOQ(nm,SGN(n),*nr);
192: } else {
193: nm = DN(n); dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
194: }
195: }
196:
197: void chsgnq(n,nr)
198: Q n,*nr;
199: {
200: Q t;
201:
202: if ( !n )
203: *nr = 0;
204: else {
205: DUPQ(n,t); SGN(t) = -SGN(t); *nr = t;
206: }
207: }
208:
209: void pwrq(n1,n,nr)
210: Q n1,n,*nr;
211: {
212: N nm,dn;
213: int sgn;
214:
215: if ( !n || UNIQ(n1) )
216: *nr = ONE;
217: else if ( !n1 )
218: *nr = 0;
219: else if ( !INT(n) ) {
220: error("can't calculate fractional power."); *nr = 0;
221: } else if ( MUNIQ(n1) )
222: *nr = BD(NM(n))[0]%2 ? n1 : ONE;
223: else if ( PL(NM(n)) > 1 ) {
224: error("exponent too big."); *nr = 0;
225: } else {
226: sgn = ((SGN(n1)>0)||EVENN(NM(n))?1:-1);
227: pwrn(NM(n1),BD(NM(n))[0],&nm);
228: if ( INT(n1) ) {
229: if ( UNIN(nm) )
230: if ( sgn == 1 )
231: *nr = ONE;
232: else
233: NTOQ(ONEN,sgn,*nr);
234: else if ( SGN(n) > 0 )
235: NTOQ(nm,sgn,*nr);
236: else
237: NDTOQ(ONEN,nm,sgn,*nr);
238: } else {
239: pwrn(DN(n1),BD(NM(n))[0],&dn);
240: if ( SGN(n) > 0 )
241: NDTOQ(nm,dn,sgn,*nr);
242: else if ( UNIN(nm) )
243: NTOQ(dn,sgn,*nr);
244: else
245: NDTOQ(dn,nm,sgn,*nr);
246: }
247: }
248: }
249:
250: int cmpq(q1,q2)
251: Q q1,q2;
252: {
253: int sgn;
254: N t,s;
255:
256: if ( !q1 )
257: if ( !q2 )
258: return ( 0 );
259: else
260: return ( -SGN(q2) );
261: else if ( !q2 )
262: return ( SGN(q1) );
263: else if ( SGN(q1) != SGN(q2) )
264: return ( SGN(q1) );
265: else if ( INT(q1) )
266: if ( INT(q2) ) {
267: sgn = cmpn(NM(q1),NM(q2));
268: if ( !sgn )
269: return ( 0 );
270: else
271: return ( SGN(q1)==sgn?1:-1 );
272: } else {
273: kmuln(NM(q1),DN(q2),&t); sgn = cmpn(t,NM(q2));
274: return ( SGN(q1) * sgn );
275: }
276: else if ( INT(q2) ) {
277: kmuln(NM(q2),DN(q1),&t); sgn = cmpn(NM(q1),t);
278: return ( SGN(q1) * sgn );
279: } else {
280: kmuln(NM(q1),DN(q2),&t); kmuln(NM(q2),DN(q1),&s); sgn = cmpn(t,s);
281: return ( SGN(q1) * sgn );
282: }
283: }
284:
285: void remq(n,m,nr)
286: Q n,m;
287: Q *nr;
288: {
289: N q,r,s;
290:
291: if ( !n ) {
292: *nr = 0;
293: return;
294: }
295: divn(NM(n),NM(m),&q,&r);
296: if ( !r )
297: *nr = 0;
298: else if ( SGN(n) > 0 )
299: NTOQ(r,1,*nr);
300: else {
301: subn(NM(m),r,&s); NTOQ(s,1,*nr);
302: }
303: }
304:
305: void mkbc(n,t)
306: int n;
307: Q *t;
308: {
309: int i;
310: N c,d;
311:
312: for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
313: c = NALLOC(1); PL(c) = 1; BD(c)[0] = n-i+1;
314: kmuln(NM(t[i-1]),c,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
315: }
316: for ( ; i <= n; i++ )
317: t[i] = t[n-i];
318: }
319:
320: void factorial(n,r)
321: int n;
322: Q *r;
323: {
324: Q t,iq,s;
325: unsigned int i,m,m0;
326: unsigned int c;
327:
328: if ( !n )
329: *r = ONE;
330: else if ( n < 0 )
331: *r = 0;
332: else {
333: for ( i = 1, t = ONE; (int)i <= n; ) {
334: for ( m0 = m = 1, c = 0; ((int)i <= n) && !c; i++ ) {
335: m0 = m; DM(m0,i,c,m)
336: }
337: if ( c ) {
338: m = m0; i--;
339: }
340: UTOQ(m,iq); mulq(t,iq,&s); t = s;
341: }
342: *r = t;
343: }
344: }
345:
346: void invl(a,mod,ar)
347: Q a,mod,*ar;
348: {
349: Q f1,f2,a1,a2,q,m,s;
350: N qn,rn;
351:
352: a1 = ONE; a2 = 0;
353: DUPQ(a,f1); SGN(f1) = 1; f2 = mod;
354:
355: while ( 1 ) {
356: divn(NM(f1),NM(f2),&qn,&rn);
357: if ( !qn )
358: q = 0;
359: else
360: NTOQ(qn,1,q);
361:
362: f1 = f2;
363: if ( !rn )
364: break;
365: else
366: NTOQ(rn,1,f2);
367:
368: mulq(a2,q,&m); subq(a1,m,&s); a1 = a2;
369: if ( !s )
370: a2 = 0;
371: else
372: remq(s,mod,&a2);
373: }
374: if ( SGN(a) < 0 )
375: chsgnq(a2,&m);
376: else
377: m = a2;
378:
379: if ( SGN(m) >= 0 )
380: *ar = m;
381: else
382: addq(m,mod,ar);
383: }
384:
385: int kara_mag=100;
386:
387: void kmuln(n1,n2,nr)
388: N n1,n2,*nr;
389: {
390: N n,t,s,m,carry;
391: int d,d1,d2,len,i,l;
392: int *r,*r0;
393:
394: if ( !n1 || !n2 ) {
395: *nr = 0; return;
396: }
397: d1 = PL(n1); d2 = PL(n2);
398: if ( (d1 < kara_mag) || (d2 < kara_mag) ) {
399: muln(n1,n2,nr); return;
400: }
401: if ( d1 < d2 ) {
402: n = n1; n1 = n2; n2 = n;
403: d = d1; d1 = d2; d2 = d;
404: }
405: if ( d2 > (d1+1)/2 ) {
406: kmulnmain(n1,n2,nr); return;
407: }
408: d = (d1/d2)+((d1%d2)!=0);
409: len = (d+1)*d2;
410: r0 = (int *)ALLOCA(len*sizeof(int));
411: bzero((char *)r0,len*sizeof(int));
412: for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) {
413: extractn(n1,i*d2,d2,&m);
414: if ( m ) {
415: kmuln(m,n2,&t);
416: addn(t,carry,&s);
417: copyn(s,d2,r);
418: extractn(s,d2,d2,&carry);
419: } else {
420: copyn(carry,d2,r);
421: carry = 0;
422: }
423: }
424: copyn(carry,d2,r);
425: for ( l = len - 1; !r0[l]; l-- );
426: l++;
427: *nr = t = NALLOC(l); PL(t) = l;
428: bcopy((char *)r0,(char *)BD(t),l*sizeof(int));
429: }
430:
431: void extractn(n,index,len,nr)
432: N n;
433: int index,len;
434: N *nr;
435: {
436: unsigned int *m;
437: int l;
438: N t;
439:
440: if ( !n ) {
441: *nr = 0; return;
442: }
443: m = BD(n)+index;
444: if ( (l = PL(n)-index) >= len ) {
445: for ( l = len - 1; (l >= 0) && !m[l]; l-- );
446: l++;
447: }
448: if ( l <= 0 ) {
449: *nr = 0; return;
450: } else {
451: *nr = t = NALLOC(l); PL(t) = l;
452: bcopy((char *)m,(char *)BD(t),l*sizeof(int));
453: }
454: }
455:
456: void copyn(n,len,p)
457: N n;
458: int len;
459: int *p;
460: {
461: if ( n )
462: bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int));
463: }
464:
465: void dupn(n,p)
466: N n;
467: N p;
468: {
469: if ( n )
470: bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int));
471: }
472:
473: void kmulnmain(n1,n2,nr)
474: N n1,n2,*nr;
475: {
476: int d1,d2,h,sgn,sgn1,sgn2,len;
477: N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
478:
479: d1 = PL(n1); d2 = PL(n2); h = (d1+1)/2;
480: extractn(n1,0,h,&n1lo); extractn(n1,h,d1-h,&n1hi);
481: extractn(n2,0,h,&n2lo); extractn(n2,h,d2-h,&n2hi);
482: kmuln(n1hi,n2hi,&hi); kmuln(n1lo,n2lo,&lo);
483: len = PL(hi)+2*h; t1 = NALLOC(len); PL(t1) = len;
484: bzero((char *)BD(t1),len*sizeof(int));
485: if ( lo )
486: bcopy((char *)BD(lo),(char *)BD(t1),PL(lo)*sizeof(int));
487: if ( hi )
488: bcopy((char *)BD(hi),(char *)(BD(t1)+2*h),PL(hi)*sizeof(int));
489:
490: addn(hi,lo,&mid1);
491: sgn1 = subn(n1hi,n1lo,&s1); sgn2 = subn(n2lo,n2hi,&s2);
492: kmuln(s1,s2,&mid2); sgn = sgn1*sgn2;
493: if ( sgn > 0 )
494: addn(mid1,mid2,&mid);
495: else
496: subn(mid1,mid2,&mid);
497: if ( mid ) {
498: len = PL(mid)+h; t2 = NALLOC(len); PL(t2) = len;
499: bzero((char *)BD(t2),len*sizeof(int));
500: bcopy((char *)BD(mid),(char *)(BD(t2)+h),PL(mid)*sizeof(int));
501: addn(t1,t2,nr);
502: } else
503: *nr = t1;
504: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>