![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Z.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 Z.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: #include "ca.h"
1.5 noro 2: #include "base.h"
1.1 noro 3: #include "inline.h"
4:
1.7 fujiwara 5: #if defined(__GNUC__)
1.12 ohara 6: #define INLINE static inline
1.14 ! fujimoto 7: #elif defined(VISUAL) || defined(__MINGW32__)
1.7 fujiwara 8: #define INLINE __inline
9: #else
10: #define INLINE
11: #endif
12:
13: INLINE void _addz(Z n1,Z n2,Z nr);
14: INLINE void _subz(Z n1,Z n2,Z nr);
15: INLINE void _mulz(Z n1,Z n2,Z nr);
1.8 noro 16: int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr);
17: int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr);
1.3 noro 18:
19: /* immediate int -> Z */
1.5 noro 20: #define UTOZ(c,n) (n)=(!((unsigned int)(c))?0:(((unsigned int)(c))<=IMM_MAX?((Z)((((unsigned int)(c))<<1)|1)):utoz((unsigned int)(c))))
21: #define STOZ(c,n) (n)=(!((int)(c))?0:(((int)(c))>=IMM_MIN&&((int)(c))<=IMM_MAX?((Z)((((int)(c))<<1)|1)):stoz((int)(c))))
1.3 noro 22: /* immediate Z ? */
23: #define IS_IMM(c) (((unsigned int)c)&1)
24: /* Z can be conver to immediate ? */
1.5 noro 25: #define IS_SZ(n) (((SL(n) == 1)||(SL(n)==-1))&&BD(n)[0]<=IMM_MAX)
1.3 noro 26: /* Z -> immediate Z */
1.5 noro 27: #define SZTOZ(n,z) (z)=(Z)(SL(n)<0?(((-BD(n)[0])<<1)|1):(((BD(n)[0])<<1)|1))
1.3 noro 28: /* Z -> immediate int */
1.5 noro 29: #define ZTOS(c) (((int)(c))>>1)
30:
31: int uniz(Z a)
32: {
33: if ( IS_IMM(a) && ZTOS(a) == 1 ) return 1;
34: else return 0;
35: }
36:
37: int cmpz(Z a,Z b)
38: {
39: int *ma,*mb;
40: int sa,sb,da,db,ca,cb,i;
41:
42: if ( !a )
43: return -sgnz(b);
44: else if ( !b )
45: return sgnz(a);
46: else {
47: sa = sgnz(a); sb = sgnz(b);
48: if ( sa > sb ) return 1;
49: else if ( sa < sb ) return -1;
50: else if ( IS_IMM(a) )
51: if ( IS_IMM(b) ) {
52: ca = ZTOS(a); cb = ZTOS(b);
53: if ( ca > cb ) return sa;
54: else if ( ca < cb ) return -sa;
55: else return 0;
56: } else
57: return -sa;
58: else if ( IS_IMM(b) )
59: return sa;
60: else {
61: da = SL(a)*sa; db = SL(b)*sa;
62: if ( da > db ) return sa;
63: else if ( da < db ) return -sa;
64: else {
65: for ( i = da-1, ma = BD(a)+i, mb = BD(b)+i; i >= 0; i-- )
66: if ( *ma > *mb ) return sa;
67: else if ( *ma < *mb ) return -sa;
68: return 0;
69: }
70: }
71: }
72: }
1.1 noro 73:
1.5 noro 74: Z stoz(int c)
1.3 noro 75: {
76: Z z;
77:
78: z = ZALLOC(1);
79: if ( c < 0 ) {
80: SL(z) = -1; BD(z)[0] = -c;
81: } else {
82: SL(z) = 1; BD(z)[0] = c;
83: }
84: return z;
85: }
86:
1.5 noro 87: Z utoz(unsigned int c)
88: {
89: Z z;
90:
91: z = ZALLOC(1);
92: SL(z) = 1; BD(z)[0] = c;
93: return z;
94: }
95:
96: Z simpz(Z n)
97: {
98: Z n1;
99:
100: if ( !n ) return 0;
101: else if ( IS_IMM(n) ) return n;
102: else if ( IS_SZ(n) ) {
103: SZTOZ(n,n1); return n1;
104: } else
105: return n;
106: }
107:
1.3 noro 108: int sgnz(Z n)
1.2 noro 109: {
110: if ( !n ) return 0;
1.5 noro 111: else if ( IS_IMM(n) ) return ZTOS(n)>0?1:-1;
1.2 noro 112: else if ( SL(n) < 0 ) return -1;
113: else return 1;
114: }
115:
1.1 noro 116: z_mag(Z n)
117: {
1.4 noro 118: int c,i;
1.3 noro 119:
120: if ( !n ) return 0;
121: else if ( IS_IMM(n) ) {
1.5 noro 122: c = ZTOS(n);
1.3 noro 123: if ( c < 0 ) c = -c;
124: for ( i = 0; c; c >>= 1, i++ );
125: return i;
126: }
127: else return n_bits((N)n);
1.1 noro 128: }
129:
130: Z qtoz(Q n)
131: {
1.3 noro 132: Z r,t;
1.4 noro 133: int c;
1.1 noro 134:
135: if ( !n ) return 0;
136: else if ( !INT(n) )
137: error("qtoz : invalid input");
138: else {
1.3 noro 139: t = (Z)NM(n);
1.5 noro 140: if ( IS_SZ(t) ) {
1.4 noro 141: c = SGN(n) < 0 ? -BD(t)[0] : BD(t)[0];
1.5 noro 142: STOZ(c,r);
1.3 noro 143: } else {
144: r = dupz((Z)t);
1.4 noro 145: if ( SGN(n) < 0 ) SL(r) = -SL(r);
1.3 noro 146: }
1.1 noro 147: return r;
148: }
149: }
150:
151: Q ztoq(Z n)
152: {
153: Q r;
154: Z nm;
1.3 noro 155: int sgn,c;
1.1 noro 156:
157: if ( !n ) return 0;
1.3 noro 158: else if ( IS_IMM(n) ) {
1.5 noro 159: c = ZTOS(n);
1.3 noro 160: STOQ(c,r);
161: return r;
162: } else {
1.1 noro 163: nm = dupz(n);
164: if ( SL(nm) < 0 ) {
165: sgn = -1;
166: SL(nm) = -SL(nm);
167: } else
168: sgn = 1;
169: NTOQ((N)nm,sgn,r);
170: return r;
171: }
172: }
173:
174: Z dupz(Z n)
175: {
1.5 noro 176: Z r;
1.1 noro 177: int sd,i;
178:
179: if ( !n ) return 0;
1.3 noro 180: else if ( IS_IMM(n) ) return n;
181: else {
182: if ( (sd = SL(n)) < 0 ) sd = -sd;
1.5 noro 183: r = ZALLOC(sd);
184: SL(r) = SL(n);
185: for ( i = 0; i < sd; i++ ) BD(r)[i] = BD(n)[i];
186: return r;
1.3 noro 187: }
1.1 noro 188: }
189:
190: Z chsgnz(Z n)
191: {
192: Z r;
1.3 noro 193: int c;
1.1 noro 194:
195: if ( !n ) return 0;
1.3 noro 196: else if ( IS_IMM(n) ) {
1.5 noro 197: c = -ZTOS(n);
198: STOZ(c,r);
1.3 noro 199: return r;
200: } else {
1.1 noro 201: r = dupz(n);
202: SL(r) = -SL(r);
203: return r;
204: }
205: }
206:
1.5 noro 207: Z absz(Z n)
208: {
209: Z r;
210: int c;
211:
212: if ( !n ) return 0;
213: else if ( sgnz(n) > 0 )
214: return n;
215: else
216: return chsgnz(n);
217: }
1.3 noro 218:
1.1 noro 219: Z addz(Z n1,Z n2)
220: {
1.3 noro 221: int d1,d2,d,c;
222: Z r,r1;
223: struct oZ t;
1.1 noro 224:
225: if ( !n1 ) return dupz(n2);
226: else if ( !n2 ) return dupz(n1);
1.3 noro 227: else if ( IS_IMM(n1) ) {
228: if ( IS_IMM(n2) ) {
1.5 noro 229: c = ZTOS(n1)+ZTOS(n2);
230: STOZ(c,r);
1.3 noro 231: } else {
1.5 noro 232: c = ZTOS(n1);
1.3 noro 233: if ( c < 0 ) {
234: t.p = -1; t.b[0] = -c;
1.1 noro 235: } else {
1.3 noro 236: t.p = 1; t.b[0] = c;
1.1 noro 237: }
1.3 noro 238: if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
239: r = ZALLOC(d2+1);
240: _addz(&t,n2,r);
1.5 noro 241: if ( !SL(r) ) r = 0;
1.3 noro 242: }
243: } else if ( IS_IMM(n2) ) {
1.5 noro 244: c = ZTOS(n2);
1.3 noro 245: if ( c < 0 ) {
246: t.p = -1; t.b[0] = -c;
1.1 noro 247: } else {
1.3 noro 248: t.p = 1; t.b[0] = c;
249: }
250: if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
251: r = ZALLOC(d1+1);
252: _addz(n1,&t,r);
1.5 noro 253: if ( !SL(r) ) r = 0;
1.3 noro 254: } else {
255: if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
256: if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
257: d = MAX(d1,d2)+1;
258: r = ZALLOC(d);
259: _addz(n1,n2,r);
260: if ( !SL(r) ) r = 0;
1.1 noro 261: }
1.5 noro 262: if ( r && !((int)r&1) && IS_SZ(r) ) {
263: SZTOZ(r,r1); r = r1;
264: }
265: return r;
1.1 noro 266: }
267:
268: Z subz(Z n1,Z n2)
269: {
1.3 noro 270: int d1,d2,d,c;
271: Z r,r1;
272: struct oZ t;
1.1 noro 273:
1.5 noro 274: if ( !n1 )
275: return chsgnz(n2);
276: else if ( !n2 ) return n1;
1.3 noro 277: else if ( IS_IMM(n1) ) {
278: if ( IS_IMM(n2) ) {
1.5 noro 279: c = ZTOS(n1)-ZTOS(n2);
280: STOZ(c,r);
1.3 noro 281: } else {
1.5 noro 282: c = ZTOS(n1);
1.3 noro 283: if ( c < 0 ) {
284: t.p = -1; t.b[0] = -c;
285: } else {
286: t.p = 1; t.b[0] = c;
287: }
288: if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
289: r = ZALLOC(d2+1);
290: _subz(&t,n2,r);
1.5 noro 291: if ( !SL(r) ) r = 0;
1.3 noro 292: }
293: } else if ( IS_IMM(n2) ) {
1.5 noro 294: c = ZTOS(n2);
1.3 noro 295: if ( c < 0 ) {
296: t.p = -1; t.b[0] = -c;
297: } else {
298: t.p = 1; t.b[0] = c;
299: }
300: if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
301: r = ZALLOC(d1+1);
302: _subz(n1,&t,r);
1.5 noro 303: if ( !SL(r) ) r = 0;
1.3 noro 304: } else {
305: if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
306: if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
307: d = MAX(d1,d2)+1;
308: r = ZALLOC(d);
309: _subz(n1,n2,r);
310: if ( !SL(r) ) r = 0;
1.1 noro 311: }
1.5 noro 312: if ( r && !((int)r&1) && IS_SZ(r) ) {
313: SZTOZ(r,r1); r = r1;
314: }
315: return r;
1.1 noro 316: }
317:
318: Z mulz(Z n1,Z n2)
319: {
1.3 noro 320: int d1,d2,sgn,i;
1.5 noro 321: int c1,c2;
1.1 noro 322: unsigned int u1,u2,u,l;
1.3 noro 323: Z r;
1.1 noro 324:
325: if ( !n1 || !n2 ) return 0;
1.3 noro 326:
327: if ( IS_IMM(n1) ) {
1.5 noro 328: c1 = ZTOS(n1);
1.3 noro 329: if ( IS_IMM(n2) ) {
1.5 noro 330: c2 = ZTOS(n2);
331: if ( c1 == 1 )
332: return n2;
333: else if ( c1 == -1 )
334: return chsgnz(n2);
335: else if ( c2 == 1 )
336: return n1;
337: else if ( c2 == -1 )
338: return chsgnz(n1);
339: else {
340: sgn = 1;
341: if ( c1 < 0 ) { c1 = -c1; sgn = -sgn; }
342: if ( c2 < 0 ) { c2 = -c2; sgn = -sgn; }
343: u1 = (unsigned int)c1; u2 = (unsigned int)c2;
344: DM(u1,u2,u,l);
345: if ( !u ) {
346: UTOZ(l,r);
347: } else {
348: r = ZALLOC(2); SL(r) = 2; BD(r)[1] = u; BD(r)[0] = l;
349: }
1.1 noro 350: }
351: } else {
1.5 noro 352: sgn = 1;
353: if ( c1 < 0 ) { c1 = -c1; sgn = -sgn; }
354: u1 = (unsigned int)c1;
1.3 noro 355: if ( (d2 = SL(n2)) < 0 ) { sgn = -sgn; d2 = -d2; }
356: r = ZALLOC(d2+1);
357: for ( i = d2; i >= 0; i-- ) BD(r)[i] = 0;
358: muln_1(BD(n2),d2,u1,BD(r));
359: SL(r) = BD(r)[d2]?d2+1:d2;
1.1 noro 360: }
1.3 noro 361: } else if ( IS_IMM(n2) ) {
1.5 noro 362: c2 = ZTOS(n2);
363: if ( c2 == 1 )
364: return n1;
365: else if ( c2 == -1 )
366: return chsgnz(n1);
367:
1.3 noro 368: sgn = 1;
1.5 noro 369: if ( c2 < 0 ) { sgn = -sgn; c2 = -c2; }
370: u2 = (unsigned int)c2;
1.3 noro 371: if ( (d1 = SL(n1)) < 0 ) { sgn = -sgn; d1 = -d1; }
372: r = ZALLOC(d1+1);
373: for ( i = d1; i >= 0; i-- ) BD(r)[i] = 0;
374: muln_1(BD(n1),d1,u2,BD(r));
375: SL(r) = BD(r)[d1]?d1+1:d1;
376: } else {
377: sgn = 1;
1.5 noro 378: if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
379: if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
1.3 noro 380: r = ZALLOC(d1+d2);
381: _mulz(n1,n2,r);
1.1 noro 382: }
1.5 noro 383: if ( sgn < 0 ) r = chsgnz(r);
1.3 noro 384: return r;
1.1 noro 385: }
386:
1.3 noro 387: /* kokokara */
1.4 noro 388: #if 0
1.1 noro 389: Z divsz(Z n1,Z n2)
390: {
1.3 noro 391: int sgn,d1,d2;
392: Z q;
1.1 noro 393:
394: if ( !n2 ) error("divsz : division by 0");
1.3 noro 395: if ( !n1 ) return 0;
396:
397: if ( IS_IMM(n1) ) {
398: if ( !IS_IMM(n2) )
399: error("divsz : cannot happen");
1.5 noro 400: c = ZTOS(n1)/ZTOS(n2);
401: STOZ(c,q);
1.3 noro 402: return q;
1.1 noro 403: }
1.3 noro 404: if ( IS_IMM(n2) ) {
405: sgn = 1;
1.5 noro 406: u2 = ZTOS(n2); if ( u2 < 0 ) { sgn = -sgn; u2 = -u2; }
1.3 noro 407: diviz(n1,u2,&q);
408: if ( sgn < 0 ) SL(q) = -SL(q);
409: return q;
410: }
411:
412: sgn = 1;
413: if ( (d2 = SL(n2)) < 0 ) { sgn = -sgn; d2 = -d2; }
414: if ( d2 == 1 ) {
415: diviz(n1,BD(u2)[0],&q);
416: if ( sgn < 0 ) SL(q) = -SL(q);
417: return q;
418: }
419: if ( (d1 = SL(n1)) < 0 ) { sgn = -sgn; d1 = -d1; }
420: if ( d1 < d2 ) error("divsz : cannot happen");
421: return q;
1.1 noro 422: }
1.5 noro 423: #endif
424:
425: /* XXX */
426: Z divz(Z n1,Z n2,Z *r)
427: {
428: int s1,s2;
429: Q t1,t2,qq,rq;
430: N qn,rn;
431:
432: if ( !n1 ) {
433: *r = 0; return 0;
434: }
435: if ( !n2 )
436: error("divz : division by 0");
437: t1 = ztoq(n1); t2 = ztoq(n2);
438: s1 = sgnz(n1); s2 = sgnz(n2);
439: /* n1 = qn*SGN(n1)*SGN(n2)*n2+SGN(n1)*rn */
440: divn(NM(t1),NM(t2),&qn,&rn);
441: NTOQ(qn,s1*s2,qq);
442: NTOQ(rn,s1,rq);
443: *r = qtoz(rq);
444: return qtoz(qq);
445: }
446:
447: Z divsz(Z n1,Z n2)
448: {
449: int s1,s2;
450: Q t1,t2,qq;
451: N qn;
452:
453: if ( !n1 )
454: return 0;
455: if ( !n2 )
456: error("divsz : division by 0");
457: t1 = ztoq(n1); t2 = ztoq(n2);
458: s1 = sgnz(n1); s2 = sgnz(n2);
459: /* n1 = qn*SGN(n1)*SGN(n2)*n2 */
460: divsn(NM(t1),NM(t2),&qn);
461: NTOQ(qn,s1*s2,qq);
462: return qtoz(qq);
463: }
464:
465: int gcdimm(int c1,int c2)
466: {
467: int r;
468:
469: if ( !c1 ) return c2;
470: else if ( !c2 ) return c1;
471: while ( 1 ) {
472: r = c1%c2;
473: if ( !r ) return c2;
474: c1 = c2; c2 = r;
475: }
476: }
1.1 noro 477:
478: Z gcdz(Z n1,Z n2)
479: {
1.5 noro 480: int c1,c2,g;
481: Z gcd,r;
482: N gn;
1.1 noro 483:
484: if ( !n1 ) return n2;
485: else if ( !n2 ) return n1;
486:
1.5 noro 487: if ( IS_IMM(n1) ) {
488: c1 = ZTOS(n1);
489: if ( c1 < 0 ) c1 = -c1;
490: if ( IS_IMM(n2) )
491: c2 = ZTOS(n2);
492: else
493: c2 = remzi(n2,c1>0?c1:-c1);
494: if ( c2 < 0 ) c2 = -c2;
495: g = gcdimm(c1,c2);
496: STOZ(g,gcd);
497: return gcd;
498: } else if ( IS_IMM(n2) ) {
499: c2 = ZTOS(n2);
500: if ( c2 < 0 ) c2 = -c2;
501: c1 = remzi(n1,c2>0?c2:-c2);
502: if ( c1 < 0 ) c1 = -c1;
503: g = gcdimm(c1,c2);
504: STOZ(g,gcd);
505: return gcd;
506: } else {
507: n1 = dupz(n1);
508: if ( SL(n1) < 0 ) SL(n1) = -SL(n1);
509: n2 = dupz(n2);
510: if ( SL(n2) < 0 ) SL(n2) = -SL(n2);
511: gcdn((N)n1,(N)n2,&gn);
512: gcd = (Z)gn;
513: if ( IS_SZ(gcd) ) {
514: SZTOZ(gcd,r); gcd = r;
515: }
516: return gcd;
517: }
1.1 noro 518: }
519:
520: int remzi(Z n,int m)
521: {
522: unsigned int *x;
523: unsigned int t,r;
1.5 noro 524: int i,c;
1.1 noro 525:
526: if ( !n ) return 0;
1.5 noro 527: if ( IS_IMM(n) ) {
528: c = ZTOS(n)%m;
529: if ( c < 0 ) c += m;
530: return c;
531: }
532:
1.1 noro 533: i = SL(n);
534: if ( i < 0 ) i = -i;
535: for ( i--, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) {
536: #if defined(sparc)
537: r = dsa(m,r,*x);
538: #else
539: DSAB(m,r,*x,t,r);
540: #endif
541: }
1.5 noro 542: if ( r && SL(n) < 0 )
543: r = m-r;
1.1 noro 544: return r;
545: }
546:
547: Z gcdz_cofactor(Z n1,Z n2,Z *c1,Z *c2)
548: {
549: Z gcd;
550:
551: gcd = gcdz(n1,n2);
552: *c1 = divsz(n1,gcd);
553: *c2 = divsz(n2,gcd);
554: return gcd;
555: }
556:
1.5 noro 557: #if 0
1.1 noro 558: Z estimate_array_gcdz(Z *b,int n)
559: {
560: int m,i,j,sd;
561: Z *a;
562: Z s,t;
563:
564: a = (Z *)ALLOCA(n*sizeof(Z));
565: for ( i = j = 0; i < n; i++ ) if ( b[i] ) a[j++] = b[i];
566: n = j;
567: if ( !n ) return 0;
568: if ( n == 1 ) return a[0];
569:
570: m = n/2;
571: for ( i = 0, s = 0; i < m; i++ ) {
572: if ( !a[i] ) continue;
573: else s = (SL(a[i])<0)?subz(s,a[i]):addz(s,a[i]);
574: }
575: for ( t = 0; i < n; i++ ) {
576: if ( !a[i] ) continue;
577: else t = (SL(a[i])<0)?subz(t,a[i]):addz(t,a[i]);
578: }
579: return gcdz(s,t);
580: }
581:
582: Z array_gcdz(Z *b,int n)
583: {
584: int m,i,j,sd;
585: Z *a;
586: Z gcd;
587:
588: a = (Z *)ALLOCA(n*sizeof(Z));
589: for ( i = j = 0; i < n; i++ ) if ( b[i] ) a[j++] = b[i];
590: n = j;
591: if ( !n ) return 0;
592: if ( n == 1 ) return a[0];
593: gcd = a[0];
594: for ( i = 1; i < n; i++ )
595: gcd = gcdz(gcd,a[i]);
596: return gcd;
597: }
1.4 noro 598: #endif
1.1 noro 599:
600: void _copyz(Z n1,Z n2)
601: {
602: int n,i;
603:
604: if ( !n1 || !SL(n1) ) SL(n2) = 0;
605: else {
606: n = SL(n2) = SL(n1);
607: if ( n < 0 ) n = -n;
608: for ( i = 0; i < n; i++ ) BD(n2)[i] = BD(n1)[i];
609: }
610: }
611:
612: void _addz(Z n1,Z n2,Z nr)
613: {
1.3 noro 614: int d1,d2;
1.1 noro 615:
616: if ( !n1 || !SL(n1) ) _copyz(n2,nr);
617: else if ( !n2 || !SL(n2) ) _copyz(n1,nr);
1.3 noro 618: else if ( (d1=SL(n1)) > 0 )
619: if ( (d2=SL(n2)) > 0 )
620: SL(nr) = _addz_main(BD(n1),d1,BD(n2),d2,BD(nr));
1.1 noro 621: else
1.3 noro 622: SL(nr) = _subz_main(BD(n1),d1,BD(n2),-d2,BD(nr));
623: else if ( (d2=SL(n2)) > 0 )
624: SL(nr) = _subz_main(BD(n2),d2,BD(n1),-d1,BD(nr));
1.1 noro 625: else
1.3 noro 626: SL(nr) = -_addz_main(BD(n1),-d1,BD(n2),-d2,BD(nr));
1.1 noro 627: }
628:
629: void _subz(Z n1,Z n2,Z nr)
630: {
1.3 noro 631: int d1,d2;
1.1 noro 632:
633: if ( !n1 || !SL(n1) ) _copyz(n2,nr);
634: else if ( !n2 || !SL(n2) ) {
635: _copyz(n1,nr);
636: SL(nr) = -SL(nr);
1.3 noro 637: } else if ( (d1=SL(n1)) > 0 )
638: if ( (d2=SL(n2)) > 0 )
639: SL(nr) = _subz_main(BD(n1),d1,BD(n2),d2,BD(nr));
1.1 noro 640: else
1.3 noro 641: SL(nr) = _addz_main(BD(n1),d1,BD(n2),-d2,BD(nr));
642: else if ( (d2=SL(n2)) > 0 )
643: SL(nr) = -_addz_main(BD(n1),-d1,BD(n2),d2,BD(nr));
1.1 noro 644: else
1.3 noro 645: SL(nr) = -_subz_main(BD(n1),-d1,BD(n2),-d2,BD(nr));
1.1 noro 646: }
647:
648: void _mulz(Z n1,Z n2,Z nr)
649: {
1.3 noro 650: int d1,d2;
1.1 noro 651: int i,d,sgn;
652: unsigned int mul;
653: unsigned int *m1,*m2;
654:
655: if ( !n1 || !SL(n1) || !n2 || !SL(n2) )
656: SL(nr) = 0;
657: else {
1.3 noro 658: d1 = SL(n1); d2 = SL(n2);
1.1 noro 659: sgn = 1;
1.3 noro 660: if ( d1 < 0 ) { sgn = -sgn; d1 = -d1; }
661: if ( d2 < 0 ) { sgn = -sgn; d2 = -d2; }
662: d = d1+d2;
1.1 noro 663: for ( i = d-1, m1 = BD(nr); i >= 0; i-- ) *m1++ = 0;
1.3 noro 664: for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < d2; i++, m2++ )
665: if ( mul = *m2 ) muln_1(m1,d1,mul,BD(nr)+i);
1.1 noro 666: SL(nr) = sgn*(BD(nr)[d-1]?d:d-1);
667: }
668: }
669:
670: int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)
671: {
672: int d,i;
673: unsigned int tmp,c;
674: unsigned int *t;
675:
676: if ( d2 > d1 ) {
677: t = m1; m1 = m2; m2 = t;
678: d = d1; d1 = d2; d2 = d;
679: }
1.11 ohara 680: #if defined(_M_IX86)
1.1 noro 681: __asm {
682: push esi
683: push edi
684: mov esi,m1
685: mov edi,m2
686: mov ebx,mr
687: mov ecx,d2
688: xor eax,eax
689: Lstart__addz:
690: mov eax,DWORD PTR [esi]
691: mov edx,DWORD PTR [edi]
692: adc eax,edx
693: mov DWORD PTR [ebx],eax
694: lea esi,DWORD PTR [esi+4]
695: lea edi,DWORD PTR [edi+4]
696: lea ebx,DWORD PTR [ebx+4]
697: dec ecx
698: jnz Lstart__addz
699: pop edi
700: pop esi
701: mov eax,0
702: adc eax,eax
703: mov c,eax
704: }
1.10 noro 705: #elif defined(i386)
1.1 noro 706: asm volatile("\
1.10 noro 707: pushl %%ebx;\
1.1 noro 708: movl %1,%%esi;\
709: movl %2,%%edi;\
710: movl %3,%%ebx;\
711: movl %4,%%ecx;\
712: testl %%eax,%%eax;\
713: Lstart__addz:;\
714: movl (%%esi),%%eax;\
715: movl (%%edi),%%edx;\
716: adcl %%edx,%%eax;\
717: movl %%eax,(%%ebx);\
718: leal 4(%%esi),%%esi;\
719: leal 4(%%edi),%%edi;\
720: leal 4(%%ebx),%%ebx;\
721: decl %%ecx;\
722: jnz Lstart__addz;\
723: movl $0,%%eax;\
724: adcl %%eax,%%eax;\
1.10 noro 725: movl %%eax,%0;\
726: popl %%ebx"\
1.1 noro 727: :"=m"(c)\
728: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
1.10 noro 729: :"eax","ecx","edx","esi","edi");
1.1 noro 730: #else
1.6 noro 731: for ( i = 0, c = 0; i < d2; i++, m1++, m2++, mr++ ) {
1.1 noro 732: tmp = *m1 + *m2;
733: if ( tmp < *m1 ) {
734: tmp += c;
735: c = 1;
736: } else {
737: tmp += c;
738: c = tmp < c ? 1 : 0;
739: }
740: *mr = tmp;
741: }
742: #endif
743: for ( i = d2, m1 += d2, mr += d2; (i < d1) && c ; i++ ) {
744: tmp = *m1++ + c;
745: c = tmp < c ? 1 : 0;
746: *mr++ = tmp;
747: }
748: for ( ; i < d1; i++ )
749: *mr++ = *m1++;
750: *mr = c;
751: return (c?d1+1:d1);
752: }
753:
754: int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)
755: {
756: int d,i,sgn;
757: unsigned int t,tmp,br;
758: unsigned int *m;
759:
760: if ( d1 > d2 ) sgn = 1;
761: else if ( d1 < d2 ) sgn = -1;
762: else {
763: for ( i = d1-1; i >= 0 && m1[i] == m2[i]; i-- );
764: if ( i < 0 ) return 0;
765: if ( m1[i] > m2[i] ) sgn = 1;
766: else if ( m1[i] < m2[i] ) sgn = -1;
767: }
768: if ( sgn < 0 ) {
769: m = m1; m1 = m2; m2 = m;
770: d = d1; d1 = d2; d2 = d;
771: }
1.11 ohara 772: #if defined(_M_IX86)
1.1 noro 773: __asm {
774: push esi
775: push edi
776: mov esi,m1
777: mov edi,m2
778: mov ebx,mr
779: mov ecx,d2
780: xor eax,eax
781: Lstart__subz:
782: mov eax,DWORD PTR [esi]
783: mov edx,DWORD PTR [edi]
784: sbb eax,edx
785: mov DWORD PTR [ebx],eax
786: lea esi,DWORD PTR [esi+4]
787: lea edi,DWORD PTR [edi+4]
788: lea ebx,DWORD PTR [ebx+4]
789: dec ecx
790: jnz Lstart__subz
791: pop edi
792: pop esi
793: mov eax,0
794: adc eax,eax
795: mov br,eax
796: }
1.10 noro 797: #elif defined(i386)
1.1 noro 798: asm volatile("\
1.10 noro 799: pushl %%ebx;\
1.1 noro 800: movl %1,%%esi;\
801: movl %2,%%edi;\
802: movl %3,%%ebx;\
803: movl %4,%%ecx;\
804: testl %%eax,%%eax;\
805: Lstart__subz:;\
806: movl (%%esi),%%eax;\
807: movl (%%edi),%%edx;\
808: sbbl %%edx,%%eax;\
809: movl %%eax,(%%ebx);\
810: leal 4(%%esi),%%esi;\
811: leal 4(%%edi),%%edi;\
812: leal 4(%%ebx),%%ebx;\
813: decl %%ecx;\
814: jnz Lstart__subz;\
815: movl $0,%%eax;\
816: adcl %%eax,%%eax;\
1.10 noro 817: movl %%eax,%0;\
818: popl %%ebx"\
1.1 noro 819: :"=m"(br)\
820: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
1.10 noro 821: :"eax","ecx","edx","esi","edi");
1.1 noro 822: #else
1.6 noro 823: for ( i = 0, br = 0; i < d2; i++, mr++ ) {
1.1 noro 824: t = *m1++;
825: tmp = *m2++ + br;
826: if ( br > 0 && !tmp ) {
827: /* tmp = 2^32 => br = 1 */
828: }else {
829: tmp = t-tmp;
830: br = tmp > t ? 1 : 0;
831: *mr = tmp;
832: }
833: }
834: #endif
835: for ( i = d2, m1 += d2, mr += d2; (i < d1) && br; i++ ) {
836: t = *m1++;
837: tmp = t - br;
838: br = tmp > t ? 1 : 0;
839: *mr++ = tmp;
840: }
841: for ( ; i < d1; i++ )
842: *mr++ = *m1++;
843: for ( i = d1-1, mr--; i >= 0 && !*mr--; i-- );
844: return sgn*(i+1);
845: }
846:
847: /* XXX */
848:
849: void printz(Z n)
850: {
1.4 noro 851: int sd,u;
1.1 noro 852:
853: if ( !n )
854: fprintf(asir_out,"0");
1.4 noro 855: else if ( IS_IMM(n) ) {
1.5 noro 856: u = ZTOS(n);
1.4 noro 857: fprintf(asir_out,"%d",u);
858: } else {
1.1 noro 859: if ( (sd = SL(n)) < 0 ) { SL(n) = -SL(n); fprintf(asir_out,"-"); }
860: printn((N)n);
861: if ( sd < 0 ) SL(n) = -SL(n);
1.2 noro 862: }
863: }
864:
865: /*
866: * Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+...
867: *
868: * t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)]
869: * where W(k,l,i) = i! * kCi * lCi
870: */
871:
872: void mkwcz(int k,int l,Z *t)
873: {
874: int i,n,up,low;
875: N nm,d,c;
876:
877: n = MIN(k,l);
878: for ( t[0] = (Z)ONEN, i = 1; i <= n; i++ ) {
879: DM(k-i+1,l-i+1,up,low);
880: if ( up ) {
881: nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = low; BD(nm)[1] = up;
882: } else {
883: nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = low;
884: }
885: kmuln((N)t[i-1],nm,&d); divin(d,i,&c); t[i] = (Z)c;
1.1 noro 886: }
887: }