in commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/MarsagliaTsangWangDiscreteSampler.java [816:880]
private static SharedStateDiscreteSampler createPoissonDistributionFromXMode(
UniformRandomProvider rng, double mean) {
// If mean >= 21.4, generate from largest p-value up, then largest down.
// The largest p-value will be at the mode (floor(mean)).
// Find p(x=mode)
final int mode = (int) mean;
// This transform is stable until mean >= 1941 where p will result in Infinity
// before the divisor i is large enough to start reducing the product (i.e. i > c).
final double c = mean * Math.exp(-mean / mode);
double p = 1.0;
for (int i = 1; i <= mode; i++) {
p *= c / i;
}
final double pMode = p;
// Find the upper limit using recursive computation of the p-value.
// Note this will exit when i overflows to negative so no check on the range
int i = mode + 1;
while (p * DOUBLE_31 >= 1) {
p *= mean / i++;
}
final int last = i - 2;
// Find the lower limit using recursive computation of the p-value.
p = pMode;
int j = -1;
for (i = mode - 1; i >= 0; i--) {
p *= (i + 1) / mean;
if (p * DOUBLE_31 < 1) {
j = i;
break;
}
}
// Probabilities are 30-bit integers, assumed denominator 2^30.
// This is the minimum sample value: prob[x - offset] = p(x)
final int offset = j + 1;
final int size = last - offset + 1;
final int[] prob = new int[size];
p = pMode;
prob[mode - offset] = toUnsignedInt30(p);
// The sum must exceed 2^30. In edges cases this is false due to round-off.
int sum = prob[mode - offset];
// From mode to upper limit
for (i = mode + 1; i <= last; i++) {
p *= mean / i;
prob[i - offset] = toUnsignedInt30(p);
sum += prob[i - offset];
}
// From mode to lower limit
p = pMode;
for (i = mode - 1; i >= offset; i--) {
p *= (i + 1) / mean;
prob[i - offset] = toUnsignedInt30(p);
sum += prob[i - offset];
}
// If the sum is < 2^30 add the remaining sum to the mode.
// If above 2^30 then the effect is truncation of the long tail of the distribution.
prob[mode - offset] += Math.max(0, INT_30 - sum);
return createSampler(rng, POISSON_NAME, prob, offset);
}