155 #if defined (HAVE_CONFIG_H)
166 #if ! defined (USE_X86_32)
167 # if defined (i386) || defined (HAVE_X86_32)
168 # define USE_X86_32 1
170 # define USE_X86_32 0
177 #define MATRIX_A 0x9908b0dfUL
178 #define UMASK 0x80000000UL
179 #define LMASK 0x7fffffffUL
180 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
181 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
195 state[0] = s & 0xffffffffUL;
196 for (j = 1; j <
MT_N; j++)
203 state[j] &= 0xffffffffUL;
219 k = (
MT_N > key_length ?
MT_N : key_length);
235 for (k =
MT_N - 1;
k; k--)
248 state[0] = 0x80000000UL;
256 uint32_t entropy[
MT_N];
260 FILE* urandom = std::fopen (
"/dev/urandom",
"rb");
265 unsigned char word[4];
266 if (std::fread (word, 4, 1, urandom) != 1)
268 entropy[n++] = word[0] + (word[1]<<8) + (word[2]<<16)
269 + (
static_cast<uint32_t
> (word[3])<<24);
282 entropy[n++] = clock ();
285 entropy[n++] = now.
usec ();
295 for (i = 0; i <
MT_N; i++)
305 for (i = 0; i <
MT_N; i++)
329 for (j =
MT_M; --j; p++)
347 y ^= (y << 7) & 0x9d2c5680UL;
348 y ^= (y << 15) & 0xefc60000UL;
349 return (y ^ (y >> 18));
355 #define randi32 randmt
360 const uint32_t lo =
randi32 ();
361 const uint32_t hi =
randi32 () & 0x1FFFFF;
362 #if defined (HAVE_X86_32)
364 uint32_t *
p = (uint32_t *)&u;
369 return ((static_cast<uint64_t> (hi) << 32) | lo);
376 const uint32_t lo =
randi32 ();
377 const uint32_t hi =
randi32 () & 0x3FFFFF;
378 #if defined (HAVE_X86_32)
380 uint32_t *
p =
static_cast<uint32_t *
> (&
u);
385 return ((static_cast<uint64_t> (hi) << 32) | lo);
393 return (static_cast<float> (
randi32 ()) + 0.5) * (1.0/4294967296.0);
403 return (a*67108864.0+b+0.4) * (1.0/9007199254740992.0);
422 #define ZIGGURAT_TABLE_SIZE 256
424 #define ZIGGURAT_NOR_R 3.6541528853610088
425 #define ZIGGURAT_NOR_INV_R 0.27366123732975828
426 #define NOR_SECTION_AREA 0.00492867323399
428 #define ZIGGURAT_EXP_R 7.69711747013104972
429 #define ZIGGURAT_EXP_INV_R 0.129918765548341586
430 #define EXP_SECTION_AREA 0.0039496598225815571993
432 #define ZIGINT uint64_t
433 #define EMANTISSA 9007199254740992.0
434 #define ERANDI randi53()
435 #define NMANTISSA EMANTISSA
436 #define NRANDI randi54()
437 #define RANDU randu53()
491 fi[255] = exp (-0.5 * x1 * x1);
502 for (i = 254; i > 0; i--)
510 fi[
i] = exp (-0.5 * x * x);
530 for (i = 254; i > 0; i--)
577 # if defined (HAVE_X86_32)
583 uint32_t *
p = (uint32_t *)&rabs;
589 p[1] = hi & 0x1FFFFF;
590 x = ( si ? -rabs : rabs ) *
wi[idx];
593 const uint64_t r =
NRANDI;
594 const int64_t rabs = r >> 1;
595 const int idx =
static_cast<int> (rabs & 0xFF);
596 const double x = ( (r & 1) ? -rabs : rabs) *
wi[idx];
598 if (rabs < static_cast<int64_t> (
ki[idx]))
618 while ( yy+yy <= xx*xx);
621 else if ((
fi[idx-1] -
fi[idx]) *
RANDU +
fi[idx] < exp (-0.5*x*x))
635 const int idx =
static_cast<int> (ri & 0xFF);
636 const double x = ri *
we[idx];
648 else if ((
fe[idx-1] -
fe[idx]) *
RANDU +
fe[idx] < exp (-x))
660 #define ZIGINT uint32_t
661 #define EMANTISSA 4294967296.0
662 #define ERANDI randi32()
663 #define NMANTISSA 2147483648.0
664 #define NRANDI randi32()
665 #define RANDU randu32()
681 ffi[255] = exp (-0.5 * x1 * x1);
692 for (i = 254; i > 0; i--)
700 ffi[
i] = exp (-0.5 * x * x);
709 ffe[255] = exp (-x1);
720 for (i = 254; i > 0; i--)
761 const uint32_t rabs = r &
LMASK;
762 const int idx =
static_cast<int> (r & 0xFF);
763 const float x =
static_cast<int32_t
> (r) *
fwi[idx];
784 while ( yy+yy <= xx*xx);
787 else if ((
ffi[idx-1] -
ffi[idx]) *
RANDU +
ffi[idx] < exp (-0.5*x*x))
801 const int idx =
static_cast<int> (ri & 0xFF);
802 const float x = ri *
fwe[idx];
824 for (i = 0; i < n; i++)
832 for (i = 0; i < n; i++)
840 for (i = 0; i < n; i++)
848 for (i = 0; i < n; i++)
856 for (i = 0; i < n; i++)
864 for (i = 0; i < n; i++)
static void next_state(void)
time_t unix_time(void) const
static void create_ziggurat_float_tables(void)
float oct_float_randu(void)
void oct_init_by_int(uint32_t s)
void oct_get_state(uint32_t *save)
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the base of natural logarithms The constant ex $e satisfies the equation log(e)
void oct_fill_float_randu(octave_idx_type n, float *p)
static uint64_t randi53(void)
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
void oct_fill_randu(octave_idx_type n, double *p)
float oct_float_randn(void)
static double randu53(void)
#define ZIGGURAT_TABLE_SIZE
static float randu32(void)
void oct_fill_randn(octave_idx_type n, double *p)
static uint32_t randmt(void)
static uint32_t state[624]
void oct_set_state(uint32_t *save)
void oct_fill_float_randn(octave_idx_type n, float *p)
static uint64_t randi54(void)
#define ZIGGURAT_NOR_INV_R
=val(i)}if ode{val(i)}occurs in table i
float oct_float_rande(void)
void oct_fill_rande(octave_idx_type n, double *p)
the element is set to zero In other the statement xample y
void oct_init_by_array(uint32_t *init_key, int key_length)
void oct_init_by_entropy(void)
static void create_ziggurat_tables(void)
void oct_fill_float_rande(octave_idx_type n, float *p)
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x