梅森旋转算法原理c语言,梅森旋转素数算法（MT199937）c语言代码

满意答案

sddzycz

2016.06.23

采纳率：42%    等级：12

已帮助：8913人

#include /*Periodparameters*/

#define N 624

#define M 397

#define MATRIX_A 0x9908b0dfUL /*constantvectora*/

#define UPPER_MASK 0x80000000UL/*mostsignificantw-rbits*/

#define LOWER_MASK 0x7fffffffUL/*leastsignificantrbits*/

static unsigned long mt[N];/*thearrayforthestatevector*/

static int mti=N+1;/*mti==N+1meansmt[N]isnotinitialized*/ /*initializesmt[N]withaseed*/

void init_genrand(unsigned long s)

{

mt[0]=s & 0xffffffffUL;

for(mti=1;mti

{

mt[mti]=

(1812433253UL*(mt[mti-1]^(mt[mti-1]>>30))+mti); /*SeeKnuthTAOCPVol2.3rdEd.P.106formultiplier.*//*Inthepreviousversions,MSBsoftheseedaffect*//*onlyMSBsofthearraymt[].*//*2002/01/09modifiedbyMakotoMatsumoto*/

mt[mti]&=0xffffffffUL;/*for>32bitmachines*/

}

} /*initializebyanarraywitharray-length*//*init_keyisthearrayforinitializingkeys*//*key_lengthisitslength*/ /*slightchangeforC++,2004/2/26*/

void init_by_array(unsigned long init_key[],int key_length)

{

int i,j,k; init_genrand(19650218UL);

i=1;j=0;

k=(N>key_length?N:key_length);

for(;k;k--)

{

mt[i]=(mt[i]^((mt[i-1]^(mt[i-1]>>30))*1664525UL))

+init_key[j]+j;/*nonlinear*/

mt[i]&=0xffffffffUL;/*forWORDSIZE>32machines*/

i++;j++; if(i>=N){mt[0]=mt[N-1];i=1;}if(j>=key_length)j=0;

}

for(k=N-1;k;k--)

{ mt[i]=(mt[i]^((mt[i-1]^(mt[i-1]>>30))*1566083941UL)) -i;/*nonlinear*/

mt[i]&=0xffffffffUL;/*forWORDSIZE>32machines*/i++;

if(i>=N){mt[0]=mt[N-1];i=1;} }

mt[0]=0x80000000UL;

/*MSBis1;assuringnon-zeroinitialarray*/ } /*generatesarandomnumberon[0,0xffffffff]-interval*/

unsigned long genrand_int32(void)

{

unsigned long y;

static unsigned long mag01[2]={0x0UL,MATRIX_A};/*mag01[x]=x*MATRIX_Aforx=0,1*/

if(mti>=N){/*generateNwordsatonetime*/ int kk;

if(mti==N+1)/*ifinit_genrand()hasnotbeencalled,*/

init_genrand(5489UL);/*adefaultinitialseedisused*/

for(kk=0;kk

{

y=(mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);mt[kk]=mt[kk+M]^(y>>1)^mag01[y&0x1UL];

} for(;kk

{ y=(mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);

mt[kk]=mt[kk+(M-N)]^(y>>1)^mag01[y&0x1UL];

}

y=(mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);

mt[N-1]=mt[M-1]^(y>>1)^mag01[y&0x1UL];mti=0; }

y=mt[mti++]; /*Tempering*/

y^=(y>>11); y^=(y<<7)&0x9d2c5680UL;y^=(y<<15)&0xefc60000UL;y^=(y>>18);

return y;

} /*generatesarandomnumberon[0,0x7fffffff]-interval*/

long genrand_int31(void)

{ return(long)(genrand_int32()>>1); } /*generatesarandomnumberon[0,1]-real-interval*/

double genrand_real1(void)

{ return genrand_int32()*(1.0/4294967295.0); }

double genrand_real2(void)

{ return genrand_int32()*(1.0/4294967296.0);/*dividedby2^32*/ } /*generatesarandomnumberon(0,1)-real-interval*/

double genrand_real3(void)

{ return(((double)genrand_int32())+0.5)*(1.0/4294967296.0);/*dividedby2^32*/ }

/*generatesarandomnumberon[0,1)with53-bitresolution*/

double genrand_res53(void)

{ unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;

return(a*67108864.0+b)*(1.0/9007199254740992.0); } /*TheserealversionsareduetoIsakuWada,2002/01/09added*/

int main(void){

int i; unsigned long init[4]={0x123,0x234,0x345,0x456},length=4;init_by_array(init,length); printf("1000outputsofgenrand_int32()\n");

for(i=0;i<1000;i++)

{

printf("%10lu",genrand_int32());

if(i%5==4)printf("\n"); }

printf("\n1000outputsofgenrand_real2()\n");

for(i=0;i<1000;i++){ printf("%10.8f",genrand_real2());

if(i%5==4)printf("\n"); }

return 0;

}

00分享举报