哪些C++随机数引擎具有 O(1) 丢弃功能

Which C++ random number engines have a O(1) discard function?

本文关键字:功能 随机数 C++ 引擎 哪些      更新时间:2023-10-16

自C++11以来,有许多std随机数引擎。他们实现的成员函数之一是跳过 z 随机生成的数字的void discard(int long long z)。此函数的复杂度在 www.cplusplus.com (http://www.cplusplus.com/reference/random/mersenne_twister_engine/discard/) 上为 O(z) 给出

但是,在 www.cppreference.com(http://en.cppreference.com/w/cpp/numeric/random/mersenne_twister_engine/discard)上有一个注释说

对于某些引擎,"快速跳跃"算法是已知的,它推进了 状态通过许多步骤(数百万的数量级)而不计算 中间状态转换。

我如何知道哪些发动机的实际丢弃成本是 O(1)?

好吧,如果你使用预先计算的跳跃点,O(1) 将适用于现有的每个 RNG。请记住,有些算法可能比 O(z) 更好,但不是 O(1) - 比如 O(log2z)。

如果我们谈论跳到任意点,事情就会变得有趣。例如,对于线性同余生成器,有已知的O(log2z)跳跃提前算法,基于F. Brown的论文,"Random Number Generation with Arbitrary Stride",Trans. Am. Nucl。Soc. (1994年11月)。代码示例在这里。

C++11标准中有LCG RNG,不确定在特定实现中跳跃的速度有多快(http://en.cppreference.com/w/cpp/numeric/random/linear_congruential_engine)

我相信 PCG 系列的 RNG 共享相同的属性

事实是,std::linear_congruential_engine<UIntType,a,c,m>::discard(unsigned long long z)绝对可以非常有效地实现(实际上非常微不足道)。它主要等同于a的幂zm的幂(对于零和非零c) - 这意味着非常简单的软件实现以O(log(z % φ(m)))UIntType乘法+模约运算执行(φ(m)=m-1- 素数和一般情况下小于m)。

^^^请注意,在某种程度上,O(log(z % φ(m)))O(1)是因为log2(z % φ(m)) < log2(m) < sizeof(UIntType)*CHAR_BIT- 尽管在实践中它更像是O(log(z))

^^^另请注意,在生成一些依赖于模板参数的预计算表(如果合适,可以使用延迟或编译时计算)后,幂的复杂性可以降低到只有几个multiplication + modulo reduction运算(如 4 或 8) - 即O(1)在任何实际意义上

。此外,对于大多数其他引擎的discard函数,可能有足够有效的算法来满足O(P(sizeof(state))*log(z))约束(P(x)- 一些低次多项式函数,最有可能是1+o(1)度或最多3+o(1)度,因为log(z) < sizeof(unsigned long long)*CHAR_BIT可以被认为是常数)。

尽管如此:

不知何故由于某种未知的原因C++标准(截至 ISO/IEC 14882:2017)不需要以更有效的方式实施discard,而不仅仅是zoperator()()对任何 PRNG 引擎的调用,包括那些绝对允许这样做的引擎。

就我个人而言,这令人困惑,毫无意义-它违反了C++语言设计的基本原则之一(以一种非常残酷的方式)- 即在性能和实际用途方面仅在C++标准中添加合理的功能。

值得注意且非常有据可查的例子:没有通过索引(std::list<T>::operator[](size_type n))访问std::list元素这样的事情,即使它就像使用迭代器begin()调用operator++()n时间一样"简单"。自然如此 - 因为执行时间O(n)在任何实际应用中都会使这个函数成为不合理的选择(AKA愚蠢的想法)。出于这个显而易见的原因,a[n]a.at(n)不是强制性序列容器要求(ISO/IEC 14882:2017 26.2.3 表 87)的一部分,而是可选序列容器操作的一部分(ISO/IEC 14882:2017 26.2.3 表 88)。

那么,为什么世界上e.discard(z)是强制性随机数引擎要求(ISO/IEC 14882:2017 29.6.1.4 表 104)的一部分,具有这种荒谬的复杂性要求 -no worse than the complexity of z consecutive calls e()- 而不是一些具有足够复杂性要求的可选操作部分条目,如O(P(sizeof(state))*log(z))

什么...?z consecutive calls e()是指数复杂性 -从什么时候开始可以?

更令人困惑的是,在我的GCC中实际上发现了这个现实世界的实现:

void discard(unsigned long long __z)
{
for (; __z != 0ULL; --__z)
(*this)(); //<--  Wait what? Are you kidding me?
}

因此,再一次像以前一样,我们别无选择,只能自己实现必要的功能......标准库C++帮助不大。

修正案:

在深入研究Mersenne Twister的设计细节时,我们发现std::mersenne_twister_engine的discard(z)也可以非常有效地实现。

template<
class UIntType,
std::size_t w, std::size_t n, std::size_t m, std::size_t r, UIntType a,
std::size_t u, UIntType d, std::size_t s, UIntType b, std::size_t t, UIntType c,
std::size_t l, UIntType f
> class mersenne_twister_engine;

即使是discard(z)的通用实现(适用于整类线性 PRNGs 模 2 - 不仅是 Mersenne Twister,还有 WELL 和许多其他)也会像O(n^ω*log(z))一样复杂 - 其中n是上面的模板参数 -w位字的状态大小,幂ω在 2 和 3 之间恒定(取决于所选的位矩阵乘法算法)。通过合理数量的模板参数相关预计算来O(n^ω)实际执行时间,可以稍微降低这种复杂性。SIMD CPU 指令(矢量 XOR 或矢量 AND)将按恒定系数提高实际性能。并行算法(例如专用硬件解决方案)将使用O(n^3)同步单比特计算(XOR 和 AND)在O(log(n))时间内计算这一点。

当然,您可能会注意到上面的n参数通常不是那么小(例如 624 表示std::mt19937或 312 表示std::mt19937_64),n-cubed 甚至更大 - 因此O(n^3)不一定是快速的实际执行。但是现代CPU(尤其是GPU)仍然可以非常快速地执行优化的实现 - 无法与z consecutive calls of operator()()荒谬的指数复杂性相提并论。

一些一般性意见:

每个现有的PRNG(以及我能想象到的每个PRNG)都可以用以下迭代方程定义:

x[n+1] = T(x[n])
Output[n] = F(x[n])

其中x[n]是经过n次迭代后的状态(一些W位序列)(所以x[0]是初始状态AKA种子),T(x)是状态迭代转换函数(将当前状态转换为下一个状态),F(x)是输出转换,它将每个状态W位序列转换为输出v位序列(Output[n]) - 通常v < W

当然,T(x)F(x)计算都很快 - 即最大时间是多项式 -O(P(W))最坏的情况。(通常,这些功能被设计得更快 - 例如O(P(v))在大多数情况下基本上是O(1),因为通常选择v作为CPU寄存器大小,并且通常可用于该大小的快速硬件优化操作)。

我的意思是字面意思 - 所有有意义的现有(和未来)PRNG 都可以以这种方式表达。

(我能想到的唯一进一步的概括是使Wv大小不变 - 即依赖于n- 即从一个迭代到下一个迭代变化。这可能没有太多的实际意义 - 我想没有人希望他们的 PRNG 内部数据无限增长并最终消耗所有 RAM 或类似的东西。即使增长非常缓慢W也可以允许非周期性PRNG设计。

所以问题是:PRNG的什么属性会使discard(z)操作快速- 即多项式 -O(P(W))- 最坏情况下的运行时间?

IMO 非常明显的答案是,如果我们可以对函数执行快速计算 -T^z(x) = T(T(...T(x)...)) - z times对于任何z,那么我们可以实现快速discard(z).

此外,不难注意到,IFT(x) = T_p(x)是一些具有一些内部固定参数的参数化变换p它是具有不同参数值的一类变换之一,对于任何允许的参数值q变换T_q(x)都可以在O(P(W))时间内快速计算。此外,如果对于任何允许的参数值pq转换,T_p(T_q(x))也属于此类转换,具有一些允许的参数r- 即T_p(T_q(x)) = T_r(x)r可以通过参数快速计算pq...假设我们定义了一个符号r=p*q- 其中*是一些二进制运算可快速计算(最多在多项式时间内) - 所以根据定义T_{p*q}(x) = T_p(T_q(x))。(你可以注意到二元运算*不一定是可交换的——即p*q的值不得与q*p相同。但是此操作在设计上是关联的- 因为T_p(T_q(T_r(x))) = T_p(T_{q*r}(x)) = T_{p*q}(T_r(x))- 因此p*(q*r)=(p*q)*r.)

^^^这种T(x)转换结构显然允许快速计算T^z(x)转换:如果已知T(x) = T_p(x)和参数p- 我们只计算q=p^z=p*p*p*...p - z times(这只是O(log(z))*的计算,可以通过预计算和/或并行执行进行优化),然后我们计算T_q(x)

^^^ 虽然这许多先决条件看起来非常特殊 - 事实上所有这些都很常见。 例如,对于一类线性变换模 2(如 Mersenne Twister 或 WELL 等),迭代变换可以表示为常量位矩阵乘以模 2 算术中的状态位向量 - 因此常数矩阵幂(在模 2 位算术中)就可以了。有了std::linear_congruential_engine它就更简单了 - 自己做数学作为一个简单的练习。基于椭圆曲线的PRNG也满足这些条件。(实际上我想知道为什么有人会在没有这个非常有用的属性的情况下设计 PRNG。 - 但那只是我。

我认为根本不存在这样的事情。 我的启发式结论是,O(1)-jump RNG 基本上是一个哈希值,以及这意味着的一切(例如,它可能根本不是"好"的 RNG)。 (请参阅@AkiSuihkonen的新答案和其中的链接。

但即使你问的是O(log z)我也没有看到在 STL 中实现。 在GCC的所有discard函数中,我能够grep都是简单的循环。

discard(unsigned long long __z)
{
for (; __z != 0ULL; --__z)
(*this)();
}

这不仅是可悲的,而且是误导性的,因为只有当有一种有效的方法来做到这一点时,discard才应该存在。

唯一不平凡的是mersenne(下图),但它仍然是O(z).

discard(unsigned long long __z)
{
while (__z > state_size - _M_p)
{
__z -= state_size - _M_p;
_M_gen_rand();
}
_M_p += __z;
}

Boost的Mersenne具有跳过功能,但仅对大于(默认值)10000000(!?)的跳过进行调用。这已经告诉我跳过在计算上非常繁重(即使它是O(log z))。 https://www.boost.org/doc/libs/1_72_0/boost/random/mersenne_twister.hpp

最后,推力显然对线性同余有一个有效的discard,但仅限于c == 0的情况下。(我不确定它是否会降低它作为 RNG 的用处。 https://thrust.github.io/doc/classthrust_1_1random_1_1linear__congruential__engine_aec05b19d2a85d02f1ff437791ea4dd68.html#aec05b19d2a85d02f1ff437791ea4dd68

所有基于计数器的随机数生成器

这些仅通过计算函数uint64_t rnd(uint64_t counter)uint16_t rnd(uint128_t counter)来工作。然后跳过功能就像

// the method itself is "randomly" generated -- known at least
// by von Neumann to typically generate poor results
struct MyRandomCBRNG {
uint64_t counter{0};
void skip(uint64_t a) { counter+=a;}
uint64_t operator()() {
uint64_t x = counter++;
// repeat multiple times if needed
x = x * 0xdeadbeefcafebabeull ^ (x >> 53) ^ (x << 11) + 0x13459876abcdfdecull;
return x;
}
};

人们甚至可以通过使用SHA-512之类的东西对由密钥连接的计数器进行哈希处理来制造加密上强大的CBRNG,更不用说Blum-Blum-Shub了。