查找64位LCG的更多独立种子值(MMIX ( Knuth))

Question

我正在使用64位LCG (MMIX (由Knuth))。它在我的代码中生成一个特定的随机数块，这些数字用来执行一些操作。我的代码工作在单一核心，我想并行工作，以减少执行时间。

在开始考虑更高级的方法之前，在这个意义上，我想简单地并行执行更多相同的代码(实际上，这些代码在一定数量的独立模拟上重复相同的任务，所以我可以简单地将模拟的数量平分到更相同的代码中，并并行运行)。

我现在唯一的问题是为每个代码找到一个种子；特别是，为了避免在不同代码中生成的数据之间出现不必要的非平凡相关性，我必须确保在不同代码中生成的随机数不重叠。要做到这一点，从第一段代码中的某个种子开始，我必须找到一种方法来找到一个值(下一个种子)非常遥远，不是绝对值，而是伪随机序列(因此，要从第一个种子到第二个种子，我需要大量的LCG步骤)。

我的第一次尝试是：

从序列中生成的两个连续数之间的LCG关系出发

所以，原则上，我可以用n= 2^40和I_0等于第一个种子的值来计算上面的关系，并从第一个随机CLG序列中得到一个新的种子间距2^40步。

问题是，这样做，我需要溢出计算a^n。事实上，对于MMIX (由Knuth) a^62，我使用无符号长int (64位)。请注意，这里唯一的问题是上面关系中的分数。如果只有求和和乘法，由于以下模块属性，我可以忽略溢出问题(实际上，我使用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，那么在下一个模块操作之前，我将进入溢出状态，达到2^128-1。那么有什么办法可以避免呢？现在，我看到的唯一解决方案是用一个循环对位执行每个乘法，就像这里所建议的那样:在乘法过程中，C代码:防止使用大型模块的模块化操作中的溢出(模块位于溢出窗口附近)还有执行模块操作的另一种方法吗？

(注意，c= 2^64，对于mod(c)操作，我没有同样的问题，因为溢出点(对于ull int变量)与模块重合)

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

很容易看出，每个项实际上是一个几何级数的和，可以用简单公式来计算。

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

(a^n - 1)项可以用模指数快速计算，但是除以(a-1)是有点棘手的，因为(a-1)和m都是偶数的(也就是不对等的)，所以不能直接计算(a-1) mod m的模乘逆。

相反，计算(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值。