libdpf/include/grotto/hexfloat.hpp

/// @file grotto/hexfloat.hpp
/// @author Ryan Henry <ryan.henry@ucalgary.ca>
/// @brief Hexadecimal floating-point formatting for IEEE doubles and wide fixed-point words.
/// @details Prints IEEE doubles and wide fixed-point words in hexadecimal.
/// @copyright Copyright (c) 2019-2023 Ryan Henry and others
/// @license Released under a GNU General Public v2.0 (GPLv2) license;
///          see [LICENSE.md](@ref license) for details.

#ifndef LIBDPF_INCLUDE_GROTTO_HEXFLOAT_HPP__
#define LIBDPF_INCLUDE_GROTTO_HEXFLOAT_HPP__

#include "hedley/hedley.h"

#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <stdexcept>
#include <string>
#include <string_view>
#include <type_traits>
#include <vector>

#include "grotto/fixedpoint.hpp"

namespace grotto
{
/// \complexity Extracts the IEEE fields and prints a fixed-size mantissa. `Θ(1)`.
/// @see grotto::from_hexfloat
/// @param x the double
/// @return the hexfloat spelling

template <class CharT = char,
          class Traits = std::char_traits<CharT>,
          class Allocator = std::allocator<CharT>>
std::basic_string<CharT, Traits, Allocator> to_hexfloat(double x)
{
    static constexpr auto base_indicator = "0x";
    if (std::isnan(x))
    {
        return "NaN";
    }

    if (std::isinf(x))
    {
        return std::string(std::signbit(x) ? "-" : "+") + "Infinity";
    }

    if (iszero(x))
    {
        return std::string(std::signbit(x) ? "-" : "+")
            + "0x0.0000000000000p+0000";
    }

    uint32_t X[2];
    std::memcpy(&X[0], &x, sizeof(X));
    const int32_t sign_mask = 0x7fffffff;
    const int32_t exponent_mask  = 0x7ff00000;
#if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_BIG_ENDIAN__)
    constexpr unsigned hi = 0;
#else
    constexpr unsigned hi = 1;
#endif
    constexpr unsigned lo = 1u - hi;

    // extract and clear the exponent
    int32_t exponent = ((X[hi] & exponent_mask) >> 20) - 0x3ff;
    if (exponent == -1023)
    {
        x *= std::exp2(1023);
        std::memcpy(&X[0], &x, sizeof(X));
        exponent += ((X[hi] & exponent_mask) >> 20) - 0x3ff;
    }
    X[hi] = (X[hi] & ~exponent_mask) ^ 0x3ff00000;

    // extract and clear the sign bit
    std::string sign_prefix = !!(X[hi] & ~sign_mask) ? "-" : "+";
    X[hi] &= sign_mask;

    uint64_t mantissa = X[lo] | (uint64_t(X[hi]) & 0xfffff) << 32;
    char buf[13+1];
    snprintf(buf, 13+1, "%013lx", mantissa);
    std::string hex_mantissa = std::string("1.") + std::string(buf);

    std::string expsign = exponent < 0 ? "p-" : "p+";
    exponent = std::abs(exponent);
    std::string expnum = std::to_string(exponent);
    std::string exponent_suffix = expsign + std::string(4 - std::min(size_t(4), expnum.length()), '0') + expnum;

    return sign_prefix + base_indicator + hex_mantissa + exponent_suffix;
}

namespace detail
{

template <typename T>
constexpr bool is_hexfloat_integer =
    std::is_integral_v<T>
    || std::is_same_v<T, simde_int128>
    || std::is_same_v<T, simde_uint128>
    || std::is_same_v<T, uint128_t>
    || std::is_same_v<T, uint256_t>;

// libstdc++ reports std::is_signed_v<__int128> as false.
template <typename T>
constexpr bool hexfloat_is_signed =
    std::is_signed_v<T> || std::is_same_v<T, simde_int128>;

template <typename T>
constexpr bool is_hexfloat_fixedpoint = false;

template <unsigned FractionalBits, typename Integral>
constexpr bool is_hexfloat_fixedpoint<fixedpoint<FractionalBits, Integral>> = true;

HEDLEY_NO_THROW
constexpr void mask_low_bits(std::uint64_t * limbs, unsigned width) noexcept
{
    if (width >= 256u)
    {
        return;
    }
    const unsigned limb = width / 64u;
    const unsigned rem = width % 64u;
    if (rem == 0u)
    {
        for (unsigned i = limb; i < 4u; ++i)
        {
            limbs[i] = 0;
        }
        return;
    }
    limbs[limb] &= (std::uint64_t{1} << rem) - 1u;
    for (unsigned i = limb + 1u; i < 4u; ++i)
    {
        limbs[i] = 0;
    }
}

HEDLEY_NO_THROW
constexpr bool bit_is_set(const std::uint64_t * limbs, unsigned bit) noexcept
{
    return ((limbs[bit / 64u] >> (bit % 64u)) & 1u) != 0u;
}

HEDLEY_NO_THROW
constexpr void negate_low_bits(std::uint64_t * limbs, unsigned width) noexcept
{
    for (unsigned i = 0; i < 4u; ++i)
    {
        limbs[i] = ~limbs[i];
    }
    std::uint64_t carry = 1;
    for (unsigned i = 0; i < 4u && carry != 0u; ++i)
    {
        const std::uint64_t old = limbs[i];
        limbs[i] = old + carry;
        carry = limbs[i] < old ? 1u : 0u;
    }
    mask_low_bits(limbs, width);
}

template <typename T>
HEDLEY_NO_THROW
constexpr void store_integer_bits(const T & value, std::uint64_t out[4]) noexcept
{
    out[0] = out[1] = out[2] = out[3] = 0;
    if constexpr (std::is_same_v<T, uint256_t>)
    {
        out[0] = value.lower().lower();
        out[1] = value.lower().upper();
        out[2] = value.upper().lower();
        out[3] = value.upper().upper();
    }
    else if constexpr (std::is_same_v<T, uint128_t>)
    {
        out[0] = value.lower();
        out[1] = value.upper();
    }
    else if constexpr (std::is_same_v<T, simde_int128> || std::is_same_v<T, simde_uint128>)
    {
        const simde_uint128 bits = static_cast<simde_uint128>(value);
        out[0] = static_cast<std::uint64_t>(bits);
        out[1] = static_cast<std::uint64_t>(bits >> 64);
    }
    else
    {
        using unsigned_same = std::make_unsigned_t<T>;
        out[0] = static_cast<std::uint64_t>(static_cast<unsigned_same>(value));
    }
}

template <typename T>
HEDLEY_NO_THROW
constexpr T load_integer_bits(const std::uint64_t * limbs) noexcept
{
    if constexpr (std::is_same_v<T, uint256_t>)
    {
        return uint256_t{uint128_t{limbs[3], limbs[2]}, uint128_t{limbs[1], limbs[0]}};
    }
    else if constexpr (std::is_same_v<T, uint128_t>)
    {
        return uint128_t{limbs[1], limbs[0]};
    }
    else if constexpr (std::is_same_v<T, simde_int128> || std::is_same_v<T, simde_uint128>)
    {
        const simde_uint128 bits = simde_uint128(limbs[0]) | (simde_uint128(limbs[1]) << 64);
        return static_cast<T>(bits);
    }
    else
    {
        using unsigned_same = std::make_unsigned_t<T>;
        return static_cast<T>(static_cast<unsigned_same>(limbs[0]));
    }
}

HEDLEY_NO_THROW
inline int highest_bit(const std::uint64_t * limbs) noexcept
{
    for (int i = 3; i >= 0; --i)
    {
        if (limbs[i] != 0u)
        {
            return i * 64 + 63 - __builtin_clzll(limbs[i]);
        }
    }
    return -1;
}

/// @brief Exact hexfloat of a two's-complement word whose binary point sits `fractional_bits`
/// below bit 0. The fraction keeps every bit below the leading 1, so a 256-bit word
/// round-trips.
/// @tparam T value type
/// @param value the value to convert or store
/// @param fractional_bits the number of fractional bits
/// @return Exact hexfloat of a two's-complement word whose binary point sits `fractional_bits`
/// below bit 0
template <typename T>
std::string format_hexfloat(const T & value, int fractional_bits)
{
    const unsigned width = static_cast<unsigned>(dpf::utils::bitlength_of_v<T>);
    std::uint64_t mag[4];
    store_integer_bits(value, mag);
    mask_low_bits(mag, width);
    const bool negative = hexfloat_is_signed<T> && width > 0u && bit_is_set(mag, width - 1u);
    if (negative)
    {
        negate_low_bits(mag, width);
    }
    const int msb = highest_bit(mag);
    if (msb < 0)
    {
        return std::string(negative ? "-" : "+") + "0x0.0p+0000";
    }

    std::string fraction;
    for (int bit = msb - 1; bit >= 0; )
    {
        unsigned nibble = 0;
        for (int k = 0; k < 4; ++k, --bit)
        {
            nibble <<= 1u;
            if (bit >= 0 && bit_is_set(mag, static_cast<unsigned>(bit)))
            {
                nibble |= 1u;
            }
        }
        fraction.push_back("0123456789abcdef"[nibble]);
    }
    while (fraction.size() > 1u && fraction.back() == '0')
    {
        fraction.pop_back();
    }
    if (fraction.empty())
    {
        fraction = "0";
    }

    const long exponent = static_cast<long>(msb) - static_cast<long>(fractional_bits);
    char expbuf[32];
    std::snprintf(expbuf, sizeof expbuf, "%c%04ld",
        exponent < 0 ? '-' : '+', std::labs(exponent));
    return std::string(negative ? "-" : "+") + "0x1." + fraction + "p" + expbuf;
}

HEDLEY_NO_THROW
inline int hex_value(char c) noexcept
{
    if (c >= '0' && c <= '9') return c - '0';
    if (c >= 'a' && c <= 'f') return c - 'a' + 10;
    if (c >= 'A' && c <= 'F') return c - 'A' + 10;
    return -1;
}

inline void shift_left_limbs(std::vector<std::uint64_t> & limbs, unsigned shift)
{
    if (shift == 0u)
    {
        return;
    }
    const std::size_t limb_shift = shift / 64u;
    const unsigned bit_shift = shift % 64u;
    std::vector<std::uint64_t> out(limbs.size() + limb_shift + 1u, 0);
    for (std::size_t i = 0; i < limbs.size(); ++i)
    {
        out[i + limb_shift] |= limbs[i] << bit_shift;
        if (bit_shift != 0u)
        {
            out[i + limb_shift + 1u] |= limbs[i] >> (64u - bit_shift);
        }
    }
    limbs.swap(out);
}

inline void shift_right_limbs(std::vector<std::uint64_t> & limbs, unsigned shift)
{
    if (shift == 0u)
    {
        return;
    }
    const std::size_t limb_shift = shift / 64u;
    const unsigned bit_shift = shift % 64u;
    if (limb_shift >= limbs.size())
    {
        limbs.assign(1, 0);
        return;
    }
    std::vector<std::uint64_t> out(limbs.size() - limb_shift, 0);
    for (std::size_t i = 0; i < out.size(); ++i)
    {
        out[i] = limbs[i + limb_shift] >> bit_shift;
        if (bit_shift != 0u && i + limb_shift + 1u < limbs.size())
        {
            out[i] |= limbs[i + limb_shift + 1u] << (64u - bit_shift);
        }
    }
    limbs.swap(out);
}

/// @brief Parse `±0x<hex>.<hex>p±<dec>` into the low `width` bits of
/// significand * 2^(exponent + fractional_bits), floored toward -inf on the
/// discarded fraction and wrapped modulo 2^width.
/// @tparam T value type
/// @param text hexfloat spelling
/// @param fractional_bits the number of fractional bits
/// @return the returned `T`
/// @throws std::invalid_argument if `empty hexfloat`
template <typename T>
T integer_from_hexfloat(std::string_view text, int fractional_bits)
{
    const unsigned width = static_cast<unsigned>(dpf::utils::bitlength_of_v<T>);
    std::size_t i = 0;
    while (i < text.size() && std::isspace(static_cast<unsigned char>(text[i])))
    {
        ++i;
    }
    std::size_t end = text.size();
    while (end > i && std::isspace(static_cast<unsigned char>(text[end - 1])))
    {
        --end;
    }
    if (i >= end)
    {
        throw std::invalid_argument("empty hexfloat");
    }

    bool negative = false;
    if (text[i] == '+' || text[i] == '-')
    {
        negative = text[i] == '-';
        ++i;
    }
    if (i + 1 < end && text[i] == '0' && (text[i + 1] == 'x' || text[i + 1] == 'X'))
    {
        i += 2;
    }

    std::vector<int> nibbles;
    bool saw_digit = false;
    bool after_point = false;
    int fraction_digits = 0;
    for (; i < end; ++i)
    {
        const char c = text[i];
        if (c == '.')
        {
            if (after_point)
            {
                throw std::invalid_argument("malformed hexfloat");
            }
            after_point = true;
            continue;
        }
        const int digit = hex_value(c);
        if (digit < 0)
        {
            break;
        }
        nibbles.push_back(digit);
        saw_digit = true;
        if (after_point)
        {
            ++fraction_digits;
        }
    }
    if (!saw_digit || i >= end || (text[i] != 'p' && text[i] != 'P'))
    {
        throw std::invalid_argument("malformed hexfloat");
    }
    ++i;
    int exp_sign = 1;
    if (i < end && (text[i] == '+' || text[i] == '-'))
    {
        exp_sign = text[i] == '-' ? -1 : 1;
        ++i;
    }
    if (i >= end || !std::isdigit(static_cast<unsigned char>(text[i])))
    {
        throw std::invalid_argument("malformed hexfloat exponent");
    }
    long exponent = 0;
    for (; i < end && std::isdigit(static_cast<unsigned char>(text[i])); ++i)
    {
        exponent = std::min(1000000000L, exponent * 10 + (text[i] - '0'));
    }
    if (i != end)
    {
        throw std::invalid_argument("malformed hexfloat");
    }
    exponent *= exp_sign;

    std::vector<std::uint64_t> sig(1, 0);
    for (int nibble : nibbles)
    {
        std::uint64_t carry = static_cast<unsigned>(nibble);
        for (std::size_t k = 0; k < sig.size(); ++k)
        {
            const unsigned __int128 prod = static_cast<unsigned __int128>(sig[k]) * 16u + carry;
            sig[k] = static_cast<std::uint64_t>(prod);
            carry = static_cast<std::uint64_t>(prod >> 64);
        }
        if (carry != 0u)
        {
            sig.push_back(carry);
        }
    }

    const long long placed = static_cast<long long>(exponent)
        - 4LL * fraction_digits
        + static_cast<long long>(fractional_bits);
    std::uint64_t bits[4] = {};
    const bool sig_zero = sig.size() == 1u && sig[0] == 0u;
    if (!sig_zero && placed < static_cast<long long>(width))
    {
        if (placed >= 0)
        {
            shift_left_limbs(sig, static_cast<unsigned>(placed));
        }
        else
        {
            const unsigned down = static_cast<unsigned>(std::min<long long>(-placed, 1000000000LL));
            shift_right_limbs(sig, down);
        }
        for (std::size_t k = 0; k < 4u && k < sig.size(); ++k)
        {
            bits[k] = sig[k];
        }
        mask_low_bits(bits, width);
        if (negative && highest_bit(bits) >= 0)
        {
            negate_low_bits(bits, width);
        }
    }
    return load_integer_bits<T>(bits);
}

}  // namespace detail
/// \complexity `format_hexfloat` walks the bits below the highest set bit, four per nibble. `Θ(width)` with `width = bitlength_of_v<T>` (at most 256). Extra space is the nibble string, `Θ(width)`.
/// @see grotto::fixedpoint

template <typename T>
std::enable_if_t<detail::is_hexfloat_integer<T>, std::string>
to_hexfloat(T value)
{
    return detail::format_hexfloat(value, 0);
}
/// \complexity Forwards the raw word to `format_hexfloat` with `FractionalBits` as the binary point. `Θ(backend width)`.
/// @see grotto::fixedpoint
/// @note Prints the mathematical value of the fixed-point word, using `integral_representation()`.

template <unsigned FractionalBits, typename Integral>
std::string to_hexfloat(fixedpoint<FractionalBits, Integral> value)
{
    return detail::format_hexfloat(value.integral_representation(), FractionalBits);
}
/// \complexity Scans the text once. Integer and fixed-point paths shift a limb vector by the exponent. Time `Θ(digits + |exponent shift|)`. Extra space is that limb vector.
/// A `double` target uses `strtod`.
/// @see grotto::to_hexfloat
/// @note Fixed-point parsing is bit-exact into `from_raw` / `make_fixed_from_integral_type`: discarded fraction bits are floored, and the word wraps modulo `2^width`.
/// @param text hexfloat spelling
/// @return the parsed value

template <typename T>
T from_hexfloat(std::string_view text)
{
    if constexpr (std::is_same_v<T, double>)
    {
        const std::string owned{text};
        char * end = nullptr;
        const double parsed = std::strtod(owned.c_str(), &end);
        if (end == owned.c_str())
        {
            throw std::invalid_argument("malformed hexfloat");
        }
        while (*end != '\0' && std::isspace(static_cast<unsigned char>(*end)))
        {
            ++end;
        }
        if (*end != '\0')
        {
            throw std::invalid_argument("malformed hexfloat");
        }
        return parsed;
    }
    else if constexpr (detail::is_hexfloat_fixedpoint<T>)
    {
        const auto bits = detail::integer_from_hexfloat<typename T::integral_type>(
            text, T::fractional_bits);
        return make_fixed_from_integral_type<T::fractional_bits, typename T::integral_type>(bits);
    }
    else
    {
        static_assert(detail::is_hexfloat_integer<T>,
            "from_hexfloat expects a floating, integral, or fixedpoint type");
        return detail::integer_from_hexfloat<T>(text, 0);
    }
}

}  // namespace grotto

#endif  // LIBDPF_INCLUDE_GROTTO_HEXFLOAT_HPP__