How to convert a number of pseudorandom bits into a statistically random floating-point value between 0 and 1?

Question

I have an implementation of the xoshiro256** PRNG algorithm in an application I'm writing in C#. This works great for generating pseudorandom values between 0 and UInt64.MaxValue, but I've hit a spot where I need a pseudorandom double-precision floating-point value between 0 inclusive and 1 exclusive (so between 0 and 0.99999999...).

I know I can use BitConverter to do a "brutal conversion" from a ulong to a double, but I'm fairly certain this will give me values that are somewhere between the number of miles in a planck length and the number of cubic millimeters in the universe, and their negatives, along with the possibility of getting infinity, negative infinity, negative 0, and NaN. That's a bit ambitious (read: entirely unsuitable) for what I'm trying to do, so I'm hoping for some way of controlling the output I get so that it can be used for dealing with percent chances of things.

I don't know enough about how IEEE floating-point values work to know for sure what bits to put where in order to get the values I'm going for. I think (and am wrong) that I can just right-shift a ulong by 12 bits (thus turning the top 52 bits into the bottom 52 bits), add 2^52 (setting the bottom bit of the exponent to 1), and then BitConverter the resulting mess into a double. So something like this:

public static double DoubleFromRand(ulong rand) {
    ulong resultUL = rand >> 12;
    resultUL += ULongPow(2UL, 52UL); // adapted from https://stackoverflow.com/a/383596/19474638
    return BitConverter.ToDouble(BitConverter.GetBytes(resultUL), 0);
}

public static ulong ULongPow(ulong x, ulong pow) {
    ulong ret = 1UL;
    while (pow != 0UL) {
        if ((pow & 1UL) == 1UL) {
            ret *= x;
        }
        x *= x;
        pow >>= 1;
    }
    return ret;
}

If this worked, I would expect that passing UInt64.MaxValue to this would give me some value very close to 1 but not quite there. What I actually get from the above algorithm is some really weird small value.

Any clue what to do here? Using C# 4.0 on Mono 6.8.0.105, on Raspberry Pi OS 64-bit, on a Raspberry Pi 4 Model B. Note that I am not interested in using "real" .NET.

(Related: How to convert a uint64_t to a double/float between 0 and 1 with maximum accuracy (C++)? This doesn't answer my question as it converts a ulong to a double with mathematical operations, which I believe doesn't guarantee the statistical randomness of the result.

Also answers relating on how to use the built-in random features of C# are not helpful as I want to use my own implementation of a PRNG. I also am not helped by answers about how to generate random 64-bit floats. I want, very specifically, to take some quantity of randomly generated bits and force-convert them into a float that will land between 0 and 1. The question is about how to do the conversion, not the generation of random things.)

Try `double DoubleFromRand(ulong rand) { return (double)(rand >> 12)/(1ul << 52); }`. — chux - Reinstate Monica, Apr 11 '23 at 21:29
Do not try to fabricate the representation of a `double`. Simply get 52 random bits as a binary numeral, convert them to a `double` (as a plain integer-to-`double` conversion), and multiply by 2^−52 to scale to [0, 1). That is what chux code above does, and it is fine for most applications. — Eric Postpischil, Apr 11 '23 at 23:01
@EricPostpischil: Actually, you can get **53** bits, due to IEEE 754's hidden bit trick giving you an extra bit of precision. — dan04, Apr 11 '23 at 23:47
@EricPostpischil Hmm, it appears that my code and chux code are doing the same thing, but chux code is obviously more portable. That might work then. — ArrayBolt3, Apr 11 '23 at 23:53

chux - Reinstate Monica · Accepted Answer · 2023-04-12T12:18:49.790

A common approach to converts "bits into a statistically random floating-point value between 0 and 1?" begins with converting such bits into a value between [1.0 ... 2.0) and then subtracting 1.0.

Floating point numbers are distributed linearly between successive powers-of-2.

With common double encoding, there are 2⁵² values between most powers of 2. So take the 64-bits of random data and reduce to 52 and form a random number [1.0... 2.0) (OP's code does this with various bit manipulations), then subtract 1.0.

This is OP's self-answer approach.

When done well, it provides a fast approach, yet has additional portability issues over the following.

Another approach is to use 52 of the 64 supplied random bits and form a value [0...2⁵²-1] and divide by 2⁵². The shift, conversion from ulong to double and division are all expected to be exact and not dependent on rounding mode. Thus we prevent statistical biases from creeping in, yet only form 2⁵² different values.

// Something like
double DoubleFromRand(ulong rand) { 
  return (double)(rand >> 12)/(1ul << 52);
}

We can even use 1 bit more @dan04 and still have "shift, conversion from ulong to double and division are all expected to be exact and not dependent on rounding mode". Using more than 53 bits loses this property.

double DoubleFromRand53(ulong rand) {
  return (double)(rand >> 11)/(1ul << 53);
}

Another approach, not shown, would use 64 of the 64 bits and form a value [0...2⁶⁴-1] and divide by 2⁶⁴. This unfortunately incurs rounding in the integer to conversion and more rounding in the division imparting an undesirable bias and range (1.0 may get returned), yet does provide closer to 2⁶⁴ different values.

As dan04 noted, we can use 53 bits. There are 2^52 representable numbers in [½, 1), and we can use 2^52 numbers with the same spacing in [0, ½), making 2^53 for [0, 1). — Eric Postpischil, Apr 12 '23 at 11:30
@EricPostpischil Yes, I recall that now too. Answer amended. — chux - Reinstate Monica, Apr 12 '23 at 12:19
I benchmarked my original solution versus this solution, and discovered that the solution in this answer is multiple times faster than my BitConverter gymnastics. So it's a full-on win in every way - it's portable, safe, and fast. — ArrayBolt3, Jul 09 '23 at 07:10

ArrayBolt3 · Answer 2 · 2023-07-09T07:11:06.183

(Edit: As pointed out in the comments in my question, this isn't a great way of doing things, but for in the event someone has a good reason to do this sort of thing, I'll leave this here. Also, this is multiple times slower than the solution by chux, so use theirs unless you have a really good reason not to.)

I think I figured it out after some experimentation.

The exponent field uses a value of 1023 to represent 0 - larger values are positive, smaller ones are negative. So a value of 1023 in the exponent allows the rest of the significand field to be interpreted "as-is", i.e. - everything after the decimal point in the number.

Doing this makes UInt64.MaxValue map to ~1.999, which is almost 2 (and it will print as 2 if you use Console.WriteLine() on it without setting the number of digits to display). UInt64.MaxValue / 2 maps to 1.5, UInt64.MaxValue / 4 maps to 1.25.

Which is great except for the fact that there's an extra 1.0 on all of these values. Thankfully it looks like I can get rid of that by simply subtracting 1.0 from the final result. This appears to be working at first glance (UInt64.MaxValue / 4 now maps to 0.25). I'm not sure if this will cause errors with smaller values, but hopefully it won't.

The final I-think-it's-working code is:

public static double DoubleFromRand(ulong rand) {
    ulong resultUL = rand >> 12;
    resultUL += ((ulong)1023 << 52);
    return BitConverter.ToDouble(BitConverter.GetBytes(resultUL), 0) - 1.0;
}

Note that this mess is probably not portable across all processors - according to Wikipedia (where I got all my info about IEEE754 floating-point math), some CPUs do really weird things with endianness when working with floating-point numbers, so depending on your CPU, this may give you seriously scrambled results. But so far, on my particular 64-bit ARM CPU, this seems to work.

How to convert a number of pseudorandom bits into a statistically random floating-point value between 0 and 1?

2 Answers2