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

【问题标题】：find more indipendent seed value for a 64 bit LCG (MMIX (by Knuth))为 64 位 LCG (MMIX (by Knuth)) 找到更多独立的种子值
【发布时间】：2021-02-12 09:05:14
【问题描述】：

我使用的是 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 位）。请注意，这里唯一的问题是上述关系中的分数。如果只有求和和乘法，由于以下模块化属性，我可以忽略溢出问题（实际上我使用 2^64 作为 c（64 位生成器））：

那么，从某个值（第一个种子）开始，如何在 LC 伪随机序列中找到第二个远离大量步长的值？

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

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

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

【问题讨论】：

标签： random parallel-processing numbers overflow random-seed

【解决方案1】：

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

从种子值 0 开始，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) 项可以通过modular exponentiation 快速计算，但除以(a-1) 有点棘手，因为(a-1) 和m 都是偶数（即不是互质数），所以我们不能直接计算(a-1) mod m的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 值。

【讨论】：

谢谢 r3mainer；有了这个，我只需要一个 128 位变量（存储 m * a1）。所以我可以继续使用一个简单的 c 函数，该函数使用 /(m * a1) 步骤迭代地执行幂幂运算以获取模数。你是在 python 中完成的，所以在 python 中你也可以处理 128 位（或者更大的）变量？
@user1172131 Python 和 Ruby 都内置了任意精度的整数运算，所以是的，您可以轻松处理超过 128 位的问题。
@user1172131 不，一点也不。 lcg_skip(2**64, 6364136223846793005, 1442695040888963407, 2**40) 立即产生结果。
@user1172131 如果你没有 Python，你可以在这里生成你自己的种子值：ideone.com/4ECnJE
@user1172131 在 Python v3.0 及更高版本中，/ 运算符为您提供浮点结果。由于双精度变量只有 53 位尾数（如果有记忆的话），它对 128 位整数没有用处。 // 整数除法运算符在这里是必不可少的。在任何情况下，您都在计算整数结果。在任何情况下，您都不希望 lcg_skip() 返回浮点值。是的，除法保证没有余数。您可以按照答案中给出的链接了解更多关于模幂运算如何工作的信息。

【解决方案2】：

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

这意味着如果你的 LCG 参数 (a, c) 满足 Hull-Dobell 定理，那么整个周期将是 2⁶⁴ 个数字，然后再重复它们，并说 Nt number pf threads you使跳转距离为 2⁶⁴ / Nt，所有线程都以相同的种子开始，并在通过 (2⁶⁴ / Nt)*threadId 初始化 LCG 后跳转，你会由于序列重叠，因此完全不受 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 中测试

【讨论】：