高斯随机数发生程序

长平狐 发布于 2012/08/13 15:44
阅读 112
收藏 0

附送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://blog.csdn.net/favormm/article/details/5707959
加载中
返回顶部
顶部