梅森旋转素数算法(MT199937)c语言代码?
#include<cstdio> /*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<N;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<N-M;kk++)
{
y=(mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);mt[kk]=mt[kk+M]^(y>>1)^mag01[y&0x1UL];
} for(;kk<N-1;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;
}