| #include "THGeneral.h" |
| #include "THRandom.h" |
| |
| #ifndef _WIN32 |
| #include <fcntl.h> |
| #include <unistd.h> |
| #endif |
| |
| /* Code for the Mersenne Twister random generator.... */ |
| #define n _MERSENNE_STATE_N |
| #define m _MERSENNE_STATE_M |
| |
| /* Creates (unseeded) new generator*/ |
| static THGenerator* THGenerator_newUnseeded() |
| { |
| THGenerator *self = THAlloc(sizeof(THGenerator)); |
| memset(self, 0, sizeof(THGenerator)); |
| self->left = 1; |
| self->seeded = 0; |
| self->normal_is_valid = 0; |
| return self; |
| } |
| |
| /* Creates new generator and makes sure it is seeded*/ |
| THGenerator* THGenerator_new() |
| { |
| THGenerator *self = THGenerator_newUnseeded(); |
| THRandom_seed(self); |
| return self; |
| } |
| |
| THGenerator* THGenerator_copy(THGenerator *self, THGenerator *from) |
| { |
| memcpy(self, from, sizeof(THGenerator)); |
| return self; |
| } |
| |
| void THGenerator_free(THGenerator *self) |
| { |
| THFree(self); |
| } |
| |
| int THGenerator_isValid(THGenerator *_generator) |
| { |
| if ((_generator->seeded == 1) && |
| (_generator->left > 0 && _generator->left <= n) && (_generator->next <= n)) |
| return 1; |
| |
| return 0; |
| } |
| |
| #ifndef _WIN32 |
| static unsigned long readURandomLong() |
| { |
| int randDev = open("/dev/urandom", O_RDONLY); |
| unsigned long randValue; |
| if (randDev < 0) { |
| THError("Unable to open /dev/urandom"); |
| } |
| ssize_t readBytes = read(randDev, &randValue, sizeof(randValue)); |
| if (readBytes < sizeof(randValue)) { |
| THError("Unable to read from /dev/urandom"); |
| } |
| close(randDev); |
| return randValue; |
| } |
| #endif // _WIN32 |
| |
| unsigned long THRandom_seed(THGenerator *_generator) |
| { |
| #ifdef _WIN32 |
| unsigned long s = (unsigned long)time(0); |
| #else |
| unsigned long s = readURandomLong(); |
| #endif |
| THRandom_manualSeed(_generator, s); |
| return s; |
| } |
| |
| /* The next 4 methods are taken from http:www.math.keio.ac.jpmatumotoemt.html |
| Here is the copyright: |
| Some minor modifications have been made to adapt to "my" C... */ |
| |
| /* |
| A C-program for MT19937, with initialization improved 2002/2/10. |
| Coded by Takuji Nishimura and Makoto Matsumoto. |
| This is a faster version by taking Shawn Cokus's optimization, |
| Matthe Bellew's simplification, Isaku Wada's double version. |
| |
| Before using, initialize the state by using init_genrand(seed) |
| or init_by_array(init_key, key_length). |
| |
| Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, |
| All rights reserved. |
| |
| Redistribution and use in source and binary forms, with or without |
| modification, are permitted provided that the following conditions |
| are met: |
| |
| 1. Redistributions of source code must retain the above copyright |
| notice, this list of conditions and the following disclaimer. |
| |
| 2. Redistributions in binary form must reproduce the above copyright |
| notice, this list of conditions and the following disclaimer in the |
| documentation and/or other materials provided with the distribution. |
| |
| 3. The names of its contributors may not be used to endorse or promote |
| products derived from this software without specific prior written |
| permission. |
| |
| THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
| "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
| LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
| A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR |
| CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
| EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
| PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
| PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
| LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
| NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
| SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
| |
| |
| Any feedback is very welcome. |
| http://www.math.keio.ac.jp/matumoto/emt.html |
| email: matumoto@math.keio.ac.jp |
| */ |
| |
| /* Macros for the Mersenne Twister random generator... */ |
| /* Period parameters */ |
| /* #define n 624 */ |
| /* #define m 397 */ |
| #define MATRIX_A 0x9908b0dfUL /* constant vector a */ |
| #define UMASK 0x80000000UL /* most significant w-r bits */ |
| #define LMASK 0x7fffffffUL /* least significant r bits */ |
| #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) |
| #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) |
| /*********************************************************** That's it. */ |
| |
| void THRandom_manualSeed(THGenerator *_generator, unsigned long the_seed_) |
| { |
| int j; |
| |
| /* This ensures reseeding resets all of the state (i.e. state for Gaussian numbers) */ |
| THGenerator *blank = THGenerator_newUnseeded(); |
| THGenerator_copy(_generator, blank); |
| THGenerator_free(blank); |
| |
| _generator->the_initial_seed = the_seed_; |
| _generator->state[0] = _generator->the_initial_seed & 0xffffffffUL; |
| for(j = 1; j < n; j++) |
| { |
| _generator->state[j] = (1812433253UL * (_generator->state[j-1] ^ (_generator->state[j-1] >> 30)) + j); |
| /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ |
| /* In the previous versions, mSBs of the seed affect */ |
| /* only mSBs of the array state[]. */ |
| /* 2002/01/09 modified by makoto matsumoto */ |
| _generator->state[j] &= 0xffffffffUL; /* for >32 bit machines */ |
| } |
| _generator->left = 1; |
| _generator->seeded = 1; |
| } |
| |
| unsigned long THRandom_initialSeed(THGenerator *_generator) |
| { |
| return _generator->the_initial_seed; |
| } |
| |
| void THRandom_nextState(THGenerator *_generator) |
| { |
| unsigned long *p = _generator->state; |
| int j; |
| |
| _generator->left = n; |
| _generator->next = 0; |
| |
| for(j = n-m+1; --j; p++) |
| *p = p[m] ^ TWIST(p[0], p[1]); |
| |
| for(j = m; --j; p++) |
| *p = p[m-n] ^ TWIST(p[0], p[1]); |
| |
| *p = p[m-n] ^ TWIST(p[0], _generator->state[0]); |
| } |
| |
| unsigned long THRandom_random(THGenerator *_generator) |
| { |
| unsigned long y; |
| |
| if (--(_generator->left) == 0) |
| THRandom_nextState(_generator); |
| y = *(_generator->state + (_generator->next)++); |
| |
| /* Tempering */ |
| y ^= (y >> 11); |
| y ^= (y << 7) & 0x9d2c5680UL; |
| y ^= (y << 15) & 0xefc60000UL; |
| y ^= (y >> 18); |
| |
| return y; |
| } |
| |
| /* generates a random number on [0,1)-double-interval */ |
| static double __uniform__(THGenerator *_generator) |
| { |
| /* divided by 2^32 */ |
| return (double)THRandom_random(_generator) * (1.0/4294967296.0); |
| } |
| |
| /********************************************************* |
| |
| Thanks *a lot* Takuji Nishimura and Makoto Matsumoto! |
| |
| Now my own code... |
| |
| *********************************************************/ |
| |
| double THRandom_uniform(THGenerator *_generator, double a, double b) |
| { |
| return(__uniform__(_generator) * (b - a) + a); |
| } |
| |
| double THRandom_normal(THGenerator *_generator, double mean, double stdv) |
| { |
| THArgCheck(stdv > 0, 2, "standard deviation must be strictly positive"); |
| |
| /* This is known as the Box-Muller method */ |
| if(!_generator->normal_is_valid) |
| { |
| _generator->normal_x = __uniform__(_generator); |
| _generator->normal_y = __uniform__(_generator); |
| _generator->normal_rho = sqrt(-2. * log(1.0-_generator->normal_y)); |
| _generator->normal_is_valid = 1; |
| } |
| else |
| _generator->normal_is_valid = 0; |
| |
| if(_generator->normal_is_valid) |
| return _generator->normal_rho*cos(2.*M_PI*_generator->normal_x)*stdv+mean; |
| else |
| return _generator->normal_rho*sin(2.*M_PI*_generator->normal_x)*stdv+mean; |
| } |
| |
| double THRandom_exponential(THGenerator *_generator, double lambda) |
| { |
| return(-1. / lambda * log(1-__uniform__(_generator))); |
| } |
| |
| double THRandom_cauchy(THGenerator *_generator, double median, double sigma) |
| { |
| return(median + sigma * tan(M_PI*(__uniform__(_generator)-0.5))); |
| } |
| |
| /* Faut etre malade pour utiliser ca. |
| M'enfin. */ |
| double THRandom_logNormal(THGenerator *_generator, double mean, double stdv) |
| { |
| THArgCheck(stdv > 0, 2, "standard deviation must be strictly positive"); |
| return(exp(THRandom_normal(_generator, mean, stdv))); |
| } |
| |
| int THRandom_geometric(THGenerator *_generator, double p) |
| { |
| THArgCheck(p > 0 && p < 1, 1, "must be > 0 and < 1"); |
| return((int)(log(1-__uniform__(_generator)) / log(p)) + 1); |
| } |
| |
| int THRandom_bernoulli(THGenerator *_generator, double p) |
| { |
| THArgCheck(p >= 0 && p <= 1, 1, "must be >= 0 and <= 1"); |
| return(__uniform__(_generator) <= p); |
| } |