pjs/js2/src/numerics.cpp

/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
*
* The contents of this file are subject to the Netscape Public
* License Version 1.1 (the "License"); you may not use this file
* except in compliance with the License. You may obtain a copy of
* the License at http://www.mozilla.org/NPL/
*
* Software distributed under the License is distributed on an "AS
* IS" basis, WITHOUT WARRANTY OF ANY KIND, either express or
* implied. See the License for the specific language governing
* rights and limitations under the License.
*
* The Original Code is the JavaScript 2 Prototype.
*
* The Initial Developer of the Original Code is Netscape
* Communications Corporation.   Portions created by Netscape are
* Copyright (C) 1998 Netscape Communications Corporation. All
* Rights Reserved.
*
* Contributor(s):
*
* Alternatively, the contents of this file may be used under the
* terms of the GNU Public License (the "GPL"), in which case the
* provisions of the GPL are applicable instead of those above.
* If you wish to allow use of your version of this file only
* under the terms of the GPL and not to allow others to use your
* version of this file under the NPL, indicate your decision by
* deleting the provisions above and replace them with the notice
* and other provisions required by the GPL.  If you do not delete
* the provisions above, a recipient may use your version of this
* file under either the NPL or the GPL.
*/

#include <cstdlib>
#include <cstring>
#include <cfloat>

#include <vector>

#include "systemtypes.h"
#include "js2value.h"

#include "numerics.h"
#include "parser.h"

#ifdef DIKDIK
#include "js2runtime.h"
#endif

#include "fdlibm_ns.h"

namespace JS = JavaScript;
using namespace JS;

#ifdef DIKDIK
using namespace JS::JS2Runtime;
#endif

//
// Portable double-precision floating point to string and back conversions
//

// ****************************************************************
//
// The author of this software is David M. Gay.
//
// Copyright (c) 1991 by Lucent Technologies.
//
// 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 strToDouble or doubleToAscii.  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.


// strToDouble for IEEE-arithmetic machines.
//
// This strToDouble returns a nearest machine number to the input decimal
// string.  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 double-precision
//      arithmetic (not IEEE double-extended).
//  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 Sudden_Underflow for IEEE-format machines without gradual
//  underflow (i.e., that flush to zero on underflow).
// #define No_leftright to omit left-right logic in fast floating-point
//  computation of doubleToAscii.
// #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
// #define ROUND_BIASED for IEEE-format with biased rounding.
// #define Inaccurate_Divide for IEEE-format with correctly rounded
//  products but inaccurate quotients, e.g., for Intel i860.
// #define NATIVE_INT64 on machines that have "long long"
//  64-bit integer types int64 and uint64.
// #define JS_THREADSAFE if the system offers preemptively scheduled
//  multiple threads.  In this case, you must provide (or suitably
//  #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
//  by FREE_DTOA_LOCK(n) for n = 0 or 1.  (The second lock, accessed
//  in pow5Mul, ensures lazy evaluation of only one copy of high
//  powers of 5; omitting this lock would introduce a small
//  probability of wasting memory, but would otherwise be harmless.)
//  You must also invoke freeDToA(s) to free the value s returned by
//  doubleToAscii.  You may do so whether or not JS_THREADSAFE is #defined.
// #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strToDouble that
//  avoids underflows on inputs whose result does not underflow.

#ifdef IS_LITTLE_ENDIAN
# define IEEE_8087
#else
# define IEEE_MC68k
#endif


// Stefan Hanske <sh990154@mail.uni-greifswald.de> reports:
//  ARM is a little endian architecture but 64 bit double words are stored
// differently: the 32 bit words are in little endian byte order, the two words
// are stored in big endian`s way.
#if defined (IEEE_8087) && !defined(__arm)
# define word0(x) ((uint32 *)&x)[1]
# define word1(x) ((uint32 *)&x)[0]
#else
# define word0(x) ((uint32 *)&x)[0]
# define word1(x) ((uint32 *)&x)[1]
#endif

// The following definition of Storeinc is appropriate for MIPS processors.
// An alternative that might be better on some machines is
// #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
#if defined(IEEE_8087)
#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
((unsigned short *)a)[0] = (unsigned short)c, a++)
#else
#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
((unsigned short *)a)[1] = (unsigned short)c, a++)
#endif

// #define P DBL_MANT_DIG
// Ten_pmax = floor(P*log(2)/log(5))
// Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16
// Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1)
// Int_max = floor(P*log(FLT_RADIX)/log(10) - 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
#define Infinite(x) (word0(x) == 0x7ff00000) // sufficient test for here
#ifndef NO_IEEE_Scale
# define Avoid_Underflow
#endif

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

#ifdef JS_THREADSAFE
static PRLock *freelist_lock;
# define ACQUIRE_DTOA_LOCK(n) PR_Lock(freelist_lock)
# define FREE_DTOA_LOCK(n) PR_Unlock(freelist_lock)
#else
# define ACQUIRE_DTOA_LOCK(n)
# define FREE_DTOA_LOCK(n)
#endif


//
// Double-precision constants
//
double JS::positiveInfinity;
double JS::negativeInfinity;
double JS::nan;
double JS::minDouble;
double JS::maxDouble;



struct InitNumerics {InitNumerics();};
static InitNumerics initNumerics;

InitNumerics::InitNumerics()
{
    word0(JS::positiveInfinity) = Exp_mask;
    word1(JS::positiveInfinity) = 0;
    word0(JS::negativeInfinity) = Exp_mask | Sign_bit;
    word1(JS::negativeInfinity) = 0;
    word0(JS::nan) = 0x7FFFFFFF;
    word1(JS::nan) = 0xFFFFFFFF;
    word0(JS::minDouble) = 0;
    word1(JS::minDouble) = 1;
    JS::maxDouble = 1.7976931348623157E+308;
}

// had to move these here since they depend upon the values
// initialized above, and we can't guarantee order other than
// lexically in a single file.
#ifdef DIKDIK
js2val JS::JS2Runtime::kUndefinedValue = JS2VAL_VOID;
js2val JS::JS2Runtime::kNaNValue = JSValue::newNumber(nan);
js2val JS::JS2Runtime::kTrueValue = JSValue::newBoolean(true);
js2val JS::JS2Runtime::kFalseValue = JSValue::newBoolean(false);
js2val JS::JS2Runtime::kNullValue = JS2VAL_NULL;
js2val JS::JS2Runtime::kNegativeZero = JSValue::newNumber(-0.0);
js2val JS::JS2Runtime::kPositiveZero = JSValue::newNumber(0.0);
js2val JS::JS2Runtime::kNegativeInfinity = JSValue::newNumber(negativeInfinity);
js2val JS::JS2Runtime::kPositiveInfinity = JSValue::newNumber(positiveInfinity);
#endif






//
// Portable double-precision floating point to string and back conversions
//


// Return the absolute difference between x and the adjacent greater-magnitude
// double number (ignoring exponent overflows).
double JS::ulp(double x)
{
    int32 L;
    double a;

    L = int32((word0(x) & Exp_mask) - (P-1)*Exp_msk1);
#ifndef Sudden_Underflow
    if (L > 0) {
#endif
        word0(a) = uint32(L);
        word1(a) = 0;
#ifndef Sudden_Underflow
    }
    else {
        L = -L >> Exp_shift;
        if (L < Exp_shift) {
            word0(a) = 0x80000u >> L;
            word1(a) = 0;
        }
        else {
            word0(a) = 0;
            L -= Exp_shift;
            word1(a) = L >= 31 ? 1u : 1u << (31 - L);
        }
    }
#endif
    return a;
}


// Return the number (0 through 32) of most significant zero bits in x.
int JS::hi0bits(uint32 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;
}


// Return the number (0 through 32) of least significant zero bits in y.
// Also shift y to the right past these 0 through 32 zeros so that y's
// least significant bit will be set unless y was originally zero.
static int lo0bits(uint32 &y)
{
    int k;
    uint32 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;
}


uint32 *JS::BigInt::freeLists[maxLgGrossSize+1];


// Allocate a BigInt with 2^lgGrossSize words.  The BigInt must not currently
// contain a number.
void JS::BigInt::allocate(uint lgGrossSize)
{
    ASSERT(lgGrossSize <= maxLgGrossSize);

    BigInt::lgGrossSize = lgGrossSize;
    negative = false;
    grossSize = 1u << lgGrossSize;
    size = 0;
    words = 0;
    ACQUIRE_DTOA_LOCK(0);
    uint32 *w = freeLists[lgGrossSize];
    if (w) {
        freeLists[lgGrossSize] = *reinterpret_cast<uint32 **>(w);
        FREE_DTOA_LOCK(0);
    } else {
        FREE_DTOA_LOCK(0);
        w = static_cast<uint32 *>(STD::malloc(v_max(uint32(grossSize * sizeof(uint32)), uint32(sizeof(uint32 *)))));
        if (!w) {
            std::bad_alloc outOfMemory;
            throw outOfMemory;
        }
    }
    words = w;
}


// Recycle this BigInt's words array, which must be non-nil.
void JS::BigInt::recycle()
{
    ASSERT(words);
    uint bucket = lgGrossSize;
    uint32 *w = words;
    ACQUIRE_DTOA_LOCK(0);
    *reinterpret_cast<uint32 **>(w) = freeLists[bucket];
    freeLists[bucket] = w;
    FREE_DTOA_LOCK(0);
}


// Copy b into this BigInt, which must be uninitialized.
void JS::BigInt::initCopy(const BigInt &b)
{
    allocate(b.lgGrossSize);
    negative = b.negative;
    size = b.size;
    std::copy(b.words, b.words+size, words);
}


// Move b into this BigInt.  The original words array of this BigInt is recycled
// and must be non-nil.  After the move b has no value.
void JS::BigInt::move(BigInt &b)
{
    recycle();
    lgGrossSize = b.lgGrossSize;
    negative = b.negative;
    grossSize = b.grossSize;
    size = b.size;
    words = b.words;
    b.words = 0;
}


