Bug Summary

File:liboctave/numeric/oct-rand.cc
Location:line 458, column 13
Description:Value stored to 'retval' is never read

Annotated Source Code

1/*
2
3Copyright (C) 2003-2013 John W. Eaton
4
5This file is part of Octave.
6
7Octave is free software; you can redistribute it and/or modify it
8under the terms of the GNU General Public License as published by the
9Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12Octave is distributed in the hope that it will be useful, but WITHOUT
13ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15for more details.
16
17You should have received a copy of the GNU General Public License
18along 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
46extern "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
76octave_rand *octave_rand::instance = 0;
77
78octave_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
87bool
88octave_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
111double
112octave_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
132static int32_t
133force_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
147void
148octave_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
174void
175octave_rand::do_reset (void)
176{
177 use_old_generators = true;
178 initialize_ranlib_generators ();
179}
180
181ColumnVector
182octave_rand::do_state (const std::string& d)
183{
184 return rand_states[d.empty () ? current_distribution : get_dist_id (d)];
185}
186
187void
188octave_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
209void
210octave_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
230std::string
231octave_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
266void
267octave_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
300void
301octave_rand::do_uniform_distribution (void)
302{
303 switch_to_generator (uniform_dist);
304
305 F77_FUNC (setcgn, SETCGN)setcgn_ (uniform_dist);
306}
307
308void
309octave_rand::do_normal_distribution (void)
310{
311 switch_to_generator (normal_dist);
312
313 F77_FUNC (setcgn, SETCGN)setcgn_ (normal_dist);
314}
315
316void
317octave_rand::do_exponential_distribution (void)
318{
319 switch_to_generator (expon_dist);
320
321 F77_FUNC (setcgn, SETCGN)setcgn_ (expon_dist);
322}
323
324void
325octave_rand::do_poisson_distribution (void)
326{
327 switch_to_generator (poisson_dist);
328
329 F77_FUNC (setcgn, SETCGN)setcgn_ (poisson_dist);
330}
331
332void
333octave_rand::do_gamma_distribution (void)
334{
335 switch_to_generator (gamma_dist);
336
337 F77_FUNC (setcgn, SETCGN)setcgn_ (gamma_dist);
338}
339
340
341double
342octave_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
422float
423octave_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
507Array<double>
508octave_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
524Array<float>
525octave_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
541NDArray
542octave_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
556FloatNDArray
557octave_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
576void
577octave_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
597void
598octave_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
623ColumnVector
624octave_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
638void
639octave_rand::save_state (void)
640{
641 rand_states[current_distribution] = get_internal_state ();;
642}
643
644int
645octave_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.
671static uint32_t
672double2uint32 (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
690void
691octave_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
707void
708octave_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
731void
732octave_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
820void
821octave_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}