Mercurial > hg > octave-lyh
diff src/rand.cc @ 1:78fd87e624cb
[project @ 1993-08-08 01:13:40 by jwe]
Initial revision
author | jwe |
---|---|
date | Sun, 08 Aug 1993 01:13:40 +0000 |
parents | |
children | d68036bcad4c |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/src/rand.cc @@ -0,0 +1,260 @@ +// tc-rand.cc -*- C++ -*- +/* + +Copyright (C) 1993 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, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#ifdef __GNUG__ +#pragma implementation +#endif + +#include "tree-const.h" +#include "f77-uscore.h" +#include "error.h" +#include "utils.h" + +// Possible distributions of random numbers. +enum rand_dist { uniform, normal }; + +// Current distribution of random numbers. +static rand_dist current_distribution = uniform; + +extern "C" +{ + int *F77_FCN (dgennor) (double*, double*, double*); + int *F77_FCN (dgenunf) (double*, double*, double*); + int *F77_FCN (setall) (int*, int*); + int *F77_FCN (getsd) (int*, int*); +} + +#ifdef WITH_DLD +tree_constant * +builtin_rand_2 (tree_constant *args, int nargin, int nargout) +{ + return rand_internal (args, nargin, nargout); +} +#endif + +static double +curr_rand_seed (void) +{ + union d2i { double d; int i[2]; }; + union d2i u; + F77_FCN (getsd) (&(u.i[0]), &(u.i[1])); + return u.d; +} + +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; +} + +static void +set_rand_seed (double val) +{ + union d2i { double d; int i[2]; }; + union d2i u; + u.d = val; + int i0 = force_to_fit_range (u.i[0], 1, 2147483563); + int i1 = force_to_fit_range (u.i[1], 1, 2147483399); + F77_FCN (setall) (&i0, &i1); +} + +static char * +curr_rand_dist (void) +{ + if (current_distribution == uniform) + return "uniform"; + else if (current_distribution == normal) + return "normal"; + else + { + panic_impossible (); + return (char *) NULL; + } +} + +tree_constant * +rand_internal (tree_constant *args, int nargin, int nargout) +{ +// Assumes that we have been given the correct number of arguments. + + tree_constant *retval = NULL_TREE_CONST; + + static int initialized = 0; + if (! initialized) + { +// Make the random number generator give us a different sequence every +// time we start octave unless we specifically set the seed. +#if 0 + int s0 = 1234567890; + int s1 = 123456789; +#else + time_t now; + struct tm *tm; + + time (&now); + tm = localtime (&now); + + int s0 = tm->tm_min * 60 + tm->tm_sec; + int s1 = (tm->tm_mday - 1) * 24 * 3600 + tm->tm_hour * 3600 + s0; +#endif + s0 = force_to_fit_range (s0, 1, 2147483563); + s1 = force_to_fit_range (s1, 1, 2147483399); + + F77_FCN (setall) (&s0, &s1); + initialized = 1; + } + + int n = 0; + int m = 0; + if (nargin == 1) + { + n = 1; + m = 1; + goto gen_matrix; + } + else if (nargin == 2) + { + switch (args[1].const_type ()) + { + case tree_constant_rep::string_constant: + char *s_arg = args[1].string_value (); + if (strcmp (s_arg, "dist") == 0) + { + retval = new tree_constant [2]; + char *s = curr_rand_dist (); + retval[0] = tree_constant (s); + } + else if (strcmp (s_arg, "seed") == 0) + { + retval = new tree_constant [2]; + double d = curr_rand_seed (); + retval[0] = tree_constant (d); + } + else if (strcmp (s_arg, "uniform") == 0) + current_distribution = uniform; + else if (strcmp (s_arg, "normal") == 0) + current_distribution = normal; + else + { + delete [] retval; + retval = NULL_TREE_CONST; + message ("rand", "unrecognized string argument"); + } + break; + case tree_constant_rep::scalar_constant: + case tree_constant_rep::complex_scalar_constant: + n = NINT (args[1].double_value ()); + m = n; + goto gen_matrix; + case tree_constant_rep::range_constant: + { + Range r = args[1].range_value (); + n = 1; + m = NINT (r.nelem ()); + } + goto gen_matrix; + case tree_constant_rep::matrix_constant: + case tree_constant_rep::complex_matrix_constant: + n = NINT (args[1].rows ()); + m = NINT (args[1].columns ()); + goto gen_matrix; + default: + panic_impossible (); + break; + } + } + else if (nargin == 3) + { + if (args[1].is_string_type () + && strcmp (args[1].string_value (), "seed") == 0) + { + double d = args[2].to_scalar (); + set_rand_seed (d); + } + else + { + n = NINT (args[1].to_scalar ()); + m = NINT (args[2].to_scalar ()); + goto gen_matrix; + } + } + + return retval; + + gen_matrix: + + if (n == 0 || m == 0) + { + retval = new tree_constant [2]; + Matrix m (0, 0); + retval[0] = tree_constant (m); + } + else if (n > 0 && m > 0) + { + retval = new tree_constant [2]; + Matrix rand_mat (n, m); + for (int j = 0; j < m; j++) + for (int i = 0; i < n; i++) + { + double d_zero = 0.0; + double d_one = 1.0; + double val; + switch (current_distribution) + { + case uniform: + F77_FCN (dgenunf) (&d_zero, &d_one, &val); + rand_mat.elem (i, j) = val; + break; + case normal: + F77_FCN (dgennor) (&d_zero, &d_one, &val); + rand_mat.elem (i, j) = val; + break; + default: + panic_impossible (); + break; + } + } + + retval[0] = tree_constant (rand_mat); + } + else + message ("rand", "invalid negative argument"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/