Random Number 1
/*
1
/*
#include ＜stdio.h>
#include ＜math.h>
#include ＜stdlib.h>
unsigned long int next = 423758209;
double rand0_1(void);
main(void)
{
int i;
for(i=0;i<=10;i++)
{
printf("%d %lf\n", i,(double)rand0_1());
}
return 0;
}
double rand0_1(void)
{
next = next*19 % 423325498;
return((double)next/(double)423325498.);
}
void srand0_1(unsigned int seed)
{
next = seed;
}
実行結果
0 0.873641 1 0.453407 2 0.614738 3 0.534233 4 0.004652 5 0.088385 6 0.679320 7 0.761296 8 0.318850 9 0.058149 10 0.104833
Random Number 2
/*
NUMERICAL RECIPE IN C P204
一様乱数(線形合同法)
*/
#include ＜stdio.h>
#include ＜stdlib.h>
#include ＜math.h>
unsigned long next = 1;
int rand(void);
int main(void)
{
int i;
for(i=0;i<=10;i++)
{
printf("%d %ld\n", i,rand());
}
return 0;
}
int rand(void)
{
next = next*1103515245 + 12345;
return (unsigned int)(next/65536) % 32768;
}
void srand(unsigned int seed)
{
next = seed;
}
実行結果
0 16838 1 38526 2 10113 3 50283 4 63819 5 38395 6 55778 7 40187 8 48980 9 4086 10 2749
/*
一様乱数(線形合同法) Part II
x_n=(A x_{n-1} + C) mod M
Mが2^n, A が8で割って余り5の数、Cが奇数という条件で、
0-M-1までの整数が周期Mで1回ずつ現れる。
A=109,C=1021, M=32768, x_0=13
*/
#include ＜stdio.h>
unsigned rndnum=13; /*乱数の初期値*/
unsigned irnd();
main()
{
int j;
for (j=0;j<10;j++){
printf("%d\n", irnd());
}
}
unsigned irnd(void) /*0-32767 の整数乱数*/
{
rndnum=(rndnum*109+1021)%32768;
return rndnum;
}
実行結果
2438 4619 12972 5945 26434 31511 27848 21797 17598 18659
Random Number 3
/* Park and Miller の乱数生成ルーチン(numerical recippes p206) */
/* #define idum = (seed) としておいて、
main の方でも忘れずにfloat rand0_1(long *);と宣言する
本文のほうでは rand0_1(&idum) と使う。
*/
#include ＜stdio.h>
#include ＜math.h>
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define MASK 123459876
unsigned long int idum = 2048;
float ran0(long *);
int main(void)
{
int i;
for(i=0;i<=10;i++)
{
printf("%d %lf\n", i,ran0(&idum));
}
return 0;
}
float ran0(long *idum)
{
long k;
float ans;
*idum ^= MASK;
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
ans=AM*(*idum);
*idum ^= MASK;
return ans;
}
実行結果
idum = 1024 0 0.250592 1 0.708165 2 0.122989 3 0.072388 4 0.621420 5 0.207334 6 0.657995 7 0.922808 8 0.631353 9 0.156461 10 0.638768 idum = 2048 0 0.226550 1 0.624079 2 0.903439 3 0.099183 4 0.960467 5 0.564713 6 0.136934 7 0.450398 8 0.832795 9 0.779885 10 0.520386
Random Number 4
/*
NUMERICAL RECIPE IN C P207-8
Park and Millerの「最小」乱数生成法とBaye-Durhamの
切り混ぜに安全機構を付けたもの。0.0と1.0の間(両端を
除く)の一様乱数を返す。初期化にはidumを負の値にして
呼び出す。それ以降は呼び出しごとにidumを変えては
ならない。RNMXはほぼ1未満の最大の浮動小数点とする。
*/
#include ＜stdio.h>
#include ＜math.h>
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
unsigned long int idum =1024;
float ran1(long *);
int main(void)
{
int i;
for(i=0;i<=10;i++)
{
printf("%d %f\n", i,ran1(&idum));
}
return 0;
}
float ran1(long *idum)
{
int j;
long k;
static long iy=0;
static long iv[NTAB];
float temp;
if(*idum <= 0 || !iy){
if(-(*idum) < 1) *idum =1;
else *idum = -(*idum);
for (j=NTAB +7;j>=0; j--){
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
if( j < NTAB) iv[j]=*idum;
}
iy=iv[0];
}
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
j=iy/NDIV;
iy=iv[j];
iv[j]= *idum;
if((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}
実行結果
0 0.415999 1 0.091965 2 0.756410 3 0.529700 4 0.930436 5 0.383502 6 0.653919 7 0.066842 8 0.722660 9 0.671149 10 0.383416
Random Number 5
/*
NUMERICAL RECIPE IN C P210
0.0と1.0の間の一様乱数を返す。乱数列を初期化、再初期化
するにはidumを任意の負の値にする。
*/
#include ＜stdio.h>
#include ＜math.h>
#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)
unsigned long int idum =1024;
float ran3(long *);
int main(void)
{
int i;
for(i=0;i<=10;i++)
{
printf("%d %f\n", i,ran3(&idum));
}
return 0;
}
float ran3(long *idum)
{
static int inext, inextp;
static long ma[56];
static int iff=0;
long mj, mk;
int i,ii,k;
if(*idum < 0 || iff == 0){
iff=1;
mj=MSEED-(*idum < 0 ? -*idum : *idum);
mj %=MBIG;
ma[55]=mj;
mk=1;
for(i=1; i<=54; i++){
ii=(21*i) % 55;
ma[ii]=mk;
mk=mj-mk;
if(mk < MZ) mk += MBIG;
mj=ma[ii];
}
for(k=1; k<4;k++)
for(i=1; i<=55;i++){
ma[i] -= ma[1+(i+30) % 55];
if(ma[i] < MZ) ma[i] += MBIG;
}
inext=0;
inextp=31;
*idum=1;
}
if(++inext == 56) inext=1;
if(++inextp == 56) inext=1;
mj=ma[inext] - ma[inextp];
if(mj < MZ) mj += MBIG;
ma[inext]=mj;
return mj*FAC;
}
実行結果
0 0.461475 1 0.633433 2 0.846511 3 0.096041 4 0.353360 5 0.450735 6 0.639113 7 0.630816 8 0.662777 9 0.684881 10 0.103685
Normal Random Numbers (正規乱数)
Random Number 1
/*
Numerical Recipe in C, p.217
正規乱数
Box-Muller Method
平均0、分散1の正規分布の乱数。
*/
#include ＜stdio.h>
#include ＜math.h>
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
unsigned long int idum =1024;
float gasdev(long *);
float ran1(long *);
int main(void)
{
int i;
double x;
for(i=0;i<=20;i++)
{
x=gasdev(&idum);
printf("%d %f\n", i, x );
}
return 0;
}
float gasdev(long *idum)
{
float ran1(long *);
static int iset=0;
static float gset;
float fac,rsq,v1,v2;
if (iset==0){
do {
v1=2.0*ran1(idum)-1.0;
v2=2.0*ran1(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;
}
}
float ran1(long *idum)
{
int j;
long k;
static long iy=0;
static long iv[NTAB];
float temp;
if(*idum <= 0 || !iy){
if(-(*idum) < 1) *idum =1;
else *idum = -(*idum);
for (j=NTAB +7;j>=0; j--){
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
if( j < NTAB) iv[j]=*idum;
}
iy=iv[0];
}
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
j=iy/NDIV;
iy=iv[j];
iv[j]= *idum;
if((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}
実行結果
0 -0.836854 1 -0.172280 2 0.187118 3 1.615440 4 -0.176774 5 0.653145 6 -0.546364 7 0.194146 8 0.925709 9 1.204321 10 1.530555 11 -1.355560 12 0.051489 13 1.020178 14 -1.226161 15 0.708497 16 0.871673 17 -0.789721 18 0.332079 19 0.205602 20 -0.169367
Random Number 2
/*正規乱数(normal random numbers)
f(x)=\frac{1}{\sqrt{2\pi}\sigma}\exp{\frac{-(x-\mu)^2}{2\sigma^2}}
expected value: \mu, variance \sigma^2
*/
#include
#include
#include
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define MASK 123459876
double normrnd(void);
double norminv(double a);
unsigned long int idum = 2048;
float ran0(long *);
main()
{
int i,numrnd;
double a;
numrnd=10;
for(i=1;i<=numrnd;i++)
{
a=normrnd();
printf("%d %f\n",i, a);
}
return ;
}
double normrnd(void)
{
double a,b;
a=ran0(&idum);
b=norminv(a);
return b;
}
double norminv(double u)
{
double a,b,w,z;
z=-log(4.0*u*(1.0-u));
a=z*(2.0611786-5.7262204/(z+11.640595));
w=sqrt(a);
if(u<0.5)b=-w;
else b=w;
return b;
}
float ran0(long *idum)
{
long k;
float ans;
*idum ^= MASK;
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
ans=AM*(*idum);
*idum ^= MASK;
return ans;
}
実行結果
1 -0.750252 2 0.316090 3 1.301970 4 -1.286779 5 1.756900 6 0.162855 7 -1.094540 8 -0.124598 9 0.965464 10 0.771815
Binomial Random Numbers (2項乱数)
Random Number 1
/*
例)
表の出る確率がp, 裏の出る確率がq=1-pであるコインを投げて、表が出ると1点、
裏が出ると0点とする。このコインをn回投げる。このときの合計点をXとすると、
取りうる値はx=0,1,2,\cdots,nですが、各点に対する確率は,
f(x)=nCxp^xq^{n-x}=nCxp^x(1-p)^{n-x}
で与えられる。
expected value: \mu=p, variance : \sigma^2 = n p(1-p)
*/
#include
#include
#include
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define MASK 123459876
int birnd(int n,float p);
int bern(float p);
unsigned long int idum = 2048;
float ran0(long *);
main()
{
int i,k,n,numrnd;
float p;
numrnd=10;
n=5;
p=0.5;
for(i=1;i<=numrnd;i++)
{
k=birnd(n,p);
printf("%d %d\n",i, k);
}
return;
}
int birnd(int n,float p)
{
int i,k;
k=0;
for(i=1;i<=n;i++)k=k+bern(p);
return k;
}
int bern(float p)
{
int j;
double a;
a=ran0(&idum);
if(a＜p)j=1;
else j=0;
return j;
}
float ran0(long *idum)
{
long k;
float ans;
*idum ^= MASK;
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
ans=AM*(*idum);
*idum ^= MASK;
return ans;
}
実行結果
1 2 2 2 3 2 4 5 5 3 6 2 7 1 8 1 9 5 10 4
Poisson Random Numbers (ポアソン乱数)
Random Number 1
/*
一定量のウランのような半減期の長い元素があったとし、一定の観測時間内に
何個の原子が崩壊するかについての分布を考える。
個々の原子が観測時間内に崩壊する確率は非常に小さいのですが、
非常に多くの原子があるため、適当な長さの観測時間を取ればその時間内
にいくつかの原子の崩壊が記録される。
このように2項分布において、対象となるnが大きいが、起こる確率(生起確率)
pが小さく両方が釣り合って n p = \lambda を満足するケースを考える。
この場合、Xの確率分布を2項分布から求めることは非常に困難ですが、
n -> \infty, p -> 0 となった極限の分布は Poissonの小数の法則から
f(x)=e^{-\lambda} \lambda^x/x!
となる。
その他にも、事故の発生件数、不良品数、突然変異数などの事象に当てはまると言われている。
expected value=variance = \lambda
*/
#include
#include
#include
#include
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define MASK 123459876
int pornd(double a);
double exprnd(double a);
unsigned long int idum = 2048;
float ran0(long *);
main()
{
int i,numrnd;
double a;
srand((unsigned)time(NULL));
numrnd=10;
for(i=1;i<=numrnd;i++){
a=pornd(3);
printf("%d %f\n",i, a);
}
return;
}
int pornd(double t)
{
double t1;
int k;
k=0;
t1=exprnd(1);
while(t1＜t){
k=k+1;
t1=t1+exprnd(1);
}
return k;
}
double exprnd(double a)
{
double b,c;
b=ran0(&idum);
c=-(1/a)*log(b);
return c;
}
float ran0(long *idum)
{
long k;
float ans;
*idum ^= MASK;
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
ans=AM*(*idum);
*idum ^= MASK;
return ans;
}
実行結果
1 1.484790 2 0.471478 3 0.101547 4 2.310794 5 0.040336 6 0.571437 7 1.988256 8 0.797625 9 0.182968 10 0.248609
Exponential Random Numbers (指数乱数)
Random Number 1
/*
例)
放射性の原子があり、この原始が崩壊するまでの時間の分布
f(x)=a e^{-ax}
F(x)=1-e^{-ax}
expected value : \mu=\frac{1}{a},
variance : \sigma^2=\frac{1}{a^2}=\mu^2
*/
#include
#include
#include
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define MASK 123459876
double exprnd(double a);
unsigned long int idum = 2048;
float ran0(long *);
main()
{
int i,numrnd;
double a;
numrnd=10;
for(i=1;i<=numrnd;i++)
{
a=exprnd(1);
printf("%d %f\n",i, a);
}
return;
}
double exprnd(double a)
{
double b,c;
b=ran0(&idum);
c=-(1/a)*log(b);
return c;
}
float ran0(long *idum)
{
long k;
float ans;
*idum ^= MASK;
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
ans=AM*(*idum);
*idum ^= MASK;
return ans;
}
実行結果
1 1.484790 2 0.471478 3 0.101547 4 2.310794 5 0.040336 6 0.571437 7 1.988256 8 0.797625 9 0.182968 10 0.248609
Test the Random Numbers
/*
一様性のテスト
*/
#include ＜stdio.h>
#define N 1000 /*乱数の発生回数*/
#define M 10 /*整数乱数の範囲*/
#define F (N/M) /*期待値*/
#define SCALE (40.0/F) /*ヒストグラムの高さ(自動スケール)*/
unsigned rndnum=13;
unsigned irnd();
double rnd();
main()
{
int i,j,rank,hist[M+1];
double e=0.0;
for (i=1;i<=M;i++)
hist[i]=0;
for (i=0;i＜N;i++){
rank=(int)(M*rnd()+1); /*1-M の乱数を1つ発生*/
hist[rank]++;
}
for (i=1;i<=M;i++){
printf("%3d:%3d", i, hist[i]);
for (j=0;j＜hist[i]*SCALE;j++) /*ヒストグラムの表示*/
printf("*");
printf("\n");
e=e+(double)(hist[i]-F)*(hist[i]-F)/F; /*\chi^2の計算*/
}
printf("chi^2=%f\n",e);
}
unsigned irnd(void) /*0-32767 の整数乱数*/
{
rndnum=(rndnum*109+1021)%32768;
return rndnum;
}
double rnd(void) /*0-1未満の実数乱数*/
{
return irnd()/32767.1;
}
Back to C Language
