ieee754_string.cpp.hpp

Go to the documentation of this file.

// fl::ieee754_string -- integer-only IEEE 754 single-precision codec.
//
// See ``src/fl/stl/ieee754_string.h`` for the design rationale and the
// embedding into FastLED #3022 / #3029.
 
#include "fl/stl/ieee754_string.h"
 
#include "fl/stl/string.h"
#include "fl/stl/charconv.h"
 
namespace fl {
 
namespace {
 
// ------------------------------------------------------------------ //
// pow10 lookup table. Generated by ci/tools/gen_pow10_table.py.
// Each entry encodes 10**k ~= mantissa * 2**binary_exp, mantissa
// normalized so the top bit (bit 63) is set.
// ------------------------------------------------------------------ //
 
static constexpr int kPow10KMin = -64;
static constexpr int kPow10KMax = 38;
static constexpr fl::size kPow10Count = 103;
 
static constexpr fl::u64 kPow10Mant[kPow10Count] = {
    0xA87FEA27A539E9A5u,  // 10**-64
    0xD29FE4B18E88640Fu,  // 10**-63
    0x83A3EEEEF9153E89u,  // 10**-62
    0xA48CEAAAB75A8E2Bu,  // 10**-61
    0xCDB02555653131B6u,  // 10**-60
    0x808E17555F3EBF12u,  // 10**-59
    0xA0B19D2AB70E6ED6u,  // 10**-58
    0xC8DE047564D20A8Cu,  // 10**-57
    0xFB158592BE068D2Fu,  // 10**-56
    0x9CED737BB6C4183Du,  // 10**-55
    0xC428D05AA4751E4Du,  // 10**-54
    0xF53304714D9265E0u,  // 10**-53
    0x993FE2C6D07B7FACu,  // 10**-52
    0xBF8FDB78849A5F97u,  // 10**-51
    0xEF73D256A5C0F77Du,  // 10**-50
    0x95A8637627989AAEu,  // 10**-49
    0xBB127C53B17EC159u,  // 10**-48
    0xE9D71B689DDE71B0u,  // 10**-47
    0x9226712162AB070Eu,  // 10**-46
    0xB6B00D69BB55C8D1u,  // 10**-45
    0xE45C10C42A2B3B06u,  // 10**-44
    0x8EB98A7A9A5B04E3u,  // 10**-43
    0xB267ED1940F1C61Cu,  // 10**-42
    0xDF01E85F912E37A3u,  // 10**-41
    0x8B61313BBABCE2C6u,  // 10**-40
    0xAE397D8AA96C1B78u,  // 10**-39
    0xD9C7DCED53C72256u,  // 10**-38
    0x881CEA14545C7575u,  // 10**-37
    0xAA242499697392D3u,  // 10**-36
    0xD4AD2DBFC3D07788u,  // 10**-35
    0x84EC3C97DA624AB5u,  // 10**-34
    0xA6274BBDD0FADD62u,  // 10**-33
    0xCFB11EAD453994BAu,  // 10**-32
    0x81CEB32C4B43FCF5u,  // 10**-31
    0xA2425FF75E14FC32u,  // 10**-30
    0xCAD2F7F5359A3B3Eu,  // 10**-29
    0xFD87B5F28300CA0Eu,  // 10**-28
    0x9E74D1B791E07E48u,  // 10**-27
    0xC612062576589DDBu,  // 10**-26
    0xF79687AED3EEC551u,  // 10**-25
    0x9ABE14CD44753B53u,  // 10**-24
    0xC16D9A0095928A27u,  // 10**-23
    0xF1C90080BAF72CB1u,  // 10**-22
    0x971DA05074DA7BEFu,  // 10**-21
    0xBCE5086492111AEBu,  // 10**-20
    0xEC1E4A7DB69561A5u,  // 10**-19
    0x9392EE8E921D5D07u,  // 10**-18
    0xB877AA3236A4B449u,  // 10**-17
    0xE69594BEC44DE15Bu,  // 10**-16
    0x901D7CF73AB0ACD9u,  // 10**-15
    0xB424DC35095CD80Fu,  // 10**-14
    0xE12E13424BB40E13u,  // 10**-13
    0x8CBCCC096F5088CCu,  // 10**-12
    0xAFEBFF0BCB24AAFFu,  // 10**-11
    0xDBE6FECEBDEDD5BFu,  // 10**-10
    0x89705F4136B4A597u,  // 10**-9
    0xABCC77118461CEFDu,  // 10**-8
    0xD6BF94D5E57A42BCu,  // 10**-7
    0x8637BD05AF6C69B6u,  // 10**-6
    0xA7C5AC471B478423u,  // 10**-5
    0xD1B71758E219652Cu,  // 10**-4
    0x83126E978D4FDF3Bu,  // 10**-3
    0xA3D70A3D70A3D70Au,  // 10**-2
    0xCCCCCCCCCCCCCCCDu,  // 10**-1
    0x8000000000000000u,  // 10**0
    0xA000000000000000u,  // 10**1
    0xC800000000000000u,  // 10**2
    0xFA00000000000000u,  // 10**3
    0x9C40000000000000u,  // 10**4
    0xC350000000000000u,  // 10**5
    0xF424000000000000u,  // 10**6
    0x9896800000000000u,  // 10**7
    0xBEBC200000000000u,  // 10**8
    0xEE6B280000000000u,  // 10**9
    0x9502F90000000000u,  // 10**10
    0xBA43B74000000000u,  // 10**11
    0xE8D4A51000000000u,  // 10**12
    0x9184E72A00000000u,  // 10**13
    0xB5E620F480000000u,  // 10**14
    0xE35FA931A0000000u,  // 10**15
    0x8E1BC9BF04000000u,  // 10**16
    0xB1A2BC2EC5000000u,  // 10**17
    0xDE0B6B3A76400000u,  // 10**18
    0x8AC7230489E80000u,  // 10**19
    0xAD78EBC5AC620000u,  // 10**20
    0xD8D726B7177A8000u,  // 10**21
    0x878678326EAC9000u,  // 10**22
    0xA968163F0A57B400u,  // 10**23
    0xD3C21BCECCEDA100u,  // 10**24
    0x84595161401484A0u,  // 10**25
    0xA56FA5B99019A5C8u,  // 10**26
    0xCECB8F27F4200F3Au,  // 10**27
    0x813F3978F8940984u,  // 10**28
    0xA18F07D736B90BE5u,  // 10**29
    0xC9F2C9CD04674EDFu,  // 10**30
    0xFC6F7C4045812296u,  // 10**31
    0x9DC5ADA82B70B59Eu,  // 10**32
    0xC5371912364CE305u,  // 10**33
    0xF684DF56C3E01BC7u,  // 10**34
    0x9A130B963A6C115Cu,  // 10**35
    0xC097CE7BC90715B3u,  // 10**36
    0xF0BDC21ABB48DB20u,  // 10**37
    0x96769950B50D88F4u,  // 10**38
};
 
static constexpr fl::i16 kPow10BExp[kPow10Count] = {
     -276,  // 10**-64
     -273,  // 10**-63
     -269,  // 10**-62
     -266,  // 10**-61
     -263,  // 10**-60
     -259,  // 10**-59
     -256,  // 10**-58
     -253,  // 10**-57
     -250,  // 10**-56
     -246,  // 10**-55
     -243,  // 10**-54
     -240,  // 10**-53
     -236,  // 10**-52
     -233,  // 10**-51
     -230,  // 10**-50
     -226,  // 10**-49
     -223,  // 10**-48
     -220,  // 10**-47
     -216,  // 10**-46
     -213,  // 10**-45
     -210,  // 10**-44
     -206,  // 10**-43
     -203,  // 10**-42
     -200,  // 10**-41
     -196,  // 10**-40
     -193,  // 10**-39
     -190,  // 10**-38
     -186,  // 10**-37
     -183,  // 10**-36
     -180,  // 10**-35
     -176,  // 10**-34
     -173,  // 10**-33
     -170,  // 10**-32
     -166,  // 10**-31
     -163,  // 10**-30
     -160,  // 10**-29
     -157,  // 10**-28
     -153,  // 10**-27
     -150,  // 10**-26
     -147,  // 10**-25
     -143,  // 10**-24
     -140,  // 10**-23
     -137,  // 10**-22
     -133,  // 10**-21
     -130,  // 10**-20
     -127,  // 10**-19
     -123,  // 10**-18
     -120,  // 10**-17
     -117,  // 10**-16
     -113,  // 10**-15
     -110,  // 10**-14
     -107,  // 10**-13
     -103,  // 10**-12
     -100,  // 10**-11
      -97,  // 10**-10
      -93,  // 10**-9
      -90,  // 10**-8
      -87,  // 10**-7
      -83,  // 10**-6
      -80,  // 10**-5
      -77,  // 10**-4
      -73,  // 10**-3
      -70,  // 10**-2
      -67,  // 10**-1
      -63,  // 10**0
      -60,  // 10**1
      -57,  // 10**2
      -54,  // 10**3
      -50,  // 10**4
      -47,  // 10**5
      -44,  // 10**6
      -40,  // 10**7
      -37,  // 10**8
      -34,  // 10**9
      -30,  // 10**10
      -27,  // 10**11
      -24,  // 10**12
      -20,  // 10**13
      -17,  // 10**14
      -14,  // 10**15
      -10,  // 10**16
       -7,  // 10**17
       -4,  // 10**18
        0,  // 10**19
        3,  // 10**20
        6,  // 10**21
       10,  // 10**22
       13,  // 10**23
       16,  // 10**24
       20,  // 10**25
       23,  // 10**26
       26,  // 10**27
       30,  // 10**28
       33,  // 10**29
       36,  // 10**30
       39,  // 10**31
       43,  // 10**32
       46,  // 10**33
       49,  // 10**34
       53,  // 10**35
       56,  // 10**36
       59,  // 10**37
       63,  // 10**38
};
 
// Top 64 bits of the 128-bit product of two u64 values, computed via four
// 32x32->64 widening multiplies. Portable to every Cortex-M target -- uses
// only integer ops, never the libgcc soft-FP cascade.
inline fl::u64 mul_hi_u64(fl::u64 a, fl::u64 b) FL_NOEXCEPT {
    const fl::u64 ah = a >> 32;
    const fl::u64 al = a & 0xFFFFFFFFull;
    const fl::u64 bh = b >> 32;
    const fl::u64 bl = b & 0xFFFFFFFFull;
    const fl::u64 ll = al * bl;
    const fl::u64 hl = ah * bl;
    const fl::u64 lh = al * bh;
    const fl::u64 hh = ah * bh;
    const fl::u64 mid = (ll >> 32) + (hl & 0xFFFFFFFFull) + (lh & 0xFFFFFFFFull);
    return hh + (hl >> 32) + (lh >> 32) + (mid >> 32);
}
 
constexpr fl::u32 kSignBitMask = 0x80000000u;
constexpr fl::u32 kInfBitsPos  = 0x7F800000u;
constexpr fl::u32 kInfBitsNeg  = 0xFF800000u;
 
} // namespace
 
u32 ieee754_parse_decimal(const char* s, fl::size len, fl::size* consumed) FL_NOEXCEPT {
    fl::size i = 0;
 
    auto fail = [&]() -> u32 {
        if (consumed) *consumed = 0;
        return 0;
    };
 
    // Skip leading whitespace.
    while (i < len && (s[i] == ' ' || s[i] == '\t' || s[i] == '\n' ||
                       s[i] == '\r' || s[i] == '\f' || s[i] == '\v')) {
        ++i;
    }
 
    // Sign.
    u32 sign_bit = 0;
    if (i < len && s[i] == '-') { sign_bit = kSignBitMask; ++i; }
    else if (i < len && s[i] == '+') { ++i; }
 
    // Mantissa: integer + optional fractional. Cap at 19 significant decimal
    // digits -- that is the max u64 can hold (10**19 < 2**64). Beyond that we
    // drop digits but keep tracking the decimal exponent so values like
    // "100000000000000000000" (1e20) still round to the right magnitude.
    fl::u64 mantissa = 0;
    int decimal_exp = 0;
    int significant_digits = 0;
    bool seen_decimal = false;
    bool seen_digit = false;
    constexpr int kMaxSignificantDigits = 19;
 
    while (i < len) {
        const char c = s[i];
        if (c >= '0' && c <= '9') {
            seen_digit = true;
            if (significant_digits < kMaxSignificantDigits) {
                // u64 max is ~1.8e19, mantissa * 10 + 9 stays below for 18 digits.
                mantissa = mantissa * 10u + static_cast<fl::u64>(c - '0');
                if (significant_digits != 0 || c != '0') {
                    ++significant_digits;
                }
                if (seen_decimal) {
                    --decimal_exp;
                }
            } else {
                // Out of mantissa precision -- track magnitude via decimal_exp only.
                if (!seen_decimal) {
                    ++decimal_exp;
                }
            }
            ++i;
        } else if (c == '.' && !seen_decimal) {
            seen_decimal = true;
            ++i;
        } else {
            break;
        }
    }
 
    // Reject inputs that contained no digit (e.g. "", ".", "-.", "+.", "e5").
    // `i == mantissa_start` only catches "" and "-"/"+" -- a lone "." passes
    // the start-position guard because the loop consumed it.
    if (!seen_digit) {
        return fail();
    }
 
    // Optional exponent.
    if (i < len && (s[i] == 'e' || s[i] == 'E')) {
        const fl::size exp_start = i + 1;
        fl::size j = exp_start;
        int exp_sign = 1;
        if (j < len && s[j] == '-') { exp_sign = -1; ++j; }
        else if (j < len && s[j] == '+') { ++j; }
        int exp_val = 0;
        const fl::size digits_begin = j;
        while (j < len && s[j] >= '0' && s[j] <= '9') {
            if (exp_val < 10000) {
                exp_val = exp_val * 10 + (s[j] - '0');
            }
            ++j;
        }
        if (j != digits_begin) {
            decimal_exp += exp_sign * exp_val;
            i = j;
        }
        // else: lone 'e' with no digits -- treat as end-of-number, leave i at 'e'.
    }
 
    if (consumed) *consumed = i;
 
    if (mantissa == 0) {
        return sign_bit;
    }
 
    // Normalize mantissa into [2**63, 2**64) so the top bit is set. The
    // shift count becomes part of the binary exponent.
    int mantissa_shift = 0;
    while ((mantissa & 0x8000000000000000ull) == 0) {
        mantissa <<= 1;
        ++mantissa_shift;
    }
    // After normalization: original_mantissa = mantissa * 2**(-mantissa_shift).
 
    // Out-of-table decimal exponent: under/overflow before we even multiply.
    // The table floor (kPow10KMin) is chosen so the parser's raw 19-digit
    // mantissa still covers IEEE 754 single-precision range -- see
    // ci/tools/gen_pow10_table.py.
    if (decimal_exp < kPow10KMin) {
        return sign_bit;  // +/- 0
    }
    if (decimal_exp > kPow10KMax) {
        return sign_bit | kInfBitsPos;  // +/- inf
    }
 
    const fl::size idx = static_cast<fl::size>(decimal_exp - kPow10KMin);
    const fl::u64 pow_mant = kPow10Mant[idx];
    const int    pow_exp  = kPow10BExp[idx];
 
    fl::u64 product_hi = mul_hi_u64(mantissa, pow_mant);
    // The full product (mantissa * pow_mant) is in [2**126, 2**128). Its top
    // 64 bits (product_hi) are in [2**62, 2**64). If bit 63 isn't set,
    // shift left by 1 to renormalize and consume one binary-exponent unit.
    int extra_shift = 0;
    if ((product_hi & 0x8000000000000000ull) == 0) {
        product_hi <<= 1;
        extra_shift = 1;
    }
 
    // Binary exponent of the renormalized mantissa: combines the input
    // mantissa shift, the pow10 binary exponent, the +64 for taking the
    // top half of the 128-bit product, and the renormalization shift.
    int bin_exp = -mantissa_shift + pow_exp + 64 - extra_shift;
 
    // IEEE 754 single-precision biased exponent: the implicit-1 form is
    // product_hi / 2**63, so the unbiased binary exponent is (bin_exp + 63).
    int biased_exp = bin_exp + 63 + 127;
 
    if (biased_exp >= 0xFF) {
        return sign_bit | kInfBitsPos;  // +/- inf
    }
    if (biased_exp <= 0) {
        // Subnormal / underflow. For embedded-grade precision we collapse
        // to +/- 0 -- the JSON-RPC use case never needs denormals, and
        // avoiding the subnormal path keeps the bit-twiddling logic compact.
        return sign_bit;
    }
 
    // Extract the 23-bit mantissa from the top 64 bits with round-half-to-even.
    //   product_hi bits [62:40] -> ieee_mant bits [22:0]
    //   product_hi bit  [39]    -> rounding bit
    //   product_hi bits [38:0]  -> sticky bits
    fl::u32 ieee_mant = static_cast<fl::u32>((product_hi >> 40) & 0x7FFFFFu);
    const fl::u64 round_bit = (product_hi >> 39) & 1u;
    const fl::u64 sticky = (product_hi & 0x7FFFFFFFFFull);  // bits [38:0]
    if (round_bit && (sticky != 0 || (ieee_mant & 1u))) {
        ++ieee_mant;
        if (ieee_mant == (1u << 23)) {
            // Mantissa overflow on rounding -- re-renormalize.
            ieee_mant = 0;
            ++biased_exp;
            if (biased_exp >= 0xFF) {
                return sign_bit | kInfBitsPos;
            }
        }
    }
 
    return sign_bit | (static_cast<fl::u32>(biased_exp) << 23) | ieee_mant;
}
 
// Append the decimal digits of `value` to `out`. Reuses the integer-only
// `fl::utoa64` helper from `fl/stl/charconv.h` -- no FP arithmetic.
static void append_u64_decimal(fl::string& out, fl::u64 value) FL_NOEXCEPT {
    char buf[24];
    const int n = fl::utoa64(value, buf, 10);
    for (int i = 0; i < n; ++i) {
        out += fl::string(1, buf[i]);
    }
}
 
fl::string ieee754_format_decimal(u32 bits, int precision) FL_NOEXCEPT {
    if (precision < 0) precision = 0;
    if (precision > 9) precision = 9;
 
    const bool neg = (bits >> 31) & 1u;
    const int  biased_exp = static_cast<int>((bits >> 23) & 0xFFu);
    const u32  mant_bits = bits & 0x7FFFFFu;
 
    // Inf / NaN.
    if (biased_exp == 0xFF) {
        if (mant_bits != 0) return fl::string("nan");
        return neg ? fl::string("-inf") : fl::string("inf");
    }
 
    fl::string s;
 
    auto append_zero_with_precision = [&]() {
        s += "0";
        if (precision > 0) {
            s += ".";
            for (int i = 0; i < precision; ++i) s += "0";
        }
    };
 
    // +/- 0 (and any subnormal -- those collapse to zero per the parser's
    // contract, so the serializer matches).
    if (biased_exp == 0) {
        if (neg) s += "-";
        append_zero_with_precision();
        return s;
    }
 
    // Normal number: value = mantissa_full * 2**bin_exp.
    // mantissa_full carries the implicit leading 1.
    const u32 mantissa_full = mant_bits | 0x800000u;       // 24 bits
    const int bin_exp_raw = biased_exp - 127 - 23;  // signed
 
    // Lift to a 64-bit normalized representation so we can multiply against
    // the shared pow10 table.  mantissa_full's bit 23 is set, so shifting left by 40
    // puts that set bit in position 63.
    const u64 mant64 = static_cast<u64>(mantissa_full) << 40;
    const int bin_exp = bin_exp_raw - 40;
 
    // Look up the normalized representation of 10**precision.
    const fl::size idx = static_cast<fl::size>(precision - kPow10KMin);
    const u64 pow_mant = kPow10Mant[idx];
    const int pow_exp = kPow10BExp[idx];
 
    // Multiply via the same widening helper the parser uses.
    u64 scaled_hi = mul_hi_u64(mant64, pow_mant);
    int scaled_bin_exp = bin_exp + pow_exp + 64;
    if ((scaled_hi & 0x8000000000000000ull) == 0) {
        scaled_hi <<= 1;
        --scaled_bin_exp;
    }
 
    // Convert (scaled_hi, scaled_bin_exp) to a u64 integer count of
    // `precision` decimal places. Overflow on the way up clamps to +/- inf;
    // underflow on the way down collapses to +/- 0 -- same contract as the
    // parser at the edges.
    u64 scaled_int;
    if (scaled_bin_exp >= 0) {
        if (scaled_bin_exp > 0) {
            // (scaled_hi << scaled_bin_exp) drops the top bit -- magnitude
            // beyond what u64 can carry. Treat as overflow.
            return neg ? fl::string("-inf") : fl::string("inf");
        }
        scaled_int = scaled_hi;
    } else {
        const int shift = -scaled_bin_exp;
        if (shift >= 64) {
            scaled_int = 0;
        } else {
            // Round-half-to-even on the shifted-off low bits.
            const u64 mask = (shift == 64) ? ~0ull : ((1ull << shift) - 1ull);
            const u64 low = scaled_hi & mask;
            scaled_int = scaled_hi >> shift;
            const u64 half = (shift == 0) ? 0 : (1ull << (shift - 1));
            if (low > half || (low == half && (scaled_int & 1ull))) {
                ++scaled_int;
            }
        }
    }
 
    // For the integer / fractional split we need the EXACT integer value of
    // 10**precision (not the normalized form). precision is in [0, 9] so this
    // fits in u32 easily -- keep it inline to avoid pulling another table.
    static constexpr u32 kPow10IntExact[] = {
        1u,           10u,           100u,         1000u,
        10000u,       100000u,       1000000u,     10000000u,
        100000000u,   1000000000u,
    };
    const u64 pow10_exact = kPow10IntExact[precision];
 
    const u64 int_part = scaled_int / pow10_exact;
    const u64 frac_part = scaled_int % pow10_exact;
 
    if (neg && (int_part != 0 || frac_part != 0)) {
        s += "-";
    }
    append_u64_decimal(s, int_part);
 
    if (precision > 0) {
        s += ".";
        // Format `frac_part` as exactly `precision` digits, left-padded with
        // zeros. utoa64 emits the minimum-digit form, so count its length and
        // prepend the difference.
        char buf[24];
        const int n = fl::utoa64(frac_part, buf, 10);
        for (int i = n; i < precision; ++i) s += "0";
        for (int i = 0; i < n; ++i) s += fl::string(1, buf[i]);
    }
 
    return s;
}
 
} // namespace fl