/* Copyright (C) 1996, 1997 John W. Eaton This file is part of Octave. Octave 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. */ // Octave modules for the Mersenne Twister (MT199337) Random Number Generator // using Richard J. Wagner's C++ implementaion MersenneTwister.h // // This file provides two Octave functions: // rand for Uniform random number // randn for Normal random number // // Based on John Eaton's rand.cc and Dirk Eddelbuettel's randmt.cc // Copyright (C) 1996, 1997 John W. Eaton // Copyright (C) 1998, 1999 Dirk Eddelbuettel // // 2001-02-10 Paul Kienzle // * use Richard J. Wagner's MersenneTwister.h // * use John Eaton's rand.cc for similar interface // * add rand("state") to obtain complete state // * renamed snorm to randn, and removed the `static' keyword from // the variables in randu after thoroughly checking that they were // reset each time the function is entered. // $Id$ #include #include #include "MersenneTwister.h" static MTRand randu; // The following routine to transform U(0,UINT_MAX) into N(0,1) is from // the GNU GPL'ed randlib library by Brown, Lovato, Russell and Venier // which is available from ftp://odin.mdacc.tmc.edu/pub/source double randn(void) /* ********************************************************************** (STANDARD-) N O R M A L DISTRIBUTION ********************************************************************** ********************************************************************** FOR DETAILS SEE: AHRENS, J.H. AND DIETER, U. EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM SAMPLING FROM THE NORMAL DISTRIBUTION. MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937. ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL' (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION) Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of SUNIF. The argument IR thus goes away. Modified by Dirk Eddelbuettel on 1 Aug 1999 to use randu() instead of RANF ********************************************************************** THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE */ { const double a[32] = { 0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904, 0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322, 0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818, 1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594, 1.862732,2.153875 }; const double d[31] = { 0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243, 0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094, 0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791, 0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039 }; const double t[31] = { 7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3, 1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2, 2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2, 4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2, 9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031 }; const double h[31] = { 3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2, 4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2, 4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2, 5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2, 8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474 }; int i; double u, s, snorm, ustar, aa, w, y, tt; u = randu(); s = 0.0; if(u > 0.5) s = 1.0; u += (u-s); u = 32.0*u; i = (int) (u); if(i == 32) i = 31; if(i == 0) goto S100; /* START CENTER */ ustar = u-(double)i; aa = *(a+i-1); S40: if(ustar <= *(t+i-1)) goto S60; w = (ustar-*(t+i-1))**(h+i-1); S50: /* EXIT (BOTH CASES) */ y = aa+w; snorm = y; if(s == 1.0) snorm = -y; return snorm; S60: /* CENTER CONTINUED */ u = randu(); w = u*(*(a+i)-aa); tt = (0.5*w+aa)*w; goto S80; S70: tt = u; ustar = randu(); S80: if(ustar > tt) goto S50; u = randu(); if(ustar >= u) goto S70; ustar = randu(); goto S40; S100: /* START TAIL */ i = 6; aa = *(a+31); goto S120; S110: aa += *(d+i-1); i += 1; S120: u += u; if(u < 1.0) goto S110; u -= 1.0; S140: w = u**(d+i-1); tt = (0.5*w+aa)*w; goto S160; S150: tt = u; S160: ustar = randu(); if(ustar > tt) goto S50; u = randu(); if(ustar >= u) goto S150; u = randu(); goto S140; } // 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 MTRand::uint32 a = randu.randInt(); randu.seed(a); retval = (double)a; } else if (s_arg == "state") { MTRand::uint32 state[randu.SAVE]; randu.save(state); RowVector a(randu.SAVE); for (int i=0; i < randu.SAVE; 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 ()) { MTRand::uint32 n = MTRand::uint32(tmp.double_value()); if (! error_state) randu.seed(n); } else if (tmp.is_matrix_type () && tmp.rows() == randu.SAVE && tmp.columns() == 1) { Array a(tmp.vector_value ()); if (! error_state) { MTRand::uint32 state[randu.SAVE]; for (int i = 0; i < randu.SAVE; i++) state[i] = MTRand::uint32(a(i)); randu.load(state); } } else error ("rand: not a state vector"); return retval; } 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); } } } } DEFUN_DLD (rand, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {} rand (@var{x})\n\ @deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\ @deftypefnx {Loadable Function} {@var{v} =} rand (\"state\", @var{x})\n\ @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 @code{eye}.\n\ \n\ You can reset the state of the random number generator using the\n\ form\n\ \n\ @example\n\ v = rand (\"state\", x)\n\ @end example\n\ \n\ where @var{x} is a scalar value. This returns the current state\n\ of the random number generator in the column vector @var{v}. If\n\ @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 @var{v}\n\ using the form\n\ \n\ @example\n\ u = rand (\"state\", v)\n\ @end example\n\ \n\ @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\ @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\ @end deftypefn\n\ @seealso{randn}\n") { octave_value_list retval; // list of return values int nargin = args.length (); // number of arguments supplied if (nargin > 2) print_usage("rand"); else if (nargin > 0 && args(0).is_string()) retval(0) = do_seed (args); else { int nr=0, nc=0; do_size (args, nr, nc); if (! error_state) { Matrix X(nr, nc); for (int c=0; c < nc; c++) for (int r=0; r < nr; r++) X(r,c) = randu.randExc(); retval(0) = X; } } return retval; } DEFUN_DLD (randn, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {} randn (@var{x})\n\ @deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\ @deftypefnx {Loadable Function} {@var{v} =} randn (\"state\", @var{x})\n\ @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 @code{rand}.\n\ \n\ @code{randn} uses Ahrens and Dieter (1973) to transform from U to N(0,1).\n\n\ @end deftypefn\n\ @seealso{rand}\n") { octave_value_list retval; // list of return values int nargin = args.length (); // number of arguments supplied if (nargin > 2) print_usage("randn"); else if (nargin > 0 && args(0).is_string()) retval(0) = do_seed (args); else { int nr=0, nc=0; do_size (args, nr, nc); if (! error_state) { Matrix X(nr, nc); for (int c=0; c < nc; c++) for (int r=0; r < nr; r++) X(r,c) = randn(); retval(0) = X; } } return retval; } /* ;;; Local Variables: *** ;;; mode: C++ *** ;;; End: *** */