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