Bug Summary

File:liboctave/numeric/randpoisson.c
Location:line 399, column 13
Description:The right operand of '<=' is a garbage value

Annotated Source Code

1/*
2
3Copyright (C) 2006-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/* Original version written by Paul Kienzle distributed as free
24 software in the in the public domain. */
25
26/* Needs the following defines:
27 * NAN: value to return for Not-A-Number
28 * RUNI: uniform generator on (0,1)
29 * RNOR: normal generator
30 * LGAMMA: log gamma function
31 * INFINITE: function to test whether a value is infinite
32 */
33
34#if defined (HAVE_CONFIG_H1)
35#include <config.h>
36#endif
37
38#include <stdio.h>
39
40#include "f77-fcn.h"
41#include "lo-error.h"
42#include "lo-ieee.h"
43#include "lo-math.h"
44#include "randmtzig.h"
45#include "randpoisson.h"
46
47#undef NANoctave_NaN
48#define NANoctave_NaN octave_NaN
49#undef INFINITElo_ieee_isinf
50#define INFINITElo_ieee_isinf lo_ieee_isinf
51#define RUNIoct_randu() oct_randu()
52#define RNORoct_randn() oct_randn()
53#define LGAMMAxlgamma xlgamma
54
55F77_RET_Tint
56F77_FUNC (dlgams, DLGAMS)dlgams_ (const double *, double *, double *);
57
58static double
59xlgamma (double x)
60{
61 double result;
62#ifdef HAVE_LGAMMA1
63 result = lgamma (x);
64#else
65 double sgngam;
66
67 if (lo_ieee_isnan (x)(sizeof (x) == sizeof (float) ? __lo_ieee_float_isnan (x) : __lo_ieee_isnan
(x))
)
68 result = x;
69 else if (x <= 0 || lo_ieee_isinf (x)(sizeof (x) == sizeof (float) ? __lo_ieee_float_isinf (x) : __lo_ieee_isinf
(x))
)
70 result = octave_Inf;
71 else
72 F77_XFCN (dlgams, DLGAMS, (&x, &result, &sgngam))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dlgams_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dlgams_ (&x, &result, &sgngam); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
73#endif
74 return result;
75}
76
77/* ---- pprsc.c from Stadloeber's winrand --- */
78
79/* flogfak(k) = ln(k!) */
80static double
81flogfak (double k)
82{
83#define C09.18938533204672742e-01 9.18938533204672742e-01
84#define C18.33333333333333333e-02 8.33333333333333333e-02
85#define C3-2.77777777777777778e-03 -2.77777777777777778e-03
86#define C57.93650793650793651e-04 7.93650793650793651e-04
87#define C7-5.95238095238095238e-04 -5.95238095238095238e-04
88
89 static double logfak[30L] =
90 {
91 0.00000000000000000, 0.00000000000000000, 0.69314718055994531,
92 1.79175946922805500, 3.17805383034794562, 4.78749174278204599,
93 6.57925121201010100, 8.52516136106541430, 10.60460290274525023,
94 12.80182748008146961, 15.10441257307551530, 17.50230784587388584,
95 19.98721449566188615, 22.55216385312342289, 25.19122118273868150,
96 27.89927138384089157, 30.67186010608067280, 33.50507345013688888,
97 36.39544520803305358, 39.33988418719949404, 42.33561646075348503,
98 45.38013889847690803, 48.47118135183522388, 51.60667556776437357,
99 54.78472939811231919, 58.00360522298051994, 61.26170176100200198,
100 64.55753862700633106, 67.88974313718153498, 71.25703896716800901
101 };
102
103 double r, rr;
104
105 if (k >= 30.0)
106 {
107 r = 1.0 / k;
108 rr = r * r;
109 return ((k + 0.5)*log (k) - k + C09.18938533204672742e-01 + r*(C18.33333333333333333e-02 + rr*(C3-2.77777777777777778e-03 + rr*(C57.93650793650793651e-04 + rr*C7-5.95238095238095238e-04))));
110 }
111 else
112 return (logfak[(int)k]);
113}
114
115
116/******************************************************************
117 * *
118 * Poisson Distribution - Patchwork Rejection/Inversion *
119 * *
120 ******************************************************************
121 * *
122 * For parameter my < 10 Tabulated Inversion is applied. *
123 * For my >= 10 Patchwork Rejection is employed: *
124 * The area below the histogram function f(x) is rearranged in *
125 * its body by certain point reflections. Within a large center *
126 * interval variates are sampled efficiently by rejection from *
127 * uniform hats. Rectangular immediate acceptance regions speed *
128 * up the generation. The remaining tails are covered by *
129 * exponential functions. *
130 * *
131 ******************************************************************
132 * *
133 * FUNCTION : - pprsc samples a random number from the Poisson *
134 * distribution with parameter my > 0. *
135 * REFERENCE : - H. Zechner (1994): Efficient sampling from *
136 * continuous and discrete unimodal distributions, *
137 * Doctoral Dissertation, 156 pp., Technical *
138 * University Graz, Austria. *
139 * SUBPROGRAM : - drand(seed) ... (0,1)-Uniform generator with *
140 * unsigned long integer *seed. *
141 * *
142 * Implemented by H. Zechner, January 1994 *
143 * Revised by F. Niederl, July 1994 *
144 * *
145 ******************************************************************/
146
147static double
148f (double k, double l_nu, double c_pm)
149{
150 return exp (k * l_nu - flogfak (k) - c_pm);
151}
152
153static double
154pprsc (double my)
155{
156 static double my_last = -1.0;
157 static double m, k2, k4, k1, k5;
158 static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
159 f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
160 double Dk, X, Y;
161 double Ds, U, V, W;
162
163 if (my != my_last)
164 { /* set-up */
165 my_last = my;
166 /* approximate deviation of reflection points k2, k4 from my - 1/2 */
167 Ds = sqrt (my + 0.25);
168
169 /* mode m, reflection points k2 and k4, and points k1 and k5, */
170 /* which delimit the centre region of h(x) */
171 m = floor (my);
172 k2 = ceil (my - 0.5 - Ds);
173 k4 = floor (my - 0.5 + Ds);
174 k1 = k2 + k2 - m + 1L;
175 k5 = k4 + k4 - m;
176
177 /* range width of the critical left and right centre region */
178 dl = (k2 - k1);
179 dr = (k5 - k4);
180
181 /* recurrence constants r(k)=p(k)/p(k-1) at k = k1, k2, k4+1, k5+1 */
182 r1 = my / k1;
183 r2 = my / k2;
184 r4 = my / (k4 + 1.0);
185 r5 = my / (k5 + 1.0);
186
187 /* reciprocal values of the scale parameters of exp. tail envelope */
188 ll = log (r1); /* expon. tail left */
189 lr = -log (r5); /* expon. tail right*/
190
191 /* Poisson constants, necessary for computing function values f(k) */
192 l_my = log (my);
193 c_pm = m * l_my - flogfak (m);
194
195 /* function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5 */
196 f2 = f (k2, l_my, c_pm);
197 f4 = f (k4, l_my, c_pm);
198 f1 = f (k1, l_my, c_pm);
199 f5 = f (k5, l_my, c_pm);
200
201 /* area of the two centre and the two exponential tail regions */
202 /* area of the two immediate acceptance regions between k2, k4 */
203 p1 = f2 * (dl + 1.0); /* immed. left */
204 p2 = f2 * dl + p1; /* centre left */
205 p3 = f4 * (dr + 1.0) + p2; /* immed. right */
206 p4 = f4 * dr + p3; /* centre right */
207 p5 = f1 / ll + p4; /* exp. tail left */
208 p6 = f5 / lr + p5; /* exp. tail right*/
209 }
210
211 for (;;)
212 {
213 /* generate uniform number U -- U(0, p6) */
214 /* case distinction corresponding to U */
215 if ((U = RUNIoct_randu() * p6) < p2)
216 { /* centre left */
217
218 /* immediate acceptance region
219 R2 = [k2, m) *[0, f2), X = k2, ... m -1 */
220 if ((V = U - p1) < 0.0) return (k2 + floor (U/f2));
221 /* immediate acceptance region
222 R1 = [k1, k2)*[0, f1), X = k1, ... k2-1 */
223 if ((W = V / dl) < f1 ) return (k1 + floor (V/f1));
224
225 /* computation of candidate X < k2, and its counterpart Y > k2 */
226 /* either squeeze-acceptance of X or acceptance-rejection of Y */
227 Dk = floor (dl * RUNIoct_randu()) + 1.0;
228 if (W <= f2 - Dk * (f2 - f2/r2))
229 { /* quick accept of */
230 return (k2 - Dk); /* X = k2 - Dk */
231 }
232 if ((V = f2 + f2 - W) < 1.0)
233 { /* quick reject of Y*/
234 Y = k2 + Dk;
235 if (V <= f2 + Dk * (1.0 - f2)/(dl + 1.0))
236 { /* quick accept of */
237 return (Y); /* Y = k2 + Dk */
238 }
239 if (V <= f (Y, l_my, c_pm)) return (Y); /* final accept of Y*/
240 }
241 X = k2 - Dk;
242 }
243 else if (U < p4)
244 { /* centre right */
245 /* immediate acceptance region
246 R3 = [m, k4+1)*[0, f4), X = m, ... k4 */
247 if ((V = U - p3) < 0.0) return (k4 - floor ((U - p2)/f4));
248 /* immediate acceptance region
249 R4 = [k4+1, k5+1)*[0, f5) */
250 if ((W = V / dr) < f5 ) return (k5 - floor (V/f5));
251
252 /* computation of candidate X > k4, and its counterpart Y < k4 */
253 /* either squeeze-acceptance of X or acceptance-rejection of Y */
254 Dk = floor (dr * RUNIoct_randu()) + 1.0;
255 if (W <= f4 - Dk * (f4 - f4*r4))
256 { /* quick accept of */
257 return (k4 + Dk); /* X = k4 + Dk */
258 }
259 if ((V = f4 + f4 - W) < 1.0)
260 { /* quick reject of Y*/
261 Y = k4 - Dk;
262 if (V <= f4 + Dk * (1.0 - f4)/ dr)
263 { /* quick accept of */
264 return (Y); /* Y = k4 - Dk */
265 }
266 if (V <= f (Y, l_my, c_pm)) return (Y); /* final accept of Y*/
267 }
268 X = k4 + Dk;
269 }
270 else
271 {
272 W = RUNIoct_randu();
273 if (U < p5)
274 { /* expon. tail left */
275 Dk = floor (1.0 - log (W)/ll);
276 if ((X = k1 - Dk) < 0L) continue; /* 0 <= X <= k1 - 1 */
277 W *= (U - p4) * ll; /* W -- U(0, h(x)) */
278 if (W <= f1 - Dk * (f1 - f1/r1))
279 return (X); /* quick accept of X*/
280 }
281 else
282 { /* expon. tail right*/
283 Dk = floor (1.0 - log (W)/lr);
284 X = k5 + Dk; /* X >= k5 + 1 */
285 W *= (U - p5) * lr; /* W -- U(0, h(x)) */
286 if (W <= f5 - Dk * (f5 - f5*r5))
287 return (X); /* quick accept of X*/
288 }
289 }
290
291 /* acceptance-rejection test of candidate X from the original area */
292 /* test, whether W <= f(k), with W = U*h(x) and U -- U(0, 1)*/
293 /* log f(X) = (X - m)*log(my) - log X! + log m! */
294 if (log (W) <= X * l_my - flogfak (X) - c_pm) return (X);
295 }
296}
297/* ---- pprsc.c end ------ */
298
299
300/* The remainder of the file is by Paul Kienzle */
301
302/* Given uniform u, find x such that CDF(L,x)==u. Return x. */
303static void
304poisson_cdf_lookup (double lambda, double *p, size_t n)
305{
306 /* Table size is predicated on the maximum value of lambda
307 * we want to store in the table, and the maximum value of
308 * returned by the uniform random number generator on [0,1).
309 * With lambda==10 and u_max = 1 - 1/(2^32+1), we
310 * have poisson_pdf(lambda,36) < 1-u_max. If instead our
311 * generator uses more bits of mantissa or returns a value
312 * in the range [0,1], then for lambda==10 we need a table
313 * size of 46 instead. For long doubles, the table size
314 * will need to be longer still. */
315#define TABLESIZE46 46
316 double t[TABLESIZE46];
317
318 /* Precompute the table for the u up to and including 0.458.
319 * We will almost certainly need it. */
320 int intlambda = (int)floor (lambda);
321 double P;
322 int tableidx;
323 size_t i = n;
324
325 t[0] = P = exp (-lambda);
326 for (tableidx = 1; tableidx <= intlambda; tableidx++)
327 {
328 P = P*lambda/(double)tableidx;
329 t[tableidx] = t[tableidx-1] + P;
330 }
331
332 while (i-- > 0)
333 {
334 double u = RUNIoct_randu();
335
336 /* If u > 0.458 we know we can jump to floor(lambda) before
337 * comparing (this observation is based on Stadlober's winrand
338 * code). For lambda >= 1, this will be a win. Lambda < 1
339 * is already fast, so adding an extra comparison is not a
340 * problem. */
341 int k = (u > 0.458 ? intlambda : 0);
342
343 /* We aren't using a for loop here because when we find the
344 * right k we want to jump to the next iteration of the
345 * outer loop, and the continue statement will only work for
346 * the inner loop. */
347 nextk:
348 if (u <= t[k])
349 {
350 p[i] = (double) k;
351 continue;
352 }
353 if (++k < tableidx) goto nextk;
354
355 /* We only need high values of the table very rarely so we
356 * don't automatically compute the entire table. */
357 while (tableidx < TABLESIZE46)
358 {
359 P = P*lambda/(double)tableidx;
360 t[tableidx] = t[tableidx-1] + P;
361 /* Make sure we converge to 1.0 just in case u is uniform
362 * on [0,1] rather than [0,1). */
363 if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0;
364 tableidx++;
365 if (u <= t[tableidx-1]) break;
366 }
367
368 /* We are assuming that the table size is big enough here.
369 * This should be true even if RUNI is returning values in
370 * the range [0,1] rather than [0,1). */
371 p[i] = (double)(tableidx-1);
372 }
373}
374
375static void
376poisson_cdf_lookup_float (double lambda, float *p, size_t n)
377{
378 double t[TABLESIZE46];
379
380 /* Precompute the table for the u up to and including 0.458.
381 * We will almost certainly need it. */
382 int intlambda = (int)floor (lambda);
383 double P;
384 int tableidx;
385 size_t i = n;
386
387 t[0] = P = exp (-lambda);
388 for (tableidx = 1; tableidx <= intlambda; tableidx++)
4
Assuming 'tableidx' is > 'intlambda'
5
Loop condition is false. Execution continues on line 394
389 {
390 P = P*lambda/(double)tableidx;
391 t[tableidx] = t[tableidx-1] + P;
392 }
393
394 while (i-- > 0)
6
Loop condition is true. Entering loop body
17
Loop condition is true. Entering loop body
395 {
396 double u = RUNIoct_randu();
397 int k = (u > 0.458 ? intlambda : 0);
7
'?' condition is true
18
'?' condition is true
398 nextk:
399 if (u <= t[k])
8
Taking false branch
11
Taking false branch
19
The right operand of '<=' is a garbage value
400 {
401 p[i] = (float) k;
402 continue;
403 }
404 if (++k < tableidx) goto nextk;
9
Taking true branch
10
Control jumps to line 399
12
Taking false branch
405
406 while (tableidx < TABLESIZE46)
13
Loop condition is true. Entering loop body
407 {
408 P = P*lambda/(double)tableidx;
409 t[tableidx] = t[tableidx-1] + P;
410 if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0;
14
Taking false branch
411 tableidx++;
412 if (u <= t[tableidx-1]) break;
15
Taking true branch
16
Execution continues on line 415
413 }
414
415 p[i] = (float)(tableidx-1);
416 }
417}
418
419/* From Press, et al., Numerical Recipes */
420static void
421poisson_rejection (double lambda, double *p, size_t n)
422{
423 double sq = sqrt (2.0*lambda);
424 double alxm = log (lambda);
425 double g = lambda*alxm - LGAMMAxlgamma(lambda+1.0);
426 size_t i;
427
428 for (i = 0; i < n; i++)
429 {
430 double y, em, t;
431 do
432 {
433 do
434 {
435 y = tan (M_PI3.14159265358979323846*RUNIoct_randu());
436 em = sq * y + lambda;
437 } while (em < 0.0);
438 em = floor (em);
439 t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g);
440 } while (RUNIoct_randu() > t);
441 p[i] = em;
442 }
443}
444
445/* From Press, et al., Numerical Recipes */
446static void
447poisson_rejection_float (double lambda, float *p, size_t n)
448{
449 double sq = sqrt (2.0*lambda);
450 double alxm = log (lambda);
451 double g = lambda*alxm - LGAMMAxlgamma(lambda+1.0);
452 size_t i;
453
454 for (i = 0; i < n; i++)
455 {
456 double y, em, t;
457 do
458 {
459 do
460 {
461 y = tan (M_PI3.14159265358979323846*RUNIoct_randu());
462 em = sq * y + lambda;
463 } while (em < 0.0);
464 em = floor (em);
465 t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g);
466 } while (RUNIoct_randu() > t);
467 p[i] = em;
468 }
469}
470
471/* The cutoff of L <= 1e8 in the following two functions before using
472 * the normal approximation is based on:
473 * > L=1e8; x=floor(linspace(0,2*L,1000));
474 * > max(abs(normal_pdf(x,L,L)-poisson_pdf(x,L)))
475 * ans = 1.1376e-28
476 * For L=1e7, the max is around 1e-9, which is within the step size of RUNI.
477 * For L>1e10 the pprsc function breaks down, as I saw from the histogram
478 * of a large sample, so 1e8 is both small enough and large enough. */
479
480/* Generate a set of poisson numbers with the same distribution */
481void
482oct_fill_randp (double L, octave_idx_type n, double *p)
483{
484 octave_idx_type i;
485 if (L < 0.0 || INFINITE(L)(sizeof (L) == sizeof (float) ? __lo_ieee_float_isinf (L) : __lo_ieee_isinf
(L))
)
486 {
487 for (i=0; i<n; i++)
488 p[i] = NANoctave_NaN;
489 }
490 else if (L <= 10.0)
491 {
492 poisson_cdf_lookup (L, p, n);
493 }
494 else if (L <= 1e8)
495 {
496 for (i=0; i<n; i++)
497 p[i] = pprsc (L);
498 }
499 else
500 {
501 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
502 const double sqrtL = sqrt (L);
503 for (i = 0; i < n; i++)
504 {
505 p[i] = floor (RNORoct_randn()*sqrtL + L + 0.5);
506 if (p[i] < 0.0)
507 p[i] = 0.0; /* will probably never happen */
508 }
509 }
510}
511
512/* Generate one poisson variate */
513double
514oct_randp (double L)
515{
516 double ret;
517 if (L < 0.0) ret = NANoctave_NaN;
518 else if (L <= 12.0)
519 {
520 /* From Press, et al. Numerical recipes */
521 double g = exp (-L);
522 int em = -1;
523 double t = 1.0;
524 do
525 {
526 ++em;
527 t *= RUNIoct_randu();
528 } while (t > g);
529 ret = em;
530 }
531 else if (L <= 1e8)
532 {
533 /* numerical recipes */
534 poisson_rejection (L, &ret, 1);
535 }
536 else if (INFINITE(L)(sizeof (L) == sizeof (float) ? __lo_ieee_float_isinf (L) : __lo_ieee_isinf
(L))
)
537 {
538 /* FIXME R uses NaN, but the normal approx. suggests that as
539 * limit should be inf. Which is correct? */
540 ret = NANoctave_NaN;
541 }
542 else
543 {
544 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
545 ret = floor (RNORoct_randn()*sqrt (L) + L + 0.5);
546 if (ret < 0.0) ret = 0.0; /* will probably never happen */
547 }
548 return ret;
549}
550
551/* Generate a set of poisson numbers with the same distribution */
552void
553oct_fill_float_randp (float FL, octave_idx_type n, float *p)
554{
555 double L = FL;
556 octave_idx_type i;
557 if (L < 0.0 || INFINITE(L)(sizeof (L) == sizeof (float) ? __lo_ieee_float_isinf (L) : __lo_ieee_isinf
(L))
)
1
Taking false branch
558 {
559 for (i=0; i<n; i++)
560 p[i] = NANoctave_NaN;
561 }
562 else if (L <= 10.0)
2
Taking true branch
563 {
564 poisson_cdf_lookup_float (L, p, n);
3
Calling 'poisson_cdf_lookup_float'
565 }
566 else if (L <= 1e8)
567 {
568 for (i=0; i<n; i++)
569 p[i] = pprsc (L);
570 }
571 else
572 {
573 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
574 const double sqrtL = sqrt (L);
575 for (i = 0; i < n; i++)
576 {
577 p[i] = floor (RNORoct_randn()*sqrtL + L + 0.5);
578 if (p[i] < 0.0)
579 p[i] = 0.0; /* will probably never happen */
580 }
581 }
582}
583
584/* Generate one poisson variate */
585float
586oct_float_randp (float FL)
587{
588 double L = FL;
589 float ret;
590 if (L < 0.0) ret = NANoctave_NaN;
591 else if (L <= 12.0)
592 {
593 /* From Press, et al. Numerical recipes */
594 double g = exp (-L);
595 int em = -1;
596 double t = 1.0;
597 do
598 {
599 ++em;
600 t *= RUNIoct_randu();
601 } while (t > g);
602 ret = em;
603 }
604 else if (L <= 1e8)
605 {
606 /* numerical recipes */
607 poisson_rejection_float (L, &ret, 1);
608 }
609 else if (INFINITE(L)(sizeof (L) == sizeof (float) ? __lo_ieee_float_isinf (L) : __lo_ieee_isinf
(L))
)
610 {
611 /* FIXME R uses NaN, but the normal approx. suggests that as
612 * limit should be inf. Which is correct? */
613 ret = NANoctave_NaN;
614 }
615 else
616 {
617 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
618 ret = floor (RNORoct_randn()*sqrt (L) + L + 0.5);
619 if (ret < 0.0) ret = 0.0; /* will probably never happen */
620 }
621 return ret;
622}