GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
randmtzig.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2006-2017 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 /*
24  A C-program for MT19937, with initialization improved 2002/2/10.
25  Coded by Takuji Nishimura and Makoto Matsumoto.
26  This is a faster version by taking Shawn Cokus's optimization,
27  Matthe Bellew's simplification, Isaku Wada's real version.
28  David Bateman added normal and exponential distributions following
29  Marsaglia and Tang's Ziggurat algorithm.
30 
31  Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
32  Copyright (C) 2004, David Bateman
33  All rights reserved.
34 
35  Redistribution and use in source and binary forms, with or without
36  modification, are permitted provided that the following conditions
37  are met:
38 
39  1. Redistributions of source code must retain the above copyright
40  notice, this list of conditions and the following disclaimer.
41 
42  2. Redistributions in binary form must reproduce the above copyright
43  notice, this list of conditions and the following disclaimer in the
44  documentation and/or other materials provided with the distribution.
45 
46  3. The names of its contributors may not be used to endorse or promote
47  products derived from this software without specific prior written
48  permission.
49 
50  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
51  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
52  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
53  A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
54  OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
55  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
56  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
57  PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
58  LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
59  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
60  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
61 
62 
63  Any feedback is very welcome.
64  http://www.math.keio.ac.jp/matumoto/emt.html
65  email: matumoto@math.keio.ac.jp
66 
67  * 2004-01-19 Paul Kienzle
68  * * comment out main
69  * add init_by_entropy, get_state, set_state
70  * * converted to allow compiling by C++ compiler
71  *
72  * 2004-01-25 David Bateman
73  * * Add Marsaglia and Tsang Ziggurat code
74  *
75  * 2004-07-13 Paul Kienzle
76  * * make into an independent library with some docs.
77  * * introduce new main and test code.
78  *
79  * 2004-07-28 Paul Kienzle & David Bateman
80  * * add -DALLBITS flag for 32 vs. 53 bits of randomness in mantissa
81  * * make the naming scheme more uniform
82  * * add -DHAVE_X86 for faster support of 53 bit mantissa on x86 arch.
83  *
84  * 2005-02-23 Paul Kienzle
85  * * fix -DHAVE_X86_32 flag and add -DUSE_X86_32=0|1 for explicit control
86  *
87  * 2006-04-01 David Bateman
88  * * convert for use in octave, declaring static functions only used
89  * here and adding oct_ to functions visible externally
90  * * inverse sense of ALLBITS
91  *
92  * 2012-05-18 David Bateman
93  * * Remove randu64 and ALLBIT option
94  * * Add the single precision generators
95  */
96 
97 /*
98  === Build instructions ===
99 
100  Compile with -DHAVE_GETTIMEOFDAY if the gettimeofday function is
101  available. This is not necessary if your architecture has
102  /dev/urandom defined.
103 
104  Uses implicit -Di386 or explicit -DHAVE_X86_32 to determine if CPU=x86.
105  You can force X86 behavior with -DUSE_X86_32=1, or suppress it with
106  -DUSE_X86_32=0. You should also consider -march=i686 or similar for
107  extra performance. Check whether -DUSE_X86_32=0 is faster on 64-bit
108  x86 architectures.
109 
110  If you want to replace the Mersenne Twister with another
111  generator then redefine randi32 appropriately.
112 
113  === Usage instructions ===
114  Before using any of the generators, initialize the state with one of
115  oct_init_by_int, oct_init_by_array or oct_init_by_entropy.
116 
117  All generators share the same state vector.
118 
119  === Mersenne Twister ===
120  void oct_init_by_int (uint32_t s) 32-bit initial state
121  void oct_init_by_array (uint32_t k[],int m) m*32-bit initial state
122  void oct_init_by_entropy (void) random initial state
123  void oct_get_state (uint32_t save[MT_N+1]) saves state in array
124  void oct_set_state (uint32_t save[MT_N+1]) restores state from array
125  static uint32_t randmt (void) returns 32-bit unsigned int
126 
127  === inline generators ===
128  static uint32_t randi32 (void) returns 32-bit unsigned int
129  static uint64_t randi53 (void) returns 53-bit unsigned int
130  static uint64_t randi54 (void) returns 54-bit unsigned int
131  static float randu32 (void) returns 32-bit uniform in (0,1)
132  static double randu53 (void) returns 53-bit uniform in (0,1)
133 
134  double oct_randu (void) returns M-bit uniform in (0,1)
135  double oct_randn (void) returns M-bit standard normal
136  double oct_rande (void) returns N-bit standard exponential
137 
138  float oct_float_randu (void) returns M-bit uniform in (0,1)
139  float oct_float_randn (void) returns M-bit standard normal
140  float oct_float_rande (void) returns N-bit standard exponential
141 
142  === Array generators ===
143  void oct_fill_randi32 (octave_idx_type, uint32_t [])
144  void oct_fill_randi64 (octave_idx_type, uint64_t [])
145 
146  void oct_fill_randu (octave_idx_type, double [])
147  void oct_fill_randn (octave_idx_type, double [])
148  void oct_fill_rande (octave_idx_type, double [])
149 
150  void oct_fill_float_randu (octave_idx_type, float [])
151  void oct_fill_float_randn (octave_idx_type, float [])
152  void oct_fill_float_rande (octave_idx_type, float [])
153 */
154 
155 #if defined (HAVE_CONFIG_H)
156 # include "config.h"
157 #endif
158 
159 #include <cstdio>
160 
161 #include "lo-math.h"
162 #include "oct-time.h"
163 #include "randmtzig.h"
164 
165 /* FIXME: may want to suppress X86 if sizeof(long) > 4 */
166 #if ! defined (USE_X86_32)
167 # if defined (i386) || defined (HAVE_X86_32)
168 # define USE_X86_32 1
169 # else
170 # define USE_X86_32 0
171 # endif
172 #endif
173 
174 /* ===== Mersenne Twister 32-bit generator ===== */
175 
176 #define MT_M 397
177 #define MATRIX_A 0x9908b0dfUL /* constant vector a */
178 #define UMASK 0x80000000UL /* most significant w-r bits */
179 #define LMASK 0x7fffffffUL /* least significant r bits */
180 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
181 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
182 
183 static uint32_t *next;
184 static uint32_t state[MT_N]; /* the array for the state vector */
185 static int left = 1;
186 static int initf = 0;
187 static int initt = 1;
188 static int inittf = 1;
189 
190 /* initializes state[MT_N] with a seed */
191 void
192 oct_init_by_int (uint32_t s)
193 {
194  int j;
195  state[0] = s & 0xffffffffUL;
196  for (j = 1; j < MT_N; j++)
197  {
198  state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j);
199  /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
200  /* In the previous versions, MSBs of the seed affect */
201  /* only MSBs of the array state[]. */
202  /* 2002/01/09 modified by Makoto Matsumoto */
203  state[j] &= 0xffffffffUL; /* for >32 bit machines */
204  }
205  left = 1;
206  initf = 1;
207 }
208 
209 /* initialize by an array with array-length */
210 /* init_key is the array for initializing keys */
211 /* key_length is its length */
212 void
213 oct_init_by_array (uint32_t *init_key, int key_length)
214 {
215  int i, j, k;
216  oct_init_by_int (19650218UL);
217  i = 1;
218  j = 0;
219  k = (MT_N > key_length ? MT_N : key_length);
220  for (; k; k--)
221  {
222  state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL))
223  + init_key[j] + j; /* non linear */
224  state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
225  i++;
226  j++;
227  if (i >= MT_N)
228  {
229  state[0] = state[MT_N-1];
230  i = 1;
231  }
232  if (j >= key_length)
233  j = 0;
234  }
235  for (k = MT_N - 1; k; k--)
236  {
237  state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL))
238  - i; /* non linear */
239  state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
240  i++;
241  if (i >= MT_N)
242  {
243  state[0] = state[MT_N-1];
244  i = 1;
245  }
246  }
247 
248  state[0] = 0x80000000UL; /* MSB is 1; assuring nonzero initial array */
249  left = 1;
250  initf = 1;
251 }
252 
253 void
255 {
256  uint32_t entropy[MT_N];
257  int n = 0;
258 
259  /* Look for entropy in /dev/urandom */
260  FILE* urandom = std::fopen ("/dev/urandom", "rb");
261  if (urandom)
262  {
263  while (n < MT_N)
264  {
265  unsigned char word[4];
266  if (std::fread (word, 4, 1, urandom) != 1)
267  break;
268  entropy[n++] = word[0] + (word[1]<<8) + (word[2]<<16)
269  + (static_cast<uint32_t> (word[3])<<24);
270  }
271  std::fclose (urandom);
272  }
273 
274  /* If there isn't enough entropy, gather some from various sources */
275 
277 
278  if (n < MT_N)
279  entropy[n++] = now.unix_time (); /* Current time in seconds */
280 
281  if (n < MT_N)
282  entropy[n++] = clock (); /* CPU time used (usec) */
283 
284  if (n < MT_N)
285  entropy[n++] = now.usec (); /* Fractional part of current time */
286 
287  /* Send all the entropy into the initial state vector */
288  oct_init_by_array (entropy,n);
289 }
290 
291 void
292 oct_set_state (uint32_t *save)
293 {
294  int i;
295  for (i = 0; i < MT_N; i++)
296  state[i] = save[i];
297  left = save[MT_N];
298  next = state + (MT_N - left + 1);
299 }
300 
301 void
302 oct_get_state (uint32_t *save)
303 {
304  int i;
305  for (i = 0; i < MT_N; i++)
306  save[i] = state[i];
307  save[MT_N] = left;
308 }
309 
310 static void
312 {
313  uint32_t *p = state;
314  int j;
315 
316  /* if init_by_int() has not been called, */
317  /* a default initial seed is used */
318  /* if (initf==0) init_by_int(5489UL); */
319  /* Or better yet, a random seed! */
320  if (initf == 0)
322 
323  left = MT_N;
324  next = state;
325 
326  for (j = MT_N - MT_M + 1; --j; p++)
327  *p = p[MT_M] ^ TWIST(p[0], p[1]);
328 
329  for (j = MT_M; --j; p++)
330  *p = p[MT_M-MT_N] ^ TWIST(p[0], p[1]);
331 
332  *p = p[MT_M-MT_N] ^ TWIST(p[0], state[0]);
333 }
334 
335 /* generates a random number on [0,0xffffffff]-interval */
336 static uint32_t
337 randmt (void)
338 {
339  uint32_t y;
340 
341  if (--left == 0)
342  next_state ();
343  y = *next++;
344 
345  /* Tempering */
346  y ^= (y >> 11);
347  y ^= (y << 7) & 0x9d2c5680UL;
348  y ^= (y << 15) & 0xefc60000UL;
349  return (y ^ (y >> 18));
350 }
351 
352 /* ===== Uniform generators ===== */
353 
354 /* Select which 32 bit generator to use */
355 #define randi32 randmt
356 
357 static uint64_t
358 randi53 (void)
359 {
360  const uint32_t lo = randi32 ();
361  const uint32_t hi = randi32 () & 0x1FFFFF;
362 #if defined (HAVE_X86_32)
363  uint64_t u;
364  uint32_t *p = (uint32_t *)&u;
365  p[0] = lo;
366  p[1] = hi;
367  return u;
368 #else
369  return ((static_cast<uint64_t> (hi) << 32) | lo);
370 #endif
371 }
372 
373 static uint64_t
374 randi54 (void)
375 {
376  const uint32_t lo = randi32 ();
377  const uint32_t hi = randi32 () & 0x3FFFFF;
378 #if defined (HAVE_X86_32)
379  uint64_t u;
380  uint32_t *p = static_cast<uint32_t *> (&u);
381  p[0] = lo;
382  p[1] = hi;
383  return u;
384 #else
385  return ((static_cast<uint64_t> (hi) << 32) | lo);
386 #endif
387 }
388 
389 /* generates a random number on (0,1)-real-interval */
390 static float
391 randu32 (void)
392 {
393  return (static_cast<float> (randi32 ()) + 0.5) * (1.0/4294967296.0);
394  /* divided by 2^32 */
395 }
396 
397 /* generates a random number on (0,1) with 53-bit resolution */
398 static double
399 randu53 (void)
400 {
401  const uint32_t a = randi32 () >> 5;
402  const uint32_t b = randi32 () >> 6;
403  return (a*67108864.0+b+0.4) * (1.0/9007199254740992.0);
404 }
405 
406 /* Determine mantissa for uniform doubles */
407 double
408 oct_randu (void)
409 {
410  return randu53 ();
411 }
412 
413 /* Determine mantissa for uniform floats */
414 float
416 {
417  return randu32 ();
418 }
419 
420 /* ===== Ziggurat normal and exponential generators ===== */
421 
422 #define ZIGGURAT_TABLE_SIZE 256
423 
424 #define ZIGGURAT_NOR_R 3.6541528853610088
425 #define ZIGGURAT_NOR_INV_R 0.27366123732975828
426 #define NOR_SECTION_AREA 0.00492867323399
427 
428 #define ZIGGURAT_EXP_R 7.69711747013104972
429 #define ZIGGURAT_EXP_INV_R 0.129918765548341586
430 #define EXP_SECTION_AREA 0.0039496598225815571993
431 
432 #define ZIGINT uint64_t
433 #define EMANTISSA 9007199254740992.0 /* 53 bit mantissa */
434 #define ERANDI randi53() /* 53 bits for mantissa */
435 #define NMANTISSA EMANTISSA
436 #define NRANDI randi54() /* 53 bits for mantissa + 1 bit sign */
437 #define RANDU randu53()
438 
443 
444 /*
445 This code is based on the paper Marsaglia and Tsang, "The ziggurat method
446 for generating random variables", Journ. Statistical Software. Code was
447 presented in this paper for a Ziggurat of 127 levels and using a 32 bit
448 integer random number generator. This version of the code, uses the
449 Mersenne Twister as the integer generator and uses 256 levels in the
450 Ziggurat. This has several advantages.
451 
452  1) As Marsaglia and Tsang themselves states, the more levels the few
453  times the expensive tail algorithm must be called
454  2) The cycle time of the generator is determined by the integer
455  generator, thus the use of a Mersenne Twister for the core random
456  generator makes this cycle extremely long.
457  3) The license on the original code was unclear, thus rewriting the code
458  from the article means we are free of copyright issues.
459  4) Compile flag for full 53-bit random mantissa.
460 
461 It should be stated that the authors made my life easier, by the fact that
462 the algorithm developed in the text of the article is for a 256 level
463 ziggurat, even if the code itself isn't...
464 
465 One modification to the algorithm developed in the article, is that it is
466 assumed that 0 <= x < Inf, and "unsigned long"s are used, thus resulting in
467 terms like 2^32 in the code. As the normal distribution is defined between
468 -Inf < x < Inf, we effectively only have 31 bit integers plus a sign. Thus
469 in Marsaglia and Tsang, terms like 2^32 become 2^31. We use NMANTISSA for
470 this term. The exponential distribution is one sided so we use the
471 full 32 bits. We use EMANTISSA for this term.
472 
473 It appears that I'm slightly slower than the code in the article, this
474 is partially due to a better generator of random integers than they
475 use. But might also be that the case of rapid return was optimized by
476 inlining the relevant code with a #define. As the basic Mersenne
477 Twister is only 25% faster than this code I suspect that the main
478 reason is just the use of the Mersenne Twister and not the inlining,
479 so I'm not going to try and optimize further.
480 */
481 
482 static void
484 {
485  int i;
486  double x, x1;
487 
488  /* Ziggurat tables for the normal distribution */
489  x1 = ZIGGURAT_NOR_R;
490  wi[255] = x1 / NMANTISSA;
491  fi[255] = exp (-0.5 * x1 * x1);
492 
493  /* Index zero is special for tail strip, where Marsaglia and Tsang
494  * defines this as
495  * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1,
496  * where v is the area of each strip of the ziggurat.
497  */
498  ki[0] = static_cast<ZIGINT> (x1 * fi[255] / NOR_SECTION_AREA * NMANTISSA);
499  wi[0] = NOR_SECTION_AREA / fi[255] / NMANTISSA;
500  fi[0] = 1.;
501 
502  for (i = 254; i > 0; i--)
503  {
504  /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
505  * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y))
506  */
507  x = sqrt (-2. * std::log (NOR_SECTION_AREA / x1 + fi[i+1]));
508  ki[i+1] = static_cast<ZIGINT> (x / x1 * NMANTISSA);
509  wi[i] = x / NMANTISSA;
510  fi[i] = exp (-0.5 * x * x);
511  x1 = x;
512  }
513 
514  ki[1] = 0;
515 
516  /* Zigurrat tables for the exponential distribution */
517  x1 = ZIGGURAT_EXP_R;
518  we[255] = x1 / EMANTISSA;
519  fe[255] = exp (-x1);
520 
521  /* Index zero is special for tail strip, where Marsaglia and Tsang
522  * defines this as
523  * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1,
524  * where v is the area of each strip of the ziggurat.
525  */
526  ke[0] = static_cast<ZIGINT> (x1 * fe[255] / EXP_SECTION_AREA * EMANTISSA);
527  we[0] = EXP_SECTION_AREA / fe[255] / EMANTISSA;
528  fe[0] = 1.;
529 
530  for (i = 254; i > 0; i--)
531  {
532  /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
533  * need inverse operator of y = exp(-x) -> x = -ln(y)
534  */
535  x = - std::log (EXP_SECTION_AREA / x1 + fe[i+1]);
536  ke[i+1] = static_cast<ZIGINT> (x / x1 * EMANTISSA);
537  we[i] = x / EMANTISSA;
538  fe[i] = exp (-x);
539  x1 = x;
540  }
541  ke[1] = 0;
542 
543  initt = 0;
544 }
545 
546 /*
547  * Here is the guts of the algorithm. As Marsaglia and Tsang state the
548  * algorithm in their paper
549  *
550  * 1) Calculate a random signed integer j and let i be the index
551  * provided by the rightmost 8-bits of j
552  * 2) Set x = j * w_i. If j < k_i return x
553  * 3) If i = 0, then return x from the tail
554  * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x
555  * 5) goto step 1
556  *
557  * Where f is the functional form of the distribution, which for a normal
558  * distribution is exp(-0.5*x*x)
559  */
560 
561 double
562 oct_randn (void)
563 {
564  if (initt)
566 
567  while (1)
568  {
569  /* The following code is specialized for 32-bit mantissa.
570  * Compared to the arbitrary mantissa code, there is a performance
571  * gain for 32-bits: PPC: 2%, MIPS: 8%, x86: 40%
572  * There is a bigger performance gain compared to using a full
573  * 53-bit mantissa: PPC: 60%, MIPS: 65%, x86: 240%
574  * Of course, different compilers and operating systems may
575  * have something to do with this.
576  */
577 # if defined (HAVE_X86_32)
578  /* 53-bit mantissa, 1-bit sign, x86 32-bit architecture */
579  double x;
580  int si,idx;
581  uint32_t lo, hi;
582  int64_t rabs;
583  uint32_t *p = (uint32_t *)&rabs;
584  lo = randi32 ();
585  idx = lo & 0xFF;
586  hi = randi32 ();
587  si = hi & UMASK;
588  p[0] = lo;
589  p[1] = hi & 0x1FFFFF;
590  x = ( si ? -rabs : rabs ) * wi[idx];
591 # else
592  /* arbitrary mantissa (selected by NRANDI, with 1 bit for sign) */
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];
597 # endif
598  if (rabs < static_cast<int64_t> (ki[idx]))
599  return x; /* 99.3% of the time we return here 1st try */
600  else if (idx == 0)
601  {
602  /* As stated in Marsaglia and Tsang
603  *
604  * For the normal tail, the method of Marsaglia[5] provides:
605  * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x,
606  * then return r+x. Except that r+x is always in the positive
607  * tail!!!! Any thing random might be used to determine the
608  * sign, but as we already have r we might as well use it
609  *
610  * [PAK] but not the bottom 8 bits, since they are all 0 here!
611  */
612  double xx, yy;
613  do
614  {
615  xx = - ZIGGURAT_NOR_INV_R * std::log (RANDU);
616  yy = - std::log (RANDU);
617  }
618  while ( yy+yy <= xx*xx);
619  return ((rabs & 0x100) ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx);
620  }
621  else if ((fi[idx-1] - fi[idx]) * RANDU + fi[idx] < exp (-0.5*x*x))
622  return x;
623  }
624 }
625 
626 double
627 oct_rande (void)
628 {
629  if (initt)
631 
632  while (1)
633  {
634  ZIGINT ri = ERANDI;
635  const int idx = static_cast<int> (ri & 0xFF);
636  const double x = ri * we[idx];
637  if (ri < ke[idx])
638  return x; /* 98.9% of the time we return here 1st try */
639  else if (idx == 0)
640  {
641  /* As stated in Marsaglia and Tsang
642  *
643  * For the exponential tail, the method of Marsaglia[5] provides:
644  * x = r - ln(U);
645  */
646  return ZIGGURAT_EXP_R - std::log (RANDU);
647  }
648  else if ((fe[idx-1] - fe[idx]) * RANDU + fe[idx] < exp (-x))
649  return x;
650  }
651 }
652 
653 #undef ZIGINT
654 #undef EMANTISSA
655 #undef ERANDI
656 #undef NMANTISSA
657 #undef NRANDI
658 #undef RANDU
659 
660 #define ZIGINT uint32_t
661 #define EMANTISSA 4294967296.0 /* 32 bit mantissa */
662 #define ERANDI randi32() /* 32 bits for mantissa */
663 #define NMANTISSA 2147483648.0 /* 31 bit mantissa */
664 #define NRANDI randi32() /* 31 bits for mantissa + 1 bit sign */
665 #define RANDU randu32()
666 
671 
672 static void
674 {
675  int i;
676  float x, x1;
677 
678  /* Ziggurat tables for the normal distribution */
679  x1 = ZIGGURAT_NOR_R;
680  fwi[255] = x1 / NMANTISSA;
681  ffi[255] = exp (-0.5 * x1 * x1);
682 
683  /* Index zero is special for tail strip, where Marsaglia and Tsang
684  * defines this as
685  * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1,
686  * where v is the area of each strip of the ziggurat.
687  */
688  fki[0] = static_cast<ZIGINT> (x1 * ffi[255] / NOR_SECTION_AREA * NMANTISSA);
689  fwi[0] = NOR_SECTION_AREA / ffi[255] / NMANTISSA;
690  ffi[0] = 1.;
691 
692  for (i = 254; i > 0; i--)
693  {
694  /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
695  * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y))
696  */
697  x = sqrt (-2. * std::log (NOR_SECTION_AREA / x1 + ffi[i+1]));
698  fki[i+1] = static_cast<ZIGINT> (x / x1 * NMANTISSA);
699  fwi[i] = x / NMANTISSA;
700  ffi[i] = exp (-0.5 * x * x);
701  x1 = x;
702  }
703 
704  fki[1] = 0;
705 
706  /* Zigurrat tables for the exponential distribution */
707  x1 = ZIGGURAT_EXP_R;
708  fwe[255] = x1 / EMANTISSA;
709  ffe[255] = exp (-x1);
710 
711  /* Index zero is special for tail strip, where Marsaglia and Tsang
712  * defines this as
713  * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1,
714  * where v is the area of each strip of the ziggurat.
715  */
716  fke[0] = static_cast<ZIGINT> (x1 * ffe[255] / EXP_SECTION_AREA * EMANTISSA);
717  fwe[0] = EXP_SECTION_AREA / ffe[255] / EMANTISSA;
718  ffe[0] = 1.;
719 
720  for (i = 254; i > 0; i--)
721  {
722  /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
723  * need inverse operator of y = exp(-x) -> x = -ln(y)
724  */
725  x = - std::log (EXP_SECTION_AREA / x1 + ffe[i+1]);
726  fke[i+1] = static_cast<ZIGINT> (x / x1 * EMANTISSA);
727  fwe[i] = x / EMANTISSA;
728  ffe[i] = exp (-x);
729  x1 = x;
730  }
731  fke[1] = 0;
732 
733  inittf = 0;
734 }
735 
736 /*
737  * Here is the guts of the algorithm. As Marsaglia and Tsang state the
738  * algorithm in their paper
739  *
740  * 1) Calculate a random signed integer j and let i be the index
741  * provided by the rightmost 8-bits of j
742  * 2) Set x = j * w_i. If j < k_i return x
743  * 3) If i = 0, then return x from the tail
744  * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x
745  * 5) goto step 1
746  *
747  * Where f is the functional form of the distribution, which for a normal
748  * distribution is exp(-0.5*x*x)
749  */
750 
751 float
753 {
754  if (inittf)
756 
757  while (1)
758  {
759  /* 32-bit mantissa */
760  const uint32_t r = randi32 ();
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];
764  if (rabs < fki[idx])
765  return x; /* 99.3% of the time we return here 1st try */
766  else if (idx == 0)
767  {
768  /* As stated in Marsaglia and Tsang
769  *
770  * For the normal tail, the method of Marsaglia[5] provides:
771  * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x,
772  * then return r+x. Except that r+x is always in the positive
773  * tail!!!! Any thing random might be used to determine the
774  * sign, but as we already have r we might as well use it
775  *
776  * [PAK] but not the bottom 8 bits, since they are all 0 here!
777  */
778  float xx, yy;
779  do
780  {
781  xx = - ZIGGURAT_NOR_INV_R * std::log (RANDU);
782  yy = - std::log (RANDU);
783  }
784  while ( yy+yy <= xx*xx);
785  return ((rabs & 0x100) ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx);
786  }
787  else if ((ffi[idx-1] - ffi[idx]) * RANDU + ffi[idx] < exp (-0.5*x*x))
788  return x;
789  }
790 }
791 
792 float
794 {
795  if (inittf)
797 
798  while (1)
799  {
800  ZIGINT ri = ERANDI;
801  const int idx = static_cast<int> (ri & 0xFF);
802  const float x = ri * fwe[idx];
803  if (ri < fke[idx])
804  return x; /* 98.9% of the time we return here 1st try */
805  else if (idx == 0)
806  {
807  /* As stated in Marsaglia and Tsang
808  *
809  * For the exponential tail, the method of Marsaglia[5] provides:
810  * x = r - ln(U);
811  */
812  return ZIGGURAT_EXP_R - std::log (RANDU);
813  }
814  else if ((ffe[idx-1] - ffe[idx]) * RANDU + ffe[idx] < exp (-x))
815  return x;
816  }
817 }
818 
819 /* Array generators */
820 void
822 {
824  for (i = 0; i < n; i++)
825  p[i] = oct_randu ();
826 }
827 
828 void
830 {
832  for (i = 0; i < n; i++)
833  p[i] = oct_randn ();
834 }
835 
836 void
838 {
840  for (i = 0; i < n; i++)
841  p[i] = oct_rande ();
842 }
843 
844 void
846 {
848  for (i = 0; i < n; i++)
849  p[i] = oct_float_randu ();
850 }
851 
852 void
854 {
856  for (i = 0; i < n; i++)
857  p[i] = oct_float_randn ();
858 }
859 
860 void
862 {
864  for (i = 0; i < n; i++)
865  p[i] = oct_float_rande ();
866 }
double oct_rande(void)
Definition: randmtzig.cc:627
double oct_randu(void)
Definition: randmtzig.cc:408
#define NMANTISSA
Definition: randmtzig.cc:663
static int left
Definition: randmtzig.cc:185
static void next_state(void)
Definition: randmtzig.cc:311
time_t unix_time(void) const
Definition: oct-time.h:95
#define UMASK
Definition: randmtzig.cc:178
#define ERANDI
Definition: randmtzig.cc:662
static void create_ziggurat_float_tables(void)
Definition: randmtzig.cc:673
float oct_float_randu(void)
Definition: randmtzig.cc:415
#define ZIGINT
Definition: randmtzig.cc:660
static int initf
Definition: randmtzig.cc:186
void oct_init_by_int(uint32_t s)
Definition: randmtzig.cc:192
void oct_get_state(uint32_t *save)
Definition: randmtzig.cc:302
fclose(in)
static float ffe[256]
Definition: randmtzig.cc:670
#define MT_N
Definition: randmtzig.h:69
octave::sys::time now
Definition: data.cc:6299
for large enough k
Definition: lu.cc:606
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)
static double fi[256]
Definition: randmtzig.cc:440
#define TWIST(u, v)
Definition: randmtzig.cc:181
u
Definition: lu.cc:138
#define NOR_SECTION_AREA
Definition: randmtzig.cc:426
s
Definition: file-io.cc:2682
static uint32_t fki[256]
Definition: randmtzig.cc:667
static uint64_t ke[256]
Definition: randmtzig.cc:441
void oct_fill_float_randu(octave_idx_type n, float *p)
Definition: randmtzig.cc:845
static uint64_t randi53(void)
Definition: randmtzig.cc:358
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
Definition: cellfun.cc:398
void oct_fill_randu(octave_idx_type n, double *p)
Definition: randmtzig.cc:821
#define randi32
Definition: randmtzig.cc:355
#define ZIGGURAT_EXP_R
Definition: randmtzig.cc:428
static uint64_t ki[256]
Definition: randmtzig.cc:439
float oct_float_randn(void)
Definition: randmtzig.cc:752
static double randu53(void)
Definition: randmtzig.cc:399
#define LMASK
Definition: randmtzig.cc:179
static uint32_t * next
Definition: randmtzig.cc:183
#define NRANDI
Definition: randmtzig.cc:664
static uint32_t fke[256]
Definition: randmtzig.cc:669
#define ZIGGURAT_TABLE_SIZE
Definition: randmtzig.cc:422
static float randu32(void)
Definition: randmtzig.cc:391
void oct_fill_randn(octave_idx_type n, double *p)
Definition: randmtzig.cc:829
static float ffi[256]
Definition: randmtzig.cc:668
static uint32_t randmt(void)
Definition: randmtzig.cc:337
#define RANDU
Definition: randmtzig.cc:665
static uint32_t state[624]
Definition: randmtzig.cc:184
#define MT_M
Definition: randmtzig.cc:176
void oct_set_state(uint32_t *save)
Definition: randmtzig.cc:292
void oct_fill_float_randn(octave_idx_type n, float *p)
Definition: randmtzig.cc:853
static float fwe[256]
Definition: randmtzig.cc:670
static uint64_t randi54(void)
Definition: randmtzig.cc:374
#define ZIGGURAT_NOR_INV_R
Definition: randmtzig.cc:425
static double wi[256]
Definition: randmtzig.cc:440
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
p
Definition: lu.cc:138
float oct_float_rande(void)
Definition: randmtzig.cc:793
static int initt
Definition: randmtzig.cc:187
void oct_fill_rande(octave_idx_type n, double *p)
Definition: randmtzig.cc:837
the element is set to zero In other the statement xample y
Definition: data.cc:5342
b
Definition: cellfun.cc:398
void oct_init_by_array(uint32_t *init_key, int key_length)
Definition: randmtzig.cc:213
void oct_init_by_entropy(void)
Definition: randmtzig.cc:254
static int inittf
Definition: randmtzig.cc:188
static void create_ziggurat_tables(void)
Definition: randmtzig.cc:483
static double fe[256]
Definition: randmtzig.cc:442
#define EMANTISSA
Definition: randmtzig.cc:661
#define EXP_SECTION_AREA
Definition: randmtzig.cc:430
static double we[256]
Definition: randmtzig.cc:442
double oct_randn(void)
Definition: randmtzig.cc:562
void oct_fill_float_rande(octave_idx_type n, float *p)
Definition: randmtzig.cc:861
#define ZIGGURAT_NOR_R
Definition: randmtzig.cc:424
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
long usec(void) const
Definition: oct-time.h:97
static float fwi[256]
Definition: randmtzig.cc:668