Skip to content

[resubmitted] Fix poisson deviate for expected values ~> LONG_MAX #34

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 3 commits into from
Closed

[resubmitted] Fix poisson deviate for expected values ~> LONG_MAX #34

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6449c2d
Fix poisson deviate for expected values ~>LONG_MAX
olsonse Dec 24, 2010
27b30de
add error check to poisson(...) for lam > LONG_MAX
olsonse Jan 10, 2011
9836e26
regenerate mtrand.c with cython 0.12.1 (same version as last generation)
olsonse Jan 10, 2011
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 33 additions & 5 deletions numpy/random/mtrand/distributions.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@
#include <math.h>
#include "distributions.h"
#include <stdio.h>
#include <limits.h>

#ifndef min
#define min(x,y) ((x<y)?x:y)
Expand Down Expand Up @@ -519,17 +520,44 @@ long rk_poisson_ptrs(rk_state *state, double lam)

}

/* lambda max represents the highest (nearly) complete poissonian distribution
* that can be (nearly) fully represented by the bits in the return value.
* The 10*sqrt(LONG_MAX) provides us with 10*sigma upwards width before we hit
* the rollover of long integers.
*
* Note that sqrt(LONG_MAX) can also be written
* sqrt( 1 << (sizeof(long)*8 - 1 ) - 1 )
* (^ -1 here is because longs are signed).
* To avoid non-constant expressions in the initializer of lambda_max (thus to
* conform with ansi c), we will approximate this as
* 0.5* ( sqrt( 1 << (sizeof(long)*8) ) + sqrt( 1 << (sizeof(long)*8-2) ) )
* which can be written as
* 0.5* ( ( 1 << (sizeof(long)*8/2) ) + ( 1 << ((sizeof(long)*8-2)/2) ) )
* or rather
* 0.5* ( ( 1 << (sizeof(long)*4) ) + ( 1 << (sizeof(long)*4-1) ) )
* This is a little larger than sqrt(LONG_MAX), but not by too much.
*/
const double lambda_max
= (double)LONG_MAX
- 10.0 * 0.5 * ( (double)( 1uL << ( sizeof(long)*4uL ) )
+ (double)( 1uL << ( sizeof(long)*4uL - 1ul ) ) );

long rk_poisson(rk_state *state, double lam)
{
if (lam >= 10)
if (lam == 0)
{
return rk_poisson_ptrs(state, lam);
return 0;
}
else if (lam == 0)
else if (lam >= lambda_max)
{
return 0;
/* we return the expected value when we run out of bits */
return (long)lam;
}
else
else if (lam >= 10)
{
return rk_poisson_ptrs(state, lam);
}
else
{
return rk_poisson_mult(state, lam);
}
Expand Down
Loading