c - 通过 6 面模具实现提高 7 面模具滚动模拟的性能

Question

在不同的StackExchange , 以下算法(基于算术编码)被提出作为生成 7 面骰子结果的有效方法，当所有给定的是 6 面骰子时:

int rand7()
{
  static double a=0, width=7;  // persistent state

  while ((int)(a+width) != (int)a)
  {
    width /= 6;
    a += (rand6()-1)*width;
  }

  int n = (int)a;
  a -= n; 
  a *= 7; width *= 7;
  return (n+1);
}

我不是真正的数学家，我会尽我所能解释这个算法是如何工作的:

在每次调用 rand7() 开始时，width 是 7^s/6^t，a 是一个非负值，a + width 位于基本情况之后的区间 [0, 7) 内。当进入while 循环时，width 是可以添加到a 的最大值。如果 floor(a + width) 与 floor(a) 不同，则随机选择 { 0, width*1/6, 宽度*1/3, 宽度*1/2, 宽度*2/3, 宽度*5/6 } 被添加到 a，指数 t 增加 1(将 width 的值减少 6 的幂)。请注意，在迭代之后，a + width 位于区间 [0, 7) 的属性保持不变。当 width 变得小于 ceil(a) - a 的差时，迭代停止。循环将更多的熵添加到 a 中，只要这样做实际上可以影响掷骰子的结果，并且直观地，这是在 [0, 7 范围内构建随机实数) 使用基数 6。离开循环后，掷骰子为 floor(a) + 1，a 减为小数部分。此时 a + width 位于区间 [0, 1) 内。为了准备下一次调用并保持不变属性，a 和 width 都被放大了 7 倍(对于 width，这会使指数 s 增加 1)。

上面解释了归纳步骤的工作原理。基本案例的分析留给感兴趣的读者作为练习。

当然，从效率的角度来看，使用浮点运算会立即弹出性能拖累(假设 rand6() 的性能已经足够并且本身无法提高)。在保持这种算术编码算法的同时，去除 float 使用的最佳方法是什么？

Answer 1

根据我发表的评论，这是该算法的定点版本。它使用无符号的 4.60(即数字的小数部分有 60 位)，这比你从 double 中得到的要多几位:

int rand7fixed() {
    static uint64_t a = 0;
    static uint64_t width = 7UL<<60;
    static const uint64_t intmask = 0xfUL<<60;

    while (((a+width)&intmask) != (a&intmask)) {
      width /= 6;
      a += (rand6()-1)*width;
    }

    int n = a >> 60;
    a &=~intmask;
    a *= 7;
    width *= 7;
    return n+1;
}

上面的结果比 OP 中的 double 版本快了大约三分之一(参见下面的基准测试结果和注释)。我怀疑浪费时间的不是浮点运算，而是从 double 到 int 的转换。

正如 R.. 指出的那样，这并不能解决偏见；它只是减少它。一种简单且相当快速的无偏算法是重复生成两个 rand6() 值 h 和 l 直到其中至少一个不为零，然后返回 h*6+l % 7 :

int rand7_2() {
  int hi, lo;
  do {
    hi = rand6() - 1;
    lo = rand6() - 1;
  } while (hi + lo);
  return (hi * 6 + lo) % 7 + 1;
}

如果您觉得需要减少对 rand6 的调用次数，您可以使用 6¹² 仅略大于 7¹¹ 的事实来生成 11 7 -一次掷骰子。为了消除偏差，仍然需要丢弃一些 12 个 6 卷的集合；丢弃集的频率将是 (6<sup>12</sup>−7<sup>11</sup>)/6<sup>12</sup>) ，或大约 11 分之一，因此平均每 7 卷需要大约 1.19 6 卷。您可以通过使用 25 个 6 卷生成 23 个 7 卷(每 7 卷 1.13 个 6 卷)来做得更好，但这不太适合 64 位算法，因此调用 rand6 的边际优势将被削弱以 128 位计算。

这是 11/12 的解决方案:

int rand7_12() {
    static int avail = 0;
    static uint32_t state = 0;
    static const uint32_t discard = 7*7*7*7*7*7*7*7*7*7*7; // 7 ** 11
    static const int out_per_round = 11;
    static const int in_per_round = 12;

    if (!avail) {
      do {
        state = rand6() - 1;
        for (int needed = in_per_round - 1; needed; --needed)
          state = state * 6 + rand6() - 1;
      }
      while (state >= discard);
      avail = out_per_round;
    }
    int rv = state % 7;
    state /= 7;
    --avail;
    return rv + 1;
}

理论上，您应该能够将比率降低到 log<sub>7</sub>6，大约为 1.086。例如，您可以通过从 972 个 6 卷生成 895 个 7 卷，丢弃大约 1600 个集合中的一个来实现平均 1.087 个 6 卷/7 卷，但是您需要 2513 位算法来保持状态。

我使用不太精确的基准测试了所有四个函数，该基准调用 rand7 700,000,000 次，然后打印结果的直方图。结果:

                                                  User time with
Algorithm       User time      rand6() calls     cycling rand6()
----------      ---------      -------------     ---------------
double          32.6 secs          760223193           13.2 secs
fixed           29.4 secs          760223194            7.9 secs
2 for 1         40.2 secs         1440004276
12 for 11       23.7 secs          840670008

上面的底层 rand6() 实现是 Gnu 标准 c++ 库的 uniform_int_distribution<>(1,6)，使用 mt19937_64(64 位 Mersenne Twister)作为 PRNG。为了更好地处理在标准库中花费的时间，我还使用一个简单的循环计数器作为伪随机数生成器来运行测试；剩下的 13.2 秒和 7.9 秒代表(粗略地)算法本身花费的时间，从中我们可以说定点算法快了大约 40%。 (很难很好地阅读组合算法，因为固定序列使分支预测更容易并减少了 rand6 调用的次数，但两者都花费了不到 5 秒。)

最后，直方图，以防有人想检查偏差(还包括使用 std::uniform_int_distribution(1,7) 运行以供引用):

Algorithm          1          2          3          4          5          6          7
---------  ---------  ---------  ---------  ---------  ---------  ---------  ---------
reference  100007522  100002456  100015800  100005923   99972185  100008908   99987206
double     100014597  100005975   99982219   99986299  100004561  100011049   99995300
fixed      100009603  100009639  100034790   99989935   99995502   99981886   99978645
2 for 1    100004476   99997766   99999521  100001382   99992802  100003868  100000185
12 for 11   99988156  100004974  100020070  100001912   99997472   99995015   99992401