How to compute the Fourier coefficients of $\exp(f(x))$ given the Fourier coefficients of $f(x)$ In an attempt to solve a certain differential equation numerically with the Galerkin method I bumped into this problem:
Given the real function $f(x) = \sum_{n = -N}^{N} a_n \, \exp(inx) $ compute the Fourier coefficients with indices between $-N$ and $N$ of the function $\exp(f(x))$.
Of course the coefficients can be computed to any required accuracy by brute force numerical methods such as Gaussian Quadrature. In terms of computation time the $O(N^2)$ operations needed by the Gaussian quadrature are however to expensive in this case. What I am looking for is therefore a closed formula or algorithm for the coefficients such that the $2N+1$ coefficients can be computed with $O(N \log N)$ operations or less.
Here is what I have tried:
Let $b_p$ denote the $p$'th Fourier coefficient of $\exp(f(x))$. For notational convinience we will consider the coefficient $b_{-p}$, which can be calculated as
$$b_{-p} = \frac{1}{2 \pi} \int_{0}^{2 \pi} \exp(f(x)) \exp(ipx) \, dx. $$
Inserting the expansion for $f$ into this expression yields
$$b_{-p} = \frac{1}{2 \pi} \int_{0}^{2 \pi} \exp\bigg( \sum_{n = -N}^{N} a_n \exp(inx) \bigg) \exp(ipx) \, dx. $$
We now rewrite this as a contour integral along the unit circle in the complex plane. One parametrization of the unit circle is $\gamma(x) = \exp(ix)$ where $0 \le x \le 2 \pi$ and by definition of a complex contour integral we therefore have 
$$b_{-p} = - \frac{i}{2 \pi} \int_{\gamma} \exp\bigg( \sum_{n = -N}^{N} a_n z^n \bigg)  \, z^{p-1} dz. $$
By the residue theorem of Cauchy this integral is easily calculated if the residues of the function $\exp\bigg( \sum_{n = -N}^{N} a_n z^n \bigg)  \, z^{p-1}$ are known but this is where I get stuck. To get the residues I have tried to rewrite the exponential as
$$ \exp\bigg( \sum_{n = -N}^{N} a_n z^n \bigg)  = \exp(a_{-N} z^{-N}) \exp(a_{-N+1} z^{-N+1}) \cdots \exp(a_{N-1} z^{N-1}) \exp(a_{N} z^{N}) $$
and then expanded each of the exponential functions as a Taylor series (the positive powers of $z$) or a Laurent series (the negative powers of $z$). This unfortunately creates a mess and I have not been able to extract a useful formula for $b_{-p}$ from it.
 A: Maybe this helps:  If $\exp\bigl(f(x)\bigr)=: g(x)$ then $g'(x)=f'(x)\>g(x)$. This immediately translates into
$$b_n={1\over n}\sum_{k+l=n} k\, a_k \,b_l\ .$$
A: In case anybody should be interested, here is my solution to the problem.


Given the Fourier coefficients of $f(x)$ we can compute the Fourier coefficients of $f(x)^n$ for $n = 0, 1, ..., K$ in $O(N log(N))$ operations (if $K$ is negligible compared to $N$) using the fast Fourier transform. We will in the following let $a_p^{(n)}$ denote the $p$'th Fourier coefficient of $f(x)^n$.


Let us now consider the $-p$'th Fourier coefficient of $\exp(f(x))$ given by 
$$ b_{-p} = \frac{1}{2 \pi} \int_{0}^{2\pi} \exp(f(x)) \exp(ipx) dx .$$
Taylor expanding the exponential function and taking the integral inside the sum gives
$$ b_{-p} = \sum_{n = 0}^{\infty} \frac{1}{n!} \frac{1}{2 \pi} \int_{0}^{2\pi} f(x)^n \exp(ipx) dx $$
The integral divided by $2 \pi$ is just the $-p$'th Fourier coefficient of $f(x)^n$ meaning that we can express $b_{-p}$ as
$$ b_{-p} = \sum_{n = 0}^{\infty} \frac{1}{n!} a_{-p}^{(n)}$$
To obtain an approximation for $b_{-p}$ I truncate this series after $K$ terms.
In my case it turns out that the series is rapidly convergent and truncating it after $K = 5$ terms gives an accurate enough solution.
