高斯随机数发生程序

  附送16bits,32bits均匀随机数发生程序

程序只产生均值为0，方差为1的随机数，要产生均值为E,方差为D的随机数，只要
随机数*D+E就可以了。

高斯随机数程序还是带参数的，参数是用来描述正态分布的一个浮点数表。在执行
程序时，先生成这个表（只做一次），而后就可以任意多次地执行高斯随机数产生
程序了。
使用C是为了保证通用性，如果有人觉得麻烦，可以用C++做个类，把这些东西都封
装进去。
另外，如果有人有兴趣，也可以把它修改成任意形式分布的连续随机数产生程序，
修改非常简单，这里就不提示了。

#include <memory.h>
#include <time.h>
#include <math.h>


#define GAUSS_TABLE_LENGTH 256

#include "random.h"

extern float Gauss_F[GAUSS_TABLE_LENGTH];
extern float Gauss_f[GAUSS_TABLE_LENGTH];
extern DWORD RandomKey,RandomKey32;

DWORD MCoef_32[2]={0xE7BD2160,0xDA3A2A9C};
WORD MCode_16[4]={0xD31C,0xC36C,0xD0A2,0xD228} ;
// two m sequence

#define INVALID_BLOCK 0xFFFF
#define BLOCK_SIZE 512

DWORD m_Seq_32(DWORD Key)
{
int i;
for(i=0;i<32;i++)
{
_asm
{
MOV EBX,Key;
SHL EBX,1
MOV EAX,Key;
MOV EDX,MCoef_32[0];
AND EAX,EDX; //select the bit for xor
MOV EDX,EAX;
SHR EAX,16;
XOR AX,DX;
XOR AH,AL; // because P only judge one byte
// so must XOR to judge the p of whole word
JP NEXT //jp equals the xor
INC EBX
NEXT: MOV Key,EBX;

}
}
return Key;
}

WORD RNG()
{
DWORD A,B;
_asm
{
_emit 0x0f
_emit 0x31
MOV A,EAX
MOV B,EDX
}

RandomKey = m_Seq_32(RandomKey^A^B);
return (WORD)RandomKey;
}

DWORD RNG32()
{
DWORD A,B;
_asm
{
_emit 0x0f
_emit 0x31
MOV A,EAX
MOV B,EDX
}

RandomKey32 = m_Seq_32(RandomKey32^A^B);
return RandomKey32;
}

void NormalTable(float *FTable, float * fTable,int Length)
/*
生成正态分布函数和正态分布概率密度函数表，表宽度为3倍方差
正态分布函数均值为0，方差为1
FTable: 正态分布函数表
fTable: 正态分布概率密度表
Length: 表的长度
*/

{
int i;
float h; /* 步长 */
float x,temp;
float C;

x=-3;
h= 6.0/Length/2;
FTable[0]=0;
C=1/sqrt(2*3.1415927);
/* 初始参数设置 */

fTable[0]= exp(-x*x/2)*C;
/* 起始点的概率密度 ，exp(-x*x/2)*C为概率密度函数 */

for(i=1;i<Length;i++)
{
x+=h;
temp = exp(-x*x/2)*C;
x+=h;
fTable[i] = exp(-x*x/2)*C;
/* 计算正态分布概率密度函数 */

FTable[i] = FTable[i-1]+(fTable[i-1]+4*temp+fTable[i])*h/3;
/* 辛普森数值积分公式计算正态分布函数 */
}
}

float NormalRNG(float * Gauss_F, int Length)
/*
生成均值为0,方差为1的高斯随机数
Gauss_F: 正态分布函数表，规格如上，可以由NormalTable函数生成
Length: 正态分布函数表的长度
返回值：均值为0,方差为1的高斯分布随机数
*/

{
float RandomNumber;
float temp,h;
int i;

h= 6.0/Length;
/* 正态分布表的步长 */


temp = (float)RNG()/65536;
/*产生一个[0,1)区间内均匀随机数 */

if(temp == 0)
return 0;

for(i=0;i<Length;i++)
{
if(temp<=Gauss_F[i]) /*计算随机数落入正态分布表的哪个区间 */
{
RandomNumber = (-3+h*i)+(temp-Gauss_F[i-1])/(Gauss_F[i]-Gauss_F[i-1])*h;
/* 进行线性差值 */
break;
}
}
return RandomNumber;
}

本文地址：http://com.8s8s.com/it/it1883.htm