// Change this BigInt's lgGrossSize to the given value without affecting the
// BigInt value contained.  This BigInt must currently have a value that will
// fit in the new lgGrossSize.
void JS::BigInt::setLgGrossSize(uint lgGrossSize)
{
    ASSERT(words);
    BigInt b(lgGrossSize);
    ASSERT(size <= b.grossSize);
    std::copy(words, words+size, b.words);
    move(b);
}


// Set the BigInt to the given integer value.
// The BigInt must not currently have a value.
void JS::BigInt::init(uint32 i)
{
    ASSERT(!words);
    allocate(1);    // Allocate two words to allow a little room for growth.
    if (i) {
        size = 1;
        words[0] = i;
    } else
        size = 0;
}


// Convert d into the form b*2^e, where b is an odd integer.  b is the BigInt
// returned in this, and e is the returned binary exponent.  Return the number
// of significant bits in b in bits.  d must be finite and nonzero.
void JS::BigInt::init(double d, int32 &e, int32 &bits)
{
    allocate(1);

    uint32 *x = words;
    uint32 z = word0(d) & Frac_mask;
    word0(d) &= 0x7fffffff;   // clear sign bit, which we ignore
    int32 de = (int32)(word0(d) >> Exp_shift);
#ifdef Sudden_Underflow
    z |= Exp_msk11;
#else
    if (de)
        z |= Exp_msk1;
#endif
    uint32 y = word1(d);
    int k;
    int i;
    if (y) {
        if ((k = lo0bits(y)) != 0) {
            x[0] = y | z << (32 - k);
            z >>= k;
        }
        else
            x[0] = y;
        i = (x[1] = z) ? 2 : 1;
    } else {
        ASSERT(z);
        k = lo0bits(z);
        x[0] = z;
        i = 1;
        k += 32;
    }
    size = uint32(i);

#ifndef Sudden_Underflow
    if (de) {
#endif
        e = de - Bias - (P-1) + k;
        bits = P - k;
#ifndef Sudden_Underflow
    } else {
        e = de - Bias - (P-1) + 1 + k;
        bits = 32*i - hi0bits(x[i-1]);
    }
#endif
}


// Let this = this*m + a.  Both a and m must be between 0 and 65535 inclusive.
void JS::BigInt::mulAdd(uint m, uint a)
{
#ifdef NATIVE_INT64
    uint64 carry, y;
#else
    uint32 carry, y;
    uint32 xi, z;
#endif

    ASSERT(m <= 0xFFFF && a <= 0xFFFF);
    uint32 sz = size;
    uint32 *x = words;
    carry = a;
    for (uint32 i = 0; i != sz; i++) {
#ifdef NATIVE_INT64
        y = *x * (uint64)m + carry;
        carry = y >> 32;
        *x++ = (uint32)(y & 0xffffffffUL);
#else
        xi = *x;
        y = (xi & 0xffff) * m + carry;
        z = (xi >> 16) * m + (y >> 16);
        carry = z >> 16;
        *x++ = (z << 16) + (y & 0xffff);
#endif
    }
    if (carry) {
        if (sz >= grossSize)
            setLgGrossSize(lgGrossSize+1);
        words[sz++] = (uint32)carry;
        size = sz;
    }
}


// Let this = this*m.
void JS::BigInt::operator*=(const BigInt &m)
{
    const BigInt *a;
    const BigInt *b;
    if (size >= m.size) {
        a = this;
        b = &m;
    } else {
        a = &m;
        b = this;
    }

    uint k = a->lgGrossSize;
    uint32 wa = a->size;
    uint32 wb = b->size;
    uint32 wc = wa + wb;
    if (wc > a->grossSize)
        k++;

    BigInt c(k);
    uint32 *xc;
    uint32 *xce;
    for (xc = c.words, xce = xc + wc; xc < xce; xc++)
        *xc = 0;
    const uint32 *xa = a->words;
    const uint32 *xae = xa + wa;
    const uint32 *xb = b->words;
    const uint32 *xbe = xb + wb;
    uint32 *xc0 = c.words;
    uint32 y;
#ifdef NATIVE_INT64
    for (; xb < xbe; xc0++) {
        if ((y = *xb++) != 0) {
            const uint32 *x = xa;
            xc = xc0;
            uint64 carry = 0;
            do {
                uint64 z = *x++ * (uint64)y + *xc + carry;
                carry = z >> 32;
                *xc++ = (uint32)(z & 0xffffffffUL);
            } while (x < xae);
            *xc = (uint32)carry;
        }
    }
#else
    for (; xb < xbe; xb++, xc0++) {
        uint32 carry;
        uint32 z, z2;
        const uint32 *x;

        if ((y = *xb & 0xffff) != 0) {
            x = xa;
            xc = xc0;
            carry = 0;
            do {
                z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
                carry = z >> 16;
                z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
                carry = z2 >> 16;
                Storeinc(xc, z2, z);
            } while (x < xae);
            *xc = carry;
        }
        if ((y = *xb >> 16) != 0) {
            x = xa;
            xc = xc0;
            carry = 0;
            z2 = *xc;
            do {
                z = (*x & 0xffff) * y + (*xc >> 16) + carry;
                carry = z >> 16;
                Storeinc(xc, z, z2);
                z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
                carry = z2 >> 16;
            } while (x < xae);
            *xc = z2;
        }
    }
#endif
    for (xc0 = c.words, xc = xc0 + wc; wc && !*--xc; --wc) ;
    c.size = wc;
    move(c);
}


// Let this = this * 2^k.  k must be nonnegative.
void JS::BigInt::pow2Mul(int32 k)
{
    ASSERT(k >= 0);

    uint32 n = uint32(k) >> 5;
    uint k1 = lgGrossSize;
    uint32 n1 = n + size + 1;
    uint32 i;
    for (i = grossSize; n1 > i; i <<= 1)
        k1++;

    BigInt b1(k1);
    uint32 *x1 = b1.words;
    for (i = 0; i < n; i++)
        *x1++ = 0;
    uint32 *x = words;
    uint32 *xe = x + size;
    if (k &= 0x1f) {
        k1 = uint(32 - k);
        uint32 z = 0;
        while (x != xe) {
            *x1++ = *x << k | z;
            z = *x++ >> k1;
        }
        if ((*x1 = z) != 0)
            ++n1;
    }
    else
        while (x != xe)
            *x1++ = *x++;
    b1.size = n1 - 1;
    move(b1);
}


// 'p5s' points to a linked list of Bigints that are powers of 625.
// This list grows on demand, and it can only grow: it won't change
// in any other way.  So if we read 'p5s' or the 'next' field of
// some BigInt on the list, and it is not NULL, we know it won't
// change to NULL or some other value.  Only when the value of
// 'p5s' or 'next' is NULL do we need to acquire the lock and add
// a new BigInt to the list.
struct PowerOf625 {
    PowerOf625 *next;
    JS::BigInt b;
};

static PowerOf625 *p5s;

#ifdef JS_THREADSAFE
static PRLock *p5s_lock;
#endif


// Let this = this * 5^k.  k must be nonnegative.
void JS::BigInt::pow5Mul(int32 k)
{
    static const uint p05[3] = {5, 25, 125};

    ASSERT(k >= 0);
    int32 i = k & 3;
    if (i)
        mulAdd(p05[i-1], 0);

    if (!(k >>= 2))
        return;
    PowerOf625 *p5 = p5s;
    if (!p5) {
        auto_ptr<PowerOf625> p(new PowerOf625);
        p->next = 0;
        p->b.init(625);
#ifdef JS_THREADSAFE
        // We take great care to not call init() and recycle() while holding
        // the lock. lock and check again
        PR_Lock(p5s_lock);
        if (!(p5 = p5s)) {
#endif
            // first time
            p5 = p.release();
            p5s = p5;
#ifdef JS_THREADSAFE
        }
        PR_Unlock(p5s_lock);
#endif
    }
    for (;;) {
        if (k & 1)
            *this *= p5->b;
        if (!(k >>= 1))
            break;
        PowerOf625 *p51 = p5->next;
        if (!p51) {
            auto_ptr<PowerOf625> p(new PowerOf625);
            p->next = 0;
            p->b = p5->b;
            p->b *= p5->b;
#ifdef JS_THREADSAFE
            PR_Lock(p5s_lock);
            if (!(p51 = p5->next)) {
#endif
                p51 = p.release();
                p5->next = p51;
#ifdef JS_THREADSAFE
            }
            PR_Unlock(p5s_lock);
#endif
        }
        p5 = p51;
    }
}


// Return -1, 0, or 1 depending on whether this<b, this==b, or this>b,
// respectively.
int JS::BigInt::cmp(const BigInt &b) const
{
    uint32 i = size;
    uint32 j = b.size;
#ifdef DEBUG
    if (i > 1 && !words[i-1])
        NOT_REACHED("cmp called with words[size-1] == 0");
    if (j > 1 && !b.words[j-1])
        NOT_REACHED("cmp called with b.words[b.size-1] == 0");
#endif
    if (i != j)
        return i<j ? -1 : 1;
    const uint32 *xa0 = words;
    const uint32 *xa = xa0 + j;
    const uint32 *xb = b.words + j;
    while (xa != xa0) {
        if (*--xa != *--xb)
            return *xa < *xb ? -1 : 1;
    }
    return 0;
}


