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