Annotation of OpenXM_contrib2/asir2000/engine/Z.c, Revision 1.16
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.15 fujimoto 7: #elif defined(VISUAL)
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: {
1.16 ! noro 33: if ( IS_IMM(a) && ZTOS(a) == 1 ) return 1;
! 34: else return 0;
1.5 noro 35: }
36:
37: int cmpz(Z a,Z b)
38: {
1.16 ! noro 39: int *ma,*mb;
! 40: int sa,sb,da,db,ca,cb,i;
1.5 noro 41:
1.16 ! noro 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: }
1.5 noro 72: }
1.1 noro 73:
1.5 noro 74: Z stoz(int c)
1.3 noro 75: {
1.16 ! noro 76: Z z;
1.3 noro 77:
1.16 ! noro 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;
1.3 noro 85: }
86:
1.5 noro 87: Z utoz(unsigned int c)
88: {
1.16 ! noro 89: Z z;
1.5 noro 90:
1.16 ! noro 91: z = ZALLOC(1);
! 92: SL(z) = 1; BD(z)[0] = c;
! 93: return z;
1.5 noro 94: }
95:
96: Z simpz(Z n)
97: {
1.16 ! noro 98: Z n1;
1.5 noro 99:
1.16 ! noro 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;
1.5 noro 106: }
107:
1.3 noro 108: int sgnz(Z n)
1.2 noro 109: {
1.16 ! noro 110: if ( !n ) return 0;
! 111: else if ( IS_IMM(n) ) return ZTOS(n)>0?1:-1;
! 112: else if ( SL(n) < 0 ) return -1;
! 113: else return 1;
1.2 noro 114: }
115:
1.1 noro 116: z_mag(Z n)
117: {
1.16 ! noro 118: int c,i;
1.3 noro 119:
1.16 ! noro 120: if ( !n ) return 0;
! 121: else if ( IS_IMM(n) ) {
! 122: c = ZTOS(n);
! 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.16 ! noro 132: Z r,t;
! 133: int c;
1.1 noro 134:
1.16 ! noro 135: if ( !n ) return 0;
! 136: else if ( !INT(n) )
! 137: error("qtoz : invalid input");
! 138: else {
! 139: t = (Z)NM(n);
! 140: if ( IS_SZ(t) ) {
! 141: c = SGN(n) < 0 ? -BD(t)[0] : BD(t)[0];
! 142: STOZ(c,r);
! 143: } else {
! 144: r = dupz((Z)t);
! 145: if ( SGN(n) < 0 ) SL(r) = -SL(r);
! 146: }
! 147: return r;
! 148: }
1.1 noro 149: }
150:
151: Q ztoq(Z n)
152: {
1.16 ! noro 153: Q r;
! 154: Z nm;
! 155: int sgn,c;
! 156:
! 157: if ( !n ) return 0;
! 158: else if ( IS_IMM(n) ) {
! 159: c = ZTOS(n);
! 160: STOQ(c,r);
! 161: return r;
! 162: } else {
! 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: }
1.1 noro 172: }
173:
174: Z dupz(Z n)
175: {
1.16 ! noro 176: Z r;
! 177: int sd,i;
1.1 noro 178:
1.16 ! noro 179: if ( !n ) return 0;
! 180: else if ( IS_IMM(n) ) return n;
! 181: else {
! 182: if ( (sd = SL(n)) < 0 ) sd = -sd;
! 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;
! 187: }
1.1 noro 188: }
189:
190: Z chsgnz(Z n)
191: {
1.16 ! noro 192: Z r;
! 193: int c;
1.1 noro 194:
1.16 ! noro 195: if ( !n ) return 0;
! 196: else if ( IS_IMM(n) ) {
! 197: c = -ZTOS(n);
! 198: STOZ(c,r);
! 199: return r;
! 200: } else {
! 201: r = dupz(n);
! 202: SL(r) = -SL(r);
! 203: return r;
! 204: }
1.1 noro 205: }
206:
1.5 noro 207: Z absz(Z n)
208: {
1.16 ! noro 209: Z r;
! 210: int c;
1.5 noro 211:
1.16 ! noro 212: if ( !n ) return 0;
! 213: else if ( sgnz(n) > 0 )
! 214: return n;
! 215: else
! 216: return chsgnz(n);
1.5 noro 217: }
1.3 noro 218:
1.1 noro 219: Z addz(Z n1,Z n2)
220: {
1.16 ! noro 221: int d1,d2,d,c;
! 222: Z r,r1;
! 223: struct oZ t;
! 224:
! 225: if ( !n1 ) return dupz(n2);
! 226: else if ( !n2 ) return dupz(n1);
! 227: else if ( IS_IMM(n1) ) {
! 228: if ( IS_IMM(n2) ) {
! 229: c = ZTOS(n1)+ZTOS(n2);
! 230: STOZ(c,r);
! 231: } else {
! 232: c = ZTOS(n1);
! 233: if ( c < 0 ) {
! 234: t.p = -1; t.b[0] = -c;
! 235: } else {
! 236: t.p = 1; t.b[0] = c;
! 237: }
! 238: if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
! 239: r = ZALLOC(d2+1);
! 240: _addz(&t,n2,r);
! 241: if ( !SL(r) ) r = 0;
! 242: }
! 243: } else if ( IS_IMM(n2) ) {
! 244: c = ZTOS(n2);
! 245: if ( c < 0 ) {
! 246: t.p = -1; t.b[0] = -c;
! 247: } else {
! 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);
! 253: if ( !SL(r) ) r = 0;
! 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;
! 261: }
! 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.16 ! noro 270: int d1,d2,d,c;
! 271: Z r,r1;
! 272: struct oZ t;
! 273:
! 274: if ( !n1 )
! 275: return chsgnz(n2);
! 276: else if ( !n2 ) return n1;
! 277: else if ( IS_IMM(n1) ) {
! 278: if ( IS_IMM(n2) ) {
! 279: c = ZTOS(n1)-ZTOS(n2);
! 280: STOZ(c,r);
! 281: } else {
! 282: c = ZTOS(n1);
! 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);
! 291: if ( !SL(r) ) r = 0;
! 292: }
! 293: } else if ( IS_IMM(n2) ) {
! 294: c = ZTOS(n2);
! 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);
! 303: if ( !SL(r) ) r = 0;
! 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;
! 311: }
! 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.16 ! noro 320: int d1,d2,sgn,i;
! 321: int c1,c2;
! 322: unsigned int u1,u2,u,l;
! 323: Z r;
! 324:
! 325: if ( !n1 || !n2 ) return 0;
! 326:
! 327: if ( IS_IMM(n1) ) {
! 328: c1 = ZTOS(n1);
! 329: if ( IS_IMM(n2) ) {
! 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: }
! 350: }
! 351: } else {
! 352: sgn = 1;
! 353: if ( c1 < 0 ) { c1 = -c1; sgn = -sgn; }
! 354: u1 = (unsigned int)c1;
! 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;
! 360: }
! 361: } else if ( IS_IMM(n2) ) {
! 362: c2 = ZTOS(n2);
! 363: if ( c2 == 1 )
! 364: return n1;
! 365: else if ( c2 == -1 )
! 366: return chsgnz(n1);
! 367:
! 368: sgn = 1;
! 369: if ( c2 < 0 ) { sgn = -sgn; c2 = -c2; }
! 370: u2 = (unsigned int)c2;
! 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;
! 378: if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
! 379: if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
! 380: r = ZALLOC(d1+d2);
! 381: _mulz(n1,n2,r);
! 382: }
! 383: if ( sgn < 0 ) r = chsgnz(r);
! 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.16 ! noro 391: int sgn,d1,d2;
! 392: Z q;
1.1 noro 393:
1.16 ! noro 394: if ( !n2 ) error("divsz : division by 0");
! 395: if ( !n1 ) return 0;
1.3 noro 396:
1.16 ! noro 397: if ( IS_IMM(n1) ) {
! 398: if ( !IS_IMM(n2) )
! 399: error("divsz : cannot happen");
! 400: c = ZTOS(n1)/ZTOS(n2);
! 401: STOZ(c,q);
! 402: return q;
! 403: }
! 404: if ( IS_IMM(n2) ) {
! 405: sgn = 1;
! 406: u2 = ZTOS(n2); if ( u2 < 0 ) { sgn = -sgn; u2 = -u2; }
! 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: {
1.16 ! noro 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);
1.5 noro 445: }
446:
447: Z divsz(Z n1,Z n2)
448: {
1.16 ! noro 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);
1.5 noro 463: }
464:
465: int gcdimm(int c1,int c2)
466: {
1.16 ! noro 467: int r;
1.5 noro 468:
1.16 ! noro 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: }
1.5 noro 476: }
1.1 noro 477:
478: Z gcdz(Z n1,Z n2)
479: {
1.16 ! noro 480: int c1,c2,g;
! 481: Z gcd,r;
! 482: N gn;
! 483:
! 484: if ( !n1 ) return n2;
! 485: else if ( !n2 ) return n1;
! 486:
! 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: {
1.16 ! noro 522: unsigned int *x;
! 523: unsigned int t,r;
! 524: int i,c;
! 525:
! 526: if ( !n ) return 0;
! 527: if ( IS_IMM(n) ) {
! 528: c = ZTOS(n)%m;
! 529: if ( c < 0 ) c += m;
! 530: return c;
! 531: }
! 532:
! 533: i = SL(n);
! 534: if ( i < 0 ) i = -i;
! 535: for ( i--, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) {
1.1 noro 536: #if defined(sparc)
1.16 ! noro 537: r = dsa(m,r,*x);
1.1 noro 538: #else
1.16 ! noro 539: DSAB(m,r,*x,t,r);
1.1 noro 540: #endif
1.16 ! noro 541: }
! 542: if ( r && SL(n) < 0 )
! 543: r = m-r;
! 544: return r;
1.1 noro 545: }
546:
547: Z gcdz_cofactor(Z n1,Z n2,Z *c1,Z *c2)
548: {
1.16 ! noro 549: Z gcd;
1.1 noro 550:
1.16 ! noro 551: gcd = gcdz(n1,n2);
! 552: *c1 = divsz(n1,gcd);
! 553: *c2 = divsz(n2,gcd);
! 554: return gcd;
1.1 noro 555: }
556:
1.5 noro 557: #if 0
1.1 noro 558: Z estimate_array_gcdz(Z *b,int n)
559: {
1.16 ! noro 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);
1.1 noro 580: }
581:
582: Z array_gcdz(Z *b,int n)
583: {
1.16 ! noro 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;
1.1 noro 597: }
1.4 noro 598: #endif
1.1 noro 599:
600: void _copyz(Z n1,Z n2)
601: {
1.16 ! noro 602: int n,i;
1.1 noro 603:
1.16 ! noro 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: }
1.1 noro 610: }
611:
612: void _addz(Z n1,Z n2,Z nr)
613: {
1.16 ! noro 614: int d1,d2;
1.1 noro 615:
1.16 ! noro 616: if ( !n1 || !SL(n1) ) _copyz(n2,nr);
! 617: else if ( !n2 || !SL(n2) ) _copyz(n1,nr);
! 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));
! 621: else
! 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));
! 625: else
! 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.16 ! noro 631: int d1,d2;
1.1 noro 632:
1.16 ! noro 633: if ( !n1 || !SL(n1) ) _copyz(n2,nr);
! 634: else if ( !n2 || !SL(n2) ) {
! 635: _copyz(n1,nr);
! 636: SL(nr) = -SL(nr);
! 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));
! 640: else
! 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));
! 644: else
! 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.16 ! noro 650: int d1,d2;
! 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 {
! 658: d1 = SL(n1); d2 = SL(n2);
! 659: sgn = 1;
! 660: if ( d1 < 0 ) { sgn = -sgn; d1 = -d1; }
! 661: if ( d2 < 0 ) { sgn = -sgn; d2 = -d2; }
! 662: d = d1+d2;
! 663: for ( i = d-1, m1 = BD(nr); i >= 0; i-- ) *m1++ = 0;
! 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);
! 666: SL(nr) = sgn*(BD(nr)[d-1]?d:d-1);
! 667: }
1.1 noro 668: }
669:
670: int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)
671: {
1.16 ! noro 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.15 fujimoto 680: #if defined(_M_IX86) && !defined(__MINGW32__)
1.16 ! 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.15 fujimoto 705: #elif defined(i386) && !defined(__MINGW32__)
1.16 ! noro 706: asm volatile("\
! 707: pushl %%ebx;\
! 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;\
! 725: movl %%eax,%0;\
! 726: popl %%ebx"\
! 727: :"=m"(c)\
! 728: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
! 729: :"eax","ecx","edx","esi","edi");
1.1 noro 730: #else
1.16 ! noro 731: for ( i = 0, c = 0; i < d2; i++, m1++, m2++, mr++ ) {
! 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: }
1.1 noro 742: #endif
1.16 ! noro 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);
1.1 noro 752: }
753:
754: int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)
755: {
1.16 ! noro 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.15 fujimoto 772: #if defined(_M_IX86) && !defined(__MINGW32__)
1.16 ! 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.15 fujimoto 797: #elif defined(i386) && !defined(__MINGW32__)
1.16 ! noro 798: asm volatile("\
! 799: pushl %%ebx;\
! 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;\
! 817: movl %%eax,%0;\
! 818: popl %%ebx"\
! 819: :"=m"(br)\
! 820: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
! 821: :"eax","ecx","edx","esi","edi");
1.1 noro 822: #else
1.16 ! noro 823: for ( i = 0, br = 0; i < d2; i++, mr++ ) {
! 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: }
1.1 noro 834: #endif
1.16 ! noro 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);
1.1 noro 845: }
846:
847: /* XXX */
848:
849: void printz(Z n)
850: {
1.16 ! noro 851: int sd,u;
1.1 noro 852:
1.16 ! noro 853: if ( !n )
! 854: fprintf(asir_out,"0");
! 855: else if ( IS_IMM(n) ) {
! 856: u = ZTOS(n);
! 857: fprintf(asir_out,"%d",u);
! 858: } else {
! 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);
! 862: }
1.2 noro 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: {
1.16 ! noro 874: int i,n,up,low;
! 875: N nm,d,c;
1.2 noro 876:
1.16 ! noro 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;
! 886: }
1.1 noro 887: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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