in commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LargeMeanPoissonSampler.java [253:316]
public int sample() {
// This will never be null. It may be a no-op delegate that returns zero.
final int y2 = smallMeanPoissonSampler.sample();
double x;
double y;
double v;
int a;
double t;
double qr;
double qa;
while (true) {
// Step 1:
final double u = rng.nextDouble();
if (u <= p1) {
// Step 2:
final double n = gaussian.sample();
x = n * sqrtLambdaPlusHalfDelta - 0.5d;
if (x > delta || x < -lambda) {
continue;
}
y = x < 0 ? Math.floor(x) : Math.ceil(x);
final double e = exponential.sample();
v = -e - 0.5 * n * n + c1;
} else {
// Step 3:
if (u > p1 + p2) {
y = lambda;
break;
}
x = delta + (twolpd / delta) * exponential.sample();
y = Math.ceil(x);
v = -exponential.sample() - delta * (x + 1) / twolpd;
}
// The Squeeze Principle
// Step 4.1:
a = x < 0 ? 1 : 0;
t = y * (y + 1) / (2 * lambda);
// Step 4.2
if (v < -t && a == 0) {
y = lambda + y;
break;
}
// Step 4.3:
qr = t * ((2 * y + 1) / (6 * lambda) - 1);
qa = qr - (t * t) / (3 * (lambda + a * (y + 1)));
// Step 4.4:
if (v < qa) {
y = lambda + y;
break;
}
// Step 4.5:
if (v > qr) {
continue;
}
// Step 4.6:
if (v < y * logLambda - getFactorialLog((int) (y + lambda)) + logLambdaFactorial) {
y = lambda + y;
break;
}
}
return (int) Math.min(y2 + (long) y, Integer.MAX_VALUE);
}