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 */
storeInc(uint32_t * p,uint16_t high,uint16_t low)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 {
BigIntWTF::BigInt166 BigInt() : sign(0) { }
167 int sign;
168
clearWTF::BigInt169 void clear()
170 {
171 sign = 0;
172 m_words.clear();
173 }
174
sizeWTF::BigInt175 size_t size() const
176 {
177 return m_words.size();
178 }
179
resizeWTF::BigInt180 void resize(size_t s)
181 {
182 m_words.resize(s);
183 }
184
wordsWTF::BigInt185 uint32_t* words()
186 {
187 return m_words.data();
188 }
189
wordsWTF::BigInt190 const uint32_t* words() const
191 {
192 return m_words.data();
193 }
194
appendWTF::BigInt195 void append(uint32_t w)
196 {
197 m_words.append(w);
198 }
199
200 Vector<uint32_t, 16> m_words;
201 };
202
multadd(BigInt & b,int m,int a)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
s2b(BigInt & b,const char * s,int nd0,int nd,uint32_t y9)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
hi0bits(uint32_t x)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
lo0bits(uint32_t * y)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
i2b(BigInt & b,int i)322 static void i2b(BigInt& b, int i)
323 {
324 b.sign = 0;
325 b.resize(1);
326 b.words()[0] = i;
327 }
328
mult(BigInt & aRef,const BigInt & bRef)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:
P5NodeWTF::P5Node420 P5Node() { }
421 BigInt val;
422 P5Node* next;
423 };
424
425 static P5Node* p5s;
426 static int p5sCount;
427
pow5mult(BigInt & b,int k)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
lshift(BigInt & b,int k)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
cmp(const BigInt & a,const BigInt & b)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
diff(BigInt & c,const BigInt & aRef,const BigInt & bRef)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
ulp(U * x)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
b2d(const BigInt & a,int * e)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
d2b(BigInt & b,U * d,int * e,int * bits)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
ratio(const BigInt & a,const BigInt & b)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
strtod(const char * s00,char ** se)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
quorem(BigInt & b,BigInt & S)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>
dtoa(DtoaBuffer result,double dd,int ndigits,bool & signOut,int & exponentOut,unsigned & precisionOut)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