規格與設計考量

std::uniform_int_distribution 的規格是：給定一個能產生 [0, 2^N) 區間整數的亂數生成器，產生一個介於 [min, max] 區間的整數。區間 [0, 2^N) 到區間 [min, max] 的映射必須保證每個 [min, max] 區間內每個整數的機率是均等的。

計算過程中 std::uniform_int_distribution 可以捨棄亂數生成器先前提供的亂數，再次從亂數生成器提取新的亂數。舉例來說，如果亂數生成器可以產生 [0, 2^32) 的整數，而 std::uniform_int_distribution 的值域是 [1, 6]，一個簡單的實作是捨棄所有不屬於 [1, 6] 區間的亂數。然而，一個良好的映射應該要儘可能避免捨棄亂數，畢竟偽亂數生成器也需要計算時間。以 Entropy Pool 為基礎的偽亂數生成器，更需等待作業系統蒐集足夠的 Entropy。

另一個考量是映射本身需要的計算。一些數值模擬程式會需要大量的亂數，因此這個映射必須越快越好。本文後兩個演算法就是以不同的方法節省餘數運算子（Modulo, %）。

OpenBSD 演算法

第一個演算法的概念非常簡單：令 s 為 max - min + 1，將 [0, 2^N) 從後面開始以 s 為單位分段，如果亂數屬於這些段落，我們就回傳 x % s。如果最前面一段的長度不足 s，我們就捨棄該區間內的亂數。另外，為了使用 uint64_t 計算 2^64 % s，這裡使用一個恆等式：2^64 % s == (2^64 - s) % s。完整程式碼如下：

class uniform_int_distribution_openbsd {
private:
  const uint64_t min_value;
  const uint64_t s;

public:
  uniform_int_distribution_openbsd(uint64_t min, uint64_t max)
      : min_value(min), s(max - min + 1) {}

  template <typename UniformRandomBitGenerator>
  uint64_t operator()(UniformRandomBitGenerator &urbg) const {
    // 2^64 % s = (2^64 - s) % s = uint64_t(-s) % s
    uint64_t discard_bound = static_cast<uint64_t>(-s) % s;
    uint64_t x = urbg();
    while (x < discard_bound) {
      x = urbg();
    }
    return x % s + min_value;
  }
};

如果亂數生成區間是 [0, 2^3) 想要映射的值域是 [0, 2]，結果對照表如下：

備註：論文作者 Daniel 將這個演算法稱為 OpenBSD 演算法，因為他可以在 OpenBSD 的原始碼庫找到這個演算法。雖然我沒有去找更早的出處，我覺得應該會有更早的實作。

OpenJDK 演算法

第二個演算法是 OpenJDK 裡面的演算法。一樣是將 [0, 2^N) 區間以 s 為單位分段。不同的是它是從前面開始分組。如果亂數屬於最後一組，我們就捨棄該亂數。與第一個演算法不同，這個方法以 x - (x % s) > 2^N - s 作為是否為最後一組的判斷式。其思路是：x % s 是計算結果必須使用的運算子。我們可以透過 x - (x % s) 得到每一個分段的開頭，如果開頭大於 2^N - s，該段沒有完整的 [0, s) 區間，因此捨棄該亂數。

這個演算法的優點是：如果 s 相對於 2^N 夠小，這個演算法有很大的機會不需要進入迴圈，進而只需要一個餘數運算子。這個演算法的缺點是：在最壞情況下這個演算法可能需要「多個」餘數運算子，餘數運算子的數量不是一個常數。與此對照，本文開頭的 OpenBSD 演算法可以保證在最差情況下永遠只需要兩個餘數運算子。

OpenJDK 演算法完整程式碼如下：

class uniform_int_distribution_openjdk {
private:
  const uint64_t min_value;
  const uint64_t s;

public:
  uniform_int_distribution_openjdk(uint64_t min, uint64_t max)
      : min_value(min), s(max - min + 1) {}

  template <typename UniformRandomBitGenerator>
  uint64_t operator()(UniformRandomBitGenerator &urbg) const {
    uint64_t x = urbg();
    uint64_t r = x % s;
    while (x - r > static_cast<uint64_t>(-s)) {  // 2^64 - s
      x = urbg();
      r = x % s;
    }
    return r + min_value;
  }
};

