org.apache.commons.math3.util.CombinatoricsUtils Java Examples

The following examples show how to use org.apache.commons.math3.util.CombinatoricsUtils. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example #1
Source File: DSCompiler.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
/** Evaluate Taylor expansion of a derivative structure.
 * @param ds array holding the derivative structure
 * @param dsOffset offset of the derivative structure in its array
 * @param delta parameters offsets (Δx, Δy, ...)
 * @return value of the Taylor expansion at x + Δx, y + Δy, ...
 * @throws MathArithmeticException if factorials becomes too large
 */
public double taylor(final double[] ds, final int dsOffset, final double ... delta)
   throws MathArithmeticException {
    double value = 0;
    for (int i = getSize() - 1; i >= 0; --i) {
        final int[] orders = getPartialDerivativeOrders(i);
        double term = ds[dsOffset + i];
        for (int k = 0; k < orders.length; ++k) {
            if (orders[k] > 0) {
                try {
                    term *= FastMath.pow(delta[k], orders[k]) /
                    CombinatoricsUtils.factorial(orders[k]);
                } catch (NotPositiveException e) {
                    // this cannot happen
                    throw new MathInternalError(e);
                }
            }
        }
        value += term;
    }
    return value;
}
 
Example #2
Source File: Binominal.java    From rapidminer-studio with GNU Affero General Public License v3.0 6 votes vote down vote up
@Override
protected double compute(double value1, double value2) {

	// special case for handling missing values
	if (Double.isNaN(value1) || Double.isNaN(value2) || Double.isInfinite(value1) || Double.isInfinite(value2)) {
		return Double.NaN;
	}

	int v1 = (int) value1;
	int v2 = (int) value2;

	if (v1 < 0 || v2 < 0) {
		throw new FunctionInputException("expression_parser.function_non_negative", getFunctionName());
	}
	// This is the common definition for the case for k > n.
	if (v2 > v1) {
		return 0;
	} else {
		return CombinatoricsUtils.binomialCoefficientDouble(v1, v2);
	}
}
 
Example #3
Source File: LocalitySensitiveHashTest.java    From oryx with Apache License 2.0 6 votes vote down vote up
private static void doTestHashesBits(double sampleRate, int numCores, int numHashes, int maxBitsDiffering) {
  LocalitySensitiveHash lsh = new LocalitySensitiveHash(sampleRate, 10, numCores);
  assertEquals(numHashes, lsh.getNumHashes());
  assertEquals(1L << numHashes, lsh.getNumPartitions());
  assertEquals(maxBitsDiffering, lsh.getMaxBitsDiffering());
  if (sampleRate >= 1.0) {
    assertEquals(lsh.getMaxBitsDiffering(), lsh.getNumHashes());
  }
  long partitionsToTry = 0;
  for (int i = 0; i <= maxBitsDiffering; i++) {
    partitionsToTry += CombinatoricsUtils.binomialCoefficient(numHashes, i);
  }
  if (numHashes < LocalitySensitiveHash.MAX_HASHES) {
    assertLessOrEqual((double) partitionsToTry / (1 << numHashes), sampleRate);
  }
}
 
Example #4
Source File: LocalitySensitiveHash.java    From oryx with Apache License 2.0 6 votes vote down vote up
/**
 * @param vector vector whose dot product with hashed vectors is to be maximized
 * @return indices of partitions containing candidates to check
 */
int[] getCandidateIndices(float[] vector) {
  int mainIndex = getIndexFor(vector);
  // Simple cases
  int numHashes = getNumHashes();
  if (numHashes == maxBitsDiffering) {
    return allIndices;
  }
  if (maxBitsDiffering == 0) {
    return new int[] { mainIndex };
  }
  // Other cases
  int howMany = 0;
  for (int i = 0; i <= maxBitsDiffering; i++) {
    howMany += (int) CombinatoricsUtils.binomialCoefficient(numHashes, i);
  }
  int[] result = new int[howMany];
  System.arraycopy(candidateIndicesPrototype, 0, result, 0, howMany);
  for (int i = 0; i < howMany; i++) {
    result[i] ^= mainIndex;
  }
  return result;
}
 
Example #5
Source File: GaussLaguerreWeightAndAbscissaFunction.java    From Strata with Apache License 2.0 6 votes vote down vote up
/**
 * {@inheritDoc}
 */
