Mercurial > octave

annotate liboctave/numeric/randmtzig.cc @ 33645:42355b7ec5d7 bytecode-interpreter tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Merge default to bytecode-interpreter.
author Nicholas R. Jankowski <jankowski.nicholas@gmail.com>
date Tue, 04 Jun 2024 11:00:11 -0400
parents 2e484f9f1f18
children
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 //
32632
2e484f9f1f18 maint: update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 32590
diff changeset
3 // Copyright (C) 2006-2024 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:
31771
21f9b34eb893 maint: Eliminate "(void)" in C++ function prototypes/declarations.
Rik <rik@octave.org>
parents: 31706
diff changeset
124 void init_mersenne_twister ()
25433
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
31771
21f9b34eb893 maint: Eliminate "(void)" in C++ function prototypes/declarations.
Rik <rik@octave.org>
parents: 31706
diff changeset
138 static uint32_t randmt () 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 ===
31771
21f9b34eb893 maint: Eliminate "(void)" in C++ function prototypes/declarations.
Rik <rik@octave.org>
parents: 31706
diff changeset
141 static uint32_t randi32 () returns 32-bit unsigned int
21f9b34eb893 maint: Eliminate "(void)" in C++ function prototypes/declarations.
Rik <rik@octave.org>
parents: 31706
diff changeset
142 static uint64_t randi53 () returns 53-bit unsigned int
21f9b34eb893 maint: Eliminate "(void)" in C++ function prototypes/declarations.
Rik <rik@octave.org>
parents: 31706
diff changeset
143 static uint64_t randi54 () returns 54-bit unsigned int
21f9b34eb893 maint: Eliminate "(void)" in C++ function prototypes/declarations.
Rik <rik@octave.org>
parents: 31706
diff changeset
144 static float randu24 () returns 24-bit uniform in (0,1)
21f9b34eb893 maint: Eliminate "(void)" in C++ function prototypes/declarations.
Rik <rik@octave.org>
parents: 31706
diff changeset
145 static double randu53 () returns 53-bit uniform in (0,1)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
146
31771
21f9b34eb893 maint: Eliminate "(void)" in C++ function prototypes/declarations.
Rik <rik@octave.org>
parents: 31706
diff changeset
147 double rand_uniform () returns M-bit uniform in (0,1)
21f9b34eb893 maint: Eliminate "(void)" in C++ function prototypes/declarations.
Rik <rik@octave.org>
parents: 31706
diff changeset
148 double rand_normal () returns M-bit standard normal
21f9b34eb893 maint: Eliminate "(void)" in C++ function prototypes/declarations.
Rik <rik@octave.org>
parents: 31706
diff changeset
149 double rand_exponential () 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>
31142
8245e773bb5b randmtzig.cc: Add missing include <ctime> (bug #62750).
Stefan Brüns <stefan.bruens@rwth-aachen.de>
parents: 30564
diff changeset
163 #include <ctime>
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
164
24856
8bb0251fcfde Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents: 24534
diff changeset
165 #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
166 #include <random>
24856
8bb0251fcfde Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents: 24534
diff changeset
167
28597
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
168 #include "oct-syscalls.h"
21928
315f4ba604c8 * randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
169 #include "oct-time.h"
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
170 #include "randmtzig.h"
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
171
21066
258c787cd9ce maint: Use "FIXME:" consistently in code base.
Rik <rik@octave.org>
parents: 20955
diff changeset
172 /* 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
173 #if ! defined (USE_X86_32)
21202
f7121e111991 maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents: 21066
diff changeset
174 # 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
175 # define USE_X86_32 1
f7121e111991 maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents: 21066
diff changeset
176 # else
f7121e111991 maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents: 21066
diff changeset
177 # define USE_X86_32 0
f7121e111991 maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents: 21066
diff changeset
178 # endif
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
179 #endif
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
180
31605
e88a07dec498 maint: Use macros to begin/end C++ namespaces.
Rik <rik@octave.org>
parents: 31142
diff changeset
181 OCTAVE_BEGIN_NAMESPACE(octave)
e88a07dec498 maint: Use macros to begin/end C++ namespaces.
Rik <rik@octave.org>
parents: 31142
diff changeset
182
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
183 /* ===== Mersenne Twister 32-bit generator ===== */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
184
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
185 #define MT_M 397
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
186 #define MATRIX_A 0x9908b0dfUL /* constant vector a */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
187 #define UMASK 0x80000000UL /* most significant w-r bits */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
188 #define LMASK 0x7fffffffUL /* least significant r bits */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
189 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
190 #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
191
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
192 static uint32_t *next;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
193 static uint32_t state[MT_N]; /* the array for the state vector */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
194 static int left = 1;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
195 static int initf = 0;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
196 static int initt = 1;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
197 static int inittf = 1;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
198
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
199 /* initializes state[MT_N] with a seed */
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
200 void
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
201 init_mersenne_twister (const uint32_t s)
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
202 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
203 int j;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
204 state[0] = s & 0xffffffffUL;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
205 for (j = 1; j < MT_N; j++)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
206 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
207 state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
208 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
209 /* In the previous versions, MSBs of the seed affect */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
210 /* only MSBs of the array state[]. */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
211 /* 2002/01/09 modified by Makoto Matsumoto */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
212 state[j] &= 0xffffffffUL; /* for >32 bit machines */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
213 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
214 left = 1;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
215 initf = 1;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
216 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
217
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
218 /* initialize by an array with array-length */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
219 /* init_key is the array for initializing keys */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
220 /* key_length is its length */
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
221 void
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
222 init_mersenne_twister (const uint32_t *init_key, const int key_length)
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
223 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
224 int i, j, k;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
225 init_mersenne_twister (19650218UL);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
226 i = 1;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
227 j = 0;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
228 k = (MT_N > key_length ? MT_N : key_length);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
229 for (; k; k--)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
230 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
231 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL))
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
232 + init_key[j] + j; /* non linear */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
233 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
234 i++;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
235 j++;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
236 if (i >= MT_N)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
237 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
238 state[0] = state[MT_N-1];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
239 i = 1;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
240 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
241 if (j >= key_length)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
242 j = 0;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
243 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
244 for (k = MT_N - 1; k; k--)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
245 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
246 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL))
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
247 - i; /* non linear */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
248 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
249 i++;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
250 if (i >= MT_N)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
251 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
252 state[0] = state[MT_N-1];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
253 i = 1;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
254 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
255 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
256
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
257 state[0] = 0x80000000UL; /* MSB is 1; assuring nonzero initial array */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
258 left = 1;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
259 initf = 1;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
260 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
261
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
262 void
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
263 init_mersenne_twister ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
264 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
265 uint32_t entropy[MT_N];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
266 int n = 0;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
267
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
268 // Gather some entropy from various sources
21928
315f4ba604c8 * randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
269
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
270 sys::time now;
21928
315f4ba604c8 * randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
271
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
272 // Current time in seconds
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
273 if (n < MT_N)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
274 entropy[n++] = now.unix_time ();
28597
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
275
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
276 // CPU time used (usec)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
277 if (n < MT_N)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
278 entropy[n++] = clock ();
28597
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
279
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
280 // Fractional part of current time
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
281 if (n < MT_N)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
282 entropy[n++] = now.usec ();
28597
a61b651914d6 Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents: 27933
diff changeset
283
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
284 // Include the PID to make sure that several processes reaching here at the
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
285 // same time use different random numbers.
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
286 if (n < MT_N)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
287 entropy[n++] = sys::getpid ();
21928
315f4ba604c8 * randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
288
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
289 if (n < MT_N)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
290 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
291 try
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
292 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
293 // The standard doesn't *guarantee* that random_device provides
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
294 // non-deterministic random numbers. So add entropy from this
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
295 // source last to make sure we gathered at least some entropy from
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
296 // the earlier sources.
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
297 std::random_device rd;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
298 std::uniform_int_distribution<uint32_t> dist;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
299 // Add 1024 bit of "true" entropy
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
300 int n_max = std::min (n + 32, MT_N);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
301 while (n < n_max)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
302 entropy[n++] = dist (rd);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
303 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
304 catch (const std::exception&)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
305 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
306 // Just ignore any exception and skip that source of entropy.
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
307 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
308 }
21928
315f4ba604c8 * randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
309
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
310 // Send all the entropy into the initial state vector
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
311 init_mersenne_twister (entropy, n);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
312 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
313
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
314 void
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
315 set_mersenne_twister_state (const uint32_t *save)
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
316 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
317 std::copy_n (save, MT_N, state);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
318 left = save[MT_N];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
319 next = state + (MT_N - left + 1);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
320 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
321
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
322 void
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
323 get_mersenne_twister_state (uint32_t *save)
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
324 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
325 std::copy_n (state, MT_N, save);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
326 save[MT_N] = left;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
327 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
328
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
329 static void
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
330 next_state ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
331 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
332 uint32_t *p = state;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
333 int j;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
334
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
335 /* if init_by_int() has not been called, */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
336 /* a default initial seed is used */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
337 /* if (initf==0) init_by_int(5489UL); */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
338 /* Or better yet, a random seed! */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
339 if (initf == 0)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
340 init_mersenne_twister ();
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
341
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
342 left = MT_N;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
343 next = state;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
344
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
345 for (j = MT_N - MT_M + 1; --j; p++)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
346 *p = p[MT_M] ^ TWIST(p[0], p[1]);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
347
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
348 for (j = MT_M; --j; p++)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
349 *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
350
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
351 *p = p[MT_M-MT_N] ^ TWIST(p[0], state[0]);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
352 }
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
353
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
354 /* generates a random number on [0,0xffffffff]-interval */
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
355 static uint32_t
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
356 randmt ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
357 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
358 uint32_t y;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
359
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
360 if (--left == 0)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
361 next_state ();
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
362 y = *next++;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
363
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
364 /* Tempering */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
365 y ^= (y >> 11);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
366 y ^= (y << 7) & 0x9d2c5680UL;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
367 y ^= (y << 15) & 0xefc60000UL;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
368 return (y ^ (y >> 18));
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
369 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
370
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
371 /* ===== Uniform generators ===== */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
372
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
373 /* Select which 32 bit generator to use */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
374 #define randi32 randmt
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
375
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
376 static uint64_t
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
377 randi53 ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
378 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
379 const uint32_t lo = randi32 ();
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
380 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
381 #if defined (HAVE_X86_32)
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
382 uint64_t u;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
383 uint32_t *p = (uint32_t *)&u;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
384 p[0] = lo;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
385 p[1] = hi;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
386 return u;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
387 #else
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
388 return ((static_cast<uint64_t> (hi) << 32) | lo);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
389 #endif
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
390 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
391
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
392 static uint64_t
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
393 randi54 ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
394 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
395 const uint32_t lo = randi32 ();
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
396 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
397 #if defined (HAVE_X86_32)
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
398 uint64_t u;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
399 uint32_t *p = static_cast<uint32_t *> (&u);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
400 p[0] = lo;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
401 p[1] = hi;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
402 return u;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
403 #else
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
404 return ((static_cast<uint64_t> (hi) << 32) | lo);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
405 #endif
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
406 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
407
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
408 /* generates a random number on (0,1)-real-interval */
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
409 static float
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
410 randu24 ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
411 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
412 uint32_t i;
27855
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
413
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
414 do
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
415 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
416 i = randi32 () & static_cast<uint32_t> (0xFFFFFF);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
417 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
418 while (i == 0);
27855
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
419
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
420 return i * (1.0f / 16777216.0f);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
421 }
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
422
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
423 /* generates a random number on (0,1) with 53-bit resolution */
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
424 static double
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
425 randu53 ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
426 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
427 int32_t a, b;
27858
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
428
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
429 do
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
430 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
431 a = randi32 () >> 5;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
432 b = randi32 () >> 6;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
433 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
434 while (a == 0 && b == 0);
27858
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
435
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
436 return (a*67108864.0 + b) * (1.0/9007199254740992.0);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
437 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
438
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
439 /* Determine mantissa for uniform doubles */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
440 template <>
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
441 OCTAVE_API double
31771
21f9b34eb893 maint: Eliminate "(void)" in C++ function prototypes/declarations.
Rik <rik@octave.org>
parents: 31706
diff changeset
442 rand_uniform<double> ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
443 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
444 return randu53 ();
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
445 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
446
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
447 /* Determine mantissa for uniform floats */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
448 template <>
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
449 OCTAVE_API float
31771
21f9b34eb893 maint: Eliminate "(void)" in C++ function prototypes/declarations.
Rik <rik@octave.org>
parents: 31706
diff changeset
450 rand_uniform<float> ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
451 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
452 return randu24 ();
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
453 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
454
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
455 /* ===== Ziggurat normal and exponential generators ===== */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
456
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
457 #define ZIGGURAT_TABLE_SIZE 256
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
458
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
459 #define ZIGGURAT_NOR_R 3.6541528853610088
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
460 #define ZIGGURAT_NOR_INV_R 0.27366123732975828
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
461 #define NOR_SECTION_AREA 0.00492867323399
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
462
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
463 #define ZIGGURAT_EXP_R 7.69711747013104972
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
464 #define ZIGGURAT_EXP_INV_R 0.129918765548341586
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
465 #define EXP_SECTION_AREA 0.0039496598225815571993
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
466
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
467 #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
468 #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
469 #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
470 #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
471 #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
472 #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
473
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
474 static ZIGINT ki[ZIGGURAT_TABLE_SIZE];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
475 static double wi[ZIGGURAT_TABLE_SIZE], fi[ZIGGURAT_TABLE_SIZE];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
476 static ZIGINT ke[ZIGGURAT_TABLE_SIZE];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
477 static double we[ZIGGURAT_TABLE_SIZE], fe[ZIGGURAT_TABLE_SIZE];
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
478
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
479 /*
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
480 This code is based on the paper Marsaglia and Tsang, "The ziggurat method
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
481 for generating random variables", Journ. Statistical Software. Code was
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
482 presented in this paper for a Ziggurat of 127 levels and using a 32 bit
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
483 integer random number generator. This version of the code, uses the
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
484 Mersenne Twister as the integer generator and uses 256 levels in the
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
485 Ziggurat. This has several advantages.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
486
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
487 1) As Marsaglia and Tsang themselves states, the more levels the few
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
488 times the expensive tail algorithm must be called
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
489 2) The cycle time of the generator is determined by the integer
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
490 generator, thus the use of a Mersenne Twister for the core random
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
491 generator makes this cycle extremely long.
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
492 3) The license on the original code was unclear, thus rewriting the code
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
493 from the article means we are free of copyright issues.
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
494 4) Compile flag for full 53-bit random mantissa.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
495
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
496 It should be stated that the authors made my life easier, by the fact that
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
497 the algorithm developed in the text of the article is for a 256 level
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
498 ziggurat, even if the code itself isn't...
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
499
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
500 One modification to the algorithm developed in the article, is that it is
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
501 assumed that 0 <= x < Inf, and "unsigned long"s are used, thus resulting in
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
502 terms like 2^32 in the code. As the normal distribution is defined between
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
503 -Inf < x < Inf, we effectively only have 31 bit integers plus a sign. Thus
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
504 in Marsaglia and Tsang, terms like 2^32 become 2^31. We use NMANTISSA for
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
505 this term. The exponential distribution is one sided so we use the
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
506 full 32 bits. We use EMANTISSA for this term.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
507
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
508 It appears that I'm slightly slower than the code in the article, this
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
509 is partially due to a better generator of random integers than they
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
510 use. But might also be that the case of rapid return was optimized by
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
511 inlining the relevant code with a #define. As the basic Mersenne
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
512 Twister is only 25% faster than this code I suspect that the main
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
513 reason is just the use of the Mersenne Twister and not the inlining,
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
514 so I'm not going to try and optimize further.
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
515 */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
516
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
517 void
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
518 create_ziggurat_tables ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
519 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
520 int i;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
521 double x, x1;
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
522
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
523 /* Ziggurat tables for the normal distribution */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
524 x1 = ZIGGURAT_NOR_R;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
525 wi[255] = x1 / NMANTISSA;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
526 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
527
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
528 /* Index zero is special for tail strip, where Marsaglia and Tsang
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
529 * defines this as
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
530 * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1,
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
531 * where v is the area of each strip of the ziggurat.
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
532 */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
533 ki[0] = static_cast<ZIGINT> (x1 * fi[255] / NOR_SECTION_AREA * NMANTISSA);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
534 wi[0] = NOR_SECTION_AREA / fi[255] / NMANTISSA;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
535 fi[0] = 1.;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
536
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
537 for (i = 254; i > 0; i--)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
538 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
539 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
540 * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y))
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
541 */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
542 x = std::sqrt (-2. * std::log (NOR_SECTION_AREA / x1 + fi[i+1]));
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
543 ki[i+1] = static_cast<ZIGINT> (x / x1 * NMANTISSA);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
544 wi[i] = x / NMANTISSA;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
545 fi[i] = exp (-0.5 * x * x);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
546 x1 = x;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
547 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
548
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
549 ki[1] = 0;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
550
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
551 /* Zigurrat tables for the exponential distribution */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
552 x1 = ZIGGURAT_EXP_R;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
553 we[255] = x1 / EMANTISSA;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
554 fe[255] = exp (-x1);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
555
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
556 /* Index zero is special for tail strip, where Marsaglia and Tsang
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
557 * defines this as
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
558 * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1,
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
559 * where v is the area of each strip of the ziggurat.
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
560 */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
561 ke[0] = static_cast<ZIGINT> (x1 * fe[255] / EXP_SECTION_AREA * EMANTISSA);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
562 we[0] = EXP_SECTION_AREA / fe[255] / EMANTISSA;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
563 fe[0] = 1.;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
564
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
565 for (i = 254; i > 0; i--)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
566 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
567 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
568 * need inverse operator of y = exp(-x) -> x = -ln(y)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
569 */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
570 x = - std::log (EXP_SECTION_AREA / x1 + fe[i+1]);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
571 ke[i+1] = static_cast<ZIGINT> (x / x1 * EMANTISSA);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
572 we[i] = x / EMANTISSA;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
573 fe[i] = exp (-x);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
574 x1 = x;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
575 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
576 ke[1] = 0;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
577
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
578 initt = 0;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
579 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
580
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
581 /*
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
582 * Here is the guts of the algorithm. As Marsaglia and Tsang state the
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
583 * algorithm in their paper
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
584 *
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
585 * 1) Calculate a random signed integer j and let i be the index
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
586 * provided by the rightmost 8-bits of j
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
587 * 2) Set x = j * w_i. If j < k_i return x
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
588 * 3) If i = 0, then return x from the tail
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
589 * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
590 * 5) goto step 1
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
591 *
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
592 * Where f is the functional form of the distribution, which for a normal
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
593 * distribution is exp(-0.5*x*x)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
594 */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
595
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
596
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
597 template <>
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
598 OCTAVE_API double
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
599 rand_normal<double> ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
600 {
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
601 if (initt)
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
602 create_ziggurat_tables ();
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
603
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
604 while (1)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
605 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
606 /* The following code is specialized for 32-bit mantissa.
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
607 * Compared to the arbitrary mantissa code, there is a performance
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
608 * gain for 32-bits: PPC: 2%, MIPS: 8%, x86: 40%
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
609 * There is a bigger performance gain compared to using a full
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
610 * 53-bit mantissa: PPC: 60%, MIPS: 65%, x86: 240%
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
611 * Of course, different compilers and operating systems may
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
612 * have something to do with this.
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
613 */
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21720
diff changeset
614 # if defined (HAVE_X86_32)
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
615 /* 53-bit mantissa, 1-bit sign, x86 32-bit architecture */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
616 double x;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
617 int si, idx;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
618 uint32_t lo, hi;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
619 int64_t rabs;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
620 uint32_t *p = (uint32_t *)&rabs;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
621 lo = randi32 ();
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
622 idx = lo & 0xFF;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
623 hi = randi32 ();
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
624 si = hi & UMASK;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
625 p[0] = lo;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
626 p[1] = hi & 0x1FFFFF;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
627 x = ( si ? -rabs : rabs ) * wi[idx];
22188
1344509a480c style fixes
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
628 # else
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
629 /* arbitrary mantissa (selected by NRANDI, with 1 bit for sign) */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
630 const uint64_t r = NRANDI;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
631 const int64_t rabs = r >> 1;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
632 const int idx = static_cast<int> (rabs & 0xFF);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
633 const double x = ( (r & 1) ? -rabs : rabs) * wi[idx];
22188
1344509a480c style fixes
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
634 # endif
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
635 if (rabs < static_cast<int64_t> (ki[idx]))
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
636 return x; /* 99.3% of the time we return here 1st try */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
637 else if (idx == 0)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
638 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
639 /* As stated in Marsaglia and Tsang
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
640 *
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
641 * For the normal tail, the method of Marsaglia[5] provides:
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
642 * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x,
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
643 * then return r+x. Except that r+x is always in the positive
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
644 * tail!!!! Any thing random might be used to determine the
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
645 * sign, but as we already have r we might as well use it
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
646 *
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
647 * [PAK] but not the bottom 8 bits, since they are all 0 here!
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
648 */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
649 double xx, yy;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
650 do
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
651 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
652 xx = - ZIGGURAT_NOR_INV_R * std::log (RANDU);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
653 yy = - std::log (RANDU);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
654 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
655 while ( yy+yy <= xx*xx);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
656 return ((rabs & 0x100) ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
657 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
658 else if ((fi[idx-1] - fi[idx]) * RANDU + fi[idx] < exp (-0.5*x*x))
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
659 return x;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
660 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
661 }
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
662
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
663 template <>
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
664 OCTAVE_API double
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
665 rand_exponential<double> ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
666 {
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
667 if (initt)
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
668 create_ziggurat_tables ();
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
669
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
670 while (1)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
671 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
672 ZIGINT ri = ERANDI;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
673 const int idx = static_cast<int> (ri & 0xFF);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
674 const double x = ri * we[idx];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
675 if (ri < ke[idx])
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
676 return x; /* 98.9% of the time we return here 1st try */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
677 else if (idx == 0)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
678 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
679 /* As stated in Marsaglia and Tsang
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
680 *
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
681 * For the exponential tail, the method of Marsaglia[5] provides:
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
682 * x = r - ln(U);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
683 */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
684 return ZIGGURAT_EXP_R - std::log (RANDU);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
685 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
686 else if ((fe[idx-1] - fe[idx]) * RANDU + fe[idx] < exp (-x))
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
687 return x;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
688 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
689 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
690
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
691 template <> OCTAVE_API void rand_uniform<double> (octave_idx_type n, double *p)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
692 {
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
693 std::generate_n (p, n, []() { return rand_uniform<double> (); });
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
694 }
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
695
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
696 template <> OCTAVE_API void rand_normal (octave_idx_type n, double *p)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
697 {
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
698 std::generate_n (p, n, []() { return rand_normal<double> (); });
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
699 }
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
700
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
701 template <> OCTAVE_API void rand_exponential (octave_idx_type n, double *p)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
702 {
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
703 std::generate_n (p, n, []() { return rand_exponential<double> (); });
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
704 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
705
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
706 #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
707 #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
708 #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
709 #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
710 #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
711 #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
712
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
713 #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
714 #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
715 #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
716 #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
717 #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
718 #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
719
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
720 static ZIGINT fki[ZIGGURAT_TABLE_SIZE];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
721 static float fwi[ZIGGURAT_TABLE_SIZE], ffi[ZIGGURAT_TABLE_SIZE];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
722 static ZIGINT fke[ZIGGURAT_TABLE_SIZE];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
723 static float fwe[ZIGGURAT_TABLE_SIZE], ffe[ZIGGURAT_TABLE_SIZE];
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
724
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
725 static void
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
726 create_ziggurat_float_tables ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
727 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
728 int i;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
729 float x, x1;
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
730
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
731 /* Ziggurat tables for the normal distribution */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
732 x1 = ZIGGURAT_NOR_R;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
733 fwi[255] = x1 / NMANTISSA;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
734 ffi[255] = exp (-0.5 * x1 * 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
735
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
736 /* Index zero is special for tail strip, where Marsaglia and Tsang
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
737 * defines this as
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
738 * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1,
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
739 * where v is the area of each strip of the ziggurat.
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
740 */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
741 fki[0] = static_cast<ZIGINT> (x1 * ffi[255] / NOR_SECTION_AREA * NMANTISSA);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
742 fwi[0] = NOR_SECTION_AREA / ffi[255] / NMANTISSA;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
743 ffi[0] = 1.;
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
744
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
745 for (i = 254; i > 0; i--)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
746 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
747 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
748 * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y))
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
749 */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
750 x = std::sqrt (-2. * std::log (NOR_SECTION_AREA / x1 + ffi[i+1]));
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
751 fki[i+1] = static_cast<ZIGINT> (x / x1 * NMANTISSA);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
752 fwi[i] = x / NMANTISSA;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
753 ffi[i] = exp (-0.5 * x * x);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
754 x1 = x;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
755 }
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
756
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
757 fki[1] = 0;
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
758
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
759 /* Zigurrat tables for the exponential distribution */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
760 x1 = ZIGGURAT_EXP_R;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
761 fwe[255] = x1 / EMANTISSA;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
762 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
763
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
764 /* Index zero is special for tail strip, where Marsaglia and Tsang
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
765 * defines this as
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
766 * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1,
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
767 * where v is the area of each strip of the ziggurat.
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
768 */
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
769 fke[0] = static_cast<ZIGINT> (x1 * ffe[255] / EXP_SECTION_AREA * EMANTISSA);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
770 fwe[0] = EXP_SECTION_AREA / ffe[255] / EMANTISSA;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
771 ffe[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
772
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
773 for (i = 254; i > 0; i--)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
774 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
775 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
776 * need inverse operator of y = exp(-x) -> x = -ln(y)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
777 */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
778 x = - std::log (EXP_SECTION_AREA / x1 + ffe[i+1]);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
779 fke[i+1] = static_cast<ZIGINT> (x / x1 * EMANTISSA);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
780 fwe[i] = x / EMANTISSA;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
781 ffe[i] = exp (-x);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
782 x1 = x;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
783 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
784 fke[1] = 0;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
785
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
786 inittf = 0;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
787 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
788
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
789 /*
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
790 * Here is the guts of the algorithm. As Marsaglia and Tsang state the
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
791 * algorithm in their paper
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
792 *
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
793 * 1) Calculate a random signed integer j and let i be the index
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
794 * provided by the rightmost 8-bits of j
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
795 * 2) Set x = j * w_i. If j < k_i return x
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
796 * 3) If i = 0, then return x from the tail
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
797 * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
798 * 5) goto step 1
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
799 *
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
800 * Where f is the functional form of the distribution, which for a normal
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
801 * distribution is exp(-0.5*x*x)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
802 */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
803
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
804 template <>
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
805 OCTAVE_API float
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
806 rand_normal<float> ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
807 {
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
808 if (inittf)
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
809 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
810
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
811 while (1)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
812 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
813 /* 32-bit mantissa */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
814 const uint32_t r = randi32 ();
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
815 const uint32_t rabs = r & LMASK;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
816 const int idx = static_cast<int> (r & 0xFF);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
817 const float x = static_cast<int32_t> (r) * fwi[idx];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
818 if (rabs < fki[idx])
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
819 return x; /* 99.3% of the time we return here 1st try */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
820 else if (idx == 0)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
821 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
822 /* As stated in Marsaglia and Tsang
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
823 *
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
824 * For the normal tail, the method of Marsaglia[5] provides:
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
825 * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x,
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
826 * then return r+x. Except that r+x is always in the positive
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
827 * tail!!!! Any thing random might be used to determine the
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
828 * sign, but as we already have r we might as well use it
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
829 *
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
830 * [PAK] but not the bottom 8 bits, since they are all 0 here!
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
831 */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
832 float xx, yy;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
833 do
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
834 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
835 xx = - ZIGGURAT_NOR_INV_R * std::log (RANDU);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
836 yy = - std::log (RANDU);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
837 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
838 while ( yy+yy <= xx*xx);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
839 return ((rabs & 0x100) ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
840 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
841 else if ((ffi[idx-1] - ffi[idx]) * RANDU + ffi[idx] < exp (-0.5*x*x))
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
842 return x;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
843 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
844 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
845
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
846 template <>
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
847 OCTAVE_API float
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
848 rand_exponential<float> ()
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
849 {
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
850 if (inittf)
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
851 create_ziggurat_float_tables ();
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
852
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
853 while (1)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
854 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
855 ZIGINT ri = ERANDI;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
856 const int idx = static_cast<int> (ri & 0xFF);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
857 const float x = ri * fwe[idx];
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
858 if (ri < fke[idx])
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
859 return x; /* 98.9% of the time we return here 1st try */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
860 else if (idx == 0)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
861 {
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
862 /* As stated in Marsaglia and Tsang
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
863 *
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
864 * For the exponential tail, the method of Marsaglia[5] provides:
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
865 * x = r - ln(U);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
866 */
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
867 return ZIGGURAT_EXP_R - std::log (RANDU);
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
868 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
869 else if ((ffe[idx-1] - ffe[idx]) * RANDU + ffe[idx] < exp (-x))
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
870 return x;
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
871 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
872 }
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
873
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
874 template <> OCTAVE_API void rand_uniform (octave_idx_type n, float *p)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
875 {
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
876 std::generate_n (p, n, []() { return rand_uniform<float> (); });
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
877 }
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
878
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
879 template <> OCTAVE_API void rand_normal (octave_idx_type n, float *p)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
880 {
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
881 std::generate_n (p, n, []() { return rand_normal<float> (); });
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
882 }
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
883
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
884 template <> OCTAVE_API void rand_exponential (octave_idx_type n, float *p)
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
885 {
32590
3c2c585965cc maint: C++ style check for liboctave/ before 9.1 release.
Rik <rik@octave.org>
parents: 32426
diff changeset
886 std::generate_n (p, n, []() { return rand_exponential<float> (); });
31607
aac27ad79be6 maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents: 31605
diff changeset
887 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
888
31605
e88a07dec498 maint: Use macros to begin/end C++ namespaces.
Rik <rik@octave.org>
parents: 31142
diff changeset
889 OCTAVE_END_NAMESPACE(octave)