arc4random_uniform 和 PCG 中的均匀分布答案

【问题标题】：Uniform distribution in arc4random_uniform and PCGarc4random_uniform 和 PCG 中的均匀分布
【发布时间】：2021-11-06 03:35:06
【问题描述】：

来自OpenBSD 的arc4random_uniform 和Melissa O'Neill 的PCG 库都具有相似的算法来生成一个无偏差的无符号整数值，最高不包括上限。

inline uint64_t
pcg_setseq_64_rxs_m_xs_64_boundedrand_r(struct pcg_state_setseq_64 * rng,
                                        uint64_t bound) {
    uint64_t threshold = -bound % bound;
    for (;;) {
        uint64_t r = pcg_setseq_64_rxs_m_xs_64_random_r(rng);
        if (r >= threshold)
            return r % bound;
    }
}

-bound % bound 不总是零吗？如果它总是为零，那为什么还要有循环和 if 语句呢？

OpenBSD 也有同样的东西。

uint32_t
arc4random_uniform(uint32_t upper_bound)
{
    uint32_t r, min;

    if (upper_bound < 2)
        return 0;

    /* 2**32 % x == (2**32 - x) % x */
    min = -upper_bound % upper_bound;

    /*
     * This could theoretically loop forever but each retry has
     * p > 0.5 (worst case, usually far better) of selecting a
     * number inside the range we need, so it should rarely need
     * to re-roll.
     */
    for (;;) {
        r = arc4random();
        if (r >= min)
            break;
    }

    return r % upper_bound;
}

Apple 版本的arc4random_uniform 有不同的版本。

u_int32_t
arc4random_uniform(u_int32_t upper_bound)
{
    u_int32_t r, min;

    if (upper_bound < 2)
        return (0);

#if (ULONG_MAX > 0xffffffffUL)
    min = 0x100000000UL % upper_bound;
#else
    /* Calculate (2**32 % upper_bound) avoiding 64-bit math */
    if (upper_bound > 0x80000000)
        min = 1 + ~upper_bound;     /* 2**32 - upper_bound */
    else {
        /* (2**32 - (x * 2)) % x == 2**32 % x when x <= 2**31 */
        min = ((0xffffffff - (upper_bound * 2)) + 1) % upper_bound;
    }
#endif

    /*
     * This could theoretically loop forever but each retry has
     * p > 0.5 (worst case, usually far better) of selecting a
     * number inside the range we need, so it should rarely need
     * to re-roll.
     */
    for (;;) {
        r = arc4random();
        if (r >= min)
            break;
    }

    return (r % upper_bound);
}

【问题讨论】：

-bound % bound 在数学上不是 -bound % bound。因为bound 是uint64_t，所以它是(2^64 - bound) % bound。它正在计算 2^64 模数的余数，以从 2^64 间隔中剥离该数量，只留下边界的倍数。
@EricPostpischil，所以给定值 3，C 中的正确表达式是 (uint64_t)-3 % 3?
表达式-bound % bound是正确的，因为bound的类型是uint64_t，所以-是模2^64执行的。使用-3，- 使用int 算术完成。因此，需要对 uint64_t 进行强制转换以将其模 2 ^ 64 包装。 -bound % bound 不需要这个演员表。（除非在 C 实现中 int 比 64 位更宽，因为那时 bound 将被提升为 int。但是，我不相信在普通使用中存在任何这样的 C 实现。）

标签： c random uniform-distribution

【解决方案1】：

因为bound 是uint64_t，所以-bound 是模2⁶⁴。结果是 2⁶⁴−bound，而不是 −bound。

然后-bound % bound计算2⁶⁴−bound模bound的余数。这等于 2⁶⁴ 模 bound 的余数。

通过将 threshold 设置为此并拒绝小于 threshold 的数字，例程将接受的间隔减少到 2⁶⁴−threshold 数字。结果是一个包含多个数字的区间，这些数字是bound 的倍数。

从在该间隔中选择的数字r，例程返回r % bound。由于区间的修整，每个残基出现的次数相等，因此结果对任何残基都没有偏差。

【讨论】：