in commons-statistics-inference/src/main/java/org/apache/commons/statistics/inference/KolmogorovSmirnovDistribution.java [333:438]
private static double pomeranz(double x, int n) {
final double t = n * x;
// Store floor(A-t) and ceil(A+t). This does not require computing A.
final int[] amt = new int[2 * n + 3];
final int[] apt = new int[2 * n + 3];
computeA(n, t, amt, apt);
// Precompute ((A[i] - A[i-1])/n)^(j-k) / (j-k)!
// A[i] - A[i-1] has 4 possible values (based on multiples of A2)
// A1 - A0 = 0 - 0 = 0
// A2 - A1 = A2 - 0 = A2
// A3 - A2 = (1 - A2) - A2 = 1 - 2 * A2
// A4 - A3 = (A2 + 1) - (1 - A2) = 2 * A2
// A5 - A4 = (1 - A2 + 1) - (A2 + 1) = 1 - 2 * A2
// A6 - A5 = (A2 + 1 + 1) - (1 - A2 + 1) = 2 * A2
// A7 - A6 = (1 - A2 + 1 + 1) - (A2 + 1 + 1) = 1 - 2 * A2
// A8 - A7 = (A2 + 1 + 1 + 1) - (1 - A2 + 1 + 1) = 2 * A2
// ...
// Ai - Ai-1 = ((i-1)/2 - A2) - (A2 + (i-2)/2) = 1 - 2 * A2 ; i = odd
// Ai - Ai-1 = (A2 + (i-1)/2) - ((i-2)/2 - A2) = 2 * A2 ; i = even
// ...
// A2n+2 - A2n+1 = n - (n - A2) = A2
// ap[][j - k] = ((A[i] - A[i-1])/n)^(j-k) / (j-k)!
// for each case: A[i] - A[i-1] in [A2, 1 - 2 * A2, 2 * A2]
// Ignore case 0 as this is not used. Factors are ap[0] = 1, else 0.
// If A2==0.5 then this is computed as a no-op due to multiplication by zero.
final int n2 = n + 2;
final double[][] ap = new double[3][n2];
final double a2 = Math.min(t - Math.floor(t), Math.ceil(t) - t);
computeAP(ap[0], a2 / n);
computeAP(ap[1], (1 - 2 * a2) / n);
computeAP(ap[2], (2 * a2) / n);
// Current and previous V
double[] vc = new double[n2];
double[] vp = new double[n2];
// Count of re-scaling
int scale = 0;
// V_1,1 = 1
vc[1] = 1;
for (int i = 2; i <= 2 * n + 2; i++) {
final double[] v = vc;
vc = vp;
vp = v;
// This is useful for following current values of vc
Arrays.fill(vc, 0);
// Select (A[i] - A[i-1]) factor
final double[] p;
if (i == 2 || i == 2 * n + 2) {
// First or last
p = ap[0];
} else {
// odd: [1] 1 - 2 * 2A
// even: [2] 2 * A2
p = ap[2 - (i & 1)];
}
// Set limits.
// j is the ultimate bound for k and should be in [1, n+1]
final int jmin = Math.max(1, amt[i] + 2);
final int jmax = Math.min(n + 1, apt[i]);
final int k1 = Math.max(1, amt[i - 1] + 2);
// All numbers will reduce in size.
// Maintain the largest close to 1.0.
// This is a change from Simard and L’Ecuyer which scaled based on the smallest.
double max = 0;
for (int j = jmin; j <= jmax; j++) {
final int k2 = Math.min(j, apt[i - 1]);
// Accurate sum.
// vp[high] is smaller
// p[high] is smaller
// Sum ascending has smaller products first.
double sum = 0;
for (int k = k1; k <= k2; k++) {
sum += vp[k] * p[j - k];
}
vc[j] = sum;
if (max < sum) {
// Note: max *may* always be the first sum: vc[jmin]
max = sum;
}
}
// Rescale if too small
if (max < P_DOWN_SCALE) {
// Only scale in current range from V
for (int j = jmin; j <= jmax; j++) {
vc[j] *= P_UP_SCALE;
}
scale -= P_SCALE_POWER;
}
}
// F_n(x) = n! V_{2n+2,n+1}
double v = vc[n + 1];
// This method is used when n < 140 where all n! are finite.
// v is below 1 so we can directly compute the result without using logs.
v *= Factorial.doubleValue(n);
// Return the CDF (rescaling as required)
return Math.scalb(v, scale);
}