@Override
public GaussianQuadratureData generate(int n) {
  ArgChecker.isTrue(n > 0);
  Pair<DoubleFunction1D, DoubleFunction1D>[] polynomials = LAGUERRE.getPolynomialsAndFirstDerivative(n, _alpha);
  Pair<DoubleFunction1D, DoubleFunction1D> pair = polynomials[n];
  DoubleFunction1D p1 = polynomials[n - 1].getFirst();
  DoubleFunction1D function = pair.getFirst();
  DoubleFunction1D derivative = pair.getSecond();
  double[] x = new double[n];
  double[] w = new double[n];
  double root = 0;
  for (int i = 0; i < n; i++) {
    root = ROOT_FINDER.getRoot(function, derivative, getInitialRootGuess(root, i, n, x));
    x[i] = root;
    w[i] =
        -GAMMA_FUNCTION.applyAsDouble(_alpha + n) / CombinatoricsUtils.factorialDouble(n) /
            (derivative.applyAsDouble(root) * p1.applyAsDouble(root));
  }
  return new GaussianQuadratureData(x, w);
}
 
Example #6
Source File: GaussJacobiWeightAndAbscissaFunction.java    From Strata with Apache License 2.0 6 votes vote down vote up
/**
 * {@inheritDoc}
 */
@Override
public GaussianQuadratureData generate(int n) {
  ArgChecker.isTrue(n > 0, "n > 0");
  Pair<DoubleFunction1D, DoubleFunction1D>[] polynomials = JACOBI.getPolynomialsAndFirstDerivative(n, _alpha, _beta);
  Pair<DoubleFunction1D, DoubleFunction1D> pair = polynomials[n];
  DoubleFunction1D previous = polynomials[n - 1].getFirst();
  DoubleFunction1D function = pair.getFirst();
  DoubleFunction1D derivative = pair.getSecond();
  double[] x = new double[n];
  double[] w = new double[n];
  double root = 0;
  for (int i = 0; i < n; i++) {
    double d = 2 * n + _c;
    root = getInitialRootGuess(root, i, n, x);
    root = ROOT_FINDER.getRoot(function, derivative, root);
    x[i] = root;
    w[i] =
        GAMMA_FUNCTION.applyAsDouble(_alpha + n) * GAMMA_FUNCTION.applyAsDouble(_beta + n) /
            CombinatoricsUtils.factorialDouble(n) / GAMMA_FUNCTION.applyAsDouble(n + _c + 1) * d *
            Math.pow(2, _c)
            / (derivative.applyAsDouble(root) * previous.applyAsDouble(root));
  }
  return new GaussianQuadratureData(x, w);
}
 
Example #7
Source File: DSCompiler.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
/** Evaluate Taylor expansion of a derivative structure.
 * @param ds array holding the derivative structure
 * @param dsOffset offset of the derivative structure in its array
 * @param delta parameters offsets (&Delta;x, &Delta;y, ...)
 * @return value of the Taylor expansion at x + &Delta;x, y + &Delta;y, ...
 * @throws MathArithmeticException if factorials becomes too large
 */
public double taylor(final double[] ds, final int dsOffset, final double ... delta)
   throws MathArithmeticException {
    double value = 0;
    for (int i = getSize() - 1; i >= 0; --i) {
        final int[] orders = getPartialDerivativeOrders(i);
        double term = ds[dsOffset + i];
        for (int k = 0; k < orders.length; ++k) {
            if (orders[k] > 0) {
                try {
                    term *= FastMath.pow(delta[k], orders[k]) /
                    CombinatoricsUtils.factorial(orders[k]);
                } catch (NotPositiveException e) {
                    // this cannot happen
                    throw new MathInternalError(e);
                }
            }
        }
        value += term;
    }
    return value;
}
 
Example #8
Source File: GlobalClusteringCoefficientTest.java    From flink with Apache License 2.0 5 votes vote down vote up
@Test
public void testWithCompleteGraph() throws Exception {
	long expectedDegree = completeGraphVertexCount - 1;
	long expectedCount = completeGraphVertexCount * CombinatoricsUtils.binomialCoefficient((int) expectedDegree, 2);

	validate(completeGraph, expectedCount, expectedCount);
}
 
