00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include<stdio.h>
00031 #include <math.h>
00032 #include "Rnd_no_gen.h"
00033
00034 #define N 624
00035 #define M 397
00036 #define MATRIX_A 0x9908b0df
00037 #define UPPER_MASK 0x80000000
00038 #define LOWER_MASK 0x7fffffff
00039
00040 #define TEMPERING_MASK_B 0x9d2c5680
00041 #define TEMPERING_MASK_C 0xefc60000
00042 #define TEMPERING_SHIFT_U(y) (y >> 11)
00043 #define TEMPERING_SHIFT_S(y) (y << 7)
00044 #define TEMPERING_SHIFT_T(y) (y << 15)
00045 #define TEMPERING_SHIFT_L(y) (y >> 18)
00046
00047 static unsigned long mt[N];
00048 static int mti=N+1;
00049
00050 void SetGenRand(unsigned long seed)
00051 {
00052 mt[0]= seed & 0xffffffff;
00053 for (mti=1; mti<N; mti++)
00054 mt[mti] = (69069 * mt[mti-1]) & 0xffffffff;
00055 }
00056
00057 double GenRand()
00058 {
00059 unsigned long y;
00060 static unsigned long mag01[2]={0x0, MATRIX_A};
00061
00062 if (mti >= N) {
00063 int kk;
00064
00065 if (mti == N+1)
00066 SetGenRand(4357);
00067
00068 for (kk=0;kk<N-M;kk++) {
00069 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
00070 mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
00071 }
00072 for (;kk<N-1;kk++) {
00073 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
00074 mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
00075 }
00076 y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
00077 mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
00078
00079 mti = 0;
00080 }
00081
00082 y = mt[mti++];
00083 y ^= TEMPERING_SHIFT_U(y);
00084 y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
00085 y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
00086 y ^= TEMPERING_SHIFT_L(y);
00087
00088 return ( (double)y / (unsigned long)0xffffffff );
00089 }
00090
00091 double GenRand(double lo, double hi)
00092 {
00093 double res;
00094 res = lo + (hi - lo)*GenRand();
00095 return res;
00096 }
00097
00098 #define IA 16807
00099 #define IM 2147483647
00100 #define AM (1.0/IM)
00101 #define IQ 127773
00102 #define IR 2836
00103 #define NTAB
00104
00105 double GaussRand()
00106 {
00107 static int iset = 0;
00108 static double gset;
00109 double fac, rsq, v1, v2;
00110
00111 if (iset == 0) {
00112 do {
00113 v1 = 2.0*GenRand() - 1.0;
00114 v2 = 2.0*GenRand() - 1.0;
00115 rsq = v1*v1 +v2*v2;
00116 } while (rsq >= 1.0 || rsq == 0.0);
00117 fac = sqrt(-2.0*log(rsq)/rsq);
00118 gset = v1*fac;
00119 iset = 1;
00120 return v2*fac;
00121 } else {
00122 iset = 0;
00123 return gset;
00124 }
00125 }
00126
00127 double GaussRand(double mean, double stdev)
00128 {
00129 double res;
00130 res = mean + stdev*GaussRand();
00131 return res;
00132 }
00133
00134
00135
00136
00137
00138