On the functional square root of $x^2+1$ There are some math quizzes like:
find a function $\phi:\mathbb{R}\rightarrow\mathbb{R}$ 
such that $\phi(\phi(x)) = f(x) \equiv x^2 + 1.$

If such $\phi$ exists (it does in this example), $\phi$ can be viewed as a "square root" of $f$ in the sense of function composition because $\phi\circ\phi = f$.  Is there a general theory on the mathematical properties of this kind of square roots? (For instance, for what $f$ will a real analytic $\phi$ exist?)
 A: (Update:  Oh sorry; I'm new here. Didn't notice the much more substantial link to the mathoverflow. But perhaps it is interesting for some newbie anyway...) 
If the function is a polynomial (or powerseries) with constant term $\ne0$ this is difficult, but sometimes a meaningful solution can be given.
If the function is as above, but the constant term is zero, then it is just the question of relatively simple manipulation of the formal powerseries/polynomial to construct a "compositional" (or "iterative") root. 
There is a very good deal in L.Comtet "Advanced Combinatorics" at pages around 130..150 (don't have it around).
Also the keywords "Bell-matrix", "Carleman-matrix" are helpful: these matrices transform the problem of funtion composition/iteration to matrix-products/matrix-powers. (The matrix-notation just implies the required formal powerseries-manipulations) . This is well established for functions $f(x)= ax + bx^2 + cx^3 +\cdots$.
If the constant term does not vanish, $f(x) = K + ax + bx^2 + cx^3 +\cdots$, then things are usually much more difficult. But there is one -I think: again well established- workaround: 
Rewrite $f(x)$ as some function $g(x+x_0) - x_0$ such that now $g(x)$ has no constant term and then apply the above mentioned procedures (solving for iterative square root and the like) to $g(x)$.
Example: denote the iterative root of a some function $f(x)$ by $f^{[1/2]}(x)$ then solve
$$f^{[1/2]}(x) = g^{[1/2]}(x-x_0)+x_0$$
where you first must find $x_0$.
[update 2] I've tried to generate that power series for $g^{[1/2]}(x)$, which has then complex coefficients. I got 
$$ g^{[1/2]}(x) \approx  \sqrt{2 x_0} x + (0.20412415 - 0.22379688 i) x^2 + (0.050024298 + 0.048042380 i) x^3 + (-0.022112213 + 0.028383580 i) x^4 + (-0.023623808 - 0.010393981 i) x^5 + (0.00074679924 - 0.021136421 i) x^6 + O(x^7)$$
where $f^{[1/2]}(x) = g^{[1/2]}(x-x_0) + x_0 $ and the fixpoint is $x_0 = \exp(i \pi /3 ) \approx 1/2 + 0.86602540 i $ and $\sqrt{2 x_0}=(1.2247449 + 0.70710678 i) $ . The range of convergence is small; however using Euler-summation (of complex orders) I could manage to reproduce $f(0)$ to $f(1)$ by applying two times the half-iterative to about 6 to 8 correct digits. [end update 2] 
For a very basic introduction you may try my small essay at 

go.helms-net.de/math/tetdocs

and look for the article "continuous functional iteration".
A: To get a closed-form solution (where possible) you can find a flow of f(x).
Let's w(x) a flow of f(x).
Then to find it we have to solve a difference equation (called Abel equation):
$$w(t+1)=f(w(t))$$
In our case it's 
$$w(t+1)=1+w(t)^2$$
or in difference form,
$$\Delta w + w - w^2-1=0$$
Unfortunately this first-order difference equation is non-linear in our case.
Supposing you somehow solved it, you will get $w_C(t)$, a function depending on a constant parameter C. Then you take C=x and t=1/2 (for square iterative root), this will be the answer.
A: Using the same approach as in this answer, it is possible to produce the functional square root by using:
$$\phi_0(x)=|x|^{\sqrt2}$$
$$\phi_{n+1}(x)=f^{-1}(\phi_n(f(x)))=\sqrt{\phi_n(x^2+1)-1}$$
which converges uniformly to a function $\phi$ satisfying $\phi(\phi(x))=f(x)$, which can be seen by noting that
$$\phi_n(\phi_n(x))=f^{-n}(|f^n(f^{-n}(|f^n(x)|^{\sqrt2}))|^{\sqrt2})=f^{-n}(f^n(x)^2)$$
is nearly equal to
$$\phi_n(\phi_n(x))\simeq f^{-n}(f^n(x)^2+1)=f(x)$$
with the error decreasing with each application of $f^{-n}$ since $|f^{-1}(a+\epsilon)-f^{-1}(a)|\simeq\epsilon/\sqrt a\ll\epsilon/2$ implies that
$$|\phi_n(\phi_n(x))-f(x)|\ll1/2^n$$
and in fact the error tends to zero much faster, since initially $a=f^n(x)^2\gg f^2(x)^{2^{n-1}}$ is insanely large.

There is the additional caveat that technically this requires us to compute iterations of $f$ since we have
$$\phi_n(x)=f^{-n}(|f^n(x)|^{\sqrt2})$$
where $f^n(x)$ is far too large to compute. We note two things at this point:


*

*The convergence is at least quadratic, so one should not need to use very large $n$ anyways.

*If one does need to use larger $n$, it is better to use an improved $\phi_0$. Expanding on what is shown in the linked question, we can use the improved$$\phi_0(x)=|x|^{\sqrt2}+\frac1{\sqrt2}|x|^{\sqrt2-2}-\frac12|x|^{-\sqrt2}$$to get faster convergence. Note that we only care about the behavior as $x\to\infty$ since $f^n(x)$ is large.
A: Introduce a new coordinate system with a fixed point of $f$ as origin, e.g. the point $\omega:=e^{\pi i/3}$. Writing $x=\omega+\xi$ with a new independent variable $\xi$  one has $\phi(\omega+\xi)=\omega +\psi(\xi)$ for a new unknown function $\psi$ with $\psi(0)=0$. This function $\psi$ satisfies in a neighbourhood of $\xi=0$ the functional equation $\psi(\psi(\xi))=2\omega\xi+\xi^2$. Now you can recursively determine the Taylor coefficients of $\psi$. If you are lucky the resulting formal series actually defines an analytical function in a neighbourhood of $\xi=0$.
A: Look also at this answer:
https://mathoverflow.net/questions/17605/how-to-solve-ffx-cosx/44727#44727
In short, the analytic solution is
$$f^{[1/2]}(x)=\phi(x)=\sum_{m=0}^{\infty} \binom {1/2}m \sum_{k=0}^m\binom mk(-1)^{m-k}f^{[k]}(x)$$
$$f^{[1/2]}(x)=\lim_{n\to\infty}\binom {1/2}n\sum_{k=0}^n\frac{1/2-n}{1/2-k}\binom nk(-1)^{n-k}f^{[k]}(x)$$
$$f^{[1/2]}(x)=\lim_{n\to\infty}\frac{\sum_{k=0}^{n} \frac{(-1)^k f^{[k]}(x)}{(1/2-k)k!(n-k)!}}{\sum_{k=0}^{n} \frac{(-1)^k }{(1/2-k) k!(n-k)!}}$$
The same way you can find not only square iterative root but iterative root of any power.
