![[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.3 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/Q.c,v 1.2 2000/04/05 08:32:17 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]);
1.3 ! noro 345: }
! 346: }
! 347:
! 348: /* mod m table */
! 349: /* XXX : should be optimized */
! 350:
! 351: void mkwcm(k,l,m,t)
! 352: int k,l,m;
! 353: int *t;
! 354: {
! 355: int i,n;
! 356: Q *s;
! 357:
! 358: n = MIN(k,l);
! 359: s = (Q *)ALLOCA((n+1)*sizeof(Q));
! 360: mkwc(k,l,s);
! 361: for ( i = 0; i <= n; i++ ) {
! 362: t[i] = rem(NM(s[i]),m);
1.2 noro 363: }
1.1 noro 364: }
365:
366: void factorial(n,r)
367: int n;
368: Q *r;
369: {
370: Q t,iq,s;
371: unsigned int i,m,m0;
372: unsigned int c;
373:
374: if ( !n )
375: *r = ONE;
376: else if ( n < 0 )
377: *r = 0;
378: else {
379: for ( i = 1, t = ONE; (int)i <= n; ) {
380: for ( m0 = m = 1, c = 0; ((int)i <= n) && !c; i++ ) {
381: m0 = m; DM(m0,i,c,m)
382: }
383: if ( c ) {
384: m = m0; i--;
385: }
386: UTOQ(m,iq); mulq(t,iq,&s); t = s;
387: }
388: *r = t;
389: }
390: }
391:
392: void invl(a,mod,ar)
393: Q a,mod,*ar;
394: {
395: Q f1,f2,a1,a2,q,m,s;
396: N qn,rn;
397:
398: a1 = ONE; a2 = 0;
399: DUPQ(a,f1); SGN(f1) = 1; f2 = mod;
400:
401: while ( 1 ) {
402: divn(NM(f1),NM(f2),&qn,&rn);
403: if ( !qn )
404: q = 0;
405: else
406: NTOQ(qn,1,q);
407:
408: f1 = f2;
409: if ( !rn )
410: break;
411: else
412: NTOQ(rn,1,f2);
413:
414: mulq(a2,q,&m); subq(a1,m,&s); a1 = a2;
415: if ( !s )
416: a2 = 0;
417: else
418: remq(s,mod,&a2);
419: }
420: if ( SGN(a) < 0 )
421: chsgnq(a2,&m);
422: else
423: m = a2;
424:
425: if ( SGN(m) >= 0 )
426: *ar = m;
427: else
428: addq(m,mod,ar);
429: }
430:
431: int kara_mag=100;
432:
433: void kmuln(n1,n2,nr)
434: N n1,n2,*nr;
435: {
436: N n,t,s,m,carry;
437: int d,d1,d2,len,i,l;
438: int *r,*r0;
439:
440: if ( !n1 || !n2 ) {
441: *nr = 0; return;
442: }
443: d1 = PL(n1); d2 = PL(n2);
444: if ( (d1 < kara_mag) || (d2 < kara_mag) ) {
445: muln(n1,n2,nr); return;
446: }
447: if ( d1 < d2 ) {
448: n = n1; n1 = n2; n2 = n;
449: d = d1; d1 = d2; d2 = d;
450: }
451: if ( d2 > (d1+1)/2 ) {
452: kmulnmain(n1,n2,nr); return;
453: }
454: d = (d1/d2)+((d1%d2)!=0);
455: len = (d+1)*d2;
456: r0 = (int *)ALLOCA(len*sizeof(int));
457: bzero((char *)r0,len*sizeof(int));
458: for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) {
459: extractn(n1,i*d2,d2,&m);
460: if ( m ) {
461: kmuln(m,n2,&t);
462: addn(t,carry,&s);
463: copyn(s,d2,r);
464: extractn(s,d2,d2,&carry);
465: } else {
466: copyn(carry,d2,r);
467: carry = 0;
468: }
469: }
470: copyn(carry,d2,r);
471: for ( l = len - 1; !r0[l]; l-- );
472: l++;
473: *nr = t = NALLOC(l); PL(t) = l;
474: bcopy((char *)r0,(char *)BD(t),l*sizeof(int));
475: }
476:
477: void extractn(n,index,len,nr)
478: N n;
479: int index,len;
480: N *nr;
481: {
482: unsigned int *m;
483: int l;
484: N t;
485:
486: if ( !n ) {
487: *nr = 0; return;
488: }
489: m = BD(n)+index;
490: if ( (l = PL(n)-index) >= len ) {
491: for ( l = len - 1; (l >= 0) && !m[l]; l-- );
492: l++;
493: }
494: if ( l <= 0 ) {
495: *nr = 0; return;
496: } else {
497: *nr = t = NALLOC(l); PL(t) = l;
498: bcopy((char *)m,(char *)BD(t),l*sizeof(int));
499: }
500: }
501:
502: void copyn(n,len,p)
503: N n;
504: int len;
505: int *p;
506: {
507: if ( n )
508: bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int));
509: }
510:
511: void dupn(n,p)
512: N n;
513: N p;
514: {
515: if ( n )
516: bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int));
517: }
518:
519: void kmulnmain(n1,n2,nr)
520: N n1,n2,*nr;
521: {
522: int d1,d2,h,sgn,sgn1,sgn2,len;
523: N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
524:
525: d1 = PL(n1); d2 = PL(n2); h = (d1+1)/2;
526: extractn(n1,0,h,&n1lo); extractn(n1,h,d1-h,&n1hi);
527: extractn(n2,0,h,&n2lo); extractn(n2,h,d2-h,&n2hi);
528: kmuln(n1hi,n2hi,&hi); kmuln(n1lo,n2lo,&lo);
529: len = PL(hi)+2*h; t1 = NALLOC(len); PL(t1) = len;
530: bzero((char *)BD(t1),len*sizeof(int));
531: if ( lo )
532: bcopy((char *)BD(lo),(char *)BD(t1),PL(lo)*sizeof(int));
533: if ( hi )
534: bcopy((char *)BD(hi),(char *)(BD(t1)+2*h),PL(hi)*sizeof(int));
535:
536: addn(hi,lo,&mid1);
537: sgn1 = subn(n1hi,n1lo,&s1); sgn2 = subn(n2lo,n2hi,&s2);
538: kmuln(s1,s2,&mid2); sgn = sgn1*sgn2;
539: if ( sgn > 0 )
540: addn(mid1,mid2,&mid);
541: else
542: subn(mid1,mid2,&mid);
543: if ( mid ) {
544: len = PL(mid)+h; t2 = NALLOC(len); PL(t2) = len;
545: bzero((char *)BD(t2),len*sizeof(int));
546: bcopy((char *)BD(mid),(char *)(BD(t2)+h),PL(mid)*sizeof(int));
547: addn(t1,t2,nr);
548: } else
549: *nr = t1;
550: }