// Initialize this BigInt to m-n.  This BigInt must not have a previous value.
void JS::BigInt::initDiff(const BigInt &m, const BigInt &n)
{
    ASSERT(!words && this != &m && this != &n);

    int i = m.cmp(n);
    if (!i)
        init(0);
    else {
        const BigInt *a;
        const BigInt *b;
        if (i < 0) {
            a = &n;
            b = &m;
        } else {
            a = &m;
            b = &n;
        }
        allocate(a->lgGrossSize);
        negative = i < 0;
        uint32 wa = a->size;
        const uint32 *xa = a->words;
        const uint32 *xae = xa + wa;
        const uint32 *xb = b->words;
        const uint32 *xbe = xb + b->size;
        uint32 *xc = words;
#ifdef NATIVE_INT64
        uint64 borrow = 0;
        uint64 y;
        while (xb < xbe) {
            y = (uint64)*xa++ - *xb++ - borrow;
            borrow = y >> 32 & 1UL;
            *xc++ = (uint32)(y & 0xffffffffUL);
        }
        while (xa < xae) {
            y = *xa++ - borrow;
            borrow = y >> 32 & 1UL;
            *xc++ = (uint32)(y & 0xffffffffUL);
        }
#else
        uint32 borrow = 0;
        uint32 y;
        uint32 z;
        while (xb < xbe) {
            y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
            borrow = (y & 0x10000) >> 16;
            z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
            borrow = (z & 0x10000) >> 16;
            Storeinc(xc, z, y);
        }
        while (xa < xae) {
            y = (*xa & 0xffff) - borrow;
            borrow = (y & 0x10000) >> 16;
            z = (*xa++ >> 16) - borrow;
            borrow = (z & 0x10000) >> 16;
            Storeinc(xc, z, y);
        }
#endif
        while (!*--xc)
            wa--;
        size = wa;
    }
}


// Return floor(this/2^k) and set this to be the remainder.
// The returned quotient must be less than 2^32.
uint32 JS::BigInt::quoRem2(int32 k)
{
    int32 n = k >> 5;
    k &= 0x1F;
    uint32 mask = (1u<<k) - 1;

    int32 w = int32(size) - n;
    if (w <= 0)
        return 0;
    ASSERT(w <= 2);
    uint32 *bx = words;
    uint32 *bxe = bx + n;
    uint32 result = *bxe >> k;
    *bxe &= mask;
    if (w == 2) {
        ASSERT(!(bxe[1] & ~mask));
        if (k)
            result |= bxe[1] << (32 - k);
    }
    n++;
    while (!*bxe && bxe != bx) {
        n--;
        bxe--;
    }
    size = uint32(n);
    return result;
}


// Return floor(this/S) and set this to be the remainder.  As added restrictions,
// this must not have more words than S, the most significant word of S must not
// start with a 1 bit, and the returned quotient must be less than 36.
int32 JS::BigInt::quoRem(const BigInt &S)
{
#ifdef NATIVE_INT64
    uint64 borrow, carry, y, ys;
#else
    uint32 borrow, carry, y, ys;
    uint32 si, z, zs;
#endif

    uint32 n = S.size;
    ASSERT(size <= n && n);
    if (size < n)
        return 0;
    const uint32 *sx = S.words;
    const uint32 *sxe = sx + --n;
    uint32 *bx = words;
    uint32 *bxe = bx + n;
    ASSERT(*sxe <= 0x7FFFFFFF);
    uint32 q = *bxe / (*sxe + 1);  // ensure q <= true quotient
    ASSERT(q < 36);
    if (q) {
        borrow = 0;
        carry = 0;
        do {
#ifdef NATIVE_INT64
            ys = *sx++ * (uint64)q + carry;
            carry = ys >> 32;
            y = *bx - (ys & 0xffffffffUL) - borrow;
            borrow = y >> 32 & 1UL;
            *bx++ = (uint32)(y & 0xffffffffUL);
#else
            si = *sx++;
            ys = (si & 0xffff) * q + carry;
            zs = (si >> 16) * q + (ys >> 16);
            carry = zs >> 16;
            y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
            borrow = (y & 0x10000) >> 16;
            z = (*bx >> 16) - (zs & 0xffff) - borrow;
            borrow = (z & 0x10000) >> 16;
            Storeinc(bx, z, y);
#endif
        } while (sx <= sxe);
        if (!*bxe) {
            bx = words;
            while (--bxe > bx && !*bxe)
                --n;
            size = n;
        }
    }
    if (cmp(S) >= 0) {
        q++;
        borrow = 0;
        carry = 0;
        bx = words;
        sx = S.words;
        do {
#ifdef NATIVE_INT64
            ys = *sx++ + carry;
            carry = ys >> 32;
            y = *bx - (ys & 0xffffffffUL) - borrow;
            borrow = y >> 32 & 1UL;
            *bx++ = (uint32)(y & 0xffffffffUL);
#else
            si = *sx++;
            ys = (si & 0xffff) + carry;
            zs = (si >> 16) + (ys >> 16);
            carry = zs >> 16;
            y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
            borrow = (y & 0x10000) >> 16;
            z = (*bx >> 16) - (zs & 0xffff) - borrow;
            borrow = (z & 0x10000) >> 16;
            Storeinc(bx, z, y);
#endif
        } while (sx <= sxe);
        bx = words;
        bxe = bx + n;
        if (!*bxe) {
            while (--bxe > bx && !*bxe)
                --n;
            size = n;
        }
    }
    return int32(q);
}


// Let this = floor(this / divisor), and return the remainder.  this must be
// nonnegative. divisor must be between 1 and 65536.
uint32 JS::BigInt::divRem(uint32 divisor)
{
    uint32 n = size;
    uint32 remainder = 0;

    ASSERT(divisor > 0 && divisor <= 65536);

    if (!n)
        return 0; // this is zero
    uint32 *bx = words;
    uint32 *bp = bx + n;
    do {
        uint32 a = *--bp;
        uint32 dividend = remainder << 16 | a >> 16;
        uint32 quotientHi = dividend / divisor;
        uint32 quotientLo;

        remainder = dividend - quotientHi*divisor;
        ASSERT(quotientHi <= 0xFFFF && remainder < divisor);
        dividend = remainder << 16 | (a & 0xFFFF);
        quotientLo = dividend / divisor;
        remainder = dividend - quotientLo*divisor;
        ASSERT(quotientLo <= 0xFFFF && remainder < divisor);
        *bp = quotientHi << 16 | quotientLo;
    } while (bp != bx);

    // Decrease the size of the number if its most significant word is now zero.
    if (bx[n-1] == 0)
        size--;
    return remainder;
}


double JS::BigInt::b2d(int32 &e) const
{
    double d;

    const uint32 *xa0 = words;
    const uint32 *xa = xa0 + size;
    ASSERT(size);
    uint32 y = *--xa;
    ASSERT(y);
    int k = hi0bits(y);
    e = 32 - k;
    if (k < Ebits) {
        word0(d) = Exp_1 | y >> (Ebits - k);
        uint32 w = xa > xa0 ? *--xa : 0;
        word1(d) = y << (32-Ebits + k) | w >> (Ebits - k);
    } else {
        uint32 z = xa > xa0 ? *--xa : 0;
        if (k -= Ebits) {
            word0(d) = Exp_1 | y << k | z >> (32 - k);
            y = xa > xa0 ? *--xa : 0;
            word1(d) = z << k | y >> (32 - k);
        }
        else {
            word0(d) = Exp_1 | y;
            word1(d) = z;
        }
    }
    return d;
}


double JS::BigInt::ratio(const BigInt &denominator) const
{
    int32 ka, kb;

    double da = b2d(ka);
    double db = denominator.b2d(kb);
    int32 k = ka - kb + 32*int32(size - denominator.size);
    if (k > 0)
        word0(da) += k*Exp_msk1;
    else
        word0(db) -= k*Exp_msk1;
    return da / db;
}


void JS::BigInt::s2b(const char *s, int32 nd0, int32 nd, uint32 y9)
{
    int32 i;
    int32 x, y;

    x = (nd + 8) / 9;
    uint k = 0;
    for (y = 1; x > y; y <<= 1, k++) ;
    ASSERT(!words);
    allocate(k);
    words[0] = y9;
    size = 1;

    i = 9;
    if (9 < nd0) {
        s += 9;
        do mulAdd(10, uint(*s++ - '0'));
        while (++i < nd0);
        s++;
    } else
        s += 10;
    for (; i < nd; i++)
        mulAdd(10, uint(*s++ - '0'));
}



static const double tens[] = {
    1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
        1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
        1e20, 1e21, 1e22
        };

static const double bigtens[] = {1e16, 1e32, 1e64, 1e128, 1e256};
static const double tinytens[] = {1e-16, 1e-32, 1e-64, 1e-128,
#ifdef Avoid_Underflow
                                  9007199254740992.e-256
#else
                                  1e-256
#endif
                                  };
// The factor of 2^53 in tinytens[4] helps us avoid setting the underflow
// flag unnecessarily.  It leads to a song and dance at the end of strToDouble.

const int32 Scale_Bit = 0x10;
const int n_bigtens = 5;


#ifdef JS_THREADSAFE
static bool initialized = false;

// hacked replica of nspr _PR_InitDtoa
static void InitDtoa(void)
{
    freelist_lock = PR_NewLock();
    p5s_lock = PR_NewLock();
    initialized = true;
}
#endif



