Fast even random float within range

Question

I did the following lcg (linear congruency generator) to draw unsigned int from 0x00000000 to 0xFFFFFFFF uniformly. Just for testing, I used the seed equal to zero and raffling the first five numbers to see what comes out.

# include <stdio.h>

unsigned long long int state = 0 ;

unsigned int rnd() {
    state *= 0x60BDC8431972EFA5ull ;
    state += 0x4D268FBC9E53A107ull ;
    return (unsigned int)( state >> 32 ) ;
}

int main(void) {
    printf( "%u\n" , rnd() ) ;
    printf( "%u\n" , rnd() ) ;
    printf( "%u\n" , rnd() ) ;
    printf( "%u\n" , rnd() ) ;
    printf( "%u\n" , rnd() ) ;
    return 0 ;
}

The result is this.

1294372796
3162674411
1090754839
3480286764
3255825657

To draw a number in single floating point format within the range [a,b] i did the following function which simply divides the drawn integer by 0xFFFFFFFFu to act as a draw of 0.0f to 1.0f and uses the result to make lerp of a to b.

float random( float a , float b ){
    float f0t1 = rnd()/(float)(0xFFFFFFFFu) ;
    return a * ( 1-f0t1 ) + f0t1 * b ;
}

Testing like this,

int main(void) {
    printf( "%f\n" , random(0,1) ) ;
    printf( "%f\n" , random(0,1) ) ;
    printf( "%f\n" , random(0,1) ) ;
    printf( "%f\n" , random(0,1) ) ;
    printf( "%f\n" , random(0,1) ) ;
    return 0 ;
}

the result is this.

0.301370
0.736368
0.253961
0.810317
0.758056

I mean, it makes a division that costs a lot. You can even convert it into multiplication, but isn’t there a faster way to calculate it? I imagine that with magic numbers you could handle binary float code efficiently and calculate, correct?

Answer 1

Surely one can exchange the Typecast of the drawn integer and the division for two integer operations and the binary code conversion. I considered the compilation optimizations and took the opportunity to make step by step to be clear. Note that printing with the same main is of the same result.

float random( float a , float b ){
    unsigned int magic = rnd() ;          // Sorteia de 0x00000000 a 0xFFFFFFFF
    magic >>= 9 ;                         // Ficou entre 0x00000000 e 0x007FFFFF
    magic += 0x3F800000 ;                 // Ficou entre 0x3F800000 e 0x3FFFFFFF
    float f1t2 = *(float*)&magic ;        // Magicamente ficou entre 1.00000000 e 1.99999988
    float f0t1 = f1t2 - 0.99999994f ;     // Ficou entre 0.00000006 e 0.99999994
    return a * ( 1-f0t1 ) + f0t1 * b ;    // Finalmente ficou entre "a" e "b"
}

The nine-bit offset to the right is to make room for the signal bits and float exponent by eliminating the least significant digits. One can use magic &= ~0x3F800000 ; also and I don’t know if it impacts performance, but it eliminates the most significant digits and the result comes out different.

As well as the addition of the "magic number" 0x3F800000 (code of 1.0f) is to fill these bits. It may be an operation | instead of a sum, I don’t know if it impacts on performance but not on outcome.

Note that after conversion is subtracted not by 0x3F800000 = 1.0f (that would make the result go from 0.00000000 and 0.99999988) but rather for 0x3F800000 - 1 = 0.99999994f to balance border approaches.

I suppose the return of a * ( 2-f1t2 ) + ( f1t2-1 ) * b overriding the need to f0t1, but I imagine that the impact on performance is the same and rounding takes control of borders.

Now notice that embed "inline" the function rnd may result not only in a decrease in the number of calls but also possibly a reduction of an integer operation, but it is not known because a shift is made before converting from 64 to 32 bits and the other later. So the following can also be done by guarantee.

float random( float a , float b ){
    state *= 0x60BDC8431972EFA5ull ;
    state += 0x4D268FBC9E53A107ull ;
    unsigned int magic = (unsigned int)( state >> (32+9) ) ;
    magic += 0x3F800000 ;
    float f1t2 = *(float*)&magic ;
    float f0t1 = f1t2 - 0.99999994f ;
    return a * ( 1-f0t1 ) + f0t1 * b ;
}

And finally realize that the reasoning holds for double floating point, after all just consider that are 32 bits drawn unless you change the lcg (can even use the 64 bits generated, as shown below), 12 bits of signal and exponent, the "magic number" equal to 1.0 is 0x3FF0000000000000 and the subtraendo 0x3FF0000000000000 - 1 magically is 0.9999999999999999.

double random( double a , double b ){
    state *= 0x60BDC8431972EFA5ull ;
    state += 0x4D268FBC9E53A107ull ;
    unsigned long long int magic = state >> 12 ;
    magic += 0x3FF0000000000000ull ;
    double f1t2 = *(double*)&magic ;
    double f0t1 = f1t2 - 0.9999999999999999 ;
    return a * ( 1-f0t1 ) + f0t1 * b ;
}

And notice that once again the printed result is the same.