Mercurial > octave

annotate liboctave/numeric/randmtzig.cc @ 27857:81ada4b7d85b

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