Calculate all possible distributions of n objects to k containers

Question

Let's say I have \$5\$ apples, thus \$n = 5\$. I have three bags (\$k = 3\$), each having the capacities of \$4\$, \$4\$ and \$2\$:

$$c = { \{4, 4, 2 \}}$$

I'd like to calculate the number of ways the \$5\$ apples can be distributed into those bags; each bag can be empty. The results for this example would be following output:

$$\{4,1,0\}\\ \{4,0,1\}\\ \{1,4,0\}\\ \{0,4,1\}\\ \{3,2,0\}\\ \{3,0,2\}\\ \{2,3,0\}\\ \{0,3,2\}\\ \{3,1,1\}\\ \{1,3,1\}\\ \{2,2,1\}\\ \{2,1,2\}\\ \{1,2,2\}\\$$

Thus there are \$13\$ possible distributions for \$n = 5\$, \$k = 3\$ and \$c = \{4, 4, 2\}\$.

I have written a simple program which can calculate this, the basic algorithm is like following:

• Generate all integer partitions for \$n\$. I'm optimizing this by limiting each partition to each number in \$c\$; as such, the partition \$\{5, 0, 0\}\$ is impossible. For the example above, this would generate following partitions:

$$\{4,1,0\}\\ \{3,2,0\}\\ \{3,1,1\}\\ \{2,2,1\}\\$$

• Permute each partition. Each partition will have \$(\frac{k!}{i!})\$ possible permutations, whereas \$i\$ describes the number of identical integers. The set \$\{3, 1, 1\}\$ would have precisely \$3\$ permutations, as \$(\frac{3!}{2!} = 3)\$:

$$\\\{3,1,1\}\\ \{1,3,1\}\\ \{1,1,3\}\\$$

• After each permutation is calculated, validate them based on \$c\$. The permutation \$\{1, 1, 3\}\$ would be invalid, as it wouldn't fit into \$c\$.
• Finally, if the permutation is valid, increase the counter.

---

I'm storing the partitions in a BlockingCollection to avoid running out of memory in case the numbers get too large (permutations will always be slower than partitions): \$n = 1000\$ could potentially yield \$2.4 \cdot 10^{31}\$ partitions!

My code works well for small numbers. Another example would yield \$48\$ possible distributions:

$$\\n = 10\\ k = 3\\ c = \{10,8,5\}$$

---

The problem of my approach is that it isn't efficient for large numbers at all. Following input could take days, if not weeks to compute:

$$\\n = 30\\ k = 20\\ c = \{20,19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1\}$$

I'm looking for ways to optimize my code, so that calculations for input like that would be feasible.

internal class Program
{
    static int n;
    static int k;
    static int[] c;

    static BlockingCollection<int[]> partitions;
    static double counter;

    static void Main(string[] args)
    {
        n = 10;
        k = 3;
        c = new[] {10, 8, 5}.OrderByDescending(x => x).ToArray();

        partitions = new BlockingCollection<int[]>(8);

        Task partitionTask = Task.Factory.StartNew(() =>
        {
            Partition(n, c[0], new List<int>());
        });
        partitionTask.ContinueWith(delegate { partitions.CompleteAdding(); });

        foreach (int[] x in partitions.GetConsumingEnumerable())
        {
            Permute(x);
        }

        Console.WriteLine($"There are {counter:n0} possible distributions.");
        Console.ReadLine();
    }

    // Generate all partitions of c
    static void Partition(int x, int limit, List<int> p)
    {
        if (x > 0)
        {
            for (int i = Math.Min(x, limit); i > 0; i--)
            {
                if (p.Count == c.Length)
                    continue;

                Partition(x - i, Math.Min(i, p.Count > c.Length - 2 ? i : c[p.Count + 1]), p.Concat(new[] { i }).ToList());
            }
        }
        else
        {
            while (p.Count < k)
                p.Add(0);

            partitions.Add(p.ToArray());
        }
    }

