Mercurial > octave

annotate liboctave/numeric/randmtzig.cc @ 30394:f3f3e3793fb5 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: style check C++ files in liboctave/ ahead of 7.1 release. * Array.cc, Array.h, CMatrix.cc, CSparse.cc, MArray.cc, MSparse.cc, MatrixType.cc, MatrixType.h, Range.cc, Range.h, Sparse.cc, Sparse.h, dMatrix.cc, dSparse.cc, fCMatrix.cc, fCMatrix.h, fMatrix.cc, fMatrix.h, idx-vector.cc, idx-vector.h, CollocWt.cc, CollocWt.h, DASPK.cc, DASRT.cc, DASSL.cc, EIG.cc, LSODE.cc, Quad.cc, aepbalance.cc, chol.cc, eigs-base.cc, gepbalance.cc, gsvd.cc, lo-amos-proto.h, lo-arpack-proto.h, lo-blas-proto.h, lo-lapack-proto.h, lo-qrupdate-proto.h, lo-specfun.cc, oct-convn.h, oct-fftw.h, oct-norm.cc, oct-spparms.cc, randmtzig.cc, sparse-chol.cc, sparse-dmsolve.cc, sparse-lu.cc, sparse-qr.cc, sparse-qr.h, svd.cc, file-ops.cc, file-stat.cc, lo-sysinfo.cc, oct-env.cc, f77-fcn.h, lo-cutils.h, lo-ieee.h, lo-regexp.cc, lo-utils.cc, oct-binmap.h, oct-inttypes.h, oct-sparse.cc, oct-sparse.h, oct-string.cc, url-transfer.cc: Style check C++ files in liboctave/ ahead of 7.1 release.
author Rik <rik@octave.org>
date Mon, 29 Nov 2021 10:03:47 -0800
parents 7854d5752dd2
children 796f54d4ddbf
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
1 ////////////////////////////////////////////////////////////////////////
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
2 //
29358
0a5b15007766 update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 27933
diff changeset
3 // Copyright (C) 2006-2021 The Octave Project Developers
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
4 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
5 // See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
6 // distribution or <https://octave.org/copyright/>.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
7 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
8 // This file is part of Octave.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
9 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
10 // Octave is free software: you can redistribute it and/or modify it
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
11 // under the terms of the GNU General Public License as published by
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
12 // the Free Software Foundation, either version 3 of the License, or
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
13 // (at your option) any later version.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
14 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
15 // Octave is distributed in the hope that it will be useful, but
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
18 // GNU General Public License for more details.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
19 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
20 // You should have received a copy of the GNU General Public License
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
21 // along with Octave; see the file COPYING. If not, see
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
22 // <https://www.gnu.org/licenses/>.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
23 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
24 ////////////////////////////////////////////////////////////////////////
7019
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6959
diff changeset
25
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
26 /*
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
27 A C-program for MT19937, with initialization improved 2002/2/10.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
28 Coded by Takuji Nishimura and Makoto Matsumoto.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
29 This is a faster version by taking Shawn Cokus's optimization,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
30 Matthe Bellew's simplification, Isaku Wada's real version.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
31 David Bateman added normal and exponential distributions following
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
32 Marsaglia and Tang's Ziggurat algorithm.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
33
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
34 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
35 Copyright (C) 2004, David Bateman
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
36 All rights reserved.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
37
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
38 Redistribution and use in source and binary forms, with or without
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
39 modification, are permitted provided that the following conditions
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
40 are met:
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
41
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
42 1. Redistributions of source code must retain the above copyright
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
43 notice, this list of conditions and the following disclaimer.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
44
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
45 2. Redistributions in binary form must reproduce the above copyright
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
46 notice, this list of conditions and the following disclaimer in the
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
47 documentation and/or other materials provided with the distribution.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
48
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
49 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
50 products derived from this software without specific prior written
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
51 permission.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
52
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
53 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
54 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
55 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
56 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
57 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
58 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
59 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
60 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
61 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
62 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
63 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
64
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
65
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
66 Any feedback is very welcome.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
67 http://www.math.keio.ac.jp/matumoto/emt.html
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
68 email: matumoto@math.keio.ac.jp
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
69
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
70 * 2004-01-19 Paul Kienzle
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
71 * * comment out main
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
72 * add init_by_entropy, get_state, set_state
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
73 * * converted to allow compiling by C++ compiler
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
74 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
75 * 2004-01-25 David Bateman
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
76 * * Add Marsaglia and Tsang Ziggurat code
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
77 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
78 * 2004-07-13 Paul Kienzle
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
79 * * make into an independent library with some docs.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
80 * * introduce new main and test code.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
81 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
82 * 2004-07-28 Paul Kienzle & David Bateman
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
83 * * add -DALLBITS flag for 32 vs. 53 bits of randomness in mantissa
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
84 * * make the naming scheme more uniform
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
85 * * add -DHAVE_X86 for faster support of 53 bit mantissa on x86 arch.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
86 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
87 * 2005-02-23 Paul Kienzle
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
88 * * 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
89 *
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
90 * 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
91 * * 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
92 * 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
93 * * 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
94 *
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
95 * 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
96 * * 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
97 * * Add the single precision generators
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
98 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
99
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
100 /*
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
101 === Build instructions ===
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
102
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
103 Compile with -DHAVE_GETTIMEOFDAY if the gettimeofday function is
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
104 available. This is not necessary if your architecture has
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
105 /dev/urandom defined.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
106
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
107 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
108 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
109 -DUSE_X86_32=0. You should also consider -march=i686 or similar for
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
110 extra performance. Check whether -DUSE_X86_32=0 is faster on 64-bit
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
111 x86 architectures.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
112
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
113 If you want to replace the Mersenne Twister with another
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
114 generator then redefine randi32 appropriately.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
115
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
116 === Usage instructions ===
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
117 Before using any of the generators, initialize the state with one of
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
118 the init_mersenne_twister functions.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
119
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
120 All generators share the same state vector.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
121
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
122 === Mersenne Twister ===
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
123 random initial state:
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
124 void init_mersenne_twister (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
125
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
126 // 32-bit initial state:
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
127 void init_mersenne_twister (uint32_t s)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
128
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
129 // m*32-bit initial state:
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
130 void init_mersenne_twister (uint32_t k[],int m)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
131
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
132 // saves state in array:
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
133 void get_mersenne_twister_state (uint32_t save[MT_N+1])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
134
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
135 // restores state from array
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
136 void set_mersenne_twister_state (uint32_t save[MT_N+1])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
137
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14846
diff changeset
138 static uint32_t randmt (void) returns 32-bit unsigned int
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
139
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
140 === inline generators ===
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14846
diff changeset
141 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
142 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
143 static uint64_t randi54 (void) returns 54-bit unsigned int
27855
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
144 static float randu24 (void) returns 24-bit uniform in (0,1)
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14846
diff changeset
145 static double randu53 (void) returns 53-bit uniform in (0,1)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
146
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
147 double rand_uniform (void) returns M-bit uniform in (0,1)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
148 double rand_normal (void) returns M-bit standard normal
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
149 double rand_exponential (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
150
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
151 === Array generators ===
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
152 void rand_uniform (octave_idx_type, double [])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
153 void rand_normal (octave_idx_type, double [])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
154 void rand_exponential (octave_idx_type, double [])
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
155 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
156
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
157 #if defined (HAVE_CONFIG_H)
21301
40de9f8f23a6 Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents: 21213
diff changeset
158 # include "config.h"
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
159 #endif
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
160
23662
bd77ab816e43 eliminate obsolete file lo-math.h
John W. Eaton <jwe@octave.org>
parents: 23649
diff changeset
161 #include <cmath>
21720
b28c26ae737c use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents: 21661
diff changeset
162 #include <cstdio>
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
163
24856
8bb0251fcfde Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents: 24534
diff changeset
164 #include <algorithm>
28597
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
165 #include <random>
24856
8bb0251fcfde Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents: 24534
diff changeset
166
28597
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
167 #include "oct-syscalls.h"
21928
315f4ba604c8 * randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
168 #include "oct-time.h"
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
169 #include "randmtzig.h"
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
170
21066
258c787cd9ce maint: Use "FIXME:" consistently in code base.
Rik <rik@octave.org>
parents: 20955
diff changeset
171 /* 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
172 #if ! defined (USE_X86_32)
21202
f7121e111991 maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents: 21066
diff changeset
173 # 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
174 # define USE_X86_32 1
f7121e111991 maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents: 21066
diff changeset
175 # else
f7121e111991 maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents: 21066
diff changeset
176 # define USE_X86_32 0
f7121e111991 maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents: 21066
diff changeset
177 # endif
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
178 #endif
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
179
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
180 namespace octave
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
181 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
182 /* ===== Mersenne Twister 32-bit generator ===== */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
183
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
184 #define MT_M 397
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
185 #define MATRIX_A 0x9908b0dfUL /* constant vector a */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
186 #define UMASK 0x80000000UL /* most significant w-r bits */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
187 #define LMASK 0x7fffffffUL /* least significant r bits */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
188 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
189 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
190
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
191 static uint32_t *next;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
192 static uint32_t state[MT_N]; /* the array for the state vector */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
193 static int left = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
194 static int initf = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
195 static int initt = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
196 static int inittf = 1;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
197
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
198 /* initializes state[MT_N] with a seed */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
199 void init_mersenne_twister (const uint32_t s)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
200 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
201 int j;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
202 state[0] = s & 0xffffffffUL;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
203 for (j = 1; j < MT_N; j++)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
204 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
205 state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
206 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
207 /* In the previous versions, MSBs of the seed affect */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
208 /* only MSBs of the array state[]. */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
209 /* 2002/01/09 modified by Makoto Matsumoto */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
210 state[j] &= 0xffffffffUL; /* for >32 bit machines */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
211 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
212 left = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
213 initf = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
214 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
215
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
216 /* initialize by an array with array-length */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
217 /* init_key is the array for initializing keys */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
218 /* key_length is its length */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
219 void init_mersenne_twister (const uint32_t *init_key, const int key_length)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
220 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
221 int i, j, k;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
222 init_mersenne_twister (19650218UL);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
223 i = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
224 j = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
225 k = (MT_N > key_length ? MT_N : key_length);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
226 for (; k; k--)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
227 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
228 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL))
27933
863ae57eee69 maint: Use Octave coding conventions in liboctave/
Rik <rik@octave.org>
parents: 27923
diff changeset
229 + init_key[j] + j; /* non linear */
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
230 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
231 i++;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
232 j++;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
233 if (i >= MT_N)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
234 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
235 state[0] = state[MT_N-1];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
236 i = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
237 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
238 if (j >= key_length)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
239 j = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
240 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
241 for (k = MT_N - 1; k; k--)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
242 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
243 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL))
27933
863ae57eee69 maint: Use Octave coding conventions in liboctave/
Rik <rik@octave.org>
parents: 27923
diff changeset
244 - i; /* non linear */
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
245 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
246 i++;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
247 if (i >= MT_N)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
248 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
249 state[0] = state[MT_N-1];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
250 i = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
251 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
252 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
253
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
254 state[0] = 0x80000000UL; /* MSB is 1; assuring nonzero initial array */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
255 left = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
256 initf = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
257 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
258
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
259 void init_mersenne_twister (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
260 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
261 uint32_t entropy[MT_N];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
262 int n = 0;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
263
28597
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
264 // Gather some entropy from various sources
21928
315f4ba604c8 * randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
265
27102
84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
266 sys::time now;
21928
315f4ba604c8 * randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
267
28597
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
268 // Current time in seconds
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
269 if (n < MT_N)
28597
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
270 entropy[n++] = now.unix_time ();
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
271
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
272 // CPU time used (usec)
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
273 if (n < MT_N)
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
274 entropy[n++] = clock ();
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
275
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
276 // Fractional part of current time
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
277 if (n < MT_N)
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
278 entropy[n++] = now.usec ();
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
279
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
280 // Include the PID to make sure that several processes reaching here at the
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
281 // same time use different random numbers.
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
282 if (n < MT_N)
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
283 entropy[n++] = sys::getpid ();
21928
315f4ba604c8 * randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
284
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
285 if (n < MT_N)
28597
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
286 {
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
287 try
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
288 {
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
289 // The standard doesn't *guarantee* that random_device provides
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
290 // non-deterministic random numbers. So add entropy from this
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
291 // source last to make sure we gathered at least some entropy from
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
292 // the earlier sources.
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
293 std::random_device rd;
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
294 std::uniform_int_distribution<uint32_t> dist;
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
295 // Add 1024 bit of "true" entropy
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
296 int n_max = std::min (n + 32, MT_N);
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
297 while (n < n_max)
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
298 entropy[n++] = dist (rd);
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
299 }
29163
8f67ad8b3103 maint: Updating naming conventions for exceptions and use const where possible.
Rik <rik@octave.org>
parents: 28597
diff changeset
300 catch (const std::exception&)
28597
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
301 {
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
302 // Just ignore any exception and skip that source of entropy.
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
303 }
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
304 }
21928
315f4ba604c8 * randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
305
28597
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
306 // Send all the entropy into the initial state vector
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
307 init_mersenne_twister (entropy, n);
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
308 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
309
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
310 void set_mersenne_twister_state (const uint32_t *save)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
311 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
312 std::copy_n (save, MT_N, state);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
313 left = save[MT_N];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
314 next = state + (MT_N - left + 1);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
315 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
316
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
317 void get_mersenne_twister_state (uint32_t *save)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
318 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
319 std::copy_n (state, MT_N, save);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
320 save[MT_N] = left;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
321 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
322
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
323 static void next_state (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
324 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
325 uint32_t *p = state;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
326 int j;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
327
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
328 /* if init_by_int() has not been called, */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
329 /* a default initial seed is used */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
330 /* if (initf==0) init_by_int(5489UL); */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
331 /* Or better yet, a random seed! */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
332 if (initf == 0)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
333 init_mersenne_twister ();
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
334
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
335 left = MT_N;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
336 next = state;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
337
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
338 for (j = MT_N - MT_M + 1; --j; p++)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
339 *p = p[MT_M] ^ TWIST(p[0], p[1]);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
340
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
341 for (j = MT_M; --j; p++)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
342 *p = p[MT_M-MT_N] ^ TWIST(p[0], p[1]);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
343
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
344 *p = p[MT_M-MT_N] ^ TWIST(p[0], state[0]);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
345 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
346
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
347 /* generates a random number on [0,0xffffffff]-interval */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
348 static uint32_t randmt (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
349 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
350 uint32_t y;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
351
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
352 if (--left == 0)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
353 next_state ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
354 y = *next++;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
355
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
356 /* Tempering */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
357 y ^= (y >> 11);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
358 y ^= (y << 7) & 0x9d2c5680UL;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
359 y ^= (y << 15) & 0xefc60000UL;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
360 return (y ^ (y >> 18));
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
361 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
362
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
363 /* ===== Uniform generators ===== */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
364
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
365 /* Select which 32 bit generator to use */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
366 #define randi32 randmt
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
367
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
368 static uint64_t randi53 (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
369 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
370 const uint32_t lo = randi32 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
371 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
372 #if defined (HAVE_X86_32)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
373 uint64_t u;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
374 uint32_t *p = (uint32_t *)&u;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
375 p[0] = lo;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
376 p[1] = hi;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
377 return u;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
378 #else
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
379 return ((static_cast<uint64_t> (hi) << 32) | lo);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
380 #endif
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
381 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
382
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
383 static uint64_t randi54 (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
384 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
385 const uint32_t lo = randi32 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
386 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
387 #if defined (HAVE_X86_32)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
388 uint64_t u;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
389 uint32_t *p = static_cast<uint32_t *> (&u);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
390 p[0] = lo;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
391 p[1] = hi;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
392 return u;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
393 #else
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
394 return ((static_cast<uint64_t> (hi) << 32) | lo);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
395 #endif
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
396 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
397
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
398 /* generates a random number on (0,1)-real-interval */
27855
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
399 static float randu24 (void)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
400 {
27855
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
401 uint32_t i;
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
402
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
403 do
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
404 {
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
405 i = randi32 () & static_cast<uint32_t> (0xFFFFFF);
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
406 }
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
407 while (i == 0);
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
408
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
409 return i * (1.0f / 16777216.0f);
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
410 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
411
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
412 /* generates a random number on (0,1) with 53-bit resolution */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
413 static double randu53 (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
414 {
27858
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
415 int32_t a, b;
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
416
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
417 do
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
418 {
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
419 a = randi32 () >> 5;
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
420 b = randi32 () >> 6;
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
421 }
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
422 while (a == 0 && b == 0);
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
423
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
424 return (a*67108864.0 + b) * (1.0/9007199254740992.0);
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
425 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
426
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
427 /* Determine mantissa for uniform doubles */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
428 template <>
29228
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents: 29163
diff changeset
429 OCTAVE_API double
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
430 rand_uniform<double> (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
431 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
432 return randu53 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
433 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
434
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
435 /* Determine mantissa for uniform floats */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
436 template <>
29228
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents: 29163
diff changeset
437 OCTAVE_API float
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
438 rand_uniform<float> (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
439 {
27855
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
440 return randu24 ();
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
441 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
442
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
443 /* ===== Ziggurat normal and exponential generators ===== */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
444
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
445 #define ZIGGURAT_TABLE_SIZE 256
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
446
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
447 #define ZIGGURAT_NOR_R 3.6541528853610088
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
448 #define ZIGGURAT_NOR_INV_R 0.27366123732975828
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
449 #define NOR_SECTION_AREA 0.00492867323399
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
450
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
451 #define ZIGGURAT_EXP_R 7.69711747013104972
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
452 #define ZIGGURAT_EXP_INV_R 0.129918765548341586
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
453 #define EXP_SECTION_AREA 0.0039496598225815571993
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
454
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
455 #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
456 #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
457 #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
458 #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
459 #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
460 #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
461
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
462 static ZIGINT ki[ZIGGURAT_TABLE_SIZE];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
463 static double wi[ZIGGURAT_TABLE_SIZE], fi[ZIGGURAT_TABLE_SIZE];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
464 static ZIGINT ke[ZIGGURAT_TABLE_SIZE];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
465 static double we[ZIGGURAT_TABLE_SIZE], fe[ZIGGURAT_TABLE_SIZE];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
466
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
467 /*
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
468 This code is based on the paper Marsaglia and Tsang, "The ziggurat method
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
469 for generating random variables", Journ. Statistical Software. Code was
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
470 presented in this paper for a Ziggurat of 127 levels and using a 32 bit
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
471 integer random number generator. This version of the code, uses the
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
472 Mersenne Twister as the integer generator and uses 256 levels in the
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
473 Ziggurat. This has several advantages.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
474
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
475 1) As Marsaglia and Tsang themselves states, the more levels the few
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
476 times the expensive tail algorithm must be called
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
477 2) The cycle time of the generator is determined by the integer
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
478 generator, thus the use of a Mersenne Twister for the core random
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
479 generator makes this cycle extremely long.
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
480 3) The license on the original code was unclear, thus rewriting the code
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
481 from the article means we are free of copyright issues.
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
482 4) Compile flag for full 53-bit random mantissa.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
483
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
484 It should be stated that the authors made my life easier, by the fact that
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
485 the algorithm developed in the text of the article is for a 256 level
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
486 ziggurat, even if the code itself isn't...
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
487
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
488 One modification to the algorithm developed in the article, is that it is
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
489 assumed that 0 <= x < Inf, and "unsigned long"s are used, thus resulting in
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
490 terms like 2^32 in the code. As the normal distribution is defined between
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
491 -Inf < x < Inf, we effectively only have 31 bit integers plus a sign. Thus
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
492 in Marsaglia and Tsang, terms like 2^32 become 2^31. We use NMANTISSA for
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
493 this term. The exponential distribution is one sided so we use the
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
494 full 32 bits. We use EMANTISSA for this term.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
495
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
496 It appears that I'm slightly slower than the code in the article, this
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
497 is partially due to a better generator of random integers than they
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
498 use. But might also be that the case of rapid return was optimized by
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
499 inlining the relevant code with a #define. As the basic Mersenne
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
500 Twister is only 25% faster than this code I suspect that the main
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
501 reason is just the use of the Mersenne Twister and not the inlining,
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
502 so I'm not going to try and optimize further.
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
503 */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
504
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
505 void create_ziggurat_tables (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
506 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
507 int i;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
508 double x, x1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
509
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
510 /* Ziggurat tables for the normal distribution */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
511 x1 = ZIGGURAT_NOR_R;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
512 wi[255] = x1 / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
513 fi[255] = exp (-0.5 * x1 * x1);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
514
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
515 /* Index zero is special for tail strip, where Marsaglia and Tsang
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
516 * defines this as
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
517 * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1,
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
518 * where v is the area of each strip of the ziggurat.
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
519 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
520 ki[0] = static_cast<ZIGINT> (x1 * fi[255] / NOR_SECTION_AREA * NMANTISSA);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
521 wi[0] = NOR_SECTION_AREA / fi[255] / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
522 fi[0] = 1.;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
523
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
524 for (i = 254; i > 0; i--)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
525 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
526 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
527 * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y))
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
528 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
529 x = std::sqrt (-2. * std::log (NOR_SECTION_AREA / x1 + fi[i+1]));
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
530 ki[i+1] = static_cast<ZIGINT> (x / x1 * NMANTISSA);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
531 wi[i] = x / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
532 fi[i] = exp (-0.5 * x * x);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
533 x1 = x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
534 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
535
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
536 ki[1] = 0;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
537
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
538 /* Zigurrat tables for the exponential distribution */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
539 x1 = ZIGGURAT_EXP_R;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
540 we[255] = x1 / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
541 fe[255] = exp (-x1);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
542
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
543 /* Index zero is special for tail strip, where Marsaglia and Tsang
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
544 * defines this as
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
545 * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1,
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
546 * where v is the area of each strip of the ziggurat.
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
547 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
548 ke[0] = static_cast<ZIGINT> (x1 * fe[255] / EXP_SECTION_AREA * EMANTISSA);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
549 we[0] = EXP_SECTION_AREA / fe[255] / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
550 fe[0] = 1.;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
551
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
552 for (i = 254; i > 0; i--)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
553 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
554 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
555 * need inverse operator of y = exp(-x) -> x = -ln(y)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
556 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
557 x = - std::log (EXP_SECTION_AREA / x1 + fe[i+1]);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
558 ke[i+1] = static_cast<ZIGINT> (x / x1 * EMANTISSA);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
559 we[i] = x / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
560 fe[i] = exp (-x);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
561 x1 = x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
562 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
563 ke[1] = 0;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
564
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
565 initt = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
566 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
567
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
568 /*
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
569 * Here is the guts of the algorithm. As Marsaglia and Tsang state the
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
570 * algorithm in their paper
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
571 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
572 * 1) Calculate a random signed integer j and let i be the index
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
573 * provided by the rightmost 8-bits of j
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
574 * 2) Set x = j * w_i. If j < k_i return x
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
575 * 3) If i = 0, then return x from the tail
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
576 * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
577 * 5) goto step 1
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
578 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
579 * Where f is the functional form of the distribution, which for a normal
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
580 * distribution is exp(-0.5*x*x)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
581 */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
582
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
583
29228
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents: 29163
diff changeset
584 template <> OCTAVE_API double rand_normal<double> (void)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
585 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
586 if (initt)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
587 create_ziggurat_tables ();
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
588
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
589 while (1)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
590 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
591 /* The following code is specialized for 32-bit mantissa.
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
592 * Compared to the arbitrary mantissa code, there is a performance
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
593 * gain for 32-bits: PPC: 2%, MIPS: 8%, x86: 40%
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
594 * There is a bigger performance gain compared to using a full
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
595 * 53-bit mantissa: PPC: 60%, MIPS: 65%, x86: 240%
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
596 * Of course, different compilers and operating systems may
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
597 * have something to do with this.
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
598 */
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21720
diff changeset
599 # if defined (HAVE_X86_32)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
600 /* 53-bit mantissa, 1-bit sign, x86 32-bit architecture */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
601 double x;
30394
f3f3e3793fb5 maint: style check C++ files in liboctave/ ahead of 7.1 release.
Rik <rik@octave.org>
parents: 29359
diff changeset
602 int si, idx;
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
603 uint32_t lo, hi;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
604 int64_t rabs;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
605 uint32_t *p = (uint32_t *)&rabs;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
606 lo = randi32 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
607 idx = lo & 0xFF;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
608 hi = randi32 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
609 si = hi & UMASK;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
610 p[0] = lo;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
611 p[1] = hi & 0x1FFFFF;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
612 x = ( si ? -rabs : rabs ) * wi[idx];
22188
1344509a480c style fixes
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
613 # else
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
614 /* arbitrary mantissa (selected by NRANDI, with 1 bit for sign) */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
615 const uint64_t r = NRANDI;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
616 const int64_t rabs = r >> 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
617 const int idx = static_cast<int> (rabs & 0xFF);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
618 const double x = ( (r & 1) ? -rabs : rabs) * wi[idx];
22188
1344509a480c style fixes
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
619 # endif
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
620 if (rabs < static_cast<int64_t> (ki[idx]))
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
621 return x; /* 99.3% of the time we return here 1st try */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
622 else if (idx == 0)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
623 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
624 /* As stated in Marsaglia and Tsang
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
625 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
626 * For the normal tail, the method of Marsaglia[5] provides:
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
627 * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x,
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
628 * then return r+x. Except that r+x is always in the positive
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
629 * tail!!!! Any thing random might be used to determine the
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
630 * sign, but as we already have r we might as well use it
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
631 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
632 * [PAK] but not the bottom 8 bits, since they are all 0 here!
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
633 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
634 double xx, yy;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
635 do
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
636 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
637 xx = - ZIGGURAT_NOR_INV_R * std::log (RANDU);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
638 yy = - std::log (RANDU);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
639 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
640 while ( yy+yy <= xx*xx);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
641 return ((rabs & 0x100) ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
642 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
643 else if ((fi[idx-1] - fi[idx]) * RANDU + fi[idx] < exp (-0.5*x*x))
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
644 return x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
645 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
646 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
647
29228
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents: 29163
diff changeset
648 template <> OCTAVE_API double rand_exponential<double> (void)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
649 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
650 if (initt)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
651 create_ziggurat_tables ();
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
652
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
653 while (1)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
654 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
655 ZIGINT ri = ERANDI;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
656 const int idx = static_cast<int> (ri & 0xFF);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
657 const double x = ri * we[idx];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
658 if (ri < ke[idx])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
659 return x; /* 98.9% of the time we return here 1st try */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
660 else if (idx == 0)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
661 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
662 /* As stated in Marsaglia and Tsang
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
663 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
664 * For the exponential tail, the method of Marsaglia[5] provides:
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
665 * x = r - ln(U);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
666 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
667 return ZIGGURAT_EXP_R - std::log (RANDU);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
668 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
669 else if ((fe[idx-1] - fe[idx]) * RANDU + fe[idx] < exp (-x))
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
670 return x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
671 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
672 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
673
29228
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents: 29163
diff changeset
674 template <> OCTAVE_API void rand_uniform<double> (octave_idx_type n, double *p)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
675 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
676 std::generate_n (p, n, [](void) { return rand_uniform<double> (); });
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
677 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
678
29228
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents: 29163
diff changeset
679 template <> OCTAVE_API void rand_normal (octave_idx_type n, double *p)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
680 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
681 std::generate_n (p, n, [](void) { return rand_normal<double> (); });
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
682 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
683
29228
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents: 29163
diff changeset
684 template <> OCTAVE_API void rand_exponential (octave_idx_type n, double *p)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
685 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
686 std::generate_n (p, n, [](void) { return rand_exponential<double> (); });
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
687 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
688
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
689 #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
690 #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
691 #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
692 #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
693 #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
694 #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
695
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
696 #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
697 #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
698 #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
699 #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
700 #define NRANDI randi32() /* 31 bits for mantissa + 1 bit sign */
27855
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
701 #define RANDU randu24()
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
702
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
703 static ZIGINT fki[ZIGGURAT_TABLE_SIZE];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
704 static float fwi[ZIGGURAT_TABLE_SIZE], ffi[ZIGGURAT_TABLE_SIZE];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
705 static ZIGINT fke[ZIGGURAT_TABLE_SIZE];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
706 static float fwe[ZIGGURAT_TABLE_SIZE], ffe[ZIGGURAT_TABLE_SIZE];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
707
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
708 static void create_ziggurat_float_tables (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
709 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
710 int i;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
711 float x, x1;
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
712
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
713 /* Ziggurat tables for the normal distribution */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
714 x1 = ZIGGURAT_NOR_R;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
715 fwi[255] = x1 / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
716 ffi[255] = exp (-0.5 * x1 * x1);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
717
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
718 /* Index zero is special for tail strip, where Marsaglia and Tsang
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
719 * defines this as
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
720 * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1,
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
721 * where v is the area of each strip of the ziggurat.
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
722 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
723 fki[0] = static_cast<ZIGINT> (x1 * ffi[255] / NOR_SECTION_AREA * NMANTISSA);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
724 fwi[0] = NOR_SECTION_AREA / ffi[255] / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
725 ffi[0] = 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
726
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
727 for (i = 254; i > 0; i--)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
728 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
729 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
730 * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y))
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
731 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
732 x = std::sqrt (-2. * std::log (NOR_SECTION_AREA / x1 + ffi[i+1]));
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
733 fki[i+1] = static_cast<ZIGINT> (x / x1 * NMANTISSA);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
734 fwi[i] = x / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
735 ffi[i] = exp (-0.5 * x * x);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
736 x1 = x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
737 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
738
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
739 fki[1] = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
740
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
741 /* Zigurrat tables for the exponential distribution */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
742 x1 = ZIGGURAT_EXP_R;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
743 fwe[255] = x1 / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
744 ffe[255] = exp (-x1);
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
745
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
746 /* Index zero is special for tail strip, where Marsaglia and Tsang
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
747 * defines this as
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
748 * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1,
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
749 * where v is the area of each strip of the ziggurat.
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
750 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
751 fke[0] = static_cast<ZIGINT> (x1 * ffe[255] / EXP_SECTION_AREA * EMANTISSA);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
752 fwe[0] = EXP_SECTION_AREA / ffe[255] / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
753 ffe[0] = 1.;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
754
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
755 for (i = 254; i > 0; i--)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
756 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
757 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
758 * need inverse operator of y = exp(-x) -> x = -ln(y)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
759 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
760 x = - std::log (EXP_SECTION_AREA / x1 + ffe[i+1]);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
761 fke[i+1] = static_cast<ZIGINT> (x / x1 * EMANTISSA);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
762 fwe[i] = x / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
763 ffe[i] = exp (-x);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
764 x1 = x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
765 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
766 fke[1] = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
767
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
768 inittf = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
769 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
770
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
771 /*
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
772 * Here is the guts of the algorithm. As Marsaglia and Tsang state the
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
773 * algorithm in their paper
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
774 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
775 * 1) Calculate a random signed integer j and let i be the index
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
776 * provided by the rightmost 8-bits of j
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
777 * 2) Set x = j * w_i. If j < k_i return x
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
778 * 3) If i = 0, then return x from the tail
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
779 * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
780 * 5) goto step 1
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
781 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
782 * Where f is the functional form of the distribution, which for a normal
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
783 * distribution is exp(-0.5*x*x)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
784 */
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
29228
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents: 29163
diff changeset
786 template <> OCTAVE_API float rand_normal<float> (void)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
787 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
788 if (inittf)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
789 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
790
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
791 while (1)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
792 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
793 /* 32-bit mantissa */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
794 const uint32_t r = randi32 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
795 const uint32_t rabs = r & LMASK;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
796 const int idx = static_cast<int> (r & 0xFF);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
797 const float x = static_cast<int32_t> (r) * fwi[idx];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
798 if (rabs < fki[idx])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
799 return x; /* 99.3% of the time we return here 1st try */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
800 else if (idx == 0)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
801 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
802 /* As stated in Marsaglia and Tsang
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
803 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
804 * For the normal tail, the method of Marsaglia[5] provides:
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
805 * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x,
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
806 * then return r+x. Except that r+x is always in the positive
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
807 * tail!!!! Any thing random might be used to determine the
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
808 * sign, but as we already have r we might as well use it
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
809 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
810 * [PAK] but not the bottom 8 bits, since they are all 0 here!
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
811 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
812 float xx, yy;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
813 do
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
814 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
815 xx = - ZIGGURAT_NOR_INV_R * std::log (RANDU);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
816 yy = - std::log (RANDU);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
817 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
818 while ( yy+yy <= xx*xx);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
819 return ((rabs & 0x100) ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
820 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
821 else if ((ffi[idx-1] - ffi[idx]) * RANDU + ffi[idx] < exp (-0.5*x*x))
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
822 return x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
823 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
824 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
825
29228
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents: 29163
diff changeset
826 template <> OCTAVE_API float rand_exponential<float> (void)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
827 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
828 if (inittf)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
829 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
830
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
831 while (1)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
832 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
833 ZIGINT ri = ERANDI;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
834 const int idx = static_cast<int> (ri & 0xFF);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
835 const float x = ri * fwe[idx];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
836 if (ri < fke[idx])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
837 return x; /* 98.9% of the time we return here 1st try */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
838 else if (idx == 0)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
839 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
840 /* As stated in Marsaglia and Tsang
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
841 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
842 * For the exponential tail, the method of Marsaglia[5] provides:
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
843 * x = r - ln(U);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
844 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
845 return ZIGGURAT_EXP_R - std::log (RANDU);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
846 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
847 else if ((ffe[idx-1] - ffe[idx]) * RANDU + ffe[idx] < exp (-x))
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
848 return x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
849 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
850 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
851
29228
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents: 29163
diff changeset
852 template <> OCTAVE_API void rand_uniform (octave_idx_type n, float *p)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
853 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
854 std::generate_n (p, n, [](void) { return rand_uniform<float> (); });
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
855 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
856
29228
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents: 29163
diff changeset
857 template <> OCTAVE_API void rand_normal (octave_idx_type n, float *p)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
858 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
859 std::generate_n (p, n, [](void) { return rand_normal<float> (); });
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
860 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
861
29228
5c14f81e0937 Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents: 29163
diff changeset
862 template <> OCTAVE_API void rand_exponential (octave_idx_type n, float *p)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
863 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
864 std::generate_n (p, n, [](void) { return rand_exponential<float> (); });
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
865 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
866 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
867