CSDN专帖系列之一： 根据某一特殊规律的概率生成随机数

CSDN专帖系列之一： 根据某一特殊规律的概率生成随机数
Author：Jeff 2006-7-7
关键字：C++ 算法 随机数 概率
环境：Window XP Professional + SP2, VC6.0

CSDN原帖：
主　　题：求一算法，随机数和概率相关
作　　者：adintr (www.adintr.com)
链 接：http://community.csdn.net/Expert/topic/4753/4753180.xml?temp=.4426691
描 述：
有一整数 base ，其值大概在 10 - 100 左右， 现在需要生成一随机整数 x，其范围在 0 - n ( 0<n<17) 之间，要求这个随机数出现的概率为 ：
base ^ x / (base ^0 + base ^ 1 + base ^ 2 + ... + base ^ n)
也就是说如果假设生成的随机数为 0 的概率为 1 点的话，则生成的随机数为 1 的概率就是 base 点，为 2 的概率就是 base * base 点，为 3 的概率就是 base ^ 3 点。

注意 base 和 n 的取值达到较大时 100^16 是超过整数的表示范围的。
效率当然是越高越好， base 和 n 的值可以预先确定下来

在帖子中楼主已经对他的问题作了解答，他的回答算法是这样的：
首先生成 n 个 [0 ，base + 1) 之间的随机数，则每一位上出现 0 的几率是 1/(base + 1), 为非 0 的几率是 base / (base + 1)。
现在这 n 个随机数全为 0， 则让最后的结果 x = 0, 这种情况出现的几率是 1 / (base + 1)^n
如果第一位为非 0， 其余全为0， 则 x = 1, 出现的几率为 base / (base + 1) ^ n
如果前两位为非 0， 其余为 0， 则 x = 2, 出现几率为 base ^ 2 / (base + 1) ^ n
同理，当前 x 位为非 0， 其余为 0，则 x = x, 出现几率就是 base ^ x / (base + 1)^n
这样就保证了每以个数出现的几率都比小一的那个数大 base 点。

关于判断位为不为0时，是不是要求末尾的0连续？当时没有想清楚，后来明白了，这是一种排列情况。虽然有很多种组合方式，但是一旦判断条件定下来，情况也就确定了，每次都必须符合这种特定的排列，是作为组合的一种。因此对于末尾的0并不需要连续。可以这样理解：
前两位为非 0， 其余为 0， 则 x = 2, 出现几率为 base ^ 2 / (base + 1) ^ n
如果第一位为非0，末尾为非0，其余位为0，出现几率为
[base/(base + 1)] * [1/(base + 1)] * ··· * [1/(base + 1)] * [base/(base + 1)]
= base ^ 2 / (base + 1) ^ n
两个几率是一样的。在程序中只处理一种情况就够了，且只能处理一种情况。

另外，当n值越大，程序就会递归得越深(因为连续随机到n个0的情况很少)。因此层次也不确
定，可能很深， 有可能栈崩掉。

后来，经过分析，对我的算法进行修改，并加以实现。呵呵，发现速度不错。*_*
经过修改后的算法：
1． 先计算(0, RAND_MAX]内base的幂M，使得M满足K = base ^ M + base ^ (M - 1) + … + base <= RAND_MAX。(因为随机数最大为RAND_MAX)
2． 计算各段区间的上限和下限。
3． 如果n > M，直接到4；否则，跳到6。
4． 生成[0，K]范围内的一个随机数S。
5． 如果S ≠ 0, 计算S位于哪块区间，将其记为H, 则生成的随机数为H + n – M – 1。否则
n ß n – M, 跳到3。
6． 计算K2 = base ^ n + base ^ (n - 1) + … + base。
7． 生成[1，K2]范围内的一个随机数S2。
8． 计算S2位于哪块下限区间内，记为H2, 则生成的随机数为H2 – 1。