// Return as a double-precision floating-point number the value represented by
// the character string str.  The string is scanned up to the first unrecognized
// character.  The character sequences 'Infinity', '+Infinity', '-Infinity', and
// 'NaN' are also recognized.
// Return a pointer to the character terminating the scan in numEnd.
// If no number can be formed, set numEnd to str and return zero.
double JS::strToDouble(const char *str, const char *&numEnd)
{
    int32 scale;
    int32 bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
        e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0;
    const char *s, *s0, *s1;
    double aadj, aadj1, adj, rv, rv0;
    int32 L;
    uint32 y, z;

#ifdef JS_THREADSAFE
    if (!initialized) InitDtoa();
#endif

    nz0 = nz = 0;
    bool negative = false;
    bool haveSign = false;
    rv = 0.;
    for (s = str;; s++)
        switch (*s) {
            case '-':
                negative = true;
                // no break
            case '+':
                haveSign = true;
                if (*++s)
                    goto break2;
                // no break
            case 0:
                s = str;
                goto ret;
            case '\t':
            case '\n':
            case '\v':
            case '\f':
            case '\r':
            case ' ':
                continue;
            default:
                goto break2;
        }
  break2:
    switch (*s) {
        case '0':
            nz0 = 1;
            while (*++s == '0') ;
            if (!*s)
                goto ret;
            break;
        case 'I':
            if (!STD::strncmp(s+1, "nfinity", 7)) {
                rv = positiveInfinity;
                s += 8;
                goto ret;
            }
            break;
        case 'N':
            if (!haveSign && !STD::strncmp(s+1, "aN", 2)) {
                rv = nan;
                s += 3;
                goto ret;
            }
    }
    s0 = s;
    y = z = 0;
    for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
        if (nd < 9)
            y = 10*y + c - '0';
        else if (nd < 16)
            z = 10*z + c - '0';
    nd0 = nd;
    if (c == '.') {
        c = *++s;
        if (!nd) {
            for (; c == '0'; c = *++s)
                nz++;
            if (c > '0' && c <= '9') {
                s0 = s;
                nf += nz;
                nz = 0;
                goto have_dig;
            }
            goto dig_done;
        }
        for (; c >= '0' && c <= '9'; c = *++s) {
          have_dig:
            nz++;
            if (c -= '0') {
                nf += nz;
                for (i = 1; i < nz; i++)
                    if (nd++ < 9)
                        y *= 10;
                    else if (nd <= DBL_DIG + 1)
                        z *= 10;
                if (nd++ < 9)
                    y = 10*y + c;
                else if (nd <= DBL_DIG + 1)
                    z = 10*z + c;
                nz = 0;
            }
        }
    }
  dig_done:
    e = 0;
    if (c == 'e' || c == 'E') {
        if (!nd && !nz && !nz0) {
            s = str;
            goto ret;
        }
        str = s;
        esign = 0;
        switch (c = *++s) {
            case '-':
                esign = 1;
            case '+':
                c = *++s;
        }
        if (c >= '0' && c <= '9') {
            while (c == '0')
                c = *++s;
            if (c > '0' && c <= '9') {
                L = c - '0';
                s1 = s;
                while ((c = *++s) >= '0' && c <= '9')
                    L = 10*L + c - '0';
                if (s - s1 > 8 || L > 19999)
                    // Avoid confusion from exponents so large that e might overflow.
                    e = 19999; // safe for 16 bit ints
                else
                    e = (int32)L;
                if (esign)
                    e = -e;
            }
            else
                e = 0;
        }
        else
            s = str;
    }
    if (!nd) {
        if (!nz && !nz0)
            s = str;
        goto ret;
    }
    e1 = e -= nf;

    // Now we have nd0 digits, starting at s0, followed by a
    // decimal point, followed by nd-nd0 digits.  The number we're
    // after is the integer represented by those digits times 10**e

    if (!nd0)
        nd0 = nd;
    k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
    rv = y;
    if (k > 9)
        rv = tens[k - 9] * rv + z;
    if (nd <= DBL_DIG && FLT_ROUNDS == 1) {
        if (!e)
            goto ret;
        if (e > 0) {
            if (e <= Ten_pmax) {
                rv *= tens[e];
                goto ret;
            }
            i = DBL_DIG - nd;
            if (e <= Ten_pmax + i) {
                // A fancier test would sometimes let us do this for larger i values.
                e -= i;
                rv *= tens[i];
                rv *= tens[e];
                goto ret;
            }
        }
#ifndef Inaccurate_Divide
        else if (e >= -Ten_pmax) {
            rv /= tens[-e];
            goto ret;
        }
#endif
    }
    e1 += nd - k;

    scale = 0;

    // Get starting approximation = rv * 10**e1

    if (e1 > 0) {
        if ((i = e1 & 15) != 0)
            rv *= tens[i];
        if (e1 &= ~15) {
            if (e1 > DBL_MAX_10_EXP) {
              ovfl:
                // Return infinity.
                rv = positiveInfinity;
                goto ret;
            }
            e1 >>= 4;
            for (j = 0; e1 > 1; j++, e1 >>= 1)
                if (e1 & 1)
                    rv *= bigtens[j];
            // The last multiplication could overflow.
            word0(rv) -= P*Exp_msk1;
            rv *= bigtens[j];
            if ((z = word0(rv) & Exp_mask) > Exp_msk1*(DBL_MAX_EXP+Bias-P))
                goto ovfl;
            if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
                // set to largest number (Can't trust DBL_MAX)
                word0(rv) = Big0;
                word1(rv) = Big1;
            }
            else
                word0(rv) += P*Exp_msk1;
        }
    }
    else if (e1 < 0) {
        e1 = -e1;
        if ((i = e1 & 15) != 0)
            rv /= tens[i];
        if (e1 &= ~15) {
            e1 >>= 4;
            if (e1 >= 1 << n_bigtens)
                goto undfl;
#ifdef Avoid_Underflow
            if (e1 & Scale_Bit)
                scale = P;
            for (j = 0; e1 > 0; j++, e1 >>= 1)
                if (e1 & 1)
                    rv *= tinytens[j];
            if (scale && (j = P + 1 - int32((word0(rv) & Exp_mask) >> Exp_shift)) > 0) {
                // scaled rv is denormal; zap j low bits
                if (j >= 32) {
                    word1(rv) = 0;
                    word0(rv) &= 0xffffffff << (j-32);
                    if (!word0(rv))
                        word0(rv) = 1;
                }
                else
                    word1(rv) &= 0xffffffff << j;
            }
#else
            for (j = 0; e1 > 1; j++, e1 >>= 1)
                if (e1 & 1)
                    rv *= tinytens[j];
            // The last multiplication could underflow.
            rv0 = rv;
            rv *= tinytens[j];
            if (!rv) {
                rv = 2.*rv0;
                rv *= tinytens[j];
#endif
                if (!rv) {
                  undfl:
                    rv = 0.;
                    goto ret;
                }
#ifndef Avoid_Underflow
                word0(rv) = Tiny0;
                word1(rv) = Tiny1;
                // The refinement below will clean this approximation up.
            }
#endif
        }
    }

    // Now the hard part -- adjusting rv to the correct value.
    // Put digits into bd: true value = bd * 10^e
    {
        BigInt bd0;
        bd0.s2b(s0, nd0, nd, y);

        for (;;) {
            BigInt bs;

            BigInt bd = bd0;
            BigInt bb;
            bb.init(rv, bbe, bbbits);    // rv = bb * 2^bbe
            bs.init(1);

            if (e >= 0) {
                bb2 = bb5 = 0;
                bd2 = bd5 = e;
            }
            else {
                bb2 = bb5 = -e;
                bd2 = bd5 = 0;
            }
            if (bbe >= 0)
                bb2 += bbe;
            else
                bd2 -= bbe;
            bs2 = bb2;
#ifdef Sudden_Underflow
            j = P + 1 - bbbits;
#else
#ifdef Avoid_Underflow
            j = bbe - scale;
#else
            j = bbe;
#endif
            i = j + bbbits - 1; // logb(rv)
            if (i < Emin)   // denormal
                j += P - Emin;
            else
                j = P + 1 - bbbits;
#endif
            bb2 += j;
            bd2 += j;
#ifdef Avoid_Underflow
            bd2 += scale;
#endif
            i = bb2 < bd2 ? bb2 : bd2;
            if (i > bs2)
                i = bs2;
            if (i > 0) {
                bb2 -= i;
                bd2 -= i;
                bs2 -= i;
            }
            if (bb5 > 0) {
                bs.pow5Mul(bb5);
                bb *= bs;
            }
            if (bb2 > 0)
                bb.pow2Mul(bb2);
            if (bd5 > 0)
                bd.pow5Mul(bd5);
            if (bd2 > 0)
                bd.pow2Mul(bd2);
            if (bs2 > 0)
                bs.pow2Mul(bs2);

            BigInt delta;
            delta.initDiff(bb, bd);
            dsign = delta.negative;
            delta.negative = false;
            i = delta.cmp(bs);
            if (i < 0) {
                // Error is less than half an ulp -- check for special case of mantissa a power of two.
                if (dsign || word1(rv) || word0(rv) & Bndry_mask
#ifdef Avoid_Underflow
                    || (word0(rv) & Exp_mask) <= Exp_msk1 + P*Exp_msk1
#else
                    || (word0(rv) & Exp_mask) <= Exp_msk1
#endif
                    ) {
#ifdef Avoid_Underflow
                    if (delta.isZero())
                        dsign = 2;
#endif
                    break;
                }
                delta.pow2Mul(Log2P);
                if (delta.cmp(bs) > 0)
                    goto drop_down;
                break;
            }
            if (i == 0) {
                // exactly half-way between
                if (dsign) {
                    if ((word0(rv) & Bndry_mask1) == Bndry_mask1 && word1(rv) == 0xffffffff) {
                        // boundary case -- increment exponent
                        word0(rv) = (word0(rv) & Exp_mask) + Exp_msk1;
                        word1(rv) = 0;
#ifdef Avoid_Underflow
                        dsign = 0;
#endif
                        break;
                    }
                }
                else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
#ifdef Avoid_Underflow
                    dsign = 2;
#endif
                  drop_down:
                    // boundary case -- decrement exponent
#ifdef Sudden_Underflow
                    L = word0(rv) & Exp_mask;
                    if (L <= Exp_msk1)
                        goto undfl;
                    L -= Exp_msk1;
#else
                    L = int32((word0(rv) & Exp_mask) - Exp_msk1);
#endif
                    word0(rv) = uint32(L) | Bndry_mask1;
                    word1(rv) = 0xffffffff;
                    break;
                }
#ifndef ROUND_BIASED
                if (!(word1(rv) & LSB))
                    break;
#endif
                if (dsign)
                    rv += ulp(rv);
#ifndef ROUND_BIASED
                else {
                    rv -= ulp(rv);
#ifndef Sudden_Underflow
                    if (!rv)
                        goto undfl;
#endif
                }
#ifdef Avoid_Underflow
                dsign = 1 - dsign;
#endif
#endif
                break;
            }
            if ((aadj = delta.ratio(bs)) <= 2.) {
                if (dsign)
                    aadj = aadj1 = 1.;
                else if (word1(rv) || word0(rv) & Bndry_mask) {
#ifndef Sudden_Underflow
                    if (word1(rv) == Tiny1 && !word0(rv))
                        goto undfl;
#endif
                    aadj = 1.;
                    aadj1 = -1.;
                }
                else {
                    // special case -- power of FLT_RADIX to be rounded down...
                    if (aadj < 2./FLT_RADIX)
                        aadj = 1./FLT_RADIX;
                    else
                        aadj *= 0.5;
                    aadj1 = -aadj;
                }
            }
            else {
                aadj *= 0.5;
                aadj1 = dsign ? aadj : -aadj;
#ifdef Check_FLT_ROUNDS
                switch (FLT_ROUNDS) {
                    case 2: // towards +infinity
                        aadj1 -= 0.5;
                        break;
                    case 0: // towards 0
                    case 3: // towards -infinity
                        aadj1 += 0.5;
                }
#else
                if (FLT_ROUNDS == 0)
                    aadj1 += 0.5;
#endif
            }
            y = word0(rv) & Exp_mask;

            // Check for overflow

            if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
                rv0 = rv;
                word0(rv) -= P*Exp_msk1;
                adj = aadj1 * ulp(rv);
                rv += adj;
                if ((word0(rv) & Exp_mask) >= Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
                    if (word0(rv0) == Big0 && word1(rv0) == Big1)
                        goto ovfl;
                    word0(rv) = Big0;
                    word1(rv) = Big1;
                    continue;
                }
                else
                    word0(rv) += P*Exp_msk1;
            }
            else {
#ifdef Sudden_Underflow
                if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
                    rv0 = rv;
                    word0(rv) += P*Exp_msk1;
                    adj = aadj1 * ulp(rv);
                    rv += adj;
                    if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
                        if (word0(rv0) == Tiny0 && word1(rv0) == Tiny1)
                            goto undfl;
                        word0(rv) = Tiny0;
                        word1(rv) = Tiny1;
                        continue;
                    } else
                        word0(rv) -= P*Exp_msk1;
                } else {
                    adj = aadj1 * ulp(rv);
                    rv += adj;
                }
