Comparing algorithms for computing binomial coefficients in Java

I have these 3 different algorithms for computing binomial coefficients (I also had the 4th recursive one, yet I discarded it since it is super slow). The first uses the factorial formula, the second optimizes it a bit, and the last is a dynamic programming algorithm that maintains a Pascal's triangle which reduces the computation to a single addition provided that the triangle is large enough (and if it is not, it is expanded rather efficiently). The formula behind the last algorithm is

$$\binom{n}{k} = \binom{n - 1}{k} + \binom{n - 1}{k - 1}.$$

You can think of the above like that $n$ selects the row of the Pascal's triangle (zero-based indexing), and $k, k - 1$ select two consecutive entries on a row (also zero-based indexing).

See what I have:

AbstractBinomialCoefficientComputer.java:

package net.coderodde.math;

import java.math.BigInteger;

/**
* This abstract class defines the API for computing binomial coefficients.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Jul 8, 2016)
*/
public abstract class AbstractBinomialCoefficientComputer {

/**
* Computes the binomial coefficient {@code n} over {@code k}.
*
* @param n the number of elements in the set.
* @param k the number of elements to choose.
* @return the number of distinct combinations when choosing {@code k} out
*         of {@code n} elements.
*/
public abstract BigInteger compute(final BigInteger n, final BigInteger k);

protected void checkArguments(final BigInteger n, final BigInteger k) {
if (n.compareTo(BigInteger.ZERO) < 0) {
throw new IllegalArgumentException("The 'n' is negative.");
}

if (k.compareTo(BigInteger.ZERO) < 0) {
throw new IllegalArgumentException("The 'k' is negative.");
}

if (k.compareTo(n) > 0) {
throw new IllegalArgumentException(
"The 'k' is larger than 'n'. (" + k + " > " + n + ").");
}
}

}


FactorialBinomialCoefficientComputer.java:

package net.coderodde.math.support;

import java.math.BigInteger;
import net.coderodde.math.AbstractBinomialCoefficientComputer;

/**
* This binomial coefficient computer computes the coefficients by means of
* factorial formula <tt>n! / (k! (n - k)!)</tt>. See
* <a href="https://en.wikipedia.org/wiki/Binomial_coefficient#Factorial_formula">Wikipedia</a>
* for more details.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Jul 8, 2016)
*/
public class FactorialBinomialCoefficientComputer
extends AbstractBinomialCoefficientComputer {

@Override
public BigInteger compute(final BigInteger n, final BigInteger k) {
checkArguments(n, k);

return factorial(n).divide(
factorial(k).multiply(
factorial(n.subtract(k))
)
);
}

static BigInteger factorial(final BigInteger number) {
BigInteger ret = BigInteger.ONE;

for (BigInteger i = BigInteger.valueOf(2L);
i.compareTo(number) <= 0;
ret = ret.multiply(i);
}

return ret;
}
}


MultiplicativeBinomialCoefficientComputer.java:

package net.coderodde.math.support;

import java.math.BigInteger;
import net.coderodde.math.AbstractBinomialCoefficientComputer;

/**
* This binomial coefficient computer computes the coefficients by means of a
* multiplicative formula described in
* <a href="https://en.wikipedia.org/wiki/Binomial_coefficient#Multiplicative_formula">Wikipedia</a>.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Jul 8, 2016)
*/
public class MultiplicativeBinomialCoefficientComputer
extends AbstractBinomialCoefficientComputer {

@Override
public BigInteger compute(final BigInteger n, final BigInteger k) {
checkArguments(n, k);

final BigInteger denominator =
FactorialBinomialCoefficientComputer.factorial(k);

BigInteger numerator = BigInteger.ONE;

i.compareTo(n) <= 0;
numerator = numerator.multiply(i);
}

return numerator.divide(denominator);
}
}


DynamicProgrammingBinomialCoefficientComputer.java:

package net.coderodde.math.support;

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;
import net.coderodde.math.AbstractBinomialCoefficientComputer;