Example #9
Source File: PascalDistribution.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
/** {@inheritDoc} */
@Override
public double logProbability(int x) {
    double ret;
    if (x < 0) {
        ret = Double.NEGATIVE_INFINITY;
    } else {
        ret = CombinatoricsUtils.binomialCoefficientLog(x +
              numberOfSuccesses - 1, numberOfSuccesses - 1) +
              logProbabilityOfSuccess * numberOfSuccesses +
              log1mProbabilityOfSuccess * x;
    }
    return ret;
}
 
Example #10
Source File: KolmogorovSmirnovTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
/**
 * Computes \(P(D_{n,m} > d)\) if {@code strict} is {@code true}; otherwise \(P(D_{n,m} \ge
 * d)\), where \(D_{n,m}\) is the 2-sample Kolmogorov-Smirnov statistic. See
 * {@link #kolmogorovSmirnovStatistic(double[], double[])} for the definition of \(D_{n,m}\).
 * <p>
 * The returned probability is exact, obtained by enumerating all partitions of {@code m + n}
 * into {@code m} and {@code n} sets, computing \(D_{n,m}\) for each partition and counting the
 * number of partitions that yield \(D_{n,m}\) values exceeding (resp. greater than or equal to)
 * {@code d}.
 * </p>
 * <p>
 * <strong>USAGE NOTE</strong>: Since this method enumerates all combinations in \({m+n} \choose
 * {n}\), it is very slow if called for large {@code m, n}. For this reason,
 * {@link #kolmogorovSmirnovTest(double[], double[])} uses this only for {@code m * n < }
 * {@value #SMALL_SAMPLE_PRODUCT}.
 * </p>
 *
 * @param d D-statistic value
 * @param n first sample size
 * @param m second sample size
 * @param strict whether or not the probability to compute is expressed as a strict inequality
 * @return probability that a randomly selected m-n partition of m + n generates \(D_{n,m}\)
 *         greater than (resp. greater than or equal to) {@code d}
 */
public double exactP(double d, int n, int m, boolean strict) {
    Iterator<int[]> combinationsIterator = CombinatoricsUtils.combinationsIterator(n + m, n);
    long tail = 0;
    final double[] nSet = new double[n];
    final double[] mSet = new double[m];
    while (combinationsIterator.hasNext()) {
        // Generate an n-set
        final int[] nSetI = combinationsIterator.next();
        // Copy the n-set to nSet and its complement to mSet
        int j = 0;
        int k = 0;
        for (int i = 0; i < n + m; i++) {
            if (j < n && nSetI[j] == i) {
                nSet[j++] = i;
            } else {
                mSet[k++] = i;
            }
        }
        final double curD = kolmogorovSmirnovStatistic(nSet, mSet);
        if (curD > d) {
            tail++;
        } else if (curD == d && !strict) {
            tail++;
        }
    }
    return (double) tail / (double) CombinatoricsUtils.binomialCoefficient(n + m, n);
}
 
Example #11
Source File: TriangleListingTest.java    From flink with Apache License 2.0 5 votes vote down vote up
@Test
public void testCompleteGraph() throws Exception {
	long expectedDegree = completeGraphVertexCount - 1;
	long expectedCount = completeGraphVertexCount * CombinatoricsUtils.binomialCoefficient((int) expectedDegree, 2) / 3;

	DataSet<Result<LongValue>> tl = completeGraph
		.run(new TriangleListing<>());

	Checksum checksum = new ChecksumHashCode<Result<LongValue>>()
		.run(tl)
		.execute();

	assertEquals(expectedCount, checksum.getCount());
}
 
Example #12
Source File: GlobalClusteringCoefficientTest.java    From flink with Apache License 2.0 5 votes vote down vote up
@Test
public void testWithCompleteGraph() throws Exception {
	long expectedDegree = completeGraphVertexCount - 1;
	long expectedCount = completeGraphVertexCount * CombinatoricsUtils.binomialCoefficient((int) expectedDegree, 2);

	validate(completeGraph, expectedCount, expectedCount);
}
 
