Mercurial > octave
annotate liboctave/numeric/randmtzig.cc @ 24856:8bb0251fcfde
Use a uint32 state vector for random number generators (bug #50256).
* oct-rand.h: Remove #include "dColVector.h" and replace with #include
"uint32NDArray.h". Replace all instances of ColumnVector with uint32NDArray.
* oct-rand.cc: Replace all instances of ColumnVector with uint32NDArray.
* oct-rand.cc (get_internal_state): Remove use of OCTAVE_LOCAL_BUFFER.
Use s.fortran_vec () to pass pointer to state vector data directly to
oct_get_state().
* oct-rand.cc (set_internal_state): Remove use of OCTAVE_LOCAL_BUFFER.
Use s.data () to pass const pointer to state vector data directly to
oct_set_state().
* oct-rand.cc (double2uint32): Remove obsolete function.
* randmtzig.h (oct_init_by_int): Change prototype to accept const value.
* randmtzig.h (oct_init_by_array): Change prototype to accept const value.
* randmtzig.h (oct_set_state): Change prototype to accept const value.
* randmtzig.cc (oct_init_by_int, oct_init_by_array, oct_set_state): Change
functions to match new prototypes. Declare temporary loop variable
just-in-time.
* randmtzig.cc (oct_get_state, oct_set_state): Use std::copy_n instead of
for loop to simplify code.
* rand.cc (Frand, Frandn, Frande, Frandg, Frandp): Change tolerances for BIST
tests which use the new RNG to be much tighter. Use '_' in large numbers for
readability. Change out-of-range seed testing to match new behavior.
| author | Rik <rik@octave.org> |
|---|---|
| date | Fri, 09 Mar 2018 16:40:12 -0800 |
| parents | 194eb4bd202b |
| children | b8ce68627441 |
| rev | line source |
|---|---|
| 7019 | 1 /* |
| 2 | |
|
23219
3ac9f9ecfae5
maint: Update copyright dates.
John W. Eaton <jwe@octave.org>
parents:
23083
diff
changeset
|
3 Copyright (C) 2006-2017 John W. Eaton |
| 7019 | 4 |
| 5 This file is part of Octave. | |
| 6 | |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23662
diff
changeset
|
7 Octave is free software: you can redistribute it and/or modify it |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
8 under the terms of the GNU General Public License as published by |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23662
diff
changeset
|
9 the Free Software Foundation, either version 3 of the License, or |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
10 (at your option) any later version. |
| 7019 | 11 |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
12 Octave is distributed in the hope that it will be useful, but |
|
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
13 WITHOUT ANY WARRANTY; without even the implied warranty of |
|
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
15 GNU General Public License for more details. |
| 7019 | 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 | |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23662
diff
changeset
|
19 <https://www.gnu.org/licenses/>. |
| 7019 | 20 |
| 21 */ | |
| 22 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
23 /* |
| 5742 | 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: | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
38 |
| 5742 | 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 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
46 3. The names of its contributors may not be used to endorse or promote |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
47 products derived from this software without specific prior written |
| 5742 | 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 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
53 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER |
| 5742 | 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 | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
86 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
87 * 2006-04-01 David Bateman |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
88 * * convert for use in octave, declaring static functions only used |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
89 * here and adding oct_ to functions visible externally |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
90 * * inverse sense of ALLBITS |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
91 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
92 * 2012-05-18 David Bateman |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
93 * * Remove randu64 and ALLBIT option |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
94 * * Add the single precision generators |
| 5742 | 95 */ |
| 96 | |
| 97 /* | |
| 98 === Build instructions === | |
| 99 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
100 Compile with -DHAVE_GETTIMEOFDAY if the gettimeofday function is |
| 5742 | 101 available. This is not necessary if your architecture has |
| 102 /dev/urandom defined. | |
| 103 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
104 Uses implicit -Di386 or explicit -DHAVE_X86_32 to determine if CPU=x86. |
|
21308
c53bfd6d8e08
maint: Use American spelling for "behavior".
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
105 You can force X86 behavior with -DUSE_X86_32=1, or suppress it with |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
106 -DUSE_X86_32=0. You should also consider -march=i686 or similar for |
| 5742 | 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 === | |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
120 void oct_init_by_int (uint32_t s) 32-bit initial state |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
121 void oct_init_by_array (uint32_t k[],int m) m*32-bit initial state |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
122 void oct_init_by_entropy (void) random initial state |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
123 void oct_get_state (uint32_t save[MT_N+1]) saves state in array |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
124 void oct_set_state (uint32_t save[MT_N+1]) restores state from array |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
125 static uint32_t randmt (void) returns 32-bit unsigned int |
| 5742 | 126 |
| 127 === inline generators === | |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
128 static uint32_t randi32 (void) returns 32-bit unsigned int |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
129 static uint64_t randi53 (void) returns 53-bit unsigned int |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
130 static uint64_t randi54 (void) returns 54-bit unsigned int |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
131 static float randu32 (void) returns 32-bit uniform in (0,1) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
132 static double randu53 (void) returns 53-bit uniform in (0,1) |
| 5742 | 133 |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
134 double oct_randu (void) returns M-bit uniform in (0,1) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
135 double oct_randn (void) returns M-bit standard normal |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
136 double oct_rande (void) returns N-bit standard exponential |
| 5742 | 137 |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
138 float oct_float_randu (void) returns M-bit uniform in (0,1) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
139 float oct_float_randn (void) returns M-bit standard normal |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
140 float oct_float_rande (void) returns N-bit standard exponential |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
141 |
| 5742 | 142 === Array generators === |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
143 void oct_fill_randi32 (octave_idx_type, uint32_t []) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
144 void oct_fill_randi64 (octave_idx_type, uint64_t []) |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
145 |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
146 void oct_fill_randu (octave_idx_type, double []) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
147 void oct_fill_randn (octave_idx_type, double []) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
148 void oct_fill_rande (octave_idx_type, double []) |
| 5742 | 149 |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
150 void oct_fill_float_randu (octave_idx_type, float []) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
151 void oct_fill_float_randn (octave_idx_type, float []) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
152 void oct_fill_float_rande (octave_idx_type, float []) |
| 5742 | 153 */ |
| 154 | |
| 155 #if defined (HAVE_CONFIG_H) | |
|
21301
40de9f8f23a6
Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents:
21213
diff
changeset
|
156 # include "config.h" |
| 5742 | 157 #endif |
| 158 | |
|
23662
bd77ab816e43
eliminate obsolete file lo-math.h
John W. Eaton <jwe@octave.org>
parents:
23649
diff
changeset
|
159 #include <cmath> |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
160 #include <cstdio> |
| 5742 | 161 |
|
24856
8bb0251fcfde
Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents:
24534
diff
changeset
|
162 #include <algorithm> |
|
8bb0251fcfde
Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents:
24534
diff
changeset
|
163 |
|
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
164 #include "oct-time.h" |
| 5742 | 165 #include "randmtzig.h" |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
166 |
|
21066
258c787cd9ce
maint: Use "FIXME:" consistently in code base.
Rik <rik@octave.org>
parents:
20955
diff
changeset
|
167 /* FIXME: may want to suppress X86 if sizeof(long) > 4 */ |
|
20791
f7084eae3318
maint: Use Octave coding conventions for #if statements.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
168 #if ! defined (USE_X86_32) |
|
21202
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
169 # if defined (i386) || defined (HAVE_X86_32) |
|
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
170 # define USE_X86_32 1 |
|
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
171 # else |
|
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
172 # define USE_X86_32 0 |
|
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
173 # endif |
| 5742 | 174 #endif |
| 175 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
176 /* ===== Mersenne Twister 32-bit generator ===== */ |
| 5742 | 177 |
| 178 #define MT_M 397 | |
| 179 #define MATRIX_A 0x9908b0dfUL /* constant vector a */ | |
| 180 #define UMASK 0x80000000UL /* most significant w-r bits */ | |
| 181 #define LMASK 0x7fffffffUL /* least significant r bits */ | |
| 182 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) | |
| 183 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) | |
| 184 | |
| 185 static uint32_t *next; | |
| 186 static uint32_t state[MT_N]; /* the array for the state vector */ | |
| 187 static int left = 1; | |
| 188 static int initf = 0; | |
| 189 static int initt = 1; | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
190 static int inittf = 1; |
| 5742 | 191 |
| 192 /* initializes state[MT_N] with a seed */ | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
193 void |
|
24856
8bb0251fcfde
Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents:
24534
diff
changeset
|
194 oct_init_by_int (const uint32_t s) |
| 5742 | 195 { |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
196 int j; |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
197 state[0] = s & 0xffffffffUL; |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
198 for (j = 1; j < MT_N; j++) |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
199 { |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
200 state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j); |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
201 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
202 /* In the previous versions, MSBs of the seed affect */ |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
203 /* only MSBs of the array state[]. */ |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
204 /* 2002/01/09 modified by Makoto Matsumoto */ |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
205 state[j] &= 0xffffffffUL; /* for >32 bit machines */ |
| 5742 | 206 } |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
207 left = 1; |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
208 initf = 1; |
| 5742 | 209 } |
| 210 | |
| 211 /* initialize by an array with array-length */ | |
| 212 /* init_key is the array for initializing keys */ | |
| 213 /* key_length is its length */ | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
214 void |
|
24856
8bb0251fcfde
Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents:
24534
diff
changeset
|
215 oct_init_by_array (const uint32_t *init_key, const int key_length) |
| 5742 | 216 { |
| 217 int i, j, k; | |
| 218 oct_init_by_int (19650218UL); | |
| 219 i = 1; | |
| 220 j = 0; | |
| 221 k = (MT_N > key_length ? MT_N : key_length); | |
| 222 for (; k; k--) | |
| 223 { | |
| 224 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL)) | |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
225 + init_key[j] + j; /* non linear */ |
| 5742 | 226 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ |
| 227 i++; | |
| 228 j++; | |
| 229 if (i >= MT_N) | |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
230 { |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
231 state[0] = state[MT_N-1]; |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
232 i = 1; |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
233 } |
| 5742 | 234 if (j >= key_length) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
235 j = 0; |
| 5742 | 236 } |
| 237 for (k = MT_N - 1; k; k--) | |
| 238 { | |
| 239 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL)) | |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
240 - i; /* non linear */ |
| 5742 | 241 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ |
| 242 i++; | |
| 243 if (i >= MT_N) | |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
244 { |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
245 state[0] = state[MT_N-1]; |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
246 i = 1; |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
247 } |
| 5742 | 248 } |
| 249 | |
|
18812
9ac2357f19bc
doc: Replace "non-zero" with "nonzero" to match existing usage.
Rik <rik@octave.org>
parents:
17769
diff
changeset
|
250 state[0] = 0x80000000UL; /* MSB is 1; assuring nonzero initial array */ |
| 5742 | 251 left = 1; |
| 252 initf = 1; | |
| 253 } | |
| 254 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
255 void |
| 5742 | 256 oct_init_by_entropy (void) |
| 257 { | |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
258 uint32_t entropy[MT_N]; |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
259 int n = 0; |
| 5742 | 260 |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
261 /* Look for entropy in /dev/urandom */ |
|
23449
c763214a8260
maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
262 FILE *urandom = std::fopen ("/dev/urandom", "rb"); |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
263 if (urandom) |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
264 { |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
265 while (n < MT_N) |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
266 { |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
267 unsigned char word[4]; |
|
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
268 if (std::fread (word, 4, 1, urandom) != 1) |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
269 break; |
|
22402
4caa7b28d183
maint: Style check C++ code in liboctave/
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
270 entropy[n++] = word[0] + (word[1]<<8) + (word[2]<<16) |
|
4caa7b28d183
maint: Style check C++ code in liboctave/
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
271 + (static_cast<uint32_t> (word[3])<<24); |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
272 } |
|
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
273 std::fclose (urandom); |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
274 } |
| 5742 | 275 |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
276 /* If there isn't enough entropy, gather some from various sources */ |
|
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
277 |
|
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
278 octave::sys::time now; |
|
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
279 |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
280 if (n < MT_N) |
|
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
281 entropy[n++] = now.unix_time (); /* Current time in seconds */ |
|
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
282 |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
283 if (n < MT_N) |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
284 entropy[n++] = clock (); /* CPU time used (usec) */ |
|
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
285 |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
286 if (n < MT_N) |
|
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
287 entropy[n++] = now.usec (); /* Fractional part of current time */ |
|
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
288 |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
289 /* Send all the entropy into the initial state vector */ |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
290 oct_init_by_array (entropy,n); |
| 5742 | 291 } |
| 292 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
293 void |
|
24856
8bb0251fcfde
Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents:
24534
diff
changeset
|
294 oct_set_state (const uint32_t *save) |
| 5742 | 295 { |
|
24856
8bb0251fcfde
Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents:
24534
diff
changeset
|
296 std::copy_n (save, MT_N, state); |
| 5742 | 297 left = save[MT_N]; |
| 298 next = state + (MT_N - left + 1); | |
| 299 } | |
| 300 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
301 void |
|
7533
ff52243af934
save state separately for each MT random number generator
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
302 oct_get_state (uint32_t *save) |
| 5742 | 303 { |
|
24856
8bb0251fcfde
Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents:
24534
diff
changeset
|
304 std::copy_n (state, MT_N, save); |
| 5742 | 305 save[MT_N] = left; |
| 306 } | |
| 307 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
308 static void |
| 5742 | 309 next_state (void) |
| 310 { | |
| 311 uint32_t *p = state; | |
| 312 int j; | |
| 313 | |
| 314 /* if init_by_int() has not been called, */ | |
| 315 /* a default initial seed is used */ | |
| 316 /* if (initf==0) init_by_int(5489UL); */ | |
| 317 /* Or better yet, a random seed! */ | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
318 if (initf == 0) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
319 oct_init_by_entropy (); |
| 5742 | 320 |
| 321 left = MT_N; | |
| 322 next = state; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
323 |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
324 for (j = MT_N - MT_M + 1; --j; p++) |
| 5742 | 325 *p = p[MT_M] ^ TWIST(p[0], p[1]); |
| 326 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
327 for (j = MT_M; --j; p++) |
| 5742 | 328 *p = p[MT_M-MT_N] ^ TWIST(p[0], p[1]); |
| 329 | |
| 330 *p = p[MT_M-MT_N] ^ TWIST(p[0], state[0]); | |
| 331 } | |
| 332 | |
| 333 /* generates a random number on [0,0xffffffff]-interval */ | |
| 5766 | 334 static uint32_t |
| 5742 | 335 randmt (void) |
| 336 { | |
|
21332
db574711df78
eliminate deprecated register storage class specifier
John W. Eaton <jwe@octave.org>
parents:
21308
diff
changeset
|
337 uint32_t y; |
| 5742 | 338 |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
339 if (--left == 0) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
340 next_state (); |
| 5742 | 341 y = *next++; |
| 342 | |
| 343 /* Tempering */ | |
| 344 y ^= (y >> 11); | |
| 345 y ^= (y << 7) & 0x9d2c5680UL; | |
| 346 y ^= (y << 15) & 0xefc60000UL; | |
| 347 return (y ^ (y >> 18)); | |
| 348 } | |
| 349 | |
| 350 /* ===== Uniform generators ===== */ | |
| 351 | |
| 352 /* Select which 32 bit generator to use */ | |
| 353 #define randi32 randmt | |
| 354 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
355 static uint64_t |
| 5742 | 356 randi53 (void) |
| 357 { | |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
358 const uint32_t lo = randi32 (); |
|
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
359 const uint32_t hi = randi32 () & 0x1FFFFF; |
|
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21720
diff
changeset
|
360 #if defined (HAVE_X86_32) |
| 5742 | 361 uint64_t u; |
| 362 uint32_t *p = (uint32_t *)&u; | |
| 363 p[0] = lo; | |
| 364 p[1] = hi; | |
| 365 return u; | |
| 366 #else | |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
367 return ((static_cast<uint64_t> (hi) << 32) | lo); |
| 5742 | 368 #endif |
| 369 } | |
| 370 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
371 static uint64_t |
| 5742 | 372 randi54 (void) |
| 373 { | |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
374 const uint32_t lo = randi32 (); |
|
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
375 const uint32_t hi = randi32 () & 0x3FFFFF; |
|
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21720
diff
changeset
|
376 #if defined (HAVE_X86_32) |
| 5742 | 377 uint64_t u; |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
378 uint32_t *p = static_cast<uint32_t *> (&u); |
| 5742 | 379 p[0] = lo; |
| 380 p[1] = hi; | |
| 381 return u; | |
| 382 #else | |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
383 return ((static_cast<uint64_t> (hi) << 32) | lo); |
| 5742 | 384 #endif |
| 385 } | |
| 386 | |
| 387 /* generates a random number on (0,1)-real-interval */ | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
388 static float |
| 5742 | 389 randu32 (void) |
| 390 { | |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
391 return (static_cast<float> (randi32 ()) + 0.5) * (1.0/4294967296.0); |
| 5742 | 392 /* divided by 2^32 */ |
| 393 } | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
394 |
| 5742 | 395 /* generates a random number on (0,1) with 53-bit resolution */ |
| 5766 | 396 static double |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
397 randu53 (void) |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
398 { |
|
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
399 const uint32_t a = randi32 () >> 5; |
|
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
400 const uint32_t b = randi32 () >> 6; |
| 6959 | 401 return (a*67108864.0+b+0.4) * (1.0/9007199254740992.0); |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
402 } |
| 5742 | 403 |
| 404 /* Determine mantissa for uniform doubles */ | |
| 405 double | |
| 406 oct_randu (void) | |
| 407 { | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
408 return randu53 (); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
409 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
410 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
411 /* Determine mantissa for uniform floats */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
412 float |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
413 oct_float_randu (void) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
414 { |
| 6959 | 415 return randu32 (); |
| 5742 | 416 } |
| 417 | |
| 418 /* ===== Ziggurat normal and exponential generators ===== */ | |
| 419 | |
| 420 #define ZIGGURAT_TABLE_SIZE 256 | |
| 421 | |
| 422 #define ZIGGURAT_NOR_R 3.6541528853610088 | |
| 423 #define ZIGGURAT_NOR_INV_R 0.27366123732975828 | |
| 424 #define NOR_SECTION_AREA 0.00492867323399 | |
| 425 | |
| 426 #define ZIGGURAT_EXP_R 7.69711747013104972 | |
| 427 #define ZIGGURAT_EXP_INV_R 0.129918765548341586 | |
| 428 #define EXP_SECTION_AREA 0.0039496598225815571993 | |
| 429 | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
430 #define ZIGINT uint64_t |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
431 #define EMANTISSA 9007199254740992.0 /* 53 bit mantissa */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
432 #define ERANDI randi53() /* 53 bits for mantissa */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
433 #define NMANTISSA EMANTISSA |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
434 #define NRANDI randi54() /* 53 bits for mantissa + 1 bit sign */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
435 #define RANDU randu53() |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
436 |
| 5742 | 437 static ZIGINT ki[ZIGGURAT_TABLE_SIZE]; |
| 438 static double wi[ZIGGURAT_TABLE_SIZE], fi[ZIGGURAT_TABLE_SIZE]; | |
| 439 static ZIGINT ke[ZIGGURAT_TABLE_SIZE]; | |
| 440 static double we[ZIGGURAT_TABLE_SIZE], fe[ZIGGURAT_TABLE_SIZE]; | |
| 441 | |
| 442 /* | |
| 443 This code is based on the paper Marsaglia and Tsang, "The ziggurat method | |
| 444 for generating random variables", Journ. Statistical Software. Code was | |
| 445 presented in this paper for a Ziggurat of 127 levels and using a 32 bit | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
446 integer random number generator. This version of the code, uses the |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
447 Mersenne Twister as the integer generator and uses 256 levels in the |
| 5742 | 448 Ziggurat. This has several advantages. |
| 449 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
450 1) As Marsaglia and Tsang themselves states, the more levels the few |
| 5742 | 451 times the expensive tail algorithm must be called |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
452 2) The cycle time of the generator is determined by the integer |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
453 generator, thus the use of a Mersenne Twister for the core random |
| 5742 | 454 generator makes this cycle extremely long. |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
455 3) The license on the original code was unclear, thus rewriting the code |
| 5742 | 456 from the article means we are free of copyright issues. |
| 457 4) Compile flag for full 53-bit random mantissa. | |
| 458 | |
| 459 It should be stated that the authors made my life easier, by the fact that | |
| 460 the algorithm developed in the text of the article is for a 256 level | |
| 461 ziggurat, even if the code itself isn't... | |
| 462 | |
| 463 One modification to the algorithm developed in the article, is that it is | |
| 464 assumed that 0 <= x < Inf, and "unsigned long"s are used, thus resulting in | |
| 465 terms like 2^32 in the code. As the normal distribution is defined between | |
| 466 -Inf < x < Inf, we effectively only have 31 bit integers plus a sign. Thus | |
| 467 in Marsaglia and Tsang, terms like 2^32 become 2^31. We use NMANTISSA for | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
468 this term. The exponential distribution is one sided so we use the |
| 5742 | 469 full 32 bits. We use EMANTISSA for this term. |
| 470 | |
| 471 It appears that I'm slightly slower than the code in the article, this | |
| 472 is partially due to a better generator of random integers than they | |
| 473 use. But might also be that the case of rapid return was optimized by | |
| 474 inlining the relevant code with a #define. As the basic Mersenne | |
| 475 Twister is only 25% faster than this code I suspect that the main | |
| 476 reason is just the use of the Mersenne Twister and not the inlining, | |
| 477 so I'm not going to try and optimize further. | |
| 478 */ | |
| 479 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
480 static void |
| 5742 | 481 create_ziggurat_tables (void) |
| 482 { | |
| 483 int i; | |
| 484 double x, x1; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
485 |
| 5742 | 486 /* Ziggurat tables for the normal distribution */ |
| 487 x1 = ZIGGURAT_NOR_R; | |
| 488 wi[255] = x1 / NMANTISSA; | |
| 489 fi[255] = exp (-0.5 * x1 * x1); | |
| 490 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
491 /* Index zero is special for tail strip, where Marsaglia and Tsang |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
492 * defines this as |
| 5742 | 493 * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1, |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
494 * where v is the area of each strip of the ziggurat. |
| 5742 | 495 */ |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
496 ki[0] = static_cast<ZIGINT> (x1 * fi[255] / NOR_SECTION_AREA * NMANTISSA); |
| 5742 | 497 wi[0] = NOR_SECTION_AREA / fi[255] / NMANTISSA; |
| 498 fi[0] = 1.; | |
| 499 | |
| 500 for (i = 254; i > 0; i--) | |
| 501 { | |
| 502 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus | |
| 503 * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y)) | |
| 504 */ | |
|
23649
aabf6cd222ac
Use sqrt from C++ std library in liboctave.
Rik <rik@octave.org>
parents:
23449
diff
changeset
|
505 x = std::sqrt (-2. * std::log (NOR_SECTION_AREA / x1 + fi[i+1])); |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
506 ki[i+1] = static_cast<ZIGINT> (x / x1 * NMANTISSA); |
| 5742 | 507 wi[i] = x / NMANTISSA; |
| 508 fi[i] = exp (-0.5 * x * x); | |
| 509 x1 = x; | |
| 510 } | |
| 511 | |
| 512 ki[1] = 0; | |
| 513 | |
| 514 /* Zigurrat tables for the exponential distribution */ | |
| 515 x1 = ZIGGURAT_EXP_R; | |
| 516 we[255] = x1 / EMANTISSA; | |
| 517 fe[255] = exp (-x1); | |
| 518 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
519 /* Index zero is special for tail strip, where Marsaglia and Tsang |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
520 * defines this as |
| 5742 | 521 * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1, |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
522 * where v is the area of each strip of the ziggurat. |
| 5742 | 523 */ |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
524 ke[0] = static_cast<ZIGINT> (x1 * fe[255] / EXP_SECTION_AREA * EMANTISSA); |
| 5742 | 525 we[0] = EXP_SECTION_AREA / fe[255] / EMANTISSA; |
| 526 fe[0] = 1.; | |
| 527 | |
| 528 for (i = 254; i > 0; i--) | |
| 529 { | |
| 530 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus | |
| 531 * need inverse operator of y = exp(-x) -> x = -ln(y) | |
| 532 */ | |
|
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
533 x = - std::log (EXP_SECTION_AREA / x1 + fe[i+1]); |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
534 ke[i+1] = static_cast<ZIGINT> (x / x1 * EMANTISSA); |
| 5742 | 535 we[i] = x / EMANTISSA; |
| 536 fe[i] = exp (-x); | |
| 537 x1 = x; | |
| 538 } | |
| 539 ke[1] = 0; | |
| 540 | |
| 541 initt = 0; | |
| 542 } | |
| 543 | |
| 544 /* | |
| 545 * Here is the guts of the algorithm. As Marsaglia and Tsang state the | |
| 546 * algorithm in their paper | |
| 547 * | |
| 548 * 1) Calculate a random signed integer j and let i be the index | |
| 549 * provided by the rightmost 8-bits of j | |
| 550 * 2) Set x = j * w_i. If j < k_i return x | |
| 551 * 3) If i = 0, then return x from the tail | |
| 552 * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x | |
| 553 * 5) goto step 1 | |
| 554 * | |
| 555 * Where f is the functional form of the distribution, which for a normal | |
| 556 * distribution is exp(-0.5*x*x) | |
| 557 */ | |
| 558 | |
| 559 double | |
| 560 oct_randn (void) | |
| 561 { | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
562 if (initt) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
563 create_ziggurat_tables (); |
| 5742 | 564 |
| 565 while (1) | |
| 566 { | |
| 567 /* The following code is specialized for 32-bit mantissa. | |
| 568 * Compared to the arbitrary mantissa code, there is a performance | |
| 569 * gain for 32-bits: PPC: 2%, MIPS: 8%, x86: 40% | |
| 570 * There is a bigger performance gain compared to using a full | |
| 571 * 53-bit mantissa: PPC: 60%, MIPS: 65%, x86: 240% | |
| 572 * Of course, different compilers and operating systems may | |
| 573 * have something to do with this. | |
| 574 */ | |
|
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21720
diff
changeset
|
575 # if defined (HAVE_X86_32) |
| 5742 | 576 /* 53-bit mantissa, 1-bit sign, x86 32-bit architecture */ |
| 577 double x; | |
| 578 int si,idx; | |
|
21332
db574711df78
eliminate deprecated register storage class specifier
John W. Eaton <jwe@octave.org>
parents:
21308
diff
changeset
|
579 uint32_t lo, hi; |
| 5742 | 580 int64_t rabs; |
| 581 uint32_t *p = (uint32_t *)&rabs; | |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
582 lo = randi32 (); |
|
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
583 idx = lo & 0xFF; |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
584 hi = randi32 (); |
|
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
585 si = hi & UMASK; |
| 5742 | 586 p[0] = lo; |
|
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
587 p[1] = hi & 0x1FFFFF; |
| 5742 | 588 x = ( si ? -rabs : rabs ) * wi[idx]; |
| 22188 | 589 # else |
| 5742 | 590 /* arbitrary mantissa (selected by NRANDI, with 1 bit for sign) */ |
| 591 const uint64_t r = NRANDI; | |
|
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
592 const int64_t rabs = r >> 1; |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
593 const int idx = static_cast<int> (rabs & 0xFF); |
|
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
594 const double x = ( (r & 1) ? -rabs : rabs) * wi[idx]; |
| 22188 | 595 # endif |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
596 if (rabs < static_cast<int64_t> (ki[idx])) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
597 return x; /* 99.3% of the time we return here 1st try */ |
| 5742 | 598 else if (idx == 0) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
599 { |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
600 /* As stated in Marsaglia and Tsang |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
601 * |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
602 * For the normal tail, the method of Marsaglia[5] provides: |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
603 * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x, |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
604 * then return r+x. Except that r+x is always in the positive |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
605 * tail!!!! Any thing random might be used to determine the |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
606 * sign, but as we already have r we might as well use it |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
607 * |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
608 * [PAK] but not the bottom 8 bits, since they are all 0 here! |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
609 */ |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
610 double xx, yy; |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
611 do |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
612 { |
|
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
613 xx = - ZIGGURAT_NOR_INV_R * std::log (RANDU); |
|
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
614 yy = - std::log (RANDU); |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
615 } |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
616 while ( yy+yy <= xx*xx); |
|
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
617 return ((rabs & 0x100) ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx); |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
618 } |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
619 else if ((fi[idx-1] - fi[idx]) * RANDU + fi[idx] < exp (-0.5*x*x)) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
620 return x; |
| 5742 | 621 } |
| 622 } | |
| 623 | |
| 624 double | |
| 625 oct_rande (void) | |
| 626 { | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
627 if (initt) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
628 create_ziggurat_tables (); |
| 5742 | 629 |
| 630 while (1) | |
| 631 { | |
| 632 ZIGINT ri = ERANDI; | |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
633 const int idx = static_cast<int> (ri & 0xFF); |
| 5742 | 634 const double x = ri * we[idx]; |
| 635 if (ri < ke[idx]) | |
|
15970
ca6202597201
Use C-style comments in C code.
Mike Miller <mtmiller@ieee.org>
parents:
15408
diff
changeset
|
636 return x; /* 98.9% of the time we return here 1st try */ |
| 5742 | 637 else if (idx == 0) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
638 { |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
639 /* As stated in Marsaglia and Tsang |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
640 * |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
641 * For the exponential tail, the method of Marsaglia[5] provides: |
| 5742 | 642 * x = r - ln(U); |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
643 */ |
|
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
644 return ZIGGURAT_EXP_R - std::log (RANDU); |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
645 } |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
646 else if ((fe[idx-1] - fe[idx]) * RANDU + fe[idx] < exp (-x)) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
647 return x; |
| 5742 | 648 } |
| 649 } | |
| 650 | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
651 #undef ZIGINT |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
652 #undef EMANTISSA |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
653 #undef ERANDI |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
654 #undef NMANTISSA |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
655 #undef NRANDI |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
656 #undef RANDU |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
657 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
658 #define ZIGINT uint32_t |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
659 #define EMANTISSA 4294967296.0 /* 32 bit mantissa */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
660 #define ERANDI randi32() /* 32 bits for mantissa */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
661 #define NMANTISSA 2147483648.0 /* 31 bit mantissa */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
662 #define NRANDI randi32() /* 31 bits for mantissa + 1 bit sign */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
663 #define RANDU randu32() |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
664 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
665 static ZIGINT fki[ZIGGURAT_TABLE_SIZE]; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
666 static float fwi[ZIGGURAT_TABLE_SIZE], ffi[ZIGGURAT_TABLE_SIZE]; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
667 static ZIGINT fke[ZIGGURAT_TABLE_SIZE]; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
668 static float fwe[ZIGGURAT_TABLE_SIZE], ffe[ZIGGURAT_TABLE_SIZE]; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
669 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
670 static void |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
671 create_ziggurat_float_tables (void) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
672 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
673 int i; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
674 float x, x1; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
675 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
676 /* Ziggurat tables for the normal distribution */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
677 x1 = ZIGGURAT_NOR_R; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
678 fwi[255] = x1 / NMANTISSA; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
679 ffi[255] = exp (-0.5 * x1 * x1); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
680 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
681 /* Index zero is special for tail strip, where Marsaglia and Tsang |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
682 * defines this as |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
683 * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1, |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
684 * where v is the area of each strip of the ziggurat. |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
685 */ |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
686 fki[0] = static_cast<ZIGINT> (x1 * ffi[255] / NOR_SECTION_AREA * NMANTISSA); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
687 fwi[0] = NOR_SECTION_AREA / ffi[255] / NMANTISSA; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
688 ffi[0] = 1.; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
689 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
690 for (i = 254; i > 0; i--) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
691 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
692 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
693 * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y)) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
694 */ |
|
23649
aabf6cd222ac
Use sqrt from C++ std library in liboctave.
Rik <rik@octave.org>
parents:
23449
diff
changeset
|
695 x = std::sqrt (-2. * std::log (NOR_SECTION_AREA / x1 + ffi[i+1])); |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
696 fki[i+1] = static_cast<ZIGINT> (x / x1 * NMANTISSA); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
697 fwi[i] = x / NMANTISSA; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
698 ffi[i] = exp (-0.5 * x * x); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
699 x1 = x; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
700 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
701 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
702 fki[1] = 0; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
703 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
704 /* Zigurrat tables for the exponential distribution */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
705 x1 = ZIGGURAT_EXP_R; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
706 fwe[255] = x1 / EMANTISSA; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
707 ffe[255] = exp (-x1); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
708 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
709 /* Index zero is special for tail strip, where Marsaglia and Tsang |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
710 * defines this as |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
711 * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1, |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
712 * where v is the area of each strip of the ziggurat. |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
713 */ |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
714 fke[0] = static_cast<ZIGINT> (x1 * ffe[255] / EXP_SECTION_AREA * EMANTISSA); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
715 fwe[0] = EXP_SECTION_AREA / ffe[255] / EMANTISSA; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
716 ffe[0] = 1.; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
717 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
718 for (i = 254; i > 0; i--) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
719 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
720 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
721 * need inverse operator of y = exp(-x) -> x = -ln(y) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
722 */ |
|
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
723 x = - std::log (EXP_SECTION_AREA / x1 + ffe[i+1]); |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
724 fke[i+1] = static_cast<ZIGINT> (x / x1 * EMANTISSA); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
725 fwe[i] = x / EMANTISSA; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
726 ffe[i] = exp (-x); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
727 x1 = x; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
728 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
729 fke[1] = 0; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
730 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
731 inittf = 0; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
732 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
733 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
734 /* |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
735 * Here is the guts of the algorithm. As Marsaglia and Tsang state the |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
736 * algorithm in their paper |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
737 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
738 * 1) Calculate a random signed integer j and let i be the index |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
739 * provided by the rightmost 8-bits of j |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
740 * 2) Set x = j * w_i. If j < k_i return x |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
741 * 3) If i = 0, then return x from the tail |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
742 * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
743 * 5) goto step 1 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
744 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
745 * Where f is the functional form of the distribution, which for a normal |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
746 * distribution is exp(-0.5*x*x) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
747 */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
748 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
749 float |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
750 oct_float_randn (void) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
751 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
752 if (inittf) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
753 create_ziggurat_float_tables (); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
754 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
755 while (1) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
756 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
757 /* 32-bit mantissa */ |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
758 const uint32_t r = randi32 (); |
|
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
759 const uint32_t rabs = r & LMASK; |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
760 const int idx = static_cast<int> (r & 0xFF); |
|
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
761 const float x = static_cast<int32_t> (r) * fwi[idx]; |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
762 if (rabs < fki[idx]) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
763 return x; /* 99.3% of the time we return here 1st try */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
764 else if (idx == 0) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
765 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
766 /* As stated in Marsaglia and Tsang |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
767 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
768 * For the normal tail, the method of Marsaglia[5] provides: |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
769 * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x, |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
770 * then return r+x. Except that r+x is always in the positive |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
771 * tail!!!! Any thing random might be used to determine the |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
772 * sign, but as we already have r we might as well use it |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
773 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
774 * [PAK] but not the bottom 8 bits, since they are all 0 here! |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
775 */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
776 float xx, yy; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
777 do |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
778 { |
|
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
779 xx = - ZIGGURAT_NOR_INV_R * std::log (RANDU); |
|
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
780 yy = - std::log (RANDU); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
781 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
782 while ( yy+yy <= xx*xx); |
|
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
783 return ((rabs & 0x100) ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
784 } |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
785 else if ((ffi[idx-1] - ffi[idx]) * RANDU + ffi[idx] < exp (-0.5*x*x)) |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
786 return x; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
787 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
788 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
789 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
790 float |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
791 oct_float_rande (void) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
792 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
793 if (inittf) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
794 create_ziggurat_float_tables (); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
795 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
796 while (1) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
797 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
798 ZIGINT ri = ERANDI; |
|
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
799 const int idx = static_cast<int> (ri & 0xFF); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
800 const float x = ri * fwe[idx]; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
801 if (ri < fke[idx]) |
|
15970
ca6202597201
Use C-style comments in C code.
Mike Miller <mtmiller@ieee.org>
parents:
15408
diff
changeset
|
802 return x; /* 98.9% of the time we return here 1st try */ |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
803 else if (idx == 0) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
804 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
805 /* As stated in Marsaglia and Tsang |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
806 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
807 * For the exponential tail, the method of Marsaglia[5] provides: |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
808 * x = r - ln(U); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
809 */ |
|
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
810 return ZIGGURAT_EXP_R - std::log (RANDU); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
811 } |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
812 else if ((ffe[idx-1] - ffe[idx]) * RANDU + ffe[idx] < exp (-x)) |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
813 return x; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
814 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
815 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
816 |
| 5742 | 817 /* Array generators */ |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
818 void |
| 5742 | 819 oct_fill_randu (octave_idx_type n, double *p) |
| 820 { | |
| 821 octave_idx_type i; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
822 for (i = 0; i < n; i++) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
823 p[i] = oct_randu (); |
| 5742 | 824 } |
| 825 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
826 void |
| 5742 | 827 oct_fill_randn (octave_idx_type n, double *p) |
| 828 { | |
| 829 octave_idx_type i; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
830 for (i = 0; i < n; i++) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
831 p[i] = oct_randn (); |
| 5742 | 832 } |
| 833 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
834 void |
| 5742 | 835 oct_fill_rande (octave_idx_type n, double *p) |
| 836 { | |
| 837 octave_idx_type i; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
838 for (i = 0; i < n; i++) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
839 p[i] = oct_rande (); |
| 5742 | 840 } |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
841 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
842 void |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
843 oct_fill_float_randu (octave_idx_type n, float *p) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
844 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
845 octave_idx_type i; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
846 for (i = 0; i < n; i++) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
847 p[i] = oct_float_randu (); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
848 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
849 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
850 void |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
851 oct_fill_float_randn (octave_idx_type n, float *p) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
852 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
853 octave_idx_type i; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
854 for (i = 0; i < n; i++) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
855 p[i] = oct_float_randn (); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
856 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
857 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
858 void |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
859 oct_fill_float_rande (octave_idx_type n, float *p) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
860 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
861 octave_idx_type i; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
862 for (i = 0; i < n; i++) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
863 p[i] = oct_float_rande (); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
864 } |