#else
                // Compute adj so that the IEEE rounding rules will correctly round rv + adj in some half-way cases.
                // If rv * ulp(rv) is denormalized (i.e., y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
                // trouble from bits lost to denormalization; example: 1.2e-307.
#ifdef Avoid_Underflow
                if (y <= P*Exp_msk1 && aadj > 1.)
#else
                    if (y <= (P-1)*Exp_msk1 && aadj > 1.)
#endif
                    {
                        aadj1 = (double)(int32)(aadj + 0.5);
                        if (!dsign)
                            aadj1 = -aadj1;
                    }
#ifdef Avoid_Underflow
                if (scale && y <= P*Exp_msk1)
                    word0(aadj1) += (P+1)*Exp_msk1 - y;
#endif
                adj = aadj1 * ulp(rv);
                rv += adj;
#endif
            }
            z = word0(rv) & Exp_mask;
#ifdef Avoid_Underflow
            if (!scale)
#endif
                if (y == z) {
                    // Can we stop now?
                    L = (int32)aadj;
                    aadj -= L;
                    // The tolerances below are conservative.
                    if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
                        if (aadj < .4999999 || aadj > .5000001)
                            break;
                    }
                    else if (aadj < .4999999/FLT_RADIX)
                        break;
                }
        }
#ifdef Avoid_Underflow
        if (scale) {
            word0(rv0) = Exp_1 - P*Exp_msk1;
            word1(rv0) = 0;
            if ((word0(rv) & Exp_mask) <= P*Exp_msk1 && word1(rv) & 1 && dsign != 2) {
                if (dsign) {
#ifdef Sudden_Underflow
                    // rv will be 0, but this would give the right result if only rv *= rv0 worked.
                    word0(rv) += P*Exp_msk1;
                    word0(rv0) = Exp_1 - 2*P*Exp_msk1;
#endif
                    rv += ulp(rv);
                }
                else
                    word1(rv) &= ~1;
            }
            rv *= rv0;
        }
#endif // Avoid_Underflow
    }
  ret:
    numEnd = s;
    return negative ? -rv : rv;
}


// A version of strToDouble that takes a char16 string that begins at str and
// ends just before strEnd.  The char16 string does not have to be
// null-terminated. Leading Unicode whitespace is skipped.
double JS::stringToDouble(const char16 *str, const char16 *strEnd, const char16 *&numEnd)
{
    const char16 *str1 = skipWhiteSpace(str, strEnd);

    CharAutoPtr cstr(new char[strEnd - str1 + 1]);
    char *q = cstr.get();
    for (const char16 *p = str1; p != strEnd; p++) {
        char16 ch = *p;
        if (uint16(ch) >> 8)
            break;
        *q++ = char(ch);
    }
    *q = '\0';

    const char *estr;
    double value = strToDouble(cstr.get(), estr);
    ptrdiff_t i = estr - cstr.get();
    numEnd = i ? str1 + i : str;
    return value;
}



class BinaryDigitReader
{
    uint base;               // Base of number; must be a power of 2
    uint digit;              // Current digit value in radix given by base
    uint digitMask;          // Mask to extract the next bit from digit
    const char16 *digits;    // Pointer to the remaining digits
    const char16 *digitsEnd; // Pointer to first non-digit

public:
    BinaryDigitReader(uint base, const char16 *digitsBegin, const char16 *digitsEnd) :
            base(base), digitMask(0), digits(digitsBegin), digitsEnd(digitsEnd) {}

    int next();
};


// Return the next binary digit from the number or -1 if done.
int BinaryDigitReader::next()
{
    if (digitMask == 0) {
        if (digits == digitsEnd)
            return -1;

        uint c = *digits++;
        if ('0' <= c && c <= '9')
            digit = c - '0';
        else if ('a' <= c && c <= 'z')
            digit = c - 'a' + 10;
        else digit = c - 'A' + 10;
        digitMask = base >> 1;
    }
    int bit = (digit & digitMask) != 0;
    digitMask >>= 1;
    return bit;
}


// Read an integer from a char16 string that begins at str and ends just before
// strEnd. The char16 string does not have to be null-terminated.  The integer
// is returned as a double, which is guaranteed to be the closest double number
// to the given input when base is 10 or a power of 2.
// May experience roundoff errors for very large numbers of a different radix.
// Return a pointer to the character just past the integer in numEnd.
// If the string does not have a number in it, set numEnd to str and return 0.
// Leading Unicode whitespace is skipped.
double JS::stringToInteger(const char16 *str, const char16 *strEnd, const char16 *&numEnd, uint base)
{
    const char16 *str1 = skipWhiteSpace(str, strEnd);

    bool negative = (*str1 == '-');
    if (negative || *str1 == '+')
        str1++;

    if ((base == 0 || base == 16) &&
        *str1 == '0' && (str1[1] == 'X' || str1[1] == 'x')) {
        // Skip past hex prefix.
        base = 16;
        str1 += 2;
    }
    if (base == 0)
        base = 10; // Default to decimal.

        // Find some prefix of the string that's a number in the given base.
    const char16 *start = str1; // Mark - if string is empty, we return 0.
    double value = 0.0;
    while (true) {
        uint digit;
        char16 c = *str1;
        if ('0' <= c && c <= '9')
            digit = uint(c) - '0';
        else if ('a' <= c && c <= 'z')
            digit = uint(c) - 'a' + 10;
        else if ('A' <= c && c <= 'Z')
            digit = uint(c) - 'A' + 10;
        else
            break;
        if (digit >= base)
            break;
        value = value*base + digit;
        str1++;
    }

    if (value >= 9007199254740992.0) {
        if (base == 10) {
            // If we're accumulating a decimal number and the number is >= 2^53, then the result from the
            // repeated multiply-add above may be inaccurate.  Call stringToDouble to get the correct answer.
            const char16 *numEnd2;
            value = stringToDouble(start, str1, numEnd2);
            ASSERT(numEnd2 == str1);

        } else if (base == 2 || base == 4 || base == 8 || base == 16 || base == 32) {
            // The number may also be inaccurate for one of these bases.  This happens if the addition in
            // value*base + digit causes a round-down to an even least significant mantissa bit when the
            // first dropped bit is a one.  If any of the following digits in the number (which haven't been
            // added in yet) are nonzero then the correct action would have been to round up instead of down.
            // An example of this occurs when reading the number 0x1000000000000081, which rounds to
            // 0x1000000000000000 instead of 0x1000000000000100.
            BinaryDigitReader bdr(base, start, str1);
            value = 0.0;

            // Skip leading zeros.
            int bit;
            do {
                bit = bdr.next();
            } while (bit == 0);

            if (bit == 1) {
                // Gather the 53 significant bits (including the leading 1)
                int bit2;
                value = 1.0;
                for (int j = 52; j; --j) {
                    bit = bdr.next();
                    if (bit < 0)
                        goto done;
                    value = value*2 + bit;
                }
                bit2 = bdr.next();      // bit2 is the 54th bit (the first dropped from the mantissa)
                if (bit2 >= 0) {
                    double factor = 2.0;
                    int sticky = 0;     // sticky is 1 if any bit beyond the 54th is 1
                    int bit3;

                    while ((bit3 = bdr.next()) >= 0) {
                        sticky |= bit3;
                        factor *= 2;
                    }
                    value += bit2 & (bit | sticky);
                    value *= factor;
                }
              done:;
            }
        }
    }
    // We don't worry about inaccurate numbers for any other base.

    if (str1 == start)
        numEnd = str;
    else {
        numEnd = str1;
        if (negative)
            value = -value;
    }
    return value;
}



