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?
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).
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)":
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")?
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).