Home | History | Annotate | Download | only in wtf
      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, 2012 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 #include "config.h"
     36 #include "dtoa.h"
     37 
     38 #include "wtf/CPU.h"
     39 #include "wtf/MathExtras.h"
     40 #include "wtf/ThreadingPrimitives.h"
     41 #include "wtf/Vector.h"
     42 
     43 #if COMPILER(MSVC)
     44 #pragma warning(disable: 4244)
     45 #pragma warning(disable: 4245)
     46 #pragma warning(disable: 4554)
     47 
     48 #if _MSC_VER == 1800
     49 // TODO(scottmg): VS2013 currently ICEs on a bunch of functions in this file.
     50 // Upstream bug fixed in next release. See http://crbug.com/288498.
     51 #pragma optimize("", off)
     52 #endif
     53 
     54 #endif
     55 
     56 namespace WTF {
     57 
     58 Mutex* s_dtoaP5Mutex;
     59 
     60 typedef union {
     61     double d;
     62     uint32_t L[2];
     63 } U;
     64 
     65 #if CPU(BIG_ENDIAN) || CPU(MIDDLE_ENDIAN)
     66 #define word0(x) (x)->L[0]
     67 #define word1(x) (x)->L[1]
     68 #else
     69 #define word0(x) (x)->L[1]
     70 #define word1(x) (x)->L[0]
     71 #endif
     72 #define dval(x) (x)->d
     73 
     74 #define Exp_shift  20
     75 #define Exp_shift1 20
     76 #define Exp_msk1    0x100000
     77 #define Exp_msk11   0x100000
     78 #define Exp_mask  0x7ff00000
     79 #define P 53
     80 #define Bias 1023
     81 #define Emin (-1022)
     82 #define Exp_1  0x3ff00000
     83 #define Exp_11 0x3ff00000
     84 #define Ebits 11
     85 #define Frac_mask  0xfffff
     86 #define Frac_mask1 0xfffff
     87 #define Ten_pmax 22
     88 #define Bletch 0x10
     89 #define Bndry_mask  0xfffff
     90 #define Bndry_mask1 0xfffff
     91 #define LSB 1
     92 #define Sign_bit 0x80000000
     93 #define Log2P 1
     94 #define Tiny0 0
     95 #define Tiny1 1
     96 #define Quick_max 14
     97 #define Int_max 14
     98 
     99 #define rounded_product(a, b) a *= b
    100 #define rounded_quotient(a, b) a /= b
    101 
    102 #define Big0 (Frac_mask1 | Exp_msk1 * (DBL_MAX_EXP + Bias - 1))
    103 #define Big1 0xffffffff
    104 
    105 #if CPU(X86_64)
    106 // FIXME: should we enable this on all 64-bit CPUs?
    107 // 64-bit emulation provided by the compiler is likely to be slower than dtoa own code on 32-bit hardware.
    108 #define USE_LONG_LONG
    109 #endif
    110 
    111 #ifndef USE_LONG_LONG
    112 /* The following definition of Storeinc is appropriate for MIPS processors.
    113  * An alternative that might be better on some machines is
    114  *  *p++ = high << 16 | low & 0xffff;
    115  */
    116 static ALWAYS_INLINE uint32_t* storeInc(uint32_t* p, uint16_t high, uint16_t low)
    117 {
    118     uint16_t* p16 = reinterpret_cast<uint16_t*>(p);
    119 #if CPU(BIG_ENDIAN)
    120     p16[0] = high;
    121     p16[1] = low;
    122 #else
    123     p16[1] = high;
    124     p16[0] = low;
    125 #endif
    126     return p + 1;
    127 }
    128 #endif
    129 
    130 struct BigInt {
    131     BigInt() : sign(0) { }
    132     int sign;
    133 
    134     void clear()
    135     {
    136         sign = 0;
    137         m_words.clear();
    138     }
    139 
    140     size_t size() const
    141     {
    142         return m_words.size();
    143     }
    144 
    145     void resize(size_t s)
    146     {
    147         m_words.resize(s);
    148     }
    149 
    150     uint32_t* words()
    151     {
    152         return m_words.data();
    153     }
    154 
    155     const uint32_t* words() const
    156     {
    157         return m_words.data();
    158     }
    159 
    160     void append(uint32_t w)
    161     {
    162         m_words.append(w);
    163     }
    164 
    165     Vector<uint32_t, 16> m_words;
    166 };
    167 
    168 static void multadd(BigInt& b, int m, int a)    /* multiply by m and add a */
    169 {
    170 #ifdef USE_LONG_LONG
    171     unsigned long long carry;
    172 #else
    173     uint32_t carry;
    174 #endif
    175 
    176     int wds = b.size();
    177     uint32_t* x = b.words();
    178     int i = 0;
    179     carry = a;
    180     do {
    181 #ifdef USE_LONG_LONG
    182         unsigned long long y = *x * (unsigned long long)m + carry;
    183         carry = y >> 32;
    184         *x++ = (uint32_t)y & 0xffffffffUL;
    185 #else
    186         uint32_t xi = *x;
    187         uint32_t y = (xi & 0xffff) * m + carry;
    188         uint32_t z = (xi >> 16) * m + (y >> 16);
    189         carry = z >> 16;
    190         *x++ = (z << 16) + (y & 0xffff);
    191 #endif
    192     } while (++i < wds);
    193 
    194     if (carry)
    195         b.append((uint32_t)carry);
    196 }
    197 
    198 static int hi0bits(uint32_t x)
    199 {
    200     int k = 0;
    201 
    202     if (!(x & 0xffff0000)) {
    203         k = 16;
    204         x <<= 16;
    205     }
    206     if (!(x & 0xff000000)) {
    207         k += 8;
    208         x <<= 8;
    209     }
    210     if (!(x & 0xf0000000)) {
    211         k += 4;
    212         x <<= 4;
    213     }
    214     if (!(x & 0xc0000000)) {
    215         k += 2;
    216         x <<= 2;
    217     }
    218     if (!(x & 0x80000000)) {
    219         k++;
    220         if (!(x & 0x40000000))
    221             return 32;
    222     }
    223     return k;
    224 }
    225 
    226 static int lo0bits(uint32_t* y)
    227 {
    228     int k;
    229     uint32_t x = *y;
    230 
    231     if (x & 7) {
    232         if (x & 1)
    233             return 0;
    234         if (x & 2) {
    235             *y = x >> 1;
    236             return 1;
    237         }
    238         *y = x >> 2;
    239         return 2;
    240     }
    241     k = 0;
    242     if (!(x & 0xffff)) {
    243         k = 16;
    244         x >>= 16;
    245     }
    246     if (!(x & 0xff)) {
    247         k += 8;
    248         x >>= 8;
    249     }
    250     if (!(x & 0xf)) {
    251         k += 4;
    252         x >>= 4;
    253     }
    254     if (!(x & 0x3)) {
    255         k += 2;
    256         x >>= 2;
    257     }
    258     if (!(x & 1)) {
    259         k++;
    260         x >>= 1;
    261         if (!x)
    262             return 32;
    263     }
    264     *y = x;
    265     return k;
    266 }
    267 
    268 static void i2b(BigInt& b, int i)
    269 {
    270     b.sign = 0;
    271     b.resize(1);
    272     b.words()[0] = i;
    273 }
    274 
    275 static void mult(BigInt& aRef, const BigInt& bRef)
    276 {
    277     const BigInt* a = &aRef;
    278     const BigInt* b = &bRef;
    279     BigInt c;
    280     int wa, wb, wc;
    281     const uint32_t* x = 0;
    282     const uint32_t* xa;
    283     const uint32_t* xb;
    284     const uint32_t* xae;
    285     const uint32_t* xbe;
    286     uint32_t* xc;
    287     uint32_t* xc0;
    288     uint32_t y;
    289 #ifdef USE_LONG_LONG
    290     unsigned long long carry, z;
    291 #else
    292     uint32_t carry, z;
    293 #endif
    294 
    295     if (a->size() < b->size()) {
    296         const BigInt* tmp = a;
    297         a = b;
    298         b = tmp;
    299     }
    300 
    301     wa = a->size();
    302     wb = b->size();
    303     wc = wa + wb;
    304     c.resize(wc);
    305 
    306     for (xc = c.words(), xa = xc + wc; xc < xa; xc++)
    307         *xc = 0;
    308     xa = a->words();
    309     xae = xa + wa;
    310     xb = b->words();
    311     xbe = xb + wb;
    312     xc0 = c.words();
    313 #ifdef USE_LONG_LONG
    314     for (; xb < xbe; xc0++) {
    315         if ((y = *xb++)) {
    316             x = xa;
    317             xc = xc0;
    318             carry = 0;
    319             do {
    320                 z = *x++ * (unsigned long long)y + *xc + carry;
    321                 carry = z >> 32;
    322                 *xc++ = (uint32_t)z & 0xffffffffUL;
    323             } while (x < xae);
    324             *xc = (uint32_t)carry;
    325         }
    326     }
    327 #else
    328     for (; xb < xbe; xb++, xc0++) {
    329         if ((y = *xb & 0xffff)) {
    330             x = xa;
    331             xc = xc0;
    332             carry = 0;
    333             do {
    334                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
    335                 carry = z >> 16;
    336                 uint32_t z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
    337                 carry = z2 >> 16;
    338                 xc = storeInc(xc, z2, z);
    339             } while (x < xae);
    340             *xc = carry;
    341         }
    342         if ((y = *xb >> 16)) {
    343             x = xa;
    344             xc = xc0;
    345             carry = 0;
    346             uint32_t z2 = *xc;
    347             do {
    348                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
    349                 carry = z >> 16;
    350                 xc = storeInc(xc, z, z2);
    351                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
    352                 carry = z2 >> 16;
    353             } while (x < xae);
    354             *xc = z2;
    355         }
    356     }
    357 #endif
    358     for (xc0 = c.words(), xc = xc0 + wc; wc > 0 && !*--xc; --wc) { }
    359     c.resize(wc);
    360     aRef = c;
    361 }
    362 
    363 struct P5Node {
    364     WTF_MAKE_NONCOPYABLE(P5Node); WTF_MAKE_FAST_ALLOCATED;
    365 public:
    366     P5Node() { }
    367     BigInt val;
    368     P5Node* next;
    369 };
    370 
    371 static P5Node* p5s;
    372 static int p5sCount;
    373 
    374 static ALWAYS_INLINE void pow5mult(BigInt& b, int k)
    375 {
    376     static int p05[3] = { 5, 25, 125 };
    377 
    378     if (int i = k & 3)
    379         multadd(b, p05[i - 1], 0);
    380 
    381     if (!(k >>= 2))
    382         return;
    383 
    384     s_dtoaP5Mutex->lock();
    385     P5Node* p5 = p5s;
    386 
    387     if (!p5) {
    388         /* first time */
    389         p5 = new P5Node;
    390         i2b(p5->val, 625);
    391         p5->next = 0;
    392         p5s = p5;
    393         p5sCount = 1;
    394     }
    395 
    396     int p5sCountLocal = p5sCount;
    397     s_dtoaP5Mutex->unlock();
    398     int p5sUsed = 0;
    399 
    400     for (;;) {
    401         if (k & 1)
    402             mult(b, p5->val);
    403 
    404         if (!(k >>= 1))
    405             break;
    406 
    407         if (++p5sUsed == p5sCountLocal) {
    408             s_dtoaP5Mutex->lock();
    409             if (p5sUsed == p5sCount) {
    410                 ASSERT(!p5->next);
    411                 p5->next = new P5Node;
    412                 p5->next->next = 0;
    413                 p5->next->val = p5->val;
    414                 mult(p5->next->val, p5->next->val);
    415                 ++p5sCount;
    416             }
    417 
    418             p5sCountLocal = p5sCount;
    419             s_dtoaP5Mutex->unlock();
    420         }
    421         p5 = p5->next;
    422     }
    423 }
    424 
    425 static ALWAYS_INLINE void lshift(BigInt& b, int k)
    426 {
    427     int n = k >> 5;
    428 
    429     int origSize = b.size();
    430     int n1 = n + origSize + 1;
    431 
    432     if (k &= 0x1f)
    433         b.resize(b.size() + n + 1);
    434     else
    435         b.resize(b.size() + n);
    436 
    437     const uint32_t* srcStart = b.words();
    438     uint32_t* dstStart = b.words();
    439     const uint32_t* src = srcStart + origSize - 1;
    440     uint32_t* dst = dstStart + n1 - 1;
    441     if (k) {
    442         uint32_t hiSubword = 0;
    443         int s = 32 - k;
    444         for (; src >= srcStart; --src) {
    445             *dst-- = hiSubword | *src >> s;
    446             hiSubword = *src << k;
    447         }
    448         *dst = hiSubword;
    449         ASSERT(dst == dstStart + n);
    450 
    451         b.resize(origSize + n + !!b.words()[n1 - 1]);
    452     }
    453     else {
    454         do {
    455             *--dst = *src--;
    456         } while (src >= srcStart);
    457     }
    458     for (dst = dstStart + n; dst != dstStart; )
    459         *--dst = 0;
    460 
    461     ASSERT(b.size() <= 1 || b.words()[b.size() - 1]);
    462 }
    463 
    464 static int cmp(const BigInt& a, const BigInt& b)
    465 {
    466     const uint32_t *xa, *xa0, *xb, *xb0;
    467     int i, j;
    468 
    469     i = a.size();
    470     j = b.size();
    471     ASSERT(i <= 1 || a.words()[i - 1]);
    472     ASSERT(j <= 1 || b.words()[j - 1]);
    473     if (i -= j)
    474         return i;
    475     xa0 = a.words();
    476     xa = xa0 + j;
    477     xb0 = b.words();
    478     xb = xb0 + j;
    479     for (;;) {
    480         if (*--xa != *--xb)
    481             return *xa < *xb ? -1 : 1;
    482         if (xa <= xa0)
    483             break;
    484     }
    485     return 0;
    486 }
    487 
    488 static ALWAYS_INLINE void diff(BigInt& c, const BigInt& aRef, const BigInt& bRef)
    489 {
    490     const BigInt* a = &aRef;
    491     const BigInt* b = &bRef;
    492     int i, wa, wb;
    493     uint32_t* xc;
    494 
    495     i = cmp(*a, *b);
    496     if (!i) {
    497         c.sign = 0;
    498         c.resize(1);
    499         c.words()[0] = 0;
    500         return;
    501     }
    502     if (i < 0) {
    503         const BigInt* tmp = a;
    504         a = b;
    505         b = tmp;
    506         i = 1;
    507     } else
    508         i = 0;
    509 
    510     wa = a->size();
    511     const uint32_t* xa = a->words();
    512     const uint32_t* xae = xa + wa;
    513     wb = b->size();
    514     const uint32_t* xb = b->words();
    515     const uint32_t* xbe = xb + wb;
    516 
    517     c.resize(wa);
    518     c.sign = i;
    519     xc = c.words();
    520 #ifdef USE_LONG_LONG
    521     unsigned long long borrow = 0;
    522     do {
    523         unsigned long long y = (unsigned long long)*xa++ - *xb++ - borrow;
    524         borrow = y >> 32 & (uint32_t)1;
    525         *xc++ = (uint32_t)y & 0xffffffffUL;
    526     } while (xb < xbe);
    527     while (xa < xae) {
    528         unsigned long long y = *xa++ - borrow;
    529         borrow = y >> 32 & (uint32_t)1;
    530         *xc++ = (uint32_t)y & 0xffffffffUL;
    531     }
    532 #else
    533     uint32_t borrow = 0;
    534     do {
    535         uint32_t y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
    536         borrow = (y & 0x10000) >> 16;
    537         uint32_t z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
    538         borrow = (z & 0x10000) >> 16;
    539         xc = storeInc(xc, z, y);
    540     } while (xb < xbe);
    541     while (xa < xae) {
    542         uint32_t y = (*xa & 0xffff) - borrow;
    543         borrow = (y & 0x10000) >> 16;
    544         uint32_t z = (*xa++ >> 16) - borrow;
    545         borrow = (z & 0x10000) >> 16;
    546         xc = storeInc(xc, z, y);
    547     }
    548 #endif
    549     while (!*--xc)
    550         wa--;
    551     c.resize(wa);
    552 }
    553 
    554 static ALWAYS_INLINE void d2b(BigInt& b, U* d, int* e, int* bits)
    555 {
    556     int de, k;
    557     uint32_t* x;
    558     uint32_t y, z;
    559     int i;
    560 #define d0 word0(d)
    561 #define d1 word1(d)
    562 
    563     b.sign = 0;
    564     b.resize(1);
    565     x = b.words();
    566 
    567     z = d0 & Frac_mask;
    568     d0 &= 0x7fffffff;    /* clear sign bit, which we ignore */
    569     if ((de = (int)(d0 >> Exp_shift)))
    570         z |= Exp_msk1;
    571     if ((y = d1)) {
    572         if ((k = lo0bits(&y))) {
    573             x[0] = y | (z << (32 - k));
    574             z >>= k;
    575         } else
    576             x[0] = y;
    577         if (z) {
    578             b.resize(2);
    579             x[1] = z;
    580         }
    581 
    582         i = b.size();
    583     } else {
    584         k = lo0bits(&z);
    585         x[0] = z;
    586         i = 1;
    587         b.resize(1);
    588         k += 32;
    589     }
    590     if (de) {
    591         *e = de - Bias - (P - 1) + k;
    592         *bits = P - k;
    593     } else {
    594         *e = 0 - Bias - (P - 1) + 1 + k;
    595         *bits = (32 * i) - hi0bits(x[i - 1]);
    596     }
    597 }
    598 #undef d0
    599 #undef d1
    600 
    601 static const double tens[] = {
    602     1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
    603     1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
    604     1e20, 1e21, 1e22
    605 };
    606 
    607 static const double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
    608 
    609 #define Scale_Bit 0x10
    610 #define n_bigtens 5
    611 
    612 static ALWAYS_INLINE int quorem(BigInt& b, BigInt& S)
    613 {
    614     size_t n;
    615     uint32_t* bx;
    616     uint32_t* bxe;
    617     uint32_t q;
    618     uint32_t* sx;
    619     uint32_t* sxe;
    620 #ifdef USE_LONG_LONG
    621     unsigned long long borrow, carry, y, ys;
    622 #else
    623     uint32_t borrow, carry, y, ys;
    624     uint32_t si, z, zs;
    625 #endif
    626     ASSERT(b.size() <= 1 || b.words()[b.size() - 1]);
    627     ASSERT(S.size() <= 1 || S.words()[S.size() - 1]);
    628 
    629     n = S.size();
    630     ASSERT_WITH_MESSAGE(b.size() <= n, "oversize b in quorem");
    631     if (b.size() < n)
    632         return 0;
    633     sx = S.words();
    634     sxe = sx + --n;
    635     bx = b.words();
    636     bxe = bx + n;
    637     q = *bxe / (*sxe + 1);    /* ensure q <= true quotient */
    638     ASSERT_WITH_MESSAGE(q <= 9, "oversized quotient in quorem");
    639     if (q) {
    640         borrow = 0;
    641         carry = 0;
    642         do {
    643 #ifdef USE_LONG_LONG
    644             ys = *sx++ * (unsigned long long)q + carry;
    645             carry = ys >> 32;
    646             y = *bx - (ys & 0xffffffffUL) - borrow;
    647             borrow = y >> 32 & (uint32_t)1;
    648             *bx++ = (uint32_t)y & 0xffffffffUL;
    649 #else
    650             si = *sx++;
    651             ys = (si & 0xffff) * q + carry;
    652             zs = (si >> 16) * q + (ys >> 16);
    653             carry = zs >> 16;
    654             y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
    655             borrow = (y & 0x10000) >> 16;
    656             z = (*bx >> 16) - (zs & 0xffff) - borrow;
    657             borrow = (z & 0x10000) >> 16;
    658             bx = storeInc(bx, z, y);
    659 #endif
    660         } while (sx <= sxe);
    661         if (!*bxe) {
    662             bx = b.words();
    663             while (--bxe > bx && !*bxe)
    664                 --n;
    665             b.resize(n);
    666         }
    667     }
    668     if (cmp(b, S) >= 0) {
    669         q++;
    670         borrow = 0;
    671         carry = 0;
    672         bx = b.words();
    673         sx = S.words();
    674         do {
    675 #ifdef USE_LONG_LONG
    676             ys = *sx++ + carry;
    677             carry = ys >> 32;
    678             y = *bx - (ys & 0xffffffffUL) - borrow;
    679             borrow = y >> 32 & (uint32_t)1;
    680             *bx++ = (uint32_t)y & 0xffffffffUL;
    681 #else
    682             si = *sx++;
    683             ys = (si & 0xffff) + carry;
    684             zs = (si >> 16) + (ys >> 16);
    685             carry = zs >> 16;
    686             y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
    687             borrow = (y & 0x10000) >> 16;
    688             z = (*bx >> 16) - (zs & 0xffff) - borrow;
    689             borrow = (z & 0x10000) >> 16;
    690             bx = storeInc(bx, z, y);
    691 #endif
    692         } while (sx <= sxe);
    693         bx = b.words();
    694         bxe = bx + n;
    695         if (!*bxe) {
    696             while (--bxe > bx && !*bxe)
    697                 --n;
    698             b.resize(n);
    699         }
    700     }
    701     return q;
    702 }
    703 
    704 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
    705  *
    706  * Inspired by "How to Print Floating-Point Numbers Accurately" by
    707  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
    708  *
    709  * Modifications:
    710  *    1. Rather than iterating, we use a simple numeric overestimate
    711  *       to determine k = floor(log10(d)).  We scale relevant
    712  *       quantities using O(log2(k)) rather than O(k) multiplications.
    713  *    2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
    714  *       try to generate digits strictly left to right.  Instead, we
    715  *       compute with fewer bits and propagate the carry if necessary
    716  *       when rounding the final digit up.  This is often faster.
    717  *    3. Under the assumption that input will be rounded nearest,
    718  *       mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
    719  *       That is, we allow equality in stopping tests when the
    720  *       round-nearest rule will give the same floating-point value
    721  *       as would satisfaction of the stopping test with strict
    722  *       inequality.
    723  *    4. We remove common factors of powers of 2 from relevant
    724  *       quantities.
    725  *    5. When converting floating-point integers less than 1e16,
    726  *       we use floating-point arithmetic rather than resorting
    727  *       to multiple-precision integers.
    728  *    6. When asked to produce fewer than 15 digits, we first try
    729  *       to get by with floating-point arithmetic; we resort to
    730  *       multiple-precision integer arithmetic only if we cannot
    731  *       guarantee that the floating-point calculation has given
    732  *       the correctly rounded result.  For k requested digits and
    733  *       "uniformly" distributed input, the probability is
    734  *       something like 10^(k-15) that we must resort to the int32_t
    735  *       calculation.
    736  *
    737  * Note: 'leftright' translates to 'generate shortest possible string'.
    738  */
    739 template<bool roundingNone, bool roundingSignificantFigures, bool roundingDecimalPlaces, bool leftright>
    740 void dtoa(DtoaBuffer result, double dd, int ndigits, bool& signOut, int& exponentOut, unsigned& precisionOut)
    741 {
    742     // Exactly one rounding mode must be specified.
    743     ASSERT(roundingNone + roundingSignificantFigures + roundingDecimalPlaces == 1);
    744     // roundingNone only allowed (only sensible?) with leftright set.
    745     ASSERT(!roundingNone || leftright);
    746 
    747     ASSERT(std::isfinite(dd));
    748 
    749     int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0,
    750         j, j1, k, k0, k_check, m2, m5, s2, s5,
    751         spec_case;
    752     int32_t L;
    753     int denorm;
    754     uint32_t x;
    755     BigInt b, delta, mlo, mhi, S;
    756     U d2, eps, u;
    757     double ds;
    758     char* s;
    759     char* s0;
    760 
    761     u.d = dd;
    762 
    763     /* Infinity or NaN */
    764     ASSERT((word0(&u) & Exp_mask) != Exp_mask);
    765 
    766     // JavaScript toString conversion treats -0 as 0.
    767     if (!dval(&u)) {
    768         signOut = false;
    769         exponentOut = 0;
    770         precisionOut = 1;
    771         result[0] = '0';
    772         result[1] = '\0';
    773         return;
    774     }
    775 
    776     if (word0(&u) & Sign_bit) {
    777         signOut = true;
    778         word0(&u) &= ~Sign_bit; // clear sign bit
    779     } else
    780         signOut = false;
    781 
    782     d2b(b, &u, &be, &bbits);
    783     if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask >> Exp_shift1)))) {
    784         dval(&d2) = dval(&u);
    785         word0(&d2) &= Frac_mask1;
    786         word0(&d2) |= Exp_11;
    787 
    788         /* log(x)    ~=~ log(1.5) + (x-1.5)/1.5
    789          * log10(x)     =  log(x) / log(10)
    790          *        ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
    791          * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
    792          *
    793          * This suggests computing an approximation k to log10(d) by
    794          *
    795          * k = (i - Bias)*0.301029995663981
    796          *    + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
    797          *
    798          * We want k to be too large rather than too small.
    799          * The error in the first-order Taylor series approximation
    800          * is in our favor, so we just round up the constant enough
    801          * to compensate for any error in the multiplication of
    802          * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
    803          * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
    804          * adding 1e-13 to the constant term more than suffices.
    805          * Hence we adjust the constant term to 0.1760912590558.
    806          * (We could get a more accurate k by invoking log10,
    807          *  but this is probably not worthwhile.)
    808          */
    809 
    810         i -= Bias;
    811         denorm = 0;
    812     } else {
    813         /* d is denormalized */
    814 
    815         i = bbits + be + (Bias + (P - 1) - 1);
    816         x = (i > 32) ? (word0(&u) << (64 - i)) | (word1(&u) >> (i - 32))
    817                 : word1(&u) << (32 - i);
    818         dval(&d2) = x;
    819         word0(&d2) -= 31 * Exp_msk1; /* adjust exponent */
    820         i -= (Bias + (P - 1) - 1) + 1;
    821         denorm = 1;
    822     }
    823     ds = (dval(&d2) - 1.5) * 0.289529654602168 + 0.1760912590558 + (i * 0.301029995663981);
    824     k = (int)ds;
    825     if (ds < 0. && ds != k)
    826         k--;    /* want k = floor(ds) */
    827     k_check = 1;
    828     if (k >= 0 && k <= Ten_pmax) {
    829         if (dval(&u) < tens[k])
    830             k--;
    831         k_check = 0;
    832     }
    833     j = bbits - i - 1;
    834     if (j >= 0) {
    835         b2 = 0;
    836         s2 = j;
    837     } else {
    838         b2 = -j;
    839         s2 = 0;
    840     }
    841     if (k >= 0) {
    842         b5 = 0;
    843         s5 = k;
    844         s2 += k;
    845     } else {
    846         b2 -= k;
    847         b5 = -k;
    848         s5 = 0;
    849     }
    850 
    851     if (roundingNone) {
    852         ilim = ilim1 = -1;
    853         i = 18;
    854         ndigits = 0;
    855     }
    856     if (roundingSignificantFigures) {
    857         if (ndigits <= 0)
    858             ndigits = 1;
    859         ilim = ilim1 = i = ndigits;
    860     }
    861     if (roundingDecimalPlaces) {
    862         i = ndigits + k + 1;
    863         ilim = i;
    864         ilim1 = i - 1;
    865         if (i <= 0)
    866             i = 1;
    867     }
    868 
    869     s = s0 = result;
    870 
    871     if (ilim >= 0 && ilim <= Quick_max) {
    872         /* Try to get by with floating-point arithmetic. */
    873 
    874         i = 0;
    875         dval(&d2) = dval(&u);
    876         k0 = k;
    877         ilim0 = ilim;
    878         ieps = 2; /* conservative */
    879         if (k > 0) {
    880             ds = tens[k & 0xf];
    881             j = k >> 4;
    882             if (j & Bletch) {
    883                 /* prevent overflows */
    884                 j &= Bletch - 1;
    885                 dval(&u) /= bigtens[n_bigtens - 1];
    886                 ieps++;
    887             }
    888             for (; j; j >>= 1, i++) {
    889                 if (j & 1) {
    890                     ieps++;
    891                     ds *= bigtens[i];
    892                 }
    893             }
    894             dval(&u) /= ds;
    895         } else if ((j1 = -k)) {
    896             dval(&u) *= tens[j1 & 0xf];
    897             for (j = j1 >> 4; j; j >>= 1, i++) {
    898                 if (j & 1) {
    899                     ieps++;
    900                     dval(&u) *= bigtens[i];
    901                 }
    902             }
    903         }
    904         if (k_check && dval(&u) < 1. && ilim > 0) {
    905             if (ilim1 <= 0)
    906                 goto fastFailed;
    907             ilim = ilim1;
    908             k--;
    909             dval(&u) *= 10.;
    910             ieps++;
    911         }
    912         dval(&eps) = (ieps * dval(&u)) + 7.;
    913         word0(&eps) -= (P - 1) * Exp_msk1;
    914         if (!ilim) {
    915             S.clear();
    916             mhi.clear();
    917             dval(&u) -= 5.;
    918             if (dval(&u) > dval(&eps))
    919                 goto oneDigit;
    920             if (dval(&u) < -dval(&eps))
    921                 goto noDigits;
    922             goto fastFailed;
    923         }
    924         if (leftright) {
    925             /* Use Steele & White method of only
    926              * generating digits needed.
    927              */
    928             dval(&eps) = (0.5 / tens[ilim - 1]) - dval(&eps);
    929             for (i = 0;;) {
    930                 L = (long int)dval(&u);
    931                 dval(&u) -= L;
    932                 *s++ = '0' + (int)L;
    933                 if (dval(&u) < dval(&eps))
    934                     goto ret;
    935                 if (1. - dval(&u) < dval(&eps))
    936                     goto bumpUp;
    937                 if (++i >= ilim)
    938                     break;
    939                 dval(&eps) *= 10.;
    940                 dval(&u) *= 10.;
    941             }
    942         } else {
    943             /* Generate ilim digits, then fix them up. */
    944             dval(&eps) *= tens[ilim - 1];
    945             for (i = 1;; i++, dval(&u) *= 10.) {
    946                 L = (int32_t)(dval(&u));
    947                 if (!(dval(&u) -= L))
    948                     ilim = i;
    949                 *s++ = '0' + (int)L;
    950                 if (i == ilim) {
    951                     if (dval(&u) > 0.5 + dval(&eps))
    952                         goto bumpUp;
    953                     if (dval(&u) < 0.5 - dval(&eps)) {
    954                         while (*--s == '0') { }
    955                         s++;
    956                         goto ret;
    957                     }
    958                     break;
    959                 }
    960             }
    961         }
    962 fastFailed:
    963         s = s0;
    964         dval(&u) = dval(&d2);
    965         k = k0;
    966         ilim = ilim0;
    967     }
    968 
    969     /* Do we have a "small" integer? */
    970 
    971     if (be >= 0 && k <= Int_max) {
    972         /* Yes. */
    973         ds = tens[k];
    974         if (ndigits < 0 && ilim <= 0) {
    975             S.clear();
    976             mhi.clear();
    977             if (ilim < 0 || dval(&u) <= 5 * ds)
    978                 goto noDigits;
    979             goto oneDigit;
    980         }
    981         for (i = 1;; i++, dval(&u) *= 10.) {
    982             L = (int32_t)(dval(&u) / ds);
    983             dval(&u) -= L * ds;
    984             *s++ = '0' + (int)L;
    985             if (!dval(&u)) {
    986                 break;
    987             }
    988             if (i == ilim) {
    989                 dval(&u) += dval(&u);
    990                 if (dval(&u) > ds || (dval(&u) == ds && (L & 1))) {
    991 bumpUp:
    992                     while (*--s == '9')
    993                         if (s == s0) {
    994                             k++;
    995                             *s = '0';
    996                             break;
    997                         }
    998                     ++*s++;
    999                 }
   1000                 break;
   1001             }
   1002         }
   1003         goto ret;
   1004     }
   1005 
   1006     m2 = b2;
   1007     m5 = b5;
   1008     mhi.clear();
   1009     mlo.clear();
   1010     if (leftright) {
   1011         i = denorm ? be + (Bias + (P - 1) - 1 + 1) : 1 + P - bbits;
   1012         b2 += i;
   1013         s2 += i;
   1014         i2b(mhi, 1);
   1015     }
   1016     if (m2 > 0 && s2 > 0) {
   1017         i = m2 < s2 ? m2 : s2;
   1018         b2 -= i;
   1019         m2 -= i;
   1020         s2 -= i;
   1021     }
   1022     if (b5 > 0) {
   1023         if (leftright) {
   1024             if (m5 > 0) {
   1025                 pow5mult(mhi, m5);
   1026                 mult(b, mhi);
   1027             }
   1028             if ((j = b5 - m5))
   1029                 pow5mult(b, j);
   1030         } else
   1031             pow5mult(b, b5);
   1032     }
   1033     i2b(S, 1);
   1034     if (s5 > 0)
   1035         pow5mult(S, s5);
   1036 
   1037     /* Check for special case that d is a normalized power of 2. */
   1038 
   1039     spec_case = 0;
   1040     if ((roundingNone || leftright) && (!word1(&u) && !(word0(&u) & Bndry_mask) && word0(&u) & (Exp_mask & ~Exp_msk1))) {
   1041         /* The special case */
   1042         b2 += Log2P;
   1043         s2 += Log2P;
   1044         spec_case = 1;
   1045     }
   1046 
   1047     /* Arrange for convenient computation of quotients:
   1048      * shift left if necessary so divisor has 4 leading 0 bits.
   1049      *
   1050      * Perhaps we should just compute leading 28 bits of S once
   1051      * and for all and pass them and a shift to quorem, so it
   1052      * can do shifts and ors to compute the numerator for q.
   1053      */
   1054     if ((i = ((s5 ? 32 - hi0bits(S.words()[S.size() - 1]) : 1) + s2) & 0x1f))
   1055         i = 32 - i;
   1056     if (i > 4) {
   1057         i -= 4;
   1058         b2 += i;
   1059         m2 += i;
   1060         s2 += i;
   1061     } else if (i < 4) {
   1062         i += 28;
   1063         b2 += i;
   1064         m2 += i;
   1065         s2 += i;
   1066     }
   1067     if (b2 > 0)
   1068         lshift(b, b2);
   1069     if (s2 > 0)
   1070         lshift(S, s2);
   1071     if (k_check) {
   1072         if (cmp(b, S) < 0) {
   1073             k--;
   1074             multadd(b, 10, 0);    /* we botched the k estimate */
   1075             if (leftright)
   1076                 multadd(mhi, 10, 0);
   1077             ilim = ilim1;
   1078         }
   1079     }
   1080     if (ilim <= 0 && roundingDecimalPlaces) {
   1081         if (ilim < 0)
   1082             goto noDigits;
   1083         multadd(S, 5, 0);
   1084         // For IEEE-754 unbiased rounding this check should be <=, such that 0.5 would flush to zero.
   1085         if (cmp(b, S) < 0)
   1086             goto noDigits;
   1087         goto oneDigit;
   1088     }
   1089     if (leftright) {
   1090         if (m2 > 0)
   1091             lshift(mhi, m2);
   1092 
   1093         /* Compute mlo -- check for special case
   1094          * that d is a normalized power of 2.
   1095          */
   1096 
   1097         mlo = mhi;
   1098         if (spec_case)
   1099             lshift(mhi, Log2P);
   1100 
   1101         for (i = 1;;i++) {
   1102             dig = quorem(b, S) + '0';
   1103             /* Do we yet have the shortest decimal string
   1104              * that will round to d?
   1105              */
   1106             j = cmp(b, mlo);
   1107             diff(delta, S, mhi);
   1108             j1 = delta.sign ? 1 : cmp(b, delta);
   1109 #ifdef DTOA_ROUND_BIASED
   1110             if (j < 0 || !j) {
   1111 #else
   1112             // FIXME: ECMA-262 specifies that equidistant results round away from
   1113             // zero, which probably means we shouldn't be on the unbiased code path
   1114             // (the (word1(&u) & 1) clause is looking highly suspicious). I haven't
   1115             // yet understood this code well enough to make the call, but we should
   1116             // probably be enabling DTOA_ROUND_BIASED. I think the interesting corner
   1117             // case to understand is probably "Math.pow(0.5, 24).toString()".
   1118             // I believe this value is interesting because I think it is precisely
   1119             // representable in binary floating point, and its decimal representation
   1120             // has a single digit that Steele & White reduction can remove, with the
   1121             // value 5 (thus equidistant from the next numbers above and below).
   1122             // We produce the correct answer using either codepath, and I don't as
   1123             // yet understand why. :-)
   1124             if (!j1 && !(word1(&u) & 1)) {
   1125                 if (dig == '9')
   1126                     goto round9up;
   1127                 if (j > 0)
   1128                     dig++;
   1129                 *s++ = dig;
   1130                 goto ret;
   1131             }
   1132             if (j < 0 || (!j && !(word1(&u) & 1))) {
   1133 #endif
   1134                 if ((b.words()[0] || b.size() > 1) && (j1 > 0)) {
   1135                     lshift(b, 1);
   1136                     j1 = cmp(b, S);
   1137                     // For IEEE-754 round-to-even, this check should be (j1 > 0 || (!j1 && (dig & 1))),
   1138                     // but ECMA-262 specifies that equidistant values (e.g. (.5).toFixed()) should
   1139                     // be rounded away from zero.
   1140                     if (j1 >= 0) {
   1141                         if (dig == '9')
   1142                             goto round9up;
   1143                         dig++;
   1144                     }
   1145                 }
   1146                 *s++ = dig;
   1147                 goto ret;
   1148             }
   1149             if (j1 > 0) {
   1150                 if (dig == '9') { /* possible if i == 1 */
   1151 round9up:
   1152                     *s++ = '9';
   1153                     goto roundoff;
   1154                 }
   1155                 *s++ = dig + 1;
   1156                 goto ret;
   1157             }
   1158             *s++ = dig;
   1159             if (i == ilim)
   1160                 break;
   1161             multadd(b, 10, 0);
   1162             multadd(mlo, 10, 0);
   1163             multadd(mhi, 10, 0);
   1164         }
   1165     } else {
   1166         for (i = 1;; i++) {
   1167             *s++ = dig = quorem(b, S) + '0';
   1168             if (!b.words()[0] && b.size() <= 1)
   1169                 goto ret;
   1170             if (i >= ilim)
   1171                 break;
   1172             multadd(b, 10, 0);
   1173         }
   1174     }
   1175 
   1176     /* Round off last digit */
   1177 
   1178     lshift(b, 1);
   1179     j = cmp(b, S);
   1180     // For IEEE-754 round-to-even, this check should be (j > 0 || (!j && (dig & 1))),
   1181     // but ECMA-262 specifies that equidistant values (e.g. (.5).toFixed()) should
   1182     // be rounded away from zero.
   1183     if (j >= 0) {
   1184 roundoff:
   1185         while (*--s == '9')
   1186             if (s == s0) {
   1187                 k++;
   1188                 *s++ = '1';
   1189                 goto ret;
   1190             }
   1191         ++*s++;
   1192     } else {
   1193         while (*--s == '0') { }
   1194         s++;
   1195     }
   1196     goto ret;
   1197 noDigits:
   1198     exponentOut = 0;
   1199     precisionOut = 1;
   1200     result[0] = '0';
   1201     result[1] = '\0';
   1202     return;
   1203 oneDigit:
   1204     *s++ = '1';
   1205     k++;
   1206     goto ret;
   1207 ret:
   1208     ASSERT(s > result);
   1209     *s = 0;
   1210     exponentOut = k;
   1211     precisionOut = s - result;
   1212 }
   1213 
   1214 void dtoa(DtoaBuffer result, double dd, bool& sign, int& exponent, unsigned& precision)
   1215 {
   1216     // flags are roundingNone, leftright.
   1217     dtoa<true, false, false, true>(result, dd, 0, sign, exponent, precision);
   1218 }
   1219 
   1220 void dtoaRoundSF(DtoaBuffer result, double dd, int ndigits, bool& sign, int& exponent, unsigned& precision)
   1221 {
   1222     // flag is roundingSignificantFigures.
   1223     dtoa<false, true, false, false>(result, dd, ndigits, sign, exponent, precision);
   1224 }
   1225 
   1226 void dtoaRoundDP(DtoaBuffer result, double dd, int ndigits, bool& sign, int& exponent, unsigned& precision)
   1227 {
   1228     // flag is roundingDecimalPlaces.
   1229     dtoa<false, false, true, false>(result, dd, ndigits, sign, exponent, precision);
   1230 }
   1231 
   1232 const char* numberToString(double d, NumberToStringBuffer buffer)
   1233 {
   1234     double_conversion::StringBuilder builder(buffer, NumberToStringBufferLength);
   1235     const double_conversion::DoubleToStringConverter& converter = double_conversion::DoubleToStringConverter::EcmaScriptConverter();
   1236     converter.ToShortest(d, &builder);
   1237     return builder.Finalize();
   1238 }
   1239 
   1240 static inline const char* formatStringTruncatingTrailingZerosIfNeeded(NumberToStringBuffer buffer, double_conversion::StringBuilder& builder)
   1241 {
   1242     size_t length = builder.position();
   1243     size_t decimalPointPosition = 0;
   1244     for (; decimalPointPosition < length; ++decimalPointPosition) {
   1245         if (buffer[decimalPointPosition] == '.')
   1246             break;
   1247     }
   1248 
   1249     // No decimal seperator found, early exit.
   1250     if (decimalPointPosition == length)
   1251         return builder.Finalize();
   1252 
   1253     size_t truncatedLength = length - 1;
   1254     for (; truncatedLength > decimalPointPosition; --truncatedLength) {
   1255         if (buffer[truncatedLength] != '0')
   1256             break;
   1257     }
   1258 
   1259     // No trailing zeros found to strip.
   1260     if (truncatedLength == length - 1)
   1261         return builder.Finalize();
   1262 
   1263     // If we removed all trailing zeros, remove the decimal point as well.
   1264     if (truncatedLength == decimalPointPosition) {
   1265         ASSERT(truncatedLength > 0);
   1266         --truncatedLength;
   1267     }
   1268 
   1269     // Truncate the StringBuilder, and return the final result.
   1270     builder.SetPosition(truncatedLength + 1);
   1271     return builder.Finalize();
   1272 }
   1273 
   1274 const char* numberToFixedPrecisionString(double d, unsigned significantFigures, NumberToStringBuffer buffer, bool truncateTrailingZeros)
   1275 {
   1276     // Mimic String::format("%.[precision]g", ...), but use dtoas rounding facilities.
   1277     // "g": Signed value printed in f or e format, whichever is more compact for the given value and precision.
   1278     // The e format is used only when the exponent of the value is less than 4 or greater than or equal to the
   1279     // precision argument. Trailing zeros are truncated, and the decimal point appears only if one or more digits follow it.
   1280     // "precision": The precision specifies the maximum number of significant digits printed.
   1281     double_conversion::StringBuilder builder(buffer, NumberToStringBufferLength);
   1282     const double_conversion::DoubleToStringConverter& converter = double_conversion::DoubleToStringConverter::EcmaScriptConverter();
   1283     converter.ToPrecision(d, significantFigures, &builder);
   1284     if (!truncateTrailingZeros)
   1285         return builder.Finalize();
   1286     return formatStringTruncatingTrailingZerosIfNeeded(buffer, builder);
   1287 }
   1288 
   1289 const char* numberToFixedWidthString(double d, unsigned decimalPlaces, NumberToStringBuffer buffer)
   1290 {
   1291     // Mimic String::format("%.[precision]f", ...), but use dtoas rounding facilities.
   1292     // "f": Signed value having the form [  ]dddd.dddd, where dddd is one or more decimal digits.
   1293     // The number of digits before the decimal point depends on the magnitude of the number, and
   1294     // the number of digits after the decimal point depends on the requested precision.
   1295     // "precision": The precision value specifies the number of digits after the decimal point.
   1296     // If a decimal point appears, at least one digit appears before it.
   1297     // The value is rounded to the appropriate number of digits.
   1298     double_conversion::StringBuilder builder(buffer, NumberToStringBufferLength);
   1299     const double_conversion::DoubleToStringConverter& converter = double_conversion::DoubleToStringConverter::EcmaScriptConverter();
   1300     converter.ToFixed(d, decimalPlaces, &builder);
   1301     return builder.Finalize();
   1302 }
   1303 
   1304 namespace Internal {
   1305 
   1306 double parseDoubleFromLongString(const UChar* string, size_t length, size_t& parsedLength)
   1307 {
   1308     Vector<LChar> conversionBuffer(length);
   1309     for (size_t i = 0; i < length; ++i)
   1310         conversionBuffer[i] = isASCII(string[i]) ? string[i] : 0;
   1311     return parseDouble(conversionBuffer.data(), length, parsedLength);
   1312 }
   1313 
   1314 } // namespace Internal
   1315 
   1316 } // namespace WTF
   1317