octave-nkf: liboctave/oct-rand.cc annotate

annotate liboctave/oct-rand.cc @ 8710:739141cde75a ss-3-1-52

fix typo in Array-f.cc

author	Jaroslav Hajek <highegg@gmail.com>
date	Mon, 09 Feb 2009 21:51:31 +0100
parents	25bc2d31e1bf
children	eb63fbe60fab

rev	line source
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	1 /*
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	2
7017 a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe] jwe parents: 7016 diff changeset	3 Copyright (C) 2003, 2005, 2006, 2007 John W. Eaton
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	4
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	5 This file is part of Octave.
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	6
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	7 Octave is free software; you can redistribute it and/or modify it
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	8 under the terms of the GNU General Public License as published by the
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6435 diff changeset	9 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6435 diff changeset	10 option) any later version.
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	11
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	12 Octave is distributed in the hope that it will be useful, but WITHOUT
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	15 for more details.
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	16
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	17 You should have received a copy of the GNU General Public License
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6435 diff changeset	18 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6435 diff changeset	19 <http://www.gnu.org/licenses/>.
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	20
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	21 */
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	22
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	23 #ifdef HAVE_CONFIG_H
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	24 #include <config.h>
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	25 #endif
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	26
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	27 #include <map>
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	28 #include <vector>
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	29
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	30 #include "f77-fcn.h"
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	31 #include "lo-ieee.h"
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	32 #include "lo-error.h"
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	33 #include "lo-mappers.h"
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	34 #include "oct-rand.h"
4415 ed8c4aaa8648 [project @ 2003-05-16 21:20:33 by jwe] jwe parents: 4308 diff changeset	35 #include "oct-time.h"
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	36 #include "data-conv.h"
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	37 #include "randmtzig.h"
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	38 #include "randpoisson.h"
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	39 #include "randgamma.h"
6435 71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	40 #include "mach-info.h"
8377 25bc2d31e1bf improve OCTAVE_LOCAL_BUFFER Jaroslav Hajek <highegg@gmail.com> parents: 8124 diff changeset	41 #include "oct-locbuf.h"
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	42
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	43 extern "C"
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	44 {
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	45 F77_RET_T
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	46 F77_FUNC (dgennor, DGENNOR) (const double&, const double&, double&);
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	47
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	48 F77_RET_T
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	49 F77_FUNC (dgenunf, DGENUNF) (const double&, const double&, double&);
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	50
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	51 F77_RET_T
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	52 F77_FUNC (dgenexp, DGENEXP) (const double&, double&);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	53
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	54 F77_RET_T
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	55 F77_FUNC (dignpoi, DIGNPOI) (const double&, double&);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	56
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	57 F77_RET_T
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	58 F77_FUNC (dgengam, DGENGAM) (const double&, const double&, double&);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	59
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	60 F77_RET_T
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4552 diff changeset	61 F77_FUNC (setall, SETALL) (const octave_idx_type&, const octave_idx_type&);
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	62
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	63 F77_RET_T
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4552 diff changeset	64 F77_FUNC (getsd, GETSD) (octave_idx_type&, octave_idx_type&);
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	65
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	66 F77_RET_T
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4552 diff changeset	67 F77_FUNC (setsd, SETSD) (const octave_idx_type&, const octave_idx_type&);
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	68
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	69 F77_RET_T
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4552 diff changeset	70 F77_FUNC (setcgn, SETCGN) (const octave_idx_type&);
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	71 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	72
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	73 octave_rand *octave_rand::instance = 0;
6326 129ef2ae319f [project @ 2007-02-18 14:01:27 by dbateman] dbateman parents: 5828 diff changeset	74
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	75 octave_rand::octave_rand (void)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	76 : current_distribution (uniform_dist), use_old_generators (false),
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	77 rand_states ()
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	78 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	79 initialize_ranlib_generators ();
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	80
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	81 initialize_mersenne_twister ();
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	82 }
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	83
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	84 bool
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	85 octave_rand::instance_ok (void)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	86 {
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	87 bool retval = true;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	88
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	89 if (! instance)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	90 instance = new octave_rand ();
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	91
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	92 if (! instance)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	93 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	94 (*current_liboctave_error_handler)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	95 ("unable to create octave_rand object!");
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	96
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	97 retval = false;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	98 }
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	99
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	100 return retval;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	101 }
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	102
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	103 double
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	104 octave_rand::do_seed (void)
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	105 {
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4552 diff changeset	106 union d2i { double d; octave_idx_type i[2]; };
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	107 union d2i u;
6435 71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	108
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	109 oct_mach_info::float_format ff = oct_mach_info::native_float_format ();
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	110
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	111 switch (ff)
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	112 {
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	113 case oct_mach_info::flt_fmt_ieee_big_endian:
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	114 F77_FUNC (getsd, GETSD) (u.i[1], u.i[0]);
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	115 break;
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	116 default:
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	117 F77_FUNC (getsd, GETSD) (u.i[0], u.i[1]);
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	118 break;
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	119 }
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	120
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	121 return u.d;
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	122 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	123
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	124 static octave_idx_type
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	125 force_to_fit_range (octave_idx_type i, octave_idx_type lo, octave_idx_type hi)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	126 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	127 assert (hi > lo && lo >= 0 && hi > lo);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	128
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	129 i = i > 0 ? i : -i;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	130
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	131 if (i < lo)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	132 i = lo;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	133 else if (i > hi)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	134 i = i % hi;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	135
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	136 return i;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	137 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	138
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	139 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	140 octave_rand::do_seed (double s)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	141 {
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	142 use_old_generators = true;
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	143
6435 71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	144 int i0, i1;
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4552 diff changeset	145 union d2i { double d; octave_idx_type i[2]; };
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	146 union d2i u;
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	147 u.d = s;
6435 71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	148
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	149 oct_mach_info::float_format ff = oct_mach_info::native_float_format ();
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	150
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	151 switch (ff)
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	152 {
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	153 case oct_mach_info::flt_fmt_ieee_big_endian:
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	154 i1 = force_to_fit_range (u.i[0], 1, 2147483563);
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	155 i0 = force_to_fit_range (u.i[1], 1, 2147483399);
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	156 break;
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	157 default:
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	158 i0 = force_to_fit_range (u.i[0], 1, 2147483563);
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	159 i1 = force_to_fit_range (u.i[1], 1, 2147483399);
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	160 break;
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	161 }
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	162
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	163 F77_FUNC (setsd, SETSD) (i0, i1);
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	164 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	165
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	166 ColumnVector
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	167 octave_rand::do_state (const std::string& d)
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	168 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	169 return rand_states[d.empty () ? current_distribution : get_dist_id (d)];
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	170 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	171
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	172 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	173 octave_rand::do_state (const ColumnVector& s, const std::string& d)
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	174 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	175 use_old_generators = false;
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	176
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	177 int old_dist = current_distribution;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	178
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	179 int new_dist = d.empty () ? current_distribution : get_dist_id (d);
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	180
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	181 ColumnVector saved_state;
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	182
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	183 if (old_dist != new_dist)
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	184 saved_state = get_internal_state ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	185
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	186 set_internal_state (s);
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	187
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	188 rand_states[new_dist] = get_internal_state ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	189
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	190 if (old_dist != new_dist)
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	191 rand_states[old_dist] = saved_state;
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	192 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	193
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	194 std::string
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	195 octave_rand::do_distribution (void)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	196 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	197 std::string retval;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	198
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	199 switch (current_distribution)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	200 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	201 case uniform_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	202 retval = "uniform";
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	203 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	204
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	205 case normal_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	206 retval = "normal";
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	207 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	208
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	209 case expon_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	210 retval = "exponential";
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	211 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	212
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	213 case poisson_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	214 retval = "poisson";
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	215 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	216
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	217 case gamma_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	218 retval = "gamma";
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	219 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	220
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	221 default:
7536 4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	222 (*current_liboctave_error_handler)
4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	223 ("rand: invalid distribution ID = %d", current_distribution);
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	224 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	225 }
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	226
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	227 return retval;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	228 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	229
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	230 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	231 octave_rand::do_distribution (const std::string& d)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	232 {
7536 4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	233 int id = get_dist_id (d);
4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	234
4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	235 switch (id)
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	236 {
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	237 case uniform_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	238 octave_rand::uniform_distribution ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	239 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	240
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	241 case normal_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	242 octave_rand::normal_distribution ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	243 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	244
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	245 case expon_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	246 octave_rand::exponential_distribution ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	247 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	248
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	249 case poisson_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	250 octave_rand::poisson_distribution ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	251 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	252
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	253 case gamma_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	254 octave_rand::gamma_distribution ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	255 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	256
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	257 default:
7536 4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	258 (*current_liboctave_error_handler)
4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	259 ("rand: invalid distribution ID = %d", id);
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	260 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	261 }
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	262 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	263
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	264 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	265 octave_rand::do_uniform_distribution (void)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	266 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	267 switch_to_generator (uniform_dist);
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	268
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	269 F77_FUNC (setcgn, SETCGN) (uniform_dist);
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	270 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	271
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	272 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	273 octave_rand::do_normal_distribution (void)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	274 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	275 switch_to_generator (normal_dist);
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	276
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	277 F77_FUNC (setcgn, SETCGN) (normal_dist);
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	278 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	279
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	280 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	281 octave_rand::do_exponential_distribution (void)
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	282 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	283 switch_to_generator (expon_dist);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	284
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	285 F77_FUNC (setcgn, SETCGN) (expon_dist);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	286 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	287
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	288 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	289 octave_rand::do_poisson_distribution (void)
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	290 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	291 switch_to_generator (poisson_dist);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	292
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	293 F77_FUNC (setcgn, SETCGN) (poisson_dist);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	294 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	295
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	296 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	297 octave_rand::do_gamma_distribution (void)
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	298 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	299 switch_to_generator (gamma_dist);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	300
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	301 F77_FUNC (setcgn, SETCGN) (gamma_dist);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	302 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	303
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	304
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	305 double
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	306 octave_rand::do_scalar (double a)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	307 {
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	308 double retval = 0.0;
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	309
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	310 if (use_old_generators)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	311 {
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	312 switch (current_distribution)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	313 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	314 case uniform_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	315 F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	316 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	317
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	318 case normal_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	319 F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	320 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	321
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	322 case expon_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	323 F77_FUNC (dgenexp, DGENEXP) (1.0, retval);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	324 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	325
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	326 case poisson_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	327 if (a < 0.0 \|\| xisnan(a) \|\| xisinf(a))
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	328 retval = octave_NaN;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	329 else
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	330 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	331 // workaround bug in ignpoi, by calling with different Mu
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	332 F77_FUNC (dignpoi, DIGNPOI) (a + 1, retval);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	333 F77_FUNC (dignpoi, DIGNPOI) (a, retval);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	334 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	335 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	336
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	337 case gamma_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	338 if (a <= 0.0 \|\| xisnan(a) \|\| xisinf(a))
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	339 retval = octave_NaN;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	340 else
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	341 F77_FUNC (dgengam, DGENGAM) (1.0, a, retval);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	342 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	343
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	344 default:
7536 4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	345 (*current_liboctave_error_handler)
4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	346 ("rand: invalid distribution ID = %d", current_distribution);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	347 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	348 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	349 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	350 else
4543 79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	351 {
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	352 switch (current_distribution)
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	353 {
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	354 case uniform_dist:
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	355 retval = oct_randu ();
4543 79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	356 break;
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	357
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	358 case normal_dist:
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	359 retval = oct_randn ();
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	360 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	361
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	362 case expon_dist:
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	363 retval = oct_rande ();
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	364 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	365
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	366 case poisson_dist:
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	367 retval = oct_randp (a);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	368 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	369
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	370 case gamma_dist:
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	371 retval = oct_randg (a);
4543 79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	372 break;
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	373
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	374 default:
7536 4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	375 (*current_liboctave_error_handler)
4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	376 ("rand: invalid distribution ID = %d", current_distribution);
4543 79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	377 break;
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	378 }
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	379
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	380 save_state ();
4543 79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	381 }
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	382
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	383 return retval;
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	384 }
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	385
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	386 Matrix
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	387 octave_rand::do_matrix (octave_idx_type n, octave_idx_type m, double a)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	388 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	389 Matrix retval;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	390
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	391 if (n >= 0 && m >= 0)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	392 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	393 retval.resize (n, m);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	394
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	395 if (n > 0 && m > 0)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	396 fill (retval.capacity(), retval.fortran_vec(), a);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	397 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	398 else
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	399 (*current_liboctave_error_handler) ("rand: invalid negative argument");
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	400
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	401 return retval;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	402 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	403
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	404 NDArray
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	405 octave_rand::do_nd_array (const dim_vector& dims, double a)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	406 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	407 NDArray retval;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	408
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	409 if (! dims.all_zero ())
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	410 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	411 retval.resize (dims);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	412
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	413 fill (retval.capacity(), retval.fortran_vec(), a);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	414 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	415
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	416 return retval;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	417 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	418
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	419 Array<double>
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	420 octave_rand::do_vector (octave_idx_type n, double a)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	421 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	422 Array<double> retval;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	423
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	424 if (n > 0)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	425 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	426 retval.resize (n);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	427
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	428 fill (retval.capacity (), retval.fortran_vec (), a);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	429 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	430 else if (n < 0)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	431 (*current_liboctave_error_handler) ("rand: invalid negative argument");
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	432
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	433 return retval;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	434 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	435
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	436 // Make the random number generator give us a different sequence every
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	437 // time we start octave unless we specifically set the seed. The
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	438 // technique used below will cycle monthly, but it it does seem to
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	439 // work ok to give fairly different seeds each time Octave starts.
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	440
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	441 void
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	442 octave_rand::initialize_ranlib_generators (void)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	443 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	444 octave_localtime tm;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	445 int stored_distribution = current_distribution;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	446 F77_FUNC (setcgn, SETCGN) (uniform_dist);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	447
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	448 int hour = tm.hour() + 1;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	449 int minute = tm.min() + 1;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	450 int second = tm.sec() + 1;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	451
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	452 octave_idx_type s0 = tm.mday() * hour * minute * second;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	453 octave_idx_type s1 = hour * minute * second;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	454
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	455 s0 = force_to_fit_range (s0, 1, 2147483563);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	456 s1 = force_to_fit_range (s1, 1, 2147483399);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	457
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	458 F77_FUNC (setall, SETALL) (s0, s1);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	459 F77_FUNC (setcgn, SETCGN) (stored_distribution);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	460 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	461
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	462 void
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	463 octave_rand::initialize_mersenne_twister (void)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	464 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	465 oct_init_by_entropy ();
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	466
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	467 ColumnVector s = get_internal_state ();
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	468
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	469 rand_states[uniform_dist] = s;
8124 d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	470
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	471 oct_init_by_entropy ();
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	472 s = get_internal_state ();
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	473 rand_states[normal_dist] = s;
8124 d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	474
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	475 oct_init_by_entropy ();
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	476 s = get_internal_state ();
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	477 rand_states[expon_dist] = s;
8124 d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	478
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	479 oct_init_by_entropy ();
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	480 s = get_internal_state ();
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	481 rand_states[poisson_dist] = s;
8124 d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	482
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	483 oct_init_by_entropy ();
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	484 s = get_internal_state ();
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	485 rand_states[gamma_dist] = s;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	486 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	487
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	488 ColumnVector
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	489 octave_rand::get_internal_state (void)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	490 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	491 ColumnVector s (MT_N + 1);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	492
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	493 OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	494
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	495 oct_get_state (tmp);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	496
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	497 for (octave_idx_type i = 0; i <= MT_N; i++)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	498 s.elem (i) = static_cast<double> (tmp [i]);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	499
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	500 return s;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	501 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	502
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	503 void
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	504 octave_rand::save_state (void)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	505 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	506 rand_states[current_distribution] = get_internal_state ();;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	507 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	508
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	509 int
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	510 octave_rand::get_dist_id (const std::string& d)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	511 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	512 int retval = unknown_dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	513
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	514 if (d == "uniform" \|\| d == "rand")
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	515 retval = uniform_dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	516 else if (d == "normal" \|\| d == "randn")
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	517 retval = normal_dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	518 else if (d == "exponential" \|\| d == "rande")
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	519 retval = expon_dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	520 else if (d == "poisson" \|\| d == "randp")
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	521 retval = poisson_dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	522 else if (d == "gamma" \|\| d == "randg")
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	523 retval = gamma_dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	524 else
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	525 (*current_liboctave_error_handler)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	526 ("rand: invalid distribution `%s'", d.c_str ());
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	527
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	528 return retval;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	529 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	530
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	531 void
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	532 octave_rand::set_internal_state (const ColumnVector& s)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	533 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	534 octave_idx_type len = s.length ();
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	535 octave_idx_type n = len < MT_N + 1 ? len : MT_N + 1;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	536
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	537 OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	538
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	539 for (octave_idx_type i = 0; i < n; i++)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	540 tmp[i] = static_cast<uint32_t> (s.elem(i));
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	541
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	542 if (len == MT_N + 1 && tmp[MT_N] <= MT_N && tmp[MT_N] > 0)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	543 oct_set_state (tmp);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	544 else
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	545 oct_init_by_array (tmp, len);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	546 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	547
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	548 void
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	549 octave_rand::switch_to_generator (int dist)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	550 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	551 if (dist != current_distribution)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	552 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	553 current_distribution = dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	554
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	555 set_internal_state (rand_states[dist]);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	556 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	557 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	558
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	559 #define MAKE_RAND(len) \
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	560 do \
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	561 { \
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	562 double val; \
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	563 for (volatile octave_idx_type i = 0; i < len; i++) \
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	564 { \
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	565 OCTAVE_QUIT; \
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	566 RAND_FUNC (val); \
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	567 v[i] = val; \
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	568 } \
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	569 } \
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	570 while (0)
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	571
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	572 void
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	573 octave_rand::fill (octave_idx_type len, double *v, double a)
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	574 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	575 if (len < 1)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	576 return;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	577
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	578 switch (current_distribution)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	579 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	580 case uniform_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	581 if (use_old_generators)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	582 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	583 #define RAND_FUNC(x) F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, x)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	584 MAKE_RAND (len);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	585 #undef RAND_FUNC
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	586 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	587 else
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	588 oct_fill_randu (len, v);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	589 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	590
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	591 case normal_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	592 if (use_old_generators)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	593 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	594 #define RAND_FUNC(x) F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, x)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	595 MAKE_RAND (len);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	596 #undef RAND_FUNC
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	597 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	598 else
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	599 oct_fill_randn (len, v);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	600 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	601
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	602 case expon_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	603 if (use_old_generators)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	604 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	605 #define RAND_FUNC(x) F77_FUNC (dgenexp, DGENEXP) (1.0, x)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	606 MAKE_RAND (len);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	607 #undef RAND_FUNC
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	608 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	609 else
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	610 oct_fill_rande (len, v);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	611 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	612
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	613 case poisson_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	614 if (use_old_generators)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	615 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	616 if (a < 0.0 \|\| xisnan(a) \|\| xisinf(a))
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	617 #define RAND_FUNC(x) x = octave_NaN;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	618 MAKE_RAND (len);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	619 #undef RAND_FUNC
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	620 else
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	621 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	622 // workaround bug in ignpoi, by calling with different Mu
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	623 double tmp;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	624 F77_FUNC (dignpoi, DIGNPOI) (a + 1, tmp);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	625 #define RAND_FUNC(x) F77_FUNC (dignpoi, DIGNPOI) (a, x)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	626 MAKE_RAND (len);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	627 #undef RAND_FUNC
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	628 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	629 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	630 else
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	631 oct_fill_randp (a, len, v);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	632 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	633
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	634 case gamma_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	635 if (use_old_generators)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	636 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	637 if (a <= 0.0 \|\| xisnan(a) \|\| xisinf(a))
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	638 #define RAND_FUNC(x) x = octave_NaN;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	639 MAKE_RAND (len);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	640 #undef RAND_FUNC
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	641 else
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	642 #define RAND_FUNC(x) F77_FUNC (dgengam, DGENGAM) (1.0, a, x)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	643 MAKE_RAND (len);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	644 #undef RAND_FUNC
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	645 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	646 else
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	647 oct_fill_randg (a, len, v);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	648 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	649
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	650 default:
7536 4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	651 (*current_liboctave_error_handler)
4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	652 ("rand: invalid distribution ID = %d", current_distribution);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	653 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	654 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	655
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	656 save_state ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	657
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	658 return;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	659 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	660
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	661 /*
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	662 ;;; Local Variables: ***
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	663 ;;; mode: C++ ***
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	664 ;;; End: ***
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	665 */

Mercurial > octave-nkf

annotate liboctave/oct-rand.cc @ 8710:739141cde75a ss-3-1-52