Mercurial > octave
annotate liboctave/numeric/randmtzig.cc @ 31605:e88a07dec498 stable
maint: Use macros to begin/end C++ namespaces.
* oct-conf-post-public.in.h: Define two macros (OCTAVE_BEGIN_NAMESPACE,
OCTAVE_END_NAMESPACE) that can be used to start/end a namespace.
* mk-opts.pl, build-env.h, build-env.in.cc, __betainc__.cc, __contourc__.cc,
__dsearchn__.cc, __eigs__.cc, __expint__.cc, __ftp__.cc, __gammainc__.cc,
__ichol__.cc, __ilu__.cc, __isprimelarge__.cc, __lin_interpn__.cc,
__magick_read__.cc, __pchip_deriv__.cc, __qp__.cc, amd.cc, auto-shlib.cc,
auto-shlib.h, balance.cc, base-text-renderer.cc, base-text-renderer.h,
besselj.cc, bitfcns.cc, bsxfun.cc, c-file-ptr-stream.cc, c-file-ptr-stream.h,
call-stack.cc, call-stack.h, ccolamd.cc, cellfun.cc, chol.cc, colamd.cc,
colloc.cc, conv2.cc, daspk.cc, dasrt.cc, dassl.cc, data.cc, data.h, debug.cc,
defaults.cc, defaults.h, defun-int.h, defun.cc, det.cc, dirfns.cc, display.cc,
display.h, dlmread.cc, dmperm.cc, dot.cc, dynamic-ld.cc, dynamic-ld.h, eig.cc,
ellipj.cc, environment.cc, environment.h, error.cc, error.h, errwarn.h,
event-manager.cc, event-manager.h, event-queue.cc, event-queue.h, fcn-info.cc,
fcn-info.h, fft.cc, fft2.cc, fftn.cc, file-io.cc, filter.cc, find.cc,
ft-text-renderer.cc, ft-text-renderer.h, gcd.cc, getgrent.cc, getpwent.cc,
getrusage.cc, givens.cc, gl-render.cc, gl-render.h, gl2ps-print.cc,
gl2ps-print.h, graphics-toolkit.cc, graphics-toolkit.h, graphics.cc,
graphics.in.h, gsvd.cc, gtk-manager.cc, gtk-manager.h, hash.cc, help.cc,
help.h, hess.cc, hex2num.cc, hook-fcn.cc, hook-fcn.h, input.cc, input.h,
interpreter-private.cc, interpreter-private.h, interpreter.cc, interpreter.h,
inv.cc, jsondecode.cc, jsonencode.cc, kron.cc, latex-text-renderer.cc,
latex-text-renderer.h, load-path.cc, load-path.h, load-save.cc, load-save.h,
lookup.cc, ls-ascii-helper.cc, ls-ascii-helper.h, ls-oct-text.cc, ls-utils.cc,
ls-utils.h, lsode.cc, lu.cc, mappers.cc, matrix_type.cc, max.cc, mex-private.h,
mex.cc, mgorth.cc, nproc.cc, oct-fstrm.cc, oct-fstrm.h, oct-hdf5-types.cc,
oct-hdf5-types.h, oct-hist.cc, oct-hist.h, oct-iostrm.cc, oct-iostrm.h,
oct-opengl.h, oct-prcstrm.cc, oct-prcstrm.h, oct-procbuf.cc, oct-procbuf.h,
oct-process.cc, oct-process.h, oct-stdstrm.h, oct-stream.cc, oct-stream.h,
oct-strstrm.cc, oct-strstrm.h, oct-tex-lexer.in.ll, oct-tex-parser.yy,
ordqz.cc, ordschur.cc, pager.cc, pager.h, pinv.cc, pow2.cc, pr-flt-fmt.cc,
pr-output.cc, procstream.cc, procstream.h, psi.cc, qr.cc, quad.cc, quadcc.cc,
qz.cc, rand.cc, rcond.cc, regexp.cc, schur.cc, settings.cc, settings.h,
sighandlers.cc, sighandlers.h, sparse-xdiv.cc, sparse-xdiv.h, sparse-xpow.cc,
sparse-xpow.h, sparse.cc, spparms.cc, sqrtm.cc, stack-frame.cc, stack-frame.h,
stream-euler.cc, strfind.cc, strfns.cc, sub2ind.cc, svd.cc, sylvester.cc,
symbfact.cc, syminfo.cc, syminfo.h, symrcm.cc, symrec.cc, symrec.h,
symscope.cc, symscope.h, symtab.cc, symtab.h, syscalls.cc, sysdep.cc, sysdep.h,
text-engine.cc, text-engine.h, text-renderer.cc, text-renderer.h, time.cc,
toplev.cc, tril.cc, tsearch.cc, typecast.cc, url-handle-manager.cc,
url-handle-manager.h, urlwrite.cc, utils.cc, utils.h, variables.cc,
variables.h, xdiv.cc, xdiv.h, xnorm.cc, xnorm.h, xpow.cc, xpow.h,
__delaunayn__.cc, __fltk_uigetfile__.cc, __glpk__.cc, __init_fltk__.cc,
__init_gnuplot__.cc, __ode15__.cc, __voronoi__.cc, audiodevinfo.cc,
audioread.cc, convhulln.cc, fftw.cc, gzip.cc, mk-build-env-features.sh,
mk-builtins.pl, cdef-class.cc, cdef-class.h, cdef-fwd.h, cdef-manager.cc,
cdef-manager.h, cdef-method.cc, cdef-method.h, cdef-object.cc, cdef-object.h,
cdef-package.cc, cdef-package.h, cdef-property.cc, cdef-property.h,
cdef-utils.cc, cdef-utils.h, ov-base.cc, ov-base.h, ov-bool-mat.cc,
ov-builtin.h, ov-cell.cc, ov-class.cc, ov-class.h, ov-classdef.cc,
ov-classdef.h, ov-complex.cc, ov-fcn-handle.cc, ov-fcn-handle.h, ov-fcn.h,
ov-java.cc, ov-java.h, ov-mex-fcn.h, ov-null-mat.cc, ov-oncleanup.cc,
ov-struct.cc, ov-typeinfo.cc, ov-typeinfo.h, ov-usr-fcn.cc, ov-usr-fcn.h,
ov.cc, ov.h, octave.cc, octave.h, mk-ops.sh, op-b-b.cc, op-b-bm.cc,
op-b-sbm.cc, op-bm-b.cc, op-bm-bm.cc, op-bm-sbm.cc, op-cdm-cdm.cc, op-cell.cc,
op-chm.cc, op-class.cc, op-cm-cm.cc, op-cm-cs.cc, op-cm-m.cc, op-cm-s.cc,
op-cm-scm.cc, op-cm-sm.cc, op-cs-cm.cc, op-cs-cs.cc, op-cs-m.cc, op-cs-s.cc,
op-cs-scm.cc, op-cs-sm.cc, op-dm-dm.cc, op-dm-scm.cc, op-dm-sm.cc,
op-dm-template.cc, op-dms-template.cc, op-fcdm-fcdm.cc, op-fcm-fcm.cc,
op-fcm-fcs.cc, op-fcm-fm.cc, op-fcm-fs.cc, op-fcn.cc, op-fcs-fcm.cc,
op-fcs-fcs.cc, op-fcs-fm.cc, op-fcs-fs.cc, op-fdm-fdm.cc, op-fm-fcm.cc,
op-fm-fcs.cc, op-fm-fm.cc, op-fm-fs.cc, op-fs-fcm.cc, op-fs-fcs.cc,
op-fs-fm.cc, op-fs-fs.cc, op-i16-i16.cc, op-i32-i32.cc, op-i64-i64.cc,
op-i8-i8.cc, op-int-concat.cc, op-m-cm.cc, op-m-cs.cc, op-m-m.cc, op-m-s.cc,
op-m-scm.cc, op-m-sm.cc, op-mi.cc, op-pm-pm.cc, op-pm-scm.cc, op-pm-sm.cc,
op-pm-template.cc, op-range.cc, op-s-cm.cc, op-s-cs.cc, op-s-m.cc, op-s-s.cc,
op-s-scm.cc, op-s-sm.cc, op-sbm-b.cc, op-sbm-bm.cc, op-sbm-sbm.cc,
op-scm-cm.cc, op-scm-cs.cc, op-scm-m.cc, op-scm-s.cc, op-scm-scm.cc,
op-scm-sm.cc, op-sm-cm.cc, op-sm-cs.cc, op-sm-m.cc, op-sm-s.cc, op-sm-scm.cc,
op-sm-sm.cc, op-str-m.cc, op-str-s.cc, op-str-str.cc, op-struct.cc,
op-ui16-ui16.cc, op-ui32-ui32.cc, op-ui64-ui64.cc, op-ui8-ui8.cc, ops.h,
anon-fcn-validator.cc, anon-fcn-validator.h, bp-table.cc, bp-table.h,
comment-list.cc, comment-list.h, filepos.h, lex.h, lex.ll, oct-lvalue.cc,
oct-lvalue.h, oct-parse.yy, parse.h, profiler.cc, profiler.h,
pt-anon-scopes.cc, pt-anon-scopes.h, pt-arg-list.cc, pt-arg-list.h,
pt-args-block.cc, pt-args-block.h, pt-array-list.cc, pt-array-list.h,
pt-assign.cc, pt-assign.h, pt-binop.cc, pt-binop.h, pt-bp.cc, pt-bp.h,
pt-cbinop.cc, pt-cbinop.h, pt-cell.cc, pt-cell.h, pt-check.cc, pt-check.h,
pt-classdef.cc, pt-classdef.h, pt-cmd.h, pt-colon.cc, pt-colon.h, pt-const.cc,
pt-const.h, pt-decl.cc, pt-decl.h, pt-eval.cc, pt-eval.h, pt-except.cc,
pt-except.h, pt-exp.cc, pt-exp.h, pt-fcn-handle.cc, pt-fcn-handle.h, pt-id.cc,
pt-id.h, pt-idx.cc, pt-idx.h, pt-jump.h, pt-loop.cc, pt-loop.h, pt-mat.cc,
pt-mat.h, pt-misc.cc, pt-misc.h, pt-pr-code.cc, pt-pr-code.h, pt-select.cc,
pt-select.h, pt-spmd.cc, pt-spmd.h, pt-stmt.cc, pt-stmt.h, pt-tm-const.cc,
pt-tm-const.h, pt-unop.cc, pt-unop.h, pt-vm-eval.cc, pt-walk.cc, pt-walk.h,
pt.cc, pt.h, token.cc, token.h, Range.cc, Range.h, idx-vector.cc, idx-vector.h,
range-fwd.h, CollocWt.cc, CollocWt.h, aepbalance.cc, aepbalance.h, chol.cc,
chol.h, gepbalance.cc, gepbalance.h, gsvd.cc, gsvd.h, hess.cc, hess.h,
lo-mappers.cc, lo-mappers.h, lo-specfun.cc, lo-specfun.h, lu.cc, lu.h,
oct-convn.cc, oct-convn.h, oct-fftw.cc, oct-fftw.h, oct-norm.cc, oct-norm.h,
oct-rand.cc, oct-rand.h, oct-spparms.cc, oct-spparms.h, qr.cc, qr.h, qrp.cc,
qrp.h, randgamma.cc, randgamma.h, randmtzig.cc, randmtzig.h, randpoisson.cc,
randpoisson.h, schur.cc, schur.h, sparse-chol.cc, sparse-chol.h, sparse-lu.cc,
sparse-lu.h, sparse-qr.cc, sparse-qr.h, svd.cc, svd.h, child-list.cc,
child-list.h, dir-ops.cc, dir-ops.h, file-ops.cc, file-ops.h, file-stat.cc,
file-stat.h, lo-sysdep.cc, lo-sysdep.h, lo-sysinfo.cc, lo-sysinfo.h,
mach-info.cc, mach-info.h, oct-env.cc, oct-env.h, oct-group.cc, oct-group.h,
oct-password.cc, oct-password.h, oct-syscalls.cc, oct-syscalls.h, oct-time.cc,
oct-time.h, oct-uname.cc, oct-uname.h, action-container.cc, action-container.h,
base-list.h, cmd-edit.cc, cmd-edit.h, cmd-hist.cc, cmd-hist.h, f77-fcn.h,
file-info.cc, file-info.h, lo-array-errwarn.cc, lo-array-errwarn.h, lo-hash.cc,
lo-hash.h, lo-ieee.h, lo-regexp.cc, lo-regexp.h, lo-utils.cc, lo-utils.h,
oct-base64.cc, oct-base64.h, oct-glob.cc, oct-glob.h, oct-inttypes.h,
oct-mutex.cc, oct-mutex.h, oct-refcount.h, oct-shlib.cc, oct-shlib.h,
oct-sparse.cc, oct-sparse.h, oct-string.h, octave-preserve-stream-state.h,
pathsearch.cc, pathsearch.h, quit.cc, quit.h, unwind-prot.cc, unwind-prot.h,
url-transfer.cc, url-transfer.h : Use new macros to begin/end C++ namespaces.
author | Rik <rik@octave.org> |
---|---|
date | Thu, 01 Dec 2022 14:23:45 -0800 |
parents | 8245e773bb5b |
children | aac27ad79be6 |
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 // |
30564
796f54d4ddbf
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
30394
diff
changeset
|
3 // Copyright (C) 2006-2022 The Octave Project Developers |
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
4 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
5 // See the file COPYRIGHT.md in the top-level directory of this |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
6 // distribution or <https://octave.org/copyright/>. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
7 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
8 // This file is part of Octave. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
9 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
10 // Octave is free software: you can redistribute it and/or modify it |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
11 // under the terms of the GNU General Public License as published by |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
12 // the Free Software Foundation, either version 3 of the License, or |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
13 // (at your option) any later version. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
14 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
15 // Octave is distributed in the hope that it will be useful, but |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
16 // WITHOUT ANY WARRANTY; without even the implied warranty of |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
18 // GNU General Public License for more details. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
19 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
20 // You should have received a copy of the GNU General Public License |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
21 // along with Octave; see the file COPYING. If not, see |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
22 // <https://www.gnu.org/licenses/>. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
23 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
24 //////////////////////////////////////////////////////////////////////// |
7019 | 25 |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
26 /* |
5742 | 27 A C-program for MT19937, with initialization improved 2002/2/10. |
28 Coded by Takuji Nishimura and Makoto Matsumoto. | |
29 This is a faster version by taking Shawn Cokus's optimization, | |
30 Matthe Bellew's simplification, Isaku Wada's real version. | |
31 David Bateman added normal and exponential distributions following | |
32 Marsaglia and Tang's Ziggurat algorithm. | |
33 | |
34 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, | |
35 Copyright (C) 2004, David Bateman | |
36 All rights reserved. | |
37 | |
38 Redistribution and use in source and binary forms, with or without | |
39 modification, are permitted provided that the following conditions | |
40 are met: | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
41 |
5742 | 42 1. Redistributions of source code must retain the above copyright |
43 notice, this list of conditions and the following disclaimer. | |
44 | |
45 2. Redistributions in binary form must reproduce the above copyright | |
46 notice, this list of conditions and the following disclaimer in the | |
47 documentation and/or other materials provided with the distribution. | |
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 | 51 permission. |
52 | |
53 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS | |
54 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT | |
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 | 57 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
58 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, | |
59 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR | |
60 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF | |
61 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING | |
62 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS | |
63 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | |
64 | |
65 | |
66 Any feedback is very welcome. | |
67 http://www.math.keio.ac.jp/matumoto/emt.html | |
68 email: matumoto@math.keio.ac.jp | |
69 | |
70 * 2004-01-19 Paul Kienzle | |
71 * * comment out main | |
72 * add init_by_entropy, get_state, set_state | |
73 * * converted to allow compiling by C++ compiler | |
74 * | |
75 * 2004-01-25 David Bateman | |
76 * * Add Marsaglia and Tsang Ziggurat code | |
77 * | |
78 * 2004-07-13 Paul Kienzle | |
79 * * make into an independent library with some docs. | |
80 * * introduce new main and test code. | |
81 * | |
82 * 2004-07-28 Paul Kienzle & David Bateman | |
83 * * add -DALLBITS flag for 32 vs. 53 bits of randomness in mantissa | |
84 * * make the naming scheme more uniform | |
85 * * add -DHAVE_X86 for faster support of 53 bit mantissa on x86 arch. | |
86 * | |
87 * 2005-02-23 Paul Kienzle | |
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 | 98 */ |
99 | |
100 /* | |
101 === Build instructions === | |
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 | 104 available. This is not necessary if your architecture has |
105 /dev/urandom defined. | |
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 | 110 extra performance. Check whether -DUSE_X86_32=0 is faster on 64-bit |
111 x86 architectures. | |
112 | |
113 If you want to replace the Mersenne Twister with another | |
114 generator then redefine randi32 appropriately. | |
115 | |
116 === Usage instructions === | |
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 | 119 |
120 All generators share the same state vector. | |
121 | |
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 | 139 |
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 | 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 | 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 | 155 */ |
156 | |
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 | 159 #endif |
160 | |
23662
bd77ab816e43
eliminate obsolete file lo-math.h
John W. Eaton <jwe@octave.org>
parents:
23649
diff
changeset
|
161 #include <cmath> |
21720
b28c26ae737c
use C++ for random number generator sources
John W. Eaton <jwe@octave.org>
parents:
21661
diff
changeset
|
162 #include <cstdio> |
31142
8245e773bb5b
randmtzig.cc: Add missing include <ctime> (bug #62750).
Stefan Brüns <stefan.bruens@rwth-aachen.de>
parents:
30564
diff
changeset
|
163 #include <ctime> |
5742 | 164 |
24856
8bb0251fcfde
Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents:
24534
diff
changeset
|
165 #include <algorithm> |
28597
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
166 #include <random> |
24856
8bb0251fcfde
Use a uint32 state vector for random number generators (bug #50256).
Rik <rik@octave.org>
parents:
24534
diff
changeset
|
167 |
28597
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
168 #include "oct-syscalls.h" |
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
169 #include "oct-time.h" |
5742 | 170 #include "randmtzig.h" |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
171 |
21066
258c787cd9ce
maint: Use "FIXME:" consistently in code base.
Rik <rik@octave.org>
parents:
20955
diff
changeset
|
172 /* FIXME: may want to suppress X86 if sizeof(long) > 4 */ |
20791
f7084eae3318
maint: Use Octave coding conventions for #if statements.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
173 #if ! defined (USE_X86_32) |
21202
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
174 # if defined (i386) || defined (HAVE_X86_32) |
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
175 # define USE_X86_32 1 |
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
176 # else |
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
177 # define USE_X86_32 0 |
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21066
diff
changeset
|
178 # endif |
5742 | 179 #endif |
180 | |
31605
e88a07dec498
maint: Use macros to begin/end C++ namespaces.
Rik <rik@octave.org>
parents:
31142
diff
changeset
|
181 OCTAVE_BEGIN_NAMESPACE(octave) |
e88a07dec498
maint: Use macros to begin/end C++ namespaces.
Rik <rik@octave.org>
parents:
31142
diff
changeset
|
182 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
183 /* ===== Mersenne Twister 32-bit generator ===== */ |
5742 | 184 |
185 #define MT_M 397 | |
186 #define MATRIX_A 0x9908b0dfUL /* constant vector a */ | |
187 #define UMASK 0x80000000UL /* most significant w-r bits */ | |
188 #define LMASK 0x7fffffffUL /* least significant r bits */ | |
189 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) | |
190 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) | |
191 | |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
192 static uint32_t *next; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
193 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
|
194 static int left = 1; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
195 static int initf = 0; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
196 static int initt = 1; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
197 static int inittf = 1; |
5742 | 198 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
199 /* 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
|
200 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
|
201 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
202 int j; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
203 state[0] = s & 0xffffffffUL; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
204 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
|
205 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
206 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
|
207 /* 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
|
208 /* 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
|
209 /* 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
|
210 /* 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
|
211 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
|
212 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
213 left = 1; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
214 initf = 1; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
215 } |
5742 | 216 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
217 /* 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
|
218 /* 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
|
219 /* 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
|
220 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
|
221 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
222 int i, j, k; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
223 init_mersenne_twister (19650218UL); |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
224 i = 1; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
225 j = 0; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
226 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
|
227 for (; k; k--) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
228 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
229 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
|
230 + 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
|
231 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
|
232 i++; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
233 j++; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
234 if (i >= MT_N) |
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 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
|
237 i = 1; |
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 if (j >= key_length) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
240 j = 0; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
241 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
242 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
|
243 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
244 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
|
245 - i; /* non linear */ |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
246 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
|
247 i++; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
248 if (i >= MT_N) |
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 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
|
251 i = 1; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
252 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
253 } |
5742 | 254 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
255 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
|
256 left = 1; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
257 initf = 1; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
258 } |
5742 | 259 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
260 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
|
261 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
262 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
|
263 int n = 0; |
5742 | 264 |
28597
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
265 // Gather some entropy from various sources |
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
266 |
27102
84ff9953faa1
where possible, eliminate octave:: namespace qualifier inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
267 sys::time now; |
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
268 |
28597
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
269 // Current time in seconds |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
270 if (n < MT_N) |
28597
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
271 entropy[n++] = now.unix_time (); |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
272 |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
273 // CPU time used (usec) |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
274 if (n < MT_N) |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
275 entropy[n++] = clock (); |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
276 |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
277 // Fractional part of current time |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
278 if (n < MT_N) |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
279 entropy[n++] = now.usec (); |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
280 |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
281 // Include the PID to make sure that several processes reaching here at the |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
282 // same time use different random numbers. |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
283 if (n < MT_N) |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
284 entropy[n++] = sys::getpid (); |
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
285 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
286 if (n < MT_N) |
28597
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
287 { |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
288 try |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
289 { |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
290 // The standard doesn't *guarantee* that random_device provides |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
291 // non-deterministic random numbers. So add entropy from this |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
292 // source last to make sure we gathered at least some entropy from |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
293 // the earlier sources. |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
294 std::random_device rd; |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
295 std::uniform_int_distribution<uint32_t> dist; |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
296 // Add 1024 bit of "true" entropy |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
297 int n_max = std::min (n + 32, MT_N); |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
298 while (n < n_max) |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
299 entropy[n++] = dist (rd); |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
300 } |
29163
8f67ad8b3103
maint: Updating naming conventions for exceptions and use const where possible.
Rik <rik@octave.org>
parents:
28597
diff
changeset
|
301 catch (const std::exception&) |
28597
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
302 { |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
303 // Just ignore any exception and skip that source of entropy. |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
304 } |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
305 } |
21928
315f4ba604c8
* randmtzig.cc: Use octave::sys::time.
John W. Eaton <jwe@octave.org>
parents:
21746
diff
changeset
|
306 |
28597
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
307 // Send all the entropy into the initial state vector |
a61b651914d6
Use C++11 functions to gather entropy from system (bug #58800).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27933
diff
changeset
|
308 init_mersenne_twister (entropy, n); |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
309 } |
5742 | 310 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
311 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
|
312 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
313 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
|
314 left = save[MT_N]; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
315 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
|
316 } |
5742 | 317 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
318 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
|
319 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
320 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
|
321 save[MT_N] = left; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
322 } |
5742 | 323 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
324 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
|
325 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
326 uint32_t *p = state; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
327 int j; |
5742 | 328 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
329 /* 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
|
330 /* 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
|
331 /* 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
|
332 /* 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
|
333 if (initf == 0) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
334 init_mersenne_twister (); |
5742 | 335 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
336 left = MT_N; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
337 next = state; |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
338 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
339 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
|
340 *p = p[MT_M] ^ TWIST(p[0], p[1]); |
5742 | 341 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
342 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
|
343 *p = p[MT_M-MT_N] ^ TWIST(p[0], p[1]); |
5742 | 344 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
345 *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
|
346 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
347 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
348 /* 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
|
349 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
|
350 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
351 uint32_t y; |
5742 | 352 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
353 if (--left == 0) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
354 next_state (); |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
355 y = *next++; |
5742 | 356 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
357 /* Tempering */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
358 y ^= (y >> 11); |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
359 y ^= (y << 7) & 0x9d2c5680UL; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
360 y ^= (y << 15) & 0xefc60000UL; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
361 return (y ^ (y >> 18)); |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
362 } |
5742 | 363 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
364 /* ===== Uniform generators ===== */ |
5742 | 365 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
366 /* Select which 32 bit generator to use */ |
5742 | 367 #define randi32 randmt |
368 | |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
369 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
|
370 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
371 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
|
372 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
|
373 #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
|
374 uint64_t u; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
375 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
|
376 p[0] = lo; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
377 p[1] = hi; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
378 return u; |
5742 | 379 #else |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
380 return ((static_cast<uint64_t> (hi) << 32) | lo); |
5742 | 381 #endif |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
382 } |
5742 | 383 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
384 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
|
385 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
386 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
|
387 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
|
388 #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
|
389 uint64_t u; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
390 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
|
391 p[0] = lo; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
392 p[1] = hi; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
393 return u; |
5742 | 394 #else |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
395 return ((static_cast<uint64_t> (hi) << 32) | lo); |
5742 | 396 #endif |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
397 } |
5742 | 398 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
399 /* 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
|
400 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
|
401 { |
27855
9405e2be91d0
Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents:
27102
diff
changeset
|
402 uint32_t i; |
9405e2be91d0
Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents:
27102
diff
changeset
|
403 |
9405e2be91d0
Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents:
27102
diff
changeset
|
404 do |
9405e2be91d0
Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents:
27102
diff
changeset
|
405 { |
9405e2be91d0
Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents:
27102
diff
changeset
|
406 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
|
407 } |
9405e2be91d0
Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents:
27102
diff
changeset
|
408 while (i == 0); |
9405e2be91d0
Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents:
27102
diff
changeset
|
409 |
9405e2be91d0
Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents:
27102
diff
changeset
|
410 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
|
411 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
412 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
413 /* 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
|
414 static double randu53 (void) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
415 { |
27858
797be8d10c22
Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents:
27857
diff
changeset
|
416 int32_t a, b; |
797be8d10c22
Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents:
27857
diff
changeset
|
417 |
797be8d10c22
Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents:
27857
diff
changeset
|
418 do |
797be8d10c22
Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents:
27857
diff
changeset
|
419 { |
797be8d10c22
Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents:
27857
diff
changeset
|
420 a = randi32 () >> 5; |
797be8d10c22
Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents:
27857
diff
changeset
|
421 b = randi32 () >> 6; |
797be8d10c22
Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents:
27857
diff
changeset
|
422 } |
797be8d10c22
Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents:
27857
diff
changeset
|
423 while (a == 0 && b == 0); |
797be8d10c22
Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents:
27857
diff
changeset
|
424 |
797be8d10c22
Create equally-spaced floating point values for rand ("double").
Rik <rik@octave.org>
parents:
27857
diff
changeset
|
425 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
|
426 } |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
427 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
428 /* 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
|
429 template <> |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29163
diff
changeset
|
430 OCTAVE_API double |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
431 rand_uniform<double> (void) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
432 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
433 return randu53 (); |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
434 } |
5742 | 435 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
436 /* 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
|
437 template <> |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29163
diff
changeset
|
438 OCTAVE_API float |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
439 rand_uniform<float> (void) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
440 { |
27855
9405e2be91d0
Restrict rand ("single") to range (0, 1) (bug #41742).
Rik <rik@octave.org>
parents:
27102
diff
changeset
|
441 return randu24 (); |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
442 } |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
443 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
444 /* ===== Ziggurat normal and exponential generators ===== */ |
5742 | 445 |
446 #define ZIGGURAT_TABLE_SIZE 256 | |
447 | |
448 #define ZIGGURAT_NOR_R 3.6541528853610088 | |
449 #define ZIGGURAT_NOR_INV_R 0.27366123732975828 | |
450 #define NOR_SECTION_AREA 0.00492867323399 | |
451 | |
452 #define ZIGGURAT_EXP_R 7.69711747013104972 | |
453 #define ZIGGURAT_EXP_INV_R 0.129918765548341586 | |
454 #define EXP_SECTION_AREA 0.0039496598225815571993 | |
455 | |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
456 #define 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
|
457 #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
|
458 #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
|
459 #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
|
460 #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
|
461 #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
|
462 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
463 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
|
464 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
|
465 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
|
466 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
|
467 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
468 /* |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
469 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
|
470 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
|
471 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
|
472 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
|
473 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
|
474 Ziggurat. This has several advantages. |
5742 | 475 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
476 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
|
477 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
|
478 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
|
479 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
|
480 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
|
481 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
|
482 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
|
483 4) Compile flag for full 53-bit random mantissa. |
5742 | 484 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
485 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
|
486 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
|
487 ziggurat, even if the code itself isn't... |
5742 | 488 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
489 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
|
490 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
|
491 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
|
492 -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
|
493 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
|
494 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
|
495 full 32 bits. We use EMANTISSA for this term. |
5742 | 496 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
497 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
|
498 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
|
499 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
|
500 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
|
501 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
|
502 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
|
503 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
|
504 */ |
5742 | 505 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
506 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
|
507 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
508 int i; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
509 double x, x1; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
510 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
511 /* 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
|
512 x1 = ZIGGURAT_NOR_R; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
513 wi[255] = x1 / NMANTISSA; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
514 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
|
515 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
516 /* 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
|
517 * defines this as |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
518 * 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
|
519 * 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
|
520 */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
521 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
|
522 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
|
523 fi[0] = 1.; |
5742 | 524 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
525 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
|
526 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
527 /* 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
|
528 * 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
|
529 */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
530 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
|
531 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
|
532 wi[i] = x / NMANTISSA; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
533 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
|
534 x1 = x; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
535 } |
5742 | 536 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
537 ki[1] = 0; |
5742 | 538 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
539 /* 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
|
540 x1 = ZIGGURAT_EXP_R; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
541 we[255] = x1 / EMANTISSA; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
542 fe[255] = exp (-x1); |
5742 | 543 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
544 /* 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
|
545 * defines this as |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
546 * 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
|
547 * 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
|
548 */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
549 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
|
550 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
|
551 fe[0] = 1.; |
5742 | 552 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
553 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
|
554 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
555 /* 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
|
556 * 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
|
557 */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
558 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
|
559 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
|
560 we[i] = x / EMANTISSA; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
561 fe[i] = exp (-x); |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
562 x1 = x; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
563 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
564 ke[1] = 0; |
5742 | 565 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
566 initt = 0; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
567 } |
5742 | 568 |
25433
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 * 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
|
571 * algorithm in their paper |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
572 * |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
573 * 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
|
574 * 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
|
575 * 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
|
576 * 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
|
577 * 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
|
578 * 5) goto step 1 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
579 * |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
580 * 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
|
581 * 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
|
582 */ |
5742 | 583 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
584 |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29163
diff
changeset
|
585 template <> OCTAVE_API double rand_normal<double> (void) |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
586 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
587 if (initt) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
588 create_ziggurat_tables (); |
5742 | 589 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
590 while (1) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
591 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
592 /* 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
|
593 * 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
|
594 * 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
|
595 * 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
|
596 * 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
|
597 * 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
|
598 * 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
|
599 */ |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21720
diff
changeset
|
600 # 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
|
601 /* 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
|
602 double x; |
30394
f3f3e3793fb5
maint: style check C++ files in liboctave/ ahead of 7.1 release.
Rik <rik@octave.org>
parents:
29359
diff
changeset
|
603 int si, idx; |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
604 uint32_t lo, hi; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
605 int64_t rabs; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
606 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
|
607 lo = randi32 (); |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
608 idx = lo & 0xFF; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
609 hi = randi32 (); |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
610 si = hi & UMASK; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
611 p[0] = lo; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
612 p[1] = hi & 0x1FFFFF; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
613 x = ( si ? -rabs : rabs ) * wi[idx]; |
22188 | 614 # else |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
615 /* 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
|
616 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
|
617 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
|
618 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
|
619 const double x = ( (r & 1) ? -rabs : rabs) * wi[idx]; |
22188 | 620 # endif |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
621 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
|
622 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
|
623 else if (idx == 0) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
624 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
625 /* 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
|
626 * |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
627 * 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
|
628 * 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
|
629 * 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
|
630 * 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
|
631 * 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
|
632 * |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
633 * [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
|
634 */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
635 double xx, yy; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
636 do |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
637 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
638 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
|
639 yy = - std::log (RANDU); |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
640 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
641 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
|
642 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
|
643 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
644 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
|
645 return x; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
646 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
647 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
648 |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29163
diff
changeset
|
649 template <> OCTAVE_API double rand_exponential<double> (void) |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
650 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
651 if (initt) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
652 create_ziggurat_tables (); |
5742 | 653 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
654 while (1) |
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 ZIGINT ri = ERANDI; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
657 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
|
658 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
|
659 if (ri < ke[idx]) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
660 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
|
661 else if (idx == 0) |
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 /* 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
|
664 * |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
665 * 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
|
666 * x = r - ln(U); |
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 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
|
669 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
670 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
|
671 return x; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
672 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
673 } |
5742 | 674 |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29163
diff
changeset
|
675 template <> OCTAVE_API void rand_uniform<double> (octave_idx_type n, double *p) |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
676 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
677 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
|
678 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
679 |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29163
diff
changeset
|
680 template <> OCTAVE_API void rand_normal (octave_idx_type n, double *p) |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
681 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
682 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
|
683 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
684 |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29163
diff
changeset
|
685 template <> OCTAVE_API void rand_exponential (octave_idx_type n, double *p) |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
686 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
687 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
|
688 } |
5742 | 689 |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
690 #undef ZIGINT |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
691 #undef EMANTISSA |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
692 #undef ERANDI |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
693 #undef NMANTISSA |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
694 #undef NRANDI |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
695 #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
|
696 |
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
697 #define 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
|
698 #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
|
699 #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
|
700 #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
|
701 #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
|
702 #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
|
703 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
704 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
|
705 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
|
706 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
|
707 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
|
708 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
709 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
|
710 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
711 int i; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
712 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
|
713 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
714 /* 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
|
715 x1 = ZIGGURAT_NOR_R; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
716 fwi[255] = x1 / NMANTISSA; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
717 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
|
718 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
719 /* 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
|
720 * defines this as |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
721 * 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
|
722 * 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
|
723 */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
724 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
|
725 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
|
726 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
|
727 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
728 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
|
729 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
730 /* 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
|
731 * 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
|
732 */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
733 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
|
734 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
|
735 fwi[i] = x / NMANTISSA; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
736 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
|
737 x1 = x; |
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 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
740 fki[1] = 0; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
741 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
742 /* 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
|
743 x1 = ZIGGURAT_EXP_R; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
744 fwe[255] = x1 / EMANTISSA; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
745 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
|
746 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
747 /* 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
|
748 * defines this as |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
749 * 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
|
750 * 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
|
751 */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
752 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
|
753 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
|
754 ffe[0] = 1.; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
755 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
756 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
|
757 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
758 /* 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
|
759 * 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
|
760 */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
761 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
|
762 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
|
763 fwe[i] = x / EMANTISSA; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
764 ffe[i] = exp (-x); |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
765 x1 = x; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
766 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
767 fke[1] = 0; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
768 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
769 inittf = 0; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
770 } |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
771 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
772 /* |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
773 * 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
|
774 * algorithm in their paper |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
775 * |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
776 * 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
|
777 * 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
|
778 * 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
|
779 * 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
|
780 * 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
|
781 * 5) goto step 1 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
782 * |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
783 * 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
|
784 * 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
|
785 */ |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
786 |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29163
diff
changeset
|
787 template <> OCTAVE_API float rand_normal<float> (void) |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
788 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
789 if (inittf) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
790 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
|
791 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
792 while (1) |
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 /* 32-bit mantissa */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
795 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
|
796 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
|
797 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
|
798 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
|
799 if (rabs < fki[idx]) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
800 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
|
801 else if (idx == 0) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
802 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
803 /* 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
|
804 * |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
805 * 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
|
806 * 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
|
807 * 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
|
808 * 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
|
809 * 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
|
810 * |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
811 * [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
|
812 */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
813 float xx, yy; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
814 do |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
815 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
816 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
|
817 yy = - std::log (RANDU); |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
818 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
819 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
|
820 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
|
821 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
822 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
|
823 return x; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
824 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
825 } |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
826 |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29163
diff
changeset
|
827 template <> OCTAVE_API float rand_exponential<float> (void) |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
828 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
829 if (inittf) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
830 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
|
831 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
832 while (1) |
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 ZIGINT ri = ERANDI; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
835 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
|
836 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
|
837 if (ri < fke[idx]) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
838 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
|
839 else if (idx == 0) |
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 /* 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
|
842 * |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
843 * 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
|
844 * x = r - ln(U); |
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 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
|
847 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
848 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
|
849 return x; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
850 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
851 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
852 |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29163
diff
changeset
|
853 template <> OCTAVE_API void rand_uniform (octave_idx_type n, float *p) |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
854 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
855 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
|
856 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
857 |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29163
diff
changeset
|
858 template <> OCTAVE_API void rand_normal (octave_idx_type n, float *p) |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
859 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
860 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
|
861 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
862 |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29163
diff
changeset
|
863 template <> OCTAVE_API void rand_exponential (octave_idx_type n, float *p) |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
864 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
865 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
|
866 } |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
867 |
31605
e88a07dec498
maint: Use macros to begin/end C++ namespaces.
Rik <rik@octave.org>
parents:
31142
diff
changeset
|
868 OCTAVE_END_NAMESPACE(octave) |
e88a07dec498
maint: Use macros to begin/end C++ namespaces.
Rik <rik@octave.org>
parents:
31142
diff
changeset
|
869 |