www.pudn.com > Genecircus20070919.rar > Random.cpp

#include "Random.h" 
 
CRandom::CRandom(void) 
{ 
} 
 
CRandom::~CRandom(void) 
{ 
} 
const double PI = 3.14159265358979324; 
 
void CRandom::moment(history data,int n,double & ave,double & adev,double & sdev,double & var,double & skew,double & curt) 
{ 
	 int j; 
	 if(n<=1){return;} 
	 double p,s=0.0; 
	 for(j =0;j= (en + 1.0)); 
            em = ::floor(em); 
            t = 1.2 * sq * (1.0 + y * y) * 
                ::exp(oldg - log_gamma(em + 1.0) - 
                      log_gamma(en - em + 1.0) + 
                      em * plog + (en - em) * pclog); 
        } while (unidev() > t); 
        bnl = em; 
    } 
    if (p != pp) 
        bnl = n - bnl; 
    return bnl; 
} 
 
double CRandom::log_gamma(double xx) 
{ 
    double x, y, tmp, ser; 
    static double cof[6] = {76.18009172947146, 
                            -86.50532032941677, 
                            24.01409824083091, 
                            -1.231739572450155, 
                            0.1208650973866179e-2, 
                            -0.5395239384953e-5}; 
    y = x = xx; 
    tmp = x + 5.5; 
    tmp -= (x + 0.5) * ::log(tmp); 
    ser = 1.000000000190015; 
    for (int j = 0; j <= 5; ++j) 
        ser += cof[j] / ++y; 
    return -tmp + ::log(2.5066282746310005 * ser / x); 
} 
 
static double r[98]; 
static int iff,ix1,ix2,ix3; 
 
double CRandom::rand_uniform() 
{ 
	int idum = -(::wxGetElapsedTime()); 
    double rm1,rm2,t; 
	int m1,m2,m3,ia1,ia2,ia3,ic1,ic2,ic3,j; 
	m1 = 259200,ia1 = 7141,ic1 = 54773,rm1 = 0.0000038580247; 
	m2 = 134456,ia2 = 8121,ic2 = 28411,rm2 = 0.0000074373773; 
	m3 = 243000,ia3 = 4561,ic3 = 51349; 
	if((idum<0)||(iff == 0)){ 
	    iff = 1; 
		ix1 = (ic1-idum)%m1; 
		ix1 = (ia1*ix1 + ic1)%m1; 
		ix2 = ix1%m2; 
		ix1 = (ia1 * ix1 + ic1)%m1; 
		ix3 = ix1 %m3; 
		for(j = 1;j<= 97;j++){ 
			ix1 = (ia1*ix1+ic1)%m1; 
			ix2 = (ia2*ix2+ic2)%m2; 
			r[j] = (double(ix1)+double(ix2)*rm2)*rm1; 
		} 
		idum = 1; 
	} 
	ix1 = (ia1*ix1+ic1)%m1; 
	ix2 = (ia2*ix2+ic2)%m2; 
	ix3 = (ia3*ix3+ic3)%m3; 
	j = 1 + int((97*ix3)/m3); 
	if((j>97)||(j<1)){ 
		return 1; 
	} 
	t = r[j]; 
	r[j] = (float(ix1)+float(ix2)*rm2)*rm1; 
	return t; 
}