回应CSDN肖舸《做程序，要“专注”和“客观”》，实验比较各离散采样算法

自从肖舸在其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: 學生大本營中"編程之美"老師轉載的本文被刪。