【编程珠玑】【第一章】生成随机数、随机取样的问题
一、利用随机数函数生成随机数
问题1（《编程珠玑》习题12.1后半段）：
给定一个rand()，可以产生从0到RAND_MAX的随机数，其中RAND_MAX很大（常见值:16位int能表示的最大整数32767），写出利用rand()生成[a,b]中任意整数的函数，其中a>=0, b<=RAND_MAX，且b-a<<RAND_MAX。
分析：
这是在编程工作最常见的随机函数的应用，在这里做一个起点再合适不过。把随机数区间的起点从0变为a，同时把一共RAND_MAX+1个数的区间缩小至只含有b-a+1个数的区间，写为 a + rand()%(b-a+1)，此时显然最大值是a+(b-a)= b。进一步地，这个b-a<<RAND_MAX的条件虽然看上去不起眼，其实很重要。
附加思考：
如果b-a和RAND_MAX很接近会发生什么情况？读者不妨先做思考，问题2的分析会做出解答。这个rand()，其实相当于《编程珠玑》提到的bigrand()。
问题2（笔试题变形）：
给定一个随机数函数rand7()，它能以等概率生成1~7这7个整数。请根据rand7()写出类似的rand5()。
分析：
如果直接像问题1中一样，把1+rand7()%5作为rand5()会有什么情况发生？这时确实能产生1~5的随机数没错，可是各个数的概率相等吗？
对于随机数2，既有可能来自于1+1%5，也有可能来自于1+6%5，显然其概率是2/7，而不是1/5，不满足rand5()等概率产生各随机数这一隐含要求。不同于问题1，问题1中一个很大的区间收缩成较小区间时，各个元素映射后的新元素概率虽然概率可能不完全一样，但却是近似相同的。
为了满足等概率的要求，可以这么做：
int rand5() { int res; do { res = rand7(); } while(t>=6); return res; }
虽然保证了1~5的概率都变成了1/5，但是有一个无法避免的缺点是，每当产生了6或者7都要抛弃，相当于这一次运行是“空转”，浪费了时间。如果对1/5这个概率不明白，可以有两种理解：每次产生6或7就被抛弃，剩下数的概率相等，必然为1/5；或者用更严密的推理：产生1~5的随机数，最终得到某一个的概率为：1/7+(1/7)*(2/7)+(1/7)*(2/7)2+...，无限项求和，结果是1/5。
问题3（笔试题原题）
给定一个随机数函数rand7()，它能以等概率生成1~7这7个整数。请根据rand7()写出类似的rand10()。
分析：
有了问题2的概率基础，把rand7()变成rand10()仅仅需要一点点思考了。
int rand10() { int t1,t2,res; do { t1 = rand7(); //t1以等概率1/5成为{1,2,3,4,5}中的一个数字 } while(t1>=6); do { t2 = rand7(); //t2以1/2概率成为奇数或偶数，因为{1,2,3,4,5,6]中3奇3偶 } while(t2==7); res = t1+5*(t2%2); //res是1~10中的任意一个数的概率都是1/10，注意到%和*具有相同的优先级，这里去掉括号结果是错的 return res; }
解释：
t1的概率是{1~5}中每个元素概率为1/5，t2是1/2概率为奇或者偶，所以(t2%2)是以1/2概率选择0或者1，因此5*(t2%2)是以1/2概率为0或者5。这样，选择1的概率是1/5 * 1/2=1/10，...，选择6的概率是1/5 * 1/2为1/10....由此可见。
问题2和问题3是对《编程珠玑》上使用范围不大的randint(i,j)生成其他范围随机数方法的解答。在掌握了问题2和问题3的解法后，你已经学会随机数区间的收缩和扩张，类似的问题迎刃而解。
问题4（《编程珠玑》习题12.1前半段）
C库函数rand()常返回15个随机位，写出bigrand()，能够返回30个随机位。
分析：
其实问题4和问题3有点像，但是不同之处在于，这次我们的视角是从位出发的，把rand()看做了将15个位每一个位以1/2概率设为0或1，从而生成0~RAND_MAX。从某种意义上来说，按这种理解方式来解这个问题更轻松一些，但是仅局限于2的幂减1这样的数值的区间，比如从0到11...11。在此基础之上把两部分的位拼接起来即可。
//《编程珠玑》的答案 //int bigrand() {return RAND_MAX *rand() +rand();} //最大值不是30个1，在另一个笔记中有叙述，怀疑有错我的答案，把先生成的部分左移15位作为高15位 long bigrand(){return ((long)RAND_MAX+1)*rand() +rand();}
扩展：拒绝采样
这是一个简单粗暴有效的方法，使用这个方法你可以不用考虑复杂的区间扩展和收缩的问题。以rand7()产生rand10()为例，构造一个二维数组并把它填成：
然后使用两次rand7()分别生成行号和列号，如果对应元素是0，抛弃重来；否则即是结果。其实这种方法依然用到了区间扩展：把7扩展成7*7，并把不符合的部分抛弃。这样横向生成一个1~7的随机数，纵坐标生成1~7的随机数，这样两次随机数便能定位一个数组元素，若该元素为1~10之间的数，则被选中，否则重新选择。这能保证1~10之间的每个数被选中的概率都是1/10。从这里可以看出这种方法其实还是有缺陷的：如果用rand7()生成rand50()，那这50个数可是填满这个表格后还有填不下的。
二、利用随机数函数产生随机事件
用随机事件表示随机事件？这个问题具体化为：使用rand()表示以M/N的概率发生的随机事件，其中M<=N，并可用作：if(事件A发生) ，其中P(A) = M/N，那么表示为：
if(rand()%N< M)
...
通过这种方式，我们就可以做出让程序“以M/N的概率执行某个命令”这样的设计了。
三、取样问题:从n个元素中随机选取m个
1.从概率角度出发
考虑整数0,1,2,...,n-1，可以用上节的方法以概率m/n选取0（推导方式略）。但是对于1，必须考虑之前0是否被选取而以(m-1)/(n-1)或m/(n-1)的概率选取1，后续就更加麻烦。好在迭代是计算机的长项，只需要把这个是否选择当前数的随机事件稍作修改即可，使之变成从r个其余的数中选择s个：
int gen(int m,int n){ int i, select = m,remaining = n; for(i=0;i<n;i++) { if(rand() % remaining <select) { printf("%d\n",i); select--; } remaining--; } return 0; }
其概率证明可见于Knuth的The Art of Computer Programming第2卷。进一步地可以优化为
int genknuth(int m,int n){ int i; for(i=0;i<n;i++) if(rand()%(n -i) < m) { printf("%d\n",i); m--; } return 0; }
《编程珠玑》提示，这种算法是“所有n的所有m元子集被选中的概率相等”，这个条件强于“所有元素被选中的概率相同”。下面是习题12.2中提到的“以等概率选择所有元素，但是有些m元子集被选中的概率比其他子集大”的算法：直接选择1个数，则这个m元集合为它本身即后续的一共m个数，可能包括回绕。
对于这种方法，总要产生n次随机数。进一步可以写为for(i=0;i<n&&m>0;i++)，但程序的平均运行时间是否变得更快需要权衡。对应地，习题12.7提供了一种稍微减少随机数产生次数的递归函数：
int randselect(int m,int n) { int r; //assert(m<=n && m>=0); if(m>0) { r = rand()%n; if(r < m) { printf("%d\n",n-1); randselect(m-1,n-1); } else randselect(m,n-1); } return 0; }
2.从集合插入出发
由于集合元素不重复，如果按等概率选择一个随机数，不在集合中就把它插入，反之直接抛弃，直到集合元素个数达到m个，同样可以满足要求，并且用C++的STL很容易实现：
void gensets(int m,int n) { set<int> S; while(S.size() < m) S.insert(rand()%n); set<int>::iterator i; for(i = S.begin();i!=S.end();++i) cout<<*i<<"\n"; }
这个算法的主要问题是，如果抛弃已存在的元素的次数过多，相当于多次产生随机数并进行集合操作，性能将明显下降。比如当n=100而m=99，取第99个元素时，算法“闭着眼睛乱猜整数，直到偶然碰上正确的那个为止”（《编程珠玑（续）》，13.1节）。虽然这种情况会在“从一般到特殊”提供解决方案，但下面的Floyd算法明显规避了产生随机数超过m次的问题。
习题12.9提供了一种基于STL集合的随机数取样方法，可以在最坏情况下也只产生m个随机数：限定当前从中取值的区间的大小，每当产生重复的随机数，就把这一次迭代时不会产生的第一个随机数拿来替换。
int genfloyd(int m,int n){ set<int> S; set<int>::iterator i; for(int j = n-m; j<n;j++) { int t = rand()%(j+1); if(S.find(t) == S.end()) S.insert(t); else S.insert(j); } for(i=S.begin();i!=S.end();++i) cout<<*i<<"\n"; }
不必基于集合而实现，我自己写了一个原理类似的算法，思想是把“用过”的元素“扔”到下一次迭代的随机数取样区间之外：
int gen(int m,int n) { int *array; int i,j; array = (int *)malloc(sizeof(int) * n); for(i=0;i<n;i++) array[i] = i; while(m>=1) { j = rand()%n; printf("%d\n",array[j]); if(j<n-1) array[j] = array[n-1]; m--; n--; } return 0; }
3、从“打乱顺序”出发
这是个来源于实际的想法：将所有n个元素打乱，取出前m个。更快的做法是，打乱前m个即可。对应的C++代码如下：
int genshuf(int m,int n){ int i,j; int *x = new int[n]; for(i = 0;i<n;i++) x[i] = i; for(i = 0;i<m;i++) { j = randint(i,n-1); //randint产生i到n-1之间的随机数 int t = x[i];x[i] = x[j];x[j] = t; } //sort(x,x+m); //sort是为了按序输出 for(i=0;i<m;i++) cout<<x[i]<<"\n"; }
把sort注释掉的这段代码，可以作为随机不重复序列产生器。类似的还有Floyd的算法P。（《编程珠玑（续）》，13.3节）
4、从一般到特殊
以上讨论的几种方式都不限定m和n的取值，只需m<=n即可。对于特殊的取值，有特殊的解决方案，以下是编程珠玑上的两例：
1.n = 106而m=n-10，这时m很大较为接近n，可以先生成10个元素，然后输出其余的元素。进行这种处理的额外代码可以提高算法的平均速度。
2.n=231而m = 107，m<<n，这时可以先生成1.1*107个元素，排序后去掉重复的（由于n很大，m中出现重复元素的概率很低），得到 107个元素的样本。
附：
（《编程珠玑（续）》习题13.5）Doug McIlory的求N个元素中取M种的第G种情况的组合的算法，原书这个算法我没有理解，也没有看到比较满意的解释，可能用到了某些我不熟悉的组合数性质。原书中没有对其进行进一步解释，并且似乎只是作者的题外话而已。如果看着有困难，这部分代码完全可以跳过。介绍这个算法的原因是用它能在获得随机数G后，直接获得第G种N个元素取M个的取法，相当于只产生了一次随机数。
int comb(int n,int m,int g,int* array){ int d = 1,t; while(m>0) { t = combination(n-d,m-1); //combination(n,m)是计算n中取m个的组合数 if(g-t<0) { m--; array[m] = d; printf("%d\n",d); } else g-=t; d++; } return 0; }
三、取样问题:从未知总数的元素中选择一个
从事先未知总数为n的对象中随机选择一个的方法。有两种常见的具体问题：
1.读取一个未知行数的文件，随机输出其中的的一行，同时最多只能缓冲一行的内容（《编程珠玑（续）》习题15.3利用了这种形式）；
2.对链表进行一趟遍历，随机输出其中的一个结点的元素，只能使用一个临时指针。
解法是，以概率1选择第一个元素存入缓存，以概率1/2用第二个元素替换掉缓存，...，直至遍历完所有元素，最后输出缓存的内容。可以分析出此时所有元素留在缓存的概率均为1/n，比如1，是1*(1/2)*(2/3)*...*((n-1)/n)。
int random_select(void){ int i,num=1; for(i=1;i<n;i++) if(rand()%i ==0) //等于0或者任何一个0~i-1之间的数，表示概率1/i。 num = i; //以1/i的概率替换num的缓存值。 return num; }
对于这个问题的常见应用：马尔科夫文本生成器里选择哈希表中某项对应链表的任意一个结点。