将base的幂作为分割点把整个线形空间(即[0，base ^ n + base ^ (n - 1) + … + base]空间)分成n + 1段，那么区间的分割状况就是：
[0，0]，
[1，base]，
[base + 1，base ^ 2 + base]，
……，
[base ^ (n - 1) + base ^ (n - 2) + … + base + 1，base ^ n + base ^ (n - 1) + … + base]
则各段区间对应的比例为：
1 / (base ^ n + base ^ (n - 1) + … + base + 1)
base / (base ^ n + base ^ (n - 1) + … + base + 1)
…….
base ^ n / (base ^ n + base ^ (n - 1) + … + base + 1)
和产生随机数的概率是一一对应的。因此，只需要随机数能够分布到整个线形空间，落在哪个区间，则生成了对应的随机数。
但是base ^ n可能很大，不只rand()函数产生不了这么大范围的值，连int32都可能表示不了。

这个时候，应该对区间进行适当的缩放，满足rand()函数能够产生收缩后的区间内的随机数，当然这个区间应该尽可能得大。计算一下区间的缩小比例，



L = base ^ M + base ^ (M – 1) + … + base + 1应该尽可能的接近RAND_MAX，且比RAND_MAX小或者相等。在允许误差的情况下，整个线形空间[0，base ^ M + … + base + 1 )各个区间对应的产生的随机数为：
[0，1)： 0 < x < n – M
[1, base + 1)： n – M + 1
……
[base ^ (M – 1) + … + base + 1, base ^ M + base ^ (M – 1) + … + base + 1)： n
当随机数为0时，递归处理余下的线形空间，直到没有余下空间为止。用图例说明如下：



下面的是代码及注释：

int GetPower(long base) // 也就是求上面所说的M。
{
int pow = -1;
long sum = 1L, old = sum;
long total = 0L;
// 这里没有计算最后的+1，因为那是一个开区间。
while (32767 > total && old <= sum) {
old = sum;
sum *= base;
total += sum;
++pow;
}

return pow; //返回base的幂。
}

long power(long x, long y) // 计算x ^ y。(几个函数名都意义重了，不管了~~~)
{
long sum = 1L;
while (y-- > 0L) {
sum *= x;
}
return sum;
}

int main(void) {
srand((unsigned)time(NULL));

// 下面的代码对下面宏的边界情况没有处理
#define NUM_BASE 10 // 即base（不要大于RAND_MAX）
#define NUM_POWER 6 // 即n（不要小于0）

unsigned long count[NUM_POWER];
memset(count, 0, sizeof(count));

int pow = GetPower(NUM_BASE); // 获得幂M
int smaller = pow > NUM_POWER ? NUM_POWER : pow;

long region[NUM_POWER + 1];
long lower_limit[NUM_POWER + 1];
memset(region, 0, sizeof(region)); //区间上限
memset(lower_limit, 0, sizeof(lower_limit)); // 区间下限
for (int i = 1; i <= smaller; ++i) {
region[i] = region[i - 1] + power(NUM_BASE, i); // 根据base计算上限
lower_limit[i] = region[i - 1] + 1L; // 下限
}

long k = 100000000L;
while (--k >= 0L) {
int n = NUM_POWER;
while (0 < n - pow) {
long seed = rand() % (region[smaller] + 1); //生成随机数[0， region[smaller] ]
int s = smaller + 1;
while (--s >= 0 && lower_limit[s] > seed); //查找.
if (s > 0) { //找到了，（下限不是0）
count[n - pow + s - 1]++;
break;
}
n -= pow; //seed为0的情况，循环处理剩下的空间。
}
if (0 >= n - pow) { //处理余下的一小段空间。
long seed = (long)rand() % region
+ 1; // [1， region
]
int s = n + 1;
while (--s >= 0 && lower_limit[s] > seed);
if (s > 0) {
count[s - 1]++;
}
}
}

for (k = 0L; k < NUM_POWER; ++k) {
printf("%d: %d/n", k, count[k]);
}

return 0;
}

第一次的测试结果：（5秒）
0: 847
1: 8344
2: 91541
3: 915915
4: 9153458
5: 89829895

第二次的测试结果：（5秒）
0: 843
1: 8316
2: 91336
3: 915133
4: 9155684
5: 89828688

可以看出规律大概是对的，只是0，1和其他数值相比，比例降低得更快。因为0和1合起来才占据了区间的一个单位，现在把0和1对应的值加起来差不多900，比较符合。
理论上来说，楼主的算法接近真实值一些，但是速度太慢。而这个算法相对来说，误差大，速度快。时间基本上都花在了最外围的循环上。