IEEE754 double to string conversion function upto 3 precision bits for denormals

f64toa.c

static uint64_t ftoa(double val, char* buf, uint32_t precision)
{
    // double IEEE-745: [sign:1] [exponent:11] [mantissa:52]
    // bias = 1023, all exponents have bias to store as positive numbers
    // normalized => exponent != 0, subnormal => exponent == 0
    uint64_t bytes = 0;
    union { double d; uint64_t u; } v = { val };
    uint64_t bits = v.u;

    int64_t exponent = (bits >> 52) & (0x7FF); // remove 52 bits and get the trailing 11 bits
    uint64_t mantissa = bits & 0xFFFFFFFFFFFFF; // get the trailing 52-bits
    bool is_denormal = false;

    if (exponent == 0) {
        if (mantissa == 0) {
            *buf++ = '0';
            *buf++ = '.';
            bytes += 2;
            for (uint32_t i = 0; i < precision; i++) {
                *buf++ = '0';
                bytes++;
            }
            return bytes;
        }
        else {
            is_denormal = true;
            exponent = 1LLU - DOUBLE_EXP_CENTER_BIAS;
        }
    } else if (exponent == 0x7FF) {
        // special cases
        if (mantissa == 0) {
            //  INF
            memcpy(buf, "inf", 3 * sizeof(char));
        }
        else {
            // NAN
            memcpy(buf, "nan", 3 * sizeof(char));
        }
        bytes += 3;
        return bytes;
    } else {
        // Normalized
        exponent = exponent - 1023LLU;  // unbiased
    }

    uint64_t int_part = 0;
    uint64_t frac_part = 0;
    uint64_t shift = 0;
    if (exponent < 0) {
        // All bits are fractional
        int_part = 0;
        shift = 52 + (-exponent);
        if (is_denormal) {
            frac_part = mantissa;
        }
        else {
            frac_part = (1ULL << 52) | mantissa;
        }
    }
    else if (exponent >= 52) {
        // All bits contribute to integer part (1 + mantissa)
        int_part = (1ULL << exponent) | (mantissa << (exponent - 52));
        frac_part = 0;
    }
    else {
        // Split between integer and fraction
        shift = 52 - exponent;
        // get top int bits (1 + mantissa)
        int_part = (1ULL << exponent) | (mantissa >> shift);
        // get lower fract part
        frac_part = mantissa & ((1ULL << shift) - 1);
    }

    bytes += itoa(int_part, buf, 10);
    buf += bytes;
    *buf++ = '.';
    bytes++;
    // convert fract to string
    for (uint32_t i = 0; i < precision; i++) {
        frac_part *= 10;
        uint32_t digit = frac_part >> shift;
        *buf++ = "0123456789"[digit % 10];
        bytes++;
        frac_part &= ((1ULL << shift) - 1);
    }

    return bytes;
}