このコードは、[0,1) の範囲で 50 個の乱数を生成します。
int main(){
int i,j,M=50;
// If the interval is not uniform
double interval_widths[3] = { 0.1, 0.11, 0.03};
double interval_widths_sum[3];
//Divide into subintervals
void Init() {
interval_widths_sum[0] = 0;
for (int i=1; i<N; i++) {
interval_widths_sum[i] = interval_widths_sum[i-1] + interval_widths[i];
}
}
//check in which interval R is
int Seek(double R) {
int i;
if (R < 0.0) return -1;
for (i = 0; i < N; i++) {
if (R >= interval_widths_sum[i]) {
break;
}
return (i);
}
}
unsigned long init[4] = {0x123, 0x234, 0x345, 0x456}, length = 4;
MTRand drand;
for (i = 0; i < M; i++) {
double x=("%10.8f ", drand());
double y=("%10.8f ", drand());
double R=("%10.8f ", drand());
cout<<"(x,y)="<< x <<","<< y<< endl;
a[i]=x*y;
double *p= &x;
double *q= &y;
double *r= &R;
double z= Seek(*r);
cout<<"(x,y)="<<*p<<","<<*q<<endl;
cout<< Init<<" "<< z<<endl;
}
}
1) [0,1] の範囲で x と y の値を生成します。x= 0.11, 0.23, ..... と y= 0.13, 0.33, ..... などとしましょう。
2) ここで a= x*y を定義します。a1= 0.11*0.13、a2 = 0.23*0.33 などのように。
3) ここで、間隔 [0,1] をサイズ [0,a1],[a1,a1+a2],......,[ai,1] の (3j-2) サブ間隔に分割したいと思います。 (i=1...3j-3)。
4) 次に、[0,1] の範囲内で数値 R を生成します。(3j-2) のどれにこの R が含まれているかを調べます。