    // Generate all permutations for a given set and validate them
    static void Permute(int[] set)
    {
        Action<int> permute = null;
        permute = start =>
        {
            if (start == set.Length)
            {
                // Validate this permutation, if valid, increase the counter
                if (!set.Where((t, i) => t > c[i]).Any())
                {
                    counter++;
                }
            }
            else
            {
                List<int> swaps = new List<int>();
                for (int i = start; i < set.Length; i++)
                {
                    if (swaps.Contains(set[i])) continue;
                    swaps.Add(set[i]);

                    Swap(set, start, i);
                    permute(start + 1);
                    Swap(set, start, i);
                }
            }
        };
        permute(0);
    }

    static void Swap(int[] set, int index1, int index2)
    {
        int temp = set[index1];
        set[index1] = set[index2];
        set[index2] = temp;
    }
}

Answer 1

The structure of the code makes it difficult to test. Suppose you had two different implementations, and wanted to test that they returned the same results? The code doesn't lend itself well to this situation.

My first attempt at speeding up the code was to use the answer to TopinFrassi's question on Mathematics.SE.

As the answer explains, the solution is given by the coefficient of \$x^n\$ in

$$\prod_{c \in C} f_c(x),$$

where \$f_c(x) = 1 + x + x^2 + \cdots + x^c\$.

First, a bare-bones implementation of Polynomial:

public struct Polynomial
{
    private readonly IReadOnlyList<BigInteger> _coefficients;

    public Polynomial(IEnumerable<BigInteger> coefficients)
    {
        if (coefficients == null)
        {
            throw new ArgumentNullException(nameof(coefficients));
        }

        _coefficients = coefficients.ToArray();
    }

    public static readonly Polynomial One = new Polynomial(new[] { BigInteger.One });

    public static Polynomial operator *(Polynomial a, Polynomial b)
    {
        var degree = a.Degree + b.Degree;
        var coefficients = new BigInteger[degree + 1];
        for (var k = 0; k <= degree; k++)
        {
            coefficients[k] = Enumerable.Range(0, k + 1)
                .Select(i => a[i] * b[k - i])
                .Sum();
        }

        return new Polynomial(coefficients);
    }

    public BigInteger this[int n] => n >= 0 && n < _coefficients.Count ? _coefficients[n] : BigInteger.Zero;

    private int Degree => _coefficients.Count - 1;
}

Some helper methods:

internal static class EnumerableExtensions
{
    public static BigInteger Sum(this IEnumerable<BigInteger> source)
    {
        return source.Aggregate(BigInteger.Zero, (x, y) => x + y);
    }

    public static Polynomial Product(this IEnumerable<Polynomial> source)
    {
        return source.Aggregate(Polynomial.One, (x, y) => x * y);
    }
}

And putting it all together:

public static BigInteger Solve(int n, IEnumerable<int> bags)
{
    return bags.Select(bag => new Polynomial(Enumerable.Repeat(BigInteger.One, bag + 1)))
               .Product()[n];
}