// doubleToAscii for IEEE arithmetic (dmg): convert double to ASCII string.
//
// Inspired by "How to Print Floating-Point Numbers Accurately" by
// Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
//
// Modifications:
//  1. Rather than iterating, we use a simple numeric overestimate
//     to determine k = floor(log10(d)).  We scale relevant
//     quantities using O(log2(k)) rather than O(k) multiplications.
//  2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
//     try to generate digits strictly left to right.  Instead, we
//     compute with fewer bits and propagate the carry if necessary
//     when rounding the final digit up.  This is often faster.
//  3. Under the assumption that input will be rounded nearest,
//     mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
//     That is, we allow equality in stopping tests when the
//     round-nearest rule will give the same floating-point value
//     as would satisfaction of the stopping test with strict
//     inequality.
//  4. We remove common factors of powers of 2 from relevant
//     quantities.
//  5. When converting floating-point integers less than 1e16,
//     we use floating-point arithmetic rather than resorting
//     to multiple-precision integers.
//  6. When asked to produce fewer than 15 digits, we first try
//     to get by with floating-point arithmetic; we resort to
//     multiple-precision integer arithmetic only if we cannot
//     guarantee that the floating-point calculation has given
//     the correctly rounded result.  For k requested digits and
//     "uniformly" distributed input, the probability is
//     something like 10^(k-15) that we must resort to the int32
//     calculation.
///

