Using memcpy to easily and portably get the underlying representation of an IEEE-754 float

floatrep.cpp

#include <cstddef>
#include <cstdint>
#include <cstring>

#include <limits>

std::uint32_t floatrep(float x) {
    std::uint32_t n;
    std::memcpy(&n, &x, sizeof(x));
    return n;
}
std::uint64_t floatrep(double x) {
    std::uint64_t n;
    std::memcpy(&n, &x, sizeof(x));
    return n;
}
float makefloat(std::uint32_t n) {
    float x;
    std::memcpy(&x, &n, sizeof(n));
    return x;
}
double makefloat(std::uint64_t n) {
    double x;
    std::memcpy(&x, &n, sizeof(n));
    return x;
}

// NOTE: this template will recursion-overflow if WIDTH is larger than all the types for which it is specialised
template <std::size_t WIDTH>
struct min_type_for_width {
    using Type = min_type_for_width<WIDTH + 1>::Type;
};
template <>
struct min_type_for_width<8> {
    using Type = std::uint8_t;
};
template <>
struct min_type_for_width<16> {
    using Type = std::uint16_t;
};
template <>
struct min_type_for_width<32> {
    using Type = std::uint32_t;
};
template <>
struct min_type_for_width<64> {
    using Type = std::uint64_t;
};

// NOTE: this is technically "sign separated from significand/exponent", not "sign+magnitude"
template <std::size_t WIDTH>
struct float_sign_magnitude {
    bool sign : 1;
    min_type_for_width<WIDTH - 1>::Type magnitude : WIDTH - 1;
};

float_sign_magnitude<32> unpack(float x) {
    auto rep = floatrep(x);
    return {
        (bool)(rep >> 31),
        (rep << 1) >> 1,
    };
}
float_sign_magnitude<64> unpack(double x) {
    auto rep = floatrep(x);
    return {
        (bool)(rep >> 63),
        (rep << 1) >> 1,
    };
}
float pack(float_sign_magnitude<32> rep) {
    return makefloat(((std::uint32_t)rep.sign << 31) | rep.magnitude);
}
double pack(float_sign_magnitude<64> rep) {
    return makefloat(((std::uint64_t)rep.sign << 63) | rep.magnitude);
}

std::uint32_t ordinal(float x) {
    const std::uint32_t HALF = (std::numeric_limits<std::uint32_t>::max() / 2) + 1;
    auto rep = unpack(x);
    if (rep.sign) { // negative
        return HALF - rep.magnitude - 1;
    } else { // positive
        return HALF + rep.magnitude;
    }
}
std::uint64_t ordinal(double x) {
    const std::uint64_t HALF = (std::numeric_limits<std::uint64_t>::max() / 2) + 1;
    auto rep = unpack(x);
    if (rep.sign) { // negative
        return HALF - rep.magnitude - 1;
    } else { // positive
        return HALF + rep.magnitude;
    }
}

// nth (ordinal->float) is just the reverse of the above two
float nth(std::uint32_t n) {
    const std::uint32_t HALF = (std::numeric_limits<std::uint32_t>::max() / 2) + 1;
    float_sign_magnitude<32> rep;
    if (n < HALF) { // negative
        rep.sign = true;
        rep.magnitude = HALF - 1 - n;
    } else { // positive
        rep.sign = false;
        rep.magnitude = n - HALF;
    }
    return pack(rep);
}
double nth(std::uint64_t n) {
    const std::uint64_t HALF = (std::numeric_limits<std::uint64_t>::max() / 2) + 1;
    float_sign_magnitude<64> rep;
    if (n < HALF) { // negative
        rep.sign = true;
        rep.magnitude = HALF - 1 - n;
    } else { // positive
        rep.sign = false;
        rep.magnitude = n - HALF;
    }
    return pack(rep);
}
// If you didn't care about portability, you certainly could extend this to long double, probably using numeric_limits
// to get the number of bits used for exponent and significand. But because the size, precision, etc of long double is
// totally platform-dependent, we don't know explicitly what unsigned type it will fit in (actually, 80-bit extended
// won't even fit in 64-bit!). It's a no from me.

// NOTE: this code sorts negative NaNs *LOWER* than negative infinity (the NaNs are sorted by their "payload" aka
// the significand). Likewise, positive NaNs sort *HIGHER* than positive infinity.
// This may seem counterintuitive at first (it's reasonable to expect that nothing is smaller than -Inf or greater than
// +Inf), but when you consider where would the NaNs sort if they need to be ordered with everything else, beyond the
// infinities is the only logical place to put them, as then they're a kind of range of "out of bounds" values.
// One can always explicitly constrain to ordinals in the range (-Inf..+Inf) if one never wants NaNs.
//
// Maybe I should write additional versions of the conversion functions which exclude NaNs from the range of ordinals
// (such that -Inf really is the 0th float and +Inf the MAXth one).

#include <iostream>

int main() {
    std::cout << "-Inf is the " << ordinal(-std::numeric_limits<float>::infinity()) << "th float value" << std::endl;
    std::cout << "-Inf is the " << ordinal(-std::numeric_limits<double>::infinity()) << "th double value" << std::endl;
    std::cout << "+Inf is the " << ordinal(+std::numeric_limits<float>::infinity()) << "th float value" << std::endl;
    std::cout << "+Inf is the " << ordinal(+std::numeric_limits<double>::infinity()) << "th double value" << std::endl;
}

sort_floats.cpp

#include <iostream>

int main() {
    // NOTE: using this range excludes the NaNs, meaning we don't cover ALL values (just all the ones that are numeric)
    auto start = ordinal(0.0f); // this is the smallest known value that Google Sheets shows without scientific notation
    auto stop = ordinal(std::numeric_limits<float>::infinity());
    const std::size_t STEPS = 999; // Google Sheets start with 1000 rows, we leave one for the column headings :)
    for (decltype(start) i = 0; i <= (STEPS - 1); i++) {
        decltype(i) n = start + ((stop - start) * ((double)i / (STEPS - 1)));
        std::cout << n << " " << nth(n) << std::endl;
    }
}