File: | liboctave/numeric/oct-rand.cc |
Location: | line 458, column 13 |
Description: | Value stored to 'retval' is never read |
1 | /* |
2 | |
3 | Copyright (C) 2003-2013 John W. Eaton |
4 | |
5 | This file is part of Octave. |
6 | |
7 | Octave is free software; you can redistribute it and/or modify it |
8 | under the terms of the GNU General Public License as published by the |
9 | Free Software Foundation; either version 3 of the License, or (at your |
10 | option) any later version. |
11 | |
12 | Octave is distributed in the hope that it will be useful, but WITHOUT |
13 | ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
14 | FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
15 | for more details. |
16 | |
17 | You should have received a copy of the GNU General Public License |
18 | along with Octave; see the file COPYING. If not, see |
19 | <http://www.gnu.org/licenses/>. |
20 | |
21 | */ |
22 | |
23 | #ifdef HAVE_CONFIG_H1 |
24 | #include <config.h> |
25 | #endif |
26 | |
27 | #include <map> |
28 | #include <vector> |
29 | |
30 | #include <stdint.h> |
31 | |
32 | #include "data-conv.h" |
33 | #include "f77-fcn.h" |
34 | #include "lo-error.h" |
35 | #include "lo-ieee.h" |
36 | #include "lo-mappers.h" |
37 | #include "mach-info.h" |
38 | #include "oct-locbuf.h" |
39 | #include "oct-rand.h" |
40 | #include "oct-time.h" |
41 | #include "randgamma.h" |
42 | #include "randmtzig.h" |
43 | #include "randpoisson.h" |
44 | #include "singleton-cleanup.h" |
45 | |
46 | extern "C" |
47 | { |
48 | F77_RET_Tint |
49 | F77_FUNC (dgennor, DGENNOR)dgennor_ (const double&, const double&, double&); |
50 | |
51 | F77_RET_Tint |
52 | F77_FUNC (dgenunf, DGENUNF)dgenunf_ (const double&, const double&, double&); |
53 | |
54 | F77_RET_Tint |
55 | F77_FUNC (dgenexp, DGENEXP)dgenexp_ (const double&, double&); |
56 | |
57 | F77_RET_Tint |
58 | F77_FUNC (dignpoi, DIGNPOI)dignpoi_ (const double&, double&); |
59 | |
60 | F77_RET_Tint |
61 | F77_FUNC (dgengam, DGENGAM)dgengam_ (const double&, const double&, double&); |
62 | |
63 | F77_RET_Tint |
64 | F77_FUNC (setall, SETALL)setall_ (const int32_t&, const int32_t&); |
65 | |
66 | F77_RET_Tint |
67 | F77_FUNC (getsd, GETSD)getsd_ (int32_t&, int32_t&); |
68 | |
69 | F77_RET_Tint |
70 | F77_FUNC (setsd, SETSD)setsd_ (const int32_t&, const int32_t&); |
71 | |
72 | F77_RET_Tint |
73 | F77_FUNC (setcgn, SETCGN)setcgn_ (const int32_t&); |
74 | } |
75 | |
76 | octave_rand *octave_rand::instance = 0; |
77 | |
78 | octave_rand::octave_rand (void) |
79 | : current_distribution (uniform_dist), use_old_generators (false), |
80 | rand_states () |
81 | { |
82 | initialize_ranlib_generators (); |
83 | |
84 | initialize_mersenne_twister (); |
85 | } |
86 | |
87 | bool |
88 | octave_rand::instance_ok (void) |
89 | { |
90 | bool retval = true; |
91 | |
92 | if (! instance) |
93 | { |
94 | instance = new octave_rand (); |
95 | |
96 | if (instance) |
97 | singleton_cleanup_list::add (cleanup_instance); |
98 | } |
99 | |
100 | if (! instance) |
101 | { |
102 | (*current_liboctave_error_handler) |
103 | ("unable to create octave_rand object!"); |
104 | |
105 | retval = false; |
106 | } |
107 | |
108 | return retval; |
109 | } |
110 | |
111 | double |
112 | octave_rand::do_seed (void) |
113 | { |
114 | union d2i { double d; int32_t i[2]; }; |
115 | union d2i u; |
116 | |
117 | oct_mach_info::float_format ff = oct_mach_info::native_float_format (); |
118 | |
119 | switch (ff) |
120 | { |
121 | case oct_mach_info::flt_fmt_ieee_big_endian: |
122 | F77_FUNC (getsd, GETSD)getsd_ (u.i[1], u.i[0]); |
123 | break; |
124 | default: |
125 | F77_FUNC (getsd, GETSD)getsd_ (u.i[0], u.i[1]); |
126 | break; |
127 | } |
128 | |
129 | return u.d; |
130 | } |
131 | |
132 | static int32_t |
133 | force_to_fit_range (int32_t i, int32_t lo, int32_t hi) |
134 | { |
135 | assert (hi > lo && lo >= 0 && hi > lo)((hi > lo && lo >= 0 && hi > lo) ? static_cast <void> (0) : __assert_fail ("hi > lo && lo >= 0 && hi > lo" , "numeric/oct-rand.cc", 135, __PRETTY_FUNCTION__)); |
136 | |
137 | i = i > 0 ? i : -i; |
138 | |
139 | if (i < lo) |
140 | i = lo; |
141 | else if (i > hi) |
142 | i = i % hi; |
143 | |
144 | return i; |
145 | } |
146 | |
147 | void |
148 | octave_rand::do_seed (double s) |
149 | { |
150 | use_old_generators = true; |
151 | |
152 | int i0, i1; |
153 | union d2i { double d; int32_t i[2]; }; |
154 | union d2i u; |
155 | u.d = s; |
156 | |
157 | oct_mach_info::float_format ff = oct_mach_info::native_float_format (); |
158 | |
159 | switch (ff) |
160 | { |
161 | case oct_mach_info::flt_fmt_ieee_big_endian: |
162 | i1 = force_to_fit_range (u.i[0], 1, 2147483563); |
163 | i0 = force_to_fit_range (u.i[1], 1, 2147483399); |
164 | break; |
165 | default: |
166 | i0 = force_to_fit_range (u.i[0], 1, 2147483563); |
167 | i1 = force_to_fit_range (u.i[1], 1, 2147483399); |
168 | break; |
169 | } |
170 | |
171 | F77_FUNC (setsd, SETSD)setsd_ (i0, i1); |
172 | } |
173 | |
174 | void |
175 | octave_rand::do_reset (void) |
176 | { |
177 | use_old_generators = true; |
178 | initialize_ranlib_generators (); |
179 | } |
180 | |
181 | ColumnVector |
182 | octave_rand::do_state (const std::string& d) |
183 | { |
184 | return rand_states[d.empty () ? current_distribution : get_dist_id (d)]; |
185 | } |
186 | |
187 | void |
188 | octave_rand::do_state (const ColumnVector& s, const std::string& d) |
189 | { |
190 | use_old_generators = false; |
191 | |
192 | int old_dist = current_distribution; |
193 | |
194 | int new_dist = d.empty () ? current_distribution : get_dist_id (d); |
195 | |
196 | ColumnVector saved_state; |
197 | |
198 | if (old_dist != new_dist) |
199 | saved_state = get_internal_state (); |
200 | |
201 | set_internal_state (s); |
202 | |
203 | rand_states[new_dist] = get_internal_state (); |
204 | |
205 | if (old_dist != new_dist) |
206 | rand_states[old_dist] = saved_state; |
207 | } |
208 | |
209 | void |
210 | octave_rand::do_reset (const std::string& d) |
211 | { |
212 | use_old_generators = false; |
213 | |
214 | int old_dist = current_distribution; |
215 | |
216 | int new_dist = d.empty () ? current_distribution : get_dist_id (d); |
217 | |
218 | ColumnVector saved_state; |
219 | |
220 | if (old_dist != new_dist) |
221 | saved_state = get_internal_state (); |
222 | |
223 | oct_init_by_entropy (); |
224 | rand_states[new_dist] = get_internal_state (); |
225 | |
226 | if (old_dist != new_dist) |
227 | rand_states[old_dist] = saved_state; |
228 | } |
229 | |
230 | std::string |
231 | octave_rand::do_distribution (void) |
232 | { |
233 | std::string retval; |
234 | |
235 | switch (current_distribution) |
236 | { |
237 | case uniform_dist: |
238 | retval = "uniform"; |
239 | break; |
240 | |
241 | case normal_dist: |
242 | retval = "normal"; |
243 | break; |
244 | |
245 | case expon_dist: |
246 | retval = "exponential"; |
247 | break; |
248 | |
249 | case poisson_dist: |
250 | retval = "poisson"; |
251 | break; |
252 | |
253 | case gamma_dist: |
254 | retval = "gamma"; |
255 | break; |
256 | |
257 | default: |
258 | (*current_liboctave_error_handler) |
259 | ("rand: invalid distribution ID = %d", current_distribution); |
260 | break; |
261 | } |
262 | |
263 | return retval; |
264 | } |
265 | |
266 | void |
267 | octave_rand::do_distribution (const std::string& d) |
268 | { |
269 | int id = get_dist_id (d); |
270 | |
271 | switch (id) |
272 | { |
273 | case uniform_dist: |
274 | octave_rand::uniform_distribution (); |
275 | break; |
276 | |
277 | case normal_dist: |
278 | octave_rand::normal_distribution (); |
279 | break; |
280 | |
281 | case expon_dist: |
282 | octave_rand::exponential_distribution (); |
283 | break; |
284 | |
285 | case poisson_dist: |
286 | octave_rand::poisson_distribution (); |
287 | break; |
288 | |
289 | case gamma_dist: |
290 | octave_rand::gamma_distribution (); |
291 | break; |
292 | |
293 | default: |
294 | (*current_liboctave_error_handler) |
295 | ("rand: invalid distribution ID = %d", id); |
296 | break; |
297 | } |
298 | } |
299 | |
300 | void |
301 | octave_rand::do_uniform_distribution (void) |
302 | { |
303 | switch_to_generator (uniform_dist); |
304 | |
305 | F77_FUNC (setcgn, SETCGN)setcgn_ (uniform_dist); |
306 | } |
307 | |
308 | void |
309 | octave_rand::do_normal_distribution (void) |
310 | { |
311 | switch_to_generator (normal_dist); |
312 | |
313 | F77_FUNC (setcgn, SETCGN)setcgn_ (normal_dist); |
314 | } |
315 | |
316 | void |
317 | octave_rand::do_exponential_distribution (void) |
318 | { |
319 | switch_to_generator (expon_dist); |
320 | |
321 | F77_FUNC (setcgn, SETCGN)setcgn_ (expon_dist); |
322 | } |
323 | |
324 | void |
325 | octave_rand::do_poisson_distribution (void) |
326 | { |
327 | switch_to_generator (poisson_dist); |
328 | |
329 | F77_FUNC (setcgn, SETCGN)setcgn_ (poisson_dist); |
330 | } |
331 | |
332 | void |
333 | octave_rand::do_gamma_distribution (void) |
334 | { |
335 | switch_to_generator (gamma_dist); |
336 | |
337 | F77_FUNC (setcgn, SETCGN)setcgn_ (gamma_dist); |
338 | } |
339 | |
340 | |
341 | double |
342 | octave_rand::do_scalar (double a) |
343 | { |
344 | double retval = 0.0; |
345 | |
346 | if (use_old_generators) |
347 | { |
348 | switch (current_distribution) |
349 | { |
350 | case uniform_dist: |
351 | F77_FUNC (dgenunf, DGENUNF)dgenunf_ (0.0, 1.0, retval); |
352 | break; |
353 | |
354 | case normal_dist: |
355 | F77_FUNC (dgennor, DGENNOR)dgennor_ (0.0, 1.0, retval); |
356 | break; |
357 | |
358 | case expon_dist: |
359 | F77_FUNC (dgenexp, DGENEXP)dgenexp_ (1.0, retval); |
360 | break; |
361 | |
362 | case poisson_dist: |
363 | if (a < 0.0 || ! xfinite (a)) |
364 | retval = octave_NaN; |
365 | else |
366 | { |
367 | // workaround bug in ignpoi, by calling with different Mu |
368 | F77_FUNC (dignpoi, DIGNPOI)dignpoi_ (a + 1, retval); |
369 | F77_FUNC (dignpoi, DIGNPOI)dignpoi_ (a, retval); |
370 | } |
371 | break; |
372 | |
373 | case gamma_dist: |
374 | if (a <= 0.0 || ! xfinite (a)) |
375 | retval = octave_NaN; |
376 | else |
377 | F77_FUNC (dgengam, DGENGAM)dgengam_ (1.0, a, retval); |
378 | break; |
379 | |
380 | default: |
381 | (*current_liboctave_error_handler) |
382 | ("rand: invalid distribution ID = %d", current_distribution); |
383 | break; |
384 | } |
385 | } |
386 | else |
387 | { |
388 | switch (current_distribution) |
389 | { |
390 | case uniform_dist: |
391 | retval = oct_randu (); |
392 | break; |
393 | |
394 | case normal_dist: |
395 | retval = oct_randn (); |
396 | break; |
397 | |
398 | case expon_dist: |
399 | retval = oct_rande (); |
400 | break; |
401 | |
402 | case poisson_dist: |
403 | retval = oct_randp (a); |
404 | break; |
405 | |
406 | case gamma_dist: |
407 | retval = oct_randg (a); |
408 | break; |
409 | |
410 | default: |
411 | (*current_liboctave_error_handler) |
412 | ("rand: invalid distribution ID = %d", current_distribution); |
413 | break; |
414 | } |
415 | |
416 | save_state (); |
417 | } |
418 | |
419 | return retval; |
420 | } |
421 | |
422 | float |
423 | octave_rand::do_float_scalar (float a) |
424 | { |
425 | float retval = 0.0; |
426 | |
427 | if (use_old_generators) |
428 | { |
429 | double da = a; |
430 | double dretval = 0.0; |
431 | switch (current_distribution) |
432 | { |
433 | case uniform_dist: |
434 | F77_FUNC (dgenunf, DGENUNF)dgenunf_ (0.0, 1.0, dretval); |
435 | break; |
436 | |
437 | case normal_dist: |
438 | F77_FUNC (dgennor, DGENNOR)dgennor_ (0.0, 1.0, dretval); |
439 | break; |
440 | |
441 | case expon_dist: |
442 | F77_FUNC (dgenexp, DGENEXP)dgenexp_ (1.0, dretval); |
443 | break; |
444 | |
445 | case poisson_dist: |
446 | if (da < 0.0 || ! xfinite (a)) |
447 | dretval = octave_NaN; |
448 | else |
449 | { |
450 | // workaround bug in ignpoi, by calling with different Mu |
451 | F77_FUNC (dignpoi, DIGNPOI)dignpoi_ (da + 1, dretval); |
452 | F77_FUNC (dignpoi, DIGNPOI)dignpoi_ (da, dretval); |
453 | } |
454 | break; |
455 | |
456 | case gamma_dist: |
457 | if (da <= 0.0 || ! xfinite (a)) |
458 | retval = octave_NaN; |
Value stored to 'retval' is never read | |
459 | else |
460 | F77_FUNC (dgengam, DGENGAM)dgengam_ (1.0, da, dretval); |
461 | break; |
462 | |
463 | default: |
464 | (*current_liboctave_error_handler) |
465 | ("rand: invalid distribution ID = %d", current_distribution); |
466 | break; |
467 | } |
468 | retval = dretval; |
469 | } |
470 | else |
471 | { |
472 | switch (current_distribution) |
473 | { |
474 | case uniform_dist: |
475 | retval = oct_float_randu (); |
476 | break; |
477 | |
478 | case normal_dist: |
479 | retval = oct_float_randn (); |
480 | break; |
481 | |
482 | case expon_dist: |
483 | retval = oct_float_rande (); |
484 | break; |
485 | |
486 | case poisson_dist: |
487 | // Keep poisson distribution in double precision for accuracy |
488 | retval = oct_randp (a); |
489 | break; |
490 | |
491 | case gamma_dist: |
492 | retval = oct_float_randg (a); |
493 | break; |
494 | |
495 | default: |
496 | (*current_liboctave_error_handler) |
497 | ("rand: invalid distribution ID = %d", current_distribution); |
498 | break; |
499 | } |
500 | |
501 | save_state (); |
502 | } |
503 | |
504 | return retval; |
505 | } |
506 | |
507 | Array<double> |
508 | octave_rand::do_vector (octave_idx_type n, double a) |
509 | { |
510 | Array<double> retval; |
511 | |
512 | if (n > 0) |
513 | { |
514 | retval.clear (n, 1); |
515 | |
516 | fill (retval.capacity (), retval.fortran_vec (), a); |
517 | } |
518 | else if (n < 0) |
519 | (*current_liboctave_error_handler) ("rand: invalid negative argument"); |
520 | |
521 | return retval; |
522 | } |
523 | |
524 | Array<float> |
525 | octave_rand::do_float_vector (octave_idx_type n, float a) |
526 | { |
527 | Array<float> retval; |
528 | |
529 | if (n > 0) |
530 | { |
531 | retval.clear (n, 1); |
532 | |
533 | fill (retval.capacity (), retval.fortran_vec (), a); |
534 | } |
535 | else if (n < 0) |
536 | (*current_liboctave_error_handler) ("rand: invalid negative argument"); |
537 | |
538 | return retval; |
539 | } |
540 | |
541 | NDArray |
542 | octave_rand::do_nd_array (const dim_vector& dims, double a) |
543 | { |
544 | NDArray retval; |
545 | |
546 | if (! dims.all_zero ()) |
547 | { |
548 | retval.clear (dims); |
549 | |
550 | fill (retval.capacity (), retval.fortran_vec (), a); |
551 | } |
552 | |
553 | return retval; |
554 | } |
555 | |
556 | FloatNDArray |
557 | octave_rand::do_float_nd_array (const dim_vector& dims, float a) |
558 | { |
559 | FloatNDArray retval; |
560 | |
561 | if (! dims.all_zero ()) |
562 | { |
563 | retval.clear (dims); |
564 | |
565 | fill (retval.capacity (), retval.fortran_vec (), a); |
566 | } |
567 | |
568 | return retval; |
569 | } |
570 | |
571 | // Make the random number generator give us a different sequence every |
572 | // time we start octave unless we specifically set the seed. The |
573 | // technique used below will cycle monthly, but it it does seem to |
574 | // work ok to give fairly different seeds each time Octave starts. |
575 | |
576 | void |
577 | octave_rand::initialize_ranlib_generators (void) |
578 | { |
579 | octave_localtime tm; |
580 | int stored_distribution = current_distribution; |
581 | F77_FUNC (setcgn, SETCGN)setcgn_ (uniform_dist); |
582 | |
583 | int hour = tm.hour () + 1; |
584 | int minute = tm.min () + 1; |
585 | int second = tm.sec () + 1; |
586 | |
587 | int32_t s0 = tm.mday () * hour * minute * second; |
588 | int32_t s1 = hour * minute * second; |
589 | |
590 | s0 = force_to_fit_range (s0, 1, 2147483563); |
591 | s1 = force_to_fit_range (s1, 1, 2147483399); |
592 | |
593 | F77_FUNC (setall, SETALL)setall_ (s0, s1); |
594 | F77_FUNC (setcgn, SETCGN)setcgn_ (stored_distribution); |
595 | } |
596 | |
597 | void |
598 | octave_rand::initialize_mersenne_twister (void) |
599 | { |
600 | oct_init_by_entropy (); |
601 | |
602 | ColumnVector s = get_internal_state (); |
603 | |
604 | rand_states[uniform_dist] = s; |
605 | |
606 | oct_init_by_entropy (); |
607 | s = get_internal_state (); |
608 | rand_states[normal_dist] = s; |
609 | |
610 | oct_init_by_entropy (); |
611 | s = get_internal_state (); |
612 | rand_states[expon_dist] = s; |
613 | |
614 | oct_init_by_entropy (); |
615 | s = get_internal_state (); |
616 | rand_states[poisson_dist] = s; |
617 | |
618 | oct_init_by_entropy (); |
619 | s = get_internal_state (); |
620 | rand_states[gamma_dist] = s; |
621 | } |
622 | |
623 | ColumnVector |
624 | octave_rand::get_internal_state (void) |
625 | { |
626 | ColumnVector s (MT_N624 + 1); |
627 | |
628 | OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1)octave_local_buffer<uint32_t> _buffer_tmp (624 + 1); uint32_t *tmp = _buffer_tmp; |
629 | |
630 | oct_get_state (tmp); |
631 | |
632 | for (octave_idx_type i = 0; i <= MT_N624; i++) |
633 | s.elem (i) = static_cast<double> (tmp[i]); |
634 | |
635 | return s; |
636 | } |
637 | |
638 | void |
639 | octave_rand::save_state (void) |
640 | { |
641 | rand_states[current_distribution] = get_internal_state ();; |
642 | } |
643 | |
644 | int |
645 | octave_rand::get_dist_id (const std::string& d) |
646 | { |
647 | int retval = unknown_dist; |
648 | |
649 | if (d == "uniform" || d == "rand") |
650 | retval = uniform_dist; |
651 | else if (d == "normal" || d == "randn") |
652 | retval = normal_dist; |
653 | else if (d == "exponential" || d == "rande") |
654 | retval = expon_dist; |
655 | else if (d == "poisson" || d == "randp") |
656 | retval = poisson_dist; |
657 | else if (d == "gamma" || d == "randg") |
658 | retval = gamma_dist; |
659 | else |
660 | (*current_liboctave_error_handler) |
661 | ("rand: invalid distribution '%s'", d.c_str ()); |
662 | |
663 | return retval; |
664 | } |
665 | |
666 | // Guarantee reproducible conversion of negative initialization values to |
667 | // random number algorithm. Note that Matlab employs slightly different rules. |
668 | // 1) Seed saturates at 2^32-1 for any value larger than that. |
669 | // 2) NaN, Inf are translated to 2^32-1. |
670 | // 3) -Inf is translated to 0. |
671 | static uint32_t |
672 | double2uint32 (double d) |
673 | { |
674 | uint32_t u; |
675 | static const double TWOUP32 = std::numeric_limits<uint32_t>::max() + 1.0; |
676 | |
677 | if (! xfinite (d)) |
678 | u = 0; |
679 | else |
680 | { |
681 | d = fmod (d, TWOUP32); |
682 | if (d < 0) |
683 | d += TWOUP32; |
684 | u = static_cast<uint32_t> (d); |
685 | } |
686 | |
687 | return u; |
688 | } |
689 | |
690 | void |
691 | octave_rand::set_internal_state (const ColumnVector& s) |
692 | { |
693 | octave_idx_type len = s.length (); |
694 | octave_idx_type n = len < MT_N624 + 1 ? len : MT_N624 + 1; |
695 | |
696 | OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1)octave_local_buffer<uint32_t> _buffer_tmp (624 + 1); uint32_t *tmp = _buffer_tmp; |
697 | |
698 | for (octave_idx_type i = 0; i < n; i++) |
699 | tmp[i] = double2uint32 (s.elem (i)); |
700 | |
701 | if (len == MT_N624 + 1 && tmp[MT_N624] <= MT_N624 && tmp[MT_N624] > 0) |
702 | oct_set_state (tmp); |
703 | else |
704 | oct_init_by_array (tmp, len); |
705 | } |
706 | |
707 | void |
708 | octave_rand::switch_to_generator (int dist) |
709 | { |
710 | if (dist != current_distribution) |
711 | { |
712 | current_distribution = dist; |
713 | |
714 | set_internal_state (rand_states[dist]); |
715 | } |
716 | } |
717 | |
718 | #define MAKE_RAND(len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0) \ |
719 | do \ |
720 | { \ |
721 | double val; \ |
722 | for (volatile octave_idx_type i = 0; i < len; i++) \ |
723 | { \ |
724 | octave_quit (); \ |
725 | RAND_FUNC (val); \ |
726 | v[i] = val; \ |
727 | } \ |
728 | } \ |
729 | while (0) |
730 | |
731 | void |
732 | octave_rand::fill (octave_idx_type len, double *v, double a) |
733 | { |
734 | if (len < 1) |
735 | return; |
736 | |
737 | switch (current_distribution) |
738 | { |
739 | case uniform_dist: |
740 | if (use_old_generators) |
741 | { |
742 | #define RAND_FUNC(x) F77_FUNC (dgenunf, DGENUNF)dgenunf_ (0.0, 1.0, x) |
743 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
744 | #undef RAND_FUNC |
745 | } |
746 | else |
747 | oct_fill_randu (len, v); |
748 | break; |
749 | |
750 | case normal_dist: |
751 | if (use_old_generators) |
752 | { |
753 | #define RAND_FUNC(x) F77_FUNC (dgennor, DGENNOR)dgennor_ (0.0, 1.0, x) |
754 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
755 | #undef RAND_FUNC |
756 | } |
757 | else |
758 | oct_fill_randn (len, v); |
759 | break; |
760 | |
761 | case expon_dist: |
762 | if (use_old_generators) |
763 | { |
764 | #define RAND_FUNC(x) F77_FUNC (dgenexp, DGENEXP)dgenexp_ (1.0, x) |
765 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
766 | #undef RAND_FUNC |
767 | } |
768 | else |
769 | oct_fill_rande (len, v); |
770 | break; |
771 | |
772 | case poisson_dist: |
773 | if (use_old_generators) |
774 | { |
775 | if (a < 0.0 || ! xfinite (a)) |
776 | #define RAND_FUNC(x) x = octave_NaN; |
777 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
778 | #undef RAND_FUNC |
779 | else |
780 | { |
781 | // workaround bug in ignpoi, by calling with different Mu |
782 | double tmp; |
783 | F77_FUNC (dignpoi, DIGNPOI)dignpoi_ (a + 1, tmp); |
784 | #define RAND_FUNC(x) F77_FUNC (dignpoi, DIGNPOI)dignpoi_ (a, x) |
785 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
786 | #undef RAND_FUNC |
787 | } |
788 | } |
789 | else |
790 | oct_fill_randp (a, len, v); |
791 | break; |
792 | |
793 | case gamma_dist: |
794 | if (use_old_generators) |
795 | { |
796 | if (a <= 0.0 || ! xfinite (a)) |
797 | #define RAND_FUNC(x) x = octave_NaN; |
798 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
799 | #undef RAND_FUNC |
800 | else |
801 | #define RAND_FUNC(x) F77_FUNC (dgengam, DGENGAM)dgengam_ (1.0, a, x) |
802 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
803 | #undef RAND_FUNC |
804 | } |
805 | else |
806 | oct_fill_randg (a, len, v); |
807 | break; |
808 | |
809 | default: |
810 | (*current_liboctave_error_handler) |
811 | ("rand: invalid distribution ID = %d", current_distribution); |
812 | break; |
813 | } |
814 | |
815 | save_state (); |
816 | |
817 | return; |
818 | } |
819 | |
820 | void |
821 | octave_rand::fill (octave_idx_type len, float *v, float a) |
822 | { |
823 | if (len < 1) |
824 | return; |
825 | |
826 | switch (current_distribution) |
827 | { |
828 | case uniform_dist: |
829 | if (use_old_generators) |
830 | { |
831 | #define RAND_FUNC(x) F77_FUNC (dgenunf, DGENUNF)dgenunf_ (0.0, 1.0, x) |
832 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
833 | #undef RAND_FUNC |
834 | } |
835 | else |
836 | oct_fill_float_randu (len, v); |
837 | break; |
838 | |
839 | case normal_dist: |
840 | if (use_old_generators) |
841 | { |
842 | #define RAND_FUNC(x) F77_FUNC (dgennor, DGENNOR)dgennor_ (0.0, 1.0, x) |
843 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
844 | #undef RAND_FUNC |
845 | } |
846 | else |
847 | oct_fill_float_randn (len, v); |
848 | break; |
849 | |
850 | case expon_dist: |
851 | if (use_old_generators) |
852 | { |
853 | #define RAND_FUNC(x) F77_FUNC (dgenexp, DGENEXP)dgenexp_ (1.0, x) |
854 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
855 | #undef RAND_FUNC |
856 | } |
857 | else |
858 | oct_fill_float_rande (len, v); |
859 | break; |
860 | |
861 | case poisson_dist: |
862 | if (use_old_generators) |
863 | { |
864 | double da = a; |
865 | if (da < 0.0 || ! xfinite (a)) |
866 | #define RAND_FUNC(x) x = octave_NaN; |
867 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
868 | #undef RAND_FUNC |
869 | else |
870 | { |
871 | // workaround bug in ignpoi, by calling with different Mu |
872 | double tmp; |
873 | F77_FUNC (dignpoi, DIGNPOI)dignpoi_ (da + 1, tmp); |
874 | #define RAND_FUNC(x) F77_FUNC (dignpoi, DIGNPOI)dignpoi_ (da, x) |
875 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
876 | #undef RAND_FUNC |
877 | } |
878 | } |
879 | else |
880 | oct_fill_float_randp (a, len, v); |
881 | break; |
882 | |
883 | case gamma_dist: |
884 | if (use_old_generators) |
885 | { |
886 | double da = a; |
887 | if (da <= 0.0 || ! xfinite (a)) |
888 | #define RAND_FUNC(x) x = octave_NaN; |
889 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
890 | #undef RAND_FUNC |
891 | else |
892 | #define RAND_FUNC(x) F77_FUNC (dgengam, DGENGAM)dgengam_ (1.0, da, x) |
893 | MAKE_RAND (len)do { double val; for (volatile octave_idx_type i = 0; i < len ; i++) { octave_quit (); RAND_FUNC (val); v[i] = val; } } while (0); |
894 | #undef RAND_FUNC |
895 | } |
896 | else |
897 | oct_fill_float_randg (a, len, v); |
898 | break; |
899 | |
900 | default: |
901 | (*current_liboctave_error_handler) |
902 | ("rand: invalid distribution ID = %d", current_distribution); |
903 | break; |
904 | } |
905 | |
906 | save_state (); |
907 | |
908 | return; |
909 | } |