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

annotate liboctave/oct-rand.cc @ 10350:12884915a8e4

merge MArray classes & improve Array interface

author	Jaroslav Hajek <highegg@gmail.com>
date	Sat, 23 Jan 2010 21:41:03 +0100
parents	6a50aa2cdf00
children	92a85ed5b86e

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
8920 eb63fbe60fab update copyright notices John W. Eaton <jwe@octave.org> parents: 8377 diff changeset	3 Copyright (C) 2003, 2005, 2006, 2007, 2008 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
10331 6a50aa2cdf00 never compile ranlib with 8-byte integers John W. Eaton <jwe@octave.org> parents: 10314 diff changeset	30 #include <stdint.h>
6a50aa2cdf00 never compile ranlib with 8-byte integers John W. Eaton <jwe@octave.org> parents: 10314 diff changeset	31
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	32 #include "f77-fcn.h"
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	33 #include "lo-ieee.h"
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	34 #include "lo-error.h"
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	35 #include "lo-mappers.h"
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	36 #include "oct-rand.h"
4415 ed8c4aaa8648 [project @ 2003-05-16 21:20:33 by jwe] jwe parents: 4308 diff changeset	37 #include "oct-time.h"
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	38 #include "data-conv.h"
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	39 #include "randmtzig.h"
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	40 #include "randpoisson.h"
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	41 #include "randgamma.h"
6435 71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	42 #include "mach-info.h"
8377 25bc2d31e1bf improve OCTAVE_LOCAL_BUFFER Jaroslav Hajek <highegg@gmail.com> parents: 8124 diff changeset	43 #include "oct-locbuf.h"
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	44
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	45 extern "C"
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	46 {
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	47 F77_RET_T
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	48 F77_FUNC (dgennor, DGENNOR) (const double&, const double&, double&);
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	49
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	50 F77_RET_T
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	51 F77_FUNC (dgenunf, DGENUNF) (const double&, const double&, double&);
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	52
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	53 F77_RET_T
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	54 F77_FUNC (dgenexp, DGENEXP) (const double&, double&);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	55
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	56 F77_RET_T
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	57 F77_FUNC (dignpoi, DIGNPOI) (const double&, double&);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	58
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	59 F77_RET_T
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	60 F77_FUNC (dgengam, DGENGAM) (const double&, const double&, double&);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	61
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	62 F77_RET_T
10331 6a50aa2cdf00 never compile ranlib with 8-byte integers John W. Eaton <jwe@octave.org> parents: 10314 diff changeset	63 F77_FUNC (setall, SETALL) (const int32_t&, const int32_t&);
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	64
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	65 F77_RET_T
10331 6a50aa2cdf00 never compile ranlib with 8-byte integers John W. Eaton <jwe@octave.org> parents: 10314 diff changeset	66 F77_FUNC (getsd, GETSD) (int32_t&, int32_t&);
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	67
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	68 F77_RET_T
10331 6a50aa2cdf00 never compile ranlib with 8-byte integers John W. Eaton <jwe@octave.org> parents: 10314 diff changeset	69 F77_FUNC (setsd, SETSD) (const int32_t&, const int32_t&);
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	70
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4543 diff changeset	71 F77_RET_T
10331 6a50aa2cdf00 never compile ranlib with 8-byte integers John W. Eaton <jwe@octave.org> parents: 10314 diff changeset	72 F77_FUNC (setcgn, SETCGN) (const int32_t&);
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	73 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: 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::instance = 0;
6326 129ef2ae319f [project @ 2007-02-18 14:01:27 by dbateman] dbateman parents: 5828 diff changeset	76
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	77 octave_rand::octave_rand (void)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	78 : 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	79 rand_states ()
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	80 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	81 initialize_ranlib_generators ();
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	82
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	83 initialize_mersenne_twister ();
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	84 }
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	85
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	86 bool
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	87 octave_rand::instance_ok (void)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	88 {
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	89 bool retval = true;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	90
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	91 if (! instance)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	92 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	93
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	94 if (! instance)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	95 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	96 (*current_liboctave_error_handler)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	97 ("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	98
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	99 retval = false;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	100 }
7533 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 return retval;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	103 }
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	104
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	105 double
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	106 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	107 {
10331 6a50aa2cdf00 never compile ranlib with 8-byte integers John W. Eaton <jwe@octave.org> parents: 10314 diff changeset	108 union d2i { double d; int32_t i[2]; };
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	109 union d2i u;
6435 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 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	112
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	113 switch (ff)
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	114 {
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	115 case oct_mach_info::flt_fmt_ieee_big_endian:
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	116 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	117 break;
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	118 default:
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	119 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	120 break;
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	121 }
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	122
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	123 return u.d;
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	124 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	125
10331 6a50aa2cdf00 never compile ranlib with 8-byte integers John W. Eaton <jwe@octave.org> parents: 10314 diff changeset	126 static int32_t
6a50aa2cdf00 never compile ranlib with 8-byte integers John W. Eaton <jwe@octave.org> parents: 10314 diff changeset	127 force_to_fit_range (int32_t i, int32_t lo, int32_t hi)
7537 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 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	130
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	131 i = i > 0 ? i : -i;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	132
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	133 if (i < lo)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	134 i = lo;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	135 else if (i > hi)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	136 i = i % hi;
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 return i;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	139 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	140
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	141 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	142 octave_rand::do_seed (double s)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	143 {
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	144 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	145
6435 71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	146 int i0, i1;
10331 6a50aa2cdf00 never compile ranlib with 8-byte integers John W. Eaton <jwe@octave.org> parents: 10314 diff changeset	147 union d2i { double d; int32_t i[2]; };
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	148 union d2i u;
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	149 u.d = s;
6435 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 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	152
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	153 switch (ff)
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	154 {
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	155 case oct_mach_info::flt_fmt_ieee_big_endian:
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	156 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	157 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	158 break;
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	159 default:
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	160 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	161 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	162 break;
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	163 }
71888a0ab040 [project @ 2007-03-23 14:17:02 by dbateman] dbateman parents: 6326 diff changeset	164
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	165 F77_FUNC (setsd, SETSD) (i0, i1);
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	166 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	167
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	168 ColumnVector
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	169 octave_rand::do_state (const std::string& d)
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	170 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	171 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	172 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	173
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	174 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	175 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	176 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	177 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	178
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	179 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	180
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	181 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	182
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	183 ColumnVector saved_state;
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	184
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	185 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	186 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	187
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	188 set_internal_state (s);
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 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	191
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	192 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	193 rand_states[old_dist] = saved_state;
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	194 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	195
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	196 std::string
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	197 octave_rand::do_distribution (void)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	198 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	199 std::string retval;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	200
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	201 switch (current_distribution)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	202 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	203 case uniform_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	204 retval = "uniform";
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	205 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	206
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	207 case normal_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	208 retval = "normal";
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	209 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	210
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	211 case expon_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	212 retval = "exponential";
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	213 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	214
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	215 case poisson_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	216 retval = "poisson";
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	217 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	218
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	219 case gamma_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	220 retval = "gamma";
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	221 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	222
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	223 default:
7536 4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	224 (*current_liboctave_error_handler)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	225 ("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	226 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	227 }
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	228
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	229 return retval;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	230 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	231
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	232 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	233 octave_rand::do_distribution (const std::string& d)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	234 {
7536 4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	235 int id = get_dist_id (d);
4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	236
4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	237 switch (id)
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	238 {
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	239 case uniform_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	240 octave_rand::uniform_distribution ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	241 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	242
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	243 case normal_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	244 octave_rand::normal_distribution ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	245 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	246
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	247 case expon_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	248 octave_rand::exponential_distribution ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	249 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	250
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	251 case poisson_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	252 octave_rand::poisson_distribution ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	253 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	254
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	255 case gamma_dist:
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	256 octave_rand::gamma_distribution ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	257 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	258
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	259 default:
7536 4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	260 (*current_liboctave_error_handler)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	261 ("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	262 break;
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	263 }
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	264 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	265
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	266 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	267 octave_rand::do_uniform_distribution (void)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	268 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	269 switch_to_generator (uniform_dist);
4308 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 F77_FUNC (setcgn, SETCGN) (uniform_dist);
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	272 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	273
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	274 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	275 octave_rand::do_normal_distribution (void)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	276 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	277 switch_to_generator (normal_dist);
4308 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 F77_FUNC (setcgn, SETCGN) (normal_dist);
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	280 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	281
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	282 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	283 octave_rand::do_exponential_distribution (void)
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	284 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	285 switch_to_generator (expon_dist);
5730 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 F77_FUNC (setcgn, SETCGN) (expon_dist);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	288 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	289
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	290 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	291 octave_rand::do_poisson_distribution (void)
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	292 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	293 switch_to_generator (poisson_dist);
5730 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 F77_FUNC (setcgn, SETCGN) (poisson_dist);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	296 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	297
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	298 void
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	299 octave_rand::do_gamma_distribution (void)
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	300 {
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	301 switch_to_generator (gamma_dist);
5730 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 F77_FUNC (setcgn, SETCGN) (gamma_dist);
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	304 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	305
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	306
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	307 double
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	308 octave_rand::do_scalar (double a)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	309 {
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	310 double retval = 0.0;
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 if (use_old_generators)
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	313 {
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	314 switch (current_distribution)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	315 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	316 case uniform_dist:
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	317 F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	318 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	319
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	320 case normal_dist:
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	321 F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	322 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	323
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	324 case expon_dist:
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	325 F77_FUNC (dgenexp, DGENEXP) (1.0, retval);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	326 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	327
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	328 case poisson_dist:
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	329 if (a < 0.0 \|\| xisnan(a) \|\| xisinf(a))
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	330 retval = octave_NaN;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	331 else
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	332 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	333 // workaround bug in ignpoi, by calling with different Mu
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	334 F77_FUNC (dignpoi, DIGNPOI) (a + 1, retval);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	335 F77_FUNC (dignpoi, DIGNPOI) (a, retval);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	336 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	337 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	338
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	339 case gamma_dist:
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	340 if (a <= 0.0 \|\| xisnan(a) \|\| xisinf(a))
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	341 retval = octave_NaN;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	342 else
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	343 F77_FUNC (dgengam, DGENGAM) (1.0, a, retval);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	344 break;
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	345
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	346 default:
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	347 (*current_liboctave_error_handler)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	348 ("rand: invalid distribution ID = %d", current_distribution);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	349 break;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	350 }
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	351 }
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	352 else
4543 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 switch (current_distribution)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	355 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	356 case uniform_dist:
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	357 retval = oct_randu ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	358 break;
4543 79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	359
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	360 case normal_dist:
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	361 retval = oct_randn ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	362 break;
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	363
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	364 case expon_dist:
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	365 retval = oct_rande ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	366 break;
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	367
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	368 case poisson_dist:
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	369 retval = oct_randp (a);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	370 break;
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	371
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	372 case gamma_dist:
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	373 retval = oct_randg (a);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	374 break;
4543 79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	375
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	376 default:
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	377 (*current_liboctave_error_handler)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	378 ("rand: invalid distribution ID = %d", current_distribution);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	379 break;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	380 }
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	381
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	382 save_state ();
4543 79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	383 }
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 return retval;
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	386 }
79df15d4470c [project @ 2003-10-18 03:53:52 by jwe] jwe parents: 4415 diff changeset	387
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	388 Matrix
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	389 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	390 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	391 Matrix retval;
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 if (n >= 0 && m >= 0)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	394 {
9647 54f45f883a53 optimize & extend randperm Jaroslav Hajek <highegg@gmail.com> parents: 8920 diff changeset	395 retval.clear (n, m);
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	396
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	397 if (n > 0 && m > 0)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	398 fill (retval.capacity(), retval.fortran_vec(), a);
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	399 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	400 else
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	401 (*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	402
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	403 return retval;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	404 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	405
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	406 NDArray
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	407 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	408 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	409 NDArray retval;
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 if (! dims.all_zero ())
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	412 {
9647 54f45f883a53 optimize & extend randperm Jaroslav Hajek <highegg@gmail.com> parents: 8920 diff changeset	413 retval.clear (dims);
7537 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 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	416 }
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 return retval;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	419 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	420
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	421 Array<double>
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	422 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	423 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	424 Array<double> retval;
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 if (n > 0)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	427 {
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10331 diff changeset	428 retval.clear (n, 1);
7537 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 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	431 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	432 else if (n < 0)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	433 (*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	434
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	435 return retval;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	436 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	437
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	438 // 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	439 // 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	440 // 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	441 // 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	442
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	443 void
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	444 octave_rand::initialize_ranlib_generators (void)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	445 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	446 octave_localtime tm;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	447 int stored_distribution = current_distribution;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	448 F77_FUNC (setcgn, SETCGN) (uniform_dist);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	449
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	450 int hour = tm.hour() + 1;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	451 int minute = tm.min() + 1;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	452 int second = tm.sec() + 1;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	453
10331 6a50aa2cdf00 never compile ranlib with 8-byte integers John W. Eaton <jwe@octave.org> parents: 10314 diff changeset	454 int32_t s0 = tm.mday() * hour * minute * second;
6a50aa2cdf00 never compile ranlib with 8-byte integers John W. Eaton <jwe@octave.org> parents: 10314 diff changeset	455 int32_t s1 = hour * minute * second;
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	456
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	457 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	458 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	459
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	460 F77_FUNC (setall, SETALL) (s0, s1);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	461 F77_FUNC (setcgn, SETCGN) (stored_distribution);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	462 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	463
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	464 void
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	465 octave_rand::initialize_mersenne_twister (void)
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 oct_init_by_entropy ();
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 ColumnVector s = get_internal_state ();
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	470
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	471 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	472
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	473 oct_init_by_entropy ();
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	474 s = get_internal_state ();
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	475 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	476
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	477 oct_init_by_entropy ();
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	478 s = get_internal_state ();
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	479 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	480
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	481 oct_init_by_entropy ();
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	482 s = get_internal_state ();
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	483 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	484
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	485 oct_init_by_entropy ();
d227d096d49e oct-rand.cc (initialize_mersenne_twister): use separate initializations for each generator Brian Gough parents: 7537 diff changeset	486 s = get_internal_state ();
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	487 rand_states[gamma_dist] = s;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	488 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	489
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	490 ColumnVector
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	491 octave_rand::get_internal_state (void)
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 ColumnVector s (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 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	496
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	497 oct_get_state (tmp);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	498
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	499 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	500 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	501
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	502 return s;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	503 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	504
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	505 void
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	506 octave_rand::save_state (void)
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 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	509 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	510
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	511 int
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	512 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	513 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	514 int retval = unknown_dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	515
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	516 if (d == "uniform" \|\| d == "rand")
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	517 retval = uniform_dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	518 else if (d == "normal" \|\| d == "randn")
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	519 retval = normal_dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	520 else if (d == "exponential" \|\| d == "rande")
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	521 retval = expon_dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	522 else if (d == "poisson" \|\| d == "randp")
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	523 retval = poisson_dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	524 else if (d == "gamma" \|\| d == "randg")
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	525 retval = gamma_dist;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	526 else
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	527 (*current_liboctave_error_handler)
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	528 ("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	529
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	530 return retval;
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	531 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	532
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	533 void
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	534 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	535 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	536 octave_idx_type len = s.length ();
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	537 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	538
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	539 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	540
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	541 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	542 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	543
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	544 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	545 oct_set_state (tmp);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	546 else
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	547 oct_init_by_array (tmp, len);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	548 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	549
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	550 void
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	551 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	552 {
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	553 if (dist != current_distribution)
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 current_distribution = 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 set_internal_state (rand_states[dist]);
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	558 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	559 }
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	560
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	561 #define MAKE_RAND(len) \
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	562 do \
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	563 { \
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	564 double val; \
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	565 for (volatile octave_idx_type i = 0; i < len; i++) \
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	566 { \
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	567 octave_quit (); \
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	568 RAND_FUNC (val); \
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	569 v[i] = val; \
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	570 } \
4308 b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	571 } \
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	572 while (0)
b738d1a02adb [project @ 2003-01-24 19:37:12 by jwe] jwe parents: diff changeset	573
7537 a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	574 void
a2950622f070 make octave_rand a proper singleton class John W. Eaton <jwe@octave.org> parents: 7536 diff changeset	575 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	576 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	577 if (len < 1)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	578 return;
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 switch (current_distribution)
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	581 {
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	582 case uniform_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	583 if (use_old_generators)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	584 {
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	585 #define RAND_FUNC(x) F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, x)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	586 MAKE_RAND (len);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	587 #undef RAND_FUNC
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	588 }
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	589 else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	590 oct_fill_randu (len, v);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	591 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	592
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	593 case normal_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	594 if (use_old_generators)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	595 {
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	596 #define RAND_FUNC(x) F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, x)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	597 MAKE_RAND (len);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	598 #undef RAND_FUNC
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	599 }
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	600 else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	601 oct_fill_randn (len, v);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	602 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	603
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	604 case expon_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	605 if (use_old_generators)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	606 {
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	607 #define RAND_FUNC(x) F77_FUNC (dgenexp, DGENEXP) (1.0, x)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	608 MAKE_RAND (len);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	609 #undef RAND_FUNC
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	610 }
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	611 else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	612 oct_fill_rande (len, v);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	613 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	614
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	615 case poisson_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	616 if (use_old_generators)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	617 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	618 if (a < 0.0 \|\| xisnan(a) \|\| xisinf(a))
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	619 #define RAND_FUNC(x) x = octave_NaN;
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	620 MAKE_RAND (len);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	621 #undef RAND_FUNC
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	622 else
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	623 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	624 // workaround bug in ignpoi, by calling with different Mu
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	625 double tmp;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	626 F77_FUNC (dignpoi, DIGNPOI) (a + 1, tmp);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	627 #define RAND_FUNC(x) F77_FUNC (dignpoi, DIGNPOI) (a, x)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	628 MAKE_RAND (len);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	629 #undef RAND_FUNC
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	630 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	631 }
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	632 else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	633 oct_fill_randp (a, len, v);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	634 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	635
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	636 case gamma_dist:
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	637 if (use_old_generators)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	638 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	639 if (a <= 0.0 \|\| xisnan(a) \|\| xisinf(a))
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	640 #define RAND_FUNC(x) x = octave_NaN;
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	641 MAKE_RAND (len);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	642 #undef RAND_FUNC
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	643 else
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	644 #define RAND_FUNC(x) F77_FUNC (dgengam, DGENGAM) (1.0, a, x)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	645 MAKE_RAND (len);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	646 #undef RAND_FUNC
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	647 }
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	648 else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	649 oct_fill_randg (a, len, v);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	650 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	651
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	652 default:
7536 4dda6fbc8ba6 oct-rand.cc: fix typo John W. Eaton <jwe@octave.org> parents: 7535 diff changeset	653 (*current_liboctave_error_handler)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	654 ("rand: invalid distribution ID = %d", current_distribution);
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	655 break;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	656 }
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	657
7533 ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	658 save_state ();
ff52243af934 save state separately for each MT random number generator John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	659
5730 109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	660 return;
109fdf7b3dcb [project @ 2006-04-03 19:18:26 by jwe] jwe parents: 5307 diff changeset	661 }

Mercurial > hg > octave-nkf

annotate liboctave/oct-rand.cc @ 10350:12884915a8e4