Example #13
Source File: TriadicCensusTest.java    From flink with Apache License 2.0 5 votes vote down vote up
@Test
public void testWithCompleteGraph() throws Exception {
	long expectedDegree = completeGraphVertexCount - 1;
	long expectedCount = completeGraphVertexCount * CombinatoricsUtils.binomialCoefficient((int) expectedDegree, 2) / 3;

	Result expectedResult = new Result(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, expectedCount);

	Result triadCensus = new TriadicCensus<LongValue, NullValue, NullValue>()
		.run(completeGraph)
		.execute();

	assertEquals(expectedResult, triadCensus);
}
 
Example #14
Source File: LocalClusteringCoefficientTest.java    From flink with Apache License 2.0 5 votes vote down vote up
@Test
public void testCompleteGraph() throws Exception {
	long degree = completeGraphVertexCount - 1;
	long triangleCount = CombinatoricsUtils.binomialCoefficient((int) degree, 2);

	validate(completeGraph, completeGraphVertexCount, degree, triangleCount);
}
 
Example #15
Source File: PolynomialsUtilsTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
@Test
public void testJacobiEvaluationAt1() {
    for (int v = 0; v < 10; ++v) {
        for (int w = 0; w < 10; ++w) {
            for (int i = 0; i < 10; ++i) {
                PolynomialFunction jacobi = PolynomialsUtils.createJacobiPolynomial(i, v, w);
                double binomial = CombinatoricsUtils.binomialCoefficient(v + i, i);
                Assert.assertTrue(Precision.equals(binomial, jacobi.value(1.0), 1));
            }
        }
    }
}
 
Example #16
Source File: TriadicCensusTest.java    From flink with Apache License 2.0 5 votes vote down vote up
@Test
public void testWithCompleteGraph() throws Exception {
	long expectedDegree = completeGraphVertexCount - 1;
	long expectedCount = completeGraphVertexCount * CombinatoricsUtils.binomialCoefficient((int) expectedDegree, 2) / 3;

	Result expectedResult = new Result(0, 0, 0, expectedCount);

	Result triadCensus = new TriadicCensus<LongValue, NullValue, NullValue>()
		.run(completeGraph)
		.execute();

	assertEquals(expectedResult, triadCensus);
}
 
Example #17
Source File: EdgeMetricsTest.java    From flink with Apache License 2.0 5 votes vote down vote up
@Test
public void testWithCompleteGraph() throws Exception {
	long expectedDegree = completeGraphVertexCount - 1;
	long expectedMaximumTriplets = CombinatoricsUtils.binomialCoefficient((int) expectedDegree, 2);
	long expectedTriplets = completeGraphVertexCount * expectedMaximumTriplets;

	Result expectedResult = new Result(expectedTriplets / 3, 2 * expectedTriplets / 3,
		expectedMaximumTriplets, expectedMaximumTriplets);

	validate(completeGraph, expectedResult);
}
 
Example #18
Source File: EdgeMetricsTest.java    From flink with Apache License 2.0 5 votes vote down vote up
@Test
public void testWithCompleteGraph() throws Exception {
	long expectedDegree = completeGraphVertexCount - 1;
	long expectedMaximumTriplets = CombinatoricsUtils.binomialCoefficient((int) expectedDegree, 2);
	long expectedTriplets = completeGraphVertexCount * expectedMaximumTriplets;

	Result expectedResult = new Result(expectedTriplets / 3, 2 * expectedTriplets / 3,
		expectedMaximumTriplets, expectedMaximumTriplets);

	validate(completeGraph, expectedResult);
}
 
Example #19
Source File: VertexMetricsTest.java    From flink with Apache License 2.0 5 votes vote down vote up
@Test
public void testWithCompleteGraph() throws Exception {
	long expectedDegree = completeGraphVertexCount - 1;
	long expectedEdges = completeGraphVertexCount * expectedDegree / 2;
	long expectedMaximumTriplets = CombinatoricsUtils.binomialCoefficient((int) expectedDegree, 2);
	long expectedTriplets = completeGraphVertexCount * expectedMaximumTriplets;

	Result expectedResult = new Result(completeGraphVertexCount, expectedEdges, expectedTriplets,
		expectedDegree, expectedMaximumTriplets);

	validate(completeGraph, false, expectedResult, expectedDegree, 1.0f);
}
 
