/****************************************************************

 *

 * The author of this software is David M. Gay.

 *

 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.

 * Copyright (C) 2002, 2005, 2006, 2007, 2008 Apple Inc. All rights reserved.

 *

 * Permission to use, copy, modify, and distribute this software for any

 * purpose without fee is hereby granted, provided that this entire notice

 * is included in all copies of any software which is or includes a copy

 * or modification of this software and in all copies of the supporting

 * documentation for such software.

 *

 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED

 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY

 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY

 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.

 *

 ***************************************************************/

/* Please send bug reports to

    David M. Gay

    Bell Laboratories, Room 2C-463

    600 Mountain Avenue

    Murray Hill, NJ 07974-0636

    U.S.A.

    dmg@bell-labs.com

 */

/* On a machine with IEEE extended-precision registers, it is

 * necessary to specify double-precision (53-bit) rounding precision

 * before invoking strtod or dtoa.  If the machine uses (the equivalent

 * of) Intel 80x87 arithmetic, the call

 *    _control87(PC_53, MCW_PC);

 * does this with many compilers.  Whether this or another call is

 * appropriate depends on the compiler; for this to work, it may be

 * necessary to #include "float.h" or another system-dependent header

 * file.

 */

/* strtod for IEEE-arithmetic machines.

 *

 * This strtod returns a nearest machine number to the input decimal

 * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are

 * broken by the IEEE round-even rule.  Otherwise ties are broken by

 * biased rounding (add half and chop).

 *

 * Inspired loosely by William D. Clinger's paper "How to Read Floating

 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].

 *

 * Modifications:

 *

 *    1. We only require IEEE.

 *    2. We get by with floating-point arithmetic in a case that

 *        Clinger missed -- when we're computing d * 10^n

 *        for a small integer d and the integer n is not too

 *        much larger than 22 (the maximum integer k for which

 *        we can represent 10^k exactly), we may be able to

 *        compute (d*10^k) * 10^(e-k) with just one roundoff.

 *    3. Rather than a bit-at-a-time adjustment of the binary

 *        result in the hard case, we use floating-point

 *        arithmetic to determine the adjustment to within

 *        one bit; only in really hard cases do we need to

 *        compute a second residual.

 *    4. Because of 3., we don't need a large table of powers of 10

 *        for ten-to-e (just some small tables, e.g. of 10^k

 *        for 0 <= k <= 22).

 */

