Café Math : Basic Power Series Algorithms

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.

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$.

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

Thank you

As always, great tricks well explained.