/**
* This binomial coefficient computer computes the coefficients by means of a
* dynamic programming algorithm that caches the Pascal's triangle long enough
* for computing the coefficient. The triangle is expanded to accommodate more
* coefficients if needed. Given that the internal Pascal's triangle is large
* enough, computing a new coefficient is reduced to a single addition.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Jul 8, 2016)
*/
public class DynamicProgrammingBinomialCoefficientComputer
extends AbstractBinomialCoefficientComputer {

private final List<List<BigInteger>> pascalsTriangle = new ArrayList<>();

public DynamicProgrammingBinomialCoefficientComputer() {

}

@Override
public BigInteger compute(final BigInteger n, final BigInteger k) {
checkArguments(n, k);

if (k.equals(BigInteger.ZERO) || k.equals(n)) {
return BigInteger.ONE;
}

checkTriangle(n);

final int rowIndex = n.intValue() - 1;
final int colIndex = k.intValue() - 1;

return pascalsTriangle.get(rowIndex).get(colIndex)
pascalsTriangle.get(rowIndex).get(colIndex + 1)
);
}

private void checkTriangle(final BigInteger n) {
final int requestedN = n.intValue();

while (pascalsTriangle.size() < requestedN + 1) {
populatePascalsTriangleRow();
}
}

private void populatePascalsTriangleRow() {
final int newRowLength = pascalsTriangle.size() + 1;
final List<BigInteger> topRow = pascalsTriangle.get(
pascalsTriangle.size() - 1);
final List<BigInteger> newRow = new ArrayList<>(newRowLength);

for (int index = 1; index < newRowLength - 1; ++index) {
}

}
}


Demo.java:

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;
import java.util.Random;
import java.util.stream.IntStream;
import net.coderodde.math.AbstractBinomialCoefficientComputer;
import net.coderodde.math.support.DynamicProgrammingBinomialCoefficientComputer;
import net.coderodde.math.support.FactorialBinomialCoefficientComputer;
import net.coderodde.math.support.MultiplicativeBinomialCoefficientComputer;

public class Demo {

private static final int MAXIMUM_N = 1000;
private static final int SIZE = 10_000;

public static void main(final String[] args) {
final long seed = System.nanoTime();
final Random random = new Random(seed);
final List<Pair<BigInteger>> data = getRandomInputData(MAXIMUM_N,
SIZE,
random);

System.out.println("Seed = " + seed);

final List<BigInteger> result1 =
profile(new FactorialBinomialCoefficientComputer(), data);

final List<BigInteger> result2 =
profile(new MultiplicativeBinomialCoefficientComputer(), data);

final List<BigInteger> result3 =
profile(new DynamicProgrammingBinomialCoefficientComputer(),
data);

System.out.println("Algorithms agree: " +
(result1.equals(result2) && result2.equals(result3)));
}

private static final List<Pair<BigInteger>>
getRandomInputData(final int maxN,
final int size,
final Random random) {
final List<Pair<BigInteger>> data = new ArrayList<>(size);
IntStream.range(0, size)
.forEach((i) -> {
});
return data;
}

private static final Pair<BigInteger> getRandomDatum(final int maxN,
final Random random) {
final int n = random.nextInt(maxN + 1);
final int k = random.nextInt(n + 1);
return new Pair<>(BigInteger.valueOf(n), BigInteger.valueOf(k));
}

private static final class Pair<E> {
public final E first;
public final E second;

public Pair(final E first, final E second) {
this.first = first;
this.second = second;
}
}

private static List<BigInteger> profile(
final AbstractBinomialCoefficientComputer computer,
final List<Pair<BigInteger>> data) {
final List<BigInteger> outputList = new ArrayList<>(data.size());

final long startTime = System.nanoTime();

for (final Pair<BigInteger> datum : data) {
}

final long endTime = System.nanoTime();

System.out.printf("%s in %.0f milliseconds.\n",
computer.getClass().getSimpleName(),
(endTime - startTime) / 1e6);

return outputList;
}
}


Performance figures

Seed = 6873321663935
FactorialBinomialCoefficientComputer in 3262 milliseconds.
MultiplicativeBinomialCoefficientComputer in 1090 milliseconds.
DynamicProgrammingBinomialCoefficientComputer in 66 milliseconds.
Algorithms agree: true

Critique request

Please, tell me anything that comes to mind.

• Note that benchmarking Java is pretty hard and measuring anything below one second says about nothing about the performance for bigger problem instances. – maaartinus Jul 10 '16 at 0:37

