自从肖舸在其CSDN博客上说“拒绝回答自由互联等网站网友的问题”，实质上不单是拒绝回答，而且还删去包括一些网友及本人对于纯粹技术探讨的评论。当然每位博主都有自由这么做，但个人认为这对于社区的交流发展有负面影响。为了探讨这个技术问题，本人唯有把回应发表于自由互联内。本文会阐述之前的论点，评论肖舸的实现，并进行了兩个实验比较不同算法、实现的优劣之处。
之前的“交流”约一个月前，肖舸于《实际中常用的一个随机数产生器（分类别概率随机）》(下简称《实》)一文中，介绍了一个宣称“0bug”的实现。而实质上这实现至少有两个bug，其中一个已被网友发现，而肖舸也随后更新，不过没有公开向该网友致谢，连网友的id或名字也没写。本人还发现另一个会做成崩溃bug，稍后再说。
因为有感《实》欠缺背后理论，而且该实现的需求不清楚、代码难读，所以本人于5日后撰写博文《用JavaScript玩转游戏编程(一)掉宝类型概率》(下简称《用》 )，从统计学解释这问题，并以JavaScript实作互动的示范程序。文中提及:
这里用了线性搜寻(linear search)，……另外，也可以用二分搜寻(binary search)，那么复杂度会降低为O(lg N)，……那么，还有没有更快的方法呢？答案是肯定的，例如别名方法(alias method)、*似方法等，有兴趣的读者可参考……当然，在N很小的情况下，线性搜寻和二分搜寻也足够。
笔者撰写本文，灵感来自这篇博文(指《实》)。其算法实际上是储存CDF的逆函数采样，利用空间和有限的CDF精确度，换取O(1)的时间复杂度。衡量N的大小、精确度、空间需求、缓存延迟后，或许该方法也能适合某些个别需求。但对于该文作者说N最大为100，二分搜寻只需最多7次迭代，因缓存问题可能二分搜寻更快。
之后，肖舸在另一篇博文《做程序，要“专注”和“客观”》中，当中有一段似乎是回应上文:
就好比我前面写的一篇博文《实》，我在文中的代码里，明明实现了O(1)的复杂度，但是就有人，为了攻击我本人和我的书《……》，专门撰文，说用其他办法，O(7)的复杂度也可以实现，我这个办法不值得提倡。
我晕，我们做算法优化，有时候，O(n)这个值，能减少1都是巨大的成功，因为程序是有循环的，循环次数是被乘数哦，这是乘法关系，这个核心算法复杂度减少1，放出去就是几千万甚至几亿的时钟开销，效率提升就是巨大的。很多时候，我做优化，都在为了减少这个1在努力。
不过，这毕竟是少数人，准确的讲，说这话的人不能算技术人员，因为针对到科学的，算法的，优化的问题上，一是一，二是二，不能带着个人感情讨论。这是技术人员，特别是程序员基本的职业道德。
这里，我希望广大程序员朋友一定要养成一个习惯，“客观”和“严谨”是程序员的基本职业修养，也是我们能在这个行业里面立足的根本，千万不要丢掉了。
因为文中说“O(7)”，我想应该是回应我写的7次迭代，如果不是， 就当本人对号入座吧。对于算法是否属于“科学”，大家可以想想。我自问只谈到技术上的意见，不知道为什么会说到“个人感情”，甚至推理至是否一个技术人员、有没有职业道德问题了。
也许是本人写得不够清楚，也许是读者有误解，无论如何，本人也在此解释一下。
首先，肖舸说“O(n)这个值”、“ 核心算法复杂度减少1” 、“O(7)”，都是不正确的说法。
研究算法的运行时间，最常见是采用Big-O表示法，例如O(1)、O(n)等。这表示法是指，算法在n接*无限大的时候，其运行时间的渐*复杂度(asymptotic complexity)的上界(upper bound)。举个例子，例如解决同一问题的两个算法﹐，它们的运行时间(例如单位是秒)为:
\begin{align*} t_A(n) &= 2n + 1 \\ t_B(n) &= 9 \end{align*}用Big-O表示法，A算法是线性速度O(n)，B算法是常数速度O(1)。也就是说n大到某个程度，算法A比算法B慢。那么n要有多大？只要联合两个函数，即可求出，当n > 4时，算法A比B慢。相反，当n<4时，算法A会比算法B快。如果有另一组算法，其方程求出来的n是少于或等于0，说明当中的其中一个算法永远比另一个优。
使用Big-O表示法的目的，是方便比较不同的算法。 Big-O会隐藏的算法实现时的系数(例如上面例子中的2、1、9)。而在实际应用中，设计或实现算法，还要考虑多个因素，例如
- 算法的应用场合: 算法所需的精确度、时间要求(通常二者须此消彼长)。如算法分为几个步骤(如本问题可分割为初始化和采样)，每个步骤的实际使用次数。另外亦要考虑输入数据的分布。
- 内存的使用: 同时间内所需的内存最大值(可用Big-O表示法)、内存存取模式(access pattern)。
- 算法实现的难度: 除了开发及维护时间，还影响程序是否健壮(robust)，例如复杂的浮点运算会做成误差。
- 硬件: 例如主频、核心数量、内存速度、特殊指令(如SIMD)，或者不同架构的可编程系统(如GPU)、设备等
- 软件: 会执行于什么操作系统、编译器、库等等，这些软件对算法是否有影响。
很多时候，解决同一问题有许多算法，并没有绝对的好坏之分(当然有时候也有些能完全取胜，有些完全无用)，每个算法都其优点缺点。程序员须要分析实际应用，去选择及实现算法。理想地，更可以实验多个算法，用实际数据去作出比较。
评论肖舸的实现这个是肖舸在《实》公布的实现代码(本人加入了#ifdef MILO_HOTFIX及GetMemory()):
#define CTonyRandomArea_TOKEN_MAX 100 //最大类型数
#define CTonyRandomArea_TOKEN_AREA_MAX 10000 //类型数组单元数，精确到小数点后两位
//输入最大100个元素的数组，每个数组表示每类占有的百分比，内部自带百分比调整。
//即如果外部输入的数字之和不是整数100，内部会根据百分比，自动调整其比例，使总和=100.0
//然后内部建立10000个单元的类型数组，根据传入的每种类型的比例，在类型数组中批量填充对应的类型值
//总之，类型数组中每种类型的数量，占据的比例正好是输入的百分比
//最后，在0~10000中取随机，然后在对应的类型数组单元中取类型值，即为返回的类型
class CTonyRandomArea
{
public:
CTonyRandomArea(double* pTokenPercentArray,char cTokenCount)
{
m_nTokenCount=cTokenCount;
if(CTonyRandomArea_TOKEN_MAX<m_nTokenCount)
m_nTokenCount=CTonyRandomArea_TOKEN_MAX;
int i=0;
for(i=0;i<m_nTokenCount;i++)
{
m_dTokenPercentArray[i]=*(pTokenPercentArray+i);
}
//动态调整内部的值
//有时候试验人员，测得几个状态出现的数字，可能懒得再计算成百分比
//程序帮忙自动计算
double dNumberCount=0;
for(i=0;i<m_nTokenCount;i++)
{
dNumberCount+=m_dTokenPercentArray[i];
}
if(100.0!=dNumberCount)
{
for(i=0;i<m_nTokenCount;i++)
{
m_dTokenPercentArray[i]/=dNumberCount;
m_dTokenPercentArray[i]*=100;
}
}
//以小数点后两位精度，开始计算在10000个总单元中，每种类型对应的数量。
for(i=0;i<m_nTokenCount;i++)
{
m_sTokenPercentArray[i]=(short)(m_dTokenPercentArray[i]*100);
}
//按比例填充类型数组
int j=0;
int nFillMin=0;
int nFillMax=0;
#ifdef MILO_HOTFIX
for(i=0;i<CTonyRandomArea_TOKEN_AREA_MAX;i++)
{
m_cTokenPercentArrayAreaUp[i]=0;
}
#else
for(i=0;i<m_nTokenCount;i++)
{
m_cTokenPercentArrayAreaUp[i]=-1;
}
#endif
for(i=0;i<m_nTokenCount;i++)
{
nFillMax=nFillMin+m_sTokenPercentArray[i];
for(j=nFillMin;j<nFillMax;j++)
{
m_cTokenPercentArrayAreaUp[j]=i;
}
nFillMin=nFillMax;
}
// m_cTokenPercentArrayAreaUp[CTonyRandomArea_TOKEN_AREA_MAX-1]=i-1;
}
~CTonyRandomArea(){}
void PrintfInfo(void)
{
int i=0;
double dNumberCount=0;
for(i=0;i<m_nTokenCount;i++)
{
dNumberCount+=m_dTokenPercentArray[i];
printf("%d = %f\n",i,m_dTokenPercentArray[i]);
}
printf("All = %f\n",dNumberCount);
//打印10000个单元的类型分布，看看排得对不对
//这段打印起来太长，一般隐掉，需要了再打印
// int nOutMax=400;
// int nInMax=25; //二者相乘为10000
// int j=0;
// for(i=0;i<nOutMax;i++)
// {
// printf("%05d - ",i*nInMax);
// for(j=0;j<nInMax;j++)
// {
// printf("%d ",m_cTokenPercentArrayAreaUp[i*nInMax+j]);
// }
// printf("\n");
// }
}
//取类型数组对应单元的值
char GetType(int nTypeIndex) //输入参数0~10000
{
if(10000<=nTypeIndex) nTypeIndex=9999;
if(0>nTypeIndex) nTypeIndex=0;
return m_cTokenPercentArrayAreaUp[nTypeIndex];
}
//真实的工作函数，利用输入的概率来产生随机type
char GetRandomType(void)
{
return GetType(GetRandomBetween(0,10000));
}
// MILO ADD {
size_t GetMemory() const {
return sizeof(*this);
}
// } MILO
private:
double m_dTokenPercentArray[CTonyRandomArea_TOKEN_MAX]; //比例数组
short m_sTokenPercentArray[CTonyRandomArea_TOKEN_MAX]; //比例数组,短整型，1~10000，相当于精确到小数点后两位
char m_nTokenCount;
char m_cTokenPercentArrayAreaUp[CTonyRandomArea_TOKEN_AREA_MAX]; //类型数组
};
////这是测试代码
//void TestCTonyRandomArea(void)
//{
// double dTokenPercentArray[100];
//
// dTokenPercentArray[0]=12.0;
// dTokenPercentArray[1]=40.0;
// dTokenPercentArray[2]=40.0;
// dTokenPercentArray[3]=7.0;
// dTokenPercentArray[4]=1.0;
//
// CTonyRandomArea Area1(dTokenPercentArray,5);
// Area1.PrintfInfo();
//
// int i=0;
// for(i=0;i<20;i++)
// {
// printf("RandType = %d\n",Area1.GetRandomType());
// }
//}
崩溃Bug
在测试这个实现时，发现GetRandomType()返回的值会>=N，而N为概率数组的大小。这样有机会做成程序崩溃。例如在代码中，把返回值作为索引，填进数组就会写到不合法的地址。
调试之后，发现程序的bug位于三个地方。第一，21-36行检查概率总和是否等于100个百分点，否则会自动进行正规化(normalization)。正规化后，浮点小数只能表示*似值。第二，在37-41行会把浮点小数乘上10000，再转做整数。因为这里是下取整，计算出来的整数总和有机会少于10000。于是m_cTokenPercentArrayAreaUp数组最后的几个元素未被填入。第三，第54-57行对m_cTokenPercentArrayAreaUp数组进行初始化，但用错循环次数，应为CTonyRandomArea_TOKEN_AREA_MAX而不是m_nTokenCount。
例如使用{ 3, 5, 7, 11, 13 } 作为输入，CTonyRandomArea_TOKEN_AREA_MAX数组就会填入前9998个元素，最后两个未被初始化。
这三个错处要完全更正会改动太多，不如重写。所以本人只是作最少改动，原整地初始化CTonyRandomArea_TOKEN_AREA_MAX数组为0。
用肉眼检验正确性？测试函数TestCTonyRandomArea()的作用就是列印20个采样。该程序员或测试人员是用肉眼观测那20个采样，看看是否大概和输入的概率分布接*么？或是记下每个返回值的频率，笔算其频率除以20后的百分比？
这个问题本质上就是能编写测试代码去自动测试的。以下的实验就会测试各个实现的特性，包括其误差。
浪费内存m_dTokenPercentArray和m_sTokenPercentArray都没必要作为成员变量。前者在PrintfInfo()函数还有一点用途，后者只在构造函数使用。
const-ness构造函数的指针没加const。使用者用的时候要const_cast。 C语言的函数都要适当地使用const指针(可参考strlen原型)，何况是C++呢。
乱数产生器不均分布以下是其的乱数产生器代码:
inline int _GetNot0(void)
{
int nRet=rand();
if(!nRet) nRet++;
return nRet;
}
inline int GetRandomBetween(int nBegin,int nEnd)
{
int n=_GetNot0();
int nBetween=0;
if(0>nBegin) nBegin=-nBegin;
if(0>nEnd) nEnd=-nEnd;
if(nBegin>nEnd)
{
nBetween=nEnd;
nEnd=nBegin;
nBegin=nBetween;
}
else if(nBegin==nEnd)
nEnd=nBegin+10;
nBetween=nEnd-nBegin;
n=n%nBetween;
n+=nBegin;
return n;
}
这个在书评里已写了很多，原来只需三句代码都可以写成这样，更把乱数的分布变得更不*均，不再多说了。
算法实现本来以为《用》一文中的说法，已足够解释。以上的情况，加上本人未试过实现别名方法(alias method)，便做了二个实验去比较及分析各种算法。实验有其局限性，只是在某个客观环境(计算机、编译器等)，对某个算法实现的测量。但本实验的结果可能可以支持我在《用》所提及的观点。
逆变换方法(线性搜寻)这个算法于《用》有详细说明，以下是其C++代码，复杂度是O(N):
template <typename T, typename RandomGenerator>
class InverseLinear {
public:
InverseLinear(const T* cdf, size_t N, RandomGenerator& random) : mCDF(cdf), mN(N), mRandom(random) {
assert(N > 0);
assert(cdf[N - 1] == (T)1.0);
}
size_t operator()() {
T y = mRandom();
for (size_t x = 0; x < mN - 1; x++)
if (y < mCDF[x])
return x;
return mN - 1;
}
size_t GetMemory() const {
return sizeof(*this) + mN * sizeof(T);
}
private:
const T* mCDF;
const size_t mN;
RandomGenerator mRandom;
};
逆变换方法(二元搜寻)
这个算法用二元搜寻代替线性搜寻，复杂度降低为O(lg N):
template <typename T, typename RandomGenerator>
class InverseBinary {
public:
InverseBinary(const T* cdf, size_t N, RandomGenerator& random) : mCDF(cdf), mN(N), mRandom(random) {
assert(N > 0);
assert(cdf[N - 1] == (T)1.0);
}
size_t operator()() {
T y = mRandom();
size_t left = 0, right = mN;
while (left < right - 1) {
size_t mid = (left + right) / 2;
assert(mid - 1 < mN);
if (mCDF[mid - 1] < y)
left = mid;
else
right = mid;
}
return left;
}
size_t GetMemory() const {
return sizeof(*this) + mN * sizeof(T);
}
private:
const T* mCDF;
const size_t mN;
RandomGenerator mRandom;
};
必须注意这个算法和一般的二元搜寻条件有些差异。
别名方法别名方法能提供O(1)的复杂度，但初始化过程比较复杂。参照[1]附录的代码实作:
template <typename T, typename RandomGenerator>
class Alias {
public:
Alias(const T* pdf, size_t N, RandomGenerator& random) : mN((T)N), mRandom(random), mAliasTable(N) {
assert(N > 0);
std::vector<size_t> smallBag, largeBag;
smallBag.reserve(N);
largeBag.reserve(N);
const T stdWeight = (T)1 / N;
for (size_t i = 0; i < N; i++) {
mAliasTable[i].prob = pdf[i];
if (pdf[i] > stdWeight)
largeBag.push_back(i);
else
smallBag.push_back(i);
}
while (!largeBag.empty() && !smallBag.empty()) {
const size_t sm = smallBag.back();
const size_t lg = largeBag.back();
smallBag.pop_back();
largeBag.pop_back();
mAliasTable[lg].prob -= stdWeight - mAliasTable[sm].prob;
mAliasTable[sm].prob *= (T)N;
mAliasTable[sm].index = lg;
if (mAliasTable[lg].prob > stdWeight)
largeBag.push_back(lg);
else
smallBag.push_back(lg);
}
for (std::vector<size_t>::iterator itr = smallBag.begin(); itr != smallBag.end(); ++itr)
mAliasTable[*itr].prob = (T)1;
for (std::vector<size_t>::iterator itr = largeBag.begin(); itr != largeBag.end(); ++itr)
mAliasTable[*itr].prob = (T)1;
}
__forceinline size_t operator()() {
T u = mRandom() * mN;
size_t i = (size_t)u;
return u - i > mAliasTable[i].prob ? mAliasTable[i].index : i;
}
size_t GetMemory() const {
return sizeof(*this) + sizeof(AliasItem) * mAliasTable.capacity();
}
private:
T mN;
RandomGenerator mRandom;
struct AliasItem {
T prob;
size_t index;
};
std::vector<AliasItem> mAliasTable;
};
伪随机数产生器
Visual C++的C运行时库(C runtime library, CRT)中，提供的rand()是线程安全的，因而有额外的开销(需要存取TLS)。而且，其精确度有限(只有15-bit)。因此，在实验中使用了两个伪乱数产生器，一个采用rand()，一个采用最普通的线性同余产生器(Linear congruential generator, LCG)。两个产生器都是产生半合区间[0,1)的匀均分布。
template <typename T>
class RandomCRT {
public:
RandomCRT(unsigned seed = 0) {
srand(seed);
}
T operator()() {
return rand() * ((T)1 / (RAND_MAX + 1));
}
};
template<typename T>
class RandomLCG {
public:
RandomLCG(unsigned seed = 0) : mSeed(seed) {}
T operator()() {
mSeed = 214013 * mSeed + 2531011;
return mSeed * ((T)1 / 4294967296);
}
private:
unsigned mSeed;
};
实验一: 方法
主要的测试方式如下。设N为4, 8, 16, ..., 256，用乱数产生一组N个元素的pdf，之后对每个算法，执行其初始化1,000次，采样1,000,000次，分别量度时间。此外，用采样来计算频表，之后用以下算式比较频表和pdf的总误差:
Error(x, h)=\sum_{0}^{N-1}{ \| h_i - P(X=x_i) \| }template <typename T>
T ComputeError(const unsigned* histogram, const T* pdf, size_t N, size_t sampleCount) {
T error = 0.0;
for (size_t i = 0; i < N; i++) {
T delta = (T)histogram[i] / sampleCount - pdf[i];
error += abs(delta);
}
return error;
}
当中h_i是采样等于i的出现的次数除以总采样次数。
有些算法分别采用float和double两种类型，亦会使用CRT和LCG两种伪乱数产生器，所有测试组合如下:
- Tony
- Linear<CRT, double>
- Binary<CRT, double>
- Alias<CRT, double>
- Linear<LCG, double>
- Binary<LCG, double>
- Alias<LCG, double>
- Linear<LCG, float>
- Binary<LCG, float>
- Alias<LCG, float>
实验环境:
- Intel i7 920 @2.67Hz
- Microsoft Windows 7 64-bit
编译器及主要编译参数:
- Microsoft Visual Studio 2008 9.0.30729.1 SP
- Platform: Win32
- Optimization: Maximize Speed (/O2)
- Inline Function Expansion: Any Suitable (/Ob2)
- Enable Intrinsic Function: Yes (/Oi)
- Favor Size of Speed: Favor Fast Code (/Ot)
- Enable C++ Exception (No)
- Runtime Library: Multi-threaded DLL (/MD)
- Buffer Security Check: No (/GS-)
- Enable Enhanced Instruction Set: Streaming SIMD Extensions 2 (/arch:SSE2)
- Floating Point Model: Fast (/fp:fast)
- Default Char Unsigned: Yes (/J)
- Enable Run-Time Type Info: No (/GR-)
宏:
- #define _SECURE_SCL 0
以下结果皆使用Google Visualization API显示，因此在RSS或转载中可能无法显示。我初次试用这API，其HTML和JavaScript代码是由C++测试程序中自动生成的。
原始数据表如下(init和sample的单位为秒，memory的单位为字节):
初始化时间纵轴是时间，单位为秒。注意横轴为对数比例。
除了Tony和Alias，其余算法的初始化都接*零。 Tony和Alias都是以线性上升。大部份情况下，Tony的初始化性能较低。
采样时间最好成绩的是Alias的LCG版本，接*常数性能，其float和double版本差异不大。 Tony亦为常数性能，比Alias慢30%左右。
一如所料，Linear和Binary的性能为O(N)和O(lg N)(因横轴为对数比例，对数在图里会接*直线，例如橙线的Binary)。值得注意的是，Linear和Binary的LCG版本在N=4时，超越Alias和Tony。这是由于这两个算法简单，系数低，所以有机会超越O(1)的对手。
需要较高性能的话，别用CRT的rand()。
准确度如果要说Tony算法的主要问题，就在于其准确度远远差于其他算法的实现(在数百倍以上)。就算使用整数的百份比，用那N=5的例子，Tony的误差约为0.08，其他算法则介乎0.001至0.002。
内存用量Tony算法所使用的内存是常量11008个字节。而其他的算法为O(N)，Alias因为要多存一个别名索引，其内存是其他除Tony以外的两倍(事实上，Alias的索引可以因N的大小而采用unsigned short或unsigned char，以减低开销)。
这里亦要留意一点，Linear和Binary其实只储存cdf的常数指针、N和Random，实际上只有12字节。 cdf数组可供多个采样器共用，所以Linear和Binary的内存开销是非常少的。
N的范围Tony算法的另一缺点是只支持至N=100。相反，其他的算法可支持较大的N，例如用Tony所需的內存，Alias的float版本可支持N=1375。
实验二: 方法上一个实验，量度各个算法的采样器在不同的N时的性能。但是在应用上，很多时候需要同时使用多个采样器。肖舸在《实》中提及:
实话跟你讲吧，这段代码是有前提的，我们要做5000万条记录，中间有20万个设备的记录，每个设备的采样频率不一样，我要并发模拟，你再想想我写这么复杂有道理没？
那么就实际测试一下。建立M个采样器(M设为1、10、...、100000)，对每个采样器轮流采样，M个采样器合共采样1百万次。而N则取Tony算法可接受的最大值100。测试的循环代码是这样的:
const size_t interleavedSample = sampleCount / M; size_t unused = 0; // prevent optimized out LARGE_INTEGER start, end; QueryPerformanceCounter(&start); for (size_t i = 0; i < interleavedSample; i++) for (size_t j = 0; j < M; j++) unused += (*samplers[j])(); QueryPerformanceCounter(&end); cout << unused << endl;
从理论上来说，用一个采样器去做一百万次采样，和用M个采样器合共采样一百万次，在效能上应该没有分别。但在实际的测试中又是否这样呢？
因为上一个实验已测量准确度，以及那些算法组合在这个情况下比较好，就只测试以下的组合:
- Tony
- Binary<LCG, float>
- Alias<LCG, float>
原始数据如下:
用图表显示采样性能(横轴为M個采样器，纵轴是采样1百万次时间，单位为秒):
可以发现，Alias在整个测试范围里，表现也是最佳的。若按理论内说，三条线应该都是*的，或接**的。为什么O(1)的Tony算法却在M=1000以后，输给O(lg N)的Binary？而且，如果用多个采样器是有额外开销的话，为何Alias和Binary不是一直随M上升而线性递增，反而会是个U形，有个甜蜜点(sweet spot)？
Tony算法除了准确度问题，最大问题是其内存开销。当同时使用多个采样器，其存取的内存超过CPU缓存空间，就会严重影响性能。这可以作为例子证实在某些情况下，我在《用》里说的假设:
但对于该文作者说N最大为100，二分搜寻只需最多7次迭代，因缓存问题可能二分搜寻更快。
至于甜蜜点，在实验之前也没预计过。我认为这也是由缓存做成的。因为缓存本身也有*行运作的特性，分散存取整个缓存空间，会比集中存取一小块缓存更有效率。观看上表，可发现甜蜜点的内存使用大约在40K-80K，这似乎刚能对上i7的64K L1缓存。但本人对硬件了解不足，还望各位网友指导。
那么在M=100,000之后，最终Tony是否又会再超越Binary呢？从理论上来说，只要两个算法都使用超过缓存大小的内存(以i7 920来说是8MB的L3)，缓存就会慢慢失去效果。可是我没有再尝试，因为在M=100,000，Tony算法已使用了超过1.1GB内存，而Binary才约41MB。回想起来，肖舸说要做20万个设备，每个有不同的采样率。如果要同时使用，Tony算法要2.2GB内存啊，而Binary才82MB，后者比前者简单、准确、快。
总结本文使用实验比较了4个离散随机变量的产生算法。实验量度的主要四种数据，可以帮助读者选择合适的算法。在于N比较大的情况下，建议使用别名方法，因为其采样性能最高、内存比Tony少、准确度比Tony高，只是初始化性能稍慢。而本人参照[1]的别名方法实现，也可以尝试进一步优化，例如用一个alloca()代替两个vector，减少配置内存的开销。
本人认为，使用类似Tony的算法，在某些情况是有用的。例如，可以设计一个采样器，其输入为整数，代表每个类别的相对比重。例如输入整数数组{ 1, 2, 3, 1, 2 }，表示概率{ 1/9, 2/9, 1/3, 1/9, 2/9 }，那么只需要9个元素的表就能实现快速的采样。这样的输入，大概十多行代码就能实现，而且准确度非常高，也不需浮点运算。缺点(或特征)是其复杂度与数组的元素总和成线性关系。这个实现留给读者尝试吧。
希望本文也能带出几个要点。在实际应用中，算法不能单靠Big-O时间复杂度来选择。以采样器为例，需同时考虑初始化和采样的时间。例如每次建立一个采样器之后，会执行50次采样，那么就可以比较各种算法的实际执行时间，也许Binary会比Alias好。另外，程序的准确性通常比性能更重要(当然也有例外，例如游戏中的粒子系统效果就不需要很准确)，但无论要求如何，也应谨慎对待准确性。除了嵌入式设备有内存限制，因现代计算机的内存速度和缓存缺失(cache miss)，编程时应尽量减少内存的开销。最后，也请注意不同算法，代码的复杂程度不尽相同。以上所说，皆是选择算化的重要考虑因素。
本人是否为一名技术人员、是否有职业道德，就不在此回应了。
有与趣的读者可下载本实验的源代码包(不保証0bug)。此外，C++ TR1有部份关于随机数的实现。
参考- [1] Burke D, Bidirectional importance sampling for illumination from environment maps. Master's Thesis, University of British Columbia, 2004.
- 2010-05-28: 肖舸在其CSDN博客上发表《请自由互联立即停止侵权行为的公开信》。本人认为本文只是适当引用肖舸的文章，为介绍、评论其文章及说明问题，并在发表时已加上原文名称、链结及原作者名字。因此，按法律并不需得到原作者同意。
- 2010-05-29: 推友@DOTIAN做了一个关于别名方法(alias method)的详细分析。本文主要在性能比较，没有介绍别名方法原理，读者可参考这文章。
- 2010-06-02: 學生大本營中"編程之美"老師轉載本文。
- 2010-06-06: 肖舸刪去其CSDN博客所有文章。此非吾所願。
- 2010-06-10: 肖舸的CSDN博客和學生大本營帪戶(被?)關閉。CSDN要求其他幾位學生大本營老師刪除有關
- 2010-06-11: 學生大本營中"編程之美"老師轉載的本文被刪。