IEEE754 doubles and uniform random distribution in C

IEEE754 floating point numbers have 3 parts: a sign (1bit), an exponent (11bits) and a mantissa (52bits).

To translate into C:

    union {
        struct {
            unsigned long long m : 52;
            unsigned long long e : 11;
            unsigned long long s : 1;
        };

        double d;
    } X;

To make a uniformly distributed number between -1 and 1, we need:

• A random mantissa m
• An exponent e starts at 1022 (which is actually -1) and has 50% chance of being 1022, 25% chance to be 1021, 12.5% chance to be 1020...
• a random sign s

This leads to the following:

#include <stdint.h>

// Please don't use C rand()
extern uint32_t decent_rand();

double uniform_rand_double_between_minus_one_and_one() {
    union {
        struct {
            unsigned long long m : 52;
            unsigned long long e : 11;
            unsigned long long s : 1;
        };

        double d;
    } X;

    uint32_t bit_divider = decent_rand();
    // Random mantissa
    X.m = decent_rand() + ((uint64_t)bit_divider<<32);
    // Base exponent
    X.e = 1022;
    
    // Recover unused entropy
    bit_divider >>= 20;
    // Infinite loop, for speed
    int i = 12;
    for(;;) {
        for(; i; --i) {
            // Consume entropy one bit at a time
            if(bit_divider & 1) {
                // 50% chance to reduce
                X.e -= 1;
                bit_divider >>= 1;
                // Handle underflow of the exponent
                if(X.e == 0) goto very_unlikely;
            } else {
                // 50% chance to be done
                bit_divider >>= 1;
                // Set the sign
                X.s = bit_divider & 1;
                return X.d;
            }
        }
        // We consumed all the entropy, we start again
        bit_divider = decent_rand();
        i = 32;
    }
    very_unlikely:
    X.s = bit_divider & 1;
    return X.d;
}

Then we can shift that space into another slice of the space, but we the sign is a hurdle here. Let's have the function without the sign:

double uniform_rand_double_between_zero_and_one_opti() {
    union {
        struct {
            unsigned long long m : 52;
            unsigned long long e : 11;
            unsigned long long s : 1;
        };

        double d;
    } X;

    uint32_t bit_divider = decent_rand();
    // Random mantissa
    X.m = decent_rand() + ((uint64_t)bit_divider<<32);
    // Base exponent
    X.e = 1022;
    
    // Recover unused entropy
    bit_divider >>= 20;
    // Infinite loop, for speed
    int i = 12;
    for(;;) {
        for(; i; --i) {
            // Consume entropy one bit at a time
            // Handle underflow of the exponent
            if(bit_divider & 1 && X.e != 0) [[unlikely]] {
                // 50% chance to reduce
                X.e -= 1;
                bit_divider >>= 1;
            } else {
                // 50% chance to be done
                goto out;
            }
        }
        // We consumed all the entropy, we start again
        bit_divider = decent_rand();
        i=32;
    }
    out:
    X.s = 0;
    return X.d;
}

    double uniform_rand_between(double a, double b) {
        return uniform_rand_double_between_zero_and_one()*(b-a) + a
    }

There you have it, a uniform random number generator for real numbers stored as IEEE754 doubles.