Uses a lot of space

I tried increasing MAXIMUM_N to 2000, and I got this error:

Exception in thread "main" java.lang.OutOfMemoryError: Java heap space


Since you are populating Pascal's triangle, you are using $O(n^2)$ space, and furthermore, each slot in the triangle is a BigNumber that is getting increasingly bigger with each row, so it's actually around $O(n^3)$ space.

Alternative suggestion

I modified your FactorialBinomialCoefficientComputer solution to keep an ArrayList of previously computed factorials. This requires around $O(n^2)$ space as opposed to $O(n^3)$ space. It isn't quite as fast as the Pascal's triangle version because it needs to do a division + multiply + subtract to compute each answer. But it uses less space and is able to handle larger values of N.

With MAXIMUM_N at 1000 and SIZE at 10000, I got these times:

Seed = 1572945124315246
FactorialBinomialCoefficientComputer in 5405 milliseconds.
MultiplicativeBinomialCoefficientComputer in 2286 milliseconds.
DynamicProgrammingBinomialCoefficientComputer in 166 milliseconds.
CachingFactorialBinomialCoefficientComputer in 437 milliseconds.

With MAXIMUM_N at 2000 and SIZE at 5000, I got these times:

Seed = 1572804057554760
FactorialBinomialCoefficientComputer in 10718 milliseconds.
MultiplicativeBinomialCoefficientComputer in 4387 milliseconds.
(DynamicProgrammingBinomialCoefficientComputer ran out of heap) CachingFactorialBinomialCoefficientComputer in 900 milliseconds.

• Can you show me the code for your space-optimized version? – coderodde Jul 9 '16 at 9:56
• @coderodde I don't have the code with me at the moment, but it was pretty straightforward. When you attempt to compute the factorial of n but your ArrayList has only m entries, you add all the entries between m and n as you compute the factorial. – JS1 Jul 9 '16 at 10:12

You don't need to keep all of Pascal's triangle. In fact you need only keep one row. Below is a C routine that I use, I hope it explains the idea.

static  int64_t*    pascals_triangle( int N)
{
int n,k;
int64_t*    C = calloc( N+1, sizeof *C); // enough for row N
for( n=0; n<=N; ++n) // run through the rows, starting with 0
{   C[n] = 1; // C(n,n) = 1
k = n; // start at the end of the row
while( --k>0)
{    // apply the recurence
C[k] += C[k-1]; // C[k] on input is C(n-1,k)
// C[k-1] is C(n-1,k-1)
// on output, C[k] is C(n,k)
}
}
return C;
}


Missed opportunity for code reuse

Compare

static BigInteger factorial(final BigInteger number) {
BigInteger ret = BigInteger.ONE;

for (BigInteger i = BigInteger.valueOf(2L);
i.compareTo(number) <= 0;
ret = ret.multiply(i);
}

return ret;
}


with

    BigInteger numerator = BigInteger.ONE;

i.compareTo(n) <= 0;
numerator = numerator.multiply(i);
}


They could easily be refactored into a method

static BigInteger prod(final BigInteger from, final BigInteger to) {
BigInteger ret = BigInteger.ONE;

for (BigInteger i = from;
i.compareTo(to) <= 0;
ret = ret.multiply(i);
}

return ret;
}


I've left the special cases for you to decide how you would want to handle them.

Unnecessary calculation

DynamicProgrammingBinomialCoefficientComputer currently admits two optimisations other than those mentioned in earlier answers.

1. The recurrence is $\binom{n}{k} = \binom{n-1}{k} + \binom{n-1}{k-1}$. At no point do you need values of $\binom{m}{r}$ for $r > k$.

2. $\binom{n}{k} = \binom{n}{n-k}$, so if you exploit the first optimisation you can also ensure that $k \le \frac n 2$.

The combination can reduce the number of values calculated from $O(n^2)$ to $O(n \min(k, n-k))$.

A different approach

You can avoid BigInteger divisions and still get the speed benefit of multiplication by first computing the prime factorisation of $\binom{n}{k}$ in ints and then multiplying them in BigInteger. One approach would be Kummer's theorem; another is to simply compute the prime factorisations of $n!$, $k!$, and $(n-k)!$ and subtract appropriately. See e.g. this answer on a sister site.