HELP!!!

[egao]

正态分布函数的C/C++算法或C/C++程序。谢！

[酒红绿叶]

double Restore_ran_lp(long *idum)
{
/*--------------------------------------------------------------------------
Long period(>2x10^18) random number generator of L'Ecuyer
with Bays-Durhan shuffle and added safeguard. Return a uniform
random deviate between 0.0 and 1.0 (exclusive of the endpoint
value). Call with idum a negative integer to initialize;
thereafter, do not alter idum between successive deviates in
a sequence. RNMX should approximate the largest floating value
that is less than 1.
--------------------------------------------------------------------------*/

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-27
#define RNMX (1.0-EPS)


int j;
long k;
static long idum2=123456789;
static long iy=0;
static long iv[32];
double temp;

if(*idum<= 0)
{
if(-(*idum)<1)
*idum=1;
else *idum=-(*idum);

for(j=NTAB+7;j>=0;j--)
{
k=(*idum)/IQ1;
*idum=IA1*(*idum-k*IQ1)-IR1*k;
if(*idum<0)
*idum +=IM1;
if(j<NTAB)
iv[j]=*idum;
}

iy = iv[0];
}

k=(*idum)/IQ1;

*idum = IA1*(*idum-k*IQ1)-IR1*k;

if( *idum<0 )
*idum += IM1;

k = idum2/IQ2;

idum2= IA2 *( *idum - k*IQ2 ) - IR2*k;

if( idum2 < 0 )
idum2 += IM2;

j = iy / NDIV;

iy = iv[j] - idum2;

iv[j] = *idum;

if( iy<1 )
iy += IMM1;

if( ( temp = AM*iy ) > RNMX )
return RNMX;
else
return temp;

#undef IM1
#undef IM2
#undef AM
#undef IMM1
#undef IA1
#undef IA2
#undef IQ1
#undef IQ2
#undef IR1
#undef IR2
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX
}


double Restore_ran_gauss(long *idum)
{
/*----------------------------------------------------------------------------
returns a normally distributed deviate with zero mean and
unit variance, using ran_lp() as the source of uniform deviates.
----------------------------------------------------------------------------*/

#define RAN Restore_ran_lp /* change this to use other random number generator */

static char iset=0;
static double gset;
double fac,rsq,v1,v2;

if(iset==0)
{
do
{
v1=2.0*RAN(idum)-1.0;//2.0*RAN(idum)-1.0;
v2=2.0*RAN(idum)-1.0;//
rsq=v1*v1+v2*v2;
}
while(rsq>=1.0 || rsq==0.0);

fac=sqrt(-2.0*log(rsq)/rsq);
gset=v1*fac;
iset=1;

return v2*fac;
}
else
{
iset=0;
return gset;
}
#undef RAN
}

不要问为什么！
我也不知道！

[egao]

我看着都晕！？
谁可以解释？！

[byry]

谁看谁精神!