Tiny float to ascii conversion (with lossless roundtrip, based on the Ryu algorithm)

nanof2a.cpp

// c++ -O3 -o nanof2a nanof2a.cpp && ./nanof2a
#include <cstdint>
#include <cstring>

namespace { // anonymous namespace to encourage inlining

struct f2a {
    char str[16];

    f2a(float f) { // example usage: printf("%s", f2a(f).str)
        // Based on https://github.com/ulfjack/ryu Copyright 2018 Ulf Adams
        // Added decimal format support to match double-conversion library output
        // The contents of this file may be used under the terms of the Apache License,
        // Version 2.0. (See copy at http://www.apache.org/licenses/LICENSE-2.0)
        // Unless required by applicable law or agreed to in writing, this software
        // is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
        // KIND, either express or implied.
        static const uint32_t POW5[10] = { 1, 5, 25, 125, 625, 3125, 15625, 78125, 390625, 1953125 }; static const uint64_t TABLE[78] = {
        576460752303423489u , 461168601842738791u , 368934881474191033u , 295147905179352826u , 472236648286964522u , 377789318629571618u ,
        302231454903657294u , 483570327845851670u , 386856262276681336u , 309485009821345069u , 495176015714152110u , 396140812571321688u ,
        316912650057057351u , 507060240091291761u , 405648192073033409u , 324518553658426727u , 519229685853482763u , 415383748682786211u ,
        332306998946228969u , 531691198313966350u , 425352958651173080u , 340282366920938464u , 544451787073501542u , 435561429658801234u ,
        348449143727040987u , 557518629963265579u , 446014903970612463u , 356811923176489971u , 570899077082383953u , 456719261665907162u ,
        365375409332725730u , 1152921504606846976u, 1441151880758558720u, 1801439850948198400u, 2251799813685248000u, 1407374883553280000u,
        1759218604441600000u, 2199023255552000000u, 1374389534720000000u, 1717986918400000000u, 2147483648000000000u, 1342177280000000000u,
        1677721600000000000u, 2097152000000000000u, 1310720000000000000u, 1638400000000000000u, 2048000000000000000u, 1280000000000000000u,
        1600000000000000000u, 2000000000000000000u, 1250000000000000000u, 1562500000000000000u, 1953125000000000000u, 1220703125000000000u,
        1525878906250000000u, 1907348632812500000u, 1192092895507812500u, 1490116119384765625u, 1862645149230957031u, 1164153218269348144u,
        1455191522836685180u, 1818989403545856475u, 2273736754432320594u, 1421085471520200371u, 1776356839400250464u, 2220446049250313080u,
        1387778780781445675u, 1734723475976807094u, 2168404344971008868u, 1355252715606880542u, 1694065894508600678u, 2117582368135750847u,
        1323488980084844279u, 1654361225106055349u, 2067951531382569187u, 1292469707114105741u, 1615587133892632177u, 2019483917365790221u };
        uint32_t bits; memcpy(&bits, &f, sizeof(float));
        char* p = str; *p = '-'; p += bits >> 31;
        const uint32_t ieeeMantissa = bits & 0x7fffff, ieeeExponent = (bits >> 23) & 255;
        if (ieeeExponent == 0   && ieeeMantissa == 0) { strcpy(p, "0");                          return; }
        if (ieeeExponent == 255                     ) { strcpy(p, ieeeMantissa ? "nan" : "inf"); return; }
        const  int32_t e2 = (ieeeExponent ? ieeeExponent : 1) - 152;
        const uint32_t m2 = (ieeeExponent ? 0x800000u : 0) | ieeeMantissa, mms = ieeeMantissa != 0 || ieeeExponent <= 1;
        const uint32_t mv = 4 * m2, mp = mv + 2, mm = mv - 1 - mms;
        uint32_t vr, vp, vm, output;
        int32_t e10, lastRemovedDigit = 0;
        bool vmIsTrailingZeros = false, vrIsTrailingZeros = false;
        if (e2 >= 0) {
            const int32_t q = e10 = (e2 * 78913) >> 18, k = 59 + ((q * 1217359) >> 19), i = -e2 + q + k;
            vr = mulShift(mv, TABLE[q], i);
            vp = mulShift(mp, TABLE[q], i);
            vm = mulShift(mm, TABLE[q], i);
            if (q != 0 && (vp - 1) / 10 <= vm / 10)
                lastRemovedDigit = mulShift(mv, TABLE[q - 1], -e2 + q + 58 + (((q - 1) * 1217359) >> 19)) % 10;
            if (q <= 9) {
                if (mv % 5 == 0)
                    vrIsTrailingZeros = mv % POW5[q] == 0;
                else if ((m2 & 1) == 0)
                    vmIsTrailingZeros = mm % POW5[q] == 0;
                else
                    vp -= mp % POW5[q] == 0;
            }
        } else {
            const uint32_t q = (uint32_t(-e2) * 732923) >> 20; e10 = q + e2;
            const int32_t i = -e2 - q, k = ((uint32_t(i) * 1217359) >> 19) - 60, j = q - k;
            vr = mulShift(mv, TABLE[i + 31], j);
            vp = mulShift(mp, TABLE[i + 31], j);
            vm = mulShift(mm, TABLE[i + 31], j);
            if (q != 0 && (vp - 1) / 10 <= vm / 10)
                lastRemovedDigit = mulShift(mv, TABLE[i + 32], q - 1 - ((((uint32_t(i + 1)) * 1217359) >> 19) - 60)) % 10;
            if (q <= 1) {
                vrIsTrailingZeros = true;
                vmIsTrailingZeros = (m2 & 1) == 0 && mms == 1;
                vp -= m2 & 1;
            } else if (q < 31) {
                vrIsTrailingZeros = (mv & ((1u << (q - 1)) - 1)) == 0;
            }
        }
        if (vmIsTrailingZeros || vrIsTrailingZeros) {
            while (vp / 10 > vm / 10) {
                vmIsTrailingZeros &= vm % 10 == 0;
                vrIsTrailingZeros &= lastRemovedDigit == 0;
                lastRemovedDigit = vr % 10;
                vr /= 10;
                vp /= 10;
                vm /= 10;
                e10++;
            }
            if (vmIsTrailingZeros) {
                while (vm % 10 == 0) {
                    vrIsTrailingZeros &= lastRemovedDigit == 0;
                    lastRemovedDigit = vr % 10;
                    vr /= 10;
                    vp /= 10;
                    vm /= 10;
                    e10++;
                }
            }
            if (vrIsTrailingZeros && lastRemovedDigit == 5 && vr % 2 == 0)
                lastRemovedDigit = 4;
            output = vr + ((vr == vm && ((m2 & 1) || !vmIsTrailingZeros)) || lastRemovedDigit >= 5);
        } else {
            while (vp / 10 > vm / 10) {
                lastRemovedDigit = vr % 10;
                vr /= 10;
                vp /= 10;
                vm /= 10;
                e10++;
            }
            output = vr + (vr == vm || lastRemovedDigit >= 5);
        }
        const int ndig = (output >= 100000000) ? 9 : (output >= 10000000) ? 8 : (output >= 1000000) ? 7 :
                         (output >= 100000)    ? 6 : (output >= 10000)    ? 5 : (output >= 1000)    ? 4 :
                         (output >= 100  )     ? 3 : (output >= 10)       ? 2 : 1, exp = e10 + ndig - 1;
        if (exp < -4 || exp > 8) { // exponential form
            for (int i = 0; i < ndig - 1; ++i, output /= 10) p[ndig - i] = '0' + output % 10;
            *p++ = '0' + output; if (ndig > 1) { *p = '.'; p += ndig; } *p++ = 'e';
            if      (exp < -9) { *p++ = '-'; *p++ = '0' - exp / 10; *p++ = '0' - exp % 10; }
            else if (exp < -4) { *p++ = '-'; *p++ = '0' - exp; }
            else if (exp >  9) {             *p++ = '0' + exp / 10; *p++ = '0' + exp % 10; }
            else               {             *p++ = '9'; }
        } else { // decimal form
            const int nint = ndig + e10, ndec = -e10; char d[9];
            for (int i = ndig; i--; output /= 10) d[i] = '0' + output % 10;
            if (nint > 0) {
                for (int i = 0; i < nint; i++) *p++ = i < ndig ? d[i] : '0';
                *p = '.'; p += ndec > 0;
                for (int i = 0; i < ndec; i++) *p++ = d[i + nint];
            } else {
                *p++ = '0'; *p++ = '.';
                for (int i = 0; i < -nint; i++) *p++ = '0';
                for (int i = 0; i <  ndig; i++) *p++ = d[i];
            }
        }
        *p = 0;
    }
private:
    static inline uint32_t mulShift(uint64_t m, uint64_t f, int s) {
        return (((m * uint32_t(f)) >> 32) + (m * uint32_t(f >> 32))) >> (s - 32);
    }
};

} // anonymous namespace

