当使用float精度时，相同的Mersenne Twister种子在C++随机数生成器中产生不同的数字。

Question

我需要运行可重复的蒙特卡罗模拟。这意味着我使用已知的种子并将其存储在我的结果中，如果我需要使用相同的随机数字运行相同的问题实例，我会使用该种子。这是常见的做法。

在研究数字精度的影响时，我遇到了以下问题：对于相同的Mersenne Twister种子，std::uniform_real_distribution<float>(-1, 1)返回的数字与std::uniform_real_distribution<double>(-1, 1)和std::uniform_real_distribution<long double>(-1, 1)不同，如下面的示例所示：

#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 );
}

结果：

$ /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
**********

如您所见，double和long double从同一数字开始（除了精度的差异），并继续产生相同的值。另一方面，float以完全不同的数字开始，其第二个数字类似于double和long double生成的第一个数字。

您的编译器是否有相同的行为？这种意外的（对我来说）差异有原因吗？

方法

回答表明没有理由期望使用不同底层精度生成的值相同。

我将采取的方法是始终以最高可能的精度生成值，并根据需要将它们转换为较低的精度（例如，float x = y，其中y是double或long double，视情况而定）。

- Escualo

---

你为什么期望它们生成相同的数字？ - Max Langhof

只要相同的种子导致相同的“float”的序列，我认为就没有问题。 - Hong Ooi

如果float和double的精度不同，那么对它们应用相同的变换可能会导致由于舍入误差而产生不同的结果。只要每种类型始终为相同的种子提供相同的序列，您的实现就是符合标准的。 - NathanOliver

但这是一种不同的输出类型。 当然，您不会期望整数上的均匀随机产生与双精度浮点数相同的值（这显然是不可能的，除非要求所有随机双精度浮点数都是整数...）。 对于float和double也是同样的情况。 - Max Langhof

@MaxLanghof 当然可以，但是令我意外的是，float和double之间的差异不仅在于精度（应用于浮点数类型的相同算法）。 - Escualo

Answer 1

每个分布都会通过从底层Mersenne Twister中获取足够数量的（伪）随机位，然后从中生成均匀分布的浮点数。

只有两种方法能够满足您对“相同算法，因此结果相同（减去精度）”的期望：

1. `std::uniform_real_distribution(-1, 1)` 的随机性仅与 `std::uniform_real_distribution(-1, 1)` 一样。更重要的是，前者与后者具有完全相同数量的可能输出。如果后者可以产生比前者更多的不同值，则需要从底层Mersenne Twister消耗更多的随机位。如果它不能 - 那么使用它的意义何在（它如何仍然“均匀”）？

2. `std::uniform_real_distribution(-1, 1)` 恰好会从底层Mersenne Twister中消耗（并大多丢弃）与 `std::uniform_real_distribution(-1, 1)` 相同数量的随机位。这将非常浪费和低效。

由于没有合理的实现会执行上述任何一种操作，`std::uniform_real_distribution(-1, 1)` 将为每个生成的数字推进底层Mersenne Twister的步骤更多。这当然会改变随机数的进展。这也解释了为什么 `long double` 和 `double` 变体相对接近：它们最初共享大部分随机位（而浮点数可能需要更少的位，因此更快地发散）。

- Max Langhof

---

非常合理。我想我搞错了（反过来）：double和long double是“相同”的事实是巧合（由于double的基础更高精度）。一般情况下，我不应该期望得到相同的结果。 - Escualo

我不确定将均匀生成的 long double 截断为 float 是否会得到均匀生成的 float，因此 2 可能无法符合要求。而 1 明显不行。 - Caleth

@Caleth 对于第二种变体，我的意思是它会丢弃那些位，而不是将它们用于扩展精度和截断。换句话说，我是指所有分布都会为了这个目的而推进它们的 RNG 状态相同的量（而不会改变它们实际使用的位数）。 - Max Langhof

Answer 2

初始化一个随机数生成器到一个特定的种子将指定它输出的随机位序列。然而，你在每种情况下使用这些位的方式不同。如果你的平台上sizeof(double)>sizeof(float)，那么std::uniform_real_distribution比std::uniform_real_distribution有更大的可能空间，因此它将需要消耗更多的随机位来生成完全均匀的分布。

第一个后果是伪随机位序列对于不同的分布类型将有不同的解释。第二个后果是每个分布在产生一个值时向下移动不同数量的位，意味着后续数字不会处于伪随机位序列中相同的位置。

解决问题的方法是始终使用相同类型的分布。如果你想比较使用低精度值和使用高精度值的结果，请只生成具有最高精度的值，并在需要时将它们截断。

- François Andrieux

Answer 3

补充@MaxLanghof的优秀答案，提供更多细节信息：

对于双倍精度浮点数，可以按照以下方式生成u64整数，并使用其中的53位来生成浮点数：

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

对于长双精度浮点数，假设采用Intel 80位格式，具有64位的尾数，它将执行大致相同的操作，获取64位，返回长双精度浮点数。

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

两种情况下都使用了64位的随机数，因此您会看到相同的值。

对于浮点数，使用32位来制作单个浮点数。

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

将消耗32位的随机数，另外32位用于下一个数字生成。由于字节顺序，第二个浮点数与第一种双精度/长双精度浮点数大致相同。

链接：http://xoshiro.di.unimi.it/

- Severin Pappadeux