00001 #ifndef _plRandom_h_
00002 #define _plRandom_h_
00003
00004
00005
00006
00007 #include <plNormalRandomTable.h>
00008
00009
00010 #define PL_USE_RAND48
00011
00012 #if defined(WIN32) || defined(_WIN32)
00013 #undef PL_USE_RAND48
00014 #endif
00015
00016 #define PL_USE_COMPILER_RND_GENARATOR
00017
00018 #ifndef PL_USE_COMPILER_RND_GENARATOR
00019 #include <random.hpp>
00020 static RUniform RU(1);
00021 static IUniform IU(1);
00022
00023
00024 inline unsigned int plRandomInt(unsigned int max)
00025 {
00026 return IU.number(0, max-1);
00027 }
00028
00029
00030 inline plFloat plRandomFloat(plFloat max)
00031 {
00032 return RU.number(0.0, max);
00033 }
00034
00035
00036 inline void plSRandom( int seed )
00037 {
00038 IU.seed(seed);
00039 RU.seed(seed);
00040 }
00041
00042 #else
00043
00044
00045 #ifdef PL_USE_RAND48
00046
00047
00048 inline unsigned int plRandomInt(unsigned int max)
00049 {
00050
00051 return (unsigned int)(max* drand48());
00052 }
00053
00054
00055 inline plFloat plRandomFloat(plFloat max)
00056 {
00057 return max*drand48();
00058 }
00059
00060
00061 inline void plSRandom( int seed )
00062 {
00063 srand48( seed );
00064 }
00065
00066
00067 #else
00068
00069 static const plFloat PL_RAND_MAX_PLUS_ONE = RAND_MAX+PL_ONE;
00070
00071
00072 inline unsigned int plRandomInt(unsigned int max)
00073 {
00074
00075 return (unsigned int)(max* (rand()/PL_RAND_MAX_PLUS_ONE));
00076 }
00077
00078
00079 inline plFloat plRandomFloat(plFloat max)
00080 {
00081 return max*rand()/PL_RAND_MAX_PLUS_ONE;
00082 }
00083
00084
00085 inline void plSRandom( int seed )
00086 {
00087 srand( seed );
00088 }
00089
00090 #endif // PL_USE_RAND48
00091
00092
00093 #endif // PL_USE_COMPILER_RND_GENARATOR
00094
00095 #ifndef PL_USE_TABLE_FOR_NORMAL_SAMPLING // do not use a table for normal random
00096 static plFloat plNormalRandom_z2;
00097 static bool plNormalRandom_use_z2 = false;
00098
00099
00100 inline plFloat plNormalRandom(plFloat m, plFloat sd)
00101 {
00102 plFloat v1, v2, w, z1;
00103
00104
00105
00106 if (plNormalRandom_use_z2) {
00107 plNormalRandom_use_z2 = false;
00108 return ( m + plNormalRandom_z2 * sd);
00109 }
00110
00111 do {
00112 v1 = PL_TWO * plRandomFloat(PL_ONE) - PL_ONE;
00113 v2 = PL_TWO * plRandomFloat(PL_ONE) - PL_ONE;
00114 w = v1 * v1 + v2 * v2;
00115 } while ( w >= PL_ONE);
00116
00117 w = sqrt( (- PL_TWO* log( w ) ) / w );
00118 z1 = v1 * w;
00119 plNormalRandom_z2 = v2 * w;
00120 plNormalRandom_use_z2 = true;
00121
00122 return( m + z1 * sd );
00123 }
00124
00125
00126 #else
00127
00128
00129
00130
00131 inline plFloat plNormalRandom(plFloat m, plFloat sd)
00132 {
00133 unsigned int ind = plRandomInt(pl_dim_gauss_table),rndsign = plRandomInt(2);
00134
00135 if(rndsign)
00136 return ((m + (sd * pl_gauss_table[ind])));
00137 else
00138 return ((m - (sd * pl_gauss_table[ind])));
00139 }
00140
00141 #endif
00142
00143
00144 inline plFloat plNormalRandom(plFloat m, plFloat sd, plFloat aLimit, plFloat bLimit)
00145 {
00146 plFloat result;
00147 do{
00148 result = plNormalRandom(m, sd);
00149 }
00150 while((result < aLimit) || (result >= bLimit));
00151 return result;
00152 }
00153
00154
00155 inline plFloat plNormalRandomUsingTable(plFloat m, plFloat sd,
00156 plFloat two_inta_minus_half,
00157 plFloat two_intb_minus_half)
00158 {
00159 plFloat u = two_inta_minus_half + plRandomFloat(two_intb_minus_half - two_inta_minus_half);
00160 int ind = int(u*pl_dim_gauss_table);
00161
00162 if(ind >= 0) return ((m + (sd * pl_gauss_table[ind])));
00163
00164 return ((m - (sd * pl_gauss_table[-ind])));
00165 }
00166
00167
00168 inline plFloat plHalfNormalRandom(plFloat m, plFloat sd)
00169 {
00170 if (sd > PL_ZERO)
00171 return (m + plFabs(plNormalRandom(PL_ZERO, sd)));
00172 else
00173 return (m - plFabs(plNormalRandom(PL_ZERO, -sd)));
00174 }
00175
00176
00177
00178
00179
00180
00181 inline unsigned long plDrawUsingRepartitionFunctionBetween(const vector <plProbValue> &repartition,
00182 unsigned long start_index,
00183 unsigned long end_index, const plProbValue& u)
00184 {
00185 if( end_index == start_index)
00186 return start_index;
00187
00188 unsigned long middle_index = (start_index + end_index) / 2;
00189
00190 if(repartition[middle_index] <= u)
00191 return plDrawUsingRepartitionFunctionBetween(repartition, middle_index+1, end_index, u);
00192
00193 return plDrawUsingRepartitionFunctionBetween(repartition, start_index, middle_index, u);
00194 }
00195
00196
00197 inline unsigned long plDrawUsingRepartitionFunction(const vector <plProbValue> &repartition )
00198 {
00199 plProbValue u = repartition[repartition.size() - 1] * plRandomFloat(PL_ONE);
00200 return plDrawUsingRepartitionFunctionBetween(repartition, 0, repartition.size() - 1, u);
00201 }
00202
00203
00204 inline unsigned long plDrawUsingRepartitionFunctionBetween(plProbValue *repartition,
00205 unsigned long start_index,
00206 unsigned long end_index, const plProbValue& u)
00207 {
00208 if( end_index == start_index)
00209 return start_index;
00210
00211 unsigned long middle_index = (start_index + end_index) / 2;
00212
00213 if(repartition[middle_index] <= u)
00214 return plDrawUsingRepartitionFunctionBetween(repartition, middle_index+1, end_index, u);
00215
00216 return plDrawUsingRepartitionFunctionBetween(repartition, start_index, middle_index, u);
00217 }
00218
00219
00220 inline unsigned long plDrawUsingRepartitionFunction( unsigned long size, plProbValue *repartition)
00221 {
00222 plProbValue u = repartition[size - 1] * plRandomFloat(PL_ONE);
00223 return plDrawUsingRepartitionFunctionBetween(repartition, 0, size - 1, u);
00224 }
00225
00226
00227
00228
00229 #endif