为 64 位 LCG (MMIX (by Knuth)) 找到更多独立种子值

Question

我使用的是 64 位 LCG（MMIX（作者 Knuth））。它在我的代码中生成了一定的随机数块，这些随机数使用它们来执行一些操作。我的代码在单核中工作，我想并行化工作以减少执行时间。
在开始考虑这种意义上更高级的方法之前，我想简单地并行执行更多相同的代码（事实上，代码在一定数量的独立模拟上重复相同的任务，所以我可以简单地将模拟数量分成更多相同的代码和运行它们并行）。
我现在唯一的问题是为每个代码找到一个种子；特别是，为了避免在不同代码中生成的数据之间出现不需要的非平凡相关性的可能性，我必须确保在各种代码中生成的随机数不会重叠。为此，从第一个代码中的某个种子开始，我必须找到一种方法来找到一个非常远的值（下一个种子），不是绝对值而是伪随机序列（因此，从从第一到第二种子，我需要大量的LCG步数。
我的第一次尝试是这样的：
从序列
中生成的2个连续数字之间的LCG关系开始
因此，原则上，我可以计算上述关系，例如，n = 2^40 和 I_0 等于第一个种子的值，并在随机 CLG 序列中获得距离 2^40 步的新种子从第一个开始。
问题是，这样做，我必须进行溢出计算 a^n。事实上，对于 MMIX（Knuth）a~2^62，我使用 unsigned long long int（64 位）。请注意，这里唯一的问题是上述关系中的分数。如果只有求和和乘法，我可以忽略由于以下 modular 属性引起的溢出问题（实际上我使用 2^64 作为 c（64 位生成器））：

那么，从某个值（第一个种子）开始，我怎样才能在 LC 伪随机序列中找到距离很大的第二个值？

[编辑]
r3mainer 解决方案非常适合 python 代码。我现在正在尝试使用无符号 __int128 变量在 c 中实现它。我只有一个问题：原则上我应该计算：

比如说，为了简单起见，我想计算：

n = 2^40 和 c(a-1)~2^126。我继续 cycle.Starting 和 temp = a，在每次迭代中我计算 temp = temp*temp，然后我计算 temp%c(a-1)。问题出在第二步（temp = temp*temp）。 temp 实际上原则上可以是 < c(a-1)~2^126 的任何数字。如果 temp 是一个大数字，比如 > 2^64，我将溢出，在下一个 module 操作之前达到 2^128 - 1。那么有没有办法避免呢？目前，我看到的唯一解决方案是按照此处的建议通过位循环执行每个乘法：c code: prevent overflow in modular operation with huge modules (modules near the overflow treshold) 在乘法期间是否有另一种方法来执行 module 运算？
（请注意，c = 2^64，使用 mod(c) 操作我没有同样的问题，因为溢出点（对于 ull int 变量）与 module 重合）

Answer 1

x[n+1] = (x[n] * a + c) % m 形式的任何 LCG 都可以非常快速地跳到任意位置。

从零种子值开始，LCG 的前几次迭代将为您提供以下序列：

x₀ = 0
x₁ = c % m
x₂ = (c(a + 1)) % m
x₃ = (c(a² + a + 1)) % m
x₄ = (c(a³ + a² + a + 1)) % m

很容易看出每一项实际上是一个几何级数的总和，可以用simple formula:

来计算

x_n = (c(a^{n-1} + a^{n-2} + ... + a + 1)) % m
    = (c * (a^n - 1) / (a - 1)) % m

(a^n - 1)项可以直接通过(a-1) mod m的modular exponentiation, but dividing by (a-1) is a bit tricky because (a-1) and m are both even (i.e., not coprime), so we can't calculate the modular multiplicative inverse快速计算。

相反，计算 (a^n-1) mod m*(a-1)，然后将结果直接（非模块化）除以 a-1。在 Python 中，计算将如下所示：

