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