FixedPoint.Stream.hpp

Go to the documentation of this file.

//
//  CeresEngine - A game development framework
//
//  Created by Rogiel Sulzbach.
//  Copyright (c) 2018-2023 Rogiel Sulzbach. All rights reserved.
//
 
// Based on https://archive.is/Rx9HR
 
#include "FixedPoint.hpp"
#include "LargeInteger.hpp"
 
#include <iostream>
 
namespace CeresEngine {
 
    // Returns the index of the most-signifcant set bit
    template<typename T> inline long findHighestBit(T value) noexcept {
        return sizeof(value) * 8 - 1 - LargeIntegerCLZHelper<T>::clz(value);
    }
 
    template<typename CharT, typename B, unsigned int F> std::basic_ostream<CharT>& operator<<(std::basic_ostream<CharT>& os, FixedPoint<B, F> x) noexcept {
        const auto uppercase = ((os.flags() & std::ios_base::uppercase) != 0);
        const auto showpoint = ((os.flags() & std::ios_base::showpoint) != 0);
        const auto adjustfield = (os.flags() & std::ios_base::adjustfield);
        const auto width = os.width();
        const auto& ctype = std::use_facet<std::ctype<CharT>>(os.getloc());
        const auto& numpunct = std::use_facet<std::numpunct<CharT>>(os.getloc());
 
        auto floatfield = (os.flags() & std::ios_base::floatfield);
        auto precision = os.precision();
        auto showTrailingZeros = true;
        auto useSignificantDigits = false;
 
        // Invalid precision? Reset to the default
        if(precision < 0) {
            precision = 6;
        }
 
        // Output buffer. Needs to be big enough for the formatted number without padding.
        // Optional prefixes (i.e. "+"/"-", decimal separator, exponent "e+/-" and/or "0x").
        constexpr auto worstCaseConstantSize = 6;
 
        // Maximum number of digits from the base type (covers integral + fractional digits)
        constexpr auto worstCaseDigitCount = std::numeric_limits<B>::digits10 + 2;
 
        // Exponent suffixes (i.e. maximum digits based on log of the base type size).
        // Needs a log10, but that isn't constexpr, so we're over-allocating on the stack. Can't hurt.
        constexpr auto worstCaseSuffixSize = std::numeric_limits<B>::digits;
 
        // Double the digit count: in the worst case the thousands grouping add a character per digit.
        using BufferType = std::array<CharT, worstCaseConstantSize + worstCaseDigitCount * 2 + worstCaseSuffixSize>;
        BufferType buffer;
 
        // Output cursor
        auto end = buffer.begin();
 
        // Keep track of the start of "internal" padding
        typename BufferType::iterator internalPad = buffer.end();
 
        // Representation of a number.
        // The value of the number is: raw / divisor * (10|2) ^ exponent
        // The base of the exponent is 2 in hexfloat mode, or 10 otherwise.
        struct NumberType {
            B raw; // raw fixed-point value
            B divisor; // the divisor indicating the place of the decimal point
            int exponent; // the exponent applied
        };
 
        // Convert a value without exponent to scientific representation
        // where the part before the decimal point is less than 10.
        const auto asScientific = [](NumberType value) {
            assert(value.exponent == 0);
            if(value.raw > 0) {
                while(value.raw / 10 >= value.divisor) {
                    value.divisor *= 10;
                    ++value.exponent;
                }
                while(value.raw < value.divisor) {
                    value.raw *= 10;
                    --value.exponent;
                }
            }
            return value;
        };
 
        NumberType value = {x.raw, B{1} << F, 0};
 
        auto base = B{10};
 
        // First write the sign
        if(value.raw < 0) {
            *end++ = ctype.widen('-');
            value.raw = -value.raw;
            internalPad = end;
        } else if(os.flags() & std::ios_base::showpos) {
            *end++ = ctype.widen('+');
            internalPad = end;
        }
        assert(value.raw >= B(0));
 
        switch(floatfield) {
        case std::ios_base::fixed | std::ios_base::scientific:
            // Hexadecimal mode: figure out the hexadecimal exponent and write "0x"
            if(value.raw > 0) {
                auto bit = findHighestBit(value.raw);
                value.exponent = bit - F; // exponent is applied to base 2
                value.divisor = B{1} << bit; // divisor is at the highest bit, ensuring it starts with "1."
                precision = (bit + 3) / 4; // precision is number of nibbles, so we show all of them
            }
            base = 16;
            showTrailingZeros = false; // Always strip trailing zeros in hexfloat mode
 
            *end++ = ctype.widen('0');
            *end++ = ctype.widen(uppercase ? 'X' : 'x');
            break;
 
        case std::ios_base::scientific:
            // Scientific mode, normalize value to scientific notation
            value = asScientific(value);
            break;
 
        case std::ios_base::fixed:
            // Fixed mode. Nothing to do.
            break;
 
        default: {
            // "auto" mode: figure out the exponent
            const NumberType scientificValue = asScientific(value);
 
            // Now `precision` indicates the number of *significant digits* (not fractional digits).
            useSignificantDigits = true;
            precision = std::max<std::streamsize>(precision, 1);
 
            if(scientificValue.exponent >= precision || scientificValue.exponent < -4) {
                // Display as scientific format
                floatfield = std::ios_base::scientific;
                value = scientificValue;
            } else {
                // Display as fixed format.
                // "showpoint" indicates whether or not we show trailing zeros
                floatfield = std::ios_base::fixed;
                showTrailingZeros = showpoint;
            }
            break;
        }
        };
 
        // If we didn't write a sign, any internal padding starts here
        // (after a potential "0x" for hexfloats).
        if(internalPad == buffer.end()) {
            internalPad = end;
        }
 
        // Separate out the integral part of the number
        B integral = value.raw / value.divisor;
        value.raw %= value.divisor;
 
        // Here we start printing the number itself
        const char* const digits = uppercase ? "0123456789ABCDEF" : "0123456789abcdef";
        const auto digitsStart = end;
 
        // Are we already printing significant digits? (yes if we're not counting significant digits)
        bool significantDigits = !useSignificantDigits;
 
        // Print the integral part
        B lastDigit = 0;
        if(integral == 0) {
            *end++ = ctype.widen('0');
            if(value.raw == 0) {
                // If the fraction is zero too, all zeros including the integral count
                // as significant digits.
                significantDigits = true;
            }
        } else {
            while(integral > 0) {
                lastDigit = integral % base;
                *end++ = ctype.widen(digits[(int)lastDigit]);
                integral /= base;
            }
            std::reverse(digitsStart, end);
            significantDigits = true;
        }
 
        if(useSignificantDigits && significantDigits) {
            // Apparently, the integral part was significant; subtract its
            // length from the remaining significant digits.
            precision -= (end - digitsStart);
        }
 
        // At this point, `value` contains only the fraction and
        // `precision` holds the number of digits to print.
        assert(value.raw < value.divisor);
        assert(precision >= 0);
 
        // Location of decimal point
        typename BufferType::iterator point = buffer.end();
 
        // Start (and length) of the trailing zeros to insert while printing
        // By tracking this to print them later instead of actually printing them now,
        // we can support large precisions with a small printing buffer.
        typename BufferType::iterator trailingZerosStart = buffer.end();
        std::streamsize trailingZerosCount = 0;
 
        if(precision > 0) {
            // Print the fractional part
            *(point = end++) = numpunct.decimal_point();
 
            for(int i = 0; i < precision; ++i) {
                if(value.raw == 0) {
                    // The rest of the digits are all zeros, mark them
                    // to be printed in this spot.
                    trailingZerosStart = end;
                    trailingZerosCount = precision - i;
                    break;
                }
 
                // Shift the divisor if we can to avoid overflow on the value
                if(value.divisor % base == 0) {
                    value.divisor /= base;
                } else {
                    value.raw *= base;
                }
                assert(value.divisor > 0);
                assert(value.raw >= 0);
                lastDigit = (value.raw / value.divisor) % base;
                value.raw %= value.divisor;
                *end++ = ctype.widen(digits[(int)lastDigit]);
 
                if(!significantDigits) {
                    // We're still finding the first significant digit
                    if(lastDigit != 0) {
                        // Found it
                        significantDigits = true;
                    } else {
                        // Not yet; increment number of digits to print
                        ++precision;
                    }
                }
            }
        } else if(showpoint) {
            // No fractional part to print, but we still want the point
            *(point = end++) = numpunct.decimal_point();
        }
 
        // Insert `ch` into the output at `position`, updating all references accordingly
        const auto insertCharacter = [&](typename BufferType::iterator position, CharT ch) {
            assert(position >= buffer.begin() && position < end);
            std::move_backward(position, end, end + 1);
            if(point != buffer.end() && position < point) {
                ++point;
            }
            if(trailingZerosStart != buffer.end() && position < trailingZerosStart) {
                ++trailingZerosStart;
            }
            ++end;
            *position = ch;
        };
 
        // Round the number: round to nearest
        bool increment = false;
        if(value.raw > value.divisor / 2) {
            // Round up
            increment = true;
        } else if(value.raw == value.divisor / 2) {
            // It's a tie (i.e. "xyzw.5"): round to even
            increment = ((lastDigit % 2) == 1);
        }
 
        if(increment) {
            auto p = end - 1;
            // Increment all digits backwards while we see "9"
            while(p >= digitsStart) {
                if(p == point) {
                    // Skip over the decimal point
                    --p;
                }
                if((*p)++ != ctype.widen('9')) {
                    break;
                }
                *p-- = ctype.widen('0');
            }
 
            if(p < digitsStart) {
                // We've incremented all the way to the start (all 9's), we need to insert the
                // carried-over 1 from incrementing the last 9.
                assert(p == digitsStart - 1);
                insertCharacter(++p, ctype.widen('1'));
 
                if(floatfield == std::ios::scientific) {
                    // We just made the integral part equal to 10, so we shift the decimal point
                    // back one place (if any) and tweak the exponent, so that we keep the integer part
                    // less than 10.
                    if(point != buffer.end()) {
                        assert(p + 2 == point);
                        std::swap(*(point - 1), *point);
                        --point;
                    }
                    ++value.exponent;
 
                    // We've introduced an extra digit so we need to strip the last digit
                    // to maintain the same precision
                    --end;
                }
            }
 
            if(useSignificantDigits && *p == ctype.widen('1') && point != buffer.end()) {
                // We've converted a leading zero to a 1 so we need to strip the last digit
                // (behind the decimal point) to maintain the same significant digit count.
                --end;
            }
        }
 
        if(point != buffer.end()) {
            if(!showTrailingZeros) {
                // Remove trailing zeros
                while(*(end - 1) == ctype.widen('0')) {
                    --end;
                }
 
                // Also clear the "trailing zeros to append during printing" range
                trailingZerosStart = buffer.end();
                trailingZerosCount = 0;
            }
 
            if(end - 1 == point && trailingZerosCount == 0 && !showpoint) {
                // Remove the decimal point, too
                --end;
            }
        }
 
        // Apply thousands grouping
        const auto& grouping = numpunct.grouping();
        if(!grouping.empty()) {
            // Step backwards from the end or decimal point, inserting the
            // thousands separator at every group interval.
            const CharT thousandsSeparator = ctype.widen(numpunct.thousands_sep());
            std::size_t group = 0;
            auto p = point != buffer.end() ? point : end;
            auto size = static_cast<int>(grouping[group]);
            while(size > 0 && size < CHAR_MAX && p - digitsStart > size) {
                p -= size;
                insertCharacter(p, thousandsSeparator);
                if(group < grouping.size() - 1) {
                    size = static_cast<int>(grouping[++group]);
                }
            }
        }
 
        // Print the exponent if required
        assert(floatfield != 0);
        if(floatfield & std::ios_base::scientific) {
            // Hexadecimal (%a/%A) or decimal (%e/%E) scientific notation
            if(floatfield & std::ios_base::fixed) {
                *end++ = ctype.widen(uppercase ? 'P' : 'p');
            } else {
                *end++ = ctype.widen(uppercase ? 'E' : 'e');
            }
 
            if(value.exponent < 0) {
                *end++ = ctype.widen('-');
                value.exponent = -value.exponent;
            } else {
                *end++ = ctype.widen('+');
            }
 
            if(floatfield == std::ios_base::scientific) {
                // In decimal scientific notation (%e/%E), the exponent is at least two digits
                if(value.exponent < 10) {
                    *end++ = ctype.widen('0');
                }
            }
 
            const auto exponent_start = end;
            if(value.exponent == 0) {
                *end++ = ctype.widen('0');
            } else
                while(value.exponent > 0) {
                    *end++ = ctype.widen(digits[value.exponent % 10]);
                    value.exponent /= 10;
                }
            std::reverse(exponent_start, end);
        }
 
        // Write character `ch` `count` times to the stream
        const auto putCharacter = [&](CharT ch, const std::streamsize count) {
            // Fill a buffer to output larger chunks
            constexpr std::streamsize chunkSize = 64;
            std::array<CharT, chunkSize> fillBuffer;
            std::fill_n(fillBuffer.begin(), std::min(count, chunkSize), ch);
 
            for(std::streamsize size, left = count; left > 0; left -= size) {
                size = std::min(chunkSize, left);
                os.rdbuf()->sputn(&fillBuffer[0], size);
            }
        };
 
        // Outputs a range of characters, making sure to output the trailing zeros range
        // if it lies in the specified range
        const auto putRange = [&](typename BufferType::const_iterator begin, typename BufferType::const_iterator end) {
            assert(end >= begin);
            if(trailingZerosStart >= begin && trailingZerosStart <= end) {
                // Print range with trailing zeros range in the middle
                assert(trailingZerosCount > 0);
                os.rdbuf()->sputn(&*begin, trailingZerosStart - begin);
                putCharacter(ctype.widen('0'), trailingZerosCount);
                os.rdbuf()->sputn(&*trailingZerosStart, end - trailingZerosStart);
            } else {
                // Print range as-is
                os.rdbuf()->sputn(&*begin, end - begin);
            }
        };
 
        // Pad the buffer if necessary.
        // Note that the length of trailing zeros is counted towards the length of the content.
        const auto contentSize = end - buffer.begin() + trailingZerosCount;
        if(contentSize >= width) {
            // Buffer needs no padding, output as-is
            putRange(buffer.begin(), end);
        } else {
            const auto padSize = width - contentSize;
            switch(adjustfield) {
            case std::ios_base::left:
                // Content is left-aligned, so output the buffer, followed by the padding
                putRange(buffer.begin(), end);
                putCharacter(os.fill(), padSize);
                break;
            case std::ios_base::internal:
                // Content is internally aligned, so output the buffer up to the "internal pad"
                // point, followed by the padding, followed by the remainder of the buffer.
                putRange(buffer.begin(), internalPad);
                putCharacter(os.fill(), padSize);
                putRange(internalPad, end);
                break;
            default:
                // Content is right-aligned, so output the padding, followed by the buffer
                putCharacter(os.fill(), padSize);
                putRange(buffer.begin(), end);
                break;
            }
        }
 
        // Width is reset after every write
        os.width(0);
 
        return os;
    }
 
