Solve set-counting using generating functions

[ajcr]

Currently molloy builds arrays of polynomial coefficients from constraints, and then multiplies them to get a final array whereupon the requested coefficient can be read off as the answer.

This is simple and and is quite fast, as long as the user does not ask for too large a count (say much higher than 1 million) and not too many constraints are put upon items in the possible collection (so we do not have to too much polynomial multiplication [= convolution] and flipping of ones and zeros in arrays).

However, rather than build arrays, it would be potentially much more efficient and useful to construct generating functions for each constraint, multiply them, and then grab the needed coefficient from the corresponding power series expansion of that function.

Using SymPy would make the first two steps relatively simple, but the last part is still difficult. SymPy has a built-in method for extracting the Taylor term of a function expansion, but it uses repeated differentiation and is incredibly slow for anything but tiny integers.

A more promising approach is to take the generating function, expand it into partial fractions and then recognise the coefficient formula for each one. For example, (1 - x)^-(k+1) expands to a series where x^n has coefficient (n+k) C n and so this coefficient can be computed very quickly.

So in summary there are four-ish steps to making this work:

• build a gf for each item given the constrainsts
• multiply these gfs and expand to partial fractions
• recognise the form of each fraction and compute coefficients using a known formula
• add/subtract the relevant coefficients from each partial fraction to get the final answer

Whether this will always work, I'm not yet sure. I have made some progress, but getting the code to work reliably is a challenge - it can be hard to properly expand the function and then recognise the appropriate formula for the fraction.