CODEWITHSUNDEEP
algorithmcalculustaylor-series
QuarterlyQuotaOfQuotes
QuarterlyQuotaOfQuotes

Reputation: 416

Efficient generation of Taylor (Maclaurin) series

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

Upvotes: 3

Views: 2123

Answers (4)

Jules
Jules

Reputation: 6346

You can do each operation on formal power series. Given power series for f,g you can find a recurrence relation for power series of f(z)+g(z), f(z)g(z), f(z)/g(z), f(g(z)), and even f^-1(z). Using these methods you can compute the power series of practically any function in polynomial time.

In special cases there are more efficient methods. If f(z) has a power series, then coefficients of the power series of f(z)/(1 - z) are simply the partial sums of the power series of f(z). So if f_n is the series for f, then the series for g(z) = f(z)/(1 - z) is given by g_n = f_n + g_(n-1).

You can extend this to division by any polynomial. The algorithm is basically the same as long division for polynomials. For example, let's compute 1/(1 - z^2). We add and subtract z^2 to the numerator to get (1 - z^2 + z^2)/(1 - z^2) = 1 + z^2/(1 - z^2). Then we add and subtract z^4 to get (z^2 - z^4 + z^4)/(1 - z^2) = z^2 + z^4/(1 - z^2). Going on like that you find 1/(1 - z^2) = 1 + z^2 + z^4 + z^6 and so on.

When you do this for a general polynomial of degree n you always have a numerator with less than n terms. You can store the coefficients of those terms in an array and use that as your state. From a state you can compute the next term in the power series and the next state in time O(n). This gives you a O(nk) time algorithm to find the first k terms in the power series of 1/p(z).

Note that computing a power series at a point z=z0 is the same as finding all derivatives at z=z0, so the two problems are equivalent. You can compute power series at a symbolic variable point to find a formula for the derivative, so there is theoretically no reason why Wolfram should be so much slower at finding n-th derivatives.

Upvotes: 0

Teepeemm
Teepeemm

Reputation: 4518

tl;dr: The coefficient of xN 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∞ xk. But we can replace x with xn, and the relation will still hold. This means that
1/((1-x5)(1-x7)(1-x11)) = (∑k=0∞x5k)(∑k=0∞x7k)(∑k=0∞x11k)

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 x40 is 4, which is from (x5·1)(x7·5)(x0·11)+(x5·0)(x7·1)(x11·3)+(x5·3)(x7·2)(x11·1)+(x5·8)(x7·0)(x11·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 xN 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.

Upvotes: 3

user58697
user58697

Reputation: 7923

Rational functions (and this particular one) need neither differentiation nor factorials. One way to calculate the series is to expand each factor into the series of its own (e.g. 1/(1 - x^5) = sum(n=[0,inf] x^(5n)), and then multiply results as polynomials.

Upvotes: 0

j_random_hacker
j_random_hacker

Reputation: 51326

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.

Upvotes: 2

Related Questions