Skip to content

Commit

Permalink
fixed files form Math #92
Browse files Browse the repository at this point in the history
  • Loading branch information
tdurieux committed Mar 7, 2017
1 parent 5a7c08c commit 224d218
Showing 1 changed file with 49 additions and 6 deletions.
55 changes: 49 additions & 6 deletions projects/Math/92/org/apache/commons/math/util/MathUtils.java
Original file line number Diff line number Diff line change
Expand Up @@ -181,30 +181,43 @@ public static long binomialCoefficient(final int n, final int k) {
if ((k == 1) || (k == n - 1)) {
return n;
}
long result = Math.round(binomialCoefficientDouble(n, k));
if (result == Long.MAX_VALUE) {
throw new ArithmeticException(
"result too large to represent in a long integer");
}
// Use symmetry for large k
if (k > n / 2)
return binomialCoefficient(n, n - k);

// We use the formula
// (n choose k) = n! / (n-k)! / k!
// (n choose k) == ((n-k+1)*...*n) / (1*...*k)
// which could be written
// (n choose k) == (n-1 choose k-1) * n / k
long result = 1;
if (n <= 61) {
// For n <= 61, the naive implementation cannot overflow.
for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
result = result * i / j;
}
} else if (n <= 66) {
// For n > 61 but n <= 66, the result cannot overflow,
// but we must take care not to overflow intermediate values.
for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
// We know that (result * i) is divisible by j,
// but (result * i) may overflow, so we split j:
// Filter out the gcd, d, so j/d and i/d are integer.
// result is divisible by (j/d) because (j/d)
// is relative prime to (i/d) and is a divisor of
// result * (i/d).
long d = gcd(i, j);
result = (result / (j / d)) * (i / d);
}
} else {
// For n > 66, a result overflow might occur, so we check
// the multiplication, taking care to not overflow
// unnecessary.
for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
long d = gcd(i, j);
result = mulAndCheck((result / (j / d)), (i / d));
}
}
return result;
}

Expand All @@ -231,9 +244,33 @@ public static long binomialCoefficient(final int n, final int k) {
* @throws IllegalArgumentException if preconditions are not met.
*/
public static double binomialCoefficientDouble(final int n, final int k) {
if (n < k) {
throw new IllegalArgumentException(
"must have n >= k for binomial coefficient (n,k)");
}
if (n < 0) {
throw new IllegalArgumentException(
"must have n >= 0 for binomial coefficient (n,k)");
}
if ((n == k) || (k == 0)) {
return 1d;
}
if ((k == 1) || (k == n - 1)) {
return n;
}
if (k > n/2) {
return binomialCoefficientDouble(n, n - k);
}
if (n < 67) {
return binomialCoefficient(n,k);
}

double result = 1d;
for (int i = 1; i <= k; i++) {
result *= (double)(n - k + i) / (double)i;
}

return Math.floor(Math.exp(binomialCoefficientLog(n, k)) + 0.5);
return Math.floor(result + 0.5);
}

/**
Expand Down Expand Up @@ -274,11 +311,17 @@ public static double binomialCoefficientLog(final int n, final int k) {
* For values small enough to do exact integer computation,
* return the log of the exact value
*/
if (n < 67) {
return Math.log(binomialCoefficient(n,k));
}

/*
* Return the log of binomialCoefficientDouble for values that will not
* overflow binomialCoefficientDouble
*/
if (n < 1030) {
return Math.log(binomialCoefficientDouble(n, k));
}

/*
* Sum logs for values that could overflow
Expand Down

0 comments on commit 224d218

Please sign in to comment.