/*

 * #define IEEE_8087 for IEEE-arithmetic machines where the least

 *    significant byte has the lowest address.

 * #define IEEE_MC68k for IEEE-arithmetic machines where the most

 *    significant byte has the lowest address.

 * #define No_leftright to omit left-right logic in fast floating-point

 *    computation of dtoa.

 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3

 *    and Honor_FLT_ROUNDS is not #defined.

 * #define Inaccurate_Divide for IEEE-format with correctly rounded

 *    products but inaccurate quotients, e.g., for Intel i860.

 * #define USE_LONG_LONG on machines that have a "long long"

 *    integer type (of >= 64 bits), and performance testing shows that

 *    it is faster than 32-bit fallback (which is often not the case

 *    on 32-bit machines). On such machines, you can #define Just_16

 *    to store 16 bits per 32-bit int32_t when doing high-precision integer

 *    arithmetic.  Whether this speeds things up or slows things down

 *    depends on the machine and the number being converted.

 * #define Bad_float_h if your system lacks a float.h or if it does not

 *    define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,

 *    FLT_RADIX, FLT_ROUNDS, and DBL_MAX.

 * #define INFNAN_CHECK on IEEE systems to cause strtod to check for

 *    Infinity and NaN (case insensitively).  On some systems (e.g.,

 *    some HP systems), it may be necessary to #define NAN_WORD0

 *    appropriately -- to the most significant word of a quiet NaN.

 *    (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)

 *    When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,

 *    strtod also accepts (case insensitively) strings of the form

 *    NaN(x), where x is a string of hexadecimal digits and spaces;

 *    if there is only one string of hexadecimal digits, it is taken

 *    for the 52 fraction bits of the resulting NaN; if there are two

 *    or more strings of hex digits, the first is for the high 20 bits,

 *    the second and subsequent for the low 32 bits, with intervening

 *    white space ignored; but if this results in none of the 52

 *    fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0

 *    and NAN_WORD1 are used instead.

 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that

 *    avoids underflows on inputs whose result does not underflow.

 *    If you #define NO_IEEE_Scale on a machine that uses IEEE-format

 *    floating-point numbers and flushes underflows to zero rather

 *    than implementing gradual underflow, then you must also #define

 *    Sudden_Underflow.

 * #define YES_ALIAS to permit aliasing certain double values with

 *    arrays of ULongs.  This leads to slightly better code with

 *    some compilers and was always used prior to 19990916, but it

 *    is not strictly legal and can cause trouble with aggressively

 *    optimizing compilers (e.g., gcc 2.95.1 under -O2).

 * #define SET_INEXACT if IEEE arithmetic is being used and extra

 *    computation should be done to set the inexact flag when the

 *    result is inexact and avoid setting inexact when the result

 *    is exact.  In this case, dtoa.c must be compiled in

 *    an environment, perhaps provided by #include "dtoa.c" in a

 *    suitable wrapper, that defines two functions,

 *        int get_inexact(void);

 *        void clear_inexact(void);

 *    such that get_inexact() returns a nonzero value if the

 *    inexact bit is already set, and clear_inexact() sets the

 *    inexact bit to 0.  When SET_INEXACT is #defined, strtod

 *    also does extra computations to set the underflow and overflow

 *    flags when appropriate (i.e., when the result is tiny and

 *    inexact or when it is a numeric value rounded to +-infinity).

 * #define NO_ERRNO if strtod should not assign errno = ERANGE when

 *    the result overflows to +-Infinity or underflows to 0.

 */

#include "config.h"

#include "dtoa.h"

#if HAVE(ERRNO_H)

#include <errno.h>

#else

#define NO_ERRNO

#endif

#include <math.h>

#include <stdint.h>

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <wtf/AlwaysInline.h>

#include <wtf/Assertions.h>

#include <wtf/FastMalloc.h>

#include <wtf/MathExtras.h>

#include <wtf/Threading.h>

#include <wtf/Vector.h>

#if COMPILER(MSVC)

#pragma warning(disable: 4244)

#pragma warning(disable: 4245)

#pragma warning(disable: 4554)

#endif

#if CPU(BIG_ENDIAN)

#define IEEE_MC68k

#elif CPU(MIDDLE_ENDIAN)

#define IEEE_ARM

#else

#define IEEE_8087

#endif

#define INFNAN_CHECK

#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(IEEE_ARM) != 1

Exactly one of IEEE_8087, IEEE_ARM or IEEE_MC68k should be defined.

#endif

