155 #if defined (HAVE_CONFIG_H)
162 #ifdef HAVE_GETTIMEOFDAY
163 #include <sys/time.h>
170 #if !defined (USE_X86_32)
171 # if defined (i386) || defined (HAVE_X86_32)
172 # define USE_X86_32 1
174 # define USE_X86_32 0
181 #define MATRIX_A 0x9908b0dfUL
182 #define UMASK 0x80000000UL
183 #define LMASK 0x7fffffffUL
184 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
185 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
199 state[0] = s & 0xffffffffUL;
200 for (j = 1; j <
MT_N; j++)
207 state[j] &= 0xffffffffUL;
223 k = (
MT_N > key_length ?
MT_N : key_length);
228 state[i] &= 0xffffffffUL;
239 for (k =
MT_N - 1; k; k--)
243 state[i] &= 0xffffffffUL;
252 state[0] = 0x80000000UL;
260 uint32_t entropy[
MT_N];
264 FILE* urandom = fopen (
"/dev/urandom",
"rb");
269 unsigned char word[4];
270 if (fread (word, 4, 1, urandom) != 1)
272 entropy[n++] = word[0]+(word[1]<<8)+(word[2]<<16)+((uint32_t)word[3]<<24);
279 entropy[n++] = time (NULL);
281 entropy[n++] = clock ();
282 #ifdef HAVE_GETTIMEOFDAY
286 if (gettimeofday (&tv, NULL) != -1)
287 entropy[n++] = tv.tv_usec;
298 for (i = 0; i <
MT_N; i++)
308 for (i = 0; i <
MT_N; i++)
332 for (j =
MT_M; --j; p++)
350 y ^= (y << 7) & 0x9d2c5680UL;
351 y ^= (y << 15) & 0xefc60000UL;
352 return (y ^ (y >> 18));
358 #define randi32 randmt
363 const uint32_t lo =
randi32 ();
364 const uint32_t hi =
randi32 ()&0x1FFFFF;
367 uint32_t *p = (uint32_t *)&u;
372 return (((uint64_t)hi<<32)|lo);
379 const uint32_t lo =
randi32 ();
380 const uint32_t hi =
randi32 ()&0x3FFFFF;
383 uint32_t *p = (uint32_t *)&u;
388 return (((uint64_t)hi<<32)|lo);
396 return ((
float)
randi32 () + 0.5) * (1.0/4294967296.0);
404 const uint32_t a =
randi32 ()>>5;
405 const uint32_t b =
randi32 ()>>6;
406 return (a*67108864.0+b+0.4) * (1.0/9007199254740992.0);
425 #define ZIGGURAT_TABLE_SIZE 256
427 #define ZIGGURAT_NOR_R 3.6541528853610088
428 #define ZIGGURAT_NOR_INV_R 0.27366123732975828
429 #define NOR_SECTION_AREA 0.00492867323399
431 #define ZIGGURAT_EXP_R 7.69711747013104972
432 #define ZIGGURAT_EXP_INV_R 0.129918765548341586
433 #define EXP_SECTION_AREA 0.0039496598225815571993
435 #define ZIGINT uint64_t
436 #define EMANTISSA 9007199254740992.0
437 #define ERANDI randi53()
438 #define NMANTISSA EMANTISSA
439 #define NRANDI randi54()
440 #define RANDU randu53()
494 fi[255] = exp (-0.5 * x1 * x1);
505 for (i = 254; i > 0; i--)
511 ki[i+1] = (
ZIGINT)(x / x1 * NMANTISSA);
513 fi[i] = exp (-0.5 * x * x);
533 for (i = 254; i > 0; i--)
539 ke[i+1] = (
ZIGINT)(x / x1 * EMANTISSA);
584 register uint32_t lo, hi;
586 uint32_t *p = (uint32_t *)&rabs;
593 x = ( si ? -rabs : rabs ) * wi[idx];
596 const uint64_t r =
NRANDI;
597 const int64_t rabs = r>>1;
598 const int idx = (
int)(rabs&0xFF);
599 const double x = ( r&1 ? -rabs : rabs) * wi[idx];
601 if (rabs < (int64_t)
ki[idx])
621 while ( yy+yy <= xx*xx);
624 else if ((fi[idx-1] - fi[idx]) *
RANDU + fi[idx] < exp (-0.5*x*x))
638 const int idx = (
int)(ri & 0xFF);
639 const double x = ri * we[idx];
651 else if ((fe[idx-1] - fe[idx]) *
RANDU + fe[idx] < exp (-x))
663 #define ZIGINT uint32_t
664 #define EMANTISSA 4294967296.0
665 #define ERANDI randi32()
666 #define NMANTISSA 2147483648.0
667 #define NRANDI randi32()
668 #define RANDU randu32()
685 ffi[255] = exp (-0.5 * x1 * x1);
696 for (i = 254; i > 0; i--)
704 ffi[i] = exp (-0.5 * x * x);
713 ffe[255] = exp (-x1);
724 for (i = 254; i > 0; i--)
765 const uint32_t rabs = r&
LMASK;
766 const int idx = (
int)(r&0xFF);
767 const float x = ((int32_t)r) * fwi[idx];
788 while ( yy+yy <= xx*xx);
791 else if ((ffi[idx-1] - ffi[idx]) *
RANDU + ffi[idx] < exp (-0.5*x*x))
805 const int idx = (
int)(ri & 0xFF);
806 const float x = ri * fwe[idx];
818 else if ((ffe[idx-1] - ffe[idx]) *
RANDU + ffe[idx] < exp (-x))
828 for (i = 0; i < n; i++)
836 for (i = 0; i < n; i++)
844 for (i = 0; i < n; i++)
852 for (i = 0; i < n; i++)
860 for (i = 0; i < n; i++)
868 for (i = 0; i < n; i++)
void oct_init_by_array(uint32_t *init_key, int key_length)
void oct_fill_rande(octave_idx_type n, double *p)
float oct_float_randu(void)
static uint32_t state[624]
void oct_fill_float_randn(octave_idx_type n, float *p)
void oct_init_by_int(uint32_t s)
void oct_fill_randn(octave_idx_type n, double *p)
void oct_fill_float_rande(octave_idx_type n, float *p)
float oct_float_rande(void)
static uint32_t randmt(void)
static float randu32(void)
void oct_init_by_entropy(void)
static double randu53(void)
void oct_fill_randu(octave_idx_type n, double *p)
static void create_ziggurat_float_tables(void)
static void next_state(void)
static uint64_t randi53(void)
void oct_fill_float_randu(octave_idx_type n, float *p)
float oct_float_randn(void)
static void create_ziggurat_tables(void)
static uint64_t randi54(void)
void oct_set_state(uint32_t *save)
#define ZIGGURAT_TABLE_SIZE
#define ZIGGURAT_NOR_INV_R
void oct_get_state(uint32_t *save)
F77_RET_T const double * x