Mercurial > octave

comparison liboctave/numeric/oct-rand.cc @ 31607:aac27ad79be6 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Re-indent code after switch to using namespace macros. * build-env.h, build-env.in.cc, Cell.h, __betainc__.cc, __eigs__.cc, __ftp__.cc, __ichol__.cc, __ilu__.cc, __isprimelarge__.cc, __magick_read__.cc, __pchip_deriv__.cc, amd.cc, base-text-renderer.cc, base-text-renderer.h, besselj.cc, bitfcns.cc, bsxfun.cc, c-file-ptr-stream.h, call-stack.cc, call-stack.h, ccolamd.cc, cellfun.cc, chol.cc, colamd.cc, dasrt.cc, data.cc, debug.cc, defaults.cc, defaults.h, det.cc, display.cc, display.h, dlmread.cc, dynamic-ld.cc, dynamic-ld.h, ellipj.cc, environment.cc, environment.h, error.cc, error.h, errwarn.h, event-manager.cc, event-manager.h, event-queue.cc, event-queue.h, fcn-info.cc, fcn-info.h, fft.cc, fft2.cc, file-io.cc, filter.cc, find.cc, ft-text-renderer.cc, ft-text-renderer.h, gcd.cc, gl-render.cc, gl-render.h, gl2ps-print.cc, gl2ps-print.h, graphics-toolkit.cc, graphics-toolkit.h, graphics.cc, gsvd.cc, gtk-manager.cc, gtk-manager.h, help.cc, help.h, hook-fcn.cc, hook-fcn.h, input.cc, input.h, interpreter-private.cc, interpreter-private.h, interpreter.cc, interpreter.h, inv.cc, jsondecode.cc, jsonencode.cc, latex-text-renderer.cc, latex-text-renderer.h, load-path.cc, load-path.h, load-save.cc, load-save.h, lookup.cc, ls-hdf5.cc, ls-mat4.cc, ls-mat5.cc, lsode.cc, lu.cc, mappers.cc, matrix_type.cc, max.cc, mex.cc, mexproto.h, mxarray.h, mxtypes.in.h, oct-errno.in.cc, oct-hdf5-types.cc, oct-hist.cc, oct-hist.h, oct-map.cc, oct-map.h, oct-opengl.h, oct-prcstrm.h, oct-process.cc, oct-process.h, oct-stdstrm.h, oct-stream.cc, oct-stream.h, oct-strstrm.h, octave-default-image.h, ordqz.cc, ordschur.cc, pager.cc, pager.h, pinv.cc, pow2.cc, pr-output.cc, psi.cc, qr.cc, quadcc.cc, rand.cc, regexp.cc, settings.cc, settings.h, sighandlers.cc, sighandlers.h, sparse-xpow.cc, sqrtm.cc, stack-frame.cc, stack-frame.h, stream-euler.cc, strfns.cc, svd.cc, syminfo.cc, syminfo.h, symrcm.cc, symrec.cc, symrec.h, symscope.cc, symscope.h, symtab.cc, symtab.h, sysdep.cc, sysdep.h, text-engine.cc, text-engine.h, text-renderer.cc, text-renderer.h, time.cc, toplev.cc, typecast.cc, url-handle-manager.cc, url-handle-manager.h, urlwrite.cc, utils.cc, utils.h, variables.cc, variables.h, xdiv.cc, __delaunayn__.cc, __init_fltk__.cc, __init_gnuplot__.cc, __ode15__.cc, __voronoi__.cc, audioread.cc, convhulln.cc, gzip.cc, cdef-class.cc, cdef-class.h, cdef-fwd.h, cdef-manager.cc, cdef-manager.h, cdef-method.cc, cdef-method.h, cdef-object.cc, cdef-object.h, cdef-package.cc, cdef-package.h, cdef-property.cc, cdef-property.h, cdef-utils.cc, cdef-utils.h, ov-base-diag.cc, ov-base-int.cc, ov-base-mat.cc, ov-base-mat.h, ov-base-scalar.cc, ov-base.cc, ov-base.h, ov-bool-mat.cc, ov-bool-mat.h, ov-bool-sparse.cc, ov-bool.cc, ov-builtin.h, ov-cell.cc, ov-ch-mat.cc, ov-class.cc, ov-class.h, ov-classdef.cc, ov-classdef.h, ov-complex.cc, ov-cx-diag.cc, ov-cx-mat.cc, ov-cx-sparse.cc, ov-dld-fcn.cc, ov-dld-fcn.h, ov-fcn-handle.cc, ov-fcn-handle.h, ov-fcn.h, ov-float.cc, ov-flt-complex.cc, ov-flt-cx-diag.cc, ov-flt-cx-mat.cc, ov-flt-re-diag.cc, ov-flt-re-mat.cc, ov-flt-re-mat.h, ov-intx.h, ov-java.cc, ov-lazy-idx.cc, ov-legacy-range.cc, ov-magic-int.cc, ov-mex-fcn.cc, ov-mex-fcn.h, ov-null-mat.cc, ov-perm.cc, ov-range.cc, ov-re-diag.cc, ov-re-mat.cc, ov-re-mat.h, ov-re-sparse.cc, ov-scalar.cc, ov-str-mat.cc, ov-struct.cc, ov-typeinfo.cc, ov-typeinfo.h, ov-usr-fcn.cc, ov-usr-fcn.h, ov.cc, ov.h, ovl.h, octave.cc, octave.h, op-b-sbm.cc, op-bm-sbm.cc, op-cs-scm.cc, op-fm-fcm.cc, op-fs-fcm.cc, op-s-scm.cc, op-scm-cs.cc, op-scm-s.cc, op-sm-cs.cc, ops.h, anon-fcn-validator.cc, anon-fcn-validator.h, bp-table.cc, bp-table.h, comment-list.cc, comment-list.h, filepos.h, lex.h, oct-lvalue.cc, oct-lvalue.h, parse.h, profiler.cc, profiler.h, pt-anon-scopes.cc, pt-anon-scopes.h, pt-arg-list.cc, pt-arg-list.h, pt-args-block.cc, pt-args-block.h, pt-array-list.cc, pt-array-list.h, pt-assign.cc, pt-assign.h, pt-binop.cc, pt-binop.h, pt-bp.cc, pt-bp.h, pt-cbinop.cc, pt-cbinop.h, pt-cell.cc, pt-cell.h, pt-check.cc, pt-check.h, pt-classdef.cc, pt-classdef.h, pt-cmd.h, pt-colon.cc, pt-colon.h, pt-const.cc, pt-const.h, pt-decl.cc, pt-decl.h, pt-eval.cc, pt-eval.h, pt-except.cc, pt-except.h, pt-exp.cc, pt-exp.h, pt-fcn-handle.cc, pt-fcn-handle.h, pt-id.cc, pt-id.h, pt-idx.cc, pt-idx.h, pt-jump.h, pt-loop.cc, pt-loop.h, pt-mat.cc, pt-mat.h, pt-misc.cc, pt-misc.h, pt-pr-code.cc, pt-pr-code.h, pt-select.cc, pt-select.h, pt-spmd.cc, pt-spmd.h, pt-stmt.cc, pt-stmt.h, pt-tm-const.cc, pt-tm-const.h, pt-unop.cc, pt-unop.h, pt-walk.cc, pt-walk.h, pt.cc, pt.h, token.cc, token.h, Range.cc, Range.h, idx-vector.cc, idx-vector.h, range-fwd.h, CollocWt.cc, CollocWt.h, aepbalance.cc, aepbalance.h, chol.cc, chol.h, gepbalance.cc, gepbalance.h, gsvd.cc, gsvd.h, hess.cc, hess.h, lo-mappers.cc, lo-mappers.h, lo-specfun.cc, lo-specfun.h, lu.cc, lu.h, oct-convn.cc, oct-convn.h, oct-fftw.cc, oct-fftw.h, oct-norm.cc, oct-norm.h, oct-rand.cc, oct-rand.h, oct-spparms.cc, oct-spparms.h, qr.cc, qr.h, qrp.cc, qrp.h, randgamma.cc, randgamma.h, randmtzig.cc, randmtzig.h, randpoisson.cc, randpoisson.h, schur.cc, schur.h, sparse-chol.cc, sparse-chol.h, sparse-lu.cc, sparse-lu.h, sparse-qr.cc, sparse-qr.h, svd.cc, svd.h, child-list.cc, child-list.h, dir-ops.cc, dir-ops.h, file-ops.cc, file-ops.h, file-stat.cc, file-stat.h, lo-sysdep.cc, lo-sysdep.h, lo-sysinfo.cc, lo-sysinfo.h, mach-info.cc, mach-info.h, oct-env.cc, oct-env.h, oct-group.cc, oct-group.h, oct-password.cc, oct-password.h, oct-syscalls.cc, oct-syscalls.h, oct-time.cc, oct-time.h, oct-uname.cc, oct-uname.h, action-container.cc, action-container.h, base-list.h, cmd-edit.cc, cmd-edit.h, cmd-hist.cc, cmd-hist.h, f77-fcn.h, file-info.cc, file-info.h, lo-array-errwarn.cc, lo-array-errwarn.h, lo-hash.cc, lo-hash.h, lo-ieee.h, lo-regexp.cc, lo-regexp.h, lo-utils.cc, lo-utils.h, oct-base64.cc, oct-base64.h, oct-glob.cc, oct-glob.h, oct-inttypes.h, oct-mutex.cc, oct-mutex.h, oct-refcount.h, oct-shlib.cc, oct-shlib.h, oct-sparse.cc, oct-sparse.h, oct-string.h, octave-preserve-stream-state.h, pathsearch.cc, pathsearch.h, quit.cc, quit.h, unwind-prot.cc, unwind-prot.h, url-transfer.cc, url-transfer.h: Re-indent code after switch to using namespace macros.
author Rik <rik@octave.org>
date Thu, 01 Dec 2022 18:02:15 -0800
parents e88a07dec498
children 597f3ee61a48
comparison
equal deleted inserted replaced
31605:e88a07dec498 31607:aac27ad79be6
46 #include "randpoisson.h" 46 #include "randpoisson.h"
47 #include "singleton-cleanup.h" 47 #include "singleton-cleanup.h"
48 48
49 OCTAVE_BEGIN_NAMESPACE(octave) 49 OCTAVE_BEGIN_NAMESPACE(octave)
50 50
51 rand *rand::m_instance = nullptr; 51 rand *rand::m_instance = nullptr;
52 52
53 rand::rand (void) 53 rand::rand (void)
54 : m_current_distribution (uniform_dist), m_use_old_generators (false), 54 : m_current_distribution (uniform_dist), m_use_old_generators (false),
55 m_rand_states () 55 m_rand_states ()
56 { 56 {
57 initialize_ranlib_generators (); 57 initialize_ranlib_generators ();
58 58
59 initialize_mersenne_twister (); 59 initialize_mersenne_twister ();
60 } 60 }
61 61
62 bool rand::instance_ok (void) 62 bool rand::instance_ok (void)
63 { 63 {
64 bool retval = true; 64 bool retval = true;
65 65
66 if (! m_instance) 66 if (! m_instance)
67 { 67 {
68 m_instance = new rand (); 68 m_instance = new rand ();
69 singleton_cleanup_list::add (cleanup_instance); 69 singleton_cleanup_list::add (cleanup_instance);
70 } 70 }
71 71
72 return retval; 72 return retval;
73 } 73 }
74 74
75 double rand::do_seed (void) 75 double rand::do_seed (void)
76 { 76 {
77 union d2i { double d; int32_t i[2]; }; 77 union d2i { double d; int32_t i[2]; };
78 union d2i u; 78 union d2i u;
79 79
80 mach_info::float_format ff = mach_info::native_float_format (); 80 mach_info::float_format ff = mach_info::native_float_format ();
81 81
82 switch (ff) 82 switch (ff)
83 { 83 {
84 case mach_info::flt_fmt_ieee_big_endian: 84 case mach_info::flt_fmt_ieee_big_endian:
85 F77_FUNC (getsd, GETSD) (u.i[1], u.i[0]); 85 F77_FUNC (getsd, GETSD) (u.i[1], u.i[0]);
86 break; 86 break;
87 87
88 default: 88 default:
89 F77_FUNC (getsd, GETSD) (u.i[0], u.i[1]); 89 F77_FUNC (getsd, GETSD) (u.i[0], u.i[1]);
90 break; 90 break;
91 } 91 }
92 92
93 return u.d; 93 return u.d;
94 } 94 }
95 95
96 static int32_t 96 static int32_t
97 force_to_fit_range (int32_t i, int32_t lo, int32_t hi) 97 force_to_fit_range (int32_t i, int32_t lo, int32_t hi)
98 { 98 {
99 assert (hi > lo && lo >= 0); 99 assert (hi > lo && lo >= 0);
100 100
101 i = (i > 0 ? i : -i); 101 i = (i > 0 ? i : -i);
102 102
103 if (i < lo) 103 if (i < lo)
104 i = lo; 104 i = lo;
105 else if (i > hi) 105 else if (i > hi)
106 i = i % hi; 106 i = i % hi;
107 107
108 return i; 108 return i;
109 } 109 }
110 110
111 void rand::do_seed (double s) 111 void rand::do_seed (double s)
112 { 112 {
113 m_use_old_generators = true; 113 m_use_old_generators = true;
114 114
115 int i0, i1; 115 int i0, i1;
116 union d2i { double d; int32_t i[2]; }; 116 union d2i { double d; int32_t i[2]; };
117 union d2i u; 117 union d2i u;
118 u.d = s; 118 u.d = s;
119 119
120 mach_info::float_format ff = mach_info::native_float_format (); 120 mach_info::float_format ff = mach_info::native_float_format ();
121 121
122 switch (ff) 122 switch (ff)
123 { 123 {
124 case mach_info::flt_fmt_ieee_big_endian: 124 case mach_info::flt_fmt_ieee_big_endian:
125 i1 = force_to_fit_range (u.i[0], 1, 2147483563); 125 i1 = force_to_fit_range (u.i[0], 1, 2147483563);
126 i0 = force_to_fit_range (u.i[1], 1, 2147483399); 126 i0 = force_to_fit_range (u.i[1], 1, 2147483399);
127 break; 127 break;
128 128
129 default: 129 default:
130 i0 = force_to_fit_range (u.i[0], 1, 2147483563); 130 i0 = force_to_fit_range (u.i[0], 1, 2147483563);
131 i1 = force_to_fit_range (u.i[1], 1, 2147483399); 131 i1 = force_to_fit_range (u.i[1], 1, 2147483399);
132 break; 132 break;
133 } 133 }
134 134
135 F77_FUNC (setsd, SETSD) (i0, i1); 135 F77_FUNC (setsd, SETSD) (i0, i1);
136 } 136 }
137 137
138 void rand::do_reset (void) 138 void rand::do_reset (void)
139 { 139 {
140 m_use_old_generators = true; 140 m_use_old_generators = true;
141 initialize_ranlib_generators (); 141 initialize_ranlib_generators ();
142 } 142 }
143 143
144 uint32NDArray rand::do_state (const std::string& d) 144 uint32NDArray rand::do_state (const std::string& d)
145 { 145 {
146 return m_rand_states[d.empty () ? m_current_distribution : get_dist_id (d)]; 146 return m_rand_states[d.empty () ? m_current_distribution : get_dist_id (d)];
147 } 147 }
148 148
149 void rand::do_state (const uint32NDArray& s, const std::string& d) 149 void rand::do_state (const uint32NDArray& s, const std::string& d)
150 { 150 {
151 m_use_old_generators = false; 151 m_use_old_generators = false;
152 152
153 int old_dist = m_current_distribution; 153 int old_dist = m_current_distribution;
154 154
155 int new_dist = (d.empty () ? m_current_distribution : get_dist_id (d)); 155 int new_dist = (d.empty () ? m_current_distribution : get_dist_id (d));
156 156
157 uint32NDArray saved_state; 157 uint32NDArray saved_state;
158 158
159 if (old_dist != new_dist) 159 if (old_dist != new_dist)
160 saved_state = get_internal_state (); 160 saved_state = get_internal_state ();
161 161
162 set_internal_state (s); 162 set_internal_state (s);
163 163
164 m_rand_states[new_dist] = get_internal_state (); 164 m_rand_states[new_dist] = get_internal_state ();
165 165
166 if (old_dist != new_dist) 166 if (old_dist != new_dist)
167 m_rand_states[old_dist] = saved_state; 167 m_rand_states[old_dist] = saved_state;
168 } 168 }
169 169
170 void rand::do_reset (const std::string& d) 170 void rand::do_reset (const std::string& d)
171 { 171 {
172 m_use_old_generators = false; 172 m_use_old_generators = false;
173 173
174 int old_dist = m_current_distribution; 174 int old_dist = m_current_distribution;
175 175
176 int new_dist = (d.empty () ? m_current_distribution : get_dist_id (d)); 176 int new_dist = (d.empty () ? m_current_distribution : get_dist_id (d));
177 177
178 uint32NDArray saved_state; 178 uint32NDArray saved_state;
179 179
180 if (old_dist != new_dist) 180 if (old_dist != new_dist)
181 saved_state = get_internal_state (); 181 saved_state = get_internal_state ();
182 182
183 init_mersenne_twister (); 183 init_mersenne_twister ();
184 m_rand_states[new_dist] = get_internal_state (); 184 m_rand_states[new_dist] = get_internal_state ();
185 185
186 if (old_dist != new_dist) 186 if (old_dist != new_dist)
187 m_rand_states[old_dist] = saved_state; 187 m_rand_states[old_dist] = saved_state;
188 } 188 }
189 189
190 std::string rand::do_distribution (void) 190 std::string rand::do_distribution (void)
191 { 191 {
192 std::string retval; 192 std::string retval;
193 193
194 switch (m_current_distribution) 194 switch (m_current_distribution)
195 { 195 {
196 case uniform_dist: 196 case uniform_dist:
197 retval = "uniform"; 197 retval = "uniform";
198 break; 198 break;
199 199
200 case normal_dist: 200 case normal_dist:
201 retval = "normal"; 201 retval = "normal";
202 break; 202 break;
203 203
204 case expon_dist: 204 case expon_dist:
205 retval = "exponential"; 205 retval = "exponential";
206 break; 206 break;
207 207
208 case poisson_dist: 208 case poisson_dist:
209 retval = "poisson"; 209 retval = "poisson";
210 break; 210 break;
211 211
212 case gamma_dist: 212 case gamma_dist:
213 retval = "gamma"; 213 retval = "gamma";
214 break; 214 break;
215 215
216 default: 216 default:
217 (*current_liboctave_error_handler) 217 (*current_liboctave_error_handler)
218 ("rand: invalid distribution ID = %d", m_current_distribution); 218 ("rand: invalid distribution ID = %d", m_current_distribution);
219 break; 219 break;
220 } 220 }
221 221
222 return retval; 222 return retval;
223 } 223 }
224 224
225 void rand::do_distribution (const std::string& d) 225 void rand::do_distribution (const std::string& d)
226 { 226 {
227 int id = get_dist_id (d); 227 int id = get_dist_id (d);
228 228
229 switch (id) 229 switch (id)
230 { 230 {
231 case uniform_dist: 231 case uniform_dist:
232 rand::uniform_distribution (); 232 rand::uniform_distribution ();
233 break; 233 break;
234 234
235 case normal_dist: 235 case normal_dist:
236 rand::normal_distribution (); 236 rand::normal_distribution ();
237 break; 237 break;
238 238
239 case expon_dist: 239 case expon_dist:
240 rand::exponential_distribution (); 240 rand::exponential_distribution ();
241 break; 241 break;
242 242
243 case poisson_dist: 243 case poisson_dist:
244 rand::poisson_distribution (); 244 rand::poisson_distribution ();
245 break; 245 break;
246 246
247 case gamma_dist: 247 case gamma_dist:
248 rand::gamma_distribution (); 248 rand::gamma_distribution ();
249 break; 249 break;
250 250
251 default: 251 default:
252 (*current_liboctave_error_handler) 252 (*current_liboctave_error_handler)
253 ("rand: invalid distribution ID = %d", id); 253 ("rand: invalid distribution ID = %d", id);
254 break; 254 break;
255 } 255 }
256 } 256 }
257 257
258 void rand::do_uniform_distribution (void) 258 void rand::do_uniform_distribution (void)
259 { 259 {
260 switch_to_generator (uniform_dist); 260 switch_to_generator (uniform_dist);
261 261
262 F77_FUNC (setcgn, SETCGN) (uniform_dist); 262 F77_FUNC (setcgn, SETCGN) (uniform_dist);
263 } 263 }
264 264
265 void rand::do_normal_distribution (void) 265 void rand::do_normal_distribution (void)
266 { 266 {
267 switch_to_generator (normal_dist); 267 switch_to_generator (normal_dist);
268 268
269 F77_FUNC (setcgn, SETCGN) (normal_dist); 269 F77_FUNC (setcgn, SETCGN) (normal_dist);
270 } 270 }
271 271
272 void rand::do_exponential_distribution (void) 272 void rand::do_exponential_distribution (void)
273 { 273 {
274 switch_to_generator (expon_dist); 274 switch_to_generator (expon_dist);
275 275
276 F77_FUNC (setcgn, SETCGN) (expon_dist); 276 F77_FUNC (setcgn, SETCGN) (expon_dist);
277 } 277 }
278 278
279 void rand::do_poisson_distribution (void) 279 void rand::do_poisson_distribution (void)
280 { 280 {
281 switch_to_generator (poisson_dist); 281 switch_to_generator (poisson_dist);
282 282
283 F77_FUNC (setcgn, SETCGN) (poisson_dist); 283 F77_FUNC (setcgn, SETCGN) (poisson_dist);
284 } 284 }
285 285
286 void rand::do_gamma_distribution (void) 286 void rand::do_gamma_distribution (void)
287 { 287 {
288 switch_to_generator (gamma_dist); 288 switch_to_generator (gamma_dist);
289 289
290 F77_FUNC (setcgn, SETCGN) (gamma_dist); 290 F77_FUNC (setcgn, SETCGN) (gamma_dist);
291 } 291 }
292 292
293 template <> 293 template <>
294 OCTAVE_API double rand::uniform<double> (void) 294 OCTAVE_API double rand::uniform<double> (void)
295 { 295 {
296 double retval; 296 double retval;
297 297
298 if (m_use_old_generators) 298 if (m_use_old_generators)
299 F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval); 299 F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval);
300 else 300 else
301 retval = rand_uniform<double> (); 301 retval = rand_uniform<double> ();
302 302
303 return retval; 303 return retval;
304 } 304 }
305 305
306 template <> 306 template <>
307 OCTAVE_API double rand::normal<double> (void) 307 OCTAVE_API double rand::normal<double> (void)
308 { 308 {
309 double retval; 309 double retval;
310 310
311 if (m_use_old_generators) 311 if (m_use_old_generators)
312 F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval); 312 F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval);
313 else 313 else
314 retval = rand_normal<double> (); 314 retval = rand_normal<double> ();
315 315
316 return retval; 316 return retval;
317 } 317 }
318 318
319 template <> 319 template <>
320 OCTAVE_API double rand::exponential<double> (void) 320 OCTAVE_API double rand::exponential<double> (void)
321 { 321 {
322 double retval; 322 double retval;
323 323
324 if (m_use_old_generators) 324 if (m_use_old_generators)
325 F77_FUNC (dgenexp, DGENEXP) (1.0, retval); 325 F77_FUNC (dgenexp, DGENEXP) (1.0, retval);
326 else 326 else
327 retval = rand_exponential<double> (); 327 retval = rand_exponential<double> ();
328 328
329 return retval; 329 return retval;
330 } 330 }
331 331
332 template <> 332 template <>
333 OCTAVE_API double rand::poisson<double> (double a) 333 OCTAVE_API double rand::poisson<double> (double a)
334 { 334 {
335 double retval; 335 double retval;
336 336
337 if (m_use_old_generators) 337 if (m_use_old_generators)
338 { 338 {
339 if (a < 0.0 || ! math::isfinite (a)) 339 if (a < 0.0 || ! math::isfinite (a))
340 retval = numeric_limits<double>::NaN (); 340 retval = numeric_limits<double>::NaN ();
341 else 341 else
342 { 342 {
343 // workaround bug in ignpoi, by calling with different Mu 343 // workaround bug in ignpoi, by calling with different Mu
344 F77_FUNC (dignpoi, DIGNPOI) (a + 1, retval); 344 F77_FUNC (dignpoi, DIGNPOI) (a + 1, retval);
345 F77_FUNC (dignpoi, DIGNPOI) (a, retval); 345 F77_FUNC (dignpoi, DIGNPOI) (a, retval);
346 } 346 }
347 } 347 }
348 else 348 else
349 retval = rand_poisson<double> (a);
350
351 return retval;
352 }
353
354 template <>
355 OCTAVE_API double rand::gamma<double> (double a)
356 {
357 double retval;
358
359 if (m_use_old_generators)
360 {
361 if (a <= 0.0 || ! math::isfinite (a))
362 retval = numeric_limits<double>::NaN ();
363 else
364 F77_FUNC (dgengam, DGENGAM) (1.0, a, retval);
365 }
366 else
367 retval = rand_gamma<double> (a);
368
369 return retval;
370 }
371
372 template <>
373 OCTAVE_API float rand::uniform<float> (void)
374 {
375 float retval;
376
377 if (m_use_old_generators)
378 F77_FUNC (fgenunf, FGENUNF) (0.0f, 1.0f, retval);
379 else
380 retval = rand_uniform<float> ();
381
382 return retval;
383 }
384
385 template <>
386 OCTAVE_API float rand::normal<float> (void)
387 {
388 float retval;
389
390 if (m_use_old_generators)
391 F77_FUNC (fgennor, FGENNOR) (0.0f, 1.0f, retval);
392 else
393 retval = rand_normal<float> ();
394
395 return retval;
396 }
397
398 template <>
399 OCTAVE_API float rand::exponential<float> (void)
400 {
401 float retval;
402
403 if (m_use_old_generators)
404 F77_FUNC (fgenexp, FGENEXP) (1.0f, retval);
405 else
406 retval = rand_exponential<float> ();
407
408 return retval;
409 }
410
411 template <>
412 OCTAVE_API float rand::poisson<float> (float a)
413 {
414 float retval;
415
416 if (m_use_old_generators)
417 {
418 if (a < 0.0f || ! math::isfinite (a))
419 retval = numeric_limits<float>::NaN ();
420 else
421 {
422 // workaround bug in ignpoi, by calling with different Mu
423 F77_FUNC (fignpoi, FIGNPOI) (a + 1, retval);
424 F77_FUNC (fignpoi, FIGNPOI) (a, retval);
425 }
426 }
427 else
428 {
429 // Keep poisson distribution in double precision for accuracy
349 retval = rand_poisson<double> (a); 430 retval = rand_poisson<double> (a);
350 431 }
351 return retval; 432
352 } 433 return retval;
353 434 }
354 template <> 435
355 OCTAVE_API double rand::gamma<double> (double a) 436 template <>
356 { 437 OCTAVE_API float rand::gamma<float> (float a)
357 double retval; 438 {
358 439 float retval;
359 if (m_use_old_generators) 440
360 { 441 if (m_use_old_generators)
361 if (a <= 0.0 || ! math::isfinite (a)) 442 {
362 retval = numeric_limits<double>::NaN (); 443 if (a <= 0.0f || ! math::isfinite (a))
363 else 444 retval = numeric_limits<float>::NaN ();
364 F77_FUNC (dgengam, DGENGAM) (1.0, a, retval); 445 else
365 } 446 F77_FUNC (fgengam, FGENGAM) (1.0f, a, retval);
366 else 447 }
367 retval = rand_gamma<double> (a); 448 else
368 449 retval = rand_gamma<float> (a);
369 return retval; 450
370 } 451 return retval;
371 452 }
372 template <> 453
373 OCTAVE_API float rand::uniform<float> (void) 454 template <typename T>
374 { 455 T rand::do_scalar (T a)
375 float retval; 456 {
376 457 T retval = 0;
377 if (m_use_old_generators) 458
378 F77_FUNC (fgenunf, FGENUNF) (0.0f, 1.0f, retval); 459 switch (m_current_distribution)
379 else 460 {
380 retval = rand_uniform<float> (); 461 case uniform_dist:
381 462 retval = uniform<T> ();
382 return retval; 463 break;
383 } 464
384 465 case normal_dist:
385 template <> 466 retval = normal<T> ();
386 OCTAVE_API float rand::normal<float> (void) 467 break;
387 { 468
388 float retval; 469 case expon_dist:
389 470 retval = exponential<T> ();
390 if (m_use_old_generators) 471 break;
391 F77_FUNC (fgennor, FGENNOR) (0.0f, 1.0f, retval); 472
392 else 473 case poisson_dist:
393 retval = rand_normal<float> (); 474 retval = poisson<T> (a);
394 475 break;
395 return retval; 476
396 } 477 case gamma_dist:
397 478 retval = gamma<T> (a);
398 template <> 479 break;
399 OCTAVE_API float rand::exponential<float> (void) 480
400 { 481 default:
401 float retval;
402
403 if (m_use_old_generators)
404 F77_FUNC (fgenexp, FGENEXP) (1.0f, retval);
405 else
406 retval = rand_exponential<float> ();
407
408 return retval;
409 }
410
411 template <>
412 OCTAVE_API float rand::poisson<float> (float a)
413 {
414 float retval;
415
416 if (m_use_old_generators)
417 {
418 if (a < 0.0f || ! math::isfinite (a))
419 retval = numeric_limits<float>::NaN ();
420 else
421 {
422 // workaround bug in ignpoi, by calling with different Mu
423 F77_FUNC (fignpoi, FIGNPOI) (a + 1, retval);
424 F77_FUNC (fignpoi, FIGNPOI) (a, retval);
425 }
426 }
427 else
428 {
429 // Keep poisson distribution in double precision for accuracy
430 retval = rand_poisson<double> (a);
431 }
432
433 return retval;
434 }
435
436 template <>
437 OCTAVE_API float rand::gamma<float> (float a)
438 {
439 float retval;
440
441 if (m_use_old_generators)
442 {
443 if (a <= 0.0f || ! math::isfinite (a))
444 retval = numeric_limits<float>::NaN ();
445 else
446 F77_FUNC (fgengam, FGENGAM) (1.0f, a, retval);
447 }
448 else
449 retval = rand_gamma<float> (a);
450
451 return retval;
452 }
453
454 template <typename T>
455 T rand::do_scalar (T a)
456 {
457 T retval = 0;
458
459 switch (m_current_distribution)
460 {
461 case uniform_dist:
462 retval = uniform<T> ();
463 break;
464
465 case normal_dist:
466 retval = normal<T> ();
467 break;
468
469 case expon_dist:
470 retval = exponential<T> ();
471 break;
472
473 case poisson_dist:
474 retval = poisson<T> (a);
475 break;
476
477 case gamma_dist:
478 retval = gamma<T> (a);
479 break;
480
481 default:
482 (*current_liboctave_error_handler)
483 ("rand: invalid distribution ID = %d", m_current_distribution);
484 break;
485 }
486
487 if (! m_use_old_generators)
488 save_state ();
489
490 return retval;
491 }
492
493 template OCTAVE_API double rand::do_scalar<double> (double);
494 template OCTAVE_API float rand::do_scalar<float> (float);
495
496 template <typename T>
497 Array<T>
498 rand::do_vector (octave_idx_type n, T a)
499 {
500 Array<T> retval;
501
502 if (n > 0)
503 {
504 retval.clear (n, 1);
505
506 fill (retval.numel (), retval.fortran_vec (), a);
507 }
508 else if (n < 0)
509 (*current_liboctave_error_handler) ("rand: invalid negative argument");
510
511 return retval;
512 }
513
514 template OCTAVE_API Array<double>
515 rand::do_vector<double> (octave_idx_type, double);
516 template OCTAVE_API Array<float>
517 rand::do_vector<float> (octave_idx_type, float);
518
519 NDArray rand::do_nd_array (const dim_vector& dims, double a)
520 {
521 NDArray retval;
522
523 if (! dims.all_zero ())
524 {
525 retval.clear (dims);
526
527 fill (retval.numel (), retval.fortran_vec (), a);
528 }
529
530 return retval;
531 }
532
533 FloatNDArray rand::do_float_nd_array (const dim_vector& dims, float a)
534 {
535 FloatNDArray retval;
536
537 if (! dims.all_zero ())
538 {
539 retval.clear (dims);
540
541 fill (retval.numel (), retval.fortran_vec (), a);
542 }
543
544 return retval;
545 }
546
547 // Make the random number generator give us a different sequence every
548 // time we start octave unless we specifically set the seed. The
549 // technique used below will cycle monthly, but it does seem to
550 // work ok to give fairly different seeds each time Octave starts.
551
552 void rand::initialize_ranlib_generators (void)
553 {
554 sys::localtime tm;
555 int stored_distribution = m_current_distribution;
556 F77_FUNC (setcgn, SETCGN) (uniform_dist);
557
558 int hour = tm.hour () + 1;
559 int minute = tm.min () + 1;
560 int second = tm.sec () + 1;
561
562 int32_t s0 = tm.mday () * hour * minute * second;
563 int32_t s1 = hour * minute * second;
564
565 s0 = force_to_fit_range (s0, 1, 2147483563);
566 s1 = force_to_fit_range (s1, 1, 2147483399);
567
568 F77_FUNC (setall, SETALL) (s0, s1);
569 F77_FUNC (setcgn, SETCGN) (stored_distribution);
570 }
571
572 void rand::initialize_mersenne_twister (void)
573 {
574 uint32NDArray s;
575
576 init_mersenne_twister ();
577 s = get_internal_state ();
578 m_rand_states[uniform_dist] = s;
579
580 init_mersenne_twister ();
581 s = get_internal_state ();
582 m_rand_states[normal_dist] = s;
583
584 init_mersenne_twister ();
585 s = get_internal_state ();
586 m_rand_states[expon_dist] = s;
587
588 init_mersenne_twister ();
589 s = get_internal_state ();
590 m_rand_states[poisson_dist] = s;
591
592 init_mersenne_twister ();
593 s = get_internal_state ();
594 m_rand_states[gamma_dist] = s;
595
596 // All of the initializations above have messed with the internal state.
597 // Restore the state of the currently selected distribution.
598 set_internal_state (m_rand_states[m_current_distribution]);
599 }
600
601 uint32NDArray rand::get_internal_state (void)
602 {
603 uint32NDArray s (dim_vector (MT_N + 1, 1));
604
605 get_mersenne_twister_state (reinterpret_cast<uint32_t *> (s.fortran_vec ()));
606
607 return s;
608 }
609
610 void rand::save_state (void)
611 {
612 m_rand_states[m_current_distribution] = get_internal_state ();;
613 }
614
615 int rand::get_dist_id (const std::string& d)
616 {
617 int retval = unknown_dist;
618
619 if (d == "uniform" || d == "rand")
620 retval = uniform_dist;
621 else if (d == "normal" || d == "randn")
622 retval = normal_dist;
623 else if (d == "exponential" || d == "rande")
624 retval = expon_dist;
625 else if (d == "poisson" || d == "randp")
626 retval = poisson_dist;
627 else if (d == "gamma" || d == "randg")
628 retval = gamma_dist;
629 else
630 (*current_liboctave_error_handler) 482 (*current_liboctave_error_handler)
631 ("rand: invalid distribution '%s'", d.c_str ()); 483 ("rand: invalid distribution ID = %d", m_current_distribution);
632 484 break;
633 return retval; 485 }
634 } 486
635 487 if (! m_use_old_generators)
636 void rand::set_internal_state (const uint32NDArray& s)
637 {
638 octave_idx_type len = s.numel ();
639
640 const uint32_t *sdata = reinterpret_cast <const uint32_t *> (s.data ());
641
642 if (len == MT_N + 1 && sdata[MT_N] <= MT_N && sdata[MT_N] > 0)
643 set_mersenne_twister_state (sdata);
644 else
645 init_mersenne_twister (sdata, len);
646 }
647
648 void rand::switch_to_generator (int dist)
649 {
650 if (dist != m_current_distribution)
651 {
652 m_current_distribution = dist;
653
654 set_internal_state (m_rand_states[dist]);
655 }
656 }
657
658 void rand::fill (octave_idx_type len, double *v, double a)
659 {
660 if (len < 1)
661 return;
662
663 switch (m_current_distribution)
664 {
665 case uniform_dist:
666 if (m_use_old_generators)
667 std::generate_n (v, len, [](void) { double x; F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, x); return x; });
668 else
669 rand_uniform<double> (len, v);
670 break;
671
672 case normal_dist:
673 if (m_use_old_generators)
674 std::generate_n (v, len, [](void) { double x; F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, x); return x; });
675 else
676 rand_normal<double> (len, v);
677 break;
678
679 case expon_dist:
680 if (m_use_old_generators)
681 std::generate_n (v, len, [](void) { double x; F77_FUNC (dgenexp, DGENEXP) (1.0, x); return x; });
682 else
683 rand_exponential<double> (len, v);
684 break;
685
686 case poisson_dist:
687 if (m_use_old_generators)
688 {
689 if (a < 0.0 || ! math::isfinite (a))
690 std::fill_n (v, len, numeric_limits<double>::NaN ());
691 else
692 {
693 // workaround bug in ignpoi, by calling with different Mu
694 double tmp;
695 F77_FUNC (dignpoi, DIGNPOI) (a + 1, tmp);
696 std::generate_n (v, len, [a](void) { double x; F77_FUNC (dignpoi, DIGNPOI) (a, x); return x; });
697 }
698 }
699 else
700 rand_poisson<double> (a, len, v);
701 break;
702
703 case gamma_dist:
704 if (m_use_old_generators)
705 {
706 if (a <= 0.0 || ! math::isfinite (a))
707 std::fill_n (v, len, numeric_limits<double>::NaN ());
708 else
709 std::generate_n (v, len, [a](void) { double x; F77_FUNC (dgengam, DGENGAM) (1.0, a, x); return x; });
710 }
711 else
712 rand_gamma<double> (a, len, v);
713 break;
714
715 default:
716 (*current_liboctave_error_handler)
717 ("rand: invalid distribution ID = %d", m_current_distribution);
718 break;
719 }
720
721 save_state (); 488 save_state ();
722 489
490 return retval;
491 }
492
493 template OCTAVE_API double rand::do_scalar<double> (double);
494 template OCTAVE_API float rand::do_scalar<float> (float);
495
496 template <typename T>
497 Array<T>
498 rand::do_vector (octave_idx_type n, T a)
499 {
500 Array<T> retval;
501
502 if (n > 0)
503 {
504 retval.clear (n, 1);
505
506 fill (retval.numel (), retval.fortran_vec (), a);
507 }
508 else if (n < 0)
509 (*current_liboctave_error_handler) ("rand: invalid negative argument");
510
511 return retval;
512 }
513
514 template OCTAVE_API Array<double>
515 rand::do_vector<double> (octave_idx_type, double);
516 template OCTAVE_API Array<float>
517 rand::do_vector<float> (octave_idx_type, float);
518
519 NDArray rand::do_nd_array (const dim_vector& dims, double a)
520 {
521 NDArray retval;
522
523 if (! dims.all_zero ())
524 {
525 retval.clear (dims);
526
527 fill (retval.numel (), retval.fortran_vec (), a);
528 }
529
530 return retval;
531 }
532
533 FloatNDArray rand::do_float_nd_array (const dim_vector& dims, float a)
534 {
535 FloatNDArray retval;
536
537 if (! dims.all_zero ())
538 {
539 retval.clear (dims);
540
541 fill (retval.numel (), retval.fortran_vec (), a);
542 }
543
544 return retval;
545 }
546
547 // Make the random number generator give us a different sequence every
548 // time we start octave unless we specifically set the seed. The
549 // technique used below will cycle monthly, but it does seem to
550 // work ok to give fairly different seeds each time Octave starts.
551
552 void rand::initialize_ranlib_generators (void)
553 {
554 sys::localtime tm;
555 int stored_distribution = m_current_distribution;
556 F77_FUNC (setcgn, SETCGN) (uniform_dist);
557
558 int hour = tm.hour () + 1;
559 int minute = tm.min () + 1;
560 int second = tm.sec () + 1;
561
562 int32_t s0 = tm.mday () * hour * minute * second;
563 int32_t s1 = hour * minute * second;
564
565 s0 = force_to_fit_range (s0, 1, 2147483563);
566 s1 = force_to_fit_range (s1, 1, 2147483399);
567
568 F77_FUNC (setall, SETALL) (s0, s1);
569 F77_FUNC (setcgn, SETCGN) (stored_distribution);
570 }
571
572 void rand::initialize_mersenne_twister (void)
573 {
574 uint32NDArray s;
575
576 init_mersenne_twister ();
577 s = get_internal_state ();
578 m_rand_states[uniform_dist] = s;
579
580 init_mersenne_twister ();
581 s = get_internal_state ();
582 m_rand_states[normal_dist] = s;
583
584 init_mersenne_twister ();
585 s = get_internal_state ();
586 m_rand_states[expon_dist] = s;
587
588 init_mersenne_twister ();
589 s = get_internal_state ();
590 m_rand_states[poisson_dist] = s;
591
592 init_mersenne_twister ();
593 s = get_internal_state ();
594 m_rand_states[gamma_dist] = s;
595
596 // All of the initializations above have messed with the internal state.
597 // Restore the state of the currently selected distribution.
598 set_internal_state (m_rand_states[m_current_distribution]);
599 }
600
601 uint32NDArray rand::get_internal_state (void)
602 {
603 uint32NDArray s (dim_vector (MT_N + 1, 1));
604
605 get_mersenne_twister_state (reinterpret_cast<uint32_t *> (s.fortran_vec ()));
606
607 return s;
608 }
609
610 void rand::save_state (void)
611 {
612 m_rand_states[m_current_distribution] = get_internal_state ();;
613 }
614
615 int rand::get_dist_id (const std::string& d)
616 {
617 int retval = unknown_dist;
618
619 if (d == "uniform" || d == "rand")
620 retval = uniform_dist;
621 else if (d == "normal" || d == "randn")
622 retval = normal_dist;
623 else if (d == "exponential" || d == "rande")
624 retval = expon_dist;
625 else if (d == "poisson" || d == "randp")
626 retval = poisson_dist;
627 else if (d == "gamma" || d == "randg")
628 retval = gamma_dist;
629 else
630 (*current_liboctave_error_handler)
631 ("rand: invalid distribution '%s'", d.c_str ());
632
633 return retval;
634 }
635
636 void rand::set_internal_state (const uint32NDArray& s)
637 {
638 octave_idx_type len = s.numel ();
639
640 const uint32_t *sdata = reinterpret_cast <const uint32_t *> (s.data ());
641
642 if (len == MT_N + 1 && sdata[MT_N] <= MT_N && sdata[MT_N] > 0)
643 set_mersenne_twister_state (sdata);
644 else
645 init_mersenne_twister (sdata, len);
646 }
647
648 void rand::switch_to_generator (int dist)
649 {
650 if (dist != m_current_distribution)
651 {
652 m_current_distribution = dist;
653
654 set_internal_state (m_rand_states[dist]);
655 }
656 }
657
658 void rand::fill (octave_idx_type len, double *v, double a)
659 {
660 if (len < 1)
723 return; 661 return;
724 } 662
725 663 switch (m_current_distribution)
726 void rand::fill (octave_idx_type len, float *v, float a) 664 {
727 { 665 case uniform_dist:
728 if (len < 1) 666 if (m_use_old_generators)
729 return; 667 std::generate_n (v, len, [](void) { double x; F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, x); return x; });
730 668 else
731 switch (m_current_distribution) 669 rand_uniform<double> (len, v);
732 { 670 break;
733 case uniform_dist: 671
734 if (m_use_old_generators) 672 case normal_dist:
735 std::generate_n (v, len, [](void) { float x; F77_FUNC (fgenunf, FGENUNF) (0.0f, 1.0f, x); return x; }); 673 if (m_use_old_generators)
736 else 674 std::generate_n (v, len, [](void) { double x; F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, x); return x; });
737 rand_uniform<float> (len, v); 675 else
738 break; 676 rand_normal<double> (len, v);
739 677 break;
740 case normal_dist: 678
741 if (m_use_old_generators) 679 case expon_dist:
742 std::generate_n (v, len, [](void) { float x; F77_FUNC (fgennor, FGENNOR) (0.0f, 1.0f, x); return x; }); 680 if (m_use_old_generators)
743 else 681 std::generate_n (v, len, [](void) { double x; F77_FUNC (dgenexp, DGENEXP) (1.0, x); return x; });
744 rand_normal<float> (len, v); 682 else
745 break; 683 rand_exponential<double> (len, v);
746 684 break;
747 case expon_dist: 685
748 if (m_use_old_generators) 686 case poisson_dist:
749 std::generate_n (v, len, [](void) { float x; F77_FUNC (fgenexp, FGENEXP) (1.0f, x); return x; }); 687 if (m_use_old_generators)
750 else 688 {
751 rand_exponential<float> (len, v); 689 if (a < 0.0 || ! math::isfinite (a))
752 break; 690 std::fill_n (v, len, numeric_limits<double>::NaN ());
753 691 else
754 case poisson_dist: 692 {
755 if (m_use_old_generators) 693 // workaround bug in ignpoi, by calling with different Mu
756 { 694 double tmp;
757 if (a < 0.0f || ! math::isfinite (a)) 695 F77_FUNC (dignpoi, DIGNPOI) (a + 1, tmp);
758 std::fill_n (v, len, numeric_limits<float>::NaN ()); 696 std::generate_n (v, len, [a](void) { double x; F77_FUNC (dignpoi, DIGNPOI) (a, x); return x; });
759 else 697 }
760 { 698 }
761 // workaround bug in ignpoi, by calling with different Mu 699 else
762 float tmp; 700 rand_poisson<double> (a, len, v);
763 F77_FUNC (fignpoi, FIGNPOI) (a + 1, tmp); 701 break;
764 std::generate_n (v, len, [a](void) { float x; F77_FUNC (fignpoi, FIGNPOI) (a, x); return x; }); 702
765 } 703 case gamma_dist:
766 } 704 if (m_use_old_generators)
767 else 705 {
768 rand_poisson<float> (a, len, v); 706 if (a <= 0.0 || ! math::isfinite (a))
769 break; 707 std::fill_n (v, len, numeric_limits<double>::NaN ());
770 708 else
771 case gamma_dist: 709 std::generate_n (v, len, [a](void) { double x; F77_FUNC (dgengam, DGENGAM) (1.0, a, x); return x; });
772 if (m_use_old_generators) 710 }
773 { 711 else
774 if (a <= 0.0f || ! math::isfinite (a)) 712 rand_gamma<double> (a, len, v);
775 std::fill_n (v, len, numeric_limits<float>::NaN ()); 713 break;
776 else 714
777 std::generate_n (v, len, [a](void) { float x; F77_FUNC (fgengam, FGENGAM) (1.0f, a, x); return x; }); 715 default:
778 } 716 (*current_liboctave_error_handler)
779 else 717 ("rand: invalid distribution ID = %d", m_current_distribution);
780 rand_gamma<float> (a, len, v); 718 break;
781 break; 719 }
782 720
783 default: 721 save_state ();
784 (*current_liboctave_error_handler) 722
785 ("rand: invalid distribution ID = %d", m_current_distribution); 723 return;
786 break; 724 }
787 } 725
788 726 void rand::fill (octave_idx_type len, float *v, float a)
789 save_state (); 727 {
790 728 if (len < 1)
791 return; 729 return;
792 } 730
731 switch (m_current_distribution)
732 {
733 case uniform_dist:
734 if (m_use_old_generators)
735 std::generate_n (v, len, [](void) { float x; F77_FUNC (fgenunf, FGENUNF) (0.0f, 1.0f, x); return x; });
736 else
737 rand_uniform<float> (len, v);
738 break;
739
740 case normal_dist:
741 if (m_use_old_generators)
742 std::generate_n (v, len, [](void) { float x; F77_FUNC (fgennor, FGENNOR) (0.0f, 1.0f, x); return x; });
743 else
744 rand_normal<float> (len, v);
745 break;
746
747 case expon_dist:
748 if (m_use_old_generators)
749 std::generate_n (v, len, [](void) { float x; F77_FUNC (fgenexp, FGENEXP) (1.0f, x); return x; });
750 else
751 rand_exponential<float> (len, v);
752 break;
753
754 case poisson_dist:
755 if (m_use_old_generators)
756 {
757 if (a < 0.0f || ! math::isfinite (a))
758 std::fill_n (v, len, numeric_limits<float>::NaN ());
759 else
760 {
761 // workaround bug in ignpoi, by calling with different Mu
762 float tmp;
763 F77_FUNC (fignpoi, FIGNPOI) (a + 1, tmp);
764 std::generate_n (v, len, [a](void) { float x; F77_FUNC (fignpoi, FIGNPOI) (a, x); return x; });
765 }
766 }
767 else
768 rand_poisson<float> (a, len, v);
769 break;
770
771 case gamma_dist:
772 if (m_use_old_generators)
773 {
774 if (a <= 0.0f || ! math::isfinite (a))
775 std::fill_n (v, len, numeric_limits<float>::NaN ());
776 else
777 std::generate_n (v, len, [a](void) { float x; F77_FUNC (fgengam, FGENGAM) (1.0f, a, x); return x; });
778 }
779 else
780 rand_gamma<float> (a, len, v);
781 break;
782
783 default:
784 (*current_liboctave_error_handler)
785 ("rand: invalid distribution ID = %d", m_current_distribution);
786 break;
787 }
788
789 save_state ();
790
791 return;
792 }
793 793
794 OCTAVE_END_NAMESPACE(octave) 794 OCTAVE_END_NAMESPACE(octave)