Café Math : Basic Power Series Algorithms
1, 1, 2, 2, 1, 8, 6, 7, 14, 27, 26, 80, 133, 170, 348, 765, 1002, 2176, 4682, ...   (A121350)
Math Blog
Phd Thesis
Financial Maths
About me

Basic Power Series Algorithms

Thu, 05 Sep 2013 22:10:07 GMT

Hi, the other day I wanted to do some formal series computations for my next research article and I was tiered of using existing computer algebra programs. I decided to implement some fast formal power series algorithms in C#, just for fun. Today, I want to share what I came up with. As always feel free to comment and propose if you find anything better.

1. Basic Algorithms

Problem. How to compute addition and substraction, \begin{align*} C(x) &= A(x) + B(x), &&& A(x) &= C(x) - B(x). \end{align*}

This is the easy part. Addition and substraction of power series are done term by term, \begin{align*} c_n &= a_n + b_n, &&& a_n &= c_n - b_n. \end{align*}

Problem. How to compute differentiation and integration, \begin{align*} B(x) = A'(x), &&& A(x) = \int B(x) dx. \end{align*}

Well easy again, term by term, \begin{align*} b_{n} &= (n + 1) a_{n + 1}, &&& a_{n + 1} &= \frac{1}{n + 1} b_{n}. \end{align*}

Problem. How to compute multiplication, $$C(x) = A(x)B(x).$$

Multiplication is done by convolution, \begin{align*} \sum_{n \ge 0} c_n z^n &= \left( \sum_{k \ge 0} a_k z^k \right)\left( \sum_{k \ge 0} b_k z^k \right) \\ &= \sum_{n \ge 0} \left( \sum_{k = 0}^n a_k b_{n - k}\right) z^n \end{align*} That is, equating like powers, $$ c_n = \sum_{k = 0}^n a_k b_{n - k} $$ There are faster tricks but this one is fine, even for decent filtrations.

Problem. How to compute division, $$A(x) = \frac{C(x)}{B(x)}.$$

Take the convolution equation for the multiplication and solve for $a_n$, $$ a_n = \frac{1}{b_0} \left(c_n - \sum_{k = 0}^{n - 1} a_k b_{n - k}\right). $$ Please note that $B(x)$ is inversible if and only if $b_0 \neq 0$. Again, there are faster tricks but this is already a decent algorithm.

Problem. How to compute $A(x) = \log\left(B(x)\right)$ fast.

We've already noted that integration and differentiation of formal series can be done term by term in a very efficient way. The best idea I've found it to write, $$ A(x) = \int \frac{B'(x)}{B(x)} dx $$ The constant of integration is $$a_0 = \log(b_0).$$ The only costly part here is the division but this method is considerably faster than considering substituting $A(x)$ in the power series expansion of the logarithm.

Problem. How to compute $B(x) = \exp\left(A(x)\right)$ fast.

Consider the folowing differential equation, $$ \frac{d}{dx} B(x) = A'(x) B(x) $$ You can now transform it in a recurrence equation on the coefficients by equating terms of like power. You get, $$ (k+1) b_{k + 1} = a_1 b_{k} + 2 a_2 b_{k - 1} + 3 a_3 b_{k - 2} + ... $$ The recurrence is initialized by $$b_0 = \exp(a_0).$$ The above recurrence relation is valid for all $k \ge 0$ as long as you agree to consider $b_k = 0$ for all negative $k$.


Alfredo de la Fuente  posted 2013-09-05 23:46:53

As always, great tricks well explained.

Samuel Vidal  posted 2013-09-09 00:57:16

Thank you