octave-nkf: liboctave/randmtzig.c annotate

annotate liboctave/randmtzig.c @ 14655:43db83eff9db

Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) * oct-rand.cc (octave_rand:do_matrix): Remove method. (octave_rand::do_float_scalar, octave_rand::do_float_nd_array, octave_rand::do_float_vector): New methods. * oct-rand.h (octave_rand:do_matrix, octave_rand::matrix): Remove methods. (octave_rand::do_float_scalar, octave_rand::do_float_nd_array, octave_rand::do_float_vector, octave_rand::float_scalar, octave_rand::do_nd_array, octave_rand::float_vector): Declare new methods. * randgamma.c (oct_fill_float_randg, oct_float_randg): New functions. * randgamma.h (oct_fill_float_randg, oct_float_randg): Declare new functions. * randpoisson.c (oct_fill_float_randp, oct_float_randp): New functions. (poisson_cdf_lookup_float, poisson_rejection_float): New static functions. * randpoisson.h (oct_fill_float_randp, oct_float_randp): Declare new functions. * randmtzig.c (randu64) : Remove static function. (ALLBITS): Remove compile flag logic. (randu32): Change return type to float. (oct_float_randu, oct_float_randn, oct_float_rande, oct_fill_float_randu, oct_fill_float_randn, oct_fill_float_rande): New functions. (create_ziggurat_float_tables): New static function. * randmtzig.h (oct_float_randu, oct_float_randn, oct_float_rande, oct_fill_float_randu, oct_fill_float_randn, oct_fill_float_rande): Declare new functions. * rand.cc (do_rand): Add logic to parse "double" or "single" as final argument.

author	David Bateman <dbateman@free.fr>
date	Sat, 19 May 2012 00:14:59 +0200
parents	72c96de7a403
children	460a3c6d8bf1

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

Mercurial > octave-nkf

annotate liboctave/randmtzig.c @ 14655:43db83eff9db