in commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/MarsagliaTsangWangDiscreteSampler.java [1049:1095]
private static SharedStateDiscreteSampler createBinomialDistributionSampler(
UniformRandomProvider rng, int trials, double probabilityOfSuccess) {
// The maximum supported value for Math.exp is approximately -744.
// This occurs when trials is large and p is close to 1.
// Handle this by using an inversion: generate j=Binomial(n,1-p), return n-j
final boolean useInversion = probabilityOfSuccess > 0.5;
final double p = useInversion ? 1 - probabilityOfSuccess : probabilityOfSuccess;
// Check if the distribution can be computed
final double p0 = Math.exp(trials * Math.log(1 - p));
if (p0 < Double.MIN_VALUE) {
throw new IllegalArgumentException("Unable to compute distribution");
}
// First find size of probability array
double t = p0;
final double h = p / (1 - p);
// Find first probability above the threshold of 2^-31
int begin = 0;
if (t * DOUBLE_31 < 1) {
// Somewhere after p(0)
// Note:
// If this loop is entered p(0) is < 2^-31.
// This has been tested at the extreme for p(0)=Double.MIN_VALUE and either
// p=0.5 or trials=2^16-1 and does not fail to find the beginning.
for (int i = 1; i <= trials; i++) {
t *= (trials + 1 - i) * h / i;
if (t * DOUBLE_31 >= 1) {
begin = i;
break;
}
}
}
// Find last probability
int end = trials;
for (int i = begin + 1; i <= trials; i++) {
t *= (trials + 1 - i) * h / i;
if (t * DOUBLE_31 < 1) {
end = i - 1;
break;
}
}
return createBinomialDistributionSamplerFromRange(rng, trials, p, useInversion,
p0, begin, end);
}