aboutsummaryrefslogtreecommitdiffstats
path: root/contrib/lua-torch/torch7/lib/TH/THRandom.c
diff options
context:
space:
mode:
Diffstat (limited to 'contrib/lua-torch/torch7/lib/TH/THRandom.c')
-rw-r--r--contrib/lua-torch/torch7/lib/TH/THRandom.c272
1 files changed, 272 insertions, 0 deletions
diff --git a/contrib/lua-torch/torch7/lib/TH/THRandom.c b/contrib/lua-torch/torch7/lib/TH/THRandom.c
new file mode 100644
index 000000000..86d721e7b
--- /dev/null
+++ b/contrib/lua-torch/torch7/lib/TH/THRandom.c
@@ -0,0 +1,272 @@
+#include "THGeneral.h"
+#include "THRandom.h"
+
+/* 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(void)
+{
+ 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(void)
+{
+ 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
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+
+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);
+}