不久前我偶然发现了这一点,正在寻找解决同一问题的方法。最近,我想出了如何在低常数 O(log(n)) 时间内做到这一点。虽然这与作者要求的 O(1) 不太匹配,但它可能足够快(使用 -O3 编译的示例运行,实现了 10 亿个任意索引随机数的性能,n 在 1 和 2 之间变化^ 48,在 55.7 秒内 - 略低于 1800 万个数字/秒)。
一、解决方案背后的理论:
一种常见的 RNG 类型是Linear Congruential Generators,基本上,它们的工作方式如下:
随机(n) = (m*随机(n-1) + b) mod p
其中 m 和 b 以及 p 是常数(有关如何选择它们的信息,请参阅 LCG 的参考资料)。由此,我们可以使用一些模运算来设计以下内容:
random(0) = seed mod p
random(1) = m*seed + b mod p
random(2) = m^2*seed + m*b + b mod p
...
random(n) = m^n*seed + b*Sum_{i = 0 to n - 1} m^i mod p
= m^n*seed + b*(m^n - 1)/(m - 1) mod p
计算上述内容可能会出现问题,因为这些数字很快就会超过数字限制。一般情况的解决方案是用 p*(m - 1) 以模计算 m^n,但是,如果我们取 b = 0(LCG 的子情况,有时称为Multiplicative congruential Generators),我们有一个更简单的解,并且只能在模 p 中进行计算。
在下文中,我使用了 RANF(由 CRAY 开发)使用的常量参数,其中 p = 2^48 和 g = 44485709377909。p 是 2 的幂的事实减少了所需的操作数量(如预期的那样) ):
#include <cassert>
#include <stdint.h>
#include <cstdlib>
class RANF{
// MCG constants and state data
static const uint64_t m = 44485709377909ULL;
static const uint64_t n = 0x0000010000000000ULL; // 2^48
static const uint64_t randMax = n - 1;
const uint64_t seed;
uint64_t state;
public:
// Constructors, which define the seed
RANF(uint64_t seed) : seed(seed), state(seed) {
assert(seed > 0 && "A seed of 0 breaks the LCG!");
}
// Gets the next random number in the sequence
inline uint64_t getNext(){
state *= m;
return state & randMax;
}
// Sets the MCG to a specific index
inline void setPosition(size_t index){
state = seed;
uint64_t mPower = m;
for (uint64_t b = 1; index; b <<= 1){
if (index & b){
state *= mPower;
index ^= b;
}
mPower *= mPower;
}
}
};
#include <cstdio>
void example(){
RANF R(1);
// Gets the number through random-access -- O(log(n))
R.setPosition(12345); // Goes to the nth random number
printf("fast nth number = %lu\n", R.getNext());
// Gets the number through standard, sequential access -- O(n)
R.setPosition(0);
for(size_t i = 0; i < 12345; i++) R.getNext();
printf("slow nth number = %lu\n", R.getNext());
}
虽然我认为作者现在已经继续前进,但希望这对其他人有用。
如果您真的关心运行时性能,使用查找表可以使上述速度提高大约 10 倍,但代价是编译时间和二进制大小(它也是 O(1) w.r.t OP 要求的所需随机索引)
在下面的版本中,我在编译时使用 c++14 constexpr 生成查找表,每秒达到 176M 任意索引随机数(这样做确实增加了大约 12 秒的额外编译时间,并且二进制大小增加了 1.5MB——如果使用部分重新编译,增加的时间可能会减少)。
class RANF{
// MCG constants and state data
static const uint64_t m = 44485709377909ULL;
static const uint64_t n = 0x0000010000000000ULL; // 2^48
static const uint64_t randMax = n - 1;
const uint64_t seed;
uint64_t state;
// Lookup table
struct lookup_t{
uint64_t v[3][65536];
constexpr lookup_t() : v() {
uint64_t mi = RANF::m;
for (size_t i = 0; i < 3; i++){
v[i][0] = 1;
uint64_t val = mi;
for (uint16_t j = 0x0001; j; j++){
v[i][j] = val;
val *= mi;
}
mi = val;
}
}
};
friend struct lookup_t;
public:
// Constructors, which define the seed
RANF(uint64_t seed) : seed(seed), state(seed) {
assert(seed > 0 && "A seed of 0 breaks the LCG!");
}
// Gets the next random number in the sequence
inline uint64_t getNext(){
state *= m;
return state & randMax;
}
// Sets the MCG to a specific index
// Note: idx.u16 indices need to be adapted for big-endian machines!
inline void setPosition(size_t index){
static constexpr auto lookup = lookup_t();
union { uint16_t u16[4]; uint64_t u64; } idx;
idx.u64 = index;
state = seed * lookup.v[0][idx.u16[0]] * lookup.v[1][idx.u16[1]] * lookup.v[2][idx.u16[2]];
}
};
基本上,它的作用是将m^0xAAAABBBBCCCC mod p 等的计算拆分为(m^0xAAAA00000000 mod p)*(m^0xBBBB0000 mod p)*(m^CCCC mod p) mod p,然后为0x0000 - 0xFFFF 范围内的每个值预先计算表,这些值可以填充@987654331 @、BBBB 或 CCCC。