How do I generate random subnormal numbers?

An answer to this question on Stack Overflow.

Question

I'd like to verify a piece of code works on subnormal numbers, so I'd like to generate a bunch of random subnormal single-precision numbers (including zero). How can I do this?

Solutions such as dividing by a large number to get a subnormal may just round to zero or, in the best case, probably won't give an even distribution.

nextafter iteration could work, but it'd be slow!

Answer

First, recall that a

With this in hand we have the following strategy:

  1. Draw 32/64 bits uniformly. We won't need all of these bits, but it'll be more convenient and faster to draw them this way. We'll just randomize as many bits as a floating-point number requires and then...
  2. Mask out the exponent bits so they are zero
  3. Convert the bit pattern into a floating-point number

A caveat is that the endianness of the exponent bit mask must match the endianness of the floating-point values. This is the case for most hardware, but you should test it if you want to be exceptionally rigorous or are working on something exotic.

All that said, we get this code:

// Compile with: clang++.par -O3 -march=native test2.cpp --std=c++20 -Wall -Wextra -pedantic -Werror
#include <concepts>
#include <iostream>
#include <random>
#include <type_traits>
template<std::floating_point T>
class uniform_subnormal_distribution {
  private:
    // float is SignBit + 8 Exponent Bits + 23 Mantissa Bits
    static constexpr uint32_t subnormal_mask32 = 0x807FFFFF;
    // double is SignBit + 11 Exponent Bits + 52 Mantissa Bits
    static constexpr uint64_t subnormal_mask64 = 0x800FFFFFFFFFFFFF;
  public:
    template<class Engine>
    T operator()(Engine& eng) const {
        if constexpr (std::is_same_v<T, float>){
            std::uniform_int_distribution<uint32_t> dist;
            // Get uniformaly distributed bits
            const uint32_t bits = dist(eng);
            // Make the exponent all zeros
            const uint32_t subnormal_bits = bits & subnormal_mask32;
            // Retrieve a floating-point value from the bits
            return std::bit_cast<float, uint32_t>(subnormal_bits);
        } else if constexpr (std::is_same_v<T, double>){
            std::uniform_int_distribution<uint64_t> dist;
            const uint64_t bits = dist(eng);
            const uint64_t subnormal_bits = bits & subnormal_mask32;
            return std::bit_cast<double, uint64_t>(subnormal_bits);
        } else {
            // can't use 'false' -- expression has to depend on a template parameter
            static_assert(!sizeof(T*), "Unsupported floating-point type");
        }
    }
};
int main(){
    std::random_device rd;
    std::mt19937 mt(rd());
    uniform_subnormal_distribution<float> dist;
    std::vector<float> res;
    for (unsigned i = 0; i < 20; i++) {
        const auto float_val = dist(mt);
        std::cout<<float_val<<std::endl;
        res.push_back(float_val);
    }
    return 0;
}