writenoise.cpp

用c/c++产生源码. 本人平时收集的. 主要用于算法仿真. Park and Miller with Bays-Durham shuffle and added safeguards

#ifndef _RANDOM_H
#define _RANDOM_H

int seed;
void setSeed(int s)
{
	seed=s;
}

/* L'Eculer with Bays-Durham shuffle and added safeguards.
Return a uniform random deviate between 0.0 and 1.0(exclusive of the endpoint values).
Call with seed a negative integer to initialize; thereafter, do not alter seed
between successive deviates in a sequence. RNMX should approaimate the largest floating
value that is less than 1.
this algorithm is ok when the required random numbers large than 100,000,000.*/

double ranLE(int &seed)
{
	const int IM1=2147483563, IM2=2147483399;
	const int IA1=40014,IA2=40692,IQ1=53668,IQ2=52774;
	const int IR1=12211,IR2=3791,NTAB=32,IMM1=IM1-1;
	const int NDIV=1+IMM1/NTAB;
	const double EPS=3.0e-16,RNMX=1.0-EPS,AM=1.0/double(IM1);
	static int seed2=123456789,iy=0;
	static int iv[NTAB];
	int j,k;
	double temp;
	
	if(seed<0) //initialize
	{
	   seed=(seed==0 ? 1 : -seed); // Be sure to prevent seed=0.
	   seed2=seed;
	   for(j=NTAB+7;j>=0;j--) // load the shuffle table(after 8 warm-ups)
	   {
	     k=seed/IQ1;
	     seed=IA1*(seed-k*IQ1)-k*IR1;
	     if(seed<0)
	      seed+=IM1;
	     if(j<NTA
	      iv[j]=seed;
	   }
	   iy=iv[0];
	}
	k=seed/IQ1; // Start here when not initializing
	seed=IA1*(seed-k*IQ1)-k*IR1; // compute seed=(IA1*seed) % IM1 without overflows by 
	if(seed<0) // Schrage's method
	   seed+=IM1;
	k=seed2/IQ2;
	seed2=IA2*(seed2-k*IQ2)-k*IR2; // Compute seed2=(IA2*seed) % IM2 likewise.
	if(seed2<0)
	   seed2+=IM2;
	j=iy/NDIV; // Will be in the range 0..NTAB-1
	iy=iv[j]-seed2; //Here seed is shuffled, seed and seeds are combined 
	iv[j]=seed; //to generate out put
	if(iy<1)
	   iy+=IMM1;
	if((temp=AM*iy)>RNMX) // Because users don't expect endpoint values.
	   return RNMX;
	else 
	   return temp;
}

/* Park and Miller with Bays-Durham shuffle and added safeguards. Returns a uniform 
* random deviate between 0.0 and 1.0(exclusive of the endpoint 
* values). Call with seed a negative integer to initialize; thereafter, do not alter seed
between successive deviates in a sequence. RNMX should approaimate the largest floating
value that is less than 1.
this algorithm is ok when the required random numbers less than 
* 100,000,000. */
double ranPM(int &seed)
{
	const int IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32;
	const int NDIV=(1+(IM-1)/NTA;
	const double EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS);
	static int iy=0;
	static int iv[NTAB];
	int j,k;
	double temp;
	
	if(seed<=0 || !iy) // initialize, be sure to prevent seed=0.
	{
	   if(-seed<1) 
	     seed=1;
	   else 
	     seed=-seed;
	   for(j=NTAB+7;j>=0;j--)
	   {
	     k=seed/IQ;
	     seed=IA*(seed-k*IQ)-IR*k;
	     if(seed < 0)
	      seed += IM;
	     if(j<NTA
	      iv[j]=seed;
	   }
	   iy=iv[0];
	}
	k=seed/IQ; // start here when not initializing.
	seed=IA*(seed-k*IQ)-IR*k;// compute seed=(IA*seed)%IM without overflows by Schrage's method.
	j=iy/NDIV; // will be in the range 0..NTAB-1
	iy=iv[j]; // output previously stored value and refill the shuffle table
	iv[j]=seed;
	if((temp=AM*iy)>RNMX) // because users don't expect endpoint values.
	   return RNMX;
	else 
	   return temp;

}

double rnorm()
{
	static int iset=0;
	static double gset;
	double fac,rsq,v1,v2;
	
	if(seed<0) // reinitialize
	   iset=0;
	cout <<iset<<endl;
	if(iset==0) // we don't have an extra deviate handy, so pick two uniform numbers
	{ // in the square extending from -1 to +1 in each direction.
	   do
	   {
	     v1=2.0*ranPM(seed)-1.0;
	     v2=2.0*ranPM(seed)-1.0;
	     rsq=v1*v1+v2*v2; // see if they are in the unit circle, and if 
	   } while(rsq>=1.0||rsq==0.0);// they are not, try again.
	   fac=sqrt(-2.0*log(rsq)/rsq);
	// Now make the Box-Muller transformation to get  
	// two normal deviates. Return one and save the other for next time.
	   gset=v1*fac; // we have an extra deviate handy, 
	   iset=1; // so unset the flat, 
	   return v2*fac; // and return it. 
	}
	else 
	{
	   iset=0;
	   return gset;
	}
}

#endif /* _RANDOM_H */

#include <iostream>
#include "Random.h"
using namespace std;

int main()
{
	setSeed(-5);
	for (int i=0;i<100;i++)
	{ 
	cout<< rnorm()<<endl;
	return 0;
}