34 #if defined (HAVE_CONFIG_H) 48 #define INFINITE lo_ieee_isinf 49 #define RUNI oct_randu() 50 #define RNOR oct_randn() 51 #define LGAMMA xlgamma 65 #define C0 9.18938533204672742e-01 66 #define C1 8.33333333333333333e-02 67 #define C3 -2.77777777777777778e-03 68 #define C5 7.93650793650793651e-04 69 #define C7 -5.95238095238095238e-04 71 static double logfak[30L] =
73 0.00000000000000000, 0.00000000000000000, 0.69314718055994531,
74 1.79175946922805500, 3.17805383034794562, 4.78749174278204599,
75 6.57925121201010100, 8.52516136106541430, 10.60460290274525023,
76 12.80182748008146961, 15.10441257307551530, 17.50230784587388584,
77 19.98721449566188615, 22.55216385312342289, 25.19122118273868150,
78 27.89927138384089157, 30.67186010608067280, 33.50507345013688888,
79 36.39544520803305358, 39.33988418719949404, 42.33561646075348503,
80 45.38013889847690803, 48.47118135183522388, 51.60667556776437357,
81 54.78472939811231919, 58.00360522298051994, 61.26170176100200198,
82 64.55753862700633106, 67.88974313718153498, 71.25703896716800901
92 + r*(
C1 + rr*(
C3 + rr*(
C5 + rr*
C7))));
95 return (logfak[static_cast<int> (
k)]);
130 f (
double k,
double l_nu,
double c_pm)
132 return exp (
k * l_nu -
flogfak (
k) - c_pm);
138 static double my_last = -1.0;
139 static double m, k2, k4, k1, k5;
140 static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
141 f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
149 Ds = std::sqrt (my + 0.25);
154 k2 =
ceil (my - 0.5 - Ds);
156 k1 = k2 + k2 - m + 1L;
166 r4 = my / (k4 + 1.0);
167 r5 = my / (k5 + 1.0);
178 f2 =
f (k2, l_my, c_pm);
179 f4 =
f (k4, l_my, c_pm);
180 f1 =
f (k1, l_my, c_pm);
181 f5 =
f (k5, l_my, c_pm);
185 p1 = f2 * (dl + 1.0);
187 p3 = f4 * (dr + 1.0) + p2;
197 if ((U =
RUNI * p6) < p2)
202 if ((
V = U - p1) < 0.0)
return (k2 +
std::floor (U/f2));
205 if ((W =
V / dl) < f1 )
return (k1 +
std::floor (
V/f1));
210 if (W <= f2 - Dk * (f2 - f2/r2))
214 if ((
V = f2 + f2 - W) < 1.0)
217 if (
V <= f2 + Dk * (1.0 - f2)/(dl + 1.0))
221 if (
V <=
f (Y, l_my, c_pm))
return (Y);
229 if ((
V = U - p3) < 0.0)
return (k4 -
std::floor ((U - p2)/f4));
232 if ((W =
V / dr) < f5 )
return (k5 -
std::floor (
V/f5));
237 if (W <= f4 - Dk * (f4 - f4*r4))
241 if ((
V = f4 + f4 - W) < 1.0)
244 if (
V <= f4 + Dk * (1.0 - f4)/ dr)
248 if (
V <=
f (Y, l_my, c_pm))
return (Y);
258 if ((X = k1 - Dk) < 0L)
continue;
260 if (W <= f1 - Dk * (f1 - f1/r1))
268 if (W <= f5 - Dk * (f5 - f5*r5))
301 int intlambda =
static_cast<int> (
std::floor (lambda));
306 t[0] = P = exp (-lambda);
307 for (tableidx = 1; tableidx <= intlambda; tableidx++)
309 P = P*lambda/
static_cast<double> (tableidx);
310 t[tableidx] =
t[tableidx-1] + P;
322 int k = (
u > 0.458 ? intlambda : 0);
331 p[
i] =
static_cast<double> (
k);
341 P = P*lambda/
static_cast<double> (tableidx);
342 t[tableidx] =
t[tableidx-1] + P;
345 if (
t[tableidx] ==
t[tableidx-1])
t[tableidx] = 1.0;
347 if (
u <=
t[tableidx-1])
break;
353 p[
i] =
static_cast<double> (tableidx-1);
364 int intlambda =
static_cast<int> (
std::floor (lambda));
369 t[0] = P = exp (-lambda);
370 for (tableidx = 1; tableidx <= intlambda; tableidx++)
372 P = P*lambda/
static_cast<double> (tableidx);
373 t[tableidx] =
t[tableidx-1] + P;
379 int k = (
u > 0.458 ? intlambda : 0);
383 p[
i] =
static_cast<float> (
k);
391 P = P*lambda/
static_cast<double> (tableidx);
392 t[tableidx] =
t[tableidx-1] + P;
393 if (
t[tableidx] ==
t[tableidx-1])
t[tableidx] = 1.0;
395 if (
u <=
t[tableidx-1])
break;
398 p[
i] =
static_cast<float> (tableidx-1);
406 double sq = std::sqrt (2.0*lambda);
408 double g = lambda*alxm -
LGAMMA(lambda+1.0);
411 for (
i = 0;
i < n;
i++)
419 em = sq *
y + lambda;
422 t = 0.9*(1.0+
y*
y)*exp (em*alxm-
flogfak (em)-g);
432 double sq = std::sqrt (2.0*lambda);
434 double g = lambda*alxm -
LGAMMA(lambda+1.0);
437 for (
i = 0;
i < n;
i++)
445 em = sq *
y + lambda;
448 t = 0.9*(1.0+
y*
y)*exp (em*alxm-
flogfak (em)-g);
485 const double sqrtL = std::sqrt (L);
486 for (
i = 0;
i < n;
i++)
529 if (ret < 0.0) ret = 0.0;
557 const double sqrtL = std::sqrt (L);
558 for (
i = 0;
i < n;
i++)
602 if (ret < 0.0) ret = 0.0;
F77_RET_T const F77_INT const F77_INT const F77_INT const F77_DBLE const F77_DBLE F77_INT F77_DBLE * V
static void poisson_cdf_lookup(double lambda, double *p, size_t n)
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the base of natural logarithms The constant ex $e satisfies the equation log(e)
std::complex< T > ceil(const std::complex< T > &x)
std::complex< T > floor(const std::complex< T > &x)
static double flogfak(double k)
static double xlgamma(double x)
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function t
static void poisson_rejection_float(double lambda, float *p, size_t n)
double oct_randp(double L)
static double pprsc(double my)
static void poisson_rejection(double lambda, double *p, size_t n)
static void poisson_cdf_lookup_float(double lambda, float *p, size_t n)
void oct_fill_float_randp(float FL, octave_idx_type n, float *p)
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol NaN(Not a Number). NaN is the result of operations which do not produce a well defined 0 result. Common operations which produce a NaN are arithmetic with infinity ex($\infty - \infty$)
the element is set to zero In other the statement xample y
static double f(double k, double l_nu, double c_pm)
float oct_float_randp(float FL)
F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE * x
void oct_fill_randp(double L, octave_idx_type n, double *p)