// Always emits at least one digit.
// If biasUp is set, then rounding in modes 2 and 3 will round away from zero
// when the number is exactly halfway between two representable values.  For example,
// rounding 2.5 to zero digits after the decimal point will return 3 and not 2.
// 2.49 will still round to 2, and 2.51 will still round to 3.
// The buffer should be at least 20 bytes for modes 0 and 1.  For the other modes,
// the buffer's size should be two greater than the maximum number of output characters expected.
// Return a pointer to the resulting string's trailing null.
static char *doubleToAscii(double d, int mode, bool biasUp, int ndigits, int *decpt, bool *negative, char *buf)
{
        /*  Arguments ndigits, decpt, negative are similar to those
            of ecvt and fcvt; trailing zeros are suppressed from
            the returned string.  If d is +-Infinity or NaN,
            then *decpt is set to 9999.

            mode:
            0 ==> shortest string that yields d when read in
            and rounded to nearest.
            1 ==> like 0, but with Steele & White stopping rule;
            e.g. with IEEE P754 arithmetic , mode 0 gives
            1e23 whereas mode 1 gives 9.999999999999999e22.
            2 ==> max(1,ndigits) significant digits.  This gives a
            return value similar to that of ecvt, except
            that trailing zeros are suppressed.
            3 ==> through ndigits past the decimal point.  This
            gives a return value similar to that from fcvt,
            except that trailing zeros are suppressed, and
            ndigits can be negative.
            4-9 should give the same return values as 2-3, i.e.,
            4 <= mode <= 9 ==> same return as mode
            2 + (mode & 1).  These modes are mainly for
            debugging; often they run slower but sometimes
            faster than modes 2-3.
            4,5,8,9 ==> left-to-right digit generation.
            6-9 ==> don't try fast floating-point estimate
            (if applicable).

            Values of mode other than 0-9 are treated as mode 0.

            Sufficient space is allocated to the return value
            to hold the suppressed trailing zeros.
*/

    int32 bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
        j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, try_quick;
    int32 L;
#ifndef Sudden_Underflow
    int32 denorm;
    uint32 x;
#endif
    double d2, ds, eps;
    char *s;

#ifdef JS_THREADSAFE
    if (!initialized) InitDtoa();
#endif

    if (word0(d) & Sign_bit) {
        // set negative for everything, including 0's and NaNs
        *negative = true;
        word0(d) &= ~Sign_bit;  // clear sign bit
    }
    else
        *negative = false;

    if ((word0(d) & Exp_mask) == Exp_mask) {
        // Infinity or NaN
        *decpt = 9999;
        strcpy(buf, !word1(d) && !(word0(d) & Frac_mask) ? "Infinity" : "NaN");
        return buf[3] ? buf + 8 : buf + 3;
    }
    if (!d) {
      no_digits:
        *decpt = 1;
        buf[0] = '0'; buf[1] = '\0';  // copy "0" to buffer
        return buf + 1;
    }

    BigInt b;
    b.init(d, be, bbits);
#ifdef Sudden_Underflow
    i = (int32)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
#else
    if ((i = (int32)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) {
#endif
        d2 = d;
        word0(d2) &= Frac_mask1;
        word0(d2) |= Exp_11;

        /* log(x)   ~=~ log(1.5) + (x-1.5)/1.5
         * log10(x)  =  log(x) / log(10)
         *      ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
         * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
         *
         * This suggests computing an approximation k to log10(d) by
         *
         * k = (i - Bias)*0.301029995663981
         *  + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
         *
         * We want k to be too large rather than too small.
         * The error in the first-order Taylor series approximation
         * is in our favor, so we just round up the constant enough
         * to compensate for any error in the multiplication of
         * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
         * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
         * adding 1e-13 to the constant term more than suffices.
         * Hence we adjust the constant term to 0.1760912590558.
         * (We could get a more accurate k by invoking log10,
         *  but this is probably not worthwhile.)
         */

        i -= Bias;
#ifndef Sudden_Underflow
        denorm = 0;
    }
    else {
        // d is denormalized

        i = bbits + be + (Bias + (P-1) - 1);
        x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32) : word1(d) << (32 - i);
        d2 = x;
        word0(d2) -= 31*Exp_msk1; // adjust exponent
        i -= (Bias + (P-1) - 1) + 1;
        denorm = 1;
    }
#endif

    // At this point d = f*2^i, where 1 <= f < 2.  d2 is an approximation of f.
    ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
    k = (int32)ds;
    if (ds < 0. && ds != k)
        k--;    // want k = floor(ds)
    k_check = 1;
    if (k >= 0 && k <= Ten_pmax) {
        if (d < tens[k])
            k--;
        k_check = 0;
    }

    // At this point floor(log10(d)) <= k <= floor(log10(d))+1.
    // If k_check is zero, we're guaranteed that k = floor(log10(d)).
    j = bbits - i - 1;
    // At this point d = b/2^j, where b is an odd integer.
    if (j >= 0) {
        b2 = 0;
        s2 = j;
    }
    else {
        b2 = -j;
        s2 = 0;
    }
    if (k >= 0) {
        b5 = 0;
        s5 = k;
        s2 += k;
    }
    else {
        b2 -= k;
        b5 = -k;
        s5 = 0;
    }

    // At this point d/10^k = (b * 2^b2 * 5^b5) / (2^s2 * 5^s5), where b is an odd integer,
    // b2 >= 0, b5 >= 0, s2 >= 0, and s5 >= 0.
    if (mode < 0 || mode > 9)
        mode = 0;
    try_quick = 1;
    if (mode > 5) {
        mode -= 4;
        try_quick = 0;
    }
    leftright = 1;
    ilim = ilim1 = 0;
    switch (mode) {
        case 0:
        case 1:
            ilim = ilim1 = -1;
            i = 18;
            ndigits = 0;
            break;
        case 2:
            leftright = 0;
            // no break
        case 4:
            if (ndigits <= 0)
                ndigits = 1;
            ilim = ilim1 = i = ndigits;
            break;
        case 3:
            leftright = 0;
            // no break
        case 5:
            i = ndigits + k + 1;
            ilim = i;
            ilim1 = i - 1;
            if (i <= 0)
                i = 1;
    }
    // ilim is the maximum number of significant digits we want, based on k and ndigits.
    // ilim1 is the maximum number of significant digits we want, based on k and ndigits,
    // when it turns out that k was computed too high by one.
    s = buf;

    if (ilim >= 0 && ilim <= Quick_max && try_quick) {

        // Try to get by with floating-point arithmetic.

        i = 0;
        d2 = d;
        k0 = k;
        ilim0 = ilim;
        ieps = 2; // conservative

        // Divide d by 10^k, keeping track of the roundoff error and avoiding overflows.
        if (k > 0) {
            ds = tens[k&0xf];
            j = k >> 4;
            if (j & Bletch) {
                // prevent overflows
                j &= Bletch - 1;
                d /= bigtens[n_bigtens-1];
                ieps++;
            }
            for (; j; j >>= 1, i++)
                if (j & 1) {
                    ieps++;
                    ds *= bigtens[i];
                }
            d /= ds;
        }
        else if ((j1 = -k) != 0) {
            d *= tens[j1 & 0xf];
            for (j = j1 >> 4; j; j >>= 1, i++)
                if (j & 1) {
                    ieps++;
                    d *= bigtens[i];
                }
        }

        // Check that k was computed correctly.
        if (k_check && d < 1. && ilim > 0) {
            if (ilim1 <= 0)
                goto fast_failed;
            ilim = ilim1;
            k--;
            d *= 10.;
            ieps++;
        }

        // eps bounds the cumulative error.
        eps = ieps*d + 7.;
        word0(eps) -= (P-1)*Exp_msk1;
        if (ilim == 0) {
            d -= 5.;
            if (d > eps)
                goto one_digit;
            if (d < -eps)
                goto no_digits;
            goto fast_failed;
        }
#ifndef No_leftright
        if (leftright) {
            // Use Steele & White method of only generating digits needed.
            eps = 0.5/tens[ilim-1] - eps;
            for (i = 0;;) {
                L = (int32)d;
                d -= L;
                *s++ = char('0' + L);
                if (d < eps)
                    goto ret1;
                if (1. - d < eps)
                    goto bump_up;
                if (++i >= ilim)
                    break;
                eps *= 10.;
                d *= 10.;
            }
        }
        else {
#endif
            // Generate ilim digits, then fix them up.
            eps *= tens[ilim-1];
            for (i = 1;; i++, d *= 10.) {
                L = (int32)d;
                d -= L;
                *s++ = char('0' + L);
                if (i == ilim) {
                    if (d > 0.5 + eps)
                        goto bump_up;
                    else if (d < 0.5 - eps) {
                        while (*--s == '0') ;
                        s++;
                        goto ret1;
                    }
                    break;
                }
            }
#ifndef No_leftright
        }
#endif
      fast_failed:
        s = buf;
        d = d2;
        k = k0;
        ilim = ilim0;
    }

    // Do we have a "small" integer?
    if (be >= 0 && k <= Int_max) {
        // Yes.
        ds = tens[k];
        if (ndigits < 0 && ilim <= 0) {
            if (ilim < 0 || d < 5*ds || (!biasUp && d == 5*ds))
                goto no_digits;
          one_digit:
            *s++ = '1';
            k++;
            goto ret1;
        }
        for (i = 1;; i++) {
            L = (int32) (d / ds);
            d -= L*ds;
#ifdef Check_FLT_ROUNDS
            // If FLT_ROUNDS == 2, L will usually be high by 1
            if (d < 0) {
                L--;
                d += ds;
            }
#endif
            *s++ = char('0' + L);
            if (i == ilim) {
                d += d;
                if ((d > ds) || (d == ds && (L & 1 || biasUp))) {
                  bump_up:
                    while (*--s == '9')
                        if (s == buf) {
                            k++;
                            *s = '0';
                            break;
                        }
                    ++*s++;
                }
                break;
            }
            if (!(d *= 10.))
                break;
        }
        goto ret1;
    }

    {
        m2 = b2;
        m5 = b5;
        BigInt mLow;    // If spec_case is false, assume that mLow == mHigh
        BigInt mHigh;

        if (leftright) {
            if (mode < 2) {
                i =
#ifndef Sudden_Underflow
                    denorm ? be + (Bias + (P-1) - 1 + 1) :
#endif
                    1 + P - bbits;
            // i is 1 plus the number of trailing zero bits in d's significand. Thus,
            // (2^m2 * 5^m5) / (2^(s2+i) * 5^s5) = (1/2 lsb of d)/10^k.
            }
            else {
                j = ilim - 1;
                if (m5 >= j)
                    m5 -= j;
                else {
                    s5 += j -= m5;
                    b5 += j;
                    m5 = 0;
                }
                if ((i = ilim) < 0) {
                    m2 -= i;
                    i = 0;
                }
                // (2^m2 * 5^m5) / (2^(s2+i) * 5^s5) = (1/2 * 10^(1-ilim))/10^k.
            }
            b2 += i;
            s2 += i;
            mHigh.init(1);
            // (mHigh * 2^m2 * 5^m5) / (2^s2 * 5^s5) = one-half of last printed (when mode >= 2) or
            // input (when mode < 2) significant digit, divided by 10^k.
        }
        // We still have d/10^k = (b * 2^b2 * 5^b5) / (2^s2 * 5^s5).  Reduce common factors in
        // b2, m2, and s2 without changing the equalities.
        if (m2 > 0 && s2 > 0) {
            i = m2 < s2 ? m2 : s2;
            b2 -= i;
            m2 -= i;
            s2 -= i;
        }

        // Fold b5 into b and m5 into mHigh.
        if (b5 > 0) {
            if (leftright) {
                if (m5 > 0) {
                    mHigh.pow5Mul(m5);
                    b *= mHigh;
                }
                if ((j = b5 - m5) != 0)
                    b.pow5Mul(j);
            }
            else
                b.pow5Mul(b5);
        }
        // Now we have d/10^k = (b * 2^b2) / (2^s2 * 5^s5) and
        // (mHigh * 2^m2) / (2^s2 * 5^s5) = one-half of last printed or input significant digit, divided by 10^k.

        BigInt S;
        S.init(1);
        if (s5 > 0)
            S.pow5Mul(s5);
        // Now we have d/10^k = (b * 2^b2) / (S * 2^s2) and
        // (mHigh * 2^m2) / (S * 2^s2) = one-half of last printed or input significant digit, divided by 10^k.

        // Check for special case that d is a normalized power of 2.
        bool spec_case = false;
        if (mode < 2) {
            if (!word1(d) && !(word0(d) & Bndry_mask)
#ifndef Sudden_Underflow
                && word0(d) & (Exp_mask & Exp_mask << 1)
#endif
                ) {
                // The special case.  Here we want to be within a quarter of the last input
                // significant digit instead of one half of it when the decimal output string's value is less than d.
                b2 += Log2P;
                s2 += Log2P;
                spec_case = true;
            }
        }

        // Arrange for convenient computation of quotients:
        // shift left if necessary so divisor has 4 leading 0 bits.
        //
        // Perhaps we should just compute leading 28 bits of S once
        // and for all and pass them and a shift to quoRem, so it
        // can do shifts and ors to compute the numerator for q.
        if ((i = ((s5 ? 32 - hi0bits(S.word(S.nWords()-1)) : 1) + s2) & 0x1f) != 0)
            i = 32 - i;
        // i is the number of leading zero bits in the most significant word of S*2^s2.
        if (i > 4) {
            i -= 4;
            b2 += i;
            m2 += i;
            s2 += i;
        }
        else if (i < 4) {
            i += 28;
            b2 += i;
            m2 += i;
            s2 += i;
        }
        // Now S*2^s2 has exactly four leading zero bits in its most significant word.
        if (b2 > 0)
            b.pow2Mul(b2);
        if (s2 > 0)
            S.pow2Mul(s2);
        // Now we have d/10^k = b/S and
        // (mHigh * 2^m2) / S = maximum acceptable error, divided by 10^k.
        if (k_check) {
            if (b.cmp(S) < 0) {
                k--;
                b.mulAdd(10, 0);  // we botched the k estimate
                if (leftright)
                    mHigh.mulAdd(10, 0);
                ilim = ilim1;
            }
        }
        // At this point 1 <= d/10^k = b/S < 10.

        if (ilim <= 0 && mode > 2) {
            // We're doing fixed-mode output and d is less than the minimum nonzero output in this mode.
            // Output either zero or the minimum nonzero output depending on which is closer to d.
            // Always emit at least one digit.  If the number appears to be zero
            // using the current mode, then emit one '0' digit and set decpt to 1.
            if (ilim < 0)
                goto no_digits;
            S.mulAdd(5, 0);
            if ((i = b.cmp(S)) < 0 || (i == 0 && !biasUp))
                goto no_digits;
            goto one_digit;
        }
        if (leftright) {
            if (m2 > 0)
                mHigh.pow2Mul(m2);

            // Compute mLow -- check for special case that d is a normalized power of 2.
            if (spec_case) {
                mLow = mHigh;
                mHigh.pow2Mul(Log2P);
            }
            // mLow/S = maximum acceptable error, divided by 10^k, if the output is less than d.
            // mHigh/S = maximum acceptable error, divided by 10^k, if the output is greater than d.

            for (i = 1;;i++) {
                dig = b.quoRem(S) + '0';
                // Do we yet have the shortest decimal string that will round to d?
                j = b.cmp(spec_case ? mLow : mHigh);        // j is b/S compared with mLow/S.
                {
                    BigInt delta;
                    delta.initDiff(S, mHigh);
                    j1 = delta.negative ? 1 : b.cmp(delta);
                }
                // j1 is b/S compared with 1 - mHigh/S.
#ifndef ROUND_BIASED
                if (j1 == 0 && !mode && !(word1(d) & 1)) {
                    if (dig == '9')
                        goto round_9_up;
                    if (j > 0)
                        dig++;
                    *s++ = (char)dig;
                    goto ret;
                }
#endif
                if ((j < 0) || (j == 0 && !mode
#ifndef ROUND_BIASED
                                && !(word1(d) & 1)
#endif
                                )) {
                    if (j1 > 0) {
                        // Either dig or dig+1 would work here as the least significant decimal digit.
                        // Use whichever would produce a decimal value closer to d.
                        b.pow2Mul(1);
                        j1 = b.cmp(S);
                        if (((j1 > 0) || (j1 == 0 && (dig & 1 || biasUp))) && (dig++ == '9'))
                            goto round_9_up;
                    }
                    *s++ = (char)dig;
                    goto ret;
                }
                if (j1 > 0) {
                    if (dig == '9') { // possible if i == 1
                      round_9_up:
                        *s++ = '9';
                        goto roundoff;
                    }
                    *s++ = char(dig + 1);
                    goto ret;
                }
                *s++ = (char)dig;
                if (i == ilim)
                    break;
                b.mulAdd(10, 0);
                if (spec_case)
                    mLow.mulAdd(10, 0);
                mHigh.mulAdd(10, 0);
            }
        }
        else
            for (i = 1;; i++) {
                *s++ = (char)(dig = b.quoRem(S) + '0');
                if (i >= ilim)
                    break;
                b.mulAdd(10, 0);
            }

        // Round off last digit
        b.pow2Mul(1);
        j = b.cmp(S);
        if ((j > 0) || (j == 0 && (dig & 1 || biasUp))) {
          roundoff:
            while (*--s == '9')
                if (s == buf) {
                    k++;
                    *s++ = '1';
                    goto ret;
                }
            ++*s++;
        }
        else {
            // Strip trailing zeros
            while (*--s == '0') ;
            s++;
        }
      ret: ;
        // S, mLow, and mHigh are destroyed at the end of this block scope.
    }
  ret1:
    *s = '\0';
    *decpt = k + 1;
    return s;
}


// Mapping of DToStrMode -> doubleToAscii mode
static const int doubleToAsciiModes[] = {
    0,   // dtosStandard
    0,   // dtosStandardExponential
    3,   // dtosFixed
    2,   // dtosExponential
    2};  // dtosPrecision