def lcg_skip(m, a, c, n):
    # Calculate nth term of LCG sequence with parameters m (modulus),
    # a (multiplier) and c (increment), assuming an initial seed of zero
    a1 = a - 1
    t = pow(a, n, m * a1) - 1
    t = (t * c // a1) % m
    return t

def test(nsteps):
    m = 2**64
    a = 6364136223846793005
    c = 1442695040888963407
    #
    print("Calculating by brute force:")
    seed = 0
    for i in range(nsteps):
        seed = (seed * a + c) % m
    print(seed)
    #
    print("Calculating by fast method:")
    # Calculate nth term by modular exponentiation
    print(lcg_skip(m, a, c, nsteps))

test(1000000)

因此，要创建具有非重叠输出序列的 LCG，您需要做的就是使用 lcg_skip() 生成的初始种子值与 n 的值相距足够远。

Answer 2

好吧，对于 LCG 来说，已知属性在 O(log₂(N)) 时间内向前和向后跳跃，其中 N 是跳跃点之间的距离，F. Brown 的论文，“任意步幅的随机数生成”，Trans。是。核素。社会。（1994 年 11 月）。

这意味着如果你有满足赫尔-多贝尔定理的LCG参数（a，c），那么整个周期将是2⁶⁴个数字在重复之前，并说Nt number pf threads 你使跳跃距离为 2⁶⁴ / Nt，并且所有线程都以相同的种子开始，并且在通过 (2⁶⁴ / Nt)*threadId，并且由于序列重叠，您将完全免受 RNG 相关性的影响。

对于在 NumPy 中实现的常见 64 位无符号模数学的最简单情况，下面的代码应该可以正常工作

import numpy as np

class LCG(object):

    UZERO: np.uint64 = np.uint64(0)
    UONE : np.uint64 = np.uint64(1)

    def __init__(self, seed: np.uint64, a: np.uint64, c: np.uint64) -> None:
        self._seed: np.uint64 = np.uint64(seed)
        self._a   : np.uint64 = np.uint64(a)
        self._c   : np.uint64 = np.uint64(c)

    def next(self) -> np.uint64:
        self._seed = self._a * self._seed + self._c
        return self._seed

    def seed(self) -> np.uint64:
        return self._seed

    def set_seed(self, seed: np.uint64) -> np.uint64:
        self._seed = seed

    def skip(self, ns: np.int64) -> None:
        """
        Signed argument - skip forward as well as backward

        The algorithm here to determine the parameters used to skip ahead is
        described in the paper F. Brown, "Random Number Generation with Arbitrary Stride,"
        Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in
        O(log2(N)) operations instead of O(N). It computes parameters
        A and C which can then be used to find x_N = A*x_0 + C mod 2^M.
        """

        nskip: np.uint64 = np.uint64(ns)

        a: np.uint64 = self._a
        c: np.uint64 = self._c

        a_next: np.uint64 = LCG.UONE
        c_next: np.uint64 = LCG.UZERO

        while nskip > LCG.UZERO:
            if (nskip & LCG.UONE) != LCG.UZERO:
                a_next = a_next * a
                c_next = c_next * a + c

            c = (a + LCG.UONE) * c
            a = a * a

            nskip = nskip >> LCG.UONE

        self._seed = a_next * self._seed + c_next


#%%
np.seterr(over='ignore')

seed = np.uint64(1)

rng64 = LCG(seed, np.uint64(6364136223846793005), np.uint64(1))

print(rng64.next())
print(rng64.next())
print(rng64.next())

#%%
rng64.skip(-3) # back by 3
print(rng64.next())
print(rng64.next())
print(rng64.next())

rng64.skip(-3) # back by 3
rng64.skip(2) # forward by 2
print(rng64.next())

在 Python 3.9.1、x64 Win 10 中测试

为 64 位 LCG (MMIX (by Knuth)) 找到更多独立种子值

find more indipendent seed value for a 64 bit LCG (MMIX (by Knuth))

random

parallel-processing

numbers

overflow

random-seed