Mercurial > octave

annotate liboctave/numeric/randmtzig.cc @ 25433:49e0447413ad

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use templates and move rand functions inside octave namespace * oct-rand.h, oct-rand.cc (rand): Move class inside octave namespace and rename from octave_rand. Use templates to eliminate duplicate code. Change all uses. * randgamma.h, randgamma.cc, randmtzig.h, randmtzig.cc, randpoisson.h, randpoisson.cc: Move functions inside Octave namespace. Use templates to eliminate duplicate code. Change all uses.
author John W. Eaton <jwe@octave.org>
date Sat, 17 Mar 2018 07:02:30 -0500
parents 6652d3823428
children 00f796120a6d
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
25054
6652d3823428 maint: Update copyright dates in all source files.
John W. Eaton <jwe@octave.org>
parents: 24857
diff changeset
3 Copyright (C) 2006-2018 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
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14846
diff changeset
141 static float randu32 (void) returns 32-bit uniform in (0,1)
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14846
diff changeset
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
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
276 octave::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 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
380 static float randu32 (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
381 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
382 return (static_cast<float> (randi32 ()) + 0.5) * (1.0/4294967296.0);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
383 /* divided by 2^32 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
384 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
385
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
386 /* 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
387 static double randu53 (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
388 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
389 const uint32_t a = randi32 () >> 5;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
390 const uint32_t b = randi32 () >> 6;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
391 return (a*67108864.0+b+0.4) * (1.0/9007199254740992.0);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
392 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
393
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
394 /* 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
395 template <>
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
396 double
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
397 rand_uniform<double> (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
398 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
399 return randu53 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
400 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
401
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
402 /* 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
403 template <>
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
404 float
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
405 rand_uniform<float> (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
406 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
407 return randu32 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
408 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
409
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
410 /* ===== Ziggurat normal and exponential generators ===== */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
411
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
412 #define ZIGGURAT_TABLE_SIZE 256
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
413
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
414 #define ZIGGURAT_NOR_R 3.6541528853610088
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
415 #define ZIGGURAT_NOR_INV_R 0.27366123732975828
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
416 #define NOR_SECTION_AREA 0.00492867323399
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
417
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
418 #define ZIGGURAT_EXP_R 7.69711747013104972
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
419 #define ZIGGURAT_EXP_INV_R 0.129918765548341586
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
420 #define EXP_SECTION_AREA 0.0039496598225815571993
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
421
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
422 #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
423 #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
424 #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
425 #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
426 #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
427 #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
428
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
429 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
430 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
431 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
432 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
433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
434 /*
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
435 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
436 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
437 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
438 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
439 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
440 Ziggurat. This has several advantages.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
441
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
442 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
443 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
444 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
445 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
446 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
447 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
448 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
449 4) Compile flag for full 53-bit random mantissa.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
450
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
451 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
452 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
453 ziggurat, even if the code itself isn't...
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
454
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
455 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
456 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
457 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
458 -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
459 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
460 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
461 full 32 bits. We use EMANTISSA for this term.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
462
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
463 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
464 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
465 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
466 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
467 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
468 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
469 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
470 */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
471
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
472 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
473 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
474 int i;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
475 double x, x1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
476
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
477 /* 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
478 x1 = ZIGGURAT_NOR_R;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
479 wi[255] = x1 / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
480 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
481
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
482 /* 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
483 * defines this as
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
484 * 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
485 * 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
486 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
487 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
488 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
489 fi[0] = 1.;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
490
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
491 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
492 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
493 /* 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
494 * 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
495 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
496 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
497 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
498 wi[i] = x / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
499 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
500 x1 = x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
501 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
502
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
503 ki[1] = 0;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
504
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
505 /* 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
506 x1 = ZIGGURAT_EXP_R;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
507 we[255] = x1 / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
508 fe[255] = exp (-x1);
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 /* 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
511 * defines this as
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
512 * 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
513 * 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
514 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
515 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
516 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
517 fe[0] = 1.;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
518
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
519 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
520 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
521 /* 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
522 * 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
523 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
524 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
525 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
526 we[i] = x / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
527 fe[i] = exp (-x);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
528 x1 = x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
529 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
530 ke[1] = 0;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
531
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
532 initt = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
533 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
534
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
535 /*
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
536 * 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
537 * algorithm in their paper
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
538 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
539 * 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
540 * 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
541 * 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
542 * 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
543 * 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
544 * 5) goto step 1
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 * 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
547 * 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
548 */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
549
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
550
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
551 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
552 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
553 if (initt)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
554 create_ziggurat_tables ();
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
555
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
556 while (1)
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 /* 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
559 * 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
560 * 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
561 * 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
562 * 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
563 * 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
564 * 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
565 */
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21720
diff changeset
566 # 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
567 /* 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
568 double x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
569 int si,idx;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
570 uint32_t lo, hi;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
571 int64_t rabs;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
572 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
573 lo = randi32 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
574 idx = lo & 0xFF;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
575 hi = randi32 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
576 si = hi & UMASK;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
577 p[0] = lo;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
578 p[1] = hi & 0x1FFFFF;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
579 x = ( si ? -rabs : rabs ) * wi[idx];
22188
1344509a480c style fixes
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
580 # else
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
581 /* 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
582 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
583 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
584 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
585 const double x = ( (r & 1) ? -rabs : rabs) * wi[idx];
22188
1344509a480c style fixes
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
586 # endif
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
587 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
588 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
589 else if (idx == 0)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
590 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
591 /* 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
592 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
593 * 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
594 * 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
595 * 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
596 * 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
597 * 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
598 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
599 * [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
600 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
601 double xx, yy;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
602 do
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
603 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
604 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
605 yy = - std::log (RANDU);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
606 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
607 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
608 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
609 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
610 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
611 return x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
612 }
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
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
615 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
616 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
617 if (initt)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
618 create_ziggurat_tables ();
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
619
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
620 while (1)
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 ZIGINT ri = ERANDI;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
623 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
624 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
625 if (ri < ke[idx])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
626 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
627 else if (idx == 0)
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 /* 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
630 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
631 * 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
632 * x = r - ln(U);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
633 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
634 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
635 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
636 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
637 return x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
638 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
639 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
640
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
641 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
642 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
643 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
644 }
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 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
647 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
648 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
649 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
650
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
651 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
652 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
653 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
654 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
655
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
656 #undef ZIGINT
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
657 #undef EMANTISSA
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
658 #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
659 #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
660 #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
661 #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
662
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
663 #define 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
664 #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
665 #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
666 #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
667 #define NRANDI randi32() /* 31 bits for mantissa + 1 bit sign */
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
668 #define RANDU randu32()
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
669
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
670 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
671 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
672 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
673 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
674
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
675 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
676 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
677 int i;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
678 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
679
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
680 /* 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
681 x1 = ZIGGURAT_NOR_R;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
682 fwi[255] = x1 / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
683 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
684
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
685 /* 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
686 * defines this as
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
687 * 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
688 * 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
689 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
690 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
691 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
692 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
693
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
694 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
695 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
696 /* 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
697 * 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
698 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
699 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
700 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
701 fwi[i] = x / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
702 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
703 x1 = x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
704 }
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 fki[1] = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
707
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
708 /* 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
709 x1 = ZIGGURAT_EXP_R;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
710 fwe[255] = x1 / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
711 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
712
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
713 /* 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
714 * defines this as
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
715 * 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
716 * 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
717 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
718 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
719 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
720 ffe[0] = 1.;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
721
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
722 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
723 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
724 /* 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
725 * 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
726 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
727 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
728 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
729 fwe[i] = x / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
730 ffe[i] = exp (-x);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
731 x1 = x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
732 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
733 fke[1] = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
734
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
735 inittf = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
736 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
737
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
738 /*
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
739 * 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
740 * algorithm in their paper
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 * 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
743 * 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
744 * 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
745 * 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
746 * 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
747 * 5) goto step 1
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 * 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
750 * 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
751 */
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
752
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
753 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
754 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
755 if (inittf)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
756 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
757
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
758 while (1)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
759 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
760 /* 32-bit mantissa */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
761 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
762 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
763 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
764 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
765 if (rabs < fki[idx])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
766 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
767 else if (idx == 0)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
768 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
769 /* 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
770 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
771 * 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
772 * 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
773 * 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
774 * 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
775 * 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
776 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
777 * [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
778 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
779 float xx, yy;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
780 do
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
781 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
782 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
783 yy = - std::log (RANDU);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
784 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
785 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
786 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
787 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
788 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
789 return x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
790 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
791 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
792
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
793 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
794 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
795 if (inittf)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
796 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
797
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
798 while (1)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
799 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
800 ZIGINT ri = ERANDI;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
801 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
802 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
803 if (ri < fke[idx])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
804 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
805 else if (idx == 0)
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 /* 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
808 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
809 * 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
810 * x = r - ln(U);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
811 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
812 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
813 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
814 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
815 return x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
816 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
817 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
818
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
819 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
820 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
821 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
822 }
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 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
825 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
826 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
827 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
828
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
829 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
830 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
831 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
832 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
833 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
834