1 /**************************************************************** 2 * 3 * The author of this software is David M. Gay. 4 * 5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. 6 * Copyright (C) 2002, 2005, 2006, 2007, 2008, 2010 Apple Inc. All rights reserved. 7 * 8 * Permission to use, copy, modify, and distribute this software for any 9 * purpose without fee is hereby granted, provided that this entire notice 10 * is included in all copies of any software which is or includes a copy 11 * or modification of this software and in all copies of the supporting 12 * documentation for such software. 13 * 14 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 15 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY 16 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 17 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 18 * 19 ***************************************************************/ 20 21 /* Please send bug reports to David M. Gay (dmg at acm dot org, 22 * with " at " changed at "@" and " dot " changed to "."). */ 23 24 /* On a machine with IEEE extended-precision registers, it is 25 * necessary to specify double-precision (53-bit) rounding precision 26 * before invoking strtod or dtoa. If the machine uses (the equivalent 27 * of) Intel 80x87 arithmetic, the call 28 * _control87(PC_53, MCW_PC); 29 * does this with many compilers. Whether this or another call is 30 * appropriate depends on the compiler; for this to work, it may be 31 * necessary to #include "float.h" or another system-dependent header 32 * file. 33 */ 34 35 /* strtod for IEEE-arithmetic machines. 36 * 37 * This strtod returns a nearest machine number to the input decimal 38 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 39 * broken by the IEEE round-even rule. Otherwise ties are broken by 40 * biased rounding (add half and chop). 41 * 42 * Inspired loosely by William D. Clinger's paper "How to Read Floating 43 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 44 * 45 * Modifications: 46 * 47 * 1. We only require IEEE double-precision arithmetic (not IEEE double-extended). 48 * 2. We get by with floating-point arithmetic in a case that 49 * Clinger missed -- when we're computing d * 10^n 50 * for a small integer d and the integer n is not too 51 * much larger than 22 (the maximum integer k for which 52 * we can represent 10^k exactly), we may be able to 53 * compute (d*10^k) * 10^(e-k) with just one roundoff. 54 * 3. Rather than a bit-at-a-time adjustment of the binary 55 * result in the hard case, we use floating-point 56 * arithmetic to determine the adjustment to within 57 * one bit; only in really hard cases do we need to 58 * compute a second residual. 59 * 4. Because of 3., we don't need a large table of powers of 10 60 * for ten-to-e (just some small tables, e.g. of 10^k 61 * for 0 <= k <= 22). 62 */ 63 64 #include "config.h" 65 #include "dtoa.h" 66 67 #if HAVE(ERRNO_H) 68 #include <errno.h> 69 #endif 70 #include <float.h> 71 #include <math.h> 72 #include <stdint.h> 73 #include <stdio.h> 74 #include <stdlib.h> 75 #include <string.h> 76 #include <wtf/AlwaysInline.h> 77 #include <wtf/Assertions.h> 78 #include <wtf/DecimalNumber.h> 79 #include <wtf/FastMalloc.h> 80 #include <wtf/MathExtras.h> 81 #include <wtf/Threading.h> 82 #include <wtf/UnusedParam.h> 83 #include <wtf/Vector.h> 84 85 #if COMPILER(MSVC) 86 #pragma warning(disable: 4244) 87 #pragma warning(disable: 4245) 88 #pragma warning(disable: 4554) 89 #endif 90 91 namespace WTF { 92 93 #if ENABLE(JSC_MULTIPLE_THREADS) 94 Mutex* s_dtoaP5Mutex; 95 #endif 96 97 typedef union { 98 double d; 99 uint32_t L[2]; 100 } U; 101 102 #if CPU(BIG_ENDIAN) || CPU(MIDDLE_ENDIAN) 103 #define word0(x) (x)->L[0] 104 #define word1(x) (x)->L[1] 105 #else 106 #define word0(x) (x)->L[1] 107 #define word1(x) (x)->L[0] 108 #endif 109 #define dval(x) (x)->d 110 111 /* The following definition of Storeinc is appropriate for MIPS processors. 112 * An alternative that might be better on some machines is 113 * *p++ = high << 16 | low & 0xffff; 114 */ 115 static ALWAYS_INLINE uint32_t* storeInc(uint32_t* p, uint16_t high, uint16_t low) 116 { 117 uint16_t* p16 = reinterpret_cast<uint16_t*>(p); 118 #if CPU(BIG_ENDIAN) 119 p16[0] = high; 120 p16[1] = low; 121 #else 122 p16[1] = high; 123 p16[0] = low; 124 #endif 125 return p + 1; 126 } 127 128 #define Exp_shift 20 129 #define Exp_shift1 20 130 #define Exp_msk1 0x100000 131 #define Exp_msk11 0x100000 132 #define Exp_mask 0x7ff00000 133 #define P 53 134 #define Bias 1023 135 #define Emin (-1022) 136 #define Exp_1 0x3ff00000 137 #define Exp_11 0x3ff00000 138 #define Ebits 11 139 #define Frac_mask 0xfffff 140 #define Frac_mask1 0xfffff 141 #define Ten_pmax 22 142 #define Bletch 0x10 143 #define Bndry_mask 0xfffff 144 #define Bndry_mask1 0xfffff 145 #define LSB 1 146 #define Sign_bit 0x80000000 147 #define Log2P 1 148 #define Tiny0 0 149 #define Tiny1 1 150 #define Quick_max 14 151 #define Int_max 14 152 153 #define rounded_product(a, b) a *= b 154 #define rounded_quotient(a, b) a /= b 155 156 #define Big0 (Frac_mask1 | Exp_msk1 * (DBL_MAX_EXP + Bias - 1)) 157 #define Big1 0xffffffff 158 159 #if CPU(PPC64) || CPU(X86_64) 160 // FIXME: should we enable this on all 64-bit CPUs? 161 // 64-bit emulation provided by the compiler is likely to be slower than dtoa own code on 32-bit hardware. 162 #define USE_LONG_LONG 163 #endif 164 165 struct BigInt { 166 BigInt() : sign(0) { } 167 int sign; 168 169 void clear() 170 { 171 sign = 0; 172 m_words.clear(); 173 } 174 175 size_t size() const 176 { 177 return m_words.size(); 178 } 179 180 void resize(size_t s) 181 { 182 m_words.resize(s); 183 } 184 185 uint32_t* words() 186 { 187 return m_words.data(); 188 } 189 190 const uint32_t* words() const 191 { 192 return m_words.data(); 193 } 194 195 void append(uint32_t w) 196 { 197 m_words.append(w); 198 } 199 200 Vector<uint32_t, 16> m_words; 201 }; 202 203 static void multadd(BigInt& b, int m, int a) /* multiply by m and add a */ 204 { 205 #ifdef USE_LONG_LONG 206 unsigned long long carry; 207 #else 208 uint32_t carry; 209 #endif 210 211 int wds = b.size(); 212 uint32_t* x = b.words(); 213 int i = 0; 214 carry = a; 215 do { 216 #ifdef USE_LONG_LONG 217 unsigned long long y = *x * (unsigned long long)m + carry; 218 carry = y >> 32; 219 *x++ = (uint32_t)y & 0xffffffffUL; 220 #else 221 uint32_t xi = *x; 222 uint32_t y = (xi & 0xffff) * m + carry; 223 uint32_t z = (xi >> 16) * m + (y >> 16); 224 carry = z >> 16; 225 *x++ = (z << 16) + (y & 0xffff); 226 #endif 227 } while (++i < wds); 228 229 if (carry) 230 b.append((uint32_t)carry); 231 } 232 233 static void s2b(BigInt& b, const char* s, int nd0, int nd, uint32_t y9) 234 { 235 b.sign = 0; 236 b.resize(1); 237 b.words()[0] = y9; 238 239 int i = 9; 240 if (9 < nd0) { 241 s += 9; 242 do { 243 multadd(b, 10, *s++ - '0'); 244 } while (++i < nd0); 245 s++; 246 } else 247 s += 10; 248 for (; i < nd; i++) 249 multadd(b, 10, *s++ - '0'); 250 } 251 252 static int hi0bits(uint32_t x) 253 { 254 int k = 0; 255 256 if (!(x & 0xffff0000)) { 257 k = 16; 258 x <<= 16; 259 } 260 if (!(x & 0xff000000)) { 261 k += 8; 262 x <<= 8; 263 } 264 if (!(x & 0xf0000000)) { 265 k += 4; 266 x <<= 4; 267 } 268 if (!(x & 0xc0000000)) { 269 k += 2; 270 x <<= 2; 271 } 272 if (!(x & 0x80000000)) { 273 k++; 274 if (!(x & 0x40000000)) 275 return 32; 276 } 277 return k; 278 } 279 280 static int lo0bits(uint32_t* y) 281 { 282 int k; 283 uint32_t x = *y; 284 285 if (x & 7) { 286 if (x & 1) 287 return 0; 288 if (x & 2) { 289 *y = x >> 1; 290 return 1; 291 } 292 *y = x >> 2; 293 return 2; 294 } 295 k = 0; 296 if (!(x & 0xffff)) { 297 k = 16; 298 x >>= 16; 299 } 300 if (!(x & 0xff)) { 301 k += 8; 302 x >>= 8; 303 } 304 if (!(x & 0xf)) { 305 k += 4; 306 x >>= 4; 307 } 308 if (!(x & 0x3)) { 309 k += 2; 310 x >>= 2; 311 } 312 if (!(x & 1)) { 313 k++; 314 x >>= 1; 315 if (!x) 316 return 32; 317 } 318 *y = x; 319 return k; 320 } 321 322 static void i2b(BigInt& b, int i) 323 { 324 b.sign = 0; 325 b.resize(1); 326 b.words()[0] = i; 327 } 328 329 static void mult(BigInt& aRef, const BigInt& bRef) 330 { 331 const BigInt* a = &aRef; 332 const BigInt* b = &bRef; 333 BigInt c; 334 int wa, wb, wc; 335 const uint32_t* x = 0; 336 const uint32_t* xa; 337 const uint32_t* xb; 338 const uint32_t* xae; 339 const uint32_t* xbe; 340 uint32_t* xc; 341 uint32_t* xc0; 342 uint32_t y; 343 #ifdef USE_LONG_LONG 344 unsigned long long carry, z; 345 #else 346 uint32_t carry, z; 347 #endif 348 349 if (a->size() < b->size()) { 350 const BigInt* tmp = a; 351 a = b; 352 b = tmp; 353 } 354 355 wa = a->size(); 356 wb = b->size(); 357 wc = wa + wb; 358 c.resize(wc); 359 360 for (xc = c.words(), xa = xc + wc; xc < xa; xc++) 361 *xc = 0; 362 xa = a->words(); 363 xae = xa + wa; 364 xb = b->words(); 365 xbe = xb + wb; 366 xc0 = c.words(); 367 #ifdef USE_LONG_LONG 368 for (; xb < xbe; xc0++) { 369 if ((y = *xb++)) { 370 x = xa; 371 xc = xc0; 372 carry = 0; 373 do { 374 z = *x++ * (unsigned long long)y + *xc + carry; 375 carry = z >> 32; 376 *xc++ = (uint32_t)z & 0xffffffffUL; 377 } while (x < xae); 378 *xc = (uint32_t)carry; 379 } 380 } 381 #else 382 for (; xb < xbe; xb++, xc0++) { 383 if ((y = *xb & 0xffff)) { 384 x = xa; 385 xc = xc0; 386 carry = 0; 387 do { 388 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 389 carry = z >> 16; 390 uint32_t z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 391 carry = z2 >> 16; 392 xc = storeInc(xc, z2, z); 393 } while (x < xae); 394 *xc = carry; 395 } 396 if ((y = *xb >> 16)) { 397 x = xa; 398 xc = xc0; 399 carry = 0; 400 uint32_t z2 = *xc; 401 do { 402 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 403 carry = z >> 16; 404 xc = storeInc(xc, z, z2); 405 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 406 carry = z2 >> 16; 407 } while (x < xae); 408 *xc = z2; 409 } 410 } 411 #endif 412 for (xc0 = c.words(), xc = xc0 + wc; wc > 0 && !*--xc; --wc) { } 413 c.resize(wc); 414 aRef = c; 415 } 416 417 struct P5Node { 418 WTF_MAKE_NONCOPYABLE(P5Node); WTF_MAKE_FAST_ALLOCATED; 419 public: 420 P5Node() { } 421 BigInt val; 422 P5Node* next; 423 }; 424 425 static P5Node* p5s; 426 static int p5sCount; 427 428 static ALWAYS_INLINE void pow5mult(BigInt& b, int k) 429 { 430 static int p05[3] = { 5, 25, 125 }; 431 432 if (int i = k & 3) 433 multadd(b, p05[i - 1], 0); 434 435 if (!(k >>= 2)) 436 return; 437 438 #if ENABLE(JSC_MULTIPLE_THREADS) 439 s_dtoaP5Mutex->lock(); 440 #endif 441 P5Node* p5 = p5s; 442 443 if (!p5) { 444 /* first time */ 445 p5 = new P5Node; 446 i2b(p5->val, 625); 447 p5->next = 0; 448 p5s = p5; 449 p5sCount = 1; 450 } 451 452 int p5sCountLocal = p5sCount; 453 #if ENABLE(JSC_MULTIPLE_THREADS) 454 s_dtoaP5Mutex->unlock(); 455 #endif 456 int p5sUsed = 0; 457 458 for (;;) { 459 if (k & 1) 460 mult(b, p5->val); 461 462 if (!(k >>= 1)) 463 break; 464 465 if (++p5sUsed == p5sCountLocal) { 466 #if ENABLE(JSC_MULTIPLE_THREADS) 467 s_dtoaP5Mutex->lock(); 468 #endif 469 if (p5sUsed == p5sCount) { 470 ASSERT(!p5->next); 471 p5->next = new P5Node; 472 p5->next->next = 0; 473 p5->next->val = p5->val; 474 mult(p5->next->val, p5->next->val); 475 ++p5sCount; 476 } 477 478 p5sCountLocal = p5sCount; 479 #if ENABLE(JSC_MULTIPLE_THREADS) 480 s_dtoaP5Mutex->unlock(); 481 #endif 482 } 483 p5 = p5->next; 484 } 485 } 486 487 static ALWAYS_INLINE void lshift(BigInt& b, int k) 488 { 489 int n = k >> 5; 490 491 int origSize = b.size(); 492 int n1 = n + origSize + 1; 493 494 if (k &= 0x1f) 495 b.resize(b.size() + n + 1); 496 else 497 b.resize(b.size() + n); 498 499 const uint32_t* srcStart = b.words(); 500 uint32_t* dstStart = b.words(); 501 const uint32_t* src = srcStart + origSize - 1; 502 uint32_t* dst = dstStart + n1 - 1; 503 if (k) { 504 uint32_t hiSubword = 0; 505 int s = 32 - k; 506 for (; src >= srcStart; --src) { 507 *dst-- = hiSubword | *src >> s; 508 hiSubword = *src << k; 509 } 510 *dst = hiSubword; 511 ASSERT(dst == dstStart + n); 512 513 b.resize(origSize + n + !!b.words()[n1 - 1]); 514 } 515 else { 516 do { 517 *--dst = *src--; 518 } while (src >= srcStart); 519 } 520 for (dst = dstStart + n; dst != dstStart; ) 521 *--dst = 0; 522 523 ASSERT(b.size() <= 1 || b.words()[b.size() - 1]); 524 } 525 526 static int cmp(const BigInt& a, const BigInt& b) 527 { 528 const uint32_t *xa, *xa0, *xb, *xb0; 529 int i, j; 530 531 i = a.size(); 532 j = b.size(); 533 ASSERT(i <= 1 || a.words()[i - 1]); 534 ASSERT(j <= 1 || b.words()[j - 1]); 535 if (i -= j) 536 return i; 537 xa0 = a.words(); 538 xa = xa0 + j; 539 xb0 = b.words(); 540 xb = xb0 + j; 541 for (;;) { 542 if (*--xa != *--xb) 543 return *xa < *xb ? -1 : 1; 544 if (xa <= xa0) 545 break; 546 } 547 return 0; 548 } 549 550 static ALWAYS_INLINE void diff(BigInt& c, const BigInt& aRef, const BigInt& bRef) 551 { 552 const BigInt* a = &aRef; 553 const BigInt* b = &bRef; 554 int i, wa, wb; 555 uint32_t* xc; 556 557 i = cmp(*a, *b); 558 if (!i) { 559 c.sign = 0; 560 c.resize(1); 561 c.words()[0] = 0; 562 return; 563 } 564 if (i < 0) { 565 const BigInt* tmp = a; 566 a = b; 567 b = tmp; 568 i = 1; 569 } else 570 i = 0; 571 572 wa = a->size(); 573 const uint32_t* xa = a->words(); 574 const uint32_t* xae = xa + wa; 575 wb = b->size(); 576 const uint32_t* xb = b->words(); 577 const uint32_t* xbe = xb + wb; 578 579 c.resize(wa); 580 c.sign = i; 581 xc = c.words(); 582 #ifdef USE_LONG_LONG 583 unsigned long long borrow = 0; 584 do { 585 unsigned long long y = (unsigned long long)*xa++ - *xb++ - borrow; 586 borrow = y >> 32 & (uint32_t)1; 587 *xc++ = (uint32_t)y & 0xffffffffUL; 588 } while (xb < xbe); 589 while (xa < xae) { 590 unsigned long long y = *xa++ - borrow; 591 borrow = y >> 32 & (uint32_t)1; 592 *xc++ = (uint32_t)y & 0xffffffffUL; 593 } 594 #else 595 uint32_t borrow = 0; 596 do { 597 uint32_t y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 598 borrow = (y & 0x10000) >> 16; 599 uint32_t z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 600 borrow = (z & 0x10000) >> 16; 601 xc = storeInc(xc, z, y); 602 } while (xb < xbe); 603 while (xa < xae) { 604 uint32_t y = (*xa & 0xffff) - borrow; 605 borrow = (y & 0x10000) >> 16; 606 uint32_t z = (*xa++ >> 16) - borrow; 607 borrow = (z & 0x10000) >> 16; 608 xc = storeInc(xc, z, y); 609 } 610 #endif 611 while (!*--xc) 612 wa--; 613 c.resize(wa); 614 } 615 616 static double ulp(U *x) 617 { 618 register int32_t L; 619 U u; 620 621 L = (word0(x) & Exp_mask) - (P - 1) * Exp_msk1; 622 word0(&u) = L; 623 word1(&u) = 0; 624 return dval(&u); 625 } 626 627 static double b2d(const BigInt& a, int* e) 628 { 629 const uint32_t* xa; 630 const uint32_t* xa0; 631 uint32_t w; 632 uint32_t y; 633 uint32_t z; 634 int k; 635 U d; 636 637 #define d0 word0(&d) 638 #define d1 word1(&d) 639 640 xa0 = a.words(); 641 xa = xa0 + a.size(); 642 y = *--xa; 643 ASSERT(y); 644 k = hi0bits(y); 645 *e = 32 - k; 646 if (k < Ebits) { 647 d0 = Exp_1 | (y >> (Ebits - k)); 648 w = xa > xa0 ? *--xa : 0; 649 d1 = (y << (32 - Ebits + k)) | (w >> (Ebits - k)); 650 goto returnD; 651 } 652 z = xa > xa0 ? *--xa : 0; 653 if (k -= Ebits) { 654 d0 = Exp_1 | (y << k) | (z >> (32 - k)); 655 y = xa > xa0 ? *--xa : 0; 656 d1 = (z << k) | (y >> (32 - k)); 657 } else { 658 d0 = Exp_1 | y; 659 d1 = z; 660 } 661 returnD: 662 #undef d0 663 #undef d1 664 return dval(&d); 665 } 666 667 static ALWAYS_INLINE void d2b(BigInt& b, U* d, int* e, int* bits) 668 { 669 int de, k; 670 uint32_t* x; 671 uint32_t y, z; 672 int i; 673 #define d0 word0(d) 674 #define d1 word1(d) 675 676 b.sign = 0; 677 b.resize(1); 678 x = b.words(); 679 680 z = d0 & Frac_mask; 681 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 682 if ((de = (int)(d0 >> Exp_shift))) 683 z |= Exp_msk1; 684 if ((y = d1)) { 685 if ((k = lo0bits(&y))) { 686 x[0] = y | (z << (32 - k)); 687 z >>= k; 688 } else 689 x[0] = y; 690 if (z) { 691 b.resize(2); 692 x[1] = z; 693 } 694 695 i = b.size(); 696 } else { 697 k = lo0bits(&z); 698 x[0] = z; 699 i = 1; 700 b.resize(1); 701 k += 32; 702 } 703 if (de) { 704 *e = de - Bias - (P - 1) + k; 705 *bits = P - k; 706 } else { 707 *e = de - Bias - (P - 1) + 1 + k; 708 *bits = (32 * i) - hi0bits(x[i - 1]); 709 } 710 } 711 #undef d0 712 #undef d1 713 714 static double ratio(const BigInt& a, const BigInt& b) 715 { 716 U da, db; 717 int k, ka, kb; 718 719 dval(&da) = b2d(a, &ka); 720 dval(&db) = b2d(b, &kb); 721 k = ka - kb + 32 * (a.size() - b.size()); 722 if (k > 0) 723 word0(&da) += k * Exp_msk1; 724 else { 725 k = -k; 726 word0(&db) += k * Exp_msk1; 727 } 728 return dval(&da) / dval(&db); 729 } 730 731 static const double tens[] = { 732 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 733 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 734 1e20, 1e21, 1e22 735 }; 736 737 static const double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 738 static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 739 9007199254740992. * 9007199254740992.e-256 740 /* = 2^106 * 1e-256 */ 741 }; 742 743 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 744 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 745 #define Scale_Bit 0x10 746 #define n_bigtens 5 747 748 double strtod(const char* s00, char** se) 749 { 750 int scale; 751 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, 752 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 753 const char *s, *s0, *s1; 754 double aadj, aadj1; 755 U aadj2, adj, rv, rv0; 756 int32_t L; 757 uint32_t y, z; 758 BigInt bb, bb1, bd, bd0, bs, delta; 759 760 sign = nz0 = nz = 0; 761 dval(&rv) = 0; 762 for (s = s00; ; s++) { 763 switch (*s) { 764 case '-': 765 sign = 1; 766 /* no break */ 767 case '+': 768 if (*++s) 769 goto break2; 770 /* no break */ 771 case 0: 772 goto ret0; 773 case '\t': 774 case '\n': 775 case '\v': 776 case '\f': 777 case '\r': 778 case ' ': 779 continue; 780 default: 781 goto break2; 782 } 783 } 784 break2: 785 if (*s == '0') { 786 nz0 = 1; 787 while (*++s == '0') { } 788 if (!*s) 789 goto ret; 790 } 791 s0 = s; 792 y = z = 0; 793 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 794 if (nd < 9) 795 y = (10 * y) + c - '0'; 796 else if (nd < 16) 797 z = (10 * z) + c - '0'; 798 nd0 = nd; 799 if (c == '.') { 800 c = *++s; 801 if (!nd) { 802 for (; c == '0'; c = *++s) 803 nz++; 804 if (c > '0' && c <= '9') { 805 s0 = s; 806 nf += nz; 807 nz = 0; 808 goto haveDig; 809 } 810 goto digDone; 811 } 812 for (; c >= '0' && c <= '9'; c = *++s) { 813 haveDig: 814 nz++; 815 if (c -= '0') { 816 nf += nz; 817 for (i = 1; i < nz; i++) 818 if (nd++ < 9) 819 y *= 10; 820 else if (nd <= DBL_DIG + 1) 821 z *= 10; 822 if (nd++ < 9) 823 y = (10 * y) + c; 824 else if (nd <= DBL_DIG + 1) 825 z = (10 * z) + c; 826 nz = 0; 827 } 828 } 829 } 830 digDone: 831 e = 0; 832 if (c == 'e' || c == 'E') { 833 if (!nd && !nz && !nz0) 834 goto ret0; 835 s00 = s; 836 esign = 0; 837 switch (c = *++s) { 838 case '-': 839 esign = 1; 840 case '+': 841 c = *++s; 842 } 843 if (c >= '0' && c <= '9') { 844 while (c == '0') 845 c = *++s; 846 if (c > '0' && c <= '9') { 847 L = c - '0'; 848 s1 = s; 849 while ((c = *++s) >= '0' && c <= '9') 850 L = (10 * L) + c - '0'; 851 if (s - s1 > 8 || L > 19999) 852 /* Avoid confusion from exponents 853 * so large that e might overflow. 854 */ 855 e = 19999; /* safe for 16 bit ints */ 856 else 857 e = (int)L; 858 if (esign) 859 e = -e; 860 } else 861 e = 0; 862 } else 863 s = s00; 864 } 865 if (!nd) { 866 if (!nz && !nz0) { 867 ret0: 868 s = s00; 869 sign = 0; 870 } 871 goto ret; 872 } 873 e1 = e -= nf; 874 875 /* Now we have nd0 digits, starting at s0, followed by a 876 * decimal point, followed by nd-nd0 digits. The number we're 877 * after is the integer represented by those digits times 878 * 10**e */ 879 880 if (!nd0) 881 nd0 = nd; 882 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 883 dval(&rv) = y; 884 if (k > 9) 885 dval(&rv) = tens[k - 9] * dval(&rv) + z; 886 if (nd <= DBL_DIG) { 887 if (!e) 888 goto ret; 889 if (e > 0) { 890 if (e <= Ten_pmax) { 891 /* rv = */ rounded_product(dval(&rv), tens[e]); 892 goto ret; 893 } 894 i = DBL_DIG - nd; 895 if (e <= Ten_pmax + i) { 896 /* A fancier test would sometimes let us do 897 * this for larger i values. 898 */ 899 e -= i; 900 dval(&rv) *= tens[i]; 901 /* rv = */ rounded_product(dval(&rv), tens[e]); 902 goto ret; 903 } 904 } else if (e >= -Ten_pmax) { 905 /* rv = */ rounded_quotient(dval(&rv), tens[-e]); 906 goto ret; 907 } 908 } 909 e1 += nd - k; 910 911 scale = 0; 912 913 /* Get starting approximation = rv * 10**e1 */ 914 915 if (e1 > 0) { 916 if ((i = e1 & 15)) 917 dval(&rv) *= tens[i]; 918 if (e1 &= ~15) { 919 if (e1 > DBL_MAX_10_EXP) { 920 ovfl: 921 #if HAVE(ERRNO_H) 922 errno = ERANGE; 923 #endif 924 /* Can't trust HUGE_VAL */ 925 word0(&rv) = Exp_mask; 926 word1(&rv) = 0; 927 goto ret; 928 } 929 e1 >>= 4; 930 for (j = 0; e1 > 1; j++, e1 >>= 1) 931 if (e1 & 1) 932 dval(&rv) *= bigtens[j]; 933 /* The last multiplication could overflow. */ 934 word0(&rv) -= P * Exp_msk1; 935 dval(&rv) *= bigtens[j]; 936 if ((z = word0(&rv) & Exp_mask) > Exp_msk1 * (DBL_MAX_EXP + Bias - P)) 937 goto ovfl; 938 if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) { 939 /* set to largest number */ 940 /* (Can't trust DBL_MAX) */ 941 word0(&rv) = Big0; 942 word1(&rv) = Big1; 943 } else 944 word0(&rv) += P * Exp_msk1; 945 } 946 } else if (e1 < 0) { 947 e1 = -e1; 948 if ((i = e1 & 15)) 949 dval(&rv) /= tens[i]; 950 if (e1 >>= 4) { 951 if (e1 >= 1 << n_bigtens) 952 goto undfl; 953 if (e1 & Scale_Bit) 954 scale = 2 * P; 955 for (j = 0; e1 > 0; j++, e1 >>= 1) 956 if (e1 & 1) 957 dval(&rv) *= tinytens[j]; 958 if (scale && (j = (2 * P) + 1 - ((word0(&rv) & Exp_mask) >> Exp_shift)) > 0) { 959 /* scaled rv is denormal; clear j low bits */ 960 if (j >= 32) { 961 word1(&rv) = 0; 962 if (j >= 53) 963 word0(&rv) = (P + 2) * Exp_msk1; 964 else 965 word0(&rv) &= 0xffffffff << (j - 32); 966 } else 967 word1(&rv) &= 0xffffffff << j; 968 } 969 if (!dval(&rv)) { 970 undfl: 971 dval(&rv) = 0.; 972 #if HAVE(ERRNO_H) 973 errno = ERANGE; 974 #endif 975 goto ret; 976 } 977 } 978 } 979 980 /* Now the hard part -- adjusting rv to the correct value.*/ 981 982 /* Put digits into bd: true value = bd * 10^e */ 983 984 s2b(bd0, s0, nd0, nd, y); 985 986 for (;;) { 987 bd = bd0; 988 d2b(bb, &rv, &bbe, &bbbits); /* rv = bb * 2^bbe */ 989 i2b(bs, 1); 990 991 if (e >= 0) { 992 bb2 = bb5 = 0; 993 bd2 = bd5 = e; 994 } else { 995 bb2 = bb5 = -e; 996 bd2 = bd5 = 0; 997 } 998 if (bbe >= 0) 999 bb2 += bbe; 1000 else 1001 bd2 -= bbe; 1002 bs2 = bb2; 1003 j = bbe - scale; 1004 i = j + bbbits - 1; /* logb(rv) */ 1005 if (i < Emin) /* denormal */ 1006 j += P - Emin; 1007 else 1008 j = P + 1 - bbbits; 1009 bb2 += j; 1010 bd2 += j; 1011 bd2 += scale; 1012 i = bb2 < bd2 ? bb2 : bd2; 1013 if (i > bs2) 1014 i = bs2; 1015 if (i > 0) { 1016 bb2 -= i; 1017 bd2 -= i; 1018 bs2 -= i; 1019 } 1020 if (bb5 > 0) { 1021 pow5mult(bs, bb5); 1022 mult(bb, bs); 1023 } 1024 if (bb2 > 0) 1025 lshift(bb, bb2); 1026 if (bd5 > 0) 1027 pow5mult(bd, bd5); 1028 if (bd2 > 0) 1029 lshift(bd, bd2); 1030 if (bs2 > 0) 1031 lshift(bs, bs2); 1032 diff(delta, bb, bd); 1033 dsign = delta.sign; 1034 delta.sign = 0; 1035 i = cmp(delta, bs); 1036 1037 if (i < 0) { 1038 /* Error is less than half an ulp -- check for 1039 * special case of mantissa a power of two. 1040 */ 1041 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask 1042 || (word0(&rv) & Exp_mask) <= (2 * P + 1) * Exp_msk1 1043 ) { 1044 break; 1045 } 1046 if (!delta.words()[0] && delta.size() <= 1) { 1047 /* exact result */ 1048 break; 1049 } 1050 lshift(delta, Log2P); 1051 if (cmp(delta, bs) > 0) 1052 goto dropDown; 1053 break; 1054 } 1055 if (!i) { 1056 /* exactly half-way between */ 1057 if (dsign) { 1058 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 1059 && word1(&rv) == ( 1060 (scale && (y = word0(&rv) & Exp_mask) <= 2 * P * Exp_msk1) 1061 ? (0xffffffff & (0xffffffff << (2 * P + 1 - (y >> Exp_shift)))) : 1062 0xffffffff)) { 1063 /*boundary case -- increment exponent*/ 1064 word0(&rv) = (word0(&rv) & Exp_mask) + Exp_msk1; 1065 word1(&rv) = 0; 1066 dsign = 0; 1067 break; 1068 } 1069 } else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) { 1070 dropDown: 1071 /* boundary case -- decrement exponent */ 1072 if (scale) { 1073 L = word0(&rv) & Exp_mask; 1074 if (L <= (2 * P + 1) * Exp_msk1) { 1075 if (L > (P + 2) * Exp_msk1) 1076 /* round even ==> */ 1077 /* accept rv */ 1078 break; 1079 /* rv = smallest denormal */ 1080 goto undfl; 1081 } 1082 } 1083 L = (word0(&rv) & Exp_mask) - Exp_msk1; 1084 word0(&rv) = L | Bndry_mask1; 1085 word1(&rv) = 0xffffffff; 1086 break; 1087 } 1088 if (!(word1(&rv) & LSB)) 1089 break; 1090 if (dsign) 1091 dval(&rv) += ulp(&rv); 1092 else { 1093 dval(&rv) -= ulp(&rv); 1094 if (!dval(&rv)) 1095 goto undfl; 1096 } 1097 dsign = 1 - dsign; 1098 break; 1099 } 1100 if ((aadj = ratio(delta, bs)) <= 2.) { 1101 if (dsign) 1102 aadj = aadj1 = 1.; 1103 else if (word1(&rv) || word0(&rv) & Bndry_mask) { 1104 if (word1(&rv) == Tiny1 && !word0(&rv)) 1105 goto undfl; 1106 aadj = 1.; 1107 aadj1 = -1.; 1108 } else { 1109 /* special case -- power of FLT_RADIX to be */ 1110 /* rounded down... */ 1111 1112 if (aadj < 2. / FLT_RADIX) 1113 aadj = 1. / FLT_RADIX; 1114 else 1115 aadj *= 0.5; 1116 aadj1 = -aadj; 1117 } 1118 } else { 1119 aadj *= 0.5; 1120 aadj1 = dsign ? aadj : -aadj; 1121 } 1122 y = word0(&rv) & Exp_mask; 1123 1124 /* Check for overflow */ 1125 1126 if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1)) { 1127 dval(&rv0) = dval(&rv); 1128 word0(&rv) -= P * Exp_msk1; 1129 adj.d = aadj1 * ulp(&rv); 1130 dval(&rv) += adj.d; 1131 if ((word0(&rv) & Exp_mask) >= Exp_msk1 * (DBL_MAX_EXP + Bias - P)) { 1132 if (word0(&rv0) == Big0 && word1(&rv0) == Big1) 1133 goto ovfl; 1134 word0(&rv) = Big0; 1135 word1(&rv) = Big1; 1136 goto cont; 1137 } 1138 word0(&rv) += P * Exp_msk1; 1139 } else { 1140 if (scale && y <= 2 * P * Exp_msk1) { 1141 if (aadj <= 0x7fffffff) { 1142 if ((z = (uint32_t)aadj) <= 0) 1143 z = 1; 1144 aadj = z; 1145 aadj1 = dsign ? aadj : -aadj; 1146 } 1147 dval(&aadj2) = aadj1; 1148 word0(&aadj2) += (2 * P + 1) * Exp_msk1 - y; 1149 aadj1 = dval(&aadj2); 1150 } 1151 adj.d = aadj1 * ulp(&rv); 1152 dval(&rv) += adj.d; 1153 } 1154 z = word0(&rv) & Exp_mask; 1155 if (!scale && y == z) { 1156 /* Can we stop now? */ 1157 L = (int32_t)aadj; 1158 aadj -= L; 1159 /* The tolerances below are conservative. */ 1160 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) { 1161 if (aadj < .4999999 || aadj > .5000001) 1162 break; 1163 } else if (aadj < .4999999 / FLT_RADIX) 1164 break; 1165 } 1166 cont: 1167 {} 1168 } 1169 if (scale) { 1170 word0(&rv0) = Exp_1 - 2 * P * Exp_msk1; 1171 word1(&rv0) = 0; 1172 dval(&rv) *= dval(&rv0); 1173 #if HAVE(ERRNO_H) 1174 /* try to avoid the bug of testing an 8087 register value */ 1175 if (!word0(&rv) && !word1(&rv)) 1176 errno = ERANGE; 1177 #endif 1178 } 1179 ret: 1180 if (se) 1181 *se = const_cast<char*>(s); 1182 return sign ? -dval(&rv) : dval(&rv); 1183 } 1184 1185 static ALWAYS_INLINE int quorem(BigInt& b, BigInt& S) 1186 { 1187 size_t n; 1188 uint32_t* bx; 1189 uint32_t* bxe; 1190 uint32_t q; 1191 uint32_t* sx; 1192 uint32_t* sxe; 1193 #ifdef USE_LONG_LONG 1194 unsigned long long borrow, carry, y, ys; 1195 #else 1196 uint32_t borrow, carry, y, ys; 1197 uint32_t si, z, zs; 1198 #endif 1199 ASSERT(b.size() <= 1 || b.words()[b.size() - 1]); 1200 ASSERT(S.size() <= 1 || S.words()[S.size() - 1]); 1201 1202 n = S.size(); 1203 ASSERT_WITH_MESSAGE(b.size() <= n, "oversize b in quorem"); 1204 if (b.size() < n) 1205 return 0; 1206 sx = S.words(); 1207 sxe = sx + --n; 1208 bx = b.words(); 1209 bxe = bx + n; 1210 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 1211 ASSERT_WITH_MESSAGE(q <= 9, "oversized quotient in quorem"); 1212 if (q) { 1213 borrow = 0; 1214 carry = 0; 1215 do { 1216 #ifdef USE_LONG_LONG 1217 ys = *sx++ * (unsigned long long)q + carry; 1218 carry = ys >> 32; 1219 y = *bx - (ys & 0xffffffffUL) - borrow; 1220 borrow = y >> 32 & (uint32_t)1; 1221 *bx++ = (uint32_t)y & 0xffffffffUL; 1222 #else 1223 si = *sx++; 1224 ys = (si & 0xffff) * q + carry; 1225 zs = (si >> 16) * q + (ys >> 16); 1226 carry = zs >> 16; 1227 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 1228 borrow = (y & 0x10000) >> 16; 1229 z = (*bx >> 16) - (zs & 0xffff) - borrow; 1230 borrow = (z & 0x10000) >> 16; 1231 bx = storeInc(bx, z, y); 1232 #endif 1233 } while (sx <= sxe); 1234 if (!*bxe) { 1235 bx = b.words(); 1236 while (--bxe > bx && !*bxe) 1237 --n; 1238 b.resize(n); 1239 } 1240 } 1241 if (cmp(b, S) >= 0) { 1242 q++; 1243 borrow = 0; 1244 carry = 0; 1245 bx = b.words(); 1246 sx = S.words(); 1247 do { 1248 #ifdef USE_LONG_LONG 1249 ys = *sx++ + carry; 1250 carry = ys >> 32; 1251 y = *bx - (ys & 0xffffffffUL) - borrow; 1252 borrow = y >> 32 & (uint32_t)1; 1253 *bx++ = (uint32_t)y & 0xffffffffUL; 1254 #else 1255 si = *sx++; 1256 ys = (si & 0xffff) + carry; 1257 zs = (si >> 16) + (ys >> 16); 1258 carry = zs >> 16; 1259 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 1260 borrow = (y & 0x10000) >> 16; 1261 z = (*bx >> 16) - (zs & 0xffff) - borrow; 1262 borrow = (z & 0x10000) >> 16; 1263 bx = storeInc(bx, z, y); 1264 #endif 1265 } while (sx <= sxe); 1266 bx = b.words(); 1267 bxe = bx + n; 1268 if (!*bxe) { 1269 while (--bxe > bx && !*bxe) 1270 --n; 1271 b.resize(n); 1272 } 1273 } 1274 return q; 1275 } 1276 1277 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 1278 * 1279 * Inspired by "How to Print Floating-Point Numbers Accurately" by 1280 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. 1281 * 1282 * Modifications: 1283 * 1. Rather than iterating, we use a simple numeric overestimate 1284 * to determine k = floor(log10(d)). We scale relevant 1285 * quantities using O(log2(k)) rather than O(k) multiplications. 1286 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 1287 * try to generate digits strictly left to right. Instead, we 1288 * compute with fewer bits and propagate the carry if necessary 1289 * when rounding the final digit up. This is often faster. 1290 * 3. Under the assumption that input will be rounded nearest, 1291 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 1292 * That is, we allow equality in stopping tests when the 1293 * round-nearest rule will give the same floating-point value 1294 * as would satisfaction of the stopping test with strict 1295 * inequality. 1296 * 4. We remove common factors of powers of 2 from relevant 1297 * quantities. 1298 * 5. When converting floating-point integers less than 1e16, 1299 * we use floating-point arithmetic rather than resorting 1300 * to multiple-precision integers. 1301 * 6. When asked to produce fewer than 15 digits, we first try 1302 * to get by with floating-point arithmetic; we resort to 1303 * multiple-precision integer arithmetic only if we cannot 1304 * guarantee that the floating-point calculation has given 1305 * the correctly rounded result. For k requested digits and 1306 * "uniformly" distributed input, the probability is 1307 * something like 10^(k-15) that we must resort to the int32_t 1308 * calculation. 1309 * 1310 * Note: 'leftright' translates to 'generate shortest possible string'. 1311 */ 1312 template<bool roundingNone, bool roundingSignificantFigures, bool roundingDecimalPlaces, bool leftright> 1313 void dtoa(DtoaBuffer result, double dd, int ndigits, bool& signOut, int& exponentOut, unsigned& precisionOut) 1314 { 1315 // Exactly one rounding mode must be specified. 1316 ASSERT(roundingNone + roundingSignificantFigures + roundingDecimalPlaces == 1); 1317 // roundingNone only allowed (only sensible?) with leftright set. 1318 ASSERT(!roundingNone || leftright); 1319 1320 ASSERT(!isnan(dd) && !isinf(dd)); 1321 1322 int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0, 1323 j, j1, k, k0, k_check, m2, m5, s2, s5, 1324 spec_case; 1325 int32_t L; 1326 int denorm; 1327 uint32_t x; 1328 BigInt b, delta, mlo, mhi, S; 1329 U d2, eps, u; 1330 double ds; 1331 char* s; 1332 char* s0; 1333 1334 u.d = dd; 1335 1336 /* Infinity or NaN */ 1337 ASSERT((word0(&u) & Exp_mask) != Exp_mask); 1338 1339 // JavaScript toString conversion treats -0 as 0. 1340 if (!dval(&u)) { 1341 signOut = false; 1342 exponentOut = 0; 1343 precisionOut = 1; 1344 result[0] = '0'; 1345 result[1] = '\0'; 1346 return; 1347 } 1348 1349 if (word0(&u) & Sign_bit) { 1350 signOut = true; 1351 word0(&u) &= ~Sign_bit; // clear sign bit 1352 } else 1353 signOut = false; 1354 1355 d2b(b, &u, &be, &bbits); 1356 if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask >> Exp_shift1)))) { 1357 dval(&d2) = dval(&u); 1358 word0(&d2) &= Frac_mask1; 1359 word0(&d2) |= Exp_11; 1360 1361 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 1362 * log10(x) = log(x) / log(10) 1363 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 1364 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 1365 * 1366 * This suggests computing an approximation k to log10(d) by 1367 * 1368 * k = (i - Bias)*0.301029995663981 1369 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 1370 * 1371 * We want k to be too large rather than too small. 1372 * The error in the first-order Taylor series approximation 1373 * is in our favor, so we just round up the constant enough 1374 * to compensate for any error in the multiplication of 1375 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 1376 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 1377 * adding 1e-13 to the constant term more than suffices. 1378 * Hence we adjust the constant term to 0.1760912590558. 1379 * (We could get a more accurate k by invoking log10, 1380 * but this is probably not worthwhile.) 1381 */ 1382 1383 i -= Bias; 1384 denorm = 0; 1385 } else { 1386 /* d is denormalized */ 1387 1388 i = bbits + be + (Bias + (P - 1) - 1); 1389 x = (i > 32) ? (word0(&u) << (64 - i)) | (word1(&u) >> (i - 32)) 1390 : word1(&u) << (32 - i); 1391 dval(&d2) = x; 1392 word0(&d2) -= 31 * Exp_msk1; /* adjust exponent */ 1393 i -= (Bias + (P - 1) - 1) + 1; 1394 denorm = 1; 1395 } 1396 ds = (dval(&d2) - 1.5) * 0.289529654602168 + 0.1760912590558 + (i * 0.301029995663981); 1397 k = (int)ds; 1398 if (ds < 0. && ds != k) 1399 k--; /* want k = floor(ds) */ 1400 k_check = 1; 1401 if (k >= 0 && k <= Ten_pmax) { 1402 if (dval(&u) < tens[k]) 1403 k--; 1404 k_check = 0; 1405 } 1406 j = bbits - i - 1; 1407 if (j >= 0) { 1408 b2 = 0; 1409 s2 = j; 1410 } else { 1411 b2 = -j; 1412 s2 = 0; 1413 } 1414 if (k >= 0) { 1415 b5 = 0; 1416 s5 = k; 1417 s2 += k; 1418 } else { 1419 b2 -= k; 1420 b5 = -k; 1421 s5 = 0; 1422 } 1423 1424 if (roundingNone) { 1425 ilim = ilim1 = -1; 1426 i = 18; 1427 ndigits = 0; 1428 } 1429 if (roundingSignificantFigures) { 1430 if (ndigits <= 0) 1431 ndigits = 1; 1432 ilim = ilim1 = i = ndigits; 1433 } 1434 if (roundingDecimalPlaces) { 1435 i = ndigits + k + 1; 1436 ilim = i; 1437 ilim1 = i - 1; 1438 if (i <= 0) 1439 i = 1; 1440 } 1441 1442 s = s0 = result; 1443 1444 if (ilim >= 0 && ilim <= Quick_max) { 1445 /* Try to get by with floating-point arithmetic. */ 1446 1447 i = 0; 1448 dval(&d2) = dval(&u); 1449 k0 = k; 1450 ilim0 = ilim; 1451 ieps = 2; /* conservative */ 1452 if (k > 0) { 1453 ds = tens[k & 0xf]; 1454 j = k >> 4; 1455 if (j & Bletch) { 1456 /* prevent overflows */ 1457 j &= Bletch - 1; 1458 dval(&u) /= bigtens[n_bigtens - 1]; 1459 ieps++; 1460 } 1461 for (; j; j >>= 1, i++) { 1462 if (j & 1) { 1463 ieps++; 1464 ds *= bigtens[i]; 1465 } 1466 } 1467 dval(&u) /= ds; 1468 } else if ((j1 = -k)) { 1469 dval(&u) *= tens[j1 & 0xf]; 1470 for (j = j1 >> 4; j; j >>= 1, i++) { 1471 if (j & 1) { 1472 ieps++; 1473 dval(&u) *= bigtens[i]; 1474 } 1475 } 1476 } 1477 if (k_check && dval(&u) < 1. && ilim > 0) { 1478 if (ilim1 <= 0) 1479 goto fastFailed; 1480 ilim = ilim1; 1481 k--; 1482 dval(&u) *= 10.; 1483 ieps++; 1484 } 1485 dval(&eps) = (ieps * dval(&u)) + 7.; 1486 word0(&eps) -= (P - 1) * Exp_msk1; 1487 if (!ilim) { 1488 S.clear(); 1489 mhi.clear(); 1490 dval(&u) -= 5.; 1491 if (dval(&u) > dval(&eps)) 1492 goto oneDigit; 1493 if (dval(&u) < -dval(&eps)) 1494 goto noDigits; 1495 goto fastFailed; 1496 } 1497 if (leftright) { 1498 /* Use Steele & White method of only 1499 * generating digits needed. 1500 */ 1501 dval(&eps) = (0.5 / tens[ilim - 1]) - dval(&eps); 1502 for (i = 0;;) { 1503 L = (long int)dval(&u); 1504 dval(&u) -= L; 1505 *s++ = '0' + (int)L; 1506 if (dval(&u) < dval(&eps)) 1507 goto ret; 1508 if (1. - dval(&u) < dval(&eps)) 1509 goto bumpUp; 1510 if (++i >= ilim) 1511 break; 1512 dval(&eps) *= 10.; 1513 dval(&u) *= 10.; 1514 } 1515 } else { 1516 /* Generate ilim digits, then fix them up. */ 1517 dval(&eps) *= tens[ilim - 1]; 1518 for (i = 1;; i++, dval(&u) *= 10.) { 1519 L = (int32_t)(dval(&u)); 1520 if (!(dval(&u) -= L)) 1521 ilim = i; 1522 *s++ = '0' + (int)L; 1523 if (i == ilim) { 1524 if (dval(&u) > 0.5 + dval(&eps)) 1525 goto bumpUp; 1526 if (dval(&u) < 0.5 - dval(&eps)) { 1527 while (*--s == '0') { } 1528 s++; 1529 goto ret; 1530 } 1531 break; 1532 } 1533 } 1534 } 1535 fastFailed: 1536 s = s0; 1537 dval(&u) = dval(&d2); 1538 k = k0; 1539 ilim = ilim0; 1540 } 1541 1542 /* Do we have a "small" integer? */ 1543 1544 if (be >= 0 && k <= Int_max) { 1545 /* Yes. */ 1546 ds = tens[k]; 1547 if (ndigits < 0 && ilim <= 0) { 1548 S.clear(); 1549 mhi.clear(); 1550 if (ilim < 0 || dval(&u) <= 5 * ds) 1551 goto noDigits; 1552 goto oneDigit; 1553 } 1554 for (i = 1;; i++, dval(&u) *= 10.) { 1555 L = (int32_t)(dval(&u) / ds); 1556 dval(&u) -= L * ds; 1557 *s++ = '0' + (int)L; 1558 if (!dval(&u)) { 1559 break; 1560 } 1561 if (i == ilim) { 1562 dval(&u) += dval(&u); 1563 if (dval(&u) > ds || (dval(&u) == ds && (L & 1))) { 1564 bumpUp: 1565 while (*--s == '9') 1566 if (s == s0) { 1567 k++; 1568 *s = '0'; 1569 break; 1570 } 1571 ++*s++; 1572 } 1573 break; 1574 } 1575 } 1576 goto ret; 1577 } 1578 1579 m2 = b2; 1580 m5 = b5; 1581 mhi.clear(); 1582 mlo.clear(); 1583 if (leftright) { 1584 i = denorm ? be + (Bias + (P - 1) - 1 + 1) : 1 + P - bbits; 1585 b2 += i; 1586 s2 += i; 1587 i2b(mhi, 1); 1588 } 1589 if (m2 > 0 && s2 > 0) { 1590 i = m2 < s2 ? m2 : s2; 1591 b2 -= i; 1592 m2 -= i; 1593 s2 -= i; 1594 } 1595 if (b5 > 0) { 1596 if (leftright) { 1597 if (m5 > 0) { 1598 pow5mult(mhi, m5); 1599 mult(b, mhi); 1600 } 1601 if ((j = b5 - m5)) 1602 pow5mult(b, j); 1603 } else 1604 pow5mult(b, b5); 1605 } 1606 i2b(S, 1); 1607 if (s5 > 0) 1608 pow5mult(S, s5); 1609 1610 /* Check for special case that d is a normalized power of 2. */ 1611 1612 spec_case = 0; 1613 if ((roundingNone || leftright) && (!word1(&u) && !(word0(&u) & Bndry_mask) && word0(&u) & (Exp_mask & ~Exp_msk1))) { 1614 /* The special case */ 1615 b2 += Log2P; 1616 s2 += Log2P; 1617 spec_case = 1; 1618 } 1619 1620 /* Arrange for convenient computation of quotients: 1621 * shift left if necessary so divisor has 4 leading 0 bits. 1622 * 1623 * Perhaps we should just compute leading 28 bits of S once 1624 * and for all and pass them and a shift to quorem, so it 1625 * can do shifts and ors to compute the numerator for q. 1626 */ 1627 if ((i = ((s5 ? 32 - hi0bits(S.words()[S.size() - 1]) : 1) + s2) & 0x1f)) 1628 i = 32 - i; 1629 if (i > 4) { 1630 i -= 4; 1631 b2 += i; 1632 m2 += i; 1633 s2 += i; 1634 } else if (i < 4) { 1635 i += 28; 1636 b2 += i; 1637 m2 += i; 1638 s2 += i; 1639 } 1640 if (b2 > 0) 1641 lshift(b, b2); 1642 if (s2 > 0) 1643 lshift(S, s2); 1644 if (k_check) { 1645 if (cmp(b, S) < 0) { 1646 k--; 1647 multadd(b, 10, 0); /* we botched the k estimate */ 1648 if (leftright) 1649 multadd(mhi, 10, 0); 1650 ilim = ilim1; 1651 } 1652 } 1653 if (ilim <= 0 && roundingDecimalPlaces) { 1654 if (ilim < 0) 1655 goto noDigits; 1656 multadd(S, 5, 0); 1657 // For IEEE-754 unbiased rounding this check should be <=, such that 0.5 would flush to zero. 1658 if (cmp(b, S) < 0) 1659 goto noDigits; 1660 goto oneDigit; 1661 } 1662 if (leftright) { 1663 if (m2 > 0) 1664 lshift(mhi, m2); 1665 1666 /* Compute mlo -- check for special case 1667 * that d is a normalized power of 2. 1668 */ 1669 1670 mlo = mhi; 1671 if (spec_case) 1672 lshift(mhi, Log2P); 1673 1674 for (i = 1;;i++) { 1675 dig = quorem(b, S) + '0'; 1676 /* Do we yet have the shortest decimal string 1677 * that will round to d? 1678 */ 1679 j = cmp(b, mlo); 1680 diff(delta, S, mhi); 1681 j1 = delta.sign ? 1 : cmp(b, delta); 1682 #ifdef DTOA_ROUND_BIASED 1683 if (j < 0 || !j) { 1684 #else 1685 // FIXME: ECMA-262 specifies that equidistant results round away from 1686 // zero, which probably means we shouldn't be on the unbiased code path 1687 // (the (word1(&u) & 1) clause is looking highly suspicious). I haven't 1688 // yet understood this code well enough to make the call, but we should 1689 // probably be enabling DTOA_ROUND_BIASED. I think the interesting corner 1690 // case to understand is probably "Math.pow(0.5, 24).toString()". 1691 // I believe this value is interesting because I think it is precisely 1692 // representable in binary floating point, and its decimal representation 1693 // has a single digit that Steele & White reduction can remove, with the 1694 // value 5 (thus equidistant from the next numbers above and below). 1695 // We produce the correct answer using either codepath, and I don't as 1696 // yet understand why. :-) 1697 if (!j1 && !(word1(&u) & 1)) { 1698 if (dig == '9') 1699 goto round9up; 1700 if (j > 0) 1701 dig++; 1702 *s++ = dig; 1703 goto ret; 1704 } 1705 if (j < 0 || (!j && !(word1(&u) & 1))) { 1706 #endif 1707 if ((b.words()[0] || b.size() > 1) && (j1 > 0)) { 1708 lshift(b, 1); 1709 j1 = cmp(b, S); 1710 // For IEEE-754 round-to-even, this check should be (j1 > 0 || (!j1 && (dig & 1))), 1711 // but ECMA-262 specifies that equidistant values (e.g. (.5).toFixed()) should 1712 // be rounded away from zero. 1713 if (j1 >= 0) { 1714 if (dig == '9') 1715 goto round9up; 1716 dig++; 1717 } 1718 } 1719 *s++ = dig; 1720 goto ret; 1721 } 1722 if (j1 > 0) { 1723 if (dig == '9') { /* possible if i == 1 */ 1724 round9up: 1725 *s++ = '9'; 1726 goto roundoff; 1727 } 1728 *s++ = dig + 1; 1729 goto ret; 1730 } 1731 *s++ = dig; 1732 if (i == ilim) 1733 break; 1734 multadd(b, 10, 0); 1735 multadd(mlo, 10, 0); 1736 multadd(mhi, 10, 0); 1737 } 1738 } else { 1739 for (i = 1;; i++) { 1740 *s++ = dig = quorem(b, S) + '0'; 1741 if (!b.words()[0] && b.size() <= 1) 1742 goto ret; 1743 if (i >= ilim) 1744 break; 1745 multadd(b, 10, 0); 1746 } 1747 } 1748 1749 /* Round off last digit */ 1750 1751 lshift(b, 1); 1752 j = cmp(b, S); 1753 // For IEEE-754 round-to-even, this check should be (j > 0 || (!j && (dig & 1))), 1754 // but ECMA-262 specifies that equidistant values (e.g. (.5).toFixed()) should 1755 // be rounded away from zero. 1756 if (j >= 0) { 1757 roundoff: 1758 while (*--s == '9') 1759 if (s == s0) { 1760 k++; 1761 *s++ = '1'; 1762 goto ret; 1763 } 1764 ++*s++; 1765 } else { 1766 while (*--s == '0') { } 1767 s++; 1768 } 1769 goto ret; 1770 noDigits: 1771 exponentOut = 0; 1772 precisionOut = 1; 1773 result[0] = '0'; 1774 result[1] = '\0'; 1775 return; 1776 oneDigit: 1777 *s++ = '1'; 1778 k++; 1779 goto ret; 1780 ret: 1781 ASSERT(s > result); 1782 *s = 0; 1783 exponentOut = k; 1784 precisionOut = s - result; 1785 } 1786 1787 void dtoa(DtoaBuffer result, double dd, bool& sign, int& exponent, unsigned& precision) 1788 { 1789 // flags are roundingNone, leftright. 1790 dtoa<true, false, false, true>(result, dd, 0, sign, exponent, precision); 1791 } 1792 1793 void dtoaRoundSF(DtoaBuffer result, double dd, int ndigits, bool& sign, int& exponent, unsigned& precision) 1794 { 1795 // flag is roundingSignificantFigures. 1796 dtoa<false, true, false, false>(result, dd, ndigits, sign, exponent, precision); 1797 } 1798 1799 void dtoaRoundDP(DtoaBuffer result, double dd, int ndigits, bool& sign, int& exponent, unsigned& precision) 1800 { 1801 // flag is roundingDecimalPlaces. 1802 dtoa<false, false, true, false>(result, dd, ndigits, sign, exponent, precision); 1803 } 1804 1805 static ALWAYS_INLINE void copyAsciiToUTF16(UChar* next, const char* src, unsigned size) 1806 { 1807 for (unsigned i = 0; i < size; ++i) 1808 *next++ = *src++; 1809 } 1810 1811 unsigned numberToString(double d, NumberToStringBuffer buffer) 1812 { 1813 // Handle NaN and Infinity. 1814 if (isnan(d) || isinf(d)) { 1815 if (isnan(d)) { 1816 copyAsciiToUTF16(buffer, "NaN", 3); 1817 return 3; 1818 } 1819 if (d > 0) { 1820 copyAsciiToUTF16(buffer, "Infinity", 8); 1821 return 8; 1822 } 1823 copyAsciiToUTF16(buffer, "-Infinity", 9); 1824 return 9; 1825 } 1826 1827 // Convert to decimal with rounding. 1828 DecimalNumber number(d); 1829 return number.exponent() >= -6 && number.exponent() < 21 1830 ? number.toStringDecimal(buffer, NumberToStringBufferLength) 1831 : number.toStringExponential(buffer, NumberToStringBufferLength); 1832 } 1833 1834 } // namespace WTF 1835
