![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Q.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to Q.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.2 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/Q.c,v 1.1.1.1 1999/12/03 07:39:08 noro Exp $ */
1.1 noro 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:
1.2 ! noro 305: /* t = [nC0 nC1 ... nCn] */
! 306:
1.1 noro 307: void mkbc(n,t)
308: int n;
309: Q *t;
310: {
311: int i;
312: N c,d;
313:
314: for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
315: c = NALLOC(1); PL(c) = 1; BD(c)[0] = n-i+1;
316: kmuln(NM(t[i-1]),c,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
317: }
318: for ( ; i <= n; i++ )
319: t[i] = t[n-i];
1.2 ! noro 320: }
! 321:
! 322: /*
! 323: * Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+...
! 324: *
! 325: * t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)]
! 326: * where W(k,l,i) = i! * kCi * lCi
! 327: */
! 328:
! 329: void mkwc(k,l,t)
! 330: int k,l;
! 331: Q *t;
! 332: {
! 333: int i,n,up,low;
! 334: N nm,d,c;
! 335:
! 336: n = MIN(k,l);
! 337: for ( t[0] = ONE, i = 1; i <= n; i++ ) {
! 338: DM(k-i+1,l-i+1,up,low);
! 339: if ( up ) {
! 340: nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = low; BD(nm)[1] = up;
! 341: } else {
! 342: nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = low;
! 343: }
! 344: kmuln(NM(t[i-1]),nm,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
! 345: }
1.1 noro 346: }
347:
348: void factorial(n,r)
349: int n;
350: Q *r;
351: {
352: Q t,iq,s;
353: unsigned int i,m,m0;
354: unsigned int c;
355:
356: if ( !n )
357: *r = ONE;
358: else if ( n < 0 )
359: *r = 0;
360: else {
361: for ( i = 1, t = ONE; (int)i <= n; ) {
362: for ( m0 = m = 1, c = 0; ((int)i <= n) && !c; i++ ) {
363: m0 = m; DM(m0,i,c,m)
364: }
365: if ( c ) {
366: m = m0; i--;
367: }
368: UTOQ(m,iq); mulq(t,iq,&s); t = s;
369: }
370: *r = t;
371: }
372: }
373:
374: void invl(a,mod,ar)
375: Q a,mod,*ar;
376: {
377: Q f1,f2,a1,a2,q,m,s;
378: N qn,rn;
379:
380: a1 = ONE; a2 = 0;
381: DUPQ(a,f1); SGN(f1) = 1; f2 = mod;
382:
383: while ( 1 ) {
384: divn(NM(f1),NM(f2),&qn,&rn);
385: if ( !qn )
386: q = 0;
387: else
388: NTOQ(qn,1,q);
389:
390: f1 = f2;
391: if ( !rn )
392: break;
393: else
394: NTOQ(rn,1,f2);
395:
396: mulq(a2,q,&m); subq(a1,m,&s); a1 = a2;
397: if ( !s )
398: a2 = 0;
399: else
400: remq(s,mod,&a2);
401: }
402: if ( SGN(a) < 0 )
403: chsgnq(a2,&m);
404: else
405: m = a2;
406:
407: if ( SGN(m) >= 0 )
408: *ar = m;
409: else
410: addq(m,mod,ar);
411: }
412:
413: int kara_mag=100;
414:
415: void kmuln(n1,n2,nr)
416: N n1,n2,*nr;
417: {
418: N n,t,s,m,carry;
419: int d,d1,d2,len,i,l;
420: int *r,*r0;
421:
422: if ( !n1 || !n2 ) {
423: *nr = 0; return;
424: }
425: d1 = PL(n1); d2 = PL(n2);
426: if ( (d1 < kara_mag) || (d2 < kara_mag) ) {
427: muln(n1,n2,nr); return;
428: }
429: if ( d1 < d2 ) {
430: n = n1; n1 = n2; n2 = n;
431: d = d1; d1 = d2; d2 = d;
432: }
433: if ( d2 > (d1+1)/2 ) {
434: kmulnmain(n1,n2,nr); return;
435: }
436: d = (d1/d2)+((d1%d2)!=0);
437: len = (d+1)*d2;
438: r0 = (int *)ALLOCA(len*sizeof(int));
439: bzero((char *)r0,len*sizeof(int));
440: for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) {
441: extractn(n1,i*d2,d2,&m);
442: if ( m ) {
443: kmuln(m,n2,&t);
444: addn(t,carry,&s);
445: copyn(s,d2,r);
446: extractn(s,d2,d2,&carry);
447: } else {
448: copyn(carry,d2,r);
449: carry = 0;
450: }
451: }
452: copyn(carry,d2,r);
453: for ( l = len - 1; !r0[l]; l-- );
454: l++;
455: *nr = t = NALLOC(l); PL(t) = l;
456: bcopy((char *)r0,(char *)BD(t),l*sizeof(int));
457: }
458:
459: void extractn(n,index,len,nr)
460: N n;
461: int index,len;
462: N *nr;
463: {
464: unsigned int *m;
465: int l;
466: N t;
467:
468: if ( !n ) {
469: *nr = 0; return;
470: }
471: m = BD(n)+index;
472: if ( (l = PL(n)-index) >= len ) {
473: for ( l = len - 1; (l >= 0) && !m[l]; l-- );
474: l++;
475: }
476: if ( l <= 0 ) {
477: *nr = 0; return;
478: } else {
479: *nr = t = NALLOC(l); PL(t) = l;
480: bcopy((char *)m,(char *)BD(t),l*sizeof(int));
481: }
482: }
483:
484: void copyn(n,len,p)
485: N n;
486: int len;
487: int *p;
488: {
489: if ( n )
490: bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int));
491: }
492:
493: void dupn(n,p)
494: N n;
495: N p;
496: {
497: if ( n )
498: bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int));
499: }
500:
501: void kmulnmain(n1,n2,nr)
502: N n1,n2,*nr;
503: {
504: int d1,d2,h,sgn,sgn1,sgn2,len;
505: N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
506:
507: d1 = PL(n1); d2 = PL(n2); h = (d1+1)/2;
508: extractn(n1,0,h,&n1lo); extractn(n1,h,d1-h,&n1hi);
509: extractn(n2,0,h,&n2lo); extractn(n2,h,d2-h,&n2hi);
510: kmuln(n1hi,n2hi,&hi); kmuln(n1lo,n2lo,&lo);
511: len = PL(hi)+2*h; t1 = NALLOC(len); PL(t1) = len;
512: bzero((char *)BD(t1),len*sizeof(int));
513: if ( lo )
514: bcopy((char *)BD(lo),(char *)BD(t1),PL(lo)*sizeof(int));
515: if ( hi )
516: bcopy((char *)BD(hi),(char *)(BD(t1)+2*h),PL(hi)*sizeof(int));
517:
518: addn(hi,lo,&mid1);
519: sgn1 = subn(n1hi,n1lo,&s1); sgn2 = subn(n2lo,n2hi,&s2);
520: kmuln(s1,s2,&mid2); sgn = sgn1*sgn2;
521: if ( sgn > 0 )
522: addn(mid1,mid2,&mid);
523: else
524: subn(mid1,mid2,&mid);
525: if ( mid ) {
526: len = PL(mid)+h; t2 = NALLOC(len); PL(t2) = len;
527: bzero((char *)BD(t2),len*sizeof(int));
528: bcopy((char *)BD(mid),(char *)(BD(t2)+h),PL(mid)*sizeof(int));
529: addn(t1,t2,nr);
530: } else
531: *nr = t1;
532: }