不失一般性,生成[a,b]范围内的随机整数问题可以简化为在[0,s)范围内生成随机整数的问题。从均匀PRNG生成有界范围上的随机整数的现代技术由以下最近的出版物代表:
Daniel Lemire,“Fast Random Integer Generation in an Interval.” ACM Trans. Model. Comput. Simul. 29, 1, Article 3 (January 2019) (ArXiv draft)
Lemire表明他的算法提供无偏结果,并受到非常快速且高质量的PRNG(例如Melissa O'Neill的PCG发生器)越来越受欢迎的推动,展示了如何快速计算结果,几乎总是避免慢速分割操作。
下面展示了他算法的ISO-C实现,可以在
randint()
中看到。这里我结合了George Marsaglia的旧版
KISS64伪随机数生成器进行演示。出于性能原因,所需的64×64→128位无符号乘法通常最好通过特定于机器的内部函数或内联汇编来实现,直接映射到适当的硬件指令。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
typedef struct Prng_T *Prng_T;
uint64_t random64 (Prng_T);
void umul64wide (uint64_t, uint64_t, uint64_t *, uint64_t *);
uint64_t randint (Prng_T prng, uint64_t s)
{
uint64_t x, h, l, t;
x = random64 (prng);
umul64wide (x, s, &h, &l);
if (l < s) {
t = (0 - s) % s;
while (l < t) {
x = random64 (prng);
umul64wide (x, s, &h, &l);
}
}
return h;
}
#define X86_INLINE_ASM (0)
void umul64wide (uint64_t a, uint64_t b, uint64_t *hi, uint64_t *lo)
{
#if X86_INLINE_ASM
uint64_t l, h;
__asm__ (
"movq %2, %%rax;\n\t"
"mulq %3;\n\t"
"movq %%rax, %0;\n\t"
"movq %%rdx, %1;\n\t"
: "=r"(l), "=r"(h)
: "r"(a), "r"(b)
: "%rax", "%rdx");
*lo = l;
*hi = h;
#else
uint64_t a_lo = (uint64_t)(uint32_t)a;
uint64_t a_hi = a >> 32;
uint64_t b_lo = (uint64_t)(uint32_t)b;
uint64_t b_hi = b >> 32;
uint64_t p0 = a_lo * b_lo;
uint64_t p1 = a_lo * b_hi;
uint64_t p2 = a_hi * b_lo;
uint64_t p3 = a_hi * b_hi;
uint32_t cy = (uint32_t)(((p0 >> 32) + (uint32_t)p1 + (uint32_t)p2) >> 32);
*lo = p0 + (p1 << 32) + (p2 << 32);
*hi = p3 + (p1 >> 32) + (p2 >> 32) + cy;
#endif
}
struct Prng_T {
uint64_t x, c, y, z, t;
};
struct Prng_T kiss64 = {1234567890987654321ULL, 123456123456123456ULL,
362436362436362436ULL, 1066149217761810ULL, 0ULL};
#define MWC64 (kiss64->t = (kiss64->x << 58) + kiss64->c, \
kiss64->c = (kiss64->x >> 6), kiss64->x += kiss64->t, \
kiss64->c += (kiss64->x < kiss64->t), kiss64->x)
#define XSH64 (kiss64->y ^= (kiss64->y << 13), kiss64->y ^= (kiss64->y >> 17), \
kiss64->y ^= (kiss64->y << 43))
#define CNG64 (kiss64->z = 6906969069ULL * kiss64->z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)
uint64_t random64 (Prng_T kiss64)
{
return KISS64;
}
int main (void)
{
int i;
Prng_T state = &kiss64;
for (i = 0; i < 1000; i++) {
printf ("%llu\n", randint (state, 10));
}
return EXIT_SUCCESS;
}