pjs/js/js2/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 oqr
// 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.


#include <cstdlib>
#include <cstring>
#include <cfloat>
#include "numerics.h"

namespace JS = JavaScript;
using namespace JavaScript;


//
// 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;

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

InitNumerics::InitNumerics()
{
	word0(positiveInfinity) = Exp_mask;
	word1(positiveInfinity) = 0;
	word0(negativeInfinity) = Exp_mask | Sign_bit;
	word1(negativeInfinity) = 0;
	word0(nan) = 0x7FFFFFFF;
	word1(nan) = 0xFFFFFFFF;
}


//
// 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(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;
	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 is the 54th bit (the first dropped from the mantissa)
				bit2 = bdr.next();
				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 = 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);
}