Mercurial > octave
annotate liboctave/numeric/randmtzig.cc @ 22022:278fc29b69ca
maint: Eliminate unnecessary double newline blocks (only whitespace changes).
* Figure.cc, annotation-dialog.cc, dialog.cc, dialog.h, files-dock-widget.cc,
find-files-dialog.cc, find-files-model.cc, file-editor-tab.cc, file-editor.cc,
find-dialog.cc, octave-qscintilla.cc, octave-txt-lexer.cc, main-window.cc,
main-window.h, octave-cmd.cc, octave-cmd.h, octave-dock-widget.cc,
octave-qt-link.cc, parser.h, webinfo.cc, settings-dialog.cc,
shortcut-manager.cc, shortcut-manager.h, welcome-wizard.cc, workspace-model.cc,
__ilu__.cc, __lin_interpn__.cc, besselj.cc, bsxfun.cc, cellfun.cc, conv2.cc,
data.cc, debug.cc, debug.h, dynamic-ld.cc, fft2.cc, file-io.cc,
ft-text-renderer.cc, gl-render.cc, gl2ps-print.cc, graphics.cc, graphics.in.h,
help.cc, jit-typeinfo.cc, jit-typeinfo.h, kron.cc, ls-ascii-helper.cc,
ls-hdf5.cc, ls-hdf5.h, ls-mat5.cc, max.cc, oct-map.cc, oct-map.h,
oct-stream.cc, octave-link.h, pr-output.cc, pt-jit.cc, qz.cc, sparse-xdiv.cc,
sparse-xpow.cc, strfind.cc, strfns.cc, symtab.h, toplev.cc, toplev.h,
typecast.cc, utils.cc, variables.cc, xdiv.h, xpow.cc, xpow.h, zfstream.cc,
__eigs__.cc, __init_fltk__.cc, audioread.cc, ccolamd.cc, chol.cc, qr.cc,
ov-base-diag.cc, ov-base-int.cc, ov-base-sparse.cc, ov-base.cc, ov-base.h,
ov-bool-mat.cc, ov-bool-mat.h, ov-bool-sparse.cc, ov-bool-sparse.h, ov-bool.cc,
ov-bool.h, ov-builtin.h, ov-cell.cc, ov-cell.h, ov-class.cc, ov-class.h,
ov-classdef.cc, ov-complex.cc, ov-complex.h, ov-cs-list.cc, ov-cs-list.h,
ov-cx-diag.cc, ov-cx-diag.h, ov-cx-mat.cc, ov-cx-mat.h, ov-cx-sparse.cc,
ov-cx-sparse.h, ov-dld-fcn.cc, ov-dld-fcn.h, ov-fcn-handle.h, ov-fcn-inline.h,
ov-float.cc, ov-float.h, ov-flt-complex.cc, ov-flt-complex.h,
ov-flt-cx-diag.cc, ov-flt-cx-diag.h, ov-flt-cx-mat.cc, ov-flt-cx-mat.h,
ov-flt-re-diag.cc, ov-flt-re-diag.h, ov-flt-re-mat.cc, ov-int16.cc,
ov-int32.cc, ov-int64.cc, ov-int8.cc, ov-intx.h, ov-lazy-idx.cc, ov-lazy-idx.h,
ov-mex-fcn.h, ov-null-mat.h, ov-perm.cc, ov-perm.h, ov-range.cc, ov-range.h,
ov-re-diag.cc, ov-re-diag.h, ov-re-mat.cc, ov-re-sparse.cc, ov-re-sparse.h,
ov-scalar.cc, ov-scalar.h, ov-str-mat.h, ov-struct.cc, ov-struct.h,
ov-type-conv.h, ov-uint16.cc, ov-uint32.cc, ov-uint64.cc, ov-uint8.cc,
ov-usr-fcn.cc, ov-usr-fcn.h, ov.h, op-cm-m.cc, op-fcm-fm.cc, op-fm-fm.cc,
op-int.h, op-m-m.cc, op-pm-template.cc, pt-binop.h, pt-select.cc, pt-stmt.cc,
pt-unop.h, Array-util.cc, Array.cc, Array.h, CColVector.cc, CMatrix.cc,
CMatrix.h, CSparse.cc, MArray.cc, MSparse.cc, MatrixType.cc, PermMatrix.cc,
Range.cc, Sparse-C.cc, Sparse.cc, dMatrix.cc, dim-vector.h, fCColVector.cc,
fCMatrix.cc, fMatrix.cc, fMatrix.h, idx-vector.cc, idx-vector.h, blaswrap.c,
f77-fcn.h, DASPK.cc, eigs-base.cc, hess.cc, lo-specfun.cc, oct-fftw.cc,
oct-norm.cc, oct-rand.cc, oct-rand.h, qr.cc, qrp.cc, randmtzig.cc,
randpoisson.cc, sparse-lu.cc, sparse-lu.h, sparse-qr.cc, Sparse-op-defs.h,
mx-inlines.cc, oct-env.h, cmd-edit.cc, f2c-main.c, lo-array-errwarn.cc,
lo-array-errwarn.h, lo-array-gripes.cc, lo-array-gripes.h, lo-ieee.h,
lo-regexp.cc, oct-binmap.h, oct-inttypes.cc, oct-inttypes.h, oct-sort.cc,
oct-sort.h, mkoctfile.in.cc:
maint: Eliminate unnecessary double newline blocks (only whitespace changes).
author | Rik <rik@octave.org> |
---|---|
date | Fri, 01 Jul 2016 19:48:35 -0700 |
parents | aab79a1885cc |
children | 1344509a480c |
rev | line source |
---|---|
7019 | 1 /* |
2 | |
19697
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
18812
diff
changeset
|
3 Copyright (C) 2006-2015 John W. Eaton |
7019 | 4 |
5 This file is part of Octave. | |
6 | |
7 Octave is free software; you can redistribute it and/or modify it | |
8 under the terms of the GNU General Public License as published by the | |
9 Free Software Foundation; either version 3 of the License, or (at your | |
10 option) any later version. | |
11 | |
12 Octave is distributed in the hope that it will be useful, but WITHOUT | |
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
15 for more details. | |
16 | |
17 You should have received a copy of the GNU General Public License | |
18 along with Octave; see the file COPYING. If not, see | |
19 <http://www.gnu.org/licenses/>. | |
20 | |
21 */ | |
22 | |
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 | |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
159 #include <cstdio> |
5742 | 160 |
7231 | 161 #include "lo-math.h" |
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
162 #include "oct-time.h" |
5742 | 163 #include "randmtzig.h" |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
164 |
21066
258c787cd9ce
maint: Use "FIXME:" consistently in code base.
Rik <rik@octave.org>
parents:
20955
diff
changeset
|
165 /* 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
|
166 #if ! defined (USE_X86_32) |
21202
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
167 # 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
|
168 # define USE_X86_32 1 |
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
169 # else |
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
170 # define USE_X86_32 0 |
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
171 # endif |
5742 | 172 #endif |
173 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
174 /* ===== Mersenne Twister 32-bit generator ===== */ |
5742 | 175 |
176 #define MT_M 397 | |
177 #define MATRIX_A 0x9908b0dfUL /* constant vector a */ | |
178 #define UMASK 0x80000000UL /* most significant w-r bits */ | |
179 #define LMASK 0x7fffffffUL /* least significant r bits */ | |
180 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) | |
181 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) | |
182 | |
183 static uint32_t *next; | |
184 static uint32_t state[MT_N]; /* the array for the state vector */ | |
185 static int left = 1; | |
186 static int initf = 0; | |
187 static int initt = 1; | |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
188 static int inittf = 1; |
5742 | 189 |
190 /* 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
|
191 void |
5742 | 192 oct_init_by_int (uint32_t s) |
193 { | |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
194 int j; |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
195 state[0] = s & 0xffffffffUL; |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
196 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
|
197 { |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
198 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
|
199 /* 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
|
200 /* 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
|
201 /* 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
|
202 /* 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
|
203 state[j] &= 0xffffffffUL; /* for >32 bit machines */ |
5742 | 204 } |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
205 left = 1; |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
206 initf = 1; |
5742 | 207 } |
208 | |
209 /* initialize by an array with array-length */ | |
210 /* init_key is the array for initializing keys */ | |
211 /* key_length is its length */ | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
212 void |
7533
ff52243af934
save state separately for each MT random number generator
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
213 oct_init_by_array (uint32_t *init_key, int key_length) |
5742 | 214 { |
215 int i, j, k; | |
216 oct_init_by_int (19650218UL); | |
217 i = 1; | |
218 j = 0; | |
219 k = (MT_N > key_length ? MT_N : key_length); | |
220 for (; k; k--) | |
221 { | |
222 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL)) | |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
223 + init_key[j] + j; /* non linear */ |
5742 | 224 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ |
225 i++; | |
226 j++; | |
227 if (i >= MT_N) | |
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
228 { |
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
229 state[0] = state[MT_N-1]; |
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
230 i = 1; |
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
231 } |
5742 | 232 if (j >= key_length) |
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
233 j = 0; |
5742 | 234 } |
235 for (k = MT_N - 1; k; k--) | |
236 { | |
237 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
|
238 - i; /* non linear */ |
5742 | 239 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ |
240 i++; | |
241 if (i >= MT_N) | |
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
242 { |
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
243 state[0] = state[MT_N-1]; |
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
244 i = 1; |
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
245 } |
5742 | 246 } |
247 | |
18812
9ac2357f19bc
doc: Replace "non-zero" with "nonzero" to match existing usage.
Rik <rik@octave.org>
parents:
17769
diff
changeset
|
248 state[0] = 0x80000000UL; /* MSB is 1; assuring nonzero initial array */ |
5742 | 249 left = 1; |
250 initf = 1; | |
251 } | |
252 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
253 void |
5742 | 254 oct_init_by_entropy (void) |
255 { | |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
256 uint32_t entropy[MT_N]; |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
257 int n = 0; |
5742 | 258 |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
259 /* Look for entropy in /dev/urandom */ |
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
260 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
|
261 if (urandom) |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
262 { |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
263 while (n < MT_N) |
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 unsigned char word[4]; |
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
266 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
|
267 break; |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
268 entropy[n++] = word[0]+(word[1]<<8)+(word[2]<<16)+(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
|
269 } |
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
270 std::fclose (urandom); |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
271 } |
5742 | 272 |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
273 /* 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
|
274 |
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
275 octave::sys::time now; |
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
276 |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
277 if (n < MT_N) |
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
278 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
|
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) |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
281 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
|
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) |
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
284 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
|
285 |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
286 /* 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
|
287 oct_init_by_array (entropy,n); |
5742 | 288 } |
289 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
290 void |
7533
ff52243af934
save state separately for each MT random number generator
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
291 oct_set_state (uint32_t *save) |
5742 | 292 { |
293 int i; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
294 for (i = 0; i < MT_N; i++) |
5742 | 295 state[i] = save[i]; |
296 left = save[MT_N]; | |
297 next = state + (MT_N - left + 1); | |
298 } | |
299 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
300 void |
7533
ff52243af934
save state separately for each MT random number generator
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
301 oct_get_state (uint32_t *save) |
5742 | 302 { |
303 int i; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
304 for (i = 0; i < MT_N; i++) |
5742 | 305 save[i] = state[i]; |
306 save[MT_N] = left; | |
307 } | |
308 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
309 static void |
5742 | 310 next_state (void) |
311 { | |
312 uint32_t *p = state; | |
313 int j; | |
314 | |
315 /* if init_by_int() has not been called, */ | |
316 /* a default initial seed is used */ | |
317 /* if (initf==0) init_by_int(5489UL); */ | |
318 /* Or better yet, a random seed! */ | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
319 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
|
320 oct_init_by_entropy (); |
5742 | 321 |
322 left = MT_N; | |
323 next = state; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
324 |
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
325 for (j = MT_N - MT_M + 1; --j; p++) |
5742 | 326 *p = p[MT_M] ^ TWIST(p[0], p[1]); |
327 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
328 for (j = MT_M; --j; p++) |
5742 | 329 *p = p[MT_M-MT_N] ^ TWIST(p[0], p[1]); |
330 | |
331 *p = p[MT_M-MT_N] ^ TWIST(p[0], state[0]); | |
332 } | |
333 | |
334 /* generates a random number on [0,0xffffffff]-interval */ | |
5766 | 335 static uint32_t |
5742 | 336 randmt (void) |
337 { | |
21332
db574711df78
eliminate deprecated register storage class specifier
John W. Eaton <jwe@octave.org>
parents:
21308
diff
changeset
|
338 uint32_t y; |
5742 | 339 |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
340 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
|
341 next_state (); |
5742 | 342 y = *next++; |
343 | |
344 /* Tempering */ | |
345 y ^= (y >> 11); | |
346 y ^= (y << 7) & 0x9d2c5680UL; | |
347 y ^= (y << 15) & 0xefc60000UL; | |
348 return (y ^ (y >> 18)); | |
349 } | |
350 | |
351 /* ===== Uniform generators ===== */ | |
352 | |
353 /* Select which 32 bit generator to use */ | |
354 #define randi32 randmt | |
355 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
356 static uint64_t |
5742 | 357 randi53 (void) |
358 { | |
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
|
359 const uint32_t lo = randi32 (); |
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
360 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
|
361 #if defined (HAVE_X86_32) |
5742 | 362 uint64_t u; |
363 uint32_t *p = (uint32_t *)&u; | |
364 p[0] = lo; | |
365 p[1] = hi; | |
366 return u; | |
367 #else | |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
368 return ((static_cast<uint64_t> (hi) << 32) | lo); |
5742 | 369 #endif |
370 } | |
371 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
372 static uint64_t |
5742 | 373 randi54 (void) |
374 { | |
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
|
375 const uint32_t lo = randi32 (); |
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
376 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
|
377 #if defined (HAVE_X86_32) |
5742 | 378 uint64_t u; |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
379 uint32_t *p = static_cast<uint32_t *> (&u); |
5742 | 380 p[0] = lo; |
381 p[1] = hi; | |
382 return u; | |
383 #else | |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
384 return ((static_cast<uint64_t> (hi) << 32) | lo); |
5742 | 385 #endif |
386 } | |
387 | |
388 /* 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
|
389 static float |
5742 | 390 randu32 (void) |
391 { | |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
392 return (static_cast<float> (randi32 ()) + 0.5) * (1.0/4294967296.0); |
5742 | 393 /* divided by 2^32 */ |
394 } | |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
395 |
5742 | 396 /* generates a random number on (0,1) with 53-bit resolution */ |
5766 | 397 static double |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
398 randu53 (void) |
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
399 { |
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
400 const uint32_t a = randi32 () >> 5; |
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
401 const uint32_t b = randi32 () >> 6; |
6959 | 402 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
|
403 } |
5742 | 404 |
405 /* Determine mantissa for uniform doubles */ | |
406 double | |
407 oct_randu (void) | |
408 { | |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
409 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
|
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 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
412 /* 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
|
413 float |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
414 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
|
415 { |
6959 | 416 return randu32 (); |
5742 | 417 } |
418 | |
419 /* ===== Ziggurat normal and exponential generators ===== */ | |
420 | |
421 #define ZIGGURAT_TABLE_SIZE 256 | |
422 | |
423 #define ZIGGURAT_NOR_R 3.6541528853610088 | |
424 #define ZIGGURAT_NOR_INV_R 0.27366123732975828 | |
425 #define NOR_SECTION_AREA 0.00492867323399 | |
426 | |
427 #define ZIGGURAT_EXP_R 7.69711747013104972 | |
428 #define ZIGGURAT_EXP_INV_R 0.129918765548341586 | |
429 #define EXP_SECTION_AREA 0.0039496598225815571993 | |
430 | |
14655
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 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
|
432 #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
|
433 #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
|
434 #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
|
435 #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
|
436 #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
|
437 |
5742 | 438 static ZIGINT ki[ZIGGURAT_TABLE_SIZE]; |
439 static double wi[ZIGGURAT_TABLE_SIZE], fi[ZIGGURAT_TABLE_SIZE]; | |
440 static ZIGINT ke[ZIGGURAT_TABLE_SIZE]; | |
441 static double we[ZIGGURAT_TABLE_SIZE], fe[ZIGGURAT_TABLE_SIZE]; | |
442 | |
443 /* | |
444 This code is based on the paper Marsaglia and Tsang, "The ziggurat method | |
445 for generating random variables", Journ. Statistical Software. Code was | |
446 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
|
447 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
|
448 Mersenne Twister as the integer generator and uses 256 levels in the |
5742 | 449 Ziggurat. This has several advantages. |
450 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
451 1) As Marsaglia and Tsang themselves states, the more levels the few |
5742 | 452 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
|
453 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
|
454 generator, thus the use of a Mersenne Twister for the core random |
5742 | 455 generator makes this cycle extremely long. |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
456 3) The license on the original code was unclear, thus rewriting the code |
5742 | 457 from the article means we are free of copyright issues. |
458 4) Compile flag for full 53-bit random mantissa. | |
459 | |
460 It should be stated that the authors made my life easier, by the fact that | |
461 the algorithm developed in the text of the article is for a 256 level | |
462 ziggurat, even if the code itself isn't... | |
463 | |
464 One modification to the algorithm developed in the article, is that it is | |
465 assumed that 0 <= x < Inf, and "unsigned long"s are used, thus resulting in | |
466 terms like 2^32 in the code. As the normal distribution is defined between | |
467 -Inf < x < Inf, we effectively only have 31 bit integers plus a sign. Thus | |
468 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
|
469 this term. The exponential distribution is one sided so we use the |
5742 | 470 full 32 bits. We use EMANTISSA for this term. |
471 | |
472 It appears that I'm slightly slower than the code in the article, this | |
473 is partially due to a better generator of random integers than they | |
474 use. But might also be that the case of rapid return was optimized by | |
475 inlining the relevant code with a #define. As the basic Mersenne | |
476 Twister is only 25% faster than this code I suspect that the main | |
477 reason is just the use of the Mersenne Twister and not the inlining, | |
478 so I'm not going to try and optimize further. | |
479 */ | |
480 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
481 static void |
5742 | 482 create_ziggurat_tables (void) |
483 { | |
484 int i; | |
485 double x, x1; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
486 |
5742 | 487 /* Ziggurat tables for the normal distribution */ |
488 x1 = ZIGGURAT_NOR_R; | |
489 wi[255] = x1 / NMANTISSA; | |
490 fi[255] = exp (-0.5 * x1 * x1); | |
491 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
492 /* 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
|
493 * defines this as |
5742 | 494 * 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
|
495 * where v is the area of each strip of the ziggurat. |
5742 | 496 */ |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
497 ki[0] = static_cast<ZIGINT> (x1 * fi[255] / NOR_SECTION_AREA * NMANTISSA); |
5742 | 498 wi[0] = NOR_SECTION_AREA / fi[255] / NMANTISSA; |
499 fi[0] = 1.; | |
500 | |
501 for (i = 254; i > 0; i--) | |
502 { | |
503 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus | |
504 * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y)) | |
505 */ | |
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
506 x = 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
|
507 ki[i+1] = static_cast<ZIGINT> (x / x1 * NMANTISSA); |
5742 | 508 wi[i] = x / NMANTISSA; |
509 fi[i] = exp (-0.5 * x * x); | |
510 x1 = x; | |
511 } | |
512 | |
513 ki[1] = 0; | |
514 | |
515 /* Zigurrat tables for the exponential distribution */ | |
516 x1 = ZIGGURAT_EXP_R; | |
517 we[255] = x1 / EMANTISSA; | |
518 fe[255] = exp (-x1); | |
519 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
520 /* 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
|
521 * defines this as |
5742 | 522 * 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
|
523 * where v is the area of each strip of the ziggurat. |
5742 | 524 */ |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
525 ke[0] = static_cast<ZIGINT> (x1 * fe[255] / EXP_SECTION_AREA * EMANTISSA); |
5742 | 526 we[0] = EXP_SECTION_AREA / fe[255] / EMANTISSA; |
527 fe[0] = 1.; | |
528 | |
529 for (i = 254; i > 0; i--) | |
530 { | |
531 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus | |
532 * need inverse operator of y = exp(-x) -> x = -ln(y) | |
533 */ | |
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
534 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
|
535 ke[i+1] = static_cast<ZIGINT> (x / x1 * EMANTISSA); |
5742 | 536 we[i] = x / EMANTISSA; |
537 fe[i] = exp (-x); | |
538 x1 = x; | |
539 } | |
540 ke[1] = 0; | |
541 | |
542 initt = 0; | |
543 } | |
544 | |
545 /* | |
546 * Here is the guts of the algorithm. As Marsaglia and Tsang state the | |
547 * algorithm in their paper | |
548 * | |
549 * 1) Calculate a random signed integer j and let i be the index | |
550 * provided by the rightmost 8-bits of j | |
551 * 2) Set x = j * w_i. If j < k_i return x | |
552 * 3) If i = 0, then return x from the tail | |
553 * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x | |
554 * 5) goto step 1 | |
555 * | |
556 * Where f is the functional form of the distribution, which for a normal | |
557 * distribution is exp(-0.5*x*x) | |
558 */ | |
559 | |
560 double | |
561 oct_randn (void) | |
562 { | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
563 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
|
564 create_ziggurat_tables (); |
5742 | 565 |
566 while (1) | |
567 { | |
568 /* The following code is specialized for 32-bit mantissa. | |
569 * Compared to the arbitrary mantissa code, there is a performance | |
570 * gain for 32-bits: PPC: 2%, MIPS: 8%, x86: 40% | |
571 * There is a bigger performance gain compared to using a full | |
572 * 53-bit mantissa: PPC: 60%, MIPS: 65%, x86: 240% | |
573 * Of course, different compilers and operating systems may | |
574 * have something to do with this. | |
575 */ | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21720
diff
changeset
|
576 # if defined (HAVE_X86_32) |
5742 | 577 /* 53-bit mantissa, 1-bit sign, x86 32-bit architecture */ |
578 double x; | |
579 int si,idx; | |
21332
db574711df78
eliminate deprecated register storage class specifier
John W. Eaton <jwe@octave.org>
parents:
21308
diff
changeset
|
580 uint32_t lo, hi; |
5742 | 581 int64_t rabs; |
582 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
|
583 lo = randi32 (); |
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
584 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
|
585 hi = randi32 (); |
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
586 si = hi & UMASK; |
5742 | 587 p[0] = lo; |
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
588 p[1] = hi & 0x1FFFFF; |
5742 | 589 x = ( si ? -rabs : rabs ) * wi[idx]; |
20955
77f5591878bf
maint: Use '! expr' rather than '!expr' to conform to coding guidelines.
Rik <rik@octave.org>
parents:
20791
diff
changeset
|
590 # else /* ! HAVE_X86_32 */ |
5742 | 591 /* arbitrary mantissa (selected by NRANDI, with 1 bit for sign) */ |
592 const uint64_t r = NRANDI; | |
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
593 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
|
594 const int idx = static_cast<int> (rabs & 0xFF); |
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
595 const double x = ( (r & 1) ? -rabs : rabs) * wi[idx]; |
20955
77f5591878bf
maint: Use '! expr' rather than '!expr' to conform to coding guidelines.
Rik <rik@octave.org>
parents:
20791
diff
changeset
|
596 # endif /* ! HAVE_X86_32 */ |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
597 if (rabs < static_cast<int64_t> (ki[idx])) |
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
598 return x; /* 99.3% of the time we return here 1st try */ |
5742 | 599 else if (idx == 0) |
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
600 { |
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
601 /* As stated in Marsaglia and Tsang |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
602 * |
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
603 * 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
|
604 * 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
|
605 * 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
|
606 * 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
|
607 * 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
|
608 * |
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
609 * [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
|
610 */ |
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
611 double xx, yy; |
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
612 do |
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
613 { |
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
614 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
|
615 yy = - std::log (RANDU); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
616 } |
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
617 while ( yy+yy <= xx*xx); |
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
618 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
|
619 } |
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
620 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
|
621 return x; |
5742 | 622 } |
623 } | |
624 | |
625 double | |
626 oct_rande (void) | |
627 { | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
628 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
|
629 create_ziggurat_tables (); |
5742 | 630 |
631 while (1) | |
632 { | |
633 ZIGINT ri = ERANDI; | |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
634 const int idx = static_cast<int> (ri & 0xFF); |
5742 | 635 const double x = ri * we[idx]; |
636 if (ri < ke[idx]) | |
15970
ca6202597201
Use C-style comments in C code.
Mike Miller <mtmiller@ieee.org>
parents:
15408
diff
changeset
|
637 return x; /* 98.9% of the time we return here 1st try */ |
5742 | 638 else if (idx == 0) |
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
639 { |
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
640 /* As stated in Marsaglia and Tsang |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
641 * |
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
642 * For the exponential tail, the method of Marsaglia[5] provides: |
5742 | 643 * x = r - ln(U); |
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
644 */ |
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
645 return ZIGGURAT_EXP_R - std::log (RANDU); |
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
646 } |
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
647 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
|
648 return x; |
5742 | 649 } |
650 } | |
651 | |
14655
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 ZIGINT |
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 EMANTISSA |
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 ERANDI |
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 NMANTISSA |
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 NRANDI |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
657 #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
|
658 |
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 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
|
660 #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
|
661 #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
|
662 #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
|
663 #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
|
664 #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
|
665 |
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 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
|
667 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
|
668 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
|
669 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
|
670 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
671 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
|
672 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
|
673 { |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
674 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
|
675 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
|
676 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
677 /* 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
|
678 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
|
679 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
|
680 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
|
681 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
682 /* 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
|
683 * 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
|
684 * 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
|
685 * 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
|
686 */ |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
687 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
|
688 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
|
689 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
|
690 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
691 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
|
692 { |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
693 /* 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
|
694 * 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
|
695 */ |
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
696 x = 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
|
697 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
|
698 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
|
699 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
|
700 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
|
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 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
703 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
|
704 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
705 /* 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
|
706 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
|
707 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
|
708 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
|
709 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
710 /* 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
|
711 * 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
|
712 * 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
|
713 * 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
|
714 */ |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
715 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
|
716 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
|
717 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
|
718 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
719 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
|
720 { |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
721 /* 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
|
722 * 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
|
723 */ |
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
724 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
|
725 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
|
726 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
|
727 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
|
728 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
|
729 } |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
730 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
|
731 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
732 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
|
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 /* |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
736 * 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
|
737 * 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
|
738 * |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
739 * 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
|
740 * 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
|
741 * 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
|
742 * 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
|
743 * 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
|
744 * 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
|
745 * |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
746 * 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
|
747 * 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
|
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 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
750 float |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
751 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
|
752 { |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
753 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
|
754 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
|
755 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
756 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
|
757 { |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
758 /* 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
|
759 const uint32_t r = randi32 (); |
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
760 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
|
761 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
|
762 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
|
763 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
|
764 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
|
765 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
|
766 { |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
767 /* 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
|
768 * |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
769 * 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
|
770 * 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
|
771 * 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
|
772 * 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
|
773 * 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
|
774 * |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
775 * [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
|
776 */ |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
777 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
|
778 do |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
779 { |
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
780 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
|
781 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
|
782 } |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
783 while ( yy+yy <= xx*xx); |
21661
1b9a36a66b01
maint: Octave coding convention cleanups.
Rik <rik@octave.org>
parents:
21332
diff
changeset
|
784 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
|
785 } |
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
786 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
|
787 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
|
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 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
791 float |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
792 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
|
793 { |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
794 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
|
795 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
|
796 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
797 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
|
798 { |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
799 ZIGINT ri = ERANDI; |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
800 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
|
801 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
|
802 if (ri < fke[idx]) |
15970
ca6202597201
Use C-style comments in C code.
Mike Miller <mtmiller@ieee.org>
parents:
15408
diff
changeset
|
803 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
|
804 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
|
805 { |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
806 /* 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
|
807 * |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
808 * 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
|
809 * 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
|
810 */ |
21942
aab79a1885cc
limit gnulib headers to liboctave/wrappers directory
John W. Eaton <jwe@octave.org>
parents:
21928
diff
changeset
|
811 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
|
812 } |
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
813 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
|
814 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
|
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 } |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
817 |
5742 | 818 /* Array generators */ |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
819 void |
5742 | 820 oct_fill_randu (octave_idx_type n, double *p) |
821 { | |
822 octave_idx_type i; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
823 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
|
824 p[i] = oct_randu (); |
5742 | 825 } |
826 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
827 void |
5742 | 828 oct_fill_randn (octave_idx_type n, double *p) |
829 { | |
830 octave_idx_type i; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
831 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
|
832 p[i] = oct_randn (); |
5742 | 833 } |
834 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
835 void |
5742 | 836 oct_fill_rande (octave_idx_type n, double *p) |
837 { | |
838 octave_idx_type i; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
839 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
|
840 p[i] = oct_rande (); |
5742 | 841 } |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
842 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
843 void |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
844 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
|
845 { |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
846 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
|
847 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
|
848 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
|
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 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
851 void |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
852 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
|
853 { |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
854 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
|
855 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
|
856 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
|
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 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
859 void |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
860 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
|
861 { |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
862 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
|
863 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
|
864 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
|
865 } |