高斯随机数发生程序

附送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 512DWORD 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，方差为1FTable: 正态分布函数表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;}