如果亂數生成區間是 [0, 2^3) 想要映射的值域是 [0, 2]，結果對照表如下：

Daniel Lemire 的演算法

最後是論文作者 Daniel Lemire 提出的演算法。Daniel 另闢蹊徑以 128-bit 整數乘法將 [0, 2^N) 映射到 [0, s * 2^N)，接著再將 128-bit 整數分為上半部分（即除以 2^64）與下半部分（即取 2^64 的餘數）。上半部分會介於 [0, s)，如果滿足其他條件，上半部分就會是回傳值。下半部分如果小於 s 我們要再做額外的判斷。

這個方法會將 [0, 2^N) 分成 s 個區間。如果 s 能夠整除 2^N，每個區間會有 floor(2^N / s) 個數值。如果不能整除，前 2^N % s 個區間會多分到 1 個數值。如果我們將每個區間的第一個數字作為被捨棄的對象，我們能比較下半部與 2^N % s 之間的關係。如果下半部小於 2^N % s，該亂數就屬於前 2^N % s 個區間的第一個數字，因此捨棄該亂數。

另外，因為 2^N % s < s，我們可以先以 lower < s 作為快篩。只有在 lower < s 的時候才執行餘數運算子。和 OpenBSD 演算法相同，這個餘數運算子的結果可以被重覆使用，因此 Daniel Lemire 演算法最差情況下只會使用一個餘數運算子。

但是這個演算法帶來一個問題：128-bit 整數乘法會比較好嗎？至截稿前 128-bit 整數乘法並不是 C/C++ 標準內建的整數型別。GCC 或 Clang 編譯器各自定義了 128-bit 整數。我們能透過 __uint128_t 做 128-bit 整數乘法。許多 64 位元計算機結構（例如：x86-64、aarch64、riscv64）都定義了計算乘法上半部結果的指令（例如：mulh），因此實務上 128-bit 整數乘法會比餘數運算子快。

完整程式碼如下：

class uniform_int_distribution_lemire {
private:
  const uint64_t min_value;
  const uint64_t s;

public:
  uniform_int_distribution_lemire(uint64_t min, uint64_t max)
      : min_value(min), s(max - min + 1) {}

  template <typename UniformRandomBitGenerator>
  uint64_t operator()(UniformRandomBitGenerator &urbg) const {
    uint64_t x = urbg();
    __uint128_t m = static_cast<__uint128_t>(x) * s;
    uint64_t lower = static_cast<uint64_t>(m);
    if (lower < s) {
      // (2^64 - s) % s = (2^64 - s) % s = uint64_t(-s) % s
      uint64_t discard_bound = static_cast<uint64_t>(-s) % s;
      while (lower < discard_bound) {
        x = urbg();
        m = static_cast<__uint128_t>(x) * s;
        lower = static_cast<uint64_t>(m);
      }
    }
    uint64_t upper = static_cast<uint64_t>(m >> 64);
    return upper + min_value;
  }
};

如果亂數生成區間是 [0, 2^3) 想要映射的值域是 [0, 2]，結果對照表如下：

後記

寫作本文的動機起源於我對 Abseil 亂數函式庫的不滿。Abseil 亂數函式庫會刻意打亂亂數種子，讓同一個程式、同一個亂數種子在每次執行時都會產生不同的亂數序列。Abseil 亂數函式庫的作者明確表示：他們不會提供任何可以關閉這個行為的 API。因為如果他們不打亂亂數種子，函式庫的使用者就會開始依賴特定序列，這會讓他們無法改進亂數函式庫。因為我就是想要固定的序列，這讓我非常不滿。然而，通過一些研究之後，我才知道 C++ 標準也沒有規定 std::uniform_int_distribution 的結果。C++ 標準只要求 std::uniform_int_distribution 的結果必須滿足均勻分佈的性質。接著我找到 Daniel Lemire 的論文，看到各種不同的最佳化。我對 Abseil 亂數函式庫的看法才由不滿轉向認同。我覺得這種「不規定結果，只規定結果的性質」的 API 設計非常特別，因此花了一點時間把想法記錄下來。

參考資料

• Daniel Lemire, Fast Random Integer Generation in an Interval
• The Abseil Random Library, Stability of Generated Sequences