This works fine for your example of \$n = 30\$, \$C = \left\{ 1, 2, 3, \ldots, 20 \right\}\$, calculating the solution \$6209623185136\$ in about 0.03s (solution cross-checked against Joe Wallis's answer).

I then translated Joe's answer to C#, only to find that it performed faster still, taking about 0.01s for the same problem.

The difference got worse as I tried bigger numbers. For example, Joe's solution takes just 0.8s to solve \$n = 30\$ and \$C = \left\{ 1, 2, 3, \ldots, 60 \right\}\$, while my solution above takes about 4.4s.

My next attempt was based on the fact that we only need one coefficient of the product. If we lazily calculate the coefficients of the polynomials, can we cut down on the time the method takes?

To delay evaluation, I changed the coefficients from BigIntegers to Lazy<BigInteger>s:

public struct LazyPolynomial
{
    private readonly IReadOnlyList<Lazy<BigInteger>> _coefficients;

    public LazyPolynomial(IEnumerable<Lazy<BigInteger>> coefficients)
    {
        if (coefficients == null)
        {
            throw new ArgumentNullException(nameof(coefficients));
        }

        _coefficients = coefficients.ToArray();
    }

    public static readonly LazyPolynomial One = new LazyPolynomial(new[] { new Lazy<BigInteger>(() => BigInteger.One) });

    public static LazyPolynomial operator *(LazyPolynomial a, LazyPolynomial b)
    {
        var degree = a.Degree + b.Degree;
        var coefficients = new Lazy<BigInteger>[degree + 1];
        for (var k = 0; k <= degree; k++)
        {
            var tmp = k;
            coefficients[k] = new Lazy<BigInteger>(() => Enumerable.Range(0, tmp + 1)
                .Select(i => a[i] * b[tmp - i])
                .Sum());
        }

        return new LazyPolynomial(coefficients);
    }

    public BigInteger this[int n] => n >= 0 && n < _coefficients.Count ? _coefficients[n].Value : BigInteger.Zero;

    private int Degree => _coefficients.Count - 1;
}

This turned out to be considerably faster than both my first attempt, and my translation of Joe's code. Some timings for \$n = 30\$ are below:

$$ \begin{array}{lrr} & C = \left\{ 1, \ldots, 20 \right\} & C = \left\{ 1, \ldots, 60 \right\} \\ \hline \text{Polynomial} & 0.032s & 4.39s \\ \text{Translated} & 0.010s & 0.80s \\ \text{LazyPolynomial} & 0.002s & 0.05s \\ \end{array}$$

---

I've now got some more accurate timings using BenchmarkDotNet:

```ini
BenchmarkDotNet=v0.7.8.0
OS=Microsoft Windows NT 6.2.9200.0
Processor=Intel(R) Core(TM) i7-4712HQ CPU @ 2.30GHz, ProcessorCount=8
HostCLR=MS.NET 4.0.30319.42000, Arch=32-bit
Type=ThirtyTwentyCompetition  Mode=Throughput  Platform=HostPlatform  Jit=HostJit  .NET=HostFramework
```

                     Method |    AvrTime |    StdDev |   op/s |
--------------------------- |----------- |---------- |------- |
 PolynomialLazyThirtyTwenty |  1.2705 ms | 0.0114 ms | 787.08 |
     PolynomialThirtyTwenty | 13.4990 ms | 0.2331 ms |  74.08 |
     TranslatedThirtyTwenty |  9.0769 ms | 0.4108 ms | 110.17 |

                    Method |       AvrTime |      StdDev |  op/s |
-------------------------- |-------------- |------------ |------ |
 PolynomialLazyThirtySixty |    12.3559 ms |   0.2563 ms | 80.93 |
     PolynomialThirtySixty | 2,610.2568 ms | 138.6060 ms |  0.38 |
     TranslatedThirtySixty |   712.4170 ms |  40.4284 ms |  1.40 |

Answer 2

I'd approach this in a different way. Looking at the problem I'd sum the different bag capacities until they are \$<= n\$. The moment the sum becomes \$> n\$ I have to stop incresing this sum as it would be useless.

Doing this kind of reasoning I'd go with a tree structure. The max level of the tree should be the number \$k\$ of bags/buckets. The spanning factor for each node in the i-th level is bags[i].capacity + 1 (from 0 to max capacity of the bag).

So, the operations of the algorithm should be the following:

1. Start with an empty tree with a single root
2. Expand the root (insert the nodes for the first level) with the possible elements of the first bag
3. If sum from root to current is > n then there is no need to proceed further on this path. If sum from root to current is == n increase counter and stop expanding this path. If sum from root to current is < n repeat point 2 with the current node instead of the root and by expanding the node with the possible elements of bag[currentLevel + 1].

Regarding the complexity of the algorithm we have that time complexity is \$O(AvgBagCapacity^{NumOfBags})\$ (that is reduced a little by breaking the computation early). Regarding the space complexity depends on the approach used to expand the tree. Personally, I'd suggest to use a depth-first approach as it reduces the space complexity to \$O(AvgBagCapacity * NumOfBags)\$ (if you don't keep the already expanded nodes in memory). If you want to furtherly reduce the space complexity you could use lazy initializations of the nodes so you use only one node at a time. By doing so, the space complexity is \$O(NumOfBags)\$.

