【问题标题】:Simple random number generator that can generate nth number in series in O(1) time可以在 O(1) 时间内连续生成第 n 个数字的简单随机数生成器
【发布时间】:2014-02-22 04:43:42
【问题描述】:

我不打算将其用于安全目的或统计分析。我需要创建一个简单的随机数生成器,用于我的计算机图形应用程序。我不想使用“随机数生成器”这个词,因为人们对它的看法非常严格,但我想不出任何其他词来描述它。

  • 它必须很快。
  • 它必须是可重复的,给定一个特定的种子。 例如:如果种子 = x,那么每次我使用种子 x 时,序列 a,b,c,d,e,f..... 都应该发生。

最重要的是,我需要能够在恒定时间内计算系列中的第 n 项。

看来,我无法使用 rand_r 或 srand() 来实现这一点,因为这些需求取决于状态,我可能需要以某种未知顺序计算第 n 个。

我查看了线性反馈移位寄存器,但它们也取决于状态。

到目前为止,我有这个:

int rand = (n * prime1 + 种子) % prime2

n = 用于表示词条在序列中的索引。例如:对于 第一项,n ==1

prime1 和 prime2 是质数,其中 素数1 >素数2

seed = 一些允许使用相同功能的数字 根据种子产生不同的系列,但相同的系列 对于给定的种子。

我无法判断这是好是坏,因为我没有使用它,但如果有更多经验的人能指出这方面的问题,或者帮助我改进它,那就太好了。 .

编辑 - 我不在乎它是否可以预测。我只是想在我的计算机图形中创建一些随机性。

【问题讨论】:

  • 以前计算数字,并缓存它们,一旦得到第n个数字,只需O(1)。
  • 我可能无法缓存它们,因为可能有数百万个数字,我需要在移动设备上每秒至少进行 400,000 次计算。此外,缓存和查找过程可能需要比实际计算本身更长的时间......
  • 您当前的算法看起来一点也不随机。你有没有尝试过一个情节?我尝试使用 P1=569、P2=359、seed=12345 进行绘图,并且该模式真的可见。只是给你一个提示。
  • 我一直在使用 279470273UL 和 4294967291UL ,并且至少 10000 次计算我没有任何重复...你是如何绘制的?
  • 使用 Excel。有了您的数字,该模式甚至更多可见。在 Excel/OpenOffice Calc 中做一个 XY 绘图,你会看到它。 X轴是N,Y轴是数值。

标签: algorithm random


【解决方案1】:

CTR mode 中使用加密块密码。第 N 个输出只是加密(N)。这不仅为您提供了所需的属性(第 N 个输出的 O(1) 计算);它还具有很强的不可预测性。

【讨论】:

  • 会和我上面的方程式一样快,还是更快?我不在乎它是否可以预测。这是为了在计算机图形中创建一些“随机性”。
  • 问题不在于重复,而在于可见模式。例如,如果您在每行 N 行中的列 out(N) % img_width 处绘制一个点,并且 img_widthprime2 具有一定的数值​​关系,那么您将开始看到一个明显的模式......
  • @John 如果您真的采用这种格式,请尝试((n+seed)*(n+seed)*prime1 + (n+seed)*prime2) % prime3。分组密码并不复杂。想象一个巨大的码本,每个可能的 16 字节序列都有一个条目,其中每个码字也是 16 字节,随机选择(但没有两个词具有相同的码字)。您只需从书中选择第 n 个代码字即可获得您的号码。分组密码是一种生成码字的算法。
  • 或者,在更实际的层面上:const int KEY = 0x12345678; int numbersGenerated = 0; int generateRandomNumber() {return encrypt(numbersGenerated++, KEY);}。这并不完全准确(对于大多数分组密码来说,密钥和加密数据将是 8、16 或 32 个字节,而不是 4 个字节),但它显示了总体思路。密钥不需要保密,因为您已经说过不需要安全性。
  • @Edmund 在 1990 年代,我使用具有 减少轮数的 DES 作为实用但高质量的 PRNG,因此这个答案给出的建议是合理的。 Philox 是一个非常快的 PRNG,与这个答案提出的方法密切相关。
【解决方案2】:

不久前我偶然发现了这一点,正在寻找解决同一问题的方法。最近,我想出了如何在低常数 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 @、BBBBCCCC

