![[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.4 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.5 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.4 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.6 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/Q.c,v 1.5 2000/08/22 05:04:04 noro Exp $
1.4 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "base.h"
52: #include "inline.h"
53:
54: void addq(n1,n2,nr)
55: Q n1,n2,*nr;
56: {
57: N nm,dn,nm1,nm2,nm3,dn0,dn1,dn2,g,g1,g0,m;
58: int sgn;
59: unsigned t,t1;
60:
61: if ( !n1 )
62: *nr = n2;
63: else if ( !n2 )
64: *nr = n1;
65: else if ( INT(n1) )
66: if ( INT(n2) ) {
67: nm1 = NM(n1); nm2 = NM(n2);
68: if ( SGN(n1) == SGN(n2) ) {
69: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
70: t1 = BD(nm1)[0]; t = t1+BD(nm2)[0];
71: if ( t < t1 ) {
72: m = NALLOC(2); PL(m) = 2; BD(m)[0] = t; BD(m)[1] = 1;
73: } else
74: UTON(t,m);
75: } else
76: addn(NM(n1),NM(n2),&m);
77: NTOQ(m,SGN(n1),*nr);
78: } else {
79: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
80: t1 = BD(nm1)[0]; t = t1-BD(nm2)[0];
81: if ( !t )
82: sgn = 0;
83: else if ( t > t1 ) {
84: sgn = -1; t = -((int)t); UTON(t,m);
85: } else {
86: sgn = 1; UTON(t,m);
87: }
88: } else
89: sgn = subn(NM(n1),NM(n2),&m);
90: if ( sgn ) {
91: if ( SGN(n1) == sgn )
92: NTOQ(m,1,*nr);
93: else
94: NTOQ(m,-1,*nr);
95: } else
96: *nr = 0;
97: }
98: } else {
99: kmuln(NM(n1),DN(n2),&m);
100: if ( SGN(n1) == SGN(n2) ) {
101: sgn = SGN(n1); addn(m,NM(n2),&nm);
102: } else
103: sgn = SGN(n1)*subn(m,NM(n2),&nm);
104: if ( sgn ) {
105: dn = DN(n2); NDTOQ(nm,dn,sgn,*nr);
106: } else
107: *nr = 0;
108: }
109: else if ( INT(n2) ) {
110: kmuln(NM(n2),DN(n1),&m);
111: if ( SGN(n1) == SGN(n2) ) {
112: sgn = SGN(n1); addn(m,NM(n1),&nm);
113: } else
114: sgn = SGN(n1)*subn(NM(n1),m,&nm);
115: if ( sgn ) {
116: dn = DN(n1); NDTOQ(nm,dn,sgn,*nr);
117: } else
118: *nr = 0;
119: } else {
120: gcdn(DN(n1),DN(n2),&g); divsn(DN(n1),g,&dn1); divsn(DN(n2),g,&dn2);
121: kmuln(NM(n1),dn2,&nm1); kmuln(NM(n2),dn1,&nm2);
122: if ( SGN(n1) == SGN(n2) ) {
123: sgn = SGN(n1); addn(nm1,nm2,&nm3);
124: } else
125: sgn = SGN(n1)*subn(nm1,nm2,&nm3);
126: if ( sgn ) {
127: gcdn(nm3,g,&g1); divsn(nm3,g1,&nm); divsn(g,g1,&g0);
128: kmuln(dn1,dn2,&dn0); kmuln(g0,dn0,&dn);
129: if ( UNIN(dn) )
130: NTOQ(nm,sgn,*nr);
131: else
132: NDTOQ(nm,dn,sgn,*nr);
133: } else
134: *nr = 0;
135: }
136: }
137:
138: void subq(n1,n2,nr)
139: Q n1,n2,*nr;
140: {
141: Q m;
142:
143: if ( !n1 )
144: if ( !n2 )
145: *nr = 0;
146: else {
147: DUPQ(n2,*nr); SGN(*nr) = -SGN(n2);
148: }
149: else if ( !n2 )
150: *nr = n1;
151: else if ( n1 == n2 )
152: *nr = 0;
153: else {
154: DUPQ(n2,m); SGN(m) = -SGN(n2); addq(n1,m,nr);
155: }
156: }
157:
158: void mulq(n1,n2,nr)
159: Q n1,n2,*nr;
160: {
161: N nm,nm1,nm2,dn,dn1,dn2,g;
162: int sgn;
163: unsigned u,l;
164:
165: if ( !n1 || !n2 ) *nr = 0;
166: else if ( INT(n1) )
167: if ( INT(n2) ) {
168: nm1 = NM(n1); nm2 = NM(n2);
169: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
170: DM(BD(NM(n1))[0],BD(NM(n2))[0],u,l)
171: if ( u ) {
172: nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = l; BD(nm)[1] = u;
173: } else {
174: nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = l;
175: }
176: } else
177: kmuln(nm1,nm2,&nm);
178: sgn = (SGN(n1)==SGN(n2)?1:-1); NTOQ(nm,sgn,*nr);
179: } else {
180: gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn);
181: kmuln(nm1,NM(n2),&nm); sgn = SGN(n1)*SGN(n2);
182: if ( UNIN(dn) )
183: NTOQ(nm,sgn,*nr);
184: else
185: NDTOQ(nm,dn,sgn,*nr);
186: }
187: else if ( INT(n2) ) {
188: gcdn(NM(n2),DN(n1),&g); divsn(NM(n2),g,&nm2); divsn(DN(n1),g,&dn);
189: kmuln(nm2,NM(n1),&nm); sgn = SGN(n1)*SGN(n2);
190: if ( UNIN(dn) )
191: NTOQ(nm,sgn,*nr);
192: else
193: NDTOQ(nm,dn,sgn,*nr);
194: } else {
195: gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn2);
196: gcdn(DN(n1),NM(n2),&g); divsn(DN(n1),g,&dn1); divsn(NM(n2),g,&nm2);
197: kmuln(nm1,nm2,&nm); kmuln(dn1,dn2,&dn); sgn = SGN(n1) * SGN(n2);
198: if ( UNIN(dn) )
199: NTOQ(nm,sgn,*nr);
200: else
201: NDTOQ(nm,dn,sgn,*nr);
202: }
203: }
204:
205: void divq(n1,n2,nq)
206: Q n1,n2,*nq;
207: {
208: Q m;
209:
210: if ( !n2 ) {
211: error("division by 0");
212: *nq = 0;
213: return;
214: } else if ( !n1 )
215: *nq = 0;
216: else if ( n1 == n2 )
217: *nq = ONE;
218: else {
219: invq(n2,&m); mulq(n1,m,nq);
1.6 ! noro 220: }
! 221: }
! 222:
! 223: void divsq(n1,n2,nq)
! 224: Q n1,n2,*nq;
! 225: {
! 226: Q m;
! 227: int sgn;
! 228: N tn;
! 229:
! 230: if ( !n2 ) {
! 231: error("division by 0");
! 232: *nq = 0;
! 233: return;
! 234: } else if ( !n1 )
! 235: *nq = 0;
! 236: else if ( n1 == n2 )
! 237: *nq = ONE;
! 238: else {
! 239: divsn(NM(n1),NM(n2),&tn);
! 240: sgn = SGN(n1)*SGN(n2);
! 241: NTOQ(tn,sgn,*nq);
1.1 noro 242: }
243: }
244:
245: void invq(n,nr)
246: Q n,*nr;
247: {
248: N nm,dn;
249:
250: if ( INT(n) )
251: if ( UNIN(NM(n)) )
252: if ( SGN(n) > 0 )
253: *nr = ONE;
254: else {
255: nm = ONEN; NTOQ(nm,SGN(n),*nr);
256: }
257: else {
258: nm = ONEN; dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
259: }
260: else if ( UNIN(NM(n)) ) {
261: nm = DN(n); NTOQ(nm,SGN(n),*nr);
262: } else {
263: nm = DN(n); dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
264: }
265: }
266:
267: void chsgnq(n,nr)
268: Q n,*nr;
269: {
270: Q t;
271:
272: if ( !n )
273: *nr = 0;
274: else {
275: DUPQ(n,t); SGN(t) = -SGN(t); *nr = t;
276: }
277: }
278:
279: void pwrq(n1,n,nr)
280: Q n1,n,*nr;
281: {
282: N nm,dn;
283: int sgn;
284:
285: if ( !n || UNIQ(n1) )
286: *nr = ONE;
287: else if ( !n1 )
288: *nr = 0;
289: else if ( !INT(n) ) {
290: error("can't calculate fractional power."); *nr = 0;
291: } else if ( MUNIQ(n1) )
292: *nr = BD(NM(n))[0]%2 ? n1 : ONE;
293: else if ( PL(NM(n)) > 1 ) {
294: error("exponent too big."); *nr = 0;
295: } else {
296: sgn = ((SGN(n1)>0)||EVENN(NM(n))?1:-1);
297: pwrn(NM(n1),BD(NM(n))[0],&nm);
298: if ( INT(n1) ) {
299: if ( UNIN(nm) )
300: if ( sgn == 1 )
301: *nr = ONE;
302: else
303: NTOQ(ONEN,sgn,*nr);
304: else if ( SGN(n) > 0 )
305: NTOQ(nm,sgn,*nr);
306: else
307: NDTOQ(ONEN,nm,sgn,*nr);
308: } else {
309: pwrn(DN(n1),BD(NM(n))[0],&dn);
310: if ( SGN(n) > 0 )
311: NDTOQ(nm,dn,sgn,*nr);
312: else if ( UNIN(nm) )
313: NTOQ(dn,sgn,*nr);
314: else
315: NDTOQ(dn,nm,sgn,*nr);
316: }
317: }
318: }
319:
320: int cmpq(q1,q2)
321: Q q1,q2;
322: {
323: int sgn;
324: N t,s;
325:
326: if ( !q1 )
327: if ( !q2 )
328: return ( 0 );
329: else
330: return ( -SGN(q2) );
331: else if ( !q2 )
332: return ( SGN(q1) );
333: else if ( SGN(q1) != SGN(q2) )
334: return ( SGN(q1) );
335: else if ( INT(q1) )
336: if ( INT(q2) ) {
337: sgn = cmpn(NM(q1),NM(q2));
338: if ( !sgn )
339: return ( 0 );
340: else
341: return ( SGN(q1)==sgn?1:-1 );
342: } else {
343: kmuln(NM(q1),DN(q2),&t); sgn = cmpn(t,NM(q2));
344: return ( SGN(q1) * sgn );
345: }
346: else if ( INT(q2) ) {
347: kmuln(NM(q2),DN(q1),&t); sgn = cmpn(NM(q1),t);
348: return ( SGN(q1) * sgn );
349: } else {
350: kmuln(NM(q1),DN(q2),&t); kmuln(NM(q2),DN(q1),&s); sgn = cmpn(t,s);
351: return ( SGN(q1) * sgn );
352: }
353: }
354:
355: void remq(n,m,nr)
356: Q n,m;
357: Q *nr;
358: {
359: N q,r,s;
360:
361: if ( !n ) {
362: *nr = 0;
363: return;
364: }
365: divn(NM(n),NM(m),&q,&r);
366: if ( !r )
367: *nr = 0;
368: else if ( SGN(n) > 0 )
369: NTOQ(r,1,*nr);
370: else {
371: subn(NM(m),r,&s); NTOQ(s,1,*nr);
372: }
373: }
374:
1.2 noro 375: /* t = [nC0 nC1 ... nCn] */
376:
1.1 noro 377: void mkbc(n,t)
378: int n;
379: Q *t;
380: {
381: int i;
382: N c,d;
383:
384: for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
385: c = NALLOC(1); PL(c) = 1; BD(c)[0] = n-i+1;
386: kmuln(NM(t[i-1]),c,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
387: }
388: for ( ; i <= n; i++ )
389: t[i] = t[n-i];
1.2 noro 390: }
391:
392: /*
393: * Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+...
394: *
395: * t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)]
396: * where W(k,l,i) = i! * kCi * lCi
397: */
398:
399: void mkwc(k,l,t)
400: int k,l;
401: Q *t;
402: {
403: int i,n,up,low;
404: N nm,d,c;
405:
406: n = MIN(k,l);
407: for ( t[0] = ONE, i = 1; i <= n; i++ ) {
408: DM(k-i+1,l-i+1,up,low);
409: if ( up ) {
410: nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = low; BD(nm)[1] = up;
411: } else {
412: nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = low;
413: }
414: kmuln(NM(t[i-1]),nm,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
1.3 noro 415: }
416: }
417:
418: /* mod m table */
419: /* XXX : should be optimized */
420:
421: void mkwcm(k,l,m,t)
422: int k,l,m;
423: int *t;
424: {
425: int i,n;
426: Q *s;
427:
428: n = MIN(k,l);
429: s = (Q *)ALLOCA((n+1)*sizeof(Q));
430: mkwc(k,l,s);
431: for ( i = 0; i <= n; i++ ) {
432: t[i] = rem(NM(s[i]),m);
1.2 noro 433: }
1.1 noro 434: }
435:
436: void factorial(n,r)
437: int n;
438: Q *r;
439: {
440: Q t,iq,s;
441: unsigned int i,m,m0;
442: unsigned int c;
443:
444: if ( !n )
445: *r = ONE;
446: else if ( n < 0 )
447: *r = 0;
448: else {
449: for ( i = 1, t = ONE; (int)i <= n; ) {
450: for ( m0 = m = 1, c = 0; ((int)i <= n) && !c; i++ ) {
451: m0 = m; DM(m0,i,c,m)
452: }
453: if ( c ) {
454: m = m0; i--;
455: }
456: UTOQ(m,iq); mulq(t,iq,&s); t = s;
457: }
458: *r = t;
459: }
460: }
461:
462: void invl(a,mod,ar)
463: Q a,mod,*ar;
464: {
465: Q f1,f2,a1,a2,q,m,s;
466: N qn,rn;
467:
468: a1 = ONE; a2 = 0;
469: DUPQ(a,f1); SGN(f1) = 1; f2 = mod;
470:
471: while ( 1 ) {
472: divn(NM(f1),NM(f2),&qn,&rn);
473: if ( !qn )
474: q = 0;
475: else
476: NTOQ(qn,1,q);
477:
478: f1 = f2;
479: if ( !rn )
480: break;
481: else
482: NTOQ(rn,1,f2);
483:
484: mulq(a2,q,&m); subq(a1,m,&s); a1 = a2;
485: if ( !s )
486: a2 = 0;
487: else
488: remq(s,mod,&a2);
489: }
490: if ( SGN(a) < 0 )
491: chsgnq(a2,&m);
492: else
493: m = a2;
494:
495: if ( SGN(m) >= 0 )
496: *ar = m;
497: else
498: addq(m,mod,ar);
499: }
500:
501: int kara_mag=100;
502:
503: void kmuln(n1,n2,nr)
504: N n1,n2,*nr;
505: {
506: N n,t,s,m,carry;
507: int d,d1,d2,len,i,l;
508: int *r,*r0;
509:
510: if ( !n1 || !n2 ) {
511: *nr = 0; return;
512: }
513: d1 = PL(n1); d2 = PL(n2);
514: if ( (d1 < kara_mag) || (d2 < kara_mag) ) {
515: muln(n1,n2,nr); return;
516: }
517: if ( d1 < d2 ) {
518: n = n1; n1 = n2; n2 = n;
519: d = d1; d1 = d2; d2 = d;
520: }
521: if ( d2 > (d1+1)/2 ) {
522: kmulnmain(n1,n2,nr); return;
523: }
524: d = (d1/d2)+((d1%d2)!=0);
525: len = (d+1)*d2;
526: r0 = (int *)ALLOCA(len*sizeof(int));
527: bzero((char *)r0,len*sizeof(int));
528: for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) {
529: extractn(n1,i*d2,d2,&m);
530: if ( m ) {
531: kmuln(m,n2,&t);
532: addn(t,carry,&s);
533: copyn(s,d2,r);
534: extractn(s,d2,d2,&carry);
535: } else {
536: copyn(carry,d2,r);
537: carry = 0;
538: }
539: }
540: copyn(carry,d2,r);
541: for ( l = len - 1; !r0[l]; l-- );
542: l++;
543: *nr = t = NALLOC(l); PL(t) = l;
544: bcopy((char *)r0,(char *)BD(t),l*sizeof(int));
545: }
546:
547: void extractn(n,index,len,nr)
548: N n;
549: int index,len;
550: N *nr;
551: {
552: unsigned int *m;
553: int l;
554: N t;
555:
556: if ( !n ) {
557: *nr = 0; return;
558: }
559: m = BD(n)+index;
560: if ( (l = PL(n)-index) >= len ) {
561: for ( l = len - 1; (l >= 0) && !m[l]; l-- );
562: l++;
563: }
564: if ( l <= 0 ) {
565: *nr = 0; return;
566: } else {
567: *nr = t = NALLOC(l); PL(t) = l;
568: bcopy((char *)m,(char *)BD(t),l*sizeof(int));
569: }
570: }
571:
572: void copyn(n,len,p)
573: N n;
574: int len;
575: int *p;
576: {
577: if ( n )
578: bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int));
579: }
580:
581: void dupn(n,p)
582: N n;
583: N p;
584: {
585: if ( n )
586: bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int));
587: }
588:
589: void kmulnmain(n1,n2,nr)
590: N n1,n2,*nr;
591: {
592: int d1,d2,h,sgn,sgn1,sgn2,len;
593: N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
594:
595: d1 = PL(n1); d2 = PL(n2); h = (d1+1)/2;
596: extractn(n1,0,h,&n1lo); extractn(n1,h,d1-h,&n1hi);
597: extractn(n2,0,h,&n2lo); extractn(n2,h,d2-h,&n2hi);
598: kmuln(n1hi,n2hi,&hi); kmuln(n1lo,n2lo,&lo);
599: len = PL(hi)+2*h; t1 = NALLOC(len); PL(t1) = len;
600: bzero((char *)BD(t1),len*sizeof(int));
601: if ( lo )
602: bcopy((char *)BD(lo),(char *)BD(t1),PL(lo)*sizeof(int));
603: if ( hi )
604: bcopy((char *)BD(hi),(char *)(BD(t1)+2*h),PL(hi)*sizeof(int));
605:
606: addn(hi,lo,&mid1);
607: sgn1 = subn(n1hi,n1lo,&s1); sgn2 = subn(n2lo,n2hi,&s2);
608: kmuln(s1,s2,&mid2); sgn = sgn1*sgn2;
609: if ( sgn > 0 )
610: addn(mid1,mid2,&mid);
611: else
612: subn(mid1,mid2,&mid);
613: if ( mid ) {
614: len = PL(mid)+h; t2 = NALLOC(len); PL(t2) = len;
615: bzero((char *)BD(t2),len*sizeof(int));
616: bcopy((char *)BD(mid),(char *)(BD(t2)+h),PL(mid)*sizeof(int));
617: addn(t1,t2,nr);
618: } else
619: *nr = t1;
620: }