1 #if defined(LIBC_SCCS) && !defined(lint) 2 static char sccsid[] = "@(#)strtod.c 5.2 (Berkeley) 02/20/93"; 3 #endif /* LIBC_SCCS and not lint */ 4 5 /**************************************************************** 6 * 7 * The author of this software is David M. Gay. 8 * 9 * Copyright (c) 1991 by AT&T. 10 * 11 * Permission to use, copy, modify, and distribute this software for any 12 * purpose without fee is hereby granted, provided that this entire notice 13 * is included in all copies of any software which is or includes a copy 14 * or modification of this software and in all copies of the supporting 15 * documentation for such software. 16 * 17 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 18 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY 19 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 20 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 21 * 22 ***************************************************************/ 23 24 /* Please send bug reports to 25 David M. Gay 26 AT&T Bell Laboratories, Room 2C-463 27 600 Mountain Avenue 28 Murray Hill, NJ 07974-2070 29 U.S.A. 30 dmg@research.att.com or research!dmg 31 */ 32 33 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines. 34 * 35 * This strtod returns a nearest machine number to the input decimal 36 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 37 * broken by the IEEE round-even rule. Otherwise ties are broken by 38 * biased rounding (add half and chop). 39 * 40 * Inspired loosely by William D. Clinger's paper "How to Read Floating 41 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 42 * 43 * Modifications: 44 * 45 * 1. We only require IEEE, IBM, or VAX double-precision 46 * arithmetic (not IEEE double-extended). 47 * 2. We get by with floating-point arithmetic in a case that 48 * Clinger missed -- when we're computing d * 10^n 49 * for a small integer d and the integer n is not too 50 * much larger than 22 (the maximum integer k for which 51 * we can represent 10^k exactly), we may be able to 52 * compute (d*10^k) * 10^(e-k) with just one roundoff. 53 * 3. Rather than a bit-at-a-time adjustment of the binary 54 * result in the hard case, we use floating-point 55 * arithmetic to determine the adjustment to within 56 * one bit; only in really hard cases do we need to 57 * compute a second residual. 58 * 4. Because of 3., we don't need a large table of powers of 10 59 * for ten-to-e (just some small tables, e.g. of 10^k 60 * for 0 <= k <= 22). 61 */ 62 63 /* 64 * #define IEEE_8087 for IEEE-arithmetic machines where the least 65 * significant byte has the lowest address. 66 * #define IEEE_MC68k for IEEE-arithmetic machines where the most 67 * significant byte has the lowest address. 68 * #define Sudden_Underflow for IEEE-format machines without gradual 69 * underflow (i.e., that flush to zero on underflow). 70 * #define IBM for IBM mainframe-style floating-point arithmetic. 71 * #define VAX for VAX-style floating-point arithmetic. 72 * #define Unsigned_Shifts if >> does treats its left operand as unsigned. 73 * #define No_leftright to omit left-right logic in fast floating-point 74 * computation of dtoa. 75 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3. 76 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 77 * that use extended-precision instructions to compute rounded 78 * products and quotients) with IBM. 79 * #define ROUND_BIASED for IEEE-format with biased rounding. 80 * #define Inaccurate_Divide for IEEE-format with correctly rounded 81 * products but inaccurate quotients, e.g., for Intel i860. 82 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision 83 * integer arithmetic. Whether this speeds things up or slows things 84 * down depends on the machine and the number being converted. 85 * #define KR_headers for old-style C function headers. 86 * #define Bad_float_h if your system lacks a float.h or if it does not 87 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, 88 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. 89 */ 90 91 #if defined(i386) || defined(mips) && defined(MIPSEL) 92 #define IEEE_8087 93 #else 94 #define IEEE_MC68k 95 #endif 96 97 #ifdef DEBUG 98 #include "stdio.h" 99 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 100 #endif 101 102 #ifdef __cplusplus 103 #include "malloc.h" 104 #include "memory.h" 105 #else 106 #ifndef KR_headers 107 #include "stdlib.h" 108 #include "string.h" 109 #else 110 #include "malloc.h" 111 #include "memory.h" 112 #endif 113 #endif 114 115 #include "errno.h" 116 #ifdef Bad_float_h 117 #undef __STDC__ 118 #ifdef IEEE_MC68k 119 #define IEEE_ARITHMETIC 120 #endif 121 #ifdef IEEE_8087 122 #define IEEE_ARITHMETIC 123 #endif 124 #ifdef IEEE_ARITHMETIC 125 #define DBL_DIG 15 126 #define DBL_MAX_10_EXP 308 127 #define DBL_MAX_EXP 1024 128 #define FLT_RADIX 2 129 #define FLT_ROUNDS 1 130 #define DBL_MAX 1.7976931348623157e+308 131 #endif 132 133 #ifdef IBM 134 #define DBL_DIG 16 135 #define DBL_MAX_10_EXP 75 136 #define DBL_MAX_EXP 63 137 #define FLT_RADIX 16 138 #define FLT_ROUNDS 0 139 #define DBL_MAX 7.2370055773322621e+75 140 #endif 141 142 #ifdef VAX 143 #define DBL_DIG 16 144 #define DBL_MAX_10_EXP 38 145 #define DBL_MAX_EXP 127 146 #define FLT_RADIX 2 147 #define FLT_ROUNDS 1 148 #define DBL_MAX 1.7014118346046923e+38 149 #endif 150 151 #ifndef LONG_MAX 152 #define LONG_MAX 2147483647 153 #endif 154 #else 155 #include "float.h" 156 #endif 157 #ifndef __MATH_H__ 158 #include "math.h" 159 #endif 160 161 #ifdef __cplusplus 162 extern "C" { 163 #endif 164 165 #ifndef CONST 166 #ifdef KR_headers 167 #define CONST /* blank */ 168 #else 169 #define CONST const 170 #endif 171 #endif 172 173 #ifdef Unsigned_Shifts 174 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000; 175 #else 176 #define Sign_Extend(a,b) /*no-op*/ 177 #endif 178 179 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 180 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. 181 #endif 182 183 #ifdef IEEE_8087 184 #define word0(x) ((unsigned long *)&x)[1] 185 #define word1(x) ((unsigned long *)&x)[0] 186 #else 187 #define word0(x) ((unsigned long *)&x)[0] 188 #define word1(x) ((unsigned long *)&x)[1] 189 #endif 190 191 /* The following definition of Storeinc is appropriate for MIPS processors. 192 * An alternative that might be better on some machines is 193 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 194 */ 195 #if defined(IEEE_8087) + defined(VAX) 196 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ 197 ((unsigned short *)a)[0] = (unsigned short)c, a++) 198 #else 199 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ 200 ((unsigned short *)a)[1] = (unsigned short)c, a++) 201 #endif 202 203 /* #define P DBL_MANT_DIG */ 204 /* Ten_pmax = floor(P*log(2)/log(5)) */ 205 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 206 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 207 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 208 209 #if defined(IEEE_8087) + defined(IEEE_MC68k) 210 #define Exp_shift 20 211 #define Exp_shift1 20 212 #define Exp_msk1 0x100000 213 #define Exp_msk11 0x100000 214 #define Exp_mask 0x7ff00000 215 #define P 53 216 #define Bias 1023 217 #define IEEE_Arith 218 #define Emin (-1022) 219 #define Exp_1 0x3ff00000 220 #define Exp_11 0x3ff00000 221 #define Ebits 11 222 #define Frac_mask 0xfffff 223 #define Frac_mask1 0xfffff 224 #define Ten_pmax 22 225 #define Bletch 0x10 226 #define Bndry_mask 0xfffff 227 #define Bndry_mask1 0xfffff 228 #define LSB 1 229 #define Sign_bit 0x80000000 230 #define Log2P 1 231 #define Tiny0 0 232 #define Tiny1 1 233 #define Quick_max 14 234 #define Int_max 14 235 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */ 236 #else 237 #undef Sudden_Underflow 238 #define Sudden_Underflow 239 #ifdef IBM 240 #define Exp_shift 24 241 #define Exp_shift1 24 242 #define Exp_msk1 0x1000000 243 #define Exp_msk11 0x1000000 244 #define Exp_mask 0x7f000000 245 #define P 14 246 #define Bias 65 247 #define Exp_1 0x41000000 248 #define Exp_11 0x41000000 249 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 250 #define Frac_mask 0xffffff 251 #define Frac_mask1 0xffffff 252 #define Bletch 4 253 #define Ten_pmax 22 254 #define Bndry_mask 0xefffff 255 #define Bndry_mask1 0xffffff 256 #define LSB 1 257 #define Sign_bit 0x80000000 258 #define Log2P 4 259 #define Tiny0 0x100000 260 #define Tiny1 0 261 #define Quick_max 14 262 #define Int_max 15 263 #else /* VAX */ 264 #define Exp_shift 23 265 #define Exp_shift1 7 266 #define Exp_msk1 0x80 267 #define Exp_msk11 0x800000 268 #define Exp_mask 0x7f80 269 #define P 56 270 #define Bias 129 271 #define Exp_1 0x40800000 272 #define Exp_11 0x4080 273 #define Ebits 8 274 #define Frac_mask 0x7fffff 275 #define Frac_mask1 0xffff007f 276 #define Ten_pmax 24 277 #define Bletch 2 278 #define Bndry_mask 0xffff007f 279 #define Bndry_mask1 0xffff007f 280 #define LSB 0x10000 281 #define Sign_bit 0x8000 282 #define Log2P 1 283 #define Tiny0 0x80 284 #define Tiny1 0 285 #define Quick_max 15 286 #define Int_max 15 287 #endif 288 #endif 289 290 #ifndef IEEE_Arith 291 #define ROUND_BIASED 292 #endif 293 294 #ifdef RND_PRODQUOT 295 #define rounded_product(a,b) a = rnd_prod(a, b) 296 #define rounded_quotient(a,b) a = rnd_quot(a, b) 297 #ifdef KR_headers 298 extern double rnd_prod(), rnd_quot(); 299 #else 300 extern double rnd_prod(double, double), rnd_quot(double, double); 301 #endif 302 #else 303 #define rounded_product(a,b) a *= b 304 #define rounded_quotient(a,b) a /= b 305 #endif 306 307 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 308 #define Big1 0xffffffff 309 310 #ifndef Just_16 311 /* When Pack_32 is not defined, we store 16 bits per 32-bit long. 312 * This makes some inner loops simpler and sometimes saves work 313 * during multiplications, but it often seems to make things slightly 314 * slower. Hence the default is now to store 32 bits per long. 315 */ 316 #ifndef Pack_32 317 #define Pack_32 318 #endif 319 #endif 320 321 #define Kmax 15 322 323 #ifdef __cplusplus 324 extern "C" double strtod(const char *s00, char **se); 325 extern "C" char *dtoa(double d, int mode, int ndigits, 326 int *decpt, int *sign, char **rve); 327 #endif 328 329 struct 330 Bigint { 331 struct Bigint *next; 332 int k, maxwds, sign, wds; 333 unsigned long x[1]; 334 }; 335 336 typedef struct Bigint Bigint; 337 338 static Bigint *freelist[Kmax+1]; 339 340 static Bigint * 341 Balloc 342 #ifdef KR_headers 343 (k) int k; 344 #else 345 (int k) 346 #endif 347 { 348 int x; 349 Bigint *rv; 350 351 if (rv = freelist[k]) { 352 freelist[k] = rv->next; 353 } else { 354 x = 1 << k; 355 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(long)); 356 rv->k = k; 357 rv->maxwds = x; 358 } 359 rv->sign = rv->wds = 0; 360 return rv; 361 } 362 363 static void 364 Bfree 365 #ifdef KR_headers 366 (v) Bigint *v; 367 #else 368 (Bigint *v) 369 #endif 370 { 371 if (v) { 372 v->next = freelist[v->k]; 373 freelist[v->k] = v; 374 } 375 } 376 377 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 378 y->wds*sizeof(long) + 2*sizeof(int)) 379 380 static Bigint * 381 multadd 382 #ifdef KR_headers 383 (b, m, a) Bigint *b; int m, a; 384 #else 385 (Bigint *b, int m, int a) /* multiply by m and add a */ 386 #endif 387 { 388 int i, wds; 389 unsigned long *x, y; 390 #ifdef Pack_32 391 unsigned long xi, z; 392 #endif 393 Bigint *b1; 394 395 wds = b->wds; 396 x = b->x; 397 i = 0; 398 do { 399 #ifdef Pack_32 400 xi = *x; 401 y = (xi & 0xffff) * m + a; 402 z = (xi >> 16) * m + (y >> 16); 403 a = (int)(z >> 16); 404 *x++ = (z << 16) + (y & 0xffff); 405 #else 406 y = *x * m + a; 407 a = (int)(y >> 16); 408 *x++ = y & 0xffff; 409 #endif 410 } while (++i < wds); 411 if (a) { 412 if (wds >= b->maxwds) { 413 b1 = Balloc(b->k+1); 414 Bcopy(b1, b); 415 Bfree(b); 416 b = b1; 417 } 418 b->x[wds++] = a; 419 b->wds = wds; 420 } 421 return b; 422 } 423 424 static Bigint * 425 s2b 426 #ifdef KR_headers 427 (s, nd0, nd, y9) CONST char *s; int nd0, nd; unsigned long y9; 428 #else 429 (CONST char *s, int nd0, int nd, unsigned long y9) 430 #endif 431 { 432 Bigint *b; 433 int i, k; 434 long x, y; 435 436 x = (nd + 8) / 9; 437 for (k = 0, y = 1; x > y; y <<= 1, k++) ; 438 #ifdef Pack_32 439 b = Balloc(k); 440 b->x[0] = y9; 441 b->wds = 1; 442 #else 443 b = Balloc(k+1); 444 b->x[0] = y9 & 0xffff; 445 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; 446 #endif 447 448 i = 9; 449 if (9 < nd0) { 450 s += 9; 451 do 452 b = multadd(b, 10, *s++ - '0'); 453 while (++i < nd0); 454 s++; 455 } else 456 s += 10; 457 for (; i < nd; i++) 458 b = multadd(b, 10, *s++ - '0'); 459 return b; 460 } 461 462 static int 463 hi0bits 464 #ifdef KR_headers 465 (x) register unsigned long x; 466 #else 467 (register unsigned long x) 468 #endif 469 { 470 register int k = 0; 471 472 if (!(x & 0xffff0000)) { 473 k = 16; 474 x <<= 16; 475 } 476 if (!(x & 0xff000000)) { 477 k += 8; 478 x <<= 8; 479 } 480 if (!(x & 0xf0000000)) { 481 k += 4; 482 x <<= 4; 483 } 484 if (!(x & 0xc0000000)) { 485 k += 2; 486 x <<= 2; 487 } 488 if (!(x & 0x80000000)) { 489 k++; 490 if (!(x & 0x40000000)) 491 return 32; 492 } 493 return k; 494 } 495 496 static int 497 lo0bits 498 #ifdef KR_headers 499 (y) unsigned long *y; 500 #else 501 (unsigned long *y) 502 #endif 503 { 504 register int k; 505 register unsigned long x = *y; 506 507 if (x & 7) { 508 if (x & 1) 509 return 0; 510 if (x & 2) { 511 *y = x >> 1; 512 return 1; 513 } 514 *y = x >> 2; 515 return 2; 516 } 517 k = 0; 518 if (!(x & 0xffff)) { 519 k = 16; 520 x >>= 16; 521 } 522 if (!(x & 0xff)) { 523 k += 8; 524 x >>= 8; 525 } 526 if (!(x & 0xf)) { 527 k += 4; 528 x >>= 4; 529 } 530 if (!(x & 0x3)) { 531 k += 2; 532 x >>= 2; 533 } 534 if (!(x & 1)) { 535 k++; 536 x >>= 1; 537 if (!x & 1) 538 return 32; 539 } 540 *y = x; 541 return k; 542 } 543 544 static Bigint * 545 i2b 546 #ifdef KR_headers 547 (i) int i; 548 #else 549 (int i) 550 #endif 551 { 552 Bigint *b; 553 554 b = Balloc(1); 555 b->x[0] = i; 556 b->wds = 1; 557 return b; 558 } 559 560 static Bigint * 561 mult 562 #ifdef KR_headers 563 (a, b) Bigint *a, *b; 564 #else 565 (Bigint *a, Bigint *b) 566 #endif 567 { 568 Bigint *c; 569 int k, wa, wb, wc; 570 unsigned long carry, y, z; 571 unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 572 #ifdef Pack_32 573 unsigned long z2; 574 #endif 575 576 if (a->wds < b->wds) { 577 c = a; 578 a = b; 579 b = c; 580 } 581 k = a->k; 582 wa = a->wds; 583 wb = b->wds; 584 wc = wa + wb; 585 if (wc > a->maxwds) 586 k++; 587 c = Balloc(k); 588 for (x = c->x, xa = x + wc; x < xa; x++) 589 *x = 0; 590 xa = a->x; 591 xae = xa + wa; 592 xb = b->x; 593 xbe = xb + wb; 594 xc0 = c->x; 595 #ifdef Pack_32 596 for (; xb < xbe; xb++, xc0++) { 597 if (y = *xb & 0xffff) { 598 x = xa; 599 xc = xc0; 600 carry = 0; 601 do { 602 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 603 carry = z >> 16; 604 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 605 carry = z2 >> 16; 606 Storeinc(xc, z2, z); 607 } while (x < xae); 608 *xc = carry; 609 } 610 if (y = *xb >> 16) { 611 x = xa; 612 xc = xc0; 613 carry = 0; 614 z2 = *xc; 615 do { 616 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 617 carry = z >> 16; 618 Storeinc(xc, z, z2); 619 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 620 carry = z2 >> 16; 621 } while (x < xae); 622 *xc = z2; 623 } 624 } 625 #else 626 for (; xb < xbe; xc0++) { 627 if (y = *xb++) { 628 x = xa; 629 xc = xc0; 630 carry = 0; 631 do { 632 z = *x++ * y + *xc + carry; 633 carry = z >> 16; 634 *xc++ = z & 0xffff; 635 } while (x < xae); 636 *xc = carry; 637 } 638 } 639 #endif 640 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 641 c->wds = wc; 642 return c; 643 } 644 645 static Bigint *p5s; 646 647 static Bigint * 648 pow5mult 649 #ifdef KR_headers 650 (b, k) Bigint *b; int k; 651 #else 652 (Bigint *b, int k) 653 #endif 654 { 655 Bigint *b1, *p5, *p51; 656 int i; 657 static int p05[3] = { 5, 25, 125 }; 658 659 if (i = k & 3) 660 b = multadd(b, p05[i-1], 0); 661 662 if (!(k >>= 2)) 663 return b; 664 if (!(p5 = p5s)) { 665 /* first time */ 666 p5 = p5s = i2b(625); 667 p5->next = 0; 668 } 669 for (;;) { 670 if (k & 1) { 671 b1 = mult(b, p5); 672 Bfree(b); 673 b = b1; 674 } 675 if (!(k >>= 1)) 676 break; 677 if (!(p51 = p5->next)) { 678 p51 = p5->next = mult(p5,p5); 679 p51->next = 0; 680 } 681 p5 = p51; 682 } 683 return b; 684 } 685 686 static Bigint * 687 lshift 688 #ifdef KR_headers 689 (b, k) Bigint *b; int k; 690 #else 691 (Bigint *b, int k) 692 #endif 693 { 694 int i, k1, n, n1; 695 Bigint *b1; 696 unsigned long *x, *x1, *xe, z; 697 698 #ifdef Pack_32 699 n = k >> 5; 700 #else 701 n = k >> 4; 702 #endif 703 k1 = b->k; 704 n1 = n + b->wds + 1; 705 for (i = b->maxwds; n1 > i; i <<= 1) 706 k1++; 707 b1 = Balloc(k1); 708 x1 = b1->x; 709 for (i = 0; i < n; i++) 710 *x1++ = 0; 711 x = b->x; 712 xe = x + b->wds; 713 #ifdef Pack_32 714 if (k &= 0x1f) { 715 k1 = 32 - k; 716 z = 0; 717 do { 718 *x1++ = *x << k | z; 719 z = *x++ >> k1; 720 } while (x < xe); 721 if (*x1 = z) 722 ++n1; 723 } 724 #else 725 if (k &= 0xf) { 726 k1 = 16 - k; 727 z = 0; 728 do { 729 *x1++ = *x << k & 0xffff | z; 730 z = *x++ >> k1; 731 } while (x < xe); 732 if (*x1 = z) 733 ++n1; 734 } 735 #endif 736 else 737 do 738 *x1++ = *x++; 739 while (x < xe); 740 b1->wds = n1 - 1; 741 Bfree(b); 742 return b1; 743 } 744 745 static int 746 cmp 747 #ifdef KR_headers 748 (a, b) Bigint *a, *b; 749 #else 750 (Bigint *a, Bigint *b) 751 #endif 752 { 753 unsigned long *xa, *xa0, *xb, *xb0; 754 int i, j; 755 756 i = a->wds; 757 j = b->wds; 758 #ifdef DEBUG 759 if (i > 1 && !a->x[i-1]) 760 Bug("cmp called with a->x[a->wds-1] == 0"); 761 if (j > 1 && !b->x[j-1]) 762 Bug("cmp called with b->x[b->wds-1] == 0"); 763 #endif 764 if (i -= j) 765 return i; 766 xa0 = a->x; 767 xa = xa0 + j; 768 xb0 = b->x; 769 xb = xb0 + j; 770 for (;;) { 771 if (*--xa != *--xb) 772 return *xa < *xb ? -1 : 1; 773 if (xa <= xa0) 774 break; 775 } 776 return 0; 777 } 778 779 static Bigint * 780 diff 781 #ifdef KR_headers 782 (a, b) Bigint *a, *b; 783 #else 784 (Bigint *a, Bigint *b) 785 #endif 786 { 787 Bigint *c; 788 int i, wa, wb; 789 long borrow, y; /* We need signed shifts here. */ 790 unsigned long *xa, *xae, *xb, *xbe, *xc; 791 #ifdef Pack_32 792 long z; 793 #endif 794 795 i = cmp(a,b); 796 if (!i) { 797 c = Balloc(0); 798 c->wds = 1; 799 c->x[0] = 0; 800 return c; 801 } 802 if (i < 0) { 803 c = a; 804 a = b; 805 b = c; 806 i = 1; 807 } else 808 i = 0; 809 c = Balloc(a->k); 810 c->sign = i; 811 wa = a->wds; 812 xa = a->x; 813 xae = xa + wa; 814 wb = b->wds; 815 xb = b->x; 816 xbe = xb + wb; 817 xc = c->x; 818 borrow = 0; 819 #ifdef Pack_32 820 do { 821 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow; 822 borrow = y >> 16; 823 Sign_Extend(borrow, y); 824 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow; 825 borrow = z >> 16; 826 Sign_Extend(borrow, z); 827 Storeinc(xc, z, y); 828 } while (xb < xbe); 829 while (xa < xae) { 830 y = (*xa & 0xffff) + borrow; 831 borrow = y >> 16; 832 Sign_Extend(borrow, y); 833 z = (*xa++ >> 16) + borrow; 834 borrow = z >> 16; 835 Sign_Extend(borrow, z); 836 Storeinc(xc, z, y); 837 } 838 #else 839 do { 840 y = *xa++ - *xb++ + borrow; 841 borrow = y >> 16; 842 Sign_Extend(borrow, y); 843 *xc++ = y & 0xffff; 844 } while (xb < xbe); 845 while (xa < xae) { 846 y = *xa++ + borrow; 847 borrow = y >> 16; 848 Sign_Extend(borrow, y); 849 *xc++ = y & 0xffff; 850 } 851 #endif 852 while (!*--xc) 853 wa--; 854 c->wds = wa; 855 return c; 856 } 857 858 static double 859 ulp 860 #ifdef KR_headers 861 (x) double x; 862 #else 863 (double x) 864 #endif 865 { 866 register long L; 867 double a; 868 869 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 870 #ifndef Sudden_Underflow 871 if (L > 0) { 872 #endif 873 #ifdef IBM 874 L |= Exp_msk1 >> 4; 875 #endif 876 word0(a) = L; 877 word1(a) = 0; 878 #ifndef Sudden_Underflow 879 } else { 880 L = -L >> Exp_shift; 881 if (L < Exp_shift) { 882 word0(a) = 0x80000 >> L; 883 word1(a) = 0; 884 } else { 885 word0(a) = 0; 886 L -= Exp_shift; 887 word1(a) = L >= 31 ? 1 : 1 << 31 - L; 888 } 889 } 890 #endif 891 return a; 892 } 893 894 static double 895 b2d 896 #ifdef KR_headers 897 (a, e) Bigint *a; int *e; 898 #else 899 (Bigint *a, int *e) 900 #endif 901 { 902 unsigned long *xa, *xa0, w, y, z; 903 int k; 904 double d; 905 #ifdef VAX 906 unsigned long d0, d1; 907 #else 908 #define d0 word0(d) 909 #define d1 word1(d) 910 #endif 911 912 xa0 = a->x; 913 xa = xa0 + a->wds; 914 y = *--xa; 915 #ifdef DEBUG 916 if (!y) Bug("zero y in b2d"); 917 #endif 918 k = hi0bits(y); 919 *e = 32 - k; 920 #ifdef Pack_32 921 if (k < Ebits) { 922 d0 = Exp_1 | y >> Ebits - k; 923 w = xa > xa0 ? *--xa : 0; 924 d1 = y << (32-Ebits) + k | w >> Ebits - k; 925 goto ret_d; 926 } 927 z = xa > xa0 ? *--xa : 0; 928 if (k -= Ebits) { 929 d0 = Exp_1 | y << k | z >> 32 - k; 930 y = xa > xa0 ? *--xa : 0; 931 d1 = z << k | y >> 32 - k; 932 } else { 933 d0 = Exp_1 | y; 934 d1 = z; 935 } 936 #else 937 if (k < Ebits + 16) { 938 z = xa > xa0 ? *--xa : 0; 939 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 940 w = xa > xa0 ? *--xa : 0; 941 y = xa > xa0 ? *--xa : 0; 942 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 943 goto ret_d; 944 } 945 z = xa > xa0 ? *--xa : 0; 946 w = xa > xa0 ? *--xa : 0; 947 k -= Ebits + 16; 948 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 949 y = xa > xa0 ? *--xa : 0; 950 d1 = w << k + 16 | y << k; 951 #endif 952 ret_d: 953 #ifdef VAX 954 word0(d) = d0 >> 16 | d0 << 16; 955 word1(d) = d1 >> 16 | d1 << 16; 956 #else 957 #undef d0 958 #undef d1 959 #endif 960 return d; 961 } 962 963 static Bigint * 964 d2b 965 #ifdef KR_headers 966 (d, e, bits) double d; int *e, *bits; 967 #else 968 (double d, int *e, int *bits) 969 #endif 970 { 971 Bigint *b; 972 int de, i, k; 973 unsigned long *x, y, z; 974 #ifdef VAX 975 unsigned long d0, d1; 976 d0 = word0(d) >> 16 | word0(d) << 16; 977 d1 = word1(d) >> 16 | word1(d) << 16; 978 #else 979 #define d0 word0(d) 980 #define d1 word1(d) 981 #endif 982 983 #ifdef Pack_32 984 b = Balloc(1); 985 #else 986 b = Balloc(2); 987 #endif 988 x = b->x; 989 990 z = d0 & Frac_mask; 991 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 992 #ifdef Sudden_Underflow 993 de = (int)(d0 >> Exp_shift); 994 #ifndef IBM 995 z |= Exp_msk11; 996 #endif 997 #else 998 if (de = (int)(d0 >> Exp_shift)) 999 z |= Exp_msk1; 1000 #endif 1001 #ifdef Pack_32 1002 if (y = d1) { 1003 if (k = lo0bits(&y)) { 1004 x[0] = y | z << 32 - k; 1005 z >>= k; 1006 } 1007 else 1008 x[0] = y; 1009 i = b->wds = (x[1] = z) ? 2 : 1; 1010 } else { 1011 #ifdef DEBUG 1012 if (!z) 1013 Bug("Zero passed to d2b"); 1014 #endif 1015 k = lo0bits(&z); 1016 x[0] = z; 1017 i = b->wds = 1; 1018 k += 32; 1019 } 1020 #else 1021 if (y = d1) { 1022 if (k = lo0bits(&y)) 1023 if (k >= 16) { 1024 x[0] = y | z << 32 - k & 0xffff; 1025 x[1] = z >> k - 16 & 0xffff; 1026 x[2] = z >> k; 1027 i = 2; 1028 } else { 1029 x[0] = y & 0xffff; 1030 x[1] = y >> 16 | z << 16 - k & 0xffff; 1031 x[2] = z >> k & 0xffff; 1032 x[3] = z >> k+16; 1033 i = 3; 1034 } 1035 else { 1036 x[0] = y & 0xffff; 1037 x[1] = y >> 16; 1038 x[2] = z & 0xffff; 1039 x[3] = z >> 16; 1040 i = 3; 1041 } 1042 } else { 1043 #ifdef DEBUG 1044 if (!z) 1045 Bug("Zero passed to d2b"); 1046 #endif 1047 k = lo0bits(&z); 1048 if (k >= 16) { 1049 x[0] = z; 1050 i = 0; 1051 } else { 1052 x[0] = z & 0xffff; 1053 x[1] = z >> 16; 1054 i = 1; 1055 } 1056 k += 32; 1057 } 1058 while (!x[i]) 1059 --i; 1060 b->wds = i + 1; 1061 #endif 1062 #ifndef Sudden_Underflow 1063 if (de) { 1064 #endif 1065 #ifdef IBM 1066 *e = (de - Bias - (P-1) << 2) + k; 1067 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 1068 #else 1069 *e = de - Bias - (P-1) + k; 1070 *bits = P - k; 1071 #endif 1072 #ifndef Sudden_Underflow 1073 } else { 1074 *e = de - Bias - (P-1) + 1 + k; 1075 #ifdef Pack_32 1076 *bits = 32*i - hi0bits(x[i-1]); 1077 #else 1078 *bits = (i+2)*16 - hi0bits(x[i]); 1079 #endif 1080 } 1081 #endif 1082 return b; 1083 } 1084 #undef d0 1085 #undef d1 1086 1087 static double 1088 ratio 1089 #ifdef KR_headers 1090 (a, b) Bigint *a, *b; 1091 #else 1092 (Bigint *a, Bigint *b) 1093 #endif 1094 { 1095 double da, db; 1096 int k, ka, kb; 1097 1098 da = b2d(a, &ka); 1099 db = b2d(b, &kb); 1100 #ifdef Pack_32 1101 k = ka - kb + 32*(a->wds - b->wds); 1102 #else 1103 k = ka - kb + 16*(a->wds - b->wds); 1104 #endif 1105 #ifdef IBM 1106 if (k > 0) { 1107 word0(da) += (k >> 2)*Exp_msk1; 1108 if (k &= 3) 1109 da *= 1 << k; 1110 } else { 1111 k = -k; 1112 word0(db) += (k >> 2)*Exp_msk1; 1113 if (k &= 3) 1114 db *= 1 << k; 1115 } 1116 #else 1117 if (k > 0) 1118 word0(da) += k*Exp_msk1; 1119 else { 1120 k = -k; 1121 word0(db) += k*Exp_msk1; 1122 } 1123 #endif 1124 return da / db; 1125 } 1126 1127 static double 1128 tens[] = { 1129 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1130 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1131 1e20, 1e21, 1e22 1132 #ifdef VAX 1133 , 1e23, 1e24 1134 #endif 1135 }; 1136 1137 static double 1138 #ifdef IEEE_Arith 1139 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 1140 static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 }; 1141 #define n_bigtens 5 1142 #else 1143 #ifdef IBM 1144 bigtens[] = { 1e16, 1e32, 1e64 }; 1145 static double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 1146 #define n_bigtens 3 1147 #else 1148 bigtens[] = { 1e16, 1e32 }; 1149 static double tinytens[] = { 1e-16, 1e-32 }; 1150 #define n_bigtens 2 1151 #endif 1152 #endif 1153 1154 double 1155 strtod 1156 #ifdef KR_headers 1157 (s00, se) CONST char *s00; char **se; 1158 #else 1159 (CONST char *s00, char **se) 1160 #endif 1161 { 1162 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, 1163 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 1164 CONST char *s, *s0, *s1; 1165 double aadj, aadj1, adj, rv, rv0; 1166 long L; 1167 unsigned long y, z; 1168 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; 1169 sign = nz0 = nz = 0; 1170 rv = 0.; 1171 for (s = s00;;s++) switch(*s) { 1172 case '-': 1173 sign = 1; 1174 /* no break */ 1175 case '+': 1176 if (*++s) 1177 goto break2; 1178 /* no break */ 1179 case 0: 1180 s = s00; 1181 goto ret; 1182 case '\t': 1183 case '\n': 1184 case '\v': 1185 case '\f': 1186 case '\r': 1187 case ' ': 1188 continue; 1189 default: 1190 goto break2; 1191 } 1192 break2: 1193 if (*s == '0') { 1194 nz0 = 1; 1195 while (*++s == '0') ; 1196 if (!*s) 1197 goto ret; 1198 } 1199 s0 = s; 1200 y = z = 0; 1201 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 1202 if (nd < 9) 1203 y = 10*y + c - '0'; 1204 else if (nd < 16) 1205 z = 10*z + c - '0'; 1206 nd0 = nd; 1207 if (c == '.') { 1208 c = *++s; 1209 if (!nd) { 1210 for (; c == '0'; c = *++s) 1211 nz++; 1212 if (c > '0' && c <= '9') { 1213 s0 = s; 1214 nf += nz; 1215 nz = 0; 1216 goto have_dig; 1217 } 1218 goto dig_done; 1219 } 1220 for (; c >= '0' && c <= '9'; c = *++s) { 1221 have_dig: 1222 nz++; 1223 if (c -= '0') { 1224 nf += nz; 1225 for (i = 1; i < nz; i++) 1226 if (nd++ < 9) 1227 y *= 10; 1228 else if (nd <= DBL_DIG + 1) 1229 z *= 10; 1230 if (nd++ < 9) 1231 y = 10*y + c; 1232 else if (nd <= DBL_DIG + 1) 1233 z = 10*z + c; 1234 nz = 0; 1235 } 1236 } 1237 } 1238 dig_done: 1239 e = 0; 1240 if (c == 'e' || c == 'E') { 1241 if (!nd && !nz && !nz0) { 1242 s = s00; 1243 goto ret; 1244 } 1245 s00 = s; 1246 esign = 0; 1247 switch(c = *++s) { 1248 case '-': 1249 esign = 1; 1250 case '+': 1251 c = *++s; 1252 } 1253 if (c >= '0' && c <= '9') { 1254 while (c == '0') 1255 c = *++s; 1256 if (c > '0' && c <= '9') { 1257 L = c - '0'; 1258 s1 = s; 1259 while ((c = *++s) >= '0' && c <= '9') 1260 L = 10*L + c - '0'; 1261 if (s - s1 > 8 || L > 19999) 1262 /* Avoid confusion from exponents 1263 * so large that e might overflow. 1264 */ 1265 e = 19999; /* safe for 16 bit ints */ 1266 else 1267 e = (int)L; 1268 if (esign) 1269 e = -e; 1270 } else 1271 e = 0; 1272 } else 1273 s = s00; 1274 } 1275 if (!nd) { 1276 if (!nz && !nz0) 1277 s = s00; 1278 goto ret; 1279 } 1280 e1 = e -= nf; 1281 1282 /* Now we have nd0 digits, starting at s0, followed by a 1283 * decimal point, followed by nd-nd0 digits. The number we're 1284 * after is the integer represented by those digits times 1285 * 10**e */ 1286 1287 if (!nd0) 1288 nd0 = nd; 1289 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 1290 rv = y; 1291 if (k > 9) 1292 rv = tens[k - 9] * rv + z; 1293 if (nd <= DBL_DIG 1294 #ifndef RND_PRODQUOT 1295 && FLT_ROUNDS == 1 1296 #endif 1297 ) { 1298 if (!e) 1299 goto ret; 1300 if (e > 0) { 1301 if (e <= Ten_pmax) { 1302 #ifdef VAX 1303 goto vax_ovfl_check; 1304 #else 1305 /* rv = */ rounded_product(rv, tens[e]); 1306 goto ret; 1307 #endif 1308 } 1309 i = DBL_DIG - nd; 1310 if (e <= Ten_pmax + i) { 1311 /* A fancier test would sometimes let us do 1312 * this for larger i values. 1313 */ 1314 e -= i; 1315 rv *= tens[i]; 1316 #ifdef VAX 1317 /* VAX exponent range is so narrow we must 1318 * worry about overflow here... 1319 */ 1320 vax_ovfl_check: 1321 word0(rv) -= P*Exp_msk1; 1322 /* rv = */ rounded_product(rv, tens[e]); 1323 if ((word0(rv) & Exp_mask) 1324 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 1325 goto ovfl; 1326 word0(rv) += P*Exp_msk1; 1327 #else 1328 /* rv = */ rounded_product(rv, tens[e]); 1329 #endif 1330 goto ret; 1331 } 1332 } 1333 #ifndef Inaccurate_Divide 1334 else if (e >= -Ten_pmax) { 1335 /* rv = */ rounded_quotient(rv, tens[-e]); 1336 goto ret; 1337 } 1338 #endif 1339 } 1340 e1 += nd - k; 1341 1342 /* Get starting approximation = rv * 10**e1 */ 1343 1344 if (e1 > 0) { 1345 if (i = e1 & 15) 1346 rv *= tens[i]; 1347 if (e1 &= ~15) { 1348 if (e1 > DBL_MAX_10_EXP) { 1349 ovfl: 1350 errno = ERANGE; 1351 #ifdef __STDC__ 1352 rv = HUGE_VAL; 1353 #else 1354 /* Can't trust HUGE_VAL */ 1355 #ifdef IEEE_Arith 1356 word0(rv) = Exp_mask; 1357 word1(rv) = 0; 1358 #else 1359 word0(rv) = Big0; 1360 word1(rv) = Big1; 1361 #endif 1362 #endif 1363 goto ret; 1364 } 1365 if (e1 >>= 4) { 1366 for (j = 0; e1 > 1; j++, e1 >>= 1) 1367 if (e1 & 1) 1368 rv *= bigtens[j]; 1369 /* The last multiplication could overflow. */ 1370 word0(rv) -= P*Exp_msk1; 1371 rv *= bigtens[j]; 1372 if ((z = word0(rv) & Exp_mask) 1373 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 1374 goto ovfl; 1375 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 1376 /* set to largest number */ 1377 /* (Can't trust DBL_MAX) */ 1378 word0(rv) = Big0; 1379 word1(rv) = Big1; 1380 } 1381 else 1382 word0(rv) += P*Exp_msk1; 1383 } 1384 } 1385 } else if (e1 < 0) { 1386 e1 = -e1; 1387 if (i = e1 & 15) 1388 rv /= tens[i]; 1389 if (e1 &= ~15) { 1390 e1 >>= 4; 1391 for (j = 0; e1 > 1; j++, e1 >>= 1) 1392 if (e1 & 1) 1393 rv *= tinytens[j]; 1394 /* The last multiplication could underflow. */ 1395 rv0 = rv; 1396 rv *= tinytens[j]; 1397 if (!rv) { 1398 rv = 2.*rv0; 1399 rv *= tinytens[j]; 1400 if (!rv) { 1401 undfl: 1402 rv = 0.; 1403 errno = ERANGE; 1404 goto ret; 1405 } 1406 word0(rv) = Tiny0; 1407 word1(rv) = Tiny1; 1408 /* The refinement below will clean 1409 * this approximation up. 1410 */ 1411 } 1412 } 1413 } 1414 1415 /* Now the hard part -- adjusting rv to the correct value.*/ 1416 1417 /* Put digits into bd: true value = bd * 10^e */ 1418 1419 bd0 = s2b(s0, nd0, nd, y); 1420 1421 for (;;) { 1422 bd = Balloc(bd0->k); 1423 Bcopy(bd, bd0); 1424 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */ 1425 bs = i2b(1); 1426 1427 if (e >= 0) { 1428 bb2 = bb5 = 0; 1429 bd2 = bd5 = e; 1430 } else { 1431 bb2 = bb5 = -e; 1432 bd2 = bd5 = 0; 1433 } 1434 if (bbe >= 0) 1435 bb2 += bbe; 1436 else 1437 bd2 -= bbe; 1438 bs2 = bb2; 1439 #ifdef Sudden_Underflow 1440 #ifdef IBM 1441 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 1442 #else 1443 j = P + 1 - bbbits; 1444 #endif 1445 #else 1446 i = bbe + bbbits - 1; /* logb(rv) */ 1447 if (i < Emin) /* denormal */ 1448 j = bbe + (P-Emin); 1449 else 1450 j = P + 1 - bbbits; 1451 #endif 1452 bb2 += j; 1453 bd2 += j; 1454 i = bb2 < bd2 ? bb2 : bd2; 1455 if (i > bs2) 1456 i = bs2; 1457 if (i > 0) { 1458 bb2 -= i; 1459 bd2 -= i; 1460 bs2 -= i; 1461 } 1462 if (bb5 > 0) { 1463 bs = pow5mult(bs, bb5); 1464 bb1 = mult(bs, bb); 1465 Bfree(bb); 1466 bb = bb1; 1467 } 1468 if (bb2 > 0) 1469 bb = lshift(bb, bb2); 1470 if (bd5 > 0) 1471 bd = pow5mult(bd, bd5); 1472 if (bd2 > 0) 1473 bd = lshift(bd, bd2); 1474 if (bs2 > 0) 1475 bs = lshift(bs, bs2); 1476 delta = diff(bb, bd); 1477 dsign = delta->sign; 1478 delta->sign = 0; 1479 i = cmp(delta, bs); 1480 if (i < 0) { 1481 /* Error is less than half an ulp -- check for 1482 * special case of mantissa a power of two. 1483 */ 1484 if (dsign || word1(rv) || word0(rv) & Bndry_mask) 1485 break; 1486 delta = lshift(delta,Log2P); 1487 if (cmp(delta, bs) > 0) 1488 goto drop_down; 1489 break; 1490 } 1491 if (i == 0) { 1492 /* exactly half-way between */ 1493 if (dsign) { 1494 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 1495 && word1(rv) == 0xffffffff) { 1496 /*boundary case -- increment exponent*/ 1497 word0(rv) = (word0(rv) & Exp_mask) 1498 + Exp_msk1 1499 #ifdef IBM 1500 | Exp_msk1 >> 4 1501 #endif 1502 ; 1503 word1(rv) = 0; 1504 break; 1505 } 1506 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 1507 drop_down: 1508 /* boundary case -- decrement exponent */ 1509 #ifdef Sudden_Underflow 1510 L = word0(rv) & Exp_mask; 1511 #ifdef IBM 1512 if (L < Exp_msk1) 1513 #else 1514 if (L <= Exp_msk1) 1515 #endif 1516 goto undfl; 1517 L -= Exp_msk1; 1518 #else 1519 L = (word0(rv) & Exp_mask) - Exp_msk1; 1520 #endif 1521 word0(rv) = L | Bndry_mask1; 1522 word1(rv) = 0xffffffff; 1523 #ifdef IBM 1524 goto cont; 1525 #else 1526 break; 1527 #endif 1528 } 1529 #ifndef ROUND_BIASED 1530 if (!(word1(rv) & LSB)) 1531 break; 1532 #endif 1533 if (dsign) 1534 rv += ulp(rv); 1535 #ifndef ROUND_BIASED 1536 else { 1537 rv -= ulp(rv); 1538 #ifndef Sudden_Underflow 1539 if (!rv) 1540 goto undfl; 1541 #endif 1542 } 1543 #endif 1544 break; 1545 } 1546 if ((aadj = ratio(delta, bs)) <= 2.) { 1547 if (dsign) 1548 aadj = aadj1 = 1.; 1549 else if (word1(rv) || word0(rv) & Bndry_mask) { 1550 #ifndef Sudden_Underflow 1551 if (word1(rv) == Tiny1 && !word0(rv)) 1552 goto undfl; 1553 #endif 1554 aadj = 1.; 1555 aadj1 = -1.; 1556 } else { 1557 /* special case -- power of FLT_RADIX to be */ 1558 /* rounded down... */ 1559 1560 if (aadj < 2./FLT_RADIX) 1561 aadj = 1./FLT_RADIX; 1562 else 1563 aadj *= 0.5; 1564 aadj1 = -aadj; 1565 } 1566 } else { 1567 aadj *= 0.5; 1568 aadj1 = dsign ? aadj : -aadj; 1569 #ifdef Check_FLT_ROUNDS 1570 switch(FLT_ROUNDS) { 1571 case 2: /* towards +infinity */ 1572 aadj1 -= 0.5; 1573 break; 1574 case 0: /* towards 0 */ 1575 case 3: /* towards -infinity */ 1576 aadj1 += 0.5; 1577 } 1578 #else 1579 if (FLT_ROUNDS == 0) 1580 aadj1 += 0.5; 1581 #endif 1582 } 1583 y = word0(rv) & Exp_mask; 1584 1585 /* Check for overflow */ 1586 1587 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 1588 rv0 = rv; 1589 word0(rv) -= P*Exp_msk1; 1590 adj = aadj1 * ulp(rv); 1591 rv += adj; 1592 if ((word0(rv) & Exp_mask) >= 1593 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 1594 if (word0(rv0) == Big0 && word1(rv0) == Big1) 1595 goto ovfl; 1596 word0(rv) = Big0; 1597 word1(rv) = Big1; 1598 goto cont; 1599 } else 1600 word0(rv) += P*Exp_msk1; 1601 } else { 1602 #ifdef Sudden_Underflow 1603 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 1604 rv0 = rv; 1605 word0(rv) += P*Exp_msk1; 1606 adj = aadj1 * ulp(rv); 1607 rv += adj; 1608 #ifdef IBM 1609 if ((word0(rv) & Exp_mask) < P*Exp_msk1) 1610 #else 1611 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 1612 #endif 1613 { 1614 if (word0(rv0) == Tiny0 1615 && word1(rv0) == Tiny1) 1616 goto undfl; 1617 word0(rv) = Tiny0; 1618 word1(rv) = Tiny1; 1619 goto cont; 1620 } else 1621 word0(rv) -= P*Exp_msk1; 1622 } else { 1623 adj = aadj1 * ulp(rv); 1624 rv += adj; 1625 } 1626 #else 1627 /* Compute adj so that the IEEE rounding rules will 1628 * correctly round rv + adj in some half-way cases. 1629 * If rv * ulp(rv) is denormalized (i.e., 1630 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 1631 * trouble from bits lost to denormalization; 1632 * example: 1.2e-307 . 1633 */ 1634 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) { 1635 aadj1 = (double)(int)(aadj + 0.5); 1636 if (!dsign) 1637 aadj1 = -aadj1; 1638 } 1639 adj = aadj1 * ulp(rv); 1640 rv += adj; 1641 #endif 1642 } 1643 z = word0(rv) & Exp_mask; 1644 if (y == z) { 1645 /* Can we stop now? */ 1646 L = aadj; 1647 aadj -= L; 1648 /* The tolerances below are conservative. */ 1649 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 1650 if (aadj < .4999999 || aadj > .5000001) 1651 break; 1652 } else if (aadj < .4999999/FLT_RADIX) 1653 break; 1654 } 1655 cont: 1656 Bfree(bb); 1657 Bfree(bd); 1658 Bfree(bs); 1659 Bfree(delta); 1660 } 1661 Bfree(bb); 1662 Bfree(bd); 1663 Bfree(bs); 1664 Bfree(bd0); 1665 Bfree(delta); 1666 ret: 1667 if (se) 1668 *se = (char *)s; 1669 return sign ? -rv : rv; 1670 } 1671 1672 static int 1673 quorem 1674 #ifdef KR_headers 1675 (b, S) Bigint *b, *S; 1676 #else 1677 (Bigint *b, Bigint *S) 1678 #endif 1679 { 1680 int n; 1681 long borrow, y; 1682 unsigned long carry, q, ys; 1683 unsigned long *bx, *bxe, *sx, *sxe; 1684 #ifdef Pack_32 1685 long z; 1686 unsigned long si, zs; 1687 #endif 1688 1689 n = S->wds; 1690 #ifdef DEBUG 1691 /*debug*/ if (b->wds > n) 1692 /*debug*/ Bug("oversize b in quorem"); 1693 #endif 1694 if (b->wds < n) 1695 return 0; 1696 sx = S->x; 1697 sxe = sx + --n; 1698 bx = b->x; 1699 bxe = bx + n; 1700 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 1701 #ifdef DEBUG 1702 /*debug*/ if (q > 9) 1703 /*debug*/ Bug("oversized quotient in quorem"); 1704 #endif 1705 if (q) { 1706 borrow = 0; 1707 carry = 0; 1708 do { 1709 #ifdef Pack_32 1710 si = *sx++; 1711 ys = (si & 0xffff) * q + carry; 1712 zs = (si >> 16) * q + (ys >> 16); 1713 carry = zs >> 16; 1714 y = (*bx & 0xffff) - (ys & 0xffff) + borrow; 1715 borrow = y >> 16; 1716 Sign_Extend(borrow, y); 1717 z = (*bx >> 16) - (zs & 0xffff) + borrow; 1718 borrow = z >> 16; 1719 Sign_Extend(borrow, z); 1720 Storeinc(bx, z, y); 1721 #else 1722 ys = *sx++ * q + carry; 1723 carry = ys >> 16; 1724 y = *bx - (ys & 0xffff) + borrow; 1725 borrow = y >> 16; 1726 Sign_Extend(borrow, y); 1727 *bx++ = y & 0xffff; 1728 #endif 1729 } while (sx <= sxe); 1730 if (!*bxe) { 1731 bx = b->x; 1732 while (--bxe > bx && !*bxe) 1733 --n; 1734 b->wds = n; 1735 } 1736 } 1737 if (cmp(b, S) >= 0) { 1738 q++; 1739 borrow = 0; 1740 carry = 0; 1741 bx = b->x; 1742 sx = S->x; 1743 do { 1744 #ifdef Pack_32 1745 si = *sx++; 1746 ys = (si & 0xffff) + carry; 1747 zs = (si >> 16) + (ys >> 16); 1748 carry = zs >> 16; 1749 y = (*bx & 0xffff) - (ys & 0xffff) + borrow; 1750 borrow = y >> 16; 1751 Sign_Extend(borrow, y); 1752 z = (*bx >> 16) - (zs & 0xffff) + borrow; 1753 borrow = z >> 16; 1754 Sign_Extend(borrow, z); 1755 Storeinc(bx, z, y); 1756 #else 1757 ys = *sx++ + carry; 1758 carry = ys >> 16; 1759 y = *bx - (ys & 0xffff) + borrow; 1760 borrow = y >> 16; 1761 Sign_Extend(borrow, y); 1762 *bx++ = y & 0xffff; 1763 #endif 1764 } while (sx <= sxe); 1765 bx = b->x; 1766 bxe = bx + n; 1767 if (!*bxe) { 1768 while (--bxe > bx && !*bxe) 1769 --n; 1770 b->wds = n; 1771 } 1772 } 1773 return q; 1774 } 1775 1776 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 1777 * 1778 * Inspired by "How to Print Floating-Point Numbers Accurately" by 1779 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101]. 1780 * 1781 * Modifications: 1782 * 1. Rather than iterating, we use a simple numeric overestimate 1783 * to determine k = floor(log10(d)). We scale relevant 1784 * quantities using O(log2(k)) rather than O(k) multiplications. 1785 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 1786 * try to generate digits strictly left to right. Instead, we 1787 * compute with fewer bits and propagate the carry if necessary 1788 * when rounding the final digit up. This is often faster. 1789 * 3. Under the assumption that input will be rounded nearest, 1790 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 1791 * That is, we allow equality in stopping tests when the 1792 * round-nearest rule will give the same floating-point value 1793 * as would satisfaction of the stopping test with strict 1794 * inequality. 1795 * 4. We remove common factors of powers of 2 from relevant 1796 * quantities. 1797 * 5. When converting floating-point integers less than 1e16, 1798 * we use floating-point arithmetic rather than resorting 1799 * to multiple-precision integers. 1800 * 6. When asked to produce fewer than 15 digits, we first try 1801 * to get by with floating-point arithmetic; we resort to 1802 * multiple-precision integer arithmetic only if we cannot 1803 * guarantee that the floating-point calculation has given 1804 * the correctly rounded result. For k requested digits and 1805 * "uniformly" distributed input, the probability is 1806 * something like 10^(k-15) that we must resort to the long 1807 * calculation. 1808 */ 1809 1810 char * 1811 __dtoa 1812 #ifdef KR_headers 1813 (d, mode, ndigits, decpt, sign, rve) 1814 double d; int mode, ndigits, *decpt, *sign; char **rve; 1815 #else 1816 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve) 1817 #endif 1818 { 1819 /* Arguments ndigits, decpt, sign are similar to those 1820 of ecvt and fcvt; trailing zeros are suppressed from 1821 the returned string. If not null, *rve is set to point 1822 to the end of the return value. If d is +-Infinity or NaN, 1823 then *decpt is set to 9999. 1824 1825 mode: 1826 0 ==> shortest string that yields d when read in 1827 and rounded to nearest. 1828 1 ==> like 0, but with Steele & White stopping rule; 1829 e.g. with IEEE P754 arithmetic , mode 0 gives 1830 1e23 whereas mode 1 gives 9.999999999999999e22. 1831 2 ==> max(1,ndigits) significant digits. This gives a 1832 return value similar to that of ecvt, except 1833 that trailing zeros are suppressed. 1834 3 ==> through ndigits past the decimal point. This 1835 gives a return value similar to that from fcvt, 1836 except that trailing zeros are suppressed, and 1837 ndigits can be negative. 1838 4-9 should give the same return values as 2-3, i.e., 1839 4 <= mode <= 9 ==> same return as mode 1840 2 + (mode & 1). These modes are mainly for 1841 debugging; often they run slower but sometimes 1842 faster than modes 2-3. 1843 4,5,8,9 ==> left-to-right digit generation. 1844 6-9 ==> don't try fast floating-point estimate 1845 (if applicable). 1846 1847 Values of mode other than 0-9 are treated as mode 0. 1848 1849 Sufficient space is allocated to the return value 1850 to hold the suppressed trailing zeros. 1851 */ 1852 1853 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, 1854 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, 1855 spec_case, try_quick; 1856 long L; 1857 #ifndef Sudden_Underflow 1858 int denorm; 1859 unsigned long x; 1860 #endif 1861 Bigint *b, *b1, *delta, *mlo, *mhi, *S; 1862 double d2, ds, eps; 1863 char *s, *s0; 1864 static Bigint *result; 1865 static int result_k; 1866 1867 if (result) { 1868 result->k = result_k; 1869 result->maxwds = 1 << result_k; 1870 Bfree(result); 1871 result = 0; 1872 } 1873 1874 if (word0(d) & Sign_bit) { 1875 /* set sign for everything, including 0's and NaNs */ 1876 *sign = 1; 1877 word0(d) &= ~Sign_bit; /* clear sign bit */ 1878 } 1879 else 1880 *sign = 0; 1881 1882 #if defined(IEEE_Arith) + defined(VAX) 1883 #ifdef IEEE_Arith 1884 if ((word0(d) & Exp_mask) == Exp_mask) 1885 #else 1886 if (word0(d) == 0x8000) 1887 #endif 1888 { 1889 /* Infinity or NaN */ 1890 *decpt = 9999; 1891 s = 1892 #ifdef IEEE_Arith 1893 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" : 1894 #endif 1895 "NaN"; 1896 if (rve) 1897 *rve = 1898 #ifdef IEEE_Arith 1899 s[3] ? s + 8 : 1900 #endif 1901 s + 3; 1902 return s; 1903 } 1904 #endif 1905 #ifdef IBM 1906 d += 0; /* normalize */ 1907 #endif 1908 if (!d) { 1909 *decpt = 1; 1910 s = "0"; 1911 if (rve) 1912 *rve = s + 1; 1913 return s; 1914 } 1915 1916 b = d2b(d, &be, &bbits); 1917 #ifdef Sudden_Underflow 1918 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); 1919 #else 1920 if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) { 1921 #endif 1922 d2 = d; 1923 word0(d2) &= Frac_mask1; 1924 word0(d2) |= Exp_11; 1925 #ifdef IBM 1926 if (j = 11 - hi0bits(word0(d2) & Frac_mask)) 1927 d2 /= 1 << j; 1928 #endif 1929 1930 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 1931 * log10(x) = log(x) / log(10) 1932 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 1933 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 1934 * 1935 * This suggests computing an approximation k to log10(d) by 1936 * 1937 * k = (i - Bias)*0.301029995663981 1938 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 1939 * 1940 * We want k to be too large rather than too small. 1941 * The error in the first-order Taylor series approximation 1942 * is in our favor, so we just round up the constant enough 1943 * to compensate for any error in the multiplication of 1944 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 1945 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 1946 * adding 1e-13 to the constant term more than suffices. 1947 * Hence we adjust the constant term to 0.1760912590558. 1948 * (We could get a more accurate k by invoking log10, 1949 * but this is probably not worthwhile.) 1950 */ 1951 1952 i -= Bias; 1953 #ifdef IBM 1954 i <<= 2; 1955 i += j; 1956 #endif 1957 #ifndef Sudden_Underflow 1958 denorm = 0; 1959 } else { 1960 /* d is denormalized */ 1961 1962 i = bbits + be + (Bias + (P-1) - 1); 1963 x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32 1964 : word1(d) << 32 - i; 1965 d2 = x; 1966 word0(d2) -= 31*Exp_msk1; /* adjust exponent */ 1967 i -= (Bias + (P-1) - 1) + 1; 1968 denorm = 1; 1969 } 1970 #endif 1971 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; 1972 k = (int)ds; 1973 if (ds < 0. && ds != k) 1974 k--; /* want k = floor(ds) */ 1975 k_check = 1; 1976 if (k >= 0 && k <= Ten_pmax) { 1977 if (d < tens[k]) 1978 k--; 1979 k_check = 0; 1980 } 1981 j = bbits - i - 1; 1982 if (j >= 0) { 1983 b2 = 0; 1984 s2 = j; 1985 } else { 1986 b2 = -j; 1987 s2 = 0; 1988 } 1989 if (k >= 0) { 1990 b5 = 0; 1991 s5 = k; 1992 s2 += k; 1993 } else { 1994 b2 -= k; 1995 b5 = -k; 1996 s5 = 0; 1997 } 1998 if (mode < 0 || mode > 9) 1999 mode = 0; 2000 try_quick = 1; 2001 if (mode > 5) { 2002 mode -= 4; 2003 try_quick = 0; 2004 } 2005 leftright = 1; 2006 switch(mode) { 2007 case 0: 2008 case 1: 2009 ilim = ilim1 = -1; 2010 i = 18; 2011 ndigits = 0; 2012 break; 2013 case 2: 2014 leftright = 0; 2015 /* no break */ 2016 case 4: 2017 if (ndigits <= 0) 2018 ndigits = 1; 2019 ilim = ilim1 = i = ndigits; 2020 break; 2021 case 3: 2022 leftright = 0; 2023 /* no break */ 2024 case 5: 2025 i = ndigits + k + 1; 2026 ilim = i; 2027 ilim1 = i - 1; 2028 if (i <= 0) 2029 i = 1; 2030 } 2031 j = sizeof(unsigned long); 2032 for (result_k = 0; sizeof(Bigint) - sizeof(unsigned long) + j < i; 2033 j <<= 1) result_k++; 2034 result = Balloc(result_k); 2035 s = s0 = (char *)result; 2036 2037 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 2038 2039 /* Try to get by with floating-point arithmetic. */ 2040 2041 i = 0; 2042 d2 = d; 2043 k0 = k; 2044 ilim0 = ilim; 2045 ieps = 2; /* conservative */ 2046 if (k > 0) { 2047 ds = tens[k&0xf]; 2048 j = k >> 4; 2049 if (j & Bletch) { 2050 /* prevent overflows */ 2051 j &= Bletch - 1; 2052 d /= bigtens[n_bigtens-1]; 2053 ieps++; 2054 } 2055 for (; j; j >>= 1, i++) 2056 if (j & 1) { 2057 ieps++; 2058 ds *= bigtens[i]; 2059 } 2060 d /= ds; 2061 } else if (j1 = -k) { 2062 d *= tens[j1 & 0xf]; 2063 for (j = j1 >> 4; j; j >>= 1, i++) 2064 if (j & 1) { 2065 ieps++; 2066 d *= bigtens[i]; 2067 } 2068 } 2069 if (k_check && d < 1. && ilim > 0) { 2070 if (ilim1 <= 0) 2071 goto fast_failed; 2072 ilim = ilim1; 2073 k--; 2074 d *= 10.; 2075 ieps++; 2076 } 2077 eps = ieps*d + 7.; 2078 word0(eps) -= (P-1)*Exp_msk1; 2079 if (ilim == 0) { 2080 S = mhi = 0; 2081 d -= 5.; 2082 if (d > eps) 2083 goto one_digit; 2084 if (d < -eps) 2085 goto no_digits; 2086 goto fast_failed; 2087 } 2088 #ifndef No_leftright 2089 if (leftright) { 2090 /* Use Steele & White method of only 2091 * generating digits needed. 2092 */ 2093 eps = 0.5/tens[ilim-1] - eps; 2094 for (i = 0;;) { 2095 L = d; 2096 d -= L; 2097 *s++ = '0' + (int)L; 2098 if (d < eps) 2099 goto ret1; 2100 if (1. - d < eps) 2101 goto bump_up; 2102 if (++i >= ilim) 2103 break; 2104 eps *= 10.; 2105 d *= 10.; 2106 } 2107 } else { 2108 #endif 2109 /* Generate ilim digits, then fix them up. */ 2110 eps *= tens[ilim-1]; 2111 for (i = 1;; i++, d *= 10.) { 2112 L = d; 2113 d -= L; 2114 *s++ = '0' + (int)L; 2115 if (i == ilim) { 2116 if (d > 0.5 + eps) 2117 goto bump_up; 2118 else if (d < 0.5 - eps) { 2119 while (*--s == '0'); 2120 s++; 2121 goto ret1; 2122 } 2123 break; 2124 } 2125 } 2126 #ifndef No_leftright 2127 } 2128 #endif 2129 fast_failed: 2130 s = s0; 2131 d = d2; 2132 k = k0; 2133 ilim = ilim0; 2134 } 2135 2136 /* Do we have a "small" integer? */ 2137 2138 if (be >= 0 && k <= Int_max) { 2139 /* Yes. */ 2140 ds = tens[k]; 2141 if (ndigits < 0 && ilim <= 0) { 2142 S = mhi = 0; 2143 if (ilim < 0 || d <= 5*ds) 2144 goto no_digits; 2145 goto one_digit; 2146 } 2147 for (i = 1;; i++) { 2148 L = d / ds; 2149 d -= L*ds; 2150 #ifdef Check_FLT_ROUNDS 2151 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 2152 if (d < 0) { 2153 L--; 2154 d += ds; 2155 } 2156 #endif 2157 *s++ = '0' + (int)L; 2158 if (i == ilim) { 2159 d += d; 2160 if (d > ds || d == ds && L & 1) { 2161 bump_up: 2162 while (*--s == '9') 2163 if (s == s0) { 2164 k++; 2165 *s = '0'; 2166 break; 2167 } 2168 ++*s++; 2169 } 2170 break; 2171 } 2172 if (!(d *= 10.)) 2173 break; 2174 } 2175 goto ret1; 2176 } 2177 2178 m2 = b2; 2179 m5 = b5; 2180 mhi = mlo = 0; 2181 if (leftright) { 2182 if (mode < 2) { 2183 i = 2184 #ifndef Sudden_Underflow 2185 denorm ? be + (Bias + (P-1) - 1 + 1) : 2186 #endif 2187 #ifdef IBM 2188 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); 2189 #else 2190 1 + P - bbits; 2191 #endif 2192 } else { 2193 j = ilim - 1; 2194 if (m5 >= j) 2195 m5 -= j; 2196 else { 2197 s5 += j -= m5; 2198 b5 += j; 2199 m5 = 0; 2200 } 2201 if ((i = ilim) < 0) { 2202 m2 -= i; 2203 i = 0; 2204 } 2205 } 2206 b2 += i; 2207 s2 += i; 2208 mhi = i2b(1); 2209 } 2210 if (m2 > 0 && s2 > 0) { 2211 i = m2 < s2 ? m2 : s2; 2212 b2 -= i; 2213 m2 -= i; 2214 s2 -= i; 2215 } 2216 if (b5 > 0) { 2217 if (leftright) { 2218 if (m5 > 0) { 2219 mhi = pow5mult(mhi, m5); 2220 b1 = mult(mhi, b); 2221 Bfree(b); 2222 b = b1; 2223 } 2224 if (j = b5 - m5) 2225 b = pow5mult(b, j); 2226 } else 2227 b = pow5mult(b, b5); 2228 } 2229 S = i2b(1); 2230 if (s5 > 0) 2231 S = pow5mult(S, s5); 2232 2233 /* Check for special case that d is a normalized power of 2. */ 2234 2235 if (mode < 2) { 2236 if (!word1(d) && !(word0(d) & Bndry_mask) 2237 #ifndef Sudden_Underflow 2238 && word0(d) & Exp_mask 2239 #endif 2240 ) { 2241 /* The special case */ 2242 b2 += Log2P; 2243 s2 += Log2P; 2244 spec_case = 1; 2245 } else 2246 spec_case = 0; 2247 } 2248 2249 /* Arrange for convenient computation of quotients: 2250 * shift left if necessary so divisor has 4 leading 0 bits. 2251 * 2252 * Perhaps we should just compute leading 28 bits of S once 2253 * and for all and pass them and a shift to quorem, so it 2254 * can do shifts and ors to compute the numerator for q. 2255 */ 2256 #ifdef Pack_32 2257 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) 2258 i = 32 - i; 2259 #else 2260 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) 2261 i = 16 - i; 2262 #endif 2263 if (i > 4) { 2264 i -= 4; 2265 b2 += i; 2266 m2 += i; 2267 s2 += i; 2268 } else if (i < 4) { 2269 i += 28; 2270 b2 += i; 2271 m2 += i; 2272 s2 += i; 2273 } 2274 if (b2 > 0) 2275 b = lshift(b, b2); 2276 if (s2 > 0) 2277 S = lshift(S, s2); 2278 if (k_check) { 2279 if (cmp(b,S) < 0) { 2280 k--; 2281 b = multadd(b, 10, 0); /* we botched the k estimate */ 2282 if (leftright) 2283 mhi = multadd(mhi, 10, 0); 2284 ilim = ilim1; 2285 } 2286 } 2287 if (ilim <= 0 && mode > 2) { 2288 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 2289 /* no digits, fcvt style */ 2290 no_digits: 2291 k = -1 - ndigits; 2292 goto ret; 2293 } 2294 one_digit: 2295 *s++ = '1'; 2296 k++; 2297 goto ret; 2298 } 2299 if (leftright) { 2300 if (m2 > 0) 2301 mhi = lshift(mhi, m2); 2302 2303 /* Compute mlo -- check for special case 2304 * that d is a normalized power of 2. 2305 */ 2306 2307 mlo = mhi; 2308 if (spec_case) { 2309 mhi = Balloc(mhi->k); 2310 Bcopy(mhi, mlo); 2311 mhi = lshift(mhi, Log2P); 2312 } 2313 2314 for (i = 1;;i++) { 2315 dig = quorem(b,S) + '0'; 2316 /* Do we yet have the shortest decimal string 2317 * that will round to d? 2318 */ 2319 j = cmp(b, mlo); 2320 delta = diff(S, mhi); 2321 j1 = delta->sign ? 1 : cmp(b, delta); 2322 Bfree(delta); 2323 #ifndef ROUND_BIASED 2324 if (j1 == 0 && !mode && !(word1(d) & 1)) { 2325 if (dig == '9') 2326 goto round_9_up; 2327 if (j > 0) 2328 dig++; 2329 *s++ = dig; 2330 goto ret; 2331 } 2332 #endif 2333 if (j < 0 || j == 0 && !mode 2334 #ifndef ROUND_BIASED 2335 && !(word1(d) & 1) 2336 #endif 2337 ) { 2338 if (j1 > 0) { 2339 b = lshift(b, 1); 2340 j1 = cmp(b, S); 2341 if ((j1 > 0 || j1 == 0 && dig & 1) 2342 && dig++ == '9') 2343 goto round_9_up; 2344 } 2345 *s++ = dig; 2346 goto ret; 2347 } 2348 if (j1 > 0) { 2349 if (dig == '9') { /* possible if i == 1 */ 2350 round_9_up: 2351 *s++ = '9'; 2352 goto roundoff; 2353 } 2354 *s++ = dig + 1; 2355 goto ret; 2356 } 2357 *s++ = dig; 2358 if (i == ilim) 2359 break; 2360 b = multadd(b, 10, 0); 2361 if (mlo == mhi) 2362 mlo = mhi = multadd(mhi, 10, 0); 2363 else { 2364 mlo = multadd(mlo, 10, 0); 2365 mhi = multadd(mhi, 10, 0); 2366 } 2367 } 2368 } else 2369 for (i = 1;; i++) { 2370 *s++ = dig = quorem(b,S) + '0'; 2371 if (i >= ilim) 2372 break; 2373 b = multadd(b, 10, 0); 2374 } 2375 2376 /* Round off last digit */ 2377 2378 b = lshift(b, 1); 2379 j = cmp(b, S); 2380 if (j > 0 || j == 0 && dig & 1) { 2381 roundoff: 2382 while (*--s == '9') 2383 if (s == s0) { 2384 k++; 2385 *s++ = '1'; 2386 goto ret; 2387 } 2388 ++*s++; 2389 } else { 2390 while (*--s == '0'); 2391 s++; 2392 } 2393 ret: 2394 Bfree(S); 2395 if (mhi) { 2396 if (mlo && mlo != mhi) 2397 Bfree(mlo); 2398 Bfree(mhi); 2399 } 2400 ret1: 2401 Bfree(b); 2402 if (s == s0) { /* don't return empty string */ 2403 *s++ = '0'; 2404 k = 0; 2405 } 2406 *s = 0; 2407 *decpt = k + 1; 2408 if (rve) 2409 *rve = s; 2410 return s0; 2411 } 2412 #ifdef __cplusplus 2413 } 2414 #endif 2415