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