namespace WTF {

#if ENABLE(JSC_MULTIPLE_THREADS)

Mutex* s_dtoaP5Mutex;

#endif

typedef union {

    double d;

    uint32_t L[2];

} U;

#ifdef YES_ALIAS

#define dval(x) x

#ifdef IEEE_8087

#define word0(x) ((uint32_t*)&x)[1]

#define word1(x) ((uint32_t*)&x)[0]

#else

#define word0(x) ((uint32_t*)&x)[0]

#define word1(x) ((uint32_t*)&x)[1]

#endif

#else

#ifdef IEEE_8087

#define word0(x) (x)->L[1]

#define word1(x) (x)->L[0]

#else

#define word0(x) (x)->L[0]

#define word1(x) (x)->L[1]

#endif

#define dval(x) (x)->d

#endif

/* The following definition of Storeinc is appropriate for MIPS processors.

 * An alternative that might be better on some machines is

 *  *p++ = high << 16 | low & 0xffff;

 */

static ALWAYS_INLINE uint32_t* storeInc(uint32_t* p, uint16_t high, uint16_t low)

{

    uint16_t* p16 = reinterpret_cast<uint16_t*>(p);

#if defined(IEEE_8087) || defined(IEEE_ARM)

    p16[1] = high;

    p16[0] = low;

#else

    p16[0] = high;

    p16[1] = low;

#endif

    return p + 1;

}

#define Exp_shift  20

#define Exp_shift1 20

#define Exp_msk1    0x100000

#define Exp_msk11   0x100000

#define Exp_mask  0x7ff00000

#define P 53

#define Bias 1023

#define Emin (-1022)

#define Exp_1  0x3ff00000

#define Exp_11 0x3ff00000

#define Ebits 11

#define Frac_mask  0xfffff

#define Frac_mask1 0xfffff

#define Ten_pmax 22

#define Bletch 0x10

#define Bndry_mask  0xfffff

#define Bndry_mask1 0xfffff

#define LSB 1

#define Sign_bit 0x80000000

#define Log2P 1

#define Tiny0 0

#define Tiny1 1

#define Quick_max 14

#define Int_max 14

#if !defined(NO_IEEE_Scale)

#undef Avoid_Underflow

#define Avoid_Underflow

#endif

#if !defined(Flt_Rounds)

#if defined(FLT_ROUNDS)

#define Flt_Rounds FLT_ROUNDS

#else

#define Flt_Rounds 1

#endif

#endif /* Flt_Rounds */

#define rounded_product(a, b) a *= b

#define rounded_quotient(a, b) a /= b

#define Big0 (Frac_mask1 | Exp_msk1 * (DBL_MAX_EXP + Bias - 1))

#define Big1 0xffffffff

// FIXME: we should remove non-Pack_32 mode since it is unused and unmaintained

#ifndef Pack_32

#define Pack_32

#endif

#if CPU(PPC64) || CPU(X86_64)

// FIXME: should we enable this on all 64-bit CPUs?

// 64-bit emulation provided by the compiler is likely to be slower than dtoa own code on 32-bit hardware.

#define USE_LONG_LONG

#endif

#ifndef USE_LONG_LONG

#ifdef Just_16

#undef Pack_32

/* When Pack_32 is not defined, we store 16 bits per 32-bit int32_t.

 * This makes some inner loops simpler and sometimes saves work

 * during multiplications, but it often seems to make things slightly

 * slower.  Hence the default is now to store 32 bits per int32_t.

 */

#endif

#endif

#define Kmax 15

struct BigInt {

    BigInt() : sign(0) { }

    int sign;

    void clear()

    {

        sign = 0;

        m_words.clear();

    }

    size_t size() const

    {

        return m_words.size();

    }

    void resize(size_t s)

    {

        m_words.resize(s);

    }

    uint32_t* words()

    {

        return m_words.data();

    }

    const uint32_t* words() const

    {

        return m_words.data();

    }

    void append(uint32_t w)

    {

        m_words.append(w);

    }

