34 #if defined (HAVE_CONFIG_H)
47 #define INFINITE lo_ieee_isinf
48 #define RUNI oct_randu()
49 #define RNOR oct_randn()
50 #define LGAMMA xlgamma
56 #if defined (HAVE_LGAMMA)
77 #define C0 9.18938533204672742e-01
78 #define C1 8.33333333333333333e-02
79 #define C3 -2.77777777777777778e-03
80 #define C5 7.93650793650793651e-04
81 #define C7 -5.95238095238095238e-04
83 static double logfak[30L] =
85 0.00000000000000000, 0.00000000000000000, 0.69314718055994531,
86 1.79175946922805500, 3.17805383034794562, 4.78749174278204599,
87 6.57925121201010100, 8.52516136106541430, 10.60460290274525023,
88 12.80182748008146961, 15.10441257307551530, 17.50230784587388584,
89 19.98721449566188615, 22.55216385312342289, 25.19122118273868150,
90 27.89927138384089157, 30.67186010608067280, 33.50507345013688888,
91 36.39544520803305358, 39.33988418719949404, 42.33561646075348503,
92 45.38013889847690803, 48.47118135183522388, 51.60667556776437357,
93 54.78472939811231919, 58.00360522298051994, 61.26170176100200198,
94 64.55753862700633106, 67.88974313718153498, 71.25703896716800901
106 return (logfak[static_cast<int> (k)]);
141 f (
double k,
double l_nu,
double c_pm)
143 return exp (k * l_nu -
flogfak (k) - c_pm);
149 static double my_last = -1.0;
150 static double m, k2, k4, k1, k5;
151 static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
152 f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
160 Ds = sqrt (my + 0.25);
165 k2 =
ceil (my - 0.5 - Ds);
167 k1 = k2 + k2 - m + 1L;
177 r4 = my / (k4 + 1.0);
178 r5 = my / (k5 + 1.0);
189 f2 =
f (k2, l_my, c_pm);
190 f4 =
f (k4, l_my, c_pm);
191 f1 =
f (k1, l_my, c_pm);
192 f5 =
f (k5, l_my, c_pm);
196 p1 = f2 * (dl + 1.0);
198 p3 = f4 * (dr + 1.0) + p2;
208 if ((U =
RUNI * p6) < p2)
213 if ((V = U - p1) < 0.0)
return (k2 +
std::floor (U/f2));
216 if ((W = V / dl) < f1 )
return (k1 +
std::floor (V/f1));
221 if (W <= f2 - Dk * (f2 - f2/r2))
225 if ((V = f2 + f2 - W) < 1.0)
228 if (V <= f2 + Dk * (1.0 - f2)/(dl + 1.0))
232 if (V <=
f (Y, l_my, c_pm))
return (Y);
240 if ((V = U - p3) < 0.0)
return (k4 -
std::floor ((U - p2)/f4));
243 if ((W = V / dr) < f5 )
return (k5 -
std::floor (V/f5));
248 if (W <= f4 - Dk * (f4 - f4*r4))
252 if ((V = f4 + f4 - W) < 1.0)
255 if (V <= f4 + Dk * (1.0 - f4)/ dr)
259 if (V <=
f (Y, l_my, c_pm))
return (Y);
269 if ((X = k1 - Dk) < 0L)
continue;
271 if (W <= f1 - Dk * (f1 - f1/r1))
279 if (W <= f5 - Dk * (f5 - f5*r5))
312 int intlambda =
static_cast<int> (
std::floor (lambda));
317 t[0] = P = exp (-lambda);
318 for (tableidx = 1; tableidx <= intlambda; tableidx++)
320 P = P*lambda/
static_cast<double> (tableidx);
321 t[tableidx] = t[tableidx-1] +
P;
333 int k = (u > 0.458 ? intlambda : 0);
342 p[
i] =
static_cast<double> (
k);
352 P = P*lambda/
static_cast<double> (tableidx);
353 t[tableidx] = t[tableidx-1] +
P;
356 if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0;
358 if (u <= t[tableidx-1])
break;
364 p[
i] =
static_cast<double> (tableidx-1);
375 int intlambda =
static_cast<int> (
std::floor (lambda));
380 t[0] = P = exp (-lambda);
381 for (tableidx = 1; tableidx <= intlambda; tableidx++)
383 P = P*lambda/
static_cast<double> (tableidx);
384 t[tableidx] = t[tableidx-1] +
P;
390 int k = (u > 0.458 ? intlambda : 0);
394 p[
i] =
static_cast<float> (
k);
402 P = P*lambda/
static_cast<double> (tableidx);
403 t[tableidx] = t[tableidx-1] +
P;
404 if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0;
406 if (u <= t[tableidx-1])
break;
409 p[
i] =
static_cast<float> (tableidx-1);
417 double sq = sqrt (2.0*lambda);
419 double g = lambda*alxm -
LGAMMA(lambda+1.0);
422 for (i = 0; i < n; i++)
430 em = sq * y + lambda;
433 t = 0.9*(1.0+y*
y)*exp (em*alxm-
flogfak (em)-g);
443 double sq = sqrt (2.0*lambda);
445 double g = lambda*alxm -
LGAMMA(lambda+1.0);
448 for (i = 0; i < n; i++)
456 em = sq * y + lambda;
459 t = 0.9*(1.0+y*
y)*exp (em*alxm-
flogfak (em)-g);
496 const double sqrtL = sqrt (L);
497 for (i = 0; i < n; i++)
540 if (ret < 0.0) ret = 0.0;
568 const double sqrtL = sqrt (L);
569 for (i = 0; i < n; i++)
613 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
subroutine dlgams(X, DLGAM, SGNGAM)
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)
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)
#define F77_XFCN(f, F, args)
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup calculate Y_a and Y _d item Given Y
double oct_randp(double L)
nd deftypefn *octave_map m
static double pprsc(double my)
OCTAVE_EXPORT octave_value_list W
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)
With real return the complex result
=val(i)}if ode{val(i)}occurs in table i
the element is set to zero In other the statement xample y
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 $)
static double f(double k, double l_nu, double c_pm)
float oct_float_randp(float FL)
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &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 F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
void oct_fill_randp(double L, octave_idx_type n, double *p)