1 /**************************************************************** 2 3 The author of this software is David M. Gay. 4 5 Copyright (C) 1998, 1999 by Lucent Technologies 6 All Rights Reserved 7 8 Permission to use, copy, modify, and distribute this software and 9 its documentation for any purpose and without fee is hereby 10 granted, provided that the above copyright notice appear in all 11 copies and that both that the copyright notice and this 12 permission notice and warranty disclaimer appear in supporting 13 documentation, and that the name of Lucent or any of its entities 14 not be used in advertising or publicity pertaining to 15 distribution of the software without specific, written prior 16 permission. 17 18 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, 19 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. 20 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY 21 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 22 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER 23 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 24 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF 25 THIS SOFTWARE. 26 27 ****************************************************************/ 28 29 /* Please send bug reports to 30 David M. Gay 31 Bell Laboratories, Room 2C-463 32 600 Mountain Avenue 33 Murray Hill, NJ 07974-0636 34 U.S.A. 35 dmg@bell-labs.com 36 */ 37 38 #include "gdtoaimp.h" 39 40 static Bigint *freelist[Kmax+1]; 41 #ifndef Omit_Private_Memory 42 #ifndef PRIVATE_MEM 43 #define PRIVATE_MEM 2304 44 #endif 45 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) 46 static double private_mem[PRIVATE_mem], *pmem_next = private_mem; 47 #endif 48 49 Bigint * 50 Balloc 51 #ifdef KR_headers 52 (k) int k; 53 #else 54 (int k) 55 #endif 56 { 57 int x; 58 Bigint *rv; 59 #ifndef Omit_Private_Memory 60 unsigned int len; 61 #endif 62 63 ACQUIRE_DTOA_LOCK(0); 64 if ( (rv = freelist[k]) !=0) { 65 freelist[k] = rv->next; 66 } 67 else { 68 x = 1 << k; 69 #ifdef Omit_Private_Memory 70 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); 71 #else 72 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) 73 /sizeof(double); 74 if (pmem_next - private_mem + len <= PRIVATE_mem) { 75 rv = (Bigint*)pmem_next; 76 pmem_next += len; 77 } 78 else 79 rv = (Bigint*)MALLOC(len*sizeof(double)); 80 #endif 81 rv->k = k; 82 rv->maxwds = x; 83 } 84 FREE_DTOA_LOCK(0); 85 rv->sign = rv->wds = 0; 86 return rv; 87 } 88 89 void 90 Bfree 91 #ifdef KR_headers 92 (v) Bigint *v; 93 #else 94 (Bigint *v) 95 #endif 96 { 97 if (v) { 98 ACQUIRE_DTOA_LOCK(0); 99 v->next = freelist[v->k]; 100 freelist[v->k] = v; 101 FREE_DTOA_LOCK(0); 102 } 103 } 104 105 int 106 lo0bits 107 #ifdef KR_headers 108 (y) ULong *y; 109 #else 110 (ULong *y) 111 #endif 112 { 113 register int k; 114 register ULong x = *y; 115 116 if (x & 7) { 117 if (x & 1) 118 return 0; 119 if (x & 2) { 120 *y = x >> 1; 121 return 1; 122 } 123 *y = x >> 2; 124 return 2; 125 } 126 k = 0; 127 if (!(x & 0xffff)) { 128 k = 16; 129 x >>= 16; 130 } 131 if (!(x & 0xff)) { 132 k += 8; 133 x >>= 8; 134 } 135 if (!(x & 0xf)) { 136 k += 4; 137 x >>= 4; 138 } 139 if (!(x & 0x3)) { 140 k += 2; 141 x >>= 2; 142 } 143 if (!(x & 1)) { 144 k++; 145 x >>= 1; 146 if (!x & 1) 147 return 32; 148 } 149 *y = x; 150 return k; 151 } 152 153 Bigint * 154 multadd 155 #ifdef KR_headers 156 (b, m, a) Bigint *b; int m, a; 157 #else 158 (Bigint *b, int m, int a) /* multiply by m and add a */ 159 #endif 160 { 161 int i, wds; 162 #ifdef ULLong 163 ULong *x; 164 ULLong carry, y; 165 #else 166 ULong carry, *x, y; 167 #ifdef Pack_32 168 ULong xi, z; 169 #endif 170 #endif 171 Bigint *b1; 172 173 wds = b->wds; 174 x = b->x; 175 i = 0; 176 carry = a; 177 do { 178 #ifdef ULLong 179 y = *x * (ULLong)m + carry; 180 carry = y >> 32; 181 *x++ = y & 0xffffffffUL; 182 #else 183 #ifdef Pack_32 184 xi = *x; 185 y = (xi & 0xffff) * m + carry; 186 z = (xi >> 16) * m + (y >> 16); 187 carry = z >> 16; 188 *x++ = (z << 16) + (y & 0xffff); 189 #else 190 y = *x * m + carry; 191 carry = y >> 16; 192 *x++ = y & 0xffff; 193 #endif 194 #endif 195 } 196 while(++i < wds); 197 if (carry) { 198 if (wds >= b->maxwds) { 199 b1 = Balloc(b->k+1); 200 Bcopy(b1, b); 201 Bfree(b); 202 b = b1; 203 } 204 b->x[wds++] = carry; 205 b->wds = wds; 206 } 207 return b; 208 } 209 210 int 211 hi0bits 212 #ifdef KR_headers 213 (x) register ULong x; 214 #else 215 (register ULong x) 216 #endif 217 { 218 register int k = 0; 219 220 if (!(x & 0xffff0000)) { 221 k = 16; 222 x <<= 16; 223 } 224 if (!(x & 0xff000000)) { 225 k += 8; 226 x <<= 8; 227 } 228 if (!(x & 0xf0000000)) { 229 k += 4; 230 x <<= 4; 231 } 232 if (!(x & 0xc0000000)) { 233 k += 2; 234 x <<= 2; 235 } 236 if (!(x & 0x80000000)) { 237 k++; 238 if (!(x & 0x40000000)) 239 return 32; 240 } 241 return k; 242 } 243 244 Bigint * 245 i2b 246 #ifdef KR_headers 247 (i) int i; 248 #else 249 (int i) 250 #endif 251 { 252 Bigint *b; 253 254 b = Balloc(1); 255 b->x[0] = i; 256 b->wds = 1; 257 return b; 258 } 259 260 Bigint * 261 mult 262 #ifdef KR_headers 263 (a, b) Bigint *a, *b; 264 #else 265 (Bigint *a, Bigint *b) 266 #endif 267 { 268 Bigint *c; 269 int k, wa, wb, wc; 270 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 271 ULong y; 272 #ifdef ULLong 273 ULLong carry, z; 274 #else 275 ULong carry, z; 276 #ifdef Pack_32 277 ULong z2; 278 #endif 279 #endif 280 281 if (a->wds < b->wds) { 282 c = a; 283 a = b; 284 b = c; 285 } 286 k = a->k; 287 wa = a->wds; 288 wb = b->wds; 289 wc = wa + wb; 290 if (wc > a->maxwds) 291 k++; 292 c = Balloc(k); 293 for(x = c->x, xa = x + wc; x < xa; x++) 294 *x = 0; 295 xa = a->x; 296 xae = xa + wa; 297 xb = b->x; 298 xbe = xb + wb; 299 xc0 = c->x; 300 #ifdef ULLong 301 for(; xb < xbe; xc0++) { 302 if ( (y = *xb++) !=0) { 303 x = xa; 304 xc = xc0; 305 carry = 0; 306 do { 307 z = *x++ * (ULLong)y + *xc + carry; 308 carry = z >> 32; 309 *xc++ = z & 0xffffffffUL; 310 } 311 while(x < xae); 312 *xc = carry; 313 } 314 } 315 #else 316 #ifdef Pack_32 317 for(; xb < xbe; xb++, xc0++) { 318 if ( (y = *xb & 0xffff) !=0) { 319 x = xa; 320 xc = xc0; 321 carry = 0; 322 do { 323 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 324 carry = z >> 16; 325 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 326 carry = z2 >> 16; 327 Storeinc(xc, z2, z); 328 } 329 while(x < xae); 330 *xc = carry; 331 } 332 if ( (y = *xb >> 16) !=0) { 333 x = xa; 334 xc = xc0; 335 carry = 0; 336 z2 = *xc; 337 do { 338 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 339 carry = z >> 16; 340 Storeinc(xc, z, z2); 341 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 342 carry = z2 >> 16; 343 } 344 while(x < xae); 345 *xc = z2; 346 } 347 } 348 #else 349 for(; xb < xbe; xc0++) { 350 if ( (y = *xb++) !=0) { 351 x = xa; 352 xc = xc0; 353 carry = 0; 354 do { 355 z = *x++ * y + *xc + carry; 356 carry = z >> 16; 357 *xc++ = z & 0xffff; 358 } 359 while(x < xae); 360 *xc = carry; 361 } 362 } 363 #endif 364 #endif 365 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 366 c->wds = wc; 367 return c; 368 } 369 370 static Bigint *p5s; 371 372 Bigint * 373 pow5mult 374 #ifdef KR_headers 375 (b, k) Bigint *b; int k; 376 #else 377 (Bigint *b, int k) 378 #endif 379 { 380 Bigint *b1, *p5, *p51; 381 int i; 382 static int p05[3] = { 5, 25, 125 }; 383 384 if ( (i = k & 3) !=0) 385 b = multadd(b, p05[i-1], 0); 386 387 if (!(k >>= 2)) 388 return b; 389 if ((p5 = p5s) == 0) { 390 /* first time */ 391 #ifdef MULTIPLE_THREADS 392 ACQUIRE_DTOA_LOCK(1); 393 if (!(p5 = p5s)) { 394 p5 = p5s = i2b(625); 395 p5->next = 0; 396 } 397 FREE_DTOA_LOCK(1); 398 #else 399 p5 = p5s = i2b(625); 400 p5->next = 0; 401 #endif 402 } 403 for(;;) { 404 if (k & 1) { 405 b1 = mult(b, p5); 406 Bfree(b); 407 b = b1; 408 } 409 if (!(k >>= 1)) 410 break; 411 if ((p51 = p5->next) == 0) { 412 #ifdef MULTIPLE_THREADS 413 ACQUIRE_DTOA_LOCK(1); 414 if (!(p51 = p5->next)) { 415 p51 = p5->next = mult(p5,p5); 416 p51->next = 0; 417 } 418 FREE_DTOA_LOCK(1); 419 #else 420 p51 = p5->next = mult(p5,p5); 421 p51->next = 0; 422 #endif 423 } 424 p5 = p51; 425 } 426 return b; 427 } 428 429 Bigint * 430 lshift 431 #ifdef KR_headers 432 (b, k) Bigint *b; int k; 433 #else 434 (Bigint *b, int k) 435 #endif 436 { 437 int i, k1, n, n1; 438 Bigint *b1; 439 ULong *x, *x1, *xe, z; 440 441 n = k >> kshift; 442 k1 = b->k; 443 n1 = n + b->wds + 1; 444 for(i = b->maxwds; n1 > i; i <<= 1) 445 k1++; 446 b1 = Balloc(k1); 447 x1 = b1->x; 448 for(i = 0; i < n; i++) 449 *x1++ = 0; 450 x = b->x; 451 xe = x + b->wds; 452 if (k &= kmask) { 453 #ifdef Pack_32 454 k1 = 32 - k; 455 z = 0; 456 do { 457 *x1++ = *x << k | z; 458 z = *x++ >> k1; 459 } 460 while(x < xe); 461 if ((*x1 = z) !=0) 462 ++n1; 463 #else 464 k1 = 16 - k; 465 z = 0; 466 do { 467 *x1++ = *x << k & 0xffff | z; 468 z = *x++ >> k1; 469 } 470 while(x < xe); 471 if (*x1 = z) 472 ++n1; 473 #endif 474 } 475 else do 476 *x1++ = *x++; 477 while(x < xe); 478 b1->wds = n1 - 1; 479 Bfree(b); 480 return b1; 481 } 482 483 int 484 cmp 485 #ifdef KR_headers 486 (a, b) Bigint *a, *b; 487 #else 488 (Bigint *a, Bigint *b) 489 #endif 490 { 491 ULong *xa, *xa0, *xb, *xb0; 492 int i, j; 493 494 i = a->wds; 495 j = b->wds; 496 #ifdef DEBUG 497 if (i > 1 && !a->x[i-1]) 498 Bug("cmp called with a->x[a->wds-1] == 0"); 499 if (j > 1 && !b->x[j-1]) 500 Bug("cmp called with b->x[b->wds-1] == 0"); 501 #endif 502 if (i -= j) 503 return i; 504 xa0 = a->x; 505 xa = xa0 + j; 506 xb0 = b->x; 507 xb = xb0 + j; 508 for(;;) { 509 if (*--xa != *--xb) 510 return *xa < *xb ? -1 : 1; 511 if (xa <= xa0) 512 break; 513 } 514 return 0; 515 } 516 517 Bigint * 518 diff 519 #ifdef KR_headers 520 (a, b) Bigint *a, *b; 521 #else 522 (Bigint *a, Bigint *b) 523 #endif 524 { 525 Bigint *c; 526 int i, wa, wb; 527 ULong *xa, *xae, *xb, *xbe, *xc; 528 #ifdef ULLong 529 ULLong borrow, y; 530 #else 531 ULong borrow, y; 532 #ifdef Pack_32 533 ULong z; 534 #endif 535 #endif 536 537 i = cmp(a,b); 538 if (!i) { 539 c = Balloc(0); 540 c->wds = 1; 541 c->x[0] = 0; 542 return c; 543 } 544 if (i < 0) { 545 c = a; 546 a = b; 547 b = c; 548 i = 1; 549 } 550 else 551 i = 0; 552 c = Balloc(a->k); 553 c->sign = i; 554 wa = a->wds; 555 xa = a->x; 556 xae = xa + wa; 557 wb = b->wds; 558 xb = b->x; 559 xbe = xb + wb; 560 xc = c->x; 561 borrow = 0; 562 #ifdef ULLong 563 do { 564 y = (ULLong)*xa++ - *xb++ - borrow; 565 borrow = y >> 32 & 1UL; 566 *xc++ = y & 0xffffffffUL; 567 } 568 while(xb < xbe); 569 while(xa < xae) { 570 y = *xa++ - borrow; 571 borrow = y >> 32 & 1UL; 572 *xc++ = y & 0xffffffffUL; 573 } 574 #else 575 #ifdef Pack_32 576 do { 577 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 578 borrow = (y & 0x10000) >> 16; 579 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 580 borrow = (z & 0x10000) >> 16; 581 Storeinc(xc, z, y); 582 } 583 while(xb < xbe); 584 while(xa < xae) { 585 y = (*xa & 0xffff) - borrow; 586 borrow = (y & 0x10000) >> 16; 587 z = (*xa++ >> 16) - borrow; 588 borrow = (z & 0x10000) >> 16; 589 Storeinc(xc, z, y); 590 } 591 #else 592 do { 593 y = *xa++ - *xb++ - borrow; 594 borrow = (y & 0x10000) >> 16; 595 *xc++ = y & 0xffff; 596 } 597 while(xb < xbe); 598 while(xa < xae) { 599 y = *xa++ - borrow; 600 borrow = (y & 0x10000) >> 16; 601 *xc++ = y & 0xffff; 602 } 603 #endif 604 #endif 605 while(!*--xc) 606 wa--; 607 c->wds = wa; 608 return c; 609 } 610 611 double 612 b2d 613 #ifdef KR_headers 614 (a, e) Bigint *a; int *e; 615 #else 616 (Bigint *a, int *e) 617 #endif 618 { 619 ULong *xa, *xa0, w, y, z; 620 int k; 621 double d; 622 #ifdef VAX 623 ULong d0, d1; 624 #else 625 #define d0 word0(d) 626 #define d1 word1(d) 627 #endif 628 629 xa0 = a->x; 630 xa = xa0 + a->wds; 631 y = *--xa; 632 #ifdef DEBUG 633 if (!y) Bug("zero y in b2d"); 634 #endif 635 k = hi0bits(y); 636 *e = 32 - k; 637 #ifdef Pack_32 638 if (k < Ebits) { 639 d0 = Exp_1 | y >> Ebits - k; 640 w = xa > xa0 ? *--xa : 0; 641 d1 = y << (32-Ebits) + k | w >> Ebits - k; 642 goto ret_d; 643 } 644 z = xa > xa0 ? *--xa : 0; 645 if (k -= Ebits) { 646 d0 = Exp_1 | y << k | z >> 32 - k; 647 y = xa > xa0 ? *--xa : 0; 648 d1 = z << k | y >> 32 - k; 649 } 650 else { 651 d0 = Exp_1 | y; 652 d1 = z; 653 } 654 #else 655 if (k < Ebits + 16) { 656 z = xa > xa0 ? *--xa : 0; 657 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 658 w = xa > xa0 ? *--xa : 0; 659 y = xa > xa0 ? *--xa : 0; 660 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 661 goto ret_d; 662 } 663 z = xa > xa0 ? *--xa : 0; 664 w = xa > xa0 ? *--xa : 0; 665 k -= Ebits + 16; 666 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 667 y = xa > xa0 ? *--xa : 0; 668 d1 = w << k + 16 | y << k; 669 #endif 670 ret_d: 671 #ifdef VAX 672 word0(d) = d0 >> 16 | d0 << 16; 673 word1(d) = d1 >> 16 | d1 << 16; 674 #endif 675 return dval(d); 676 } 677 #undef d0 678 #undef d1 679 680 Bigint * 681 d2b 682 #ifdef KR_headers 683 (d, e, bits) double d; int *e, *bits; 684 #else 685 (double d, int *e, int *bits) 686 #endif 687 { 688 Bigint *b; 689 int de, i, k; 690 ULong *x, y, z; 691 #ifdef VAX 692 ULong d0, d1; 693 d0 = word0(d) >> 16 | word0(d) << 16; 694 d1 = word1(d) >> 16 | word1(d) << 16; 695 #else 696 #define d0 word0(d) 697 #define d1 word1(d) 698 #endif 699 700 #ifdef Pack_32 701 b = Balloc(1); 702 #else 703 b = Balloc(2); 704 #endif 705 x = b->x; 706 707 z = d0 & Frac_mask; 708 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 709 #ifdef Sudden_Underflow 710 de = (int)(d0 >> Exp_shift); 711 #ifndef IBM 712 z |= Exp_msk11; 713 #endif 714 #else 715 if ( (de = (int)(d0 >> Exp_shift)) !=0) 716 z |= Exp_msk1; 717 #endif 718 #ifdef Pack_32 719 if ( (y = d1) !=0) { 720 if ( (k = lo0bits(&y)) !=0) { 721 x[0] = y | z << 32 - k; 722 z >>= k; 723 } 724 else 725 x[0] = y; 726 i = b->wds = (x[1] = z) !=0 ? 2 : 1; 727 } 728 else { 729 #ifdef DEBUG 730 if (!z) 731 Bug("Zero passed to d2b"); 732 #endif 733 k = lo0bits(&z); 734 x[0] = z; 735 i = b->wds = 1; 736 k += 32; 737 } 738 #else 739 if ( (y = d1) !=0) { 740 if ( (k = lo0bits(&y)) !=0) 741 if (k >= 16) { 742 x[0] = y | z << 32 - k & 0xffff; 743 x[1] = z >> k - 16 & 0xffff; 744 x[2] = z >> k; 745 i = 2; 746 } 747 else { 748 x[0] = y & 0xffff; 749 x[1] = y >> 16 | z << 16 - k & 0xffff; 750 x[2] = z >> k & 0xffff; 751 x[3] = z >> k+16; 752 i = 3; 753 } 754 else { 755 x[0] = y & 0xffff; 756 x[1] = y >> 16; 757 x[2] = z & 0xffff; 758 x[3] = z >> 16; 759 i = 3; 760 } 761 } 762 else { 763 #ifdef DEBUG 764 if (!z) 765 Bug("Zero passed to d2b"); 766 #endif 767 k = lo0bits(&z); 768 if (k >= 16) { 769 x[0] = z; 770 i = 0; 771 } 772 else { 773 x[0] = z & 0xffff; 774 x[1] = z >> 16; 775 i = 1; 776 } 777 k += 32; 778 } 779 while(!x[i]) 780 --i; 781 b->wds = i + 1; 782 #endif 783 #ifndef Sudden_Underflow 784 if (de) { 785 #endif 786 #ifdef IBM 787 *e = (de - Bias - (P-1) << 2) + k; 788 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 789 #else 790 *e = de - Bias - (P-1) + k; 791 *bits = P - k; 792 #endif 793 #ifndef Sudden_Underflow 794 } 795 else { 796 *e = de - Bias - (P-1) + 1 + k; 797 #ifdef Pack_32 798 *bits = 32*i - hi0bits(x[i-1]); 799 #else 800 *bits = (i+2)*16 - hi0bits(x[i]); 801 #endif 802 } 803 #endif 804 return b; 805 } 806 #undef d0 807 #undef d1 808 809 CONST double 810 #ifdef IEEE_Arith 811 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 812 CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 813 }; 814 #else 815 #ifdef IBM 816 bigtens[] = { 1e16, 1e32, 1e64 }; 817 CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 818 #else 819 bigtens[] = { 1e16, 1e32 }; 820 CONST double tinytens[] = { 1e-16, 1e-32 }; 821 #endif 822 #endif 823 824 CONST double 825 tens[] = { 826 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 827 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 828 1e20, 1e21, 1e22 829 #ifdef VAX 830 , 1e23, 1e24 831 #endif 832 }; 833 834 char * 835 #ifdef KR_headers 836 strcp_D2A(a, b) char *a; char *b; 837 #else 838 strcp_D2A(char *a, CONST char *b) 839 #endif 840 { 841 while(*a = *b++) 842 a++; 843 return a; 844 } 845 846 #ifdef NO_STRING_H 847 848 Char * 849 #ifdef KR_headers 850 memcpy_D2A(a, b, len) Char *a; Char *b; size_t len; 851 #else 852 memcpy_D2A(void *a1, void *b1, size_t len) 853 #endif 854 { 855 register char *a = (char*)a1, *ae = a + len; 856 register char *b = (char*)b1, *a0 = a; 857 while(a < ae) 858 *a++ = *b++; 859 return a0; 860 } 861 862 #endif /* NO_STRING_H */ 863