Example #20
Source File: PascalDistribution.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
/** {@inheritDoc} */
public double probability(int x) {
    double ret;
    if (x < 0) {
        ret = 0.0;
    } else {
        ret = CombinatoricsUtils.binomialCoefficientDouble(x +
              numberOfSuccesses - 1, numberOfSuccesses - 1) *
              FastMath.pow(probabilityOfSuccess, numberOfSuccesses) *
              FastMath.pow(1.0 - probabilityOfSuccess, x);
    }
    return ret;
}
 
Example #21
Source File: ApacheCommonsCombinationGenerator.java    From tutorials with MIT License 5 votes vote down vote up
/** 
 * Print all combinations of r elements from a set
 * @param n - number of elements in set
 * @param r - number of elements in selection
 */
public static void generate(int n, int r) {
    Iterator<int[]> iterator = CombinatoricsUtils.combinationsIterator(n, r);
    while (iterator.hasNext()) {
        final int[] combination = iterator.next();
        System.out.println(Arrays.toString(combination));
    }
}
 
Example #22
Source File: DerivativeStructureTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
@Override
@Test
public void testLog() {
    double[] epsilon = new double[] { 1.0e-16, 1.0e-16, 3.0e-14, 7.0e-13, 3.0e-11 };
    for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
        for (double x = 0.1; x < 1.2; x += 0.001) {
            DerivativeStructure log = new DerivativeStructure(1, maxOrder, 0, x).log();
            Assert.assertEquals(FastMath.log(x), log.getValue(), epsilon[0]);
            for (int n = 1; n <= maxOrder; ++n) {
                double refDer = -CombinatoricsUtils.factorial(n - 1) / FastMath.pow(-x, n);
                Assert.assertEquals(refDer, log.getPartialDerivative(n), epsilon[n]);
            }
        }
    }
}
 
Example #23
Source File: DerivativeStructureTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
@Test
public void testReciprocal() {
    for (double x = 0.1; x < 1.2; x += 0.1) {
        DerivativeStructure r = new DerivativeStructure(1, 6, 0, x).reciprocal();
        Assert.assertEquals(1 / x, r.getValue(), 1.0e-15);
        for (int i = 1; i < r.getOrder(); ++i) {
            double expected = ArithmeticUtils.pow(-1, i) * CombinatoricsUtils.factorial(i) /
                              FastMath.pow(x, i + 1);
            Assert.assertEquals(expected, r.getPartialDerivative(i), 1.0e-15 * FastMath.abs(expected));
        }
    }
}
 
Example #24
Source File: DSCompilerTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
@Test
public void testSize() {
    for (int i = 0; i < 6; ++i) {
        for (int j = 0; j < 6; ++j) {
            long expected = CombinatoricsUtils.binomialCoefficient(i + j, i);
            Assert.assertEquals(expected, DSCompiler.getCompiler(i, j).getSize());
            Assert.assertEquals(expected, DSCompiler.getCompiler(j, i).getSize());
        }
    }
}
 
Example #25
Source File: PascalDistribution.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
/** {@inheritDoc} */
public double probability(int x) {
    double ret;
    if (x < 0) {
        ret = 0.0;
    } else {
        ret = CombinatoricsUtils.binomialCoefficientDouble(x +
              numberOfSuccesses - 1, numberOfSuccesses - 1) *
              FastMath.pow(probabilityOfSuccess, numberOfSuccesses) *
              FastMath.pow(1.0 - probabilityOfSuccess, x);
    }
    return ret;
}
 
Example #26
Source File: BetaBinomialDistribution.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
/**
 * @param k number of successes.  Must be positive or zero.
 * @return the value of the pdf at k
 */
@Override
public double logProbability(int k) {
    ParamUtils.isPositiveOrZero(k, "Number of successes must be greater than or equal to zero.");

    // nchoosek * beta(k+alpha, n-k+beta)/ beta(alpha, beta)
    // binomialcoefficient is "n choose k"
    return k > n ? Double.NEGATIVE_INFINITY :
            CombinatoricsUtils.binomialCoefficientLog(n,k) + Beta.logBeta(k+alpha, n - k + beta) - Beta.logBeta(alpha, beta);
}
 
