in commons-statistics-inference/src/main/java/org/apache/commons/statistics/inference/KolmogorovSmirnovDistribution.java [260:321]
private static RealSquareMatrix createH(double x, int n) {
// MATH-437:
// This is *not* (int) (n * x) + 1.
// This is only ever called when 1/n < x < 1 - 1/n.
// => h cannot be >= 1 when using ceil. h can be 0 if nx is integral.
final int k = (int) Math.ceil(n * x);
final double h = k - n * x;
final int m = 2 * k - 1;
final double[] data = new double[m * m];
// Start by filling everything with either 0 or 1.
for (int i = 0; i < m; ++i) {
// h[i][j] = i - j + 1 < 0 ? 0 : 1
// => h[i][j<=i+1] = 1
final int jend = Math.min(m - 1, i + 1);
for (int j = i * m; j <= i * m + jend; j++) {
data[j] = 1;
}
}
// Setting up power-array to avoid calculating the same value twice:
// hp[0] = h^1, ..., hp[m-1] = h^m
final double[] hp = new double[m];
hp[0] = h;
for (int i = 1; i < m; ++i) {
// Avoid compound rounding errors using h * hp[i - 1]
// with Math.pow as it is within 1 ulp of the exact result
hp[i] = Math.pow(h, i + 1);
}
// First column and last row has special values (each other reversed).
for (int i = 0; i < m; ++i) {
data[i * m] -= hp[i];
data[(m - 1) * m + i] -= hp[m - i - 1];
}
// [1] states: "For 1/2 < h < 1 the bottom left element of the matrix should be
// (1 - 2*h^m + (2h - 1)^m )/m!"
if (2 * h - 1 > 0) {
data[(m - 1) * m] += Math.pow(2 * h - 1, m);
}
// Aside from the first column and last row, the (i, j)-th element is 1/(i - j + 1)! if i -
// j + 1 >= 0, else 0. 1's and 0's are already put, so only division with (i - j + 1)! is
// needed in the elements that have 1's. Note that i - j + 1 > 0 <=> i + 1 > j instead of
// j'ing all the way to m. Also note that we can use pre-computed factorials given
// the limits where this method is called.
for (int i = 0; i < m; ++i) {
final int im = i * m;
for (int j = 0; j < i + 1; ++j) {
// Here (i - j + 1 > 0)
// Divide by (i - j + 1)!
// Note: This method is used when:
// n <= 140; nxx < 0.754693
// n <= 100000; n x^1.5 < 1.4
// max m ~ 2nx ~ (1.4/1e5)^(2/3) * 2e5 = 116
// Use a tabulated factorial
data[im + j] /= Factorial.doubleValue(i - j + 1);
}
}
return SquareMatrixSupport.create(m, data);
}