Mercurial > hg > octave-terminal
changeset 4308:b738d1a02adb
[project @ 2003-01-24 19:37:12 by jwe]
author | jwe |
---|---|
date | Fri, 24 Jan 2003 19:37:12 +0000 |
parents | fd034cd46aea |
children | a9560cebae6e |
files | liboctave/oct-rand.cc liboctave/oct-rand.h |
diffstat | 2 files changed, 371 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/liboctave/oct-rand.cc @@ -0,0 +1,304 @@ +/* + +Copyright (C) 2003 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. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "f77-fcn.h" +#include "lo-error.h" +#include "oct-rand.h" + +// Possible distributions of random numbers. This was handled with an +// enum, but unwind_protecting that doesn't work so well. +#define uniform_dist 1 +#define normal_dist 2 + +// Current distribution of random numbers. +static int current_distribution = uniform_dist; + +// Has the seed been set yet? +static bool initialized = false; + +extern "C" +{ + int F77_FUNC (dgennor, DGENNOR) (const double&, const double&, + double&); + + int F77_FUNC (dgenunf, DGENUNF) (const double&, const double&, + double&); + + int F77_FUNC (setall, SETALL) (const int&, const int&); + + int F77_FUNC (getsd, GETSD) (int&, int&); + + int F77_FUNC (setsd, SETSD) (const int&, const int&); + + int F77_FUNC (setcgn, SETCGN) (const int&); +} + +static int +force_to_fit_range (int i, int lo, int hi) +{ + assert (hi > lo && lo >= 0 && hi > lo); + + i = i > 0 ? i : -i; + + if (i < lo) + i = lo; + else if (i > hi) + i = i % hi; + + return i; +} + +// Make the random number generator give us a different sequence every +// time we start octave unless we specifically set the seed. The +// technique used below will cycle monthly, but it it does seem to +// work ok to give fairly different seeds each time Octave starts. + +static void +do_initialization (void) +{ + time_t now; + struct tm *tm; + + time (&now); + tm = localtime (&now); + + int hour = tm->tm_hour + 1; + int minute = tm->tm_min + 1; + int second = tm->tm_sec + 1; + + int s0 = tm->tm_mday * hour * minute * second; + int s1 = hour * minute * second; + + s0 = force_to_fit_range (s0, 1, 2147483563); + s1 = force_to_fit_range (s1, 1, 2147483399); + + F77_FUNC (setall, SETALL) (s0, s1); + + initialized = true; +} + +static inline void +maybe_initialize (void) +{ + if (! initialized) + do_initialization (); +} + +double +octave_rand::seed (void) +{ + maybe_initialize (); + + union d2i { double d; int i[2]; }; + union d2i u; + F77_FUNC (getsd, GETSD) (u.i[0], u.i[1]); + return u.d; +} + +void +octave_rand::seed (double s) +{ + maybe_initialize (); + + union d2i { double d; int i[2]; }; + union d2i u; + u.d = s; + int i0 = force_to_fit_range (u.i[0], 1, 2147483563); + int i1 = force_to_fit_range (u.i[1], 1, 2147483399); + F77_FUNC (setsd, SETSD) (i0, i1); +} + +std::string +octave_rand::distribution (void) +{ + maybe_initialize (); + + if (current_distribution == uniform_dist) + return "uniform"; + else if (current_distribution == normal_dist) + return "normal"; + else + { + abort (); + return ""; + } +} + +void +octave_rand::distribution (const std::string& d) +{ + if (d == "uniform") + octave_rand::uniform_distribution (); + else if (d == "normal") + octave_rand::normal_distribution (); + else + (*current_liboctave_error_handler) ("rand: invalid distribution"); +} + +void +octave_rand::uniform_distribution (void) +{ + maybe_initialize (); + + current_distribution = uniform_dist; + + F77_FUNC (setcgn, SETCGN) (uniform_dist); +} + +void +octave_rand::normal_distribution (void) +{ + maybe_initialize (); + + current_distribution = normal_dist; + + F77_FUNC (setcgn, SETCGN) (normal_dist); +} + +double +octave_rand::scalar (void) +{ + maybe_initialize (); + + double retval = 0.0; + + switch (current_distribution) + { + case uniform_dist: + F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval); + break; + + case normal_dist: + F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval); + break; + + default: + abort (); + break; + } + + return retval; +} + +#define MAKE_RAND_MAT(mat, nr, nc, f, F) \ + do \ + { \ + double val; \ + for (volatile int j = 0; j < nc; j++) \ + for (volatile int i = 0; i < nr; i++) \ + { \ + OCTAVE_QUIT; \ + F77_FUNC (f, F) (0.0, 1.0, val); \ + mat(i,j) = val; \ + } \ + } \ + while (0) + +Matrix +octave_rand::matrix (int n, int m) +{ + maybe_initialize (); + + Matrix retval; + + if (n >= 0 && m >= 0) + { + retval.resize (n, m); + + if (n > 0 && m > 0) + { + switch (current_distribution) + { + case uniform_dist: + MAKE_RAND_MAT (retval, n, m, dgenunf, DGENUNF); + break; + + case normal_dist: + MAKE_RAND_MAT (retval, n, m, dgennor, DGENNOR); + break; + + default: + abort (); + break; + } + } + } + else + (*current_liboctave_error_handler) ("rand: invalid negative argument"); + + return retval; +} + +#define MAKE_RAND_ARRAY(array, n, f, F) \ + do \ + { \ + double val; \ + for (volatile int i = 0; i < n; i++) \ + { \ + OCTAVE_QUIT; \ + F77_FUNC (f, F) (0.0, 1.0, val); \ + array(i) = val; \ + } \ + } \ + while (0) + +Array<double> +octave_rand::vector (int n) +{ + maybe_initialize (); + + Array<double> retval; + + if (n > 0) + { + retval.resize (n); + + switch (current_distribution) + { + case uniform_dist: + MAKE_RAND_ARRAY (retval, n, dgenunf, DGENUNF); + break; + + case normal_dist: + MAKE_RAND_ARRAY (retval, n, dgennor, DGENNOR); + break; + + default: + abort (); + break; + } + } + else if (n < 0) + (*current_liboctave_error_handler) ("rand: invalid negative argument"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/oct-rand.h @@ -0,0 +1,67 @@ +/* + +Copyright (C) 2003 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. + +*/ + +#if !defined (octave_rand_h) +#define octave_rand_h 1 + +#include <string> + +#include "dMatrix.h" + +struct +octave_rand +{ + // Return the current seed. + static double seed (void); + + // Set the seed. + static void seed (double s); + + // Return the current distribution. + static std::string distribution (void); + + // Set the current distribution. May be either "uniform" (the + // default) or "normal". + static void distribution (const std::string& d); + + static void uniform_distribution (void); + + static void normal_distribution (void); + + // Return the next number from the sequence. + static double scalar (void); + + // Return a matrix of numbers from the sequence, filled in column + // major order. + static Matrix matrix (int r, int c); + + // Return an array of numbers from the sequence. + static Array<double> vector (int n); +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/