    template<typename CharT, class Traits, typename B, unsigned int F>
    std::basic_istream<CharT, Traits>& operator>>(std::basic_istream<CharT, Traits>& is, FixedPoint<B, F>& x) {
        const typename std::basic_istream<CharT, Traits>::sentry sentry(is);
        if(!sentry) {
            return is;
        }
 
        const auto& ctype = std::use_facet<std::ctype<CharT>>(is.getloc());
        const auto& numpunct = std::use_facet<std::numpunct<CharT>>(is.getloc());
 
        bool thousandsSeparatorAllowed = false;
        const bool supportsThousandsSeparators = !numpunct.grouping().empty();
 
        const auto& isValidCharacter = [](const char ch) {
            // Note: allowing ['p', 'i', 'n', 't', 'y'] is technically in violation of the spec (we are emulating std::num_get),
            // but otherwise we cannot parse hexfloats and "infinity". This is a known issue with the spec (LWG #2381).
            return std::isxdigit(ch) || ch == 'x' || ch == 'X' || ch == 'p' || ch == 'P' || ch == 'i' || ch == 'I' || ch == 'n' || ch == 'N' || ch == 't' ||
                ch == 'T' || ch == 'y' || ch == 'Y' || ch == '-' || ch == '+';
        };
 
        const auto& peek = [&]() {
            for(;;) {
                auto ch = is.rdbuf()->sgetc();
                if(ch == Traits::eof()) {
                    is.setstate(std::ios::eofbit);
                    return '\0';
                }
                if(ch == numpunct.decimal_point()) {
                    return '.';
                }
                if(ch == numpunct.thousands_sep()) {
                    if(!supportsThousandsSeparators || !thousandsSeparatorAllowed) {
                        return '\0';
                    }
                    // Ignore valid thousands separators
                    is.rdbuf()->sbumpc();
                    continue;
                }
                auto res = ctype.narrow(ch, 0);
                if(!isValidCharacter(res)) {
                    // Invalid character: end input
                    return '\0';
                }
                return res;
            }
        };
 
        const auto& bump = [&]() {
            is.rdbuf()->sbumpc();
        };
 
        const auto& next = [&]() {
            bump();
            return peek();
        };
 
        bool negate = false;
        auto ch = peek();
        if(ch == '-') {
            negate = true;
            ch = next();
        } else if(ch == '+') {
            ch = next();
        }
 
        const char infinity[] = "infinity";
        // Must be "inf" or "infinity"
        int index = 0;
        while(index < 8 && ch == infinity[index]) {
            ++index;
            ch = next();
        }
 
        if(index > 0) {
            if(index == 3 || index == 8) {
                x = negate ? std::numeric_limits<FixedPoint<B, F>>::min() : std::numeric_limits<FixedPoint<B, F>>::max();
            } else {
                is.setstate(std::ios::failbit);
            }
            return is;
        }
 
        char exponentChar = 'e';
        int base = 10;
 
        constexpr auto noFraction = std::numeric_limits<std::size_t>::max();
        std::size_t fractionStart = noFraction;
        std::vector<unsigned char> significand;
 
        if(ch == '0') {
            ch = next();
            if(ch == 'x' || ch == 'X') {
                // Hexfloat
                exponentChar = 'p';
                base = 16;
                ch = next();
            } else {
                significand.push_back(0);
            }
        }
 
        // Parse the significand
        thousandsSeparatorAllowed = true;
        for(;; ch = next()) {
            if(ch == '.') {
                if(fractionStart != noFraction) {
                    // Double decimal point. Stop parsing.
                    break;
                }
                fractionStart = significand.size();
                thousandsSeparatorAllowed = false;
            } else {
                unsigned char val = base;
                if(ch >= '0' && ch <= '9') {
                    val = ch - '0';
                } else if(ch >= 'a' && ch <= 'f') {
                    val = ch - 'a' + 10;
                } else if(ch >= 'A' && ch <= 'F') {
                    val = ch - 'A' + 10;
                }
                if(val < 0 || val >= base) {
                    break;
                }
                significand.push_back(val);
            }
        }
        if(significand.empty()) {
            // We need a significand
            is.setstate(std::ios::failbit);
            return is;
        }
        thousandsSeparatorAllowed = false;
 
        if(fractionStart == noFraction) {
            // If we haven't seen a fraction yet, place it at the end of the significand
            fractionStart = significand.size();
        }
 
        // Parse the exponent
        bool exponentOverflow = false;
        std::size_t exponent = 0;
        bool exponentNegate = false;
        if(std::tolower(ch) == exponentChar) {
            ch = next();
            if(ch == '-') {
                exponentNegate = true;
                ch = next();
            } else if(ch == '+') {
                ch = next();
            }
 
            bool parsed = false;
            while(std::isdigit(ch)) {
                if(exponent <= std::numeric_limits<int>::max() / 10) {
                    exponent = exponent * 10 + (ch - '0');
                } else {
                    exponentOverflow = true;
                }
                parsed = true;
                ch = next();
            }
            if(!parsed) {
                // If the exponent character is given, the exponent value may not be empty
                is.setstate(std::ios::failbit);
                return is;
            }
        }
 
        // We've parsed all we need. Construct the value.
        if(exponentOverflow) {
            // Absolute exponent is too large
            if(std::all_of(significand.begin(), significand.end(), [](const unsigned char x) { return x == 0; })) {
                // Significand is zero. Exponent doesn't matter.
                x = FixedPoint<B, F>(B(0));
            } else if(exponentNegate) {
                // A huge negative exponent approaches 0.
                x = FixedPoint<B, F>(B(0));
            } else {
                // A huge positive exponent approaches infinity.
                x = std::numeric_limits<FixedPoint<B, F>>::max();
            }
            return is;
        }
 
        // Shift the fraction offset according to exponent
        {
            const auto exponentMultiplier = (base == 10) ? 1 : 4;
            if(exponentNegate) {
                const auto adjust = std::min(exponent / exponentMultiplier, fractionStart);
                fractionStart -= adjust;
                exponent -= adjust * exponentMultiplier;
            } else {
                const auto adjust = std::min(exponent / exponentMultiplier, significand.size() - fractionStart);
                fractionStart += adjust;
                exponent -= adjust * exponentMultiplier;
            }
        }
 
        constexpr auto isSigned = std::is_signed<B>::value;
        constexpr auto intBits = sizeof(B) * 8 - F - (isSigned ? 1 : 0);
        constexpr auto maxInt = (B{1} << intBits) - 1;
        constexpr auto maxFraction = (B{1} << F) - 1;
        constexpr auto maxValue = (B{1} << sizeof(B) * 8) - 1;
 
        // Parse the integer part
        B integer = 0;
        for(std::size_t i = 0; i < fractionStart; ++i) {
            if(integer > maxInt / base) {
                // Overflow
                x = negate ? std::numeric_limits<FixedPoint<B, F>>::min() : std::numeric_limits<FixedPoint<B, F>>::max();
                return is;
            }
            assert(significand[i] < base);
            integer = integer * base + significand[i];
        }
 
        // Parse the fractional part
        B fraction = 0;
        B divisor = 1;
        for(std::size_t i = fractionStart; i < significand.size(); ++i) {
            assert(significand[i] < base);
            if(divisor > maxFraction / base) {
                // We're done
                break;
            }
            fraction = fraction * base + significand[i];
            divisor *= base;
        }
 
        // Construct the value from the parsed parts
        B rawValue = (integer << F) + (fraction << F) / divisor;
 
        // Apply remaining exponent
        if(exponentChar == 'p') {
            // Base-2 exponent
            if(exponentNegate) {
                rawValue >>= exponent;
            } else {
                rawValue <<= exponent;
            }
        } else {
            // Base-10 exponent
            if(exponentNegate) {
                B remainder = 0;
                for(std::size_t e = 0; e < exponent; ++e) {
                    remainder = rawValue % 10;
                    rawValue /= 10;
                }
                rawValue += remainder / 5;
            } else {
                for(std::size_t e = 0; e < exponent; ++e) {
                    if(rawValue > maxValue / 10) {
                        // Overflow
                        x = negate ? std::numeric_limits<FixedPoint<B, F>>::min() : std::numeric_limits<FixedPoint<B, F>>::max();
                        return is;
                    }
                    rawValue *= 10;
                }
            }
        }
        x = FixedPoint<B, F>(static_cast<B>(negate ? -rawValue : rawValue));
 
        return is;
    }
 
} // namespace CeresEngine