Example #27
Source File: DerivativeStructureTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
@Override
@Test
public void testLog() {
    double[] epsilon = new double[] { 1.0e-16, 1.0e-16, 3.0e-14, 7.0e-13, 3.0e-11 };
    for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
        for (double x = 0.1; x < 1.2; x += 0.001) {
            DerivativeStructure log = new DerivativeStructure(1, maxOrder, 0, x).log();
            Assert.assertEquals(FastMath.log(x), log.getValue(), epsilon[0]);
            for (int n = 1; n <= maxOrder; ++n) {
                double refDer = -CombinatoricsUtils.factorial(n - 1) / FastMath.pow(-x, n);
                Assert.assertEquals(refDer, log.getPartialDerivative(n), epsilon[n]);
            }
        }
    }
}
 
Example #28
Source File: DerivativeStructureTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
@Test
public void testReciprocal() {
    for (double x = 0.1; x < 1.2; x += 0.1) {
        DerivativeStructure r = new DerivativeStructure(1, 6, 0, x).reciprocal();
        Assert.assertEquals(1 / x, r.getValue(), 1.0e-15);
        for (int i = 1; i < r.getOrder(); ++i) {
            double expected = ArithmeticUtils.pow(-1, i) * CombinatoricsUtils.factorial(i) /
                              FastMath.pow(x, i + 1);
            Assert.assertEquals(expected, r.getPartialDerivative(i), 1.0e-15 * FastMath.abs(expected));
        }
    }
}
 
Example #29
Source File: DSCompilerTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
@Test
public void testSize() {
    for (int i = 0; i < 6; ++i) {
        for (int j = 0; j < 6; ++j) {
            long expected = CombinatoricsUtils.binomialCoefficient(i + j, i);
            Assert.assertEquals(expected, DSCompiler.getCompiler(i, j).getSize());
            Assert.assertEquals(expected, DSCompiler.getCompiler(j, i).getSize());
        }
    }
}
 
Example #30
Source File: KolmogorovSmirnovTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
/**
 * Computes \(P(D_{n,m} > d)\) if {@code strict} is {@code true}; otherwise \(P(D_{n,m} \ge
 * d)\), where \(D_{n,m}\) is the 2-sample Kolmogorov-Smirnov statistic. See
 * {@link #kolmogorovSmirnovStatistic(double[], double[])} for the definition of \(D_{n,m}\).
 * <p>
 * The returned probability is exact, obtained by enumerating all partitions of {@code m + n}
 * into {@code m} and {@code n} sets, computing \(D_{n,m}\) for each partition and counting the
 * number of partitions that yield \(D_{n,m}\) values exceeding (resp. greater than or equal to)
 * {@code d}.
 * </p>
 * <p>
 * <strong>USAGE NOTE</strong>: Since this method enumerates all combinations in \({m+n} \choose
 * {n}\), it is very slow if called for large {@code m, n}. For this reason,
 * {@link #kolmogorovSmirnovTest(double[], double[])} uses this only for {@code m * n < }
 * {@value #SMALL_SAMPLE_PRODUCT}.
 * </p>
 *
 * @param d D-statistic value
 * @param n first sample size
 * @param m second sample size
 * @param strict whether or not the probability to compute is expressed as a strict inequality
 * @return probability that a randomly selected m-n partition of m + n generates \(D_{n,m}\)
 *         greater than (resp. greater than or equal to) {@code d}
 */
public double exactP(double d, int n, int m, boolean strict) {
    Iterator<int[]> combinationsIterator = CombinatoricsUtils.combinationsIterator(n + m, n);
    long tail = 0;
    final double[] nSet = new double[n];
    final double[] mSet = new double[m];
    while (combinationsIterator.hasNext()) {
        // Generate an n-set
        final int[] nSetI = combinationsIterator.next();
        // Copy the n-set to nSet and its complement to mSet
        int j = 0;
        int k = 0;
        for (int i = 0; i < n + m; i++) {
            if (j < n && nSetI[j] == i) {
                nSet[j++] = i;
            } else {
                mSet[k++] = i;
            }
        }
        final double curD = kolmogorovSmirnovStatistic(nSet, mSet);
        if (curD > d) {
            tail++;
        } else if (curD == d && !strict) {
            tail++;
        }
    }
    return (double) tail / (double) CombinatoricsUtils.binomialCoefficient(n + m, n);
}