int js_dtoa()

in src/couch_quickjs/quickjs/dtoa.c [1108:1318]


int js_dtoa(char *buf, double d, int radix, int n_digits, int flags,
            JSDTOATempMem *tmp_mem)
{
    uint64_t a, m, *mptr = tmp_mem->mem;
    int e, sgn, l, E, P, i, E_max, radix1, radix_shift;
    char *q;
    mpb_t *tmp1, *mant_max;
    int fmt = flags & JS_DTOA_FORMAT_MASK;

    tmp1 = dtoa_malloc(&mptr, sizeof(mpb_t) + sizeof(limb_t) * DBIGNUM_LEN_MAX);
    mant_max = dtoa_malloc(&mptr, sizeof(mpb_t) + sizeof(limb_t) * MANT_LEN_MAX);
    assert((mptr - tmp_mem->mem) <= sizeof(JSDTOATempMem) / sizeof(mptr[0]));

    radix_shift = ctz32(radix);
    radix1 = radix >> radix_shift;
    a = float64_as_uint64(d);
    sgn = a >> 63;
    e = (a >> 52) & 0x7ff;
    m = a & (((uint64_t)1 << 52) - 1);
    q = buf;
    if (e == 0x7ff) {
        if (m == 0) {
            if (sgn)
                *q++ = '-';
            memcpy(q, "Infinity", 8);
            q += 8;
        } else {
            memcpy(q, "NaN", 3);
            q += 3;
        }
        goto done;
    } else if (e == 0) {
        if (m == 0) {
            tmp1->len = 1;
            tmp1->tab[0] = 0;
            E = 1;
            if (fmt == JS_DTOA_FORMAT_FREE)
                P = 1;
            else if (fmt == JS_DTOA_FORMAT_FRAC)
                P = n_digits + 1;
            else
                P = n_digits;
            /* "-0" is displayed as "0" if JS_DTOA_MINUS_ZERO is not present */
            if (sgn && (flags & JS_DTOA_MINUS_ZERO))
                *q++ = '-';
            goto output;
        }
        /* denormal number: convert to a normal number */
        l = clz64(m) - 11;
        e -= l - 1;
        m <<= l;
    } else {
        m |= (uint64_t)1 << 52;
    }
    if (sgn)
        *q++ = '-';
    /* remove the bias */
    e -= 1022;
    /* d = 2^(e-53)*m */
    //    printf("m=0x%016" PRIx64 " e=%d\n", m, e);
#ifdef USE_FAST_INT
    if (fmt == JS_DTOA_FORMAT_FREE &&
        e >= 1 && e <= 53 &&
        (m & (((uint64_t)1 << (53 - e)) - 1)) == 0 &&
        (flags & JS_DTOA_EXP_MASK) != JS_DTOA_EXP_ENABLED) {
        m >>= 53 - e;
        /* 'm' is never zero */
        q += u64toa_radix(q, m, radix);
        goto done;
    }
#endif
    
    /* this choice of E implies F=round(x*B^(P-E) is such as: 
       B^(P-1) <= F < 2.B^P. */
    E = 1 + mul_log2_radix(e - 1, radix);
    
    if (fmt == JS_DTOA_FORMAT_FREE) {
        int P_max, E0, e1, E_found, P_found;
        uint64_t m1, mant_found, mant, mant_max1;
        /* P_max is guarranteed to work by construction */
        P_max = dtoa_max_digits_table[radix - 2];
        E0 = E;
        E_found = 0;
        P_found = 0;
        mant_found = 0;
        /* find the minimum number of digits by successive tries */
        P = P_max; /* P_max is guarateed to work */
        for(;;) {
            /* mant_max always fits on 64 bits */
            mant_max1 = pow_ui(radix, P);
            /* compute the mantissa in base B */
            E = E0;
            for(;;) {
                /* XXX: add inexact flag */
                mul_pow_round(tmp1, m, e - 53, radix1, radix_shift, P - E, JS_RNDN);
                mant = mpb_get_u64(tmp1);
                if (mant < mant_max1)
                    break;
                E++; /* at most one iteration is possible */
            }
            /* remove useless trailing zero digits */
            while ((mant % radix) == 0) {
                mant /= radix;
                P--;
            }
            /* garanteed to work for P = P_max */
            if (P_found == 0)
                goto prec_found;
            /* convert back to base 2 */
            mpb_set_u64(tmp1, mant);
            m1 = mul_pow_round_to_d(&e1, tmp1, radix1, radix_shift, E - P, JS_RNDN);
            //            printf("P=%2d: m=0x%016" PRIx64 " e=%d m1=0x%016" PRIx64 " e1=%d\n", P, m, e, m1, e1);
            /* Note: (m, e) is never zero here, so the exponent for m1
               = 0 does not matter */
            if (m1 == m && e1 == e) {
            prec_found:
                P_found = P;
                E_found = E;
                mant_found = mant;
                if (P == 1)
                    break;
                P--; /* try lower exponent */
            } else {
                break;
            }
        }
        P = P_found;
        E = E_found;
        mpb_set_u64(tmp1, mant_found);
#ifdef JS_DTOA_DUMP_STATS
        if (radix == 10) {
            out_len_count[P - 1]++;
        }
#endif        
    } else if (fmt == JS_DTOA_FORMAT_FRAC) {
        int len;

        assert(n_digits >= 0 && n_digits <= JS_DTOA_MAX_DIGITS);
        /* P = max_int(E, 1) + n_digits; */
        /* frac is rounded using RNDNA */
        mul_pow_round(tmp1, m, e - 53, radix1, radix_shift, n_digits, JS_RNDNA);

        /* we add one extra digit on the left and remove it if needed
           to avoid testing if the result is < radix^P */
        len = output_digits(q, tmp1, radix, max_int(E + 1, 1) + n_digits,
                            max_int(E + 1, 1));
        if (q[0] == '0' && len >= 2 && q[1] != '.') {
            len--;
            memmove(q, q + 1, len);
        }
        q += len;
        goto done;
    } else {
        int pow_shift;
        assert(n_digits >= 1 && n_digits <= JS_DTOA_MAX_DIGITS);
        P = n_digits;
        /* mant_max = radix^P */
        mant_max->len = 1;
        mant_max->tab[0] = 1;
        pow_shift = mul_pow(mant_max, radix1, radix_shift, P, FALSE, 0);
        mpb_shr_round(mant_max, pow_shift, JS_RNDZ);
        
        for(;;) {
            /* fixed and frac are rounded using RNDNA */
            mul_pow_round(tmp1, m, e - 53, radix1, radix_shift, P - E, JS_RNDNA);
            if (mpb_cmp(tmp1, mant_max) < 0)
                break;
            E++; /* at most one iteration is possible */
        }
    }
 output:
    if (fmt == JS_DTOA_FORMAT_FIXED)
        E_max = n_digits;
    else
        E_max = dtoa_max_digits_table[radix - 2] + 4;
    if ((flags & JS_DTOA_EXP_MASK) == JS_DTOA_EXP_ENABLED ||
        ((flags & JS_DTOA_EXP_MASK) == JS_DTOA_EXP_AUTO && (E <= -6 || E > E_max))) {
        q += output_digits(q, tmp1, radix, P, 1);
        E--;
        if (radix == 10) {
            *q++ = 'e';
        } else if (radix1 == 1 && radix_shift <= 4) {
            E *= radix_shift;
            *q++ = 'p';
        } else {
            *q++ = '@';
        }
        if (E < 0) {
            *q++ = '-';
            E = -E;
        } else {
            *q++ = '+';
        }
        q += u32toa(q, E);
    } else if (E <= 0) {
        *q++ = '0';
        *q++ = '.';
        for(i = 0; i < -E; i++)
            *q++ = '0';
        q += output_digits(q, tmp1, radix, P, P);
    } else {
        q += output_digits(q, tmp1, radix, P, min_int(P, E));
        for(i = 0; i < E - P; i++)
            *q++ = '0';
    }
 done:
    *q = '\0';
    dtoa_free(mant_max);
    dtoa_free(tmp1);
    return q - buf;
}