For the sake of clarity I wrote some code which does what I described. Obviously it has a lot of room for improvement (starting with the naming :P ). I also included some tests (the first two are the same as two of your examples) so you can do comparisons and play with it.

public class Node
{
    public int[] Buckets { get; set; }
    public int Element { get; set; }
    public int ParentSum { get; set; }
    public int NodeSum { get { return Element + ParentSum; } }
    public int Level { get; set; }

    public IEnumerable<Node> GenerateChildrenForNextLevel()
    {
        if(Level < Buckets.Length - 1)
        {
            var thisNode = this;

            for(int i = 0; i <= Buckets[Level + 1]; i++)
            {
                yield return new Node()
                {
                    Element = i,
                    ParentSum = thisNode.NodeSum,
                    Level = thisNode.Level + 1, 
                    Buckets = thisNode.Buckets
                };
            }
        }

        yield break;
    }

    public int GetNumberOfPermutationsForSum(int sum)
    {
        int numberOfPermutations = 0;

        if(NodeSum == sum)
        {
            numberOfPermutations = 1;
        }
        else if(NodeSum < sum)
        {
            var children = GenerateChildrenForNextLevel();

            foreach(var child in children)
            {
                numberOfPermutations += child.GetNumberOfPermutationsForSum(sum);
            }
        }

        return numberOfPermutations;
    }
}

public class Tree
{
    public int[] Buckets { get; set; }
    public Node Root { get; set; }
    public int TotalSum { get; private set; }

    public Tree(int totalSum, int[] buckets)
    {
        Root = new Node()
        {
            Element = 0,
            ParentSum = 0,
            Level = -1,
            Buckets = buckets
        };
        TotalSum = totalSum;
        Buckets = buckets;
    }

    public int CalculateNumberOfPermutations()
    {
        return Root.GetNumberOfPermutationsForSum(TotalSum);
    }
}

class Program
{
    static void Main()
    {
        int[] buckets1 = new int[] { 4, 4, 2 };
        int[] buckets2 = new int[] { 10, 8, 5};
        int[] buckets3 = new int[] { 20, 19, 18, 17, 16, 
                                 15, 14 };
        Tree tree1 = new Tree(5, buckets1);
        Tree tree2 = new Tree(10, buckets2);
        Tree tree3 = new Tree(30, buckets3);
        int numberOfPermutations1, numberOfPermutations2, numberOfPermutations3;
        Stopwatch sw = new Stopwatch();

        sw.Start();
        numberOfPermutations1 = tree1.CalculateNumberOfPermutations();
        sw.Stop();

        Console.WriteLine("The number of permutations for the first example is {0} and was calculated in {1}",
                      numberOfPermutations1,
                      sw.Elapsed);

        sw.Start();
        numberOfPermutations2 = tree2.CalculateNumberOfPermutations();
        sw.Stop();

        Console.WriteLine("The number of permutations for the second example is {0} and was calculated in {1}",
                      numberOfPermutations2,
                      sw.Elapsed);

        sw.Start();
        numberOfPermutations3 = tree3.CalculateNumberOfPermutations();
        sw.Stop();

        Console.WriteLine("The number of permutations for the third example is {0} and was calculated in {1}",
                      numberOfPermutations3,
                      sw.Elapsed);
    }
}

Let me know if anything's unclear.