C++ random yields different numbers for same Mersenne Twister seed when using float precision

833 Views Asked by At

I need to run reproducible Monte Carlo runs. That means I use a known seed that I store with my results, and use that seed if I need to run the same problem instance using the same random numbers. This is common practice.

While investigating the effects of numeric precision, I ran into the following issue: For the same Mersenne Twister seed, std::uniform_real_distribution<float>(-1, 1) returns different numbers than std::uniform_real_distribution<double>(-1, 1) and std::uniform_real_distribution<long double>(-1, 1), as the following example shows:

#include <iomanip>
#include <iostream>
#include <random>

template < typename T >
void numbers( int seed ) {
  std::mt19937                        gen( seed );
  std::uniform_real_distribution< T > dis( -1, 1 );
  auto p = std::numeric_limits< T >::max_digits10;
  std::cout << std::setprecision( p ) << std::scientific << std::setw( p + 7 )
            << dis( gen ) << "\n"
            << std::setw( p + 7 ) << dis( gen ) << "\n"
            << std::setw( p + 7 ) << dis( gen ) << "\n"
            << "**********\n";
}

int main() {
  int seed = 123;
  numbers< float >( seed );
  numbers< double >( seed );
  numbers< long double >( seed );
}

Result:

$ /usr/bin/clang++ -v
Apple LLVM version 10.0.0 (clang-1000.11.45.5)
Target: x86_64-apple-darwin18.2.0
Thread model: posix
InstalledDir: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin

$ /usr/bin/clang++ bug.cpp -std=c++17
$ ./a.out 
 3.929383755e-01
 4.259105921e-01
-4.277213216e-01
**********
 4.25910643160561708e-01
-1.43058149942132062e-01
 3.81769702875451866e-01
**********
 4.259106431605616525145e-01
-1.430581499421320209545e-01
 3.817697028754518623166e-01
**********

As you can see, double and long double both start around at the same number (save precision differences) and continue yielding the same values. On the other hand, float starts off with a completely different number, and its second number is similar to the first number produced by double and long double.

Do you see the same behavior in your compiler? Is there a reason for this unexpected (to me) discrepancy?

Approach

The responses make it clear that there's no reason to expect that values generated with different underlying precision will be the same.

The approach that I'll take to generate reproducible runs will be to always generate values at the highest precision possible, and cast them to lower precision on demand (e.g., float x = y, where y is double or long double, as the case may be).

3

There are 3 best solutions below

3
On BEST ANSWER

Each distribution will generate floating point numbers by grabbing a sufficient number of (pseudo)random bits from the underlying Mersenne Twister and then producing uniformly distributed floating point numbers from it.

There are only two ways that an implementation could fulfill your expectation of "same algorithm, therefore same results (minus precision)":

  1. std::uniform_real_distribution<long double>(-1, 1) is only as random as std::uniform_real_distribution<float>(-1, 1). More to the point, the former has exactly as many possible outputs as the latter. If the latter can produce more different values than the former, then it needs to consume more bits of randomness from the underlying Mersenne Twister. If it cannot - well, what's the point of using it (and how would it still be "uniform")?

  2. std::uniform_real_distribution<float>(-1, 1) consumes (and mostly discards) exactly as many bits of randomness from the underlying Mersenne Twister as std::uniform_real_distribution<long double>(-1, 1). That would be very wasteful and inefficient.

Since no sane implementation will do either of the above, std::uniform_real_distribution<long double>(-1, 1) will advance the underlying Mersenne Twister by more steps than std::uniform_real_distribution<float>(-1, 1) for each generated number. That will of course change the progression of the random numbers. This also explains why the long double and double variant are relatively close together: They share most of their random bits initially (whereas float likely requires a lot fewer bits and thus diverges quicker).

0
On

Initializing a random number generator to a specific seed will specify the sequence of random bits it puts out. However, you aren't using those bits the same way in each case. A std::uniform_real_distribution<double> has a larger possibility space than std::uniform_real_distribution<float> (assuming sizeof(double) > sizeof(float) on your platform) so it will need to consume a larger quantity of random bits to generate a fully uniform distribution.

The first consequence is that the pseudo-random sequence of bits will have a different interpretation for different distribution types. The second consequence is each distribution moves a different number of bits down the pseudo-random sequence whenever is produces a value, meaning following numbers won't be at the same point in the pseudo-random bit sequence.

The solution to your problem is to always use the same type of distribution. If you want to compare the result of using lower precision values with using higher precision values, only generate values with the highest precision and truncate them down when needed.

0
On

Just to add to excellent @MaxLanghof answer with more specifics:

For double code would do something like this - generate u64 integer, and use 53bits of it to make float, along the lines

double r = (u64 >> 11) * (1.0 / (uint64_t(1) << 53));

For long double, assuming Intel 80bits format, with 64bit mantissa, it will do about the same, get 64bits, return back long double.

long double r = u64 * (1.0 / (uint64_t(1) << 64)); // pseudocode

64bit of randomness are consumed in both cases, thus you see same values.

In case of a float, 32bits are used to make single float

float r = (u32 >> 8) * (1.0f / (uint32_t(1) << 24));

32bits of randomness are consumed, and another 32bits are used for next number, which together with endianness makes second float about the same as first double/long double.

Link: http://xoshiro.di.unimi.it/