// Demo / unit-test:
#include <cmath>
#include <cstdlib>
#include <cstdio>
#include <limits>

int main() {
    float a[] = {
        0.0f, -0.0f, 1.0f, 0.1f, 0.12f, 0.123f, 0.1234f, 0.12345f, 0.123456f, 0.1234567f,
        0.12345678f, 0.123456789f, 3.0540412e5f, 8.0990312e3f, 3.4028235e38f, 1e-45f,
        2.4414062e-4f, 2.4414062e-3f, 4.3945312e-3f, 6.3476562e-3f, 6.7108864e17f,
        std::numeric_limits<float>::infinity(),
        -std::numeric_limits<float>::infinity(),
        std::numeric_limits<float>::quiet_NaN(),
        std::numeric_limits<float>::signaling_NaN(),
    };
    for (int i = 0, n = sizeof(a) / sizeof(float); i < n; i++)
        printf("%-14.9g -> %s\n", a[i], f2a(a[i]).str);
    printf("checking all positive bit patterns:\n");
    for (uint32_t i = 0; i < 0x7f800000; i++) {
        float f; memcpy(&f, &i, sizeof(float));
        float p = strtof(f2a(f).str, nullptr);
        if (f != p)
            return printf("float %.9g parsed back to %.9g\n", f, p);
        if ((i & 0x3FFFFF) == 0)
            printf("\r%5.1f%%", 100.0 * i / 0x7f800000), fflush(stdout);
    }
    printf("\rall ok!\n");
    return 0;
}