152 #if defined (HAVE_CONFIG_H) 165 #if ! defined (USE_X86_32) 166 # if defined (i386) || defined (HAVE_X86_32) 167 # define USE_X86_32 1 169 # define USE_X86_32 0 176 #define MATRIX_A 0x9908b0dfUL 177 #define UMASK 0x80000000UL 178 #define LMASK 0x7fffffffUL 179 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) 180 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) 194 state[0] =
s & 0xffffffffUL;
195 for (j = 1; j <
MT_N; j++)
202 state[j] &= 0xffffffffUL;
218 k = (
MT_N > key_length ?
MT_N : key_length);
247 state[0] = 0x80000000UL;
255 uint32_t entropy[
MT_N];
259 FILE *urandom = std::fopen (
"/dev/urandom",
"rb");
264 unsigned char word[4];
265 if (std::fread (word, 4, 1, urandom) != 1)
267 entropy[n++] = word[0] + (word[1]<<8) + (word[2]<<16)
268 + (
static_cast<uint32_t
> (word[3])<<24);
281 entropy[n++] = clock ();
324 for (j =
MT_M; --j;
p++)
342 y ^= (
y << 7) & 0x9d2c5680UL;
343 y ^= (
y << 15) & 0xefc60000UL;
344 return (
y ^ (
y >> 18));
350 #define randi32 randmt 355 const uint32_t lo =
randi32 ();
356 const uint32_t hi =
randi32 () & 0x1FFFFF;
357 #if defined (HAVE_X86_32) 359 uint32_t *
p = (uint32_t *)&
u;
364 return ((static_cast<uint64_t> (hi) << 32) | lo);
371 const uint32_t lo =
randi32 ();
372 const uint32_t hi =
randi32 () & 0x3FFFFF;
373 #if defined (HAVE_X86_32) 375 uint32_t *
p =
static_cast<uint32_t *
> (&
u);
380 return ((static_cast<uint64_t> (hi) << 32) | lo);
388 return (static_cast<float> (
randi32 ()) + 0.5) * (1.0/4294967296.0);
398 return (
a*67108864.0+
b+0.4) * (1.0/9007199254740992.0);
417 #define ZIGGURAT_TABLE_SIZE 256 419 #define ZIGGURAT_NOR_R 3.6541528853610088 420 #define ZIGGURAT_NOR_INV_R 0.27366123732975828 421 #define NOR_SECTION_AREA 0.00492867323399 423 #define ZIGGURAT_EXP_R 7.69711747013104972 424 #define ZIGGURAT_EXP_INV_R 0.129918765548341586 425 #define EXP_SECTION_AREA 0.0039496598225815571993 427 #define ZIGINT uint64_t 428 #define EMANTISSA 9007199254740992.0 429 #define ERANDI randi53() 430 #define NMANTISSA EMANTISSA 431 #define NRANDI randi54() 432 #define RANDU randu53() 486 fi[255] = exp (-0.5 * x1 * x1);
497 for (
i = 254;
i > 0;
i--)
505 fi[
i] = exp (-0.5 *
x *
x);
525 for (
i = 254;
i > 0;
i--)
572 # if defined (HAVE_X86_32) 578 uint32_t *
p = (uint32_t *)&rabs;
584 p[1] = hi & 0x1FFFFF;
585 x = ( si ? -rabs : rabs ) *
wi[idx];
588 const uint64_t r =
NRANDI;
589 const int64_t rabs = r >> 1;
590 const int idx =
static_cast<int> (rabs & 0xFF);
591 const double x = ( (r & 1) ? -rabs : rabs) *
wi[idx];
593 if (rabs < static_cast<int64_t> (
ki[idx]))
613 while ( yy+yy <= xx*xx);
616 else if ((
fi[idx-1] -
fi[idx]) *
RANDU +
fi[idx] < exp (-0.5*
x*
x))
630 const int idx =
static_cast<int> (ri & 0xFF);
631 const double x = ri *
we[idx];
643 else if ((
fe[idx-1] -
fe[idx]) *
RANDU +
fe[idx] < exp (-
x))
655 #define ZIGINT uint32_t 656 #define EMANTISSA 4294967296.0 657 #define ERANDI randi32() 658 #define NMANTISSA 2147483648.0 659 #define NRANDI randi32() 660 #define RANDU randu32() 676 ffi[255] = exp (-0.5 * x1 * x1);
687 for (
i = 254;
i > 0;
i--)
695 ffi[
i] = exp (-0.5 *
x *
x);
704 ffe[255] = exp (-x1);
715 for (
i = 254;
i > 0;
i--)
756 const uint32_t rabs = r &
LMASK;
757 const int idx =
static_cast<int> (r & 0xFF);
758 const float x =
static_cast<int32_t
> (r) *
fwi[idx];
779 while ( yy+yy <= xx*xx);
796 const int idx =
static_cast<int> (ri & 0xFF);
797 const float x = ri *
fwe[idx];
static void next_state(void)
static void create_ziggurat_float_tables(void)
float oct_float_randu(void)
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)
void oct_set_state(const uint32_t *save)
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)
OCTAVE_EXPORT octave_value_list or both For fclose
static double randu53(void)
#define ZIGGURAT_TABLE_SIZE
static float randu32(void)
void oct_fill_randn(octave_idx_type n, double *p)
time_t unix_time(void) const
static uint32_t randmt(void)
static uint32_t state[624]
void oct_fill_float_randn(octave_idx_type n, float *p)
void oct_init_by_int(const uint32_t s)
static uint64_t randi54(void)
#define ZIGGURAT_NOR_INV_R
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_entropy(void)
static void create_ziggurat_tables(void)
void oct_fill_float_rande(octave_idx_type n, float *p)
void oct_init_by_array(const uint32_t *init_key, const int key_length)
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 const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE * x