Efficient generation of Taylor (Maclaurin) series

Question

Consider function
y=1/((1-x^5)(1-x^7)(1-x^11))

WolframAlpha computes first 1000 elements of the MacLaurin series expansion in a few seconds:
https://www.wolframalpha.com/input/?i=maclaurin+series+1%2F%28%281-x%5E5%29%281-x%5E7%29%281-x%5E11%29%29

Out of curiosity I wrote a very naive java program to do the same using BigInteger for polynomial coefficients. In pseudocode it would be something like:

BigInt next=1;
BigInt factorial=1;
while(true){
   function=function.differentiate();
   factorial*=++next;
   print("Next coefficient is: " + function(0)/factorial);
}

This program crashes with java.lang.outofmemory exception after computing first seven, or so, coefficients, because numerator and denominator of the fraction become enormously long polynomials. Suppose my code is inefficient, but still it does not seem like Wolfram is using the same technique they show you if the first year calculus class.
The question is: what does Wolfram use?

For comparison Wolfram takes quite a bit more time to just compute tenth derivative of the same function than it takes to get the first 1000 terms of polynomial, which, if done naively, would require differentiating the function 1000 times.
https://www.wolframalpha.com/input/?i=tenth+derivative+1%2F%28%281-x%5E5%29%281-x%5E7%29%281-x%5E11%29%29

Answer 1

tl;dr: The coefficient of x^N is the number of ways that N can be partitioned using only 5, 7, and 11.

I’m not sure how Wolfram does it, but for this function, it is possible to more efficiently compute the coefficients (using techniques you would see at the end of your first year in calculus). As a power series, 1/(1-x)=∑_k=0^∞ x^k. But we can replace x with xⁿ, and the relation will still hold. This means that
1/((1-x⁵)(1-x⁷)(1-x¹¹)) = (∑_k=0^∞x^5k)(∑_k=0^∞x^7k)(∑_k=0^∞x^11k)

Multiplying this out would be a pain. But all of the coefficients are 1, so we only need to look at the exponents, which add together. For example, Wolfram shows that the coefficient of x⁴⁰ is 4, which is from (x^5·1)(x^7·5)(x^0·11)+(x^5·0)(x^7·1)(x^11·3)+(x^5·3)(x^7·2)(x^11·1)+(x^5·8)(x^7·0)(x^11·0).

But if we only need to add the exponents, then we don’t need to care about the coefficients or the variable x. In the end, the coefficient of x^N is the number of ways that N can be written as a sum of 5s, 7s, and 11s. This is a restricted version of the partition problem, but the same ideas would still hold. In particular, a dynamic programming approach would be able to calculate coefficients in linear time and space.

Answer 2

Not sure about the fraction's numerator, but I can see why its denominator is growing much too fast:

factorial*=factorial+1;

is not how you calculate a factorial. That more than squares the "factorial" value on the denominator with each iteration! So you will get 1, 2, 6, 42, 1806, 3263442... By contrast, factorials go 1, 2, 6, 24, 120, 720...

To calculate the factorial incrementally, maintain a loop counter, and multiply factorial by that each time.