From octave-maintainers-request at bevo dot che dot wisc dot edu Sat Jan 24 18:25:48 2004 Subject: Re: randn benchmarks From: Paul Kienzle To: octave-maintainers at bevo dot che dot wisc dot edu Cc: "Dmitri A. Sergatskov" Date: Sat, 24 Jan 2004 19:24:36 -0500 --Apple-Mail-21-46312021 Content-Transfer-Encoding: 7bit Content-Type: text/plain; charset=US-ASCII; format=flowed On Jan 24, 2004, at 3:32 PM, Paul Kienzle wrote: > Comparing apples to apples (randmt to rand instead > of randn), I find that octave-forge's rand is faster on > both IRIX and Linux. I can squeeze another 25% out > of it by directly accessing the array values. I will need > to do this anyway for N-d array support. I will post this > code in a bit (it is on another machine which just stopped > responding). The attached code uses the randn generator suggested earlier by David (not the ziggurat algorithm): while(1) { const double u = randu(); const double v = 1.7156 * (randu() - 0.5); const double x = u - 0.449871; const double y = fabs(v) + 0.386595; const double q = x*x + y*(0.196*y - 0.25472*x); if (q < 0.27597) return v/u; if ((q <= 0.27846) && (v*v < - 4. * log(u) *u*u)) return v/u; } Anyone know if the constants are valid for doubles? They should be if I understand the algorithm (although there may be a slight loss in efficiency over the ideal constants). The time for rand(1000) on my machine is now 0.17 rather than 0.25, and randn(1000) is 0.60 rather than 2.5. It also creates n-d arrays. It is not well tested. I should also take the opportunity to add in the fast poisson generator that I have, but that will have to wait until next pass. Paul Kienzle pkienzle at users dot sf dot net --Apple-Mail-21-46312021 Content-Transfer-Encoding: 7bit Content-Type: application/octet-stream; x-unix-mode=0644; name="rand.cc" Content-Disposition: attachment; filename=rand.cc /* Copyright (C) 2004 Paul Kienzle This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. Octave is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with Octave; see the file COPYING. If not, write to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. // Based on John Eaton's rand.cc and Dirk Eddelbuettel's randmt.cc and the Mersenne Twister // Copyright (C) 1996, 1997 John W. Eaton // Copyright (C) 1998, 1999 Dirk Eddelbuettel // Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, // // $Id: rand.cc,v 1.1.1.1 2001/10/10 19:54:49 pkienzle Exp $ */ #include #include #include #ifdef HAVE_GETTIMEOFDAY #include #endif /* A C-program for MT19937, with initialization improved 2002/2/10. Coded by Takuji Nishimura and Makoto Matsumoto. This is a faster version by taking Shawn Cokus's optimization, Matthe Bellew's simplification, Isaku Wada's real version. Before using, initialize the state by using init_genrand(seed) or init_by_array(init_key, key_length). Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. The names of its contributors may not be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. Any feedback is very welcome. http://www.math.keio.ac.jp/matumoto/emt.html email: matumoto at math dot keio dot ac dot jp * 2004-01-19 Paul Kienzle * * comment out main * * add init_by_entropy, get_state, set_state * * converted to allow compiling by C++ compiler */ /* Period parameters */ #define N 624 #define M 397 #define MATRIX_A 0x9908b0dfUL /* constant vector a */ #define UMASK 0x80000000UL /* most significant w-r bits */ #define LMASK 0x7fffffffUL /* least significant r bits */ #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) static unsigned long state[N]; /* the array for the state vector */ static int left = 1; static int initf = 0; static unsigned long *next; /* initializes state[N] with a seed */ void init_genrand(unsigned long s) { int j; state[0]= s & 0xffffffffUL; for (j=1; j> 30)) + j); /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ /* In the previous versions, MSBs of the seed affect */ /* only MSBs of the array state[]. */ /* 2002/01/09 modified by Makoto Matsumoto */ state[j] &= 0xffffffffUL; /* for >32 bit machines */ } left = 1; initf = 1; } /* initialize by an array with array-length */ /* init_key is the array for initializing keys */ /* key_length is its length */ void init_by_array(unsigned long init_key[], int key_length) { int i, j, k; init_genrand(19650218UL); i=1; j=0; k = (N>key_length ? N : key_length); for (; k; k--) { state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL)) + init_key[j] + j; /* non linear */ state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ i++; j++; if (i>=N) { state[0] = state[N-1]; i=1; } if (j>=key_length) j=0; } for (k=N-1; k; k--) { state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL)) - i; /* non linear */ state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ i++; if (i>=N) { state[0] = state[N-1]; i=1; } } state[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ left = 1; initf = 1; } void init_by_entropy(void) { unsigned long entropy[N]; int n = 0; /* Look for entropy in /dev/urandom */ FILE* urandom =fopen("/dev/urandom", "rb"); if (urandom) { while (n> 11); y ^= (y << 7) & 0x9d2c5680UL; y ^= (y << 15) & 0xefc60000UL; y ^= (y >> 18); return y; } /* generates a random number on (0,1)-real-interval */ inline double randu(void) { return ((double)randi() + 0.5) * (1.0/4294967296.0); /* divided by 2^32 */ } /* generates a random number on [0,1) with 53-bit resolution*/ inline double rand53(void) { unsigned long a=randi()>>5, b=randi()>>6; return(a*67108864.0+b)*(1.0/9007199254740992.0); } /* These real versions are due to Isaku Wada, 2002/01/09 added */ inline double randn(void) { while(1) { const double u = randu(); const double v = 1.7156 * (randu() - 0.5); const double x = u - 0.449871; const double y = fabs(v) + 0.386595; const double q = x*x + y*(0.196*y - 0.25472*x); if (q < 0.27597) return v/u; if ((q <= 0.27846) && (v*v < - 4. * log(u) *u*u)) return v/u; } } // Octave interface starts here static octave_value do_seed (octave_value_list args) { octave_value retval; // Check if they said the magic words std::string s_arg = args(0).string_value (); if (s_arg == "seed") { // If they ask for the current "seed", then reseed with the next // available random number unsigned long a = randi(); init_genrand(a); retval = (double)a; } else if (s_arg == "state") { unsigned long state[N+1]; get_state(state); RowVector a(N+1); for (int i=0; i < N+1; i++) a(i) = state[i]; retval = a; } else { error ("rand: unrecognized string argument"); return retval; } // Check if just getting state if (args.length() == 1) return retval; // Set the state from either a scalar or a previously returned state vector octave_value tmp = args(1); if (tmp.is_scalar_type ()) { unsigned long n = (unsigned long)(tmp.double_value()); if (! error_state) init_genrand(n); } else if (tmp.is_matrix_type () && tmp.rows() == N+1 && tmp.columns() == 1) { Array a(tmp.vector_value ()); if (! error_state) { unsigned long state[N+1]; for (int i = 0; i < N+1; i++) state[i] = (unsigned long)(a(i)); set_state(state); } } else error ("rand: not a state vector"); return retval; } #ifdef HAVE_ND_ARRAYS static void do_size(const char *fcn, octave_value_list args, dim_vector& dims) { int nargin = args.length(); if (nargin == 1) { get_dimensions(args(0),fcn,dims); } else { dims.resize (nargin); for (int i = 0; i < nargin; i++) { dims(i) = args(i).is_empty () ? 0 : args(i).nint_value (); if (error_state) return; } } int ndim = dims.length(); while (ndim > 2 && dims(ndim-1) == 1) ndim--; dims.resize (ndim); check_dimensions (dims, fcn); } #else static void do_size(octave_value_list args, int& nr, int& nc) { int nargin = args.length(); if (nargin == 0) { nr = nc = 1; } else if (nargin == 1) { octave_value tmp = args(0); if (tmp.is_scalar_type ()) { double dval = tmp.double_value (); if (xisnan (dval)) { error ("rand: NaN is invalid a matrix dimension"); } else { nr = nc = NINT (tmp.double_value ()); } } else if (tmp.is_range ()) { Range rng = tmp.range_value (); nr = 1; nc = rng.nelem (); } else if (tmp.is_matrix_type ()) { // XXX FIXME XXX -- this should probably use the function // from data.cc. Matrix a = args(0).matrix_value (); if (error_state) return; nr = a.rows (); nc = a.columns (); if (nr == 1 && nc == 2) { nr = NINT (a (0, 0)); nc = NINT (a (0, 1)); } else if (nr == 2 && nc == 1) { nr = NINT (a (0, 0)); nc = NINT (a (1, 0)); } else warning ("rand (A): use rand (size (A)) instead"); } else { gripe_wrong_type_arg ("rand", tmp); } } else if (nargin == 2) { double rval = args(0).double_value (); double cval = args(1).double_value (); if (! error_state) { if (xisnan (rval) || xisnan (cval)) { error ("rand: NaN is invalid as a matrix dimension"); } else { nr = NINT (rval); nc = NINT (cval); } } } } #endif void fill_randu(int n, double *p) { for (int i=0; i < n; i++) p[i] = randu(); } DEFUN_DLD (rand, args, nargout, "-*- texinfo -*-\n\ at deftypefn {Loadable Function} {} rand (@var{x})\n\ at deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\ at deftypefnx {Loadable Function} {@var{v} =} rand (\"state\", @var{x})\n\ at deftypefnx {Loadable Function} {@var{s} =} rand (\"seed\", @var{x})\n\ Return a matrix with random elements uniformly distributed on the\n\ semi-open interval [0, 1). The arguments are handled the same as the\n\ arguments for at code{eye} dot \n\ \n\ You can reset the state of the random number generator using the\n\ form\n\ \n\ at example\n\ v = rand (\"state\", x)\n\ at end example\n\ \n\ where at var{x} is a scalar value. This returns the current state\n\ of the random number generator in the column vector at var{v} dot If\n\ at var{x} is not given, then the state is returned but not changed.\n\ Later, you can restore the random number generator to the state at var{v}\n\ using the form\n\ \n\ at example\n\ u = rand (\"state\", v)\n\ at end example\n\ \n\ at noindent\n\ If instead of \"state\" you use \"seed\" to query the random\n\ number generator, then the state will be collapsed from roughly\n\ 20000 bits down to 32 bits.\n\ \n\ at code{rand} uses the Mersenne Twister, with a period of 2^19937-1.\n\ Do NOT use for CRYPTOGRAPHY without securely hashing several returned\n\ values together, otherwise the generator state can be learned after\n\ reading 624 consecutive values.\n\ \n\ M. Matsumoto and T. Nishimura, ``Mersenne Twister: A 623-dimensionally\n\ equidistributed uniform pseudorandom number generator'', ACM Trans. on\n\ Modeling and Computer Simulation Vol. 8, No. 1, Januray pp.3-30 1998\n\ \n\ http://www.math.keio.ac.jp/~matumoto/emt.html\n\ at end deftypefn\n\ at seealso{randn}\n") { octave_value_list retval; // list of return values int nargin = args.length (); // number of arguments supplied if (nargin > 0 && args(0).is_string()) retval(0) = do_seed (args); else { #ifdef HAVE_ND_ARRAYS dim_vector dims; do_size ("rand", args, dims); if (error_state) return retval; int ndim = dims.length(); switch (ndim) { case 0: { double v; fill_randu(1,&v); retval(0) = v; } break; case 1: case 2: { Matrix X(dims(0),dims(ndim==1?0:1)); fill_randu(X.capacity(),X.fortran_vec()); retval(0) = X; } break; default: { NDArray Xn(dims); fill_randu(Xn.capacity(),Xn.fortran_vec()); retval(0) = Xn; } break; } #else int nr=0, nc=0; do_size (args, nr, nc); if (! error_state) { Matrix X(nr, nc); fill_randu(nr*nc,X.fortran_vec()); retval(0) = X; } #endif } return retval; } void fill_randn(int n, double *p) { for (int i=0; i < n; i++) p[i] = randn(); } DEFUN_DLD (randn, args, nargout, "-*- texinfo -*-\n\ at deftypefn {Loadable Function} {} randn (@var{x})\n\ at deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\ at deftypefnx {Loadable Function} {@var{v} =} randn (\"state\", @var{x})\n\ at deftypefnx {Loadable Function} {@var{s} =} randn (\"seed\", @var{x})\n\ Return a matrix with normally distributed random elements. The\n\ arguments are handled the same as the arguments for at code{rand} dot \n\ \n\ at code{randn} uses Ahrens and Dieter (1973) to transform from U to N(0,1).\n\n\ at end deftypefn\n\ at seealso{rand}\n") { octave_value_list retval; // list of return values int nargin = args.length (); // number of arguments supplied if (nargin > 0 && args(0).is_string()) retval(0) = do_seed (args); else { #ifdef HAVE_ND_ARRAYS dim_vector dims; do_size ("randn", args, dims); if (error_state) return retval; int ndim = dims.length(); switch (ndim) { case 0: { double v; fill_randn(1,&v); retval(0) = v; } break; case 1: case 2: { Matrix X(dims(0),dims(ndim==1?0:1)); fill_randn(X.capacity(),X.fortran_vec()); retval(0) = X; } break; default: { NDArray Xn(dims); fill_randn(Xn.capacity(),Xn.fortran_vec()); retval(0) = Xn; } break; } #else int nr=0, nc=0; do_size (args, nr, nc); if (! error_state) { Matrix X(nr, nc); fill_randn(nr*nc,X.fortran_vec()); retval(0) = X; } #endif } return retval; } /* ;;; Local Variables: *** ;;; mode: C++ *** ;;; End: *** */ --Apple-Mail-21-46312021--