【讨论】:

    【解决方案3】:

    正常意义上的RNG,有f(n) = S(f(n-1))这样的序列模式

    由于计算方便,它们在某些点也失去了精度(如 % mod),因此不可能将序列扩展为像 X(n) = f(n) = 仅具有 n 的平凡函数这样的函数。

    这意味着充其量你有 O(n) 的时间。


    因此,要针对O(1),你需要放弃f(n) = S(f(n-1))的想法,直接指定一个简单的公式,以便直接计算第N个数不知道第 (N-1) 个;这也使种子变得毫无意义。

    所以,你最终得到了一个简单的代数函数,而不是一个序列。例如:

    int my_rand(int n) { return 42; } // Don't laugh!
    int my_rand(int n) { 3*n*n + 2*n + 7; }
    

    如果您想对生成的模式(如分布)施加更多约束,这将成为一个复杂的数学问题。


    但是,对于您最初的目标,如果您想要的是恒定速度获取伪随机数,我建议使用传统 RNG 预先生成它并使用查找表访问。

    编辑:我注意到您对很多数字的表大小感到担忧,但是您可能会引入一些混合模型,例如包含 N 个条目的表,然后执行 f(k) = g( tbl[k%n] , k),这至少在 N 个连续序列中提供了良好的分布。

    【讨论】:

    • 我认为种子仍然可以发挥作用。例如:如果 7 是你的种子,并且你将它乘以 n 或它会的东西。与使用不同的种子相比,您会得到不同的序列。
    • seed 不同于函数的参数。种子直接分配用于产生下一个序列的“最后一个数字”。
    • 第一段的结论是错误的。举个反例,看看 Andre Harder 使用 O (log n) 中的线性同余生成器的答案。
    【解决方案4】:

    这演示了一个作为散列计数器实现的 PRNG。这似乎重复了 R. 的建议(在 CTR 模式下使用块密码作为流密码),但为此,我避免使用加密安全原语:为了提高执行速度并且因为安全性不是所需的功能。

    如果我们试图根据您的要求创建一个安全的流密码,即任何发出的序列都是可重复的,只要知道它的索引...

    ...然后我们可以选择一个安全的散列算法(如 SHA256)和一个有很多位的计数器(可能 2048 -> 序列每 2^2048 个生成的数字重复一次而不重新播种)。

    然而,我在这里展示的版本使用 Bob Jenkins 著名的散列函数(简单且快速,但不安全)以及 64 位计数器(与我的系统上的整数一样大,无需自定义递增)代码)。

    main 中的代码表明,在初始化后知道 RNG 的计数器(种子)允许重复 PRNG 序列,只要我们知道在重复点之前生成了多少值。

    实际上,如果您知道输出序列中任意点的计数器值,您将能够检索在该点之前生成的所有值,以及之后将生成的所有值。这仅涉及在与输出序列中的已知点关联的参考计数器值中添加或减去序数差异。

    修改这个类以用作测试框架应该很容易——您可以插入其他哈希函数并更改计数器的大小,以查看对速度以及生成值分布的影响(我所做的唯一一致性分析是在 main() 打印的十六进制数字的屏幕中寻找模式)。

    #include <iostream>
    #include <iomanip>
    #include <ctime>
    
    using namespace std;
    
    class CHashedCounterRng {
        static unsigned JenkinsHash(const void *input, unsigned len) {
            unsigned hash = 0;
            for(unsigned i=0; i<len; ++i) {
                hash += static_cast<const unsigned char*>(input)[i];
                hash += hash << 10;
                hash ^= hash >> 6;
            }
            hash += hash << 3;
            hash ^= hash >> 11;
            hash += hash << 15;
            return hash;
        }
    
        unsigned long long m_counter;
    
        void IncrementCounter() { ++m_counter; }
    
    public:
        unsigned long long GetSeed() const {
            return m_counter; 
        }
        void SetSeed(unsigned long long new_seed) {
            m_counter = new_seed; 
        }
        unsigned int operator ()() {
            // the next random number is generated here
            const auto r = JenkinsHash(&m_counter, sizeof(m_counter));
            IncrementCounter();
            return r;
        }
    
        // the default coontructor uses time() 
        // to seed the counter
        CHashedCounterRng() : m_counter(time(0)) {}
    
        // you can supply a predetermined seed here, 
        // or after construction with SetSeed(seed)
        CHashedCounterRng(unsigned long long seed) : m_counter(seed) {}
    };
    
    int main() {
        CHashedCounterRng rng;
        // time()'s high bits change very slowly, so look at low digits
        // if you want to verify that the seed is different between runs
        const auto stored_counter = rng.GetSeed();
        cout << "initial seed: " << stored_counter << endl;
        for(int i=0; i<20; ++i) {
            for(int j=0; j<8; ++j) {
                const unsigned x = rng();
                cout << setfill('0') << setw(8) << hex << x << ' ';
            }
            cout << endl;
        }
        cout << endl;
    
        cout << "The last line again:" << endl;
        rng.SetSeed(stored_counter + 19 * 8);
        for(int j=0; j<8; ++j) {
            const unsigned x = rng();
            cout << setfill('0') << setw(8) << hex << x  << ' ';
        }
        cout << endl << endl;
        return 0;
    }
    

    【讨论】:

    • 我注意到您正在检查重复输出作为随机性的度量;请记住,真正的随机数生成器将始终发出重复值,因为其输出流中的每一位都独立于其他每一位。事实上,经过一定数量的没有重复输出的观察后,越来越多的证据表明您正在查看的流不包含随机值。
    猜你喜欢
    • 1970-01-01
    • 2015-11-17
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2015-04-08
    • 2018-05-08
    • 2019-06-26
    相关资源
    最近更新 更多