// Convert value according to the given mode and return a pointer to the resulting ASCII string.
// The result is held somewhere in buffer, but not necessarily at the beginning.  The size of
// buffer is given in bufferSize, and must be at least as large as given by dtosStandardBufferSize
// or dtosVariableBufferSize, whichever is appropriate.
char *JS::doubleToStr(char *buffer, size_t DEBUG_ONLY(bufferSize), double value, DToStrMode mode, int precision)
{
    int decPt;                  // Position of decimal point relative to first digit returned by doubleToAscii
    bool negative;              // True if the sign bit was set in value
    int nDigits;                // Number of significand digits returned by doubleToAscii
    char *numBegin = buffer+2;  // Pointer to the digits returned by doubleToAscii; the +2 leaves space for
                                // the sign and/or decimal point
    char *numEnd;               // Pointer past the digits returned by doubleToAscii

    ASSERT(bufferSize >= (size_t)(mode <= dtosStandardExponential ? dtosStandardBufferSize : dtosVariableBufferSize(precision)));

    if (mode == dtosFixed && (value >= 1e21 || value <= -1e21))
        mode = dtosStandard;    // Change mode here rather than below because the buffer may not be large enough to hold a large integer.

    numEnd = doubleToAscii(value, doubleToAsciiModes[mode], mode >= dtosFixed, precision, &decPt, &negative, numBegin);
    nDigits = numEnd - numBegin;

    // If Infinity, -Infinity, or NaN, return the string regardless of the mode.
    if (decPt != 9999) {
        bool exponentialNotation = false;
        int minNDigits = 0;         // Minimum number of significand digits required by mode and precision
        char *p;
        char *q;

        switch (mode) {
            case dtosStandard:
                if (decPt < -5 || decPt > 21)
                    exponentialNotation = true;
                else
                    minNDigits = decPt;
                break;

            case dtosFixed:
                if (precision >= 0)
                    minNDigits = decPt + precision;
                else
                    minNDigits = decPt;
                break;

            case dtosExponential:
                ASSERT(precision > 0);
                minNDigits = precision;
                // Fall through
            case dtosStandardExponential:
                exponentialNotation = true;
                break;

            case dtosPrecision:
                ASSERT(precision > 0);
                minNDigits = precision;
                if (decPt < -5 || decPt > precision)
                    exponentialNotation = true;
                break;
        }

        // If the number has fewer than minNDigits, pad it with zeros at the end
        if (nDigits < minNDigits) {
            p = numBegin + minNDigits;
            nDigits = minNDigits;
            do {
                *numEnd++ = '0';
            } while (numEnd != p);
            *numEnd = '\0';
        }

        if (exponentialNotation) {
            // Insert a decimal point if more than one significand digit
            if (nDigits != 1) {
                numBegin--;
                numBegin[0] = numBegin[1];
                numBegin[1] = '.';
            }
            sprintf(numEnd, "e%+d", decPt-1);
        } else if (decPt != nDigits) {
            // Some kind of a fraction in fixed notation
            ASSERT(decPt <= nDigits);
            if (decPt > 0) {
                // dd...dd . dd...dd
                p = --numBegin;
                do {
                    *p = p[1];
                    p++;
                } while (--decPt);
                *p = '.';
            } else {
                // 0 . 00...00dd...dd
                p = numEnd;
                numEnd += 1 - decPt;
                q = numEnd;
                ASSERT(numEnd < buffer + bufferSize);
                *numEnd = '\0';
                while (p != numBegin)
                    *--q = *--p;
                for (p = numBegin + 1; p != q; p++)
                    *p = '0';
                *numBegin = '.';
                *--numBegin = '0';
            }
        }
    }

    // If negative and neither -0.0 nor NaN, output a leading '-'.
    if (negative &&
        !(word0(value) == Sign_bit && word1(value) == 0) &&
        !((word0(value) & Exp_mask) == Exp_mask &&
          (word1(value) || (word0(value) & Frac_mask)))) {
        *--numBegin = '-';
    }
    return numBegin;
}


inline char baseDigit(uint32 digit)
{
    if (digit >= 10)
        return char('a' - 10 + digit);
    else
        return char('0' + digit);
}


// Convert value to a string in the given base.  The integral part of value will be printed exactly
// in that base, regardless of how large it is, because there is no exponential notation for non-base-ten
// numbers.  The fractional part will be rounded to as few digits as possible while still preserving
// the round-trip property (analogous to that of printing decimal numbers).  In other words, if one were
// to read the resulting string in via a hypothetical base-number-reading routine that rounds to the nearest
// IEEE double (and to an even significand if there are two equally near doubles), then the result would
// equal value (except for -0.0, which converts to "0", and NaN, which is not equal to itself).
//
// Store the result in the given buffer, which must have at least dtobasesBufferSize bytes.
// Return the number of characters stored.
size_t JS::doubleToBaseStr(char *buffer, double value, uint base)
{
    ASSERT(base >= 2 && base <= 36);

    char *p = buffer;       // Pointer to current position in the buffer
    if (value < 0.0
#ifdef _WIN32
        && !((word0(value) & Exp_mask) == Exp_mask && ((word0(value) & Frac_mask) || word1(value))) // Visual C++ doesn't know how to compare against NaN
#endif
        ) {
        *p++ = '-';
        value = -value;
    }

    // Check for Infinity and NaN
    if ((word0(value) & Exp_mask) == Exp_mask) {
        strcpy(p, !word1(value) && !(word0(value) & Frac_mask) ? "Infinity" : "NaN");
        return strlen(buffer);
    }

    // Output the integer part of value with the digits in reverse order.
    char *pInt = p;                 // Pointer to the beginning of the integer part of the string
    double valueInt = fd::floor(value); // value truncated to an integer
    uint32 digit;
    if (valueInt <= 4294967295.0) {
        uint32 n = (uint32)valueInt;
        if (n)
            do {
                uint32 m = n / base;
                digit = n - m*base;
                n = m;
                ASSERT(digit < base);
                *p++ = baseDigit(digit);
            } while (n);
        else *p++ = '0';
    } else {
        int32 e;
        int32 bits;  // Number of significant bits in valueInt; not used.
        BigInt b;
        b.init(valueInt, e, bits);
        b.pow2Mul(e);
        do {
            digit = b.divRem(base);
            ASSERT(digit < base);
            *p++ = baseDigit(digit);
        } while (!b.isZero());
    }
    // Reverse the digits of the integer part of value.
    char *q = p-1;
    while (q > pInt) {
        char ch = *pInt;
        *pInt++ = *q;
        *q-- = ch;
    }

    double valueFrac = value - valueInt;        // The fractional part of value
    if (valueFrac != 0.0) {
        // We have a fraction.
        *p++ = '.';

        int32 e, bbits;
        BigInt b;
        b.init(valueFrac, e, bbits);
        ASSERT(e < 0);
        // At this point valueFrac = b * 2^e.  e must be less than zero because 0 < valueFrac < 1.

        int32 s2 = -int32(word0(value) >> Exp_shift1 & Exp_mask>>Exp_shift1);
#ifndef Sudden_Underflow
        if (!s2)
            s2 = -1;
#endif
        s2 += Bias + P;
        // 1/2^s2 = (nextDouble(value) - value)/2
        ASSERT(-s2 < e);

        BigInt mLow;
        BigInt mHigh;
        bool useMHigh = false;  // If false, assume that mHigh == mLow

        mLow.init(1);
        if (!word1(value) && !(word0(value) & Bndry_mask)
#ifndef Sudden_Underflow
            && word0(value) & (Exp_mask & Exp_mask << 1)
#endif
            ) {
            // The special case.  Here we want to be within a quarter of the last input
            // significant digit instead of one half of it when the output string's value is less than value.
            s2 += Log2P;
            useMHigh = true;
            mHigh.init(1u<<Log2P);
        }
        b.pow2Mul(e + s2);

        BigInt s;
        s.init(1);
        s.pow2Mul(s2);

        // At this point we have the following:
        //   s = 2^s2;
        //   1 > valueFrac = b/2^s2 > 0;
        //   (value - prevDouble(value))/2 = mLow/2^s2;
        //   (nextDouble(value) - value)/2 = mHigh/2^s2.

        bool done = false;
        do {
            int32 j, j1;

            b.mulAdd(base, 0);
            digit = b.quoRem2(s2);
            mLow.mulAdd(base, 0);
            if (useMHigh)
                mHigh.mulAdd(base, 0);

            // Do we yet have the shortest string that will round to value?
            j = b.cmp(mLow);
            // j is b/2^s2 compared with mLow/2^s2.
            {
                BigInt delta;
                delta.initDiff(s, useMHigh ? mHigh : mLow);
                j1 = delta.negative ? 1 : b.cmp(delta);
            }
            // j1 is b/2^s2 compared with 1 - mHigh/2^s2.

#ifndef ROUND_BIASED
            if (j1 == 0 && !(word1(value) & 1)) {
                if (j > 0)
                    digit++;
                done = true;
            } else
#endif
                if (j < 0 || (j == 0
#ifndef ROUND_BIASED
                              && !(word1(value) & 1)
#endif
                              )) {
                    if (j1 > 0) {
                        // Either dig or dig+1 would work here as the least significant digit.
                        // Use whichever would produce an output value closer to value.
                        b.pow2Mul(1);
                        j1 = b.cmp(s);
                        if (j1 > 0)
                            // The even test (|| (j1 == 0 && (digit & 1))) is not here because it messes up odd base output
                            // such as 3.5 in base 3.
                            digit++;
                    }
                    done = true;
                } else if (j1 > 0) {
                    digit++;
                    done = true;
                }
            ASSERT(digit < base);
            *p++ = baseDigit(digit);
        } while (!done);
    }
    ASSERT(p < buffer + dtobasesBufferSize);
    *p = '\0';
    return size_t(p - buffer);
}


// A version of doubleToStr that appends to the end of String dst.
// precision should not exceed 101.
void JS::appendDouble(String &dst, double value, DToStrMode mode, int precision)
{
    char buffer[dtosVariableBufferSize(101)];
    ASSERT(uint(precision) <= 101);

    dst += doubleToStr(buffer, sizeof buffer, value, mode, precision);
}


// A version of doubleToStr that prints to Formatter f.
// precision should not exceed 101.
void JS::printDouble(Formatter &f, double value, DToStrMode mode, int precision)
{
    char buffer[dtosVariableBufferSize(101)];
    ASSERT(uint(precision) <= 101);

    f << doubleToStr(buffer, sizeof buffer, value, mode, precision);
}