![[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"
! 2: #include "inline.h"
! 3:
! 4: typedef struct oZ {
! 5: int p;
! 6: unsigned int b[1];
! 7: } *Z;
! 8:
! 9: #define ZALLOC(d) ((Z)MALLOC_ATOMIC(TRUESIZE(oZ,(d)-1,int)))
! 10: #define SL(n) ((n)->p)
! 11:
! 12: Z qtoz(Q n);
! 13: Q ztoq(Z n);
! 14: Z chsgnz(Z n);
! 15: Z dupz(Z n);
! 16: Z addz(Z n1,Z n2);
! 17: Z subz(Z n1,Z n2);
! 18: Z mulz(Z n1,Z n2);
! 19: Z divsz(Z n1,Z n2);
! 20: Z divz(Z n1,Z n2,Z *rem);
! 21: Z gcdz(Z n1,Z n2);
! 22: Z gcdz_cofactor(Z n1,Z n2,Z *c1,Z *c2);
! 23: Z estimate_array_gcdz(Z *a,int n);
! 24: Z array_gcdz(Z *a,int n);
! 25: int remzi(Z n,int m);
! 26: inline void _addz(Z n1,Z n2,Z nr);
! 27: inline void _subz(Z n1,Z n2,Z nr);
! 28: inline void _mulz(Z n1,Z n2,Z nr);
! 29: inline int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr);
! 30: inline int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr);
! 31:
! 32: z_mag(Z n)
! 33: {
! 34: return n_bits((N)n);
! 35: }
! 36:
! 37: Z qtoz(Q n)
! 38: {
! 39: Z r;
! 40:
! 41: if ( !n ) return 0;
! 42: else if ( !INT(n) )
! 43: error("qtoz : invalid input");
! 44: else {
! 45: r = dupz((Z)NM(n));
! 46: if ( SGN(n) < 0 ) SL(r) = -SL(r);
! 47: return r;
! 48: }
! 49: }
! 50:
! 51: Q ztoq(Z n)
! 52: {
! 53: Q r;
! 54: Z nm;
! 55: int sgn;
! 56:
! 57: if ( !n ) return 0;
! 58: else {
! 59: nm = dupz(n);
! 60: if ( SL(nm) < 0 ) {
! 61: sgn = -1;
! 62: SL(nm) = -SL(nm);
! 63: } else
! 64: sgn = 1;
! 65: NTOQ((N)nm,sgn,r);
! 66: return r;
! 67: }
! 68: }
! 69:
! 70: Z dupz(Z n)
! 71: {
! 72: Z nr;
! 73: int sd,i;
! 74:
! 75: if ( !n ) return 0;
! 76: if ( (sd = SL(n)) < 0 ) sd = -sd;
! 77: nr = ZALLOC(sd);
! 78: SL(nr) = SL(n);
! 79: for ( i = 0; i < sd; i++ ) BD(nr)[i] = BD(n)[i];
! 80: return nr;
! 81: }
! 82:
! 83: Z chsgnz(Z n)
! 84: {
! 85: Z r;
! 86:
! 87: if ( !n ) return 0;
! 88: else {
! 89: r = dupz(n);
! 90: SL(r) = -SL(r);
! 91: return r;
! 92: }
! 93: }
! 94:
! 95: Z addz(Z n1,Z n2)
! 96: {
! 97: unsigned int u1,u2,t;
! 98: int sd1,sd2,d,d1,d2;
! 99: Z nr;
! 100:
! 101: if ( !n1 ) return dupz(n2);
! 102: else if ( !n2 ) return dupz(n1);
! 103: else {
! 104: d1 = ((sd1 = SL(n1))< 0)?-sd1:sd1;
! 105: d2 = ((sd2 = SL(n2))< 0)?-sd2:sd2;
! 106: if ( d1 == 1 && d2 == 1 ) {
! 107: u1 = BD(n1)[0]; u2 = BD(n2)[0];
! 108: if ( sd1 > 0 ) {
! 109: if ( sd2 > 0 ) {
! 110: t = u1+u2;
! 111: if ( t < u1 ) {
! 112: nr=ZALLOC(2); SL(nr) = 2; BD(nr)[1] = 1;
! 113: } else {
! 114: nr=ZALLOC(1); SL(nr) = 1;
! 115: }
! 116: BD(nr)[0] = t;
! 117: } else {
! 118: if ( u1 == u2 ) nr = 0;
! 119: else if ( u1 > u2 ) {
! 120: nr=ZALLOC(1); SL(nr) = 1; BD(nr)[0] = u1-u2;
! 121: } else {
! 122: nr=ZALLOC(1); SL(nr) = -1; BD(nr)[0] = u2-u1;
! 123: }
! 124: }
! 125: } else {
! 126: if ( sd2 > 0 ) {
! 127: if ( u2 == u1 ) nr = 0;
! 128: else if ( u2 > u1 ) {
! 129: nr=ZALLOC(1); SL(nr) = 1; BD(nr)[0] = u2-u1;
! 130: } else {
! 131: nr=ZALLOC(1); SL(nr) = -1; BD(nr)[0] = u1-u2;
! 132: }
! 133: } else {
! 134: t = u1+u2;
! 135: if ( t < u1 ) {
! 136: nr=ZALLOC(2); SL(nr) = -2; BD(nr)[1] = 1;
! 137: } else {
! 138: nr=ZALLOC(1); SL(nr) = -1;
! 139: }
! 140: BD(nr)[0] = t;
! 141: }
! 142: }
! 143: } else {
! 144: d = MAX(d1,d2)+1;
! 145: nr = ZALLOC(d);
! 146: _addz(n1,n2,nr);
! 147: if ( !SL(nr) ) nr = 0;
! 148: }
! 149: return nr;
! 150: }
! 151: }
! 152:
! 153: Z subz(Z n1,Z n2)
! 154: {
! 155: int sd1,sd2,d;
! 156: Z nr;
! 157:
! 158: if ( !n1 )
! 159: if ( !n2 ) return 0;
! 160: else {
! 161: nr = dupz(n2);
! 162: SL(nr) = -SL(nr);
! 163: }
! 164: else if ( !n2 ) return dupz(n1);
! 165: else {
! 166: if ( (sd1 = SL(n1)) < 0 ) sd1 = -sd1;
! 167: if ( (sd2 = SL(n2)) < 0 ) sd2 = -sd2;
! 168: d = MAX(sd1,sd2)+1;
! 169: nr = ZALLOC(d);
! 170: _subz(n1,n2,nr);
! 171: if ( !SL(nr) ) nr = 0;
! 172: return nr;
! 173: }
! 174: }
! 175:
! 176: Z mulz(Z n1,Z n2)
! 177: {
! 178: int sd1,sd2,d1,d2;
! 179: unsigned int u1,u2,u,l;
! 180: Z nr;
! 181:
! 182: if ( !n1 || !n2 ) return 0;
! 183: else {
! 184: d1 = ((sd1 = SL(n1))< 0)?-sd1:sd1;
! 185: d2 = ((sd2 = SL(n2))< 0)?-sd2:sd2;
! 186: if ( d1 == 1 && d2 == 1 ) {
! 187: u1 = BD(n1)[0]; u2 = BD(n2)[0];
! 188: DM(u1,u2,u,l);
! 189: if ( u ) {
! 190: nr = ZALLOC(2); SL(nr) = sd1*sd2*2; BD(nr)[1] = u;
! 191: } else {
! 192: nr = ZALLOC(1); SL(nr) = sd1*sd2;
! 193: }
! 194: BD(nr)[0] = l;
! 195: } else {
! 196: nr = ZALLOC(d1+d2);
! 197: _mulz(n1,n2,nr);
! 198: }
! 199: return nr;
! 200: }
! 201: }
! 202:
! 203: /* XXX */
! 204: Z divz(Z n1,Z n2,Z *rem)
! 205: {
! 206: int sgn,sd1,sd2;
! 207: N q,r;
! 208:
! 209: if ( !n2 ) error("divz : division by 0");
! 210: else if ( !n1 ) {
! 211: *rem = 0; return 0;
! 212: } else {
! 213: sd2 = SL(n2);
! 214: if ( sd2 == 1 && BD(n2)[0] == 1 ) {
! 215: *rem = 0; return n1;
! 216: } else if ( sd2 == -1 && BD(n2)[0] == 1 ) {
! 217: *rem = 0; return chsgnz(n1);
! 218: }
! 219:
! 220: sd1 = SL(n1);
! 221: n1 = dupz(n1);
! 222: if ( SL(n1) < 0 ) SL(n1) = -SL(n1);
! 223: divn((N)n1,(N)n2,&q,&r);
! 224: if ( q && ((sd1>0&&sd2<0)||(sd1<0&&sd2>0)) ) SL((Z)q) = -SL((Z)q);
! 225: if ( r && sd1 < 0 ) SL((Z)r) = -SL((Z)r);
! 226: *rem = (Z)r;
! 227: return (Z)q;
! 228: }
! 229: }
! 230:
! 231: Z divsz(Z n1,Z n2)
! 232: {
! 233: int sgn,sd1,sd2;
! 234: N q;
! 235:
! 236: if ( !n2 ) error("divsz : division by 0");
! 237: else if ( !n1 ) return 0;
! 238: else {
! 239: sd2 = SL(n2);
! 240: if ( sd2 == 1 && BD(n2)[0] == 1 ) return n1;
! 241: else if ( sd2 == -1 && BD(n2)[0] == 1 ) return chsgnz(n1);
! 242:
! 243: sd1 = SL(n1);
! 244: n1 = dupz(n1);
! 245: if ( SL(n1) < 0 ) SL(n1) = -SL(n1);
! 246: divsn((N)n1,(N)n2,&q);
! 247: if ( q && ((sd1>0&&sd2<0)||(sd1<0&&sd2>0)) ) SL((Z)q) = -SL((Z)q);
! 248: return (Z)q;
! 249: }
! 250: }
! 251:
! 252: Z gcdz(Z n1,Z n2)
! 253: {
! 254: int sd1,sd2;
! 255: N gcd;
! 256:
! 257: if ( !n1 ) return n2;
! 258: else if ( !n2 ) return n1;
! 259:
! 260: n1 = dupz(n1);
! 261: if ( SL(n1) < 0 ) SL(n1) = -SL(n1);
! 262: n2 = dupz(n2);
! 263: if ( SL(n2) < 0 ) SL(n2) = -SL(n2);
! 264: gcdn((N)n1,(N)n2,&gcd);
! 265: return (Z)gcd;
! 266: }
! 267:
! 268: int remzi(Z n,int m)
! 269: {
! 270: unsigned int *x;
! 271: unsigned int t,r;
! 272: int i;
! 273:
! 274: if ( !n ) return 0;
! 275: i = SL(n);
! 276: if ( i < 0 ) i = -i;
! 277: for ( i--, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) {
! 278: #if defined(sparc)
! 279: r = dsa(m,r,*x);
! 280: #else
! 281: DSAB(m,r,*x,t,r);
! 282: #endif
! 283: }
! 284: return r;
! 285: }
! 286:
! 287: Z gcdz_cofactor(Z n1,Z n2,Z *c1,Z *c2)
! 288: {
! 289: Z gcd;
! 290:
! 291: gcd = gcdz(n1,n2);
! 292: *c1 = divsz(n1,gcd);
! 293: *c2 = divsz(n2,gcd);
! 294: return gcd;
! 295: }
! 296:
! 297: Z estimate_array_gcdz(Z *b,int n)
! 298: {
! 299: int m,i,j,sd;
! 300: Z *a;
! 301: Z s,t;
! 302:
! 303: a = (Z *)ALLOCA(n*sizeof(Z));
! 304: for ( i = j = 0; i < n; i++ ) if ( b[i] ) a[j++] = b[i];
! 305: n = j;
! 306: if ( !n ) return 0;
! 307: if ( n == 1 ) return a[0];
! 308:
! 309: m = n/2;
! 310: for ( i = 0, s = 0; i < m; i++ ) {
! 311: if ( !a[i] ) continue;
! 312: else s = (SL(a[i])<0)?subz(s,a[i]):addz(s,a[i]);
! 313: }
! 314: for ( t = 0; i < n; i++ ) {
! 315: if ( !a[i] ) continue;
! 316: else t = (SL(a[i])<0)?subz(t,a[i]):addz(t,a[i]);
! 317: }
! 318: return gcdz(s,t);
! 319: }
! 320:
! 321: Z array_gcdz(Z *b,int n)
! 322: {
! 323: int m,i,j,sd;
! 324: Z *a;
! 325: Z gcd;
! 326:
! 327: a = (Z *)ALLOCA(n*sizeof(Z));
! 328: for ( i = j = 0; i < n; i++ ) if ( b[i] ) a[j++] = b[i];
! 329: n = j;
! 330: if ( !n ) return 0;
! 331: if ( n == 1 ) return a[0];
! 332: gcd = a[0];
! 333: for ( i = 1; i < n; i++ )
! 334: gcd = gcdz(gcd,a[i]);
! 335: return gcd;
! 336: }
! 337:
! 338: void _copyz(Z n1,Z n2)
! 339: {
! 340: int n,i;
! 341:
! 342: if ( !n1 || !SL(n1) ) SL(n2) = 0;
! 343: else {
! 344: n = SL(n2) = SL(n1);
! 345: if ( n < 0 ) n = -n;
! 346: for ( i = 0; i < n; i++ ) BD(n2)[i] = BD(n1)[i];
! 347: }
! 348: }
! 349:
! 350: void _addz(Z n1,Z n2,Z nr)
! 351: {
! 352: int sd1,sd2;
! 353:
! 354: if ( !n1 || !SL(n1) ) _copyz(n2,nr);
! 355: else if ( !n2 || !SL(n2) ) _copyz(n1,nr);
! 356: else if ( (sd1=SL(n1)) > 0 )
! 357: if ( (sd2=SL(n2)) > 0 )
! 358: SL(nr) = _addz_main(BD(n1),sd1,BD(n2),sd2,BD(nr));
! 359: else
! 360: SL(nr) = _subz_main(BD(n1),sd1,BD(n2),-sd2,BD(nr));
! 361: else if ( (sd2=SL(n2)) > 0 )
! 362: SL(nr) = _subz_main(BD(n2),sd2,BD(n1),-sd1,BD(nr));
! 363: else
! 364: SL(nr) = -_addz_main(BD(n1),-sd1,BD(n2),-sd2,BD(nr));
! 365: }
! 366:
! 367: void _subz(Z n1,Z n2,Z nr)
! 368: {
! 369: int sd1,sd2;
! 370:
! 371: if ( !n1 || !SL(n1) ) _copyz(n2,nr);
! 372: else if ( !n2 || !SL(n2) ) {
! 373: _copyz(n1,nr);
! 374: SL(nr) = -SL(nr);
! 375: } else if ( (sd1=SL(n1)) > 0 )
! 376: if ( (sd2=SL(n2)) > 0 )
! 377: SL(nr) = _subz_main(BD(n1),sd1,BD(n2),sd2,BD(nr));
! 378: else
! 379: SL(nr) = _addz_main(BD(n1),sd1,BD(n2),-sd2,BD(nr));
! 380: else if ( (sd2=SL(n2)) > 0 )
! 381: SL(nr) = -_addz_main(BD(n1),-sd1,BD(n2),sd2,BD(nr));
! 382: else
! 383: SL(nr) = -_subz_main(BD(n1),-sd1,BD(n2),-sd2,BD(nr));
! 384: }
! 385:
! 386: void _mulz(Z n1,Z n2,Z nr)
! 387: {
! 388: int sd1,sd2;
! 389: int i,d,sgn;
! 390: unsigned int mul;
! 391: unsigned int *m1,*m2;
! 392:
! 393: if ( !n1 || !SL(n1) || !n2 || !SL(n2) )
! 394: SL(nr) = 0;
! 395: else {
! 396: sd1 = SL(n1); sd2 = SL(n2);
! 397: sgn = 1;
! 398: if ( sd1 < 0 ) { sgn = -sgn; sd1 = -sd1; }
! 399: if ( sd2 < 0 ) { sgn = -sgn; sd2 = -sd2; }
! 400: d = sd1+sd2;
! 401: for ( i = d-1, m1 = BD(nr); i >= 0; i-- ) *m1++ = 0;
! 402: for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < sd2; i++, m2++ )
! 403: if ( mul = *m2 ) muln_1(m1,sd1,mul,BD(nr)+i);
! 404: SL(nr) = sgn*(BD(nr)[d-1]?d:d-1);
! 405: }
! 406: }
! 407:
! 408: int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)
! 409: {
! 410: int d,i;
! 411: unsigned int tmp,c;
! 412: unsigned int *t;
! 413:
! 414: if ( d2 > d1 ) {
! 415: t = m1; m1 = m2; m2 = t;
! 416: d = d1; d1 = d2; d2 = d;
! 417: }
! 418: #if defined(VISUAL)
! 419: __asm {
! 420: push esi
! 421: push edi
! 422: mov esi,m1
! 423: mov edi,m2
! 424: mov ebx,mr
! 425: mov ecx,d2
! 426: xor eax,eax
! 427: Lstart__addz:
! 428: mov eax,DWORD PTR [esi]
! 429: mov edx,DWORD PTR [edi]
! 430: adc eax,edx
! 431: mov DWORD PTR [ebx],eax
! 432: lea esi,DWORD PTR [esi+4]
! 433: lea edi,DWORD PTR [edi+4]
! 434: lea ebx,DWORD PTR [ebx+4]
! 435: dec ecx
! 436: jnz Lstart__addz
! 437: pop edi
! 438: pop esi
! 439: mov eax,0
! 440: adc eax,eax
! 441: mov c,eax
! 442: }
! 443: #elif defined(i386)
! 444: asm volatile("\
! 445: movl %1,%%esi;\
! 446: movl %2,%%edi;\
! 447: movl %3,%%ebx;\
! 448: movl %4,%%ecx;\
! 449: testl %%eax,%%eax;\
! 450: Lstart__addz:;\
! 451: movl (%%esi),%%eax;\
! 452: movl (%%edi),%%edx;\
! 453: adcl %%edx,%%eax;\
! 454: movl %%eax,(%%ebx);\
! 455: leal 4(%%esi),%%esi;\
! 456: leal 4(%%edi),%%edi;\
! 457: leal 4(%%ebx),%%ebx;\
! 458: decl %%ecx;\
! 459: jnz Lstart__addz;\
! 460: movl $0,%%eax;\
! 461: adcl %%eax,%%eax;\
! 462: movl %%eax,%0"\
! 463: :"=m"(c)\
! 464: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
! 465: :"eax","ebx","ecx","edx","esi","edi");
! 466: #else
! 467: for ( i = 0, c = 0, mr = BD(nr); i < d2; i++, m1++, m2++, mr++ ) {
! 468: tmp = *m1 + *m2;
! 469: if ( tmp < *m1 ) {
! 470: tmp += c;
! 471: c = 1;
! 472: } else {
! 473: tmp += c;
! 474: c = tmp < c ? 1 : 0;
! 475: }
! 476: *mr = tmp;
! 477: }
! 478: #endif
! 479: for ( i = d2, m1 += d2, mr += d2; (i < d1) && c ; i++ ) {
! 480: tmp = *m1++ + c;
! 481: c = tmp < c ? 1 : 0;
! 482: *mr++ = tmp;
! 483: }
! 484: for ( ; i < d1; i++ )
! 485: *mr++ = *m1++;
! 486: *mr = c;
! 487: return (c?d1+1:d1);
! 488: }
! 489:
! 490: int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)
! 491: {
! 492: int d,i,sgn;
! 493: unsigned int t,tmp,br;
! 494: unsigned int *m;
! 495:
! 496: if ( d1 > d2 ) sgn = 1;
! 497: else if ( d1 < d2 ) sgn = -1;
! 498: else {
! 499: for ( i = d1-1; i >= 0 && m1[i] == m2[i]; i-- );
! 500: if ( i < 0 ) return 0;
! 501: if ( m1[i] > m2[i] ) sgn = 1;
! 502: else if ( m1[i] < m2[i] ) sgn = -1;
! 503: }
! 504: if ( sgn < 0 ) {
! 505: m = m1; m1 = m2; m2 = m;
! 506: d = d1; d1 = d2; d2 = d;
! 507: }
! 508: #if defined(VISUAL)
! 509: __asm {
! 510: push esi
! 511: push edi
! 512: mov esi,m1
! 513: mov edi,m2
! 514: mov ebx,mr
! 515: mov ecx,d2
! 516: xor eax,eax
! 517: Lstart__subz:
! 518: mov eax,DWORD PTR [esi]
! 519: mov edx,DWORD PTR [edi]
! 520: sbb eax,edx
! 521: mov DWORD PTR [ebx],eax
! 522: lea esi,DWORD PTR [esi+4]
! 523: lea edi,DWORD PTR [edi+4]
! 524: lea ebx,DWORD PTR [ebx+4]
! 525: dec ecx
! 526: jnz Lstart__subz
! 527: pop edi
! 528: pop esi
! 529: mov eax,0
! 530: adc eax,eax
! 531: mov br,eax
! 532: }
! 533: #elif defined(i386)
! 534: asm volatile("\
! 535: movl %1,%%esi;\
! 536: movl %2,%%edi;\
! 537: movl %3,%%ebx;\
! 538: movl %4,%%ecx;\
! 539: testl %%eax,%%eax;\
! 540: Lstart__subz:;\
! 541: movl (%%esi),%%eax;\
! 542: movl (%%edi),%%edx;\
! 543: sbbl %%edx,%%eax;\
! 544: movl %%eax,(%%ebx);\
! 545: leal 4(%%esi),%%esi;\
! 546: leal 4(%%edi),%%edi;\
! 547: leal 4(%%ebx),%%ebx;\
! 548: decl %%ecx;\
! 549: jnz Lstart__subz;\
! 550: movl $0,%%eax;\
! 551: adcl %%eax,%%eax;\
! 552: movl %%eax,%0"\
! 553: :"=m"(br)\
! 554: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
! 555: :"eax","ebx","ecx","edx","esi","edi");
! 556: #else
! 557: for ( i = 0, br = 0, mr = BD(nr); i < d2; i++, mr++ ) {
! 558: t = *m1++;
! 559: tmp = *m2++ + br;
! 560: if ( br > 0 && !tmp ) {
! 561: /* tmp = 2^32 => br = 1 */
! 562: }else {
! 563: tmp = t-tmp;
! 564: br = tmp > t ? 1 : 0;
! 565: *mr = tmp;
! 566: }
! 567: }
! 568: #endif
! 569: for ( i = d2, m1 += d2, mr += d2; (i < d1) && br; i++ ) {
! 570: t = *m1++;
! 571: tmp = t - br;
! 572: br = tmp > t ? 1 : 0;
! 573: *mr++ = tmp;
! 574: }
! 575: for ( ; i < d1; i++ )
! 576: *mr++ = *m1++;
! 577: for ( i = d1-1, mr--; i >= 0 && !*mr--; i-- );
! 578: return sgn*(i+1);
! 579: }
! 580:
! 581: /* XXX */
! 582:
! 583: void printz(Z n)
! 584: {
! 585: int sd;
! 586:
! 587: if ( !n )
! 588: fprintf(asir_out,"0");
! 589: else {
! 590: if ( (sd = SL(n)) < 0 ) { SL(n) = -SL(n); fprintf(asir_out,"-"); }
! 591: printn((N)n);
! 592: if ( sd < 0 ) SL(n) = -SL(n);
! 593: }
! 594: }