Comparing algorithms for computing binomial coefficients in Java

Question

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; 
                i = i.add(BigInteger.ONE)) {
            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;

        for (BigInteger i = n.subtract(k).add(BigInteger.ONE); 
                i.compareTo(n) <= 0; 
                i = i.add(BigInteger.ONE)) {
            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() {
        pascalsTriangle.add(new ArrayList<>());
        pascalsTriangle.add(new ArrayList<>());

        pascalsTriangle.get(0).add(BigInteger.ONE);
        pascalsTriangle.get(1).add(BigInteger.ONE);
        pascalsTriangle.get(1).add(BigInteger.ONE);
    }

    @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)
                .add(
               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);

        newRow.add(BigInteger.ONE);

        for (int index = 1; index < newRowLength - 1; ++index) {
            newRow.add(topRow.get(index - 1).add(topRow.get(index))); 
        }

        newRow.add(BigInteger.ONE);

        pascalsTriangle.add(newRow);
    }
}

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) -> { 
                     data.add(getRandomDatum(maxN, random)); 
                 });
        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) {
            outputList.add(computer.compute(datum.first, datum.second));
        }

        final long endTime = System.nanoTime();

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

        return outputList;
    }
}

Performance figures

I had these 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.

Answer 1

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.

Answer 2

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;
}

Answer 3

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; 
            i = i.add(BigInteger.ONE)) {
        ret = ret.multiply(i);
    }

    return ret;
}

with

    BigInteger numerator = BigInteger.ONE;

    for (BigInteger i = n.subtract(k).add(BigInteger.ONE); 
            i.compareTo(n) <= 0; 
            i = i.add(BigInteger.ONE)) {
        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; 
            i = i.add(BigInteger.ONE)) {
        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.