Somnium Somnium - 1 month ago 14
C++ Question

Best algorithm for series expansion of Rational function

I need to code function in C++ which efficiently finds coefficients of Taylor Series of given rational function (P(x) / Q(x)).

Function parameters will be power of polynomials (equal in nominator and denominator), two arrays with coefficients of polynomials and number of terms in expansion.

My idea was following.
Consider identity

P(x) / Q(x) = R(x) + ...


Where
R(x)
is a polynomial with number of terms equal to number of coefficients I need to find. Then I can multiply both sides with
Q(x)
and get

P(x) = R(x) * Q(x)

R(x) * Q(x) - P(x) = 0


Therefore, all coefficients should be zero. This is system of equations which have O(n^3) algorithm to solve. O(n^3) is not that fast as I wanted.

Is there any faster algorithm?

I know that coefficients of series are satisfying linear recurrence relation.
This makes me think that O(n) algorithm is possible.

Answer

The algorithm that I'm about to describe is justified mathematically by formal power series. Every function with a Taylor series has a formal power series. The converse is not true, but if we do arithmetic on functions with Taylor series and get a function with a Taylor series, then we can do the same arithmetic with formal power series and get the same answer.

The long division algorithm for formal power series is like the long division algorithm that you may have learned in school. I'll demonstrate it on the example (1 + 2 x)/(1 - x - x^2), which has coefficients equal to the Lucas numbers.

The denominator must have a nonzero constant term. We start by writing the numerator, which is the first residual.

             --------
1 - x - x^2 ) 1 + 2 x

[ We divide the residual's lowest-order term (1) by the denominator's constant term (1) and put the quotient up top.

              1
             --------
1 - x - x^2 ) 1 + 2 x

Now we multiply 1 - x - x^2 by 1 and subtract it from the current residual.

              1
             --------
1 - x - x^2 ) 1 + 2 x
              1 -   x - x^2
              -------------
                  3 x + x^2

Do it again.

              1 + 3 x
             --------
1 - x - x^2 ) 1 + 2 x
              1 -   x -   x^2
              ---------------
                  3 x +   x^2
                  3 x - 3 x^2 - 3 x^3
                  -------------------
                        4 x^2 + 3 x^3

And again.

              1 + 3 x + 4 x^2
             ----------------
1 - x - x^2 ) 1 + 2 x
              1 -   x -   x^2
              ---------------
                  3 x +   x^2
                  3 x - 3 x^2 - 3 x^3
                  -------------------
                        4 x^2 + 3 x^3
                        4 x^2 - 4 x^3 - 4 x^4
                        ---------------------
                                7 x^3 + 4 x^4

And again.

              1 + 3 x + 4 x^2 + 7 x^3
             ------------------------
1 - x - x^2 ) 1 + 2 x
              1 -   x -   x^2
              ---------------
                  3 x +   x^2
                  3 x - 3 x^2 - 3 x^3
                  -------------------
                        4 x^2 + 3 x^3
                        4 x^2 - 4 x^3 - 4 x^4
                        ---------------------
                                7 x^3 + 4 x^4
                                7 x^3 - 7 x^4 - 7 x^4
                                ---------------------
                                       11 x^4 + 7 x^5

The individual divisions were kind of boring because I used a divisor with a leading 1, but if I had used, say, 2 - 2 x - 2 x^2, then all of the coefficients in the quotient would be divided by 2.

Comments