Mercurial > octave

annotate liboctave/numeric/randmtzig.cc @ 27933:863ae57eee69

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use Octave coding conventions in liboctave/ * Array-C.cc, Array-d.cc, Array-f.cc, MSparse.cc, Sparse.cc, dim-vector.h, xerbla.cc, aepbalance.cc, eigs-base.cc, gepbalance.cc, oct-fftw.cc, randmtzig.cc, mx-inlines.cc, lo-sysdep.cc, base-list.h, cmd-edit.h, lo-regexp.cc, oct-atomic.h, oct-binmap.h, oct-inttypes.cc, oct-inttypes.h, quit.h, url-transfer.cc: Use Octave coding conventions in liboctave.
author Rik <rik@octave.org>
date Sat, 11 Jan 2020 12:53:20 -0800
parents bd51beb6205e
children a61b651914d6 0a5b15007766
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
1 ////////////////////////////////////////////////////////////////////////
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
2 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
3 // Copyright (C) 2006-2020 The Octave Project Developers
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
4 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
5 // See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
6 // distribution or <https://octave.org/copyright/>.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
7 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
8 // This file is part of Octave.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
9 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
10 // Octave is free software: you can redistribute it and/or modify it
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
11 // under the terms of the GNU General Public License as published by
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
12 // the Free Software Foundation, either version 3 of the License, or
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
13 // (at your option) any later version.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
14 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
15 // Octave is distributed in the hope that it will be useful, but
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
18 // GNU General Public License for more details.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
19 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
20 // You should have received a copy of the GNU General Public License
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
21 // along with Octave; see the file COPYING. If not, see
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
22 // <https://www.gnu.org/licenses/>.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
23 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
24 ////////////////////////////////////////////////////////////////////////
7019
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6959
diff changeset
25
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
26 /*
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
27 A C-program for MT19937, with initialization improved 2002/2/10.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
28 Coded by Takuji Nishimura and Makoto Matsumoto.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
29 This is a faster version by taking Shawn Cokus's optimization,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
30 Matthe Bellew's simplification, Isaku Wada's real version.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
31 David Bateman added normal and exponential distributions following
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
32 Marsaglia and Tang's Ziggurat algorithm.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
33
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
34 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
35 Copyright (C) 2004, David Bateman
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
36 All rights reserved.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
37
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
38 Redistribution and use in source and binary forms, with or without
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
39 modification, are permitted provided that the following conditions
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
40 are met:
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
41
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
42 1. Redistributions of source code must retain the above copyright
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
43 notice, this list of conditions and the following disclaimer.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
44
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
45 2. Redistributions in binary form must reproduce the above copyright
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
46 notice, this list of conditions and the following disclaimer in the
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
47 documentation and/or other materials provided with the distribution.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
48
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
49 3. The names of its contributors may not be used to endorse or promote
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
50 products derived from this software without specific prior written
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
51 permission.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
52
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
53 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
54 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
55 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
56 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
57 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
58 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
59 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
60 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
61 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
62 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
63 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
64
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
65
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
66 Any feedback is very welcome.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
67 http://www.math.keio.ac.jp/matumoto/emt.html
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
68 email: matumoto@math.keio.ac.jp
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
69
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
70 * 2004-01-19 Paul Kienzle
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
71 * * comment out main
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
72 * add init_by_entropy, get_state, set_state
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
73 * * converted to allow compiling by C++ compiler
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
74 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
75 * 2004-01-25 David Bateman
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
76 * * Add Marsaglia and Tsang Ziggurat code
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
77 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
78 * 2004-07-13 Paul Kienzle
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
79 * * make into an independent library with some docs.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
80 * * introduce new main and test code.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
81 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
82 * 2004-07-28 Paul Kienzle & David Bateman
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
83 * * add -DALLBITS flag for 32 vs. 53 bits of randomness in mantissa
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
84 * * make the naming scheme more uniform
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
85 * * add -DHAVE_X86 for faster support of 53 bit mantissa on x86 arch.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
86 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
87 * 2005-02-23 Paul Kienzle
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
88 * * fix -DHAVE_X86_32 flag and add -DUSE_X86_32=0|1 for explicit control
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
89 *
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
90 * 2006-04-01 David Bateman
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
91 * * convert for use in octave, declaring static functions only used
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
92 * here and adding oct_ to functions visible externally
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
93 * * inverse sense of ALLBITS
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
94 *
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
95 * 2012-05-18 David Bateman
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
96 * * Remove randu64 and ALLBIT option
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
97 * * Add the single precision generators
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
98 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
99
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
100 /*
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
101 === Build instructions ===
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
102
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
103 Compile with -DHAVE_GETTIMEOFDAY if the gettimeofday function is
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
104 available. This is not necessary if your architecture has
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
105 /dev/urandom defined.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
106
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
107 Uses implicit -Di386 or explicit -DHAVE_X86_32 to determine if CPU=x86.
21308
c53bfd6d8e08 maint: Use American spelling for "behavior".
Rik <rik@octave.org>
parents: 21301
diff changeset
108 You can force X86 behavior with -DUSE_X86_32=1, or suppress it with
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
109 -DUSE_X86_32=0. You should also consider -march=i686 or similar for
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
110 extra performance. Check whether -DUSE_X86_32=0 is faster on 64-bit
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
111 x86 architectures.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
112
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
113 If you want to replace the Mersenne Twister with another
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
114 generator then redefine randi32 appropriately.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
115
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
116 === Usage instructions ===
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
117 Before using any of the generators, initialize the state with one of
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
118 the init_mersenne_twister functions.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
119
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
120 All generators share the same state vector.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
121
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
122 === Mersenne Twister ===
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
123 random initial state:
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
124 void init_mersenne_twister (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
125
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
126 // 32-bit initial state:
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
127 void init_mersenne_twister (uint32_t s)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
128
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
129 // m*32-bit initial state:
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
130 void init_mersenne_twister (uint32_t k[],int m)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
131
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
132 // saves state in array:
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
133 void get_mersenne_twister_state (uint32_t save[MT_N+1])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
134
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
135 // restores state from array
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
136 void set_mersenne_twister_state (uint32_t save[MT_N+1])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
137
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14846
diff changeset
138 static uint32_t randmt (void) returns 32-bit unsigned int
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
139
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
140 === inline generators ===
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14846
diff changeset
141 static uint32_t randi32 (void) returns 32-bit unsigned int
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14846
diff changeset
142 static uint64_t randi53 (void) returns 53-bit unsigned int
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14846
diff changeset
143 static uint64_t randi54 (void) returns 54-bit unsigned int
27855
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
144 static float randu24 (void) returns 24-bit uniform in (0,1)
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14846
diff changeset
145 static double randu53 (void) returns 53-bit uniform in (0,1)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
146
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
147 double rand_uniform (void) returns M-bit uniform in (0,1)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
148 double rand_normal (void) returns M-bit standard normal
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
149 double rand_exponential (void) returns N-bit standard exponential
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
150
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
151 === Array generators ===
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
152 void rand_uniform (octave_idx_type, double [])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
153 void rand_normal (octave_idx_type, double [])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
154 void rand_exponential (octave_idx_type, double [])
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
155 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
156
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
157 #if defined (HAVE_CONFIG_H)
21301
40de9f8f23a6 Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents: 21213
diff changeset
158 # include "config.h"
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
159 #endif
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
160
23662
bd77ab816e43 eliminate obsolete file lo-math.h
John W. Eaton <jwe@octave.org>
parents: 23649
diff changeset
161 #include <cmath>
21720
b28c26ae737c use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents: 21661
diff changeset
162 #include <cstdio>
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
163
24856
8bb0251fcfde Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents: 24534
diff changeset
164 #include <algorithm>
8bb0251fcfde Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents: 24534
diff changeset
165
21928
315f4ba604c8 * randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents: 21746
diff changeset
166 #include "oct-time.h"
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
167 #include "randmtzig.h"
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
168
21066
258c787cd9ce maint: Use "FIXME:" consistently in code base.
Rik <rik@octave.org>
parents: 20955
diff changeset
169 /* 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
170 #if ! defined (USE_X86_32)
21202
f7121e111991 maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents: 21066
diff changeset
171 # 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
172 # define USE_X86_32 1
f7121e111991 maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents: 21066
diff changeset
173 # else
f7121e111991 maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents: 21066
diff changeset
174 # define USE_X86_32 0
f7121e111991 maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents: 21066
diff changeset
175 # endif
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
176 #endif
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
177
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
178 namespace octave
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
179 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
180 /* ===== Mersenne Twister 32-bit generator ===== */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
181
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
182 #define MT_M 397
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
183 #define MATRIX_A 0x9908b0dfUL /* constant vector a */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
184 #define UMASK 0x80000000UL /* most significant w-r bits */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
185 #define LMASK 0x7fffffffUL /* least significant r bits */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
186 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
187 #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
188
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
189 static uint32_t *next;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
190 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
191 static int left = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
192 static int initf = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
193 static int initt = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
194 static int inittf = 1;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
195
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
196 /* 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
197 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
198 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
199 int j;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
200 state[0] = s & 0xffffffffUL;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
201 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
202 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
203 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
204 /* 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
205 /* 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
206 /* 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
207 /* 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
208 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
209 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
210 left = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
211 initf = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
212 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
213
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
214 /* 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
215 /* 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
216 /* 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
217 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
218 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
219 int i, j, k;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
220 init_mersenne_twister (19650218UL);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
221 i = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
222 j = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
223 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
224 for (; k; k--)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
225 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
226 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL))
27933
863ae57eee69 maint: Use Octave coding conventions in liboctave/
Rik <rik@octave.org>
parents: 27923
diff changeset
227 + init_key[j] + j; /* non linear */
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
228 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
229 i++;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
230 j++;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
231 if (i >= MT_N)
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 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
234 i = 1;
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 if (j >= key_length)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
237 j = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
238 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
239 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
240 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
241 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL))
27933
863ae57eee69 maint: Use Octave coding conventions in liboctave/
Rik <rik@octave.org>
parents: 27923
diff changeset
242 - i; /* non linear */
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
243 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
244 i++;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
245 if (i >= MT_N)
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 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
248 i = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
249 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
250 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
251
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
252 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
253 left = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
254 initf = 1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
255 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
256
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
257 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
258 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
259 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
260 int n = 0;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
261
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
262 /* 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
263 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
264 if (urandom)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
265 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
266 while (n < MT_N)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
267 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
268 unsigned char word[4];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
269 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
270 break;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
271 entropy[n++] = word[0] + (word[1]<<8) + (word[2]<<16)
27933
863ae57eee69 maint: Use Octave coding conventions in liboctave/
Rik <rik@octave.org>
parents: 27923
diff changeset
272 + (static_cast<uint32_t> (word[3])<<24);
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
273 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
274 std::fclose (urandom);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
275 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
276
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
277 /* 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
278
27102
84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
279 sys::time now;
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++] = 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
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++] = clock (); /* CPU time used (usec) */
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 if (n < MT_N)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
288 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
289
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
290 /* 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
291 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
292 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
293
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
294 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
295 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
296 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
297 left = save[MT_N];
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
298 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
299 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
300
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
301 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
302 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
303 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
304 save[MT_N] = left;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
305 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
306
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
307 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
308 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
309 uint32_t *p = state;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
310 int j;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
311
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
312 /* 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
313 /* 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
314 /* 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
315 /* 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
316 if (initf == 0)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
317 init_mersenne_twister ();
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
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 left = MT_N;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
320 next = state;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
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_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
323 *p = p[MT_M] ^ 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 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
326 *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
327
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
328 *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
329 }
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 /* 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
332 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
333 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
334 uint32_t y;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
335
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
336 if (--left == 0)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
337 next_state ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
338 y = *next++;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
339
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
340 /* Tempering */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
341 y ^= (y >> 11);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
342 y ^= (y << 7) & 0x9d2c5680UL;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
343 y ^= (y << 15) & 0xefc60000UL;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
344 return (y ^ (y >> 18));
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
345 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
346
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
347 /* ===== Uniform generators ===== */
5742
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 /* Select which 32 bit generator to use */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
350 #define randi32 randmt
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
351
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
352 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
353 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
354 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
355 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
356 #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
357 uint64_t u;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
358 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
359 p[0] = lo;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
360 p[1] = hi;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
361 return u;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
362 #else
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
363 return ((static_cast<uint64_t> (hi) << 32) | lo);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
364 #endif
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
365 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
366
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
367 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
368 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
369 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
370 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
371 #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
372 uint64_t u;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
373 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
374 p[0] = lo;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
375 p[1] = hi;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
376 return u;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
377 #else
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
378 return ((static_cast<uint64_t> (hi) << 32) | lo);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
379 #endif
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
380 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
381
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
382 /* generates a random number on (0,1)-real-interval */
27855
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
383 static float randu24 (void)
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
384 {
27855
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
385 uint32_t i;
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
386
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
387 do
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
388 {
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
389 i = randi32 () & static_cast<uint32_t> (0xFFFFFF);
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
390 }
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
391 while (i == 0);
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
392
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
393 return i * (1.0f / 16777216.0f);
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
394 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
395
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
396 /* 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
397 static double randu53 (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
398 {
27858
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
399 int32_t a, b;
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
400
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
401 do
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
402 {
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
403 a = randi32 () >> 5;
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
404 b = randi32 () >> 6;
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
405 }
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
406 while (a == 0 && b == 0);
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
407
797be8d10c22 Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents: 27857
diff changeset
408 return (a*67108864.0 + b) * (1.0/9007199254740992.0);
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
409 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
410
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
411 /* 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
412 template <>
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
413 double
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
414 rand_uniform<double> (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
415 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
416 return randu53 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
417 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
418
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
419 /* 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
420 template <>
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
421 float
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
422 rand_uniform<float> (void)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
423 {
27855
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
424 return randu24 ();
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
425 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
426
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
427 /* ===== Ziggurat normal and exponential generators ===== */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
428
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
429 #define ZIGGURAT_TABLE_SIZE 256
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
430
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
431 #define ZIGGURAT_NOR_R 3.6541528853610088
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
432 #define ZIGGURAT_NOR_INV_R 0.27366123732975828
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
433 #define NOR_SECTION_AREA 0.00492867323399
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
434
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
435 #define ZIGGURAT_EXP_R 7.69711747013104972
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
436 #define ZIGGURAT_EXP_INV_R 0.129918765548341586
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
437 #define EXP_SECTION_AREA 0.0039496598225815571993
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
438
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
439 #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
440 #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
441 #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
442 #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
443 #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
444 #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
445
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
446 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
447 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
448 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
449 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
450
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
451 /*
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
452 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
453 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
454 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
455 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
456 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
457 Ziggurat. This has several advantages.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
458
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
459 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
460 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
461 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
462 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
463 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
464 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
465 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
466 4) Compile flag for full 53-bit random mantissa.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
467
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
468 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
469 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
470 ziggurat, even if the code itself isn't...
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 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
473 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
474 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
475 -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
476 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
477 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
478 full 32 bits. We use EMANTISSA for this term.
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
479
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
480 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
481 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
482 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
483 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
484 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
485 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
486 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
487 */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
488
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
489 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
490 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
491 int i;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
492 double x, x1;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
493
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
494 /* 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
495 x1 = ZIGGURAT_NOR_R;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
496 wi[255] = x1 / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
497 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
498
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
499 /* 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
500 * defines this as
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
501 * 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
502 * 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
503 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
504 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
505 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
506 fi[0] = 1.;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
507
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
508 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
509 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
510 /* 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
511 * 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
512 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
513 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
514 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
515 wi[i] = x / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
516 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
517 x1 = x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
518 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
519
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
520 ki[1] = 0;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
521
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
522 /* 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
523 x1 = ZIGGURAT_EXP_R;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
524 we[255] = x1 / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
525 fe[255] = exp (-x1);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
526
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
527 /* 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
528 * defines this as
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
529 * 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
530 * 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
531 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
532 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
533 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
534 fe[0] = 1.;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
535
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
536 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
537 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
538 /* 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
539 * 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
540 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
541 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
542 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
543 we[i] = x / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
544 fe[i] = exp (-x);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
545 x1 = x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
546 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
547 ke[1] = 0;
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
548
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
549 initt = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
550 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
551
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
552 /*
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
553 * 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
554 * algorithm in their paper
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
555 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
556 * 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
557 * 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
558 * 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
559 * 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
560 * 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
561 * 5) goto step 1
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
562 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
563 * 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
564 * 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
565 */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
566
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
567
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
568 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
569 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
570 if (initt)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
571 create_ziggurat_tables ();
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
572
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
573 while (1)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
574 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
575 /* 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
576 * 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
577 * 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
578 * 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
579 * 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
580 * 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
581 * 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
582 */
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21720
diff changeset
583 # 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
584 /* 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
585 double x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
586 int si,idx;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
587 uint32_t lo, hi;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
588 int64_t rabs;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
589 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
590 lo = randi32 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
591 idx = lo & 0xFF;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
592 hi = randi32 ();
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
593 si = hi & UMASK;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
594 p[0] = lo;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
595 p[1] = hi & 0x1FFFFF;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
596 x = ( si ? -rabs : rabs ) * wi[idx];
22188
1344509a480c style fixes
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
597 # else
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
598 /* 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
599 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
600 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
601 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
602 const double x = ( (r & 1) ? -rabs : rabs) * wi[idx];
22188
1344509a480c style fixes
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
603 # endif
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
604 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
605 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
606 else if (idx == 0)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
607 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
608 /* 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
609 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
610 * 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
611 * 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
612 * 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
613 * 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
614 * 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
615 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
616 * [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
617 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
618 double xx, yy;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
619 do
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
620 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
621 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
622 yy = - std::log (RANDU);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
623 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
624 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
625 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
626 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
627 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
628 return x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
629 }
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
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
632 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
633 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
634 if (initt)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
635 create_ziggurat_tables ();
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
636
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
637 while (1)
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 ZIGINT ri = ERANDI;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
640 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
641 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
642 if (ri < ke[idx])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
643 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
644 else if (idx == 0)
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 /* 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
647 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
648 * 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
649 * x = r - ln(U);
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 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
652 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
653 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
654 return x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
655 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
656 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
657
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
658 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
659 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
660 std::generate_n (p, n, [](void) { return rand_uniform<double> (); });
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
661 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
662
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
663 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
664 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
665 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
666 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
667
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
668 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
669 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
670 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
671 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
672
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
673 #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
674 #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
675 #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
676 #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
677 #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
678 #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
679
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
680 #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
681 #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
682 #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
683 #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
684 #define NRANDI randi32() /* 31 bits for mantissa + 1 bit sign */
27855
9405e2be91d0 Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents: 27102
diff changeset
685 #define RANDU randu24()
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
686
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
687 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
688 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
689 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
690 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
691
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
692 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
693 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
694 int i;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
695 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
696
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
697 /* 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
698 x1 = ZIGGURAT_NOR_R;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
699 fwi[255] = x1 / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
700 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
701
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
702 /* 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
703 * defines this as
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
704 * 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
705 * 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
706 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
707 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
708 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
709 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
710
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
711 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
712 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
713 /* 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
714 * 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
715 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
716 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
717 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
718 fwi[i] = x / NMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
719 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
720 x1 = x;
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
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
723 fki[1] = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
724
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
725 /* 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
726 x1 = ZIGGURAT_EXP_R;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
727 fwe[255] = x1 / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
728 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
729
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
730 /* 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
731 * defines this as
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
732 * 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
733 * 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
734 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
735 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
736 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
737 ffe[0] = 1.;
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 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
740 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
741 /* 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
742 * 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
743 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
744 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
745 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
746 fwe[i] = x / EMANTISSA;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
747 ffe[i] = exp (-x);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
748 x1 = x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
749 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
750 fke[1] = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
751
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
752 inittf = 0;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
753 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
754
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
755 /*
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
756 * 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
757 * algorithm in their paper
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
758 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
759 * 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
760 * 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
761 * 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
762 * 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
763 * 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
764 * 5) goto step 1
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
765 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
766 * 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
767 * 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
768 */
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
769
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
770 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
771 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
772 if (inittf)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
773 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
774
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
775 while (1)
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 /* 32-bit mantissa */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
778 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
779 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
780 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
781 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
782 if (rabs < fki[idx])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
783 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
784 else if (idx == 0)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
785 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
786 /* 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
787 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
788 * 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
789 * 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
790 * 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
791 * 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
792 * 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
793 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
794 * [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
795 */
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
796 float xx, yy;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
797 do
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
798 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
799 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
800 yy = - std::log (RANDU);
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
801 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
802 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
803 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
804 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
805 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
806 return x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
807 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
808 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
809
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
810 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
811 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
812 if (inittf)
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
813 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
814
25433
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
815 while (1)
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 ZIGINT ri = ERANDI;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
818 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
819 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
820 if (ri < fke[idx])
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
821 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
822 else if (idx == 0)
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 /* 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
825 *
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
826 * 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
827 * x = r - ln(U);
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 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
830 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
831 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
832 return x;
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
833 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
834 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
835
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
836 template <> void rand_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
837 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
838 std::generate_n (p, n, [](void) { return rand_uniform<float> (); });
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
839 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
840
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
841 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
842 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
843 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
844 }
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
845
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
846 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
847 {
49e0447413ad use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 25054
diff changeset
848 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
849 }
14655
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
850 }
43db83eff9db Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents: 14138
diff changeset
851