Annotation of OpenXM_contrib2/asir2000/engine/N.c, Revision 1.1
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/N.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #include "ca.h"
! 3: #include "base.h"
! 4:
! 5: void addn();
! 6: int subn();
! 7: void _bshiftn();
! 8: void dupn();
! 9:
! 10: #if defined(VISUAL) || defined(i386)
! 11: void addn(n1,n2,nr)
! 12: N n1,n2,*nr;
! 13: {
! 14: unsigned int *m1,*m2,*mr;
! 15: unsigned int c;
! 16: N r;
! 17: int i,d1,d2;
! 18: unsigned int tmp;
! 19:
! 20: if ( !n1 )
! 21: COPY(n2,*nr);
! 22: else if ( !n2 )
! 23: COPY(n1,*nr);
! 24: else {
! 25: if ( PL(n1) > PL(n2) ) {
! 26: d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
! 27: } else {
! 28: d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
! 29: }
! 30: *nr = r = NALLOC(d1 + 1); INITRC(r); mr = BD(r);
! 31:
! 32: #if defined(VISUAL)
! 33: __asm {
! 34: push esi
! 35: push edi
! 36: mov esi,m1
! 37: mov edi,m2
! 38: mov ebx,mr
! 39: mov ecx,d2
! 40: xor eax,eax
! 41: Lstart_addn:
! 42: mov eax,DWORD PTR [esi]
! 43: mov edx,DWORD PTR [edi]
! 44: adc eax,edx
! 45: mov DWORD PTR [ebx],eax
! 46: lea esi,DWORD PTR [esi+4]
! 47: lea edi,DWORD PTR [edi+4]
! 48: lea ebx,DWORD PTR [ebx+4]
! 49: dec ecx
! 50: jnz Lstart_addn
! 51: pop edi
! 52: pop esi
! 53: mov eax,0
! 54: adc eax,eax
! 55: mov c,eax
! 56: }
! 57: #else
! 58: asm volatile("\
! 59: movl %1,%%esi;\
! 60: movl %2,%%edi;\
! 61: movl %3,%%ebx;\
! 62: movl %4,%%ecx;\
! 63: testl %%eax,%%eax;\
! 64: Lstart_addn:;\
! 65: movl (%%esi),%%eax;\
! 66: movl (%%edi),%%edx;\
! 67: adcl %%edx,%%eax;\
! 68: movl %%eax,(%%ebx);\
! 69: leal 4(%%esi),%%esi;\
! 70: leal 4(%%edi),%%edi;\
! 71: leal 4(%%ebx),%%ebx;\
! 72: decl %%ecx;\
! 73: jnz Lstart_addn;\
! 74: movl $0,%%eax;\
! 75: adcl %%eax,%%eax;\
! 76: movl %%eax,%0"\
! 77: :"=m"(c)\
! 78: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
! 79: :"eax","ebx","ecx","edx","esi","edi");
! 80: #endif
! 81: for ( i = d2, m1 += d2, mr += d2; (i < d1) && c ; i++ ) {
! 82: tmp = *m1++ + c;
! 83: c = tmp < c ? 1 : 0;
! 84: *mr++ = tmp;
! 85: }
! 86: for ( ; i < d1; i++ )
! 87: *mr++ = *m1++;
! 88: *mr = c;
! 89: PL(r) = (c?d1+1:d1);
! 90: }
! 91: }
! 92:
! 93: int subn(n1,n2,nr)
! 94: N n1,n2,*nr;
! 95: {
! 96: N r;
! 97: unsigned int *m1,*m2,*mr,br;
! 98: unsigned int tmp,t;
! 99: int d1,d2,sgn,i;
! 100:
! 101: if ( !n1 ) {
! 102: if ( n2 ) {
! 103: COPY(n2,*nr);
! 104: return -1;
! 105: } else {
! 106: *nr = 0;
! 107: return 0;
! 108: }
! 109: } else if ( !n2 ) {
! 110: COPY(n1,*nr);
! 111: return 1;
! 112: } else {
! 113: d1 = PL(n1); d2 = PL(n2);
! 114: m1 = BD(n1); m2 = BD(n2);
! 115: if ( (d1 = PL(n1)) > (d2 = PL(n2)) )
! 116: sgn = 1;
! 117: else if ( d1 < d2 ) {
! 118: d1 = PL(n2); d2 = PL(n1);
! 119: m1 = BD(n2); m2 = BD(n1);
! 120: sgn = -1;
! 121: } else {
! 122: for ( i = d1-1; i >= 0 && m1[i] == m2[i]; i-- );
! 123: if ( i < 0 ) {
! 124: *nr = 0;
! 125: return 0;
! 126: }
! 127: d1 = d2 = i+1;
! 128: if ( m1[i] > m2[i] )
! 129: sgn = 1;
! 130: else {
! 131: m1 = BD(n2); m2 = BD(n1);
! 132: sgn = -1;
! 133: }
! 134: }
! 135: *nr = r = NALLOC(d1); INITRC(r); mr = BD(r);
! 136:
! 137: #if defined(VISUAL)
! 138: __asm {
! 139: push esi
! 140: push edi
! 141: mov esi,m1
! 142: mov edi,m2
! 143: mov ebx,mr
! 144: mov ecx,d2
! 145: xor eax,eax
! 146: Lstart_subn:
! 147: mov eax,DWORD PTR [esi]
! 148: mov edx,DWORD PTR [edi]
! 149: sbb eax,edx
! 150: mov DWORD PTR [ebx],eax
! 151: lea esi,DWORD PTR [esi+4]
! 152: lea edi,DWORD PTR [edi+4]
! 153: lea ebx,DWORD PTR [ebx+4]
! 154: dec ecx
! 155: jnz Lstart_subn
! 156: pop edi
! 157: pop esi
! 158: mov eax,0
! 159: adc eax,eax
! 160: mov br,eax
! 161: }
! 162: #else
! 163: asm volatile("\
! 164: movl %1,%%esi;\
! 165: movl %2,%%edi;\
! 166: movl %3,%%ebx;\
! 167: movl %4,%%ecx;\
! 168: testl %%eax,%%eax;\
! 169: Lstart_subn:;\
! 170: movl (%%esi),%%eax;\
! 171: movl (%%edi),%%edx;\
! 172: sbbl %%edx,%%eax;\
! 173: movl %%eax,(%%ebx);\
! 174: leal 4(%%esi),%%esi;\
! 175: leal 4(%%edi),%%edi;\
! 176: leal 4(%%ebx),%%ebx;\
! 177: decl %%ecx;\
! 178: jnz Lstart_subn;\
! 179: movl $0,%%eax;\
! 180: adcl %%eax,%%eax;\
! 181: movl %%eax,%0"\
! 182: :"=m"(br)\
! 183: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
! 184: :"eax","ebx","ecx","edx","esi","edi");
! 185: #endif
! 186: for ( i = d2, m1 += d2, mr += d2; (i < d1) && br; i++ ) {
! 187: t = *m1++;
! 188: tmp = t - br;
! 189: br = tmp > t ? 1 : 0;
! 190: *mr++ = tmp;
! 191: }
! 192: for ( ; i < d1; i++ )
! 193: *mr++ = *m1++;
! 194: for ( i = d1-1, mr--; i >= 0 && !*mr--; i-- );
! 195: PL(r) = i + 1;
! 196: return sgn;
! 197: }
! 198: }
! 199:
! 200: void _addn(n1,n2,nr)
! 201: N n1,n2,nr;
! 202: {
! 203: unsigned int *m1,*m2,*mr;
! 204: unsigned int c;
! 205: int i,d1,d2;
! 206: unsigned int tmp;
! 207:
! 208: if ( !n1 || !PL(n1) )
! 209: dupn(n2,nr);
! 210: else if ( !n2 || !PL(n2) )
! 211: dupn(n1,nr);
! 212: else {
! 213: if ( PL(n1) > PL(n2) ) {
! 214: d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
! 215: } else {
! 216: d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
! 217: }
! 218: mr = BD(nr);
! 219:
! 220: #if defined(VISUAL)
! 221: __asm {
! 222: push esi
! 223: push edi
! 224: mov esi,m1
! 225: mov edi,m2
! 226: mov ebx,mr
! 227: mov ecx,d2
! 228: xor eax,eax
! 229: Lstart__addn:
! 230: mov eax,DWORD PTR [esi]
! 231: mov edx,DWORD PTR [edi]
! 232: adc eax,edx
! 233: mov DWORD PTR [ebx],eax
! 234: lea esi,DWORD PTR [esi+4]
! 235: lea edi,DWORD PTR [edi+4]
! 236: lea ebx,DWORD PTR [ebx+4]
! 237: dec ecx
! 238: jnz Lstart__addn
! 239: pop edi
! 240: pop esi
! 241: mov eax,0
! 242: adc eax,eax
! 243: mov c,eax
! 244: }
! 245: #else
! 246: asm volatile("\
! 247: movl %1,%%esi;\
! 248: movl %2,%%edi;\
! 249: movl %3,%%ebx;\
! 250: movl %4,%%ecx;\
! 251: testl %%eax,%%eax;\
! 252: Lstart__addn:;\
! 253: movl (%%esi),%%eax;\
! 254: movl (%%edi),%%edx;\
! 255: adcl %%edx,%%eax;\
! 256: movl %%eax,(%%ebx);\
! 257: leal 4(%%esi),%%esi;\
! 258: leal 4(%%edi),%%edi;\
! 259: leal 4(%%ebx),%%ebx;\
! 260: decl %%ecx;\
! 261: jnz Lstart__addn;\
! 262: movl $0,%%eax;\
! 263: adcl %%eax,%%eax;\
! 264: movl %%eax,%0"\
! 265: :"=m"(c)\
! 266: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
! 267: :"eax","ebx","ecx","edx","esi","edi");
! 268: #endif
! 269: for ( i = d2, m1 += d2, mr += d2; (i < d1) && c ; i++ ) {
! 270: tmp = *m1++ + c;
! 271: c = tmp < c ? 1 : 0;
! 272: *mr++ = tmp;
! 273: }
! 274: for ( ; i < d1; i++ )
! 275: *mr++ = *m1++;
! 276: *mr = c;
! 277: PL(nr) = (c?d1+1:d1);
! 278: }
! 279: }
! 280:
! 281: int _subn(n1,n2,nr)
! 282: N n1,n2,nr;
! 283: {
! 284: unsigned int *m1,*m2,*mr,br;
! 285: unsigned int tmp,t;
! 286: int d1,d2,sgn,i;
! 287:
! 288: if ( !n1 || !PL(n1) ) {
! 289: if ( n2 && PL(n2) ) {
! 290: dupn(n2,nr);
! 291: return -1;
! 292: } else {
! 293: PL(nr) = 0;
! 294: return 0;
! 295: }
! 296: } else if ( !n2 || !PL(n2) ) {
! 297: dupn(n1,nr);
! 298: return 1;
! 299: } else {
! 300: d1 = PL(n1); d2 = PL(n2);
! 301: m1 = BD(n1); m2 = BD(n2);
! 302: if ( (d1 = PL(n1)) > (d2 = PL(n2)) )
! 303: sgn = 1;
! 304: else if ( d1 < d2 ) {
! 305: d1 = PL(n2); d2 = PL(n1);
! 306: m1 = BD(n2); m2 = BD(n1);
! 307: sgn = -1;
! 308: } else {
! 309: for ( i = d1-1; i >= 0 && m1[i] == m2[i]; i-- );
! 310: if ( i < 0 ) {
! 311: PL(nr) = 0;
! 312: return 0;
! 313: }
! 314: d1 = d2 = i+1;
! 315: if ( m1[i] > m2[i] )
! 316: sgn = 1;
! 317: else {
! 318: m1 = BD(n2); m2 = BD(n1);
! 319: sgn = -1;
! 320: }
! 321: }
! 322: mr = BD(nr);
! 323:
! 324: #if defined(VISUAL)
! 325: __asm {
! 326: push esi
! 327: push edi
! 328: mov esi,m1
! 329: mov edi,m2
! 330: mov ebx,mr
! 331: mov ecx,d2
! 332: xor eax,eax
! 333: Lstart__subn:
! 334: mov eax,DWORD PTR [esi]
! 335: mov edx,DWORD PTR [edi]
! 336: sbb eax,edx
! 337: mov DWORD PTR [ebx],eax
! 338: lea esi,DWORD PTR [esi+4]
! 339: lea edi,DWORD PTR [edi+4]
! 340: lea ebx,DWORD PTR [ebx+4]
! 341: dec ecx
! 342: jnz Lstart__subn
! 343: pop edi
! 344: pop esi
! 345: mov eax,0
! 346: adc eax,eax
! 347: mov br,eax
! 348: }
! 349: #else
! 350: asm volatile("\
! 351: movl %1,%%esi;\
! 352: movl %2,%%edi;\
! 353: movl %3,%%ebx;\
! 354: movl %4,%%ecx;\
! 355: testl %%eax,%%eax;\
! 356: Lstart__subn:;\
! 357: movl (%%esi),%%eax;\
! 358: movl (%%edi),%%edx;\
! 359: sbbl %%edx,%%eax;\
! 360: movl %%eax,(%%ebx);\
! 361: leal 4(%%esi),%%esi;\
! 362: leal 4(%%edi),%%edi;\
! 363: leal 4(%%ebx),%%ebx;\
! 364: decl %%ecx;\
! 365: jnz Lstart__subn;\
! 366: movl $0,%%eax;\
! 367: adcl %%eax,%%eax;\
! 368: movl %%eax,%0"\
! 369: :"=m"(br)\
! 370: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
! 371: :"eax","ebx","ecx","edx","esi","edi");
! 372: #endif
! 373: for ( i = d2, m1 += d2, mr += d2; (i < d1) && br; i++ ) {
! 374: t = *m1++;
! 375: tmp = t - br;
! 376: br = tmp > t ? 1 : 0;
! 377: *mr++ = tmp;
! 378: }
! 379: for ( ; i < d1; i++ )
! 380: *mr++ = *m1++;
! 381: for ( i = d1-1, mr--; i >= 0 && !*mr--; i-- );
! 382: PL(nr) = i + 1;
! 383: return sgn;
! 384: }
! 385: }
! 386: #else
! 387:
! 388: void addn(n1,n2,nr)
! 389: N n1,n2,*nr;
! 390: {
! 391: unsigned int *m1,*m2,*mr,i,c;
! 392: N r;
! 393: int d1,d2;
! 394: unsigned int tmp;
! 395:
! 396: if ( !n1 )
! 397: COPY(n2,*nr);
! 398: else if ( !n2 )
! 399: COPY(n1,*nr);
! 400: else {
! 401: if ( PL(n1) > PL(n2) ) {
! 402: d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
! 403: } else {
! 404: d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
! 405: }
! 406: *nr = r = NALLOC(d1 + 1); INITRC(r);
! 407: for ( i = 0, c = 0, mr = BD(r); i < d2; i++, m1++, m2++, mr++ ) {
! 408: tmp = *m1 + *m2;
! 409: if ( tmp < *m1 ) {
! 410: tmp += c;
! 411: c = 1;
! 412: } else {
! 413: tmp += c;
! 414: c = tmp < c ? 1 : 0;
! 415: }
! 416: *mr = tmp;
! 417: }
! 418: for ( ; (i < d1) && c ; i++, m1++, mr++ ) {
! 419: tmp = *m1 + c;
! 420: c = tmp < c ? 1 : 0;
! 421: *mr = tmp;
! 422: }
! 423: for ( ; i < d1; i++ )
! 424: *mr++ = *m1++;
! 425: *mr = c;
! 426: PL(r) = (c?d1+1:d1);
! 427: }
! 428: }
! 429:
! 430: int subn(n1,n2,nr)
! 431: N n1,n2,*nr;
! 432: {
! 433: N r;
! 434: unsigned int *m1,*m2,*mr,i,br;
! 435: L tmp;
! 436: int d1,d2,nz,sgn;
! 437:
! 438: if ( !n1 ) {
! 439: if ( n2 ) {
! 440: COPY(n2,*nr);
! 441: return -1;
! 442: } else {
! 443: *nr = 0;
! 444: return 0;
! 445: }
! 446: } else if ( !n2 ) {
! 447: COPY(n1,*nr);
! 448: return 1;
! 449: } else {
! 450: switch ( cmpn(n1,n2) ) {
! 451: case 1:
! 452: d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
! 453: sgn = 1; break;
! 454: case -1:
! 455: d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
! 456: sgn = -1; break;
! 457: case 0:
! 458: default:
! 459: *nr = 0; return ( 0 ); break;
! 460: }
! 461: *nr = r = NALLOC(d1); INITRC(r);
! 462: for ( i = 0, br = 0, nz = -1, mr = BD(r);
! 463: i < d2; i++ ) {
! 464: if ( (tmp = (L)*m1++ - (L)*m2++ - (L)br) && ( tmp > -LBASE ) )
! 465: nz = i;
! 466: if ( tmp < 0 ) {
! 467: br = 1; *mr++ = (unsigned int)(tmp + LBASE);
! 468: } else {
! 469: br = 0; *mr++ = (unsigned int)tmp;
! 470: }
! 471: }
! 472: for ( ; (i < d1) && br; i++ ) {
! 473: if ( (tmp = (L)*m1++ - (L)br) && ( tmp > -LBASE ) )
! 474: nz = i;
! 475: if ( tmp < 0 ) {
! 476: br = 1; *mr++ = (unsigned int)(tmp + LBASE);
! 477: } else {
! 478: br = 0; *mr++ = (unsigned int)tmp;
! 479: }
! 480: }
! 481: for ( ; i < d1; i++ )
! 482: if ( *mr++ = *m1++ )
! 483: nz = i;
! 484: PL(r) = nz + 1;
! 485: return sgn;
! 486: }
! 487: }
! 488:
! 489: void _addn(n1,n2,nr)
! 490: N n1,n2,nr;
! 491: {
! 492: unsigned int *m1,*m2,*mr,i,c;
! 493: int d1,d2;
! 494: unsigned int tmp;
! 495:
! 496: if ( !n1 || !PL(n1) )
! 497: dupn(n2,nr);
! 498: else if ( !n2 || !PL(n2) )
! 499: dupn(n1,nr);
! 500: else {
! 501: if ( PL(n1) > PL(n2) ) {
! 502: d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
! 503: } else {
! 504: d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
! 505: }
! 506: for ( i = 0, c = 0, mr = BD(nr); i < d2; i++, m1++, m2++, mr++ ) {
! 507: tmp = *m1 + *m2;
! 508: if ( tmp < *m1 ) {
! 509: tmp += c;
! 510: c = 1;
! 511: } else {
! 512: tmp += c;
! 513: c = tmp < c ? 1 : 0;
! 514: }
! 515: *mr = tmp;
! 516: }
! 517: for ( ; (i < d1) && c ; i++, m1++, mr++ ) {
! 518: tmp = *m1 + c;
! 519: c = tmp < c ? 1 : 0;
! 520: *mr = tmp;
! 521: }
! 522: for ( ; i < d1; i++ )
! 523: *mr++ = *m1++;
! 524: *mr = c;
! 525: PL(nr) = (c?d1+1:d1);
! 526: }
! 527: }
! 528:
! 529: int _subn(n1,n2,nr)
! 530: N n1,n2,nr;
! 531: {
! 532: N r;
! 533: unsigned int *m1,*m2,*mr,i,br;
! 534: L tmp;
! 535: int d1,d2,nz,sgn;
! 536:
! 537: if ( !n1 || !PL(n1) ) {
! 538: if ( n2 && PL(n2) ) {
! 539: dupn(n2,nr);
! 540: return -1;
! 541: } else {
! 542: PL(nr) = 0;
! 543: return 0;
! 544: }
! 545: } else if ( !n2 || !PL(n2) ) {
! 546: dupn(n1,nr);
! 547: return 1;
! 548: } else {
! 549: switch ( cmpn(n1,n2) ) {
! 550: case 1:
! 551: d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
! 552: sgn = 1; break;
! 553: case -1:
! 554: d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
! 555: sgn = -1; break;
! 556: case 0:
! 557: default:
! 558: PL(nr) = 0; return ( 0 ); break;
! 559: }
! 560: for ( i = 0, br = 0, nz = -1, mr = BD(nr);
! 561: i < d2; i++ ) {
! 562: if ( (tmp = (L)*m1++ - (L)*m2++ - (L)br) && ( tmp > -LBASE ) )
! 563: nz = i;
! 564: if ( tmp < 0 ) {
! 565: br = 1; *mr++ = (unsigned int)(tmp + LBASE);
! 566: } else {
! 567: br = 0; *mr++ = (unsigned int)tmp;
! 568: }
! 569: }
! 570: for ( ; (i < d1) && br; i++ ) {
! 571: if ( (tmp = (L)*m1++ - (L)br) && ( tmp > -LBASE ) )
! 572: nz = i;
! 573: if ( tmp < 0 ) {
! 574: br = 1; *mr++ = (unsigned int)(tmp + LBASE);
! 575: } else {
! 576: br = 0; *mr++ = (unsigned int)tmp;
! 577: }
! 578: }
! 579: for ( ; i < d1; i++ )
! 580: if ( *mr++ = *m1++ )
! 581: nz = i;
! 582: PL(nr) = nz + 1;
! 583: return sgn;
! 584: }
! 585: }
! 586: #endif
! 587:
! 588: /* a2 += a1; n2 >= n1 */
! 589:
! 590: void addarray_to(a1,n1,a2,n2)
! 591: unsigned int *a1;
! 592: int n1;
! 593: unsigned int *a2;
! 594: int n2;
! 595: {
! 596: int i;
! 597: unsigned int c,tmp;
! 598:
! 599: for ( i = 0, c = 0; i < n2; i++, a1++, a2++ ) {
! 600: tmp = *a1 + *a2;
! 601: if ( tmp < *a1 ) {
! 602: tmp += c;
! 603: c = 1;
! 604: } else {
! 605: tmp += c;
! 606: c = tmp < c ? 1 : 0;
! 607: }
! 608: *a2 = tmp;
! 609: }
! 610: for ( ; (i < n2) && c ; i++, a2++ ) {
! 611: tmp = *a2 + c;
! 612: c = tmp < c ? 1 : 0;
! 613: *a2 = tmp;
! 614: }
! 615: if ( i == n2 && c )
! 616: *a2 = c;
! 617: }
! 618:
! 619: void pwrn(n,e,nr)
! 620: N n,*nr;
! 621: int e;
! 622: {
! 623: N nw,nw1;
! 624:
! 625: if ( e == 1 ) {
! 626: COPY(n,*nr);
! 627: return;
! 628: }
! 629: pwrn(n,e / 2,&nw);
! 630: if ( e % 2 == 0 )
! 631: kmuln(nw,nw,nr);
! 632: else {
! 633: kmuln(nw,nw,&nw1); kmuln(nw1,n,nr); FREEN(nw1);
! 634: }
! 635: FREEN(nw);
! 636: }
! 637:
! 638: extern int igcd_algorithm;
! 639:
! 640: void gcdEuclidn();
! 641:
! 642: void gcdn(n1,n2,nr)
! 643: N n1,n2,*nr;
! 644: {
! 645: N m1,m2,g;
! 646:
! 647: if ( !igcd_algorithm )
! 648: gcdEuclidn(n1,n2,nr);
! 649: else {
! 650: if ( !n1 )
! 651: *nr = n2;
! 652: else if ( !n2 )
! 653: *nr = n1;
! 654: else {
! 655: n32ton27(n1,&m1);
! 656: n32ton27(n2,&m2);
! 657: gcdBinary_27n(m1,m2,&g);
! 658: n27ton32(g,nr);
! 659: }
! 660: }
! 661: }
! 662:
! 663: void gcdEuclidn(n1,n2,nr)
! 664: N n1,n2,*nr;
! 665: {
! 666: N m1,m2,q,r;
! 667: unsigned int i1,i2,ir;
! 668:
! 669: if ( !n1 )
! 670: COPY(n2,*nr);
! 671: else if ( !n2 )
! 672: COPY(n1,*nr);
! 673: else {
! 674: if ( PL(n1) > PL(n2) ) {
! 675: COPY(n1,m1); COPY(n2,m2);
! 676: } else {
! 677: COPY(n1,m2); COPY(n2,m1);
! 678: }
! 679: while ( PL(m1) > 1 ) {
! 680: divn(m1,m2,&q,&r); FREEN(m1); FREEN(q);
! 681: if ( !r ) {
! 682: *nr = m2;
! 683: return;
! 684: } else {
! 685: m1 = m2; m2 = r;
! 686: }
! 687: }
! 688: for ( i1 = BD(m1)[0], i2 = BD(m2)[0]; ir = i1 % i2; i1 = i2, i2 = ir );
! 689: if ( i2 == 1 )
! 690: COPY(ONEN,*nr);
! 691: else {
! 692: *nr = r = NALLOC(1); INITRC(r); PL(r) = 1; BD(r)[0] = i2;
! 693: }
! 694: }
! 695: }
! 696:
! 697: int cmpn(n1,n2)
! 698: N n1,n2;
! 699: {
! 700: int i;
! 701: unsigned int *m1,*m2;
! 702:
! 703: if ( !n1 )
! 704: if ( !n2 )
! 705: return 0;
! 706: else
! 707: return -1;
! 708: else if ( !n2 )
! 709: return 1;
! 710: else if ( PL(n1) > PL(n2) )
! 711: return 1;
! 712: else if ( PL(n1) < PL(n2) )
! 713: return -1;
! 714: else {
! 715: for ( i = PL(n1)-1, m1 = BD(n1)+i, m2 = BD(n2)+i;
! 716: i >= 0; i--, m1--, m2-- )
! 717: if ( *m1 > *m2 )
! 718: return 1;
! 719: else if ( *m1 < *m2 )
! 720: return -1;
! 721: return 0;
! 722: }
! 723: }
! 724:
! 725: void bshiftn(n,b,r)
! 726: N n;
! 727: int b;
! 728: N *r;
! 729: {
! 730: int w,l,nl,i,j;
! 731: N z;
! 732: unsigned int msw;
! 733: unsigned int *p,*pz;
! 734:
! 735: if ( b == 0 ) {
! 736: COPY(n,*r); return;
! 737: }
! 738: if ( b > 0 ) { /* >> */
! 739: w = b / BSH; l = PL(n)-w;
! 740: if ( l <= 0 ) {
! 741: *r = 0; return;
! 742: }
! 743: b %= BSH; p = BD(n)+w;
! 744: if ( !b ) {
! 745: *r = z = NALLOC(l); INITRC(z); PL(z) = l;
! 746: bcopy(p,BD(z),l*sizeof(unsigned int));
! 747: return;
! 748: }
! 749: msw = p[l-1];
! 750: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
! 751: if ( b >= i ) {
! 752: l--;
! 753: if ( !l ) {
! 754: *r = 0; return;
! 755: }
! 756: *r = z = NALLOC(l); INITRC(z); PL(z) = l; pz = BD(z);
! 757: for ( j = 0; j < l; j++, p++ )
! 758: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
! 759: } else {
! 760: *r = z = NALLOC(l); INITRC(z); PL(z) = l; pz = BD(z);
! 761: for ( j = 1; j < l; j++, p++ )
! 762: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
! 763: *pz = *p>>b;
! 764: }
! 765: } else { /* << */
! 766: b = -b;
! 767: w = b / BSH; b %= BSH; l = PL(n); p = BD(n);
! 768: if ( !b ) {
! 769: nl = l+w; *r = z = NALLOC(nl); INITRC(z); PL(z) = nl;
! 770: bzero((char *)BD(z),w*sizeof(unsigned int));
! 771: bcopy(p,BD(z)+w,l*sizeof(unsigned int));
! 772: return;
! 773: }
! 774: msw = p[l-1];
! 775: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
! 776: if ( b + i > BSH ) {
! 777: nl = l+w+1;
! 778: *r = z = NALLOC(nl); INITRC(z); PL(z) = nl; pz = BD(z)+w;
! 779: bzero((char *)BD(z),w*sizeof(unsigned int));
! 780: *pz++ = *p++<<b;
! 781: for ( j = 1; j < l; j++, p++ )
! 782: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
! 783: *pz = *(p-1)>>(BSH-b);
! 784: } else {
! 785: nl = l+w;
! 786: *r = z = NALLOC(nl); INITRC(z); PL(z) = nl; pz = BD(z)+w;
! 787: bzero((char *)BD(z),w*sizeof(unsigned int));
! 788: *pz++ = *p++<<b;
! 789: for ( j = 1; j < l; j++, p++ )
! 790: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
! 791: }
! 792: }
! 793: }
! 794:
! 795: #if 0
! 796: void _bshiftn(n,b,z)
! 797: N n;
! 798: int b;
! 799: N z;
! 800: {
! 801: int w,l,nl,i,j;
! 802: unsigned int msw;
! 803: unsigned int *p,*pz;
! 804:
! 805: if ( b == 0 ) {
! 806: copyn(n,PL(n),BD(z)); PL(z) = PL(n); return;
! 807: }
! 808: if ( b > 0 ) { /* >> */
! 809: w = b / BSH; l = PL(n)-w;
! 810: if ( l <= 0 ) {
! 811: PL(z) = 0; return;
! 812: }
! 813: b %= BSH; p = BD(n)+w;
! 814: if ( !b ) {
! 815: PL(z) = l;
! 816: bcopy(p,BD(z),l*sizeof(unsigned int));
! 817: return;
! 818: }
! 819: msw = p[l-1];
! 820: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
! 821: if ( b >= i ) {
! 822: l--;
! 823: if ( !l ) {
! 824: PL(z) = 0; return;
! 825: }
! 826: PL(z) = l; pz = BD(z);
! 827: for ( j = 0; j < l; j++, p++ )
! 828: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
! 829: } else {
! 830: PL(z) = l; pz = BD(z);
! 831: for ( j = 1; j < l; j++, p++ )
! 832: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
! 833: *pz = *p>>b;
! 834: }
! 835: } else { /* << */
! 836: b = -b;
! 837: w = b / BSH; b %= BSH; l = PL(n); p = BD(n);
! 838: if ( !b ) {
! 839: nl = l+w; PL(z) = nl;
! 840: bzero((char *)BD(z),w*sizeof(unsigned int));
! 841: bcopy(p,BD(z)+w,l*sizeof(unsigned int));
! 842: return;
! 843: }
! 844: msw = p[l-1];
! 845: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
! 846: if ( b + i > BSH ) {
! 847: nl = l+w+1;
! 848: PL(z) = nl; pz = BD(z)+w;
! 849: bzero((char *)BD(z),w*sizeof(unsigned int));
! 850: *pz++ = *p++<<b;
! 851: for ( j = 1; j < l; j++, p++ )
! 852: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
! 853: *pz = *(p-1)>>(BSH-b);
! 854: } else {
! 855: nl = l+w;
! 856: PL(z) = nl; pz = BD(z)+w;
! 857: bzero((char *)BD(z),w*sizeof(unsigned int));
! 858: *pz++ = *p++<<b;
! 859: for ( j = 1; j < l; j++, p++ )
! 860: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
! 861: }
! 862: }
! 863: }
! 864: #endif
! 865:
! 866: void shiftn(n,w,r)
! 867: N n;
! 868: int w;
! 869: N *r;
! 870: {
! 871: int l,nl;
! 872: N z;
! 873:
! 874: if ( w == 0 )
! 875: COPY(n,*r);
! 876: else if ( w > 0 ) { /* >> */
! 877: l = PL(n)-w;
! 878: if ( l <= 0 )
! 879: *r = 0;
! 880: else {
! 881: *r = z = NALLOC(l); INITRC(z); PL(z) = l;
! 882: bcopy(BD(n)+w,BD(z),l*sizeof(unsigned int));
! 883: }
! 884: } else { /* << */
! 885: w = -w;
! 886: l = PL(n); nl = l+w;
! 887: *r = z = NALLOC(nl); INITRC(z); PL(z) = nl;
! 888: bzero((char *)BD(z),w*sizeof(unsigned int));
! 889: bcopy(BD(n),BD(z)+w,l*sizeof(unsigned int));
! 890: }
! 891: }
! 892:
! 893: void randomn(bits,r)
! 894: int bits;
! 895: N *r;
! 896: {
! 897: int l,i;
! 898: unsigned int *tb;
! 899: N t;
! 900:
! 901: l = (bits+31)>>5; /* word length */
! 902: *r = t = NALLOC(l);
! 903: tb = BD(t);
! 904: for ( i = 0; i < l; i++ )
! 905: tb[i] = mt_genrand();
! 906: if ( bits&31 )
! 907: tb[l-1] &= (1<<(bits&31))-1;
! 908: for ( i = l-1; i >= 0 && !tb[i]; i-- );
! 909: if ( i < 0 )
! 910: *r = 0;
! 911: else
! 912: PL(t) = i+1;
! 913: }
! 914:
! 915: void freen(n)
! 916: N n;
! 917: {
! 918: if ( n && (n != ONEN) )
! 919: free(n);
! 920: }
! 921:
! 922: int n_bits(n)
! 923: N n;
! 924: {
! 925: unsigned int l,i,t;
! 926:
! 927: if ( !n )
! 928: return 0;
! 929: l = PL(n); t = BD(n)[l-1];
! 930: for ( i = 0; t; t>>=1, i++);
! 931: return i + (l-1)*BSH;
! 932: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>