in commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratSampler.java [1078:1160]
private double edgeSample(long xx) {
// Expected frequency = 0.0117188
// Drop the sign bit to create u:
long u1 = xx & MAX_INT64;
// Extract the sign bit for use later
// Use 2 - 1 or 0 - 1
final double signBit = ((xx >>> 62) & 0x2) - 1.0;
final int j = selectRegion();
// Four kinds of overhangs:
// j = 0 : Sample from tail
// 0 < j < J_INFLECTION : Overhang is concave; only sample from Lower-Left triangle
// j = J_INFLECTION : Must sample from entire overhang rectangle
// j > J_INFLECTION : Overhangs are convex; implicitly accept point in Lower-Left triangle
//
// Conditional statements are arranged such that the more likely outcomes are first.
double x;
if (j > J_INFLECTION) {
// Convex overhang
// Expected frequency: 0.00892899
// Observed loop repeat frequency: 0.389804
for (;;) {
x = interpolate(X, j, u1);
// u2 = u1 + (u2 - u1) = u1 + uDistance
final long uDistance = randomInt63() - u1;
if (uDistance >= 0) {
// Lower-left triangle
break;
}
if (uDistance >= CONVEX_E_MAX &&
// Within maximum distance of f(x) from the triangle hypotenuse.
// Frequency (per upper-right triangle): 0.431497
// Reject frequency: 0.489630
interpolate(Y, j, u1 + uDistance) < Math.exp(-0.5 * x * x)) {
break;
}
// uDistance < E_MAX (upper-right triangle) or rejected as above the curve
u1 = randomInt63();
}
} else if (j < J_INFLECTION) {
if (j == 0) {
// Tail
// Expected frequency: 0.000276902
// Note: Although less frequent than the next branch, j == 0 is a subset of
// j < J_INFLECTION and must be first.
// Observed loop repeat frequency: 0.0634786
do {
x = ONE_OVER_X_0 * exponential.sample();
} while (exponential.sample() < 0.5 * x * x);
x += X_0;
} else {
// Concave overhang
// Expected frequency: 0.00249694
// Observed loop repeat frequency: 0.0123784
for (;;) {
// Create a second uniform deviate (as u1 is recycled).
final long u = randomInt63();
// If u2 < u1 then reflect in the hypotenuse by swapping u1 and u2.
// Use conditional ternary to avoid a 50/50 branch statement to swap the pair.
final long u2 = u1 < u ? u : u1;
u1 = u1 < u ? u1 : u;
x = interpolate(X, j, u1);
if (u2 - u1 > CONCAVE_E_MAX ||
interpolate(Y, j, u2) < Math.exp(-0.5 * x * x)) {
break;
}
u1 = randomInt63();
}
}
} else {
// Inflection point
// Expected frequency: 0.000015914
// Observed loop repeat frequency: 0.500213
for (;;) {
x = interpolate(X, j, u1);
if (interpolate(Y, j, randomInt63()) < Math.exp(-0.5 * x * x)) {
break;
}
u1 = randomInt63();
}
}
return signBit * x;
}