    Vector<uint32_t, 16> m_words;

};

static void multadd(BigInt& b, int m, int a)    /* multiply by m and add a */

{

#ifdef USE_LONG_LONG

    unsigned long long carry;

#else

    uint32_t carry;

#endif

    int wds = b.size();

    uint32_t* x = b.words();

    int i = 0;

    carry = a;

    do {

#ifdef USE_LONG_LONG

        unsigned long long y = *x * (unsigned long long)m + carry;

        carry = y >> 32;

        *x++ = (uint32_t)y & 0xffffffffUL;

#else

#ifdef Pack_32

        uint32_t xi = *x;

        uint32_t y = (xi & 0xffff) * m + carry;

        uint32_t z = (xi >> 16) * m + (y >> 16);

        carry = z >> 16;

        *x++ = (z << 16) + (y & 0xffff);

#else

        uint32_t y = *x * m + carry;

        carry = y >> 16;

        *x++ = y & 0xffff;

#endif

#endif

    } while (++i < wds);

    if (carry)

        b.append((uint32_t)carry);

}

static void s2b(BigInt& b, const char* s, int nd0, int nd, uint32_t y9)

{

    int k;

    int32_t y;

    int32_t x = (nd + 8) / 9;

    for (k = 0, y = 1; x > y; y <<= 1, k++) { }

#ifdef Pack_32

    b.sign = 0;

    b.resize(1);

    b.words()[0] = y9;

#else

    b.sign = 0;

    b.resize((b->x[1] = y9 >> 16) ? 2 : 1);

    b.words()[0] = y9 & 0xffff;

#endif

    int i = 9;

    if (9 < nd0) {

        s += 9;

        do {

            multadd(b, 10, *s++ - '0');

        } while (++i < nd0);

        s++;

    } else

        s += 10;

    for (; i < nd; i++)

        multadd(b, 10, *s++ - '0');

}

static int hi0bits(uint32_t x)

{

    int k = 0;

    if (!(x & 0xffff0000)) {

        k = 16;

        x <<= 16;

    }

    if (!(x & 0xff000000)) {

        k += 8;

        x <<= 8;

    }

    if (!(x & 0xf0000000)) {

        k += 4;

        x <<= 4;

    }

    if (!(x & 0xc0000000)) {

        k += 2;

        x <<= 2;

    }

    if (!(x & 0x80000000)) {

        k++;

        if (!(x & 0x40000000))

            return 32;

    }

    return k;

}

static int lo0bits(uint32_t* y)

{

    int k;

    uint32_t x = *y;

    if (x & 7) {

        if (x & 1)

            return 0;

        if (x & 2) {

            *y = x >> 1;

            return 1;

        }

        *y = x >> 2;

        return 2;

    }

    k = 0;

    if (!(x & 0xffff)) {

        k = 16;

        x >>= 16;

    }

    if (!(x & 0xff)) {

        k += 8;

        x >>= 8;

    }

    if (!(x & 0xf)) {

        k += 4;

        x >>= 4;

    }

    if (!(x & 0x3)) {

        k += 2;

        x >>= 2;

    }

    if (!(x & 1)) {

        k++;

        x >>= 1;

        if (!x & 1)

            return 32;

    }

    *y = x;

    return k;

}

static void i2b(BigInt& b, int i)

{

    b.sign = 0;

    b.resize(1);

    b.words()[0] = i;

}

static void mult(BigInt& aRef, const BigInt& bRef)

{

    const BigInt* a = &aRef;

    const BigInt* b = &bRef;

    BigInt c;

    int wa, wb, wc;

    const uint32_t* x = 0;

    const uint32_t* xa;

    const uint32_t* xb;

    const uint32_t* xae;

    const uint32_t* xbe;

    uint32_t* xc;

    uint32_t* xc0;

    uint32_t y;

#ifdef USE_LONG_LONG

    unsigned long long carry, z;

#else

    uint32_t carry, z;

#endif

    if (a->size() < b->size()) {

        const BigInt* tmp = a;

        a = b;

        b = tmp;

    }

    wa = a->size();

    wb = b->size();

    wc = wa + wb;

    c.resize(wc);

    for (xc = c.words(), xa = xc + wc; xc < xa; xc++)

        *xc = 0;

    xa = a->words();

    xae = xa + wa;

    xb = b->words();

    xbe = xb + wb;

    xc0 = c.words();

#ifdef USE_LONG_LONG

    for (; xb < xbe; xc0++) {

        if ((y = *xb++)) {

            x = xa;

            xc = xc0;

            carry = 0;

            do {

                z = *x++ * (unsigned long long)y + *xc + carry;

                carry = z >> 32;

                *xc++ = (uint32_t)z & 0xffffffffUL;

            } while (x < xae);

            *xc = (uint32_t)carry;

        }

    }

#else

#ifdef Pack_32

    for (; xb < xbe; xb++, xc0++) {