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