A formal power series such that $f(f(x))=x$ Let
$$f(x)=\sum_{n=1}^\infty a_n x^n$$
be a formal power series with no constant term such that $f(f(x))=x$. We find that
$$f(f(x))=a_1^2x+(a_1a_2+a_1^2a_2)x^2+(a_1a_3+2a_1a_2^2+a_3a_3)x^3+\dots$$
so $a_1^2=1$. If $a_1=1$, you need all the other terms to be zero, however if $a_1=-1$ we get a family of nontrivial solutions. Let $a_2=a$, and requiring the higher coefficients of $f(f(x))$ to be zero we can find $a_3=2a^2$, $a_4=-\frac{11}2a^3$, $a_5=\frac{11}2a^4$, $a_6=\frac{105}4a^5$, $a_7=-\frac{279}2 a^6$...
Is there a closed form for these numbers?
 A: Warning: this answer is wrong. See Qiaochu Yuan's comment and answer.
The coefficients are uniquely determined recursively, hence they are those of the involutive series
$$-\frac x{1+ax}=-x+ax^2-a^2x^3+a^3x^4-\dots.$$
A: A result of Marshall Cohen posted in the arxiv answers the question. Their result is:
If $n$ is a positive integer and $\omega$ is a primitive $n$-th root of unity in the field $\mathbb F$, then for every infinite sequence $\{a_n\}_{1<k\neq nj+1}$ of elements of $\mathbb F$ there exists a unique sequence $\{a_{nj+1}\}_{n\in\mathbb N}$ such that the formal power series
$$f(z)=\omega z+\sum_{h=2}^\infty a_h z^h$$
has order $n$ as a power series over $\mathbb F$.
The question asked for the case $n=2$ and $\mathbb F=\mathbb C$, and made a specific choice for the $\{a_n\}_{1<k\neq nj+1}$.
A: I don't think your approach will work. If $f$ is subject to $f \circ f = \mathrm{id}$, then every bijective map $g$ induces another solution $h := g^{-1} \circ f \circ g$:
$$h \circ h = g^{-1} \circ f \circ g \circ g^{-1} \circ f \circ g = g^{-1} \circ f \circ f \circ g = g^{-1} \circ g = \mathrm{id}.$$
A: The claim that the coefficients are unique given the quadratic term is incorrect. Before going into that, here's some background on Anne's answer. We consider the simpler question of how to determine solutions to $f(f(x)) = x$ where $f$ is a Mobius transformation $f(x) = \frac{ax + b}{cx + d}$. Over a field $K$ the group of Mobius transformations is isomorphic to the projective general linear group $PGL_2(K)$, with the isomorphism given by sending a $2 \times 2$ matrix $\left[ \begin{array}{cc} a & b \\ c & d \end{array} \right]$ to $\frac{ax + b}{cx + d}$. So the problem reduces to finding matrices $M$ squaring to a scalar multiple of the identity.
So write $M^2 = c$ for some scalar $c$, so that the eigenvalues of $M$ are $\pm \sqrt{c}$ (over an algebraically closed field, say $\mathbb{C}$). Since we only need to work up to scale we might as well divide $M$ by $c$, or equivalently assume WLOG that $c = 1$, so $M^2 = 1$ and $M$ has eigenvalues $\pm 1$. If $M$ is conjugate to a nontrivial Jordan block then it squares to another nontrivial Jordan block so can't square to $1$; hence $M$ is diagonalizable. If $M$ has eigenvalues $\{ 1, 1 \}$ or $\{ -1, -1 \}$ then it's a scalar multiple of the identity; otherwise its eigenvalues are $\{ 1, -1 \}$. This implies that $\text{tr}(M) = 0$ and $\det(M) = -1$, so $M$ has the form $\left[ \begin{array}{cc} a & b \\ c & -a \end{array} \right]$ where $\det(M) = -a^2 - bc = -1$. This gives $bc = 1 - a^2$.
Now we add the constraint that as a formal power series $f(x)$ should have no constant term. This means $f(0) = 0$ which gives $b = 0$. Then $a = \pm 1$ and we can take $a = 1$ WLOG, which gives $M = \left[ \begin{array}{cc} 1 & 0 \\ c & -1 \end{array} \right]$, hence
$$f(x) = \frac{x}{cx - 1} = -x - cx^2 - c^2 x^3 - \dots $$
as in Anne's answer (with $a = -c$).

On the other hand Klaus's answer implies that the resulting power series cannot be unique, since we can conjugate by an invertible (with respect to composition) formal power series. With a little more effort (showing that we can even arrange for the first $n$ coefficients of this series past $x$ to vanish) we can show that the first $n$ coefficients never uniquely determine the others. So you did something funny with your calculations but I'm not sure what.
Generally, it's known that every formal power series of finite order (with respect to composition) is conjugate to $f(x) = \zeta x$ for $\zeta$ some root of unity (above we have $\zeta = -1$, and in general we probably need to work over an algebraically closed field). This implies:

Classification: If $f(x)$ is a formal power series satisfying $f(0) = 0$ and $f(f(x)) = x$, then either $f(x) = x$ or $\boxed{ f(x) = g(-g^{-1}(x)) }$ where $g(x) = x + \dots$.

Given $g(x)$, the coefficients of $g^{-1}(x)$ can be computed using Lagrange inversion. This gives solutions depending on an infinite number of parameters, namely the higher coefficients $g_i$ of $g(x)$, and the first $n$ coefficients of $f$ depend only on the first $n$ coefficients of $g$. For example, if $g(x) = x + g_2 x^2 + \dots$ then $g^{-1}(x) = x - g_2 x^2 + \dots$ which gives $f(x) = -x + 2 g_2 x^2 + \dots$.
To give a relatively explicit example, take $g(x) = x - ax^2$. Then $g^{-1}(x) = \frac{1 - \sqrt{1 - 4ax}}{2a}$ by the quadratic formula (we have to take the minus sign so that $g^{-1}(0) = 0$); this is a version of the generating function of the Catalan numbers. If $y - ay^2 = x$ then $g(-y) = -y - ay^2 = x - 2y$, which gives
$$\boxed{ \begin{align*} f(x) &= x - \frac{1 - \sqrt{1 - 4ax}}{a}  \\
 &= -x - 2ax^2 - 4a^2 x^3 - 10a^3 x^4 - \dots \\
 &= -x - \sum_{n=2}^{\infty} \frac{2}{n} {2n-2 \choose n-1} a^{n-1} x^n. \end{align*} }$$
