Finding f such that f(f(x))=g(x) given g Suppose $g(x)$ is a smooth increasing function defined for $x \ge 0$ such that $g(x) \ge x$ for all $x$. Does there exist a function $f$ with similar properties such that $f(f(x))=g(x)$ for all $x \ge 0$? (You can interpret "similar" as widely as you'd like - smoothness would be great, but even continuity would be nice.)
I asked the question given these conditions on $g$ since it seems reasonable that they would produce a positive answer. However, I'm just as interested in the same question for more general classes of $g$. For example, suppose we only assume $g$ is continuous, or even measurable - can we find an $f$ with the same properties? And let's suppose we relax the requirement $g(x) \ge x$, etc (I included that because it helps ensure the existence of a set-theoretic $f$).
Under the given conditions, how many such $f$ exist?
 A: Assume that $f(x)$ and $g(x)$ are analytic and that $f$ has a fixed point, without loss of generality $f(0)=0$ and therefore $f(f(0))=g(0)=0$. The Taylors series of $g^n(x)$ can be constructed.
 using Bell polynomials of iterated functions and then the iterate can be evaluated at $1/2$. 
Taylors series of iterated functions
Let
$$f(x)=g^{1/2}(x)$$
then
\begin{eqnarray*}
g^n(x) &=& g_1^n x+ \frac{1}{2} g_2(\sum_{k_1=0}^{n-1}g_1^{2n-k_1-2}) x^2  \\
&+&             \frac{1}{6}[g_3(\sum_{k_1=0}^{n-1}g_1^{3n-2k_1-3})
               +3g_2^2 (\sum_{k_1=0}^{n-1}
                       \sum_{k_2=0}^{n-k_1-2}
                          g_1^{3n-2k_1-k_2-5}) ] x^3 + \cdots
\end{eqnarray*}  
Proof.
Let $h(x)=g^{n-1}(x)$, then
\begin{eqnarray*}
 Dg(h(z))&=&g'(h(x))h'(x)\\
         &=&g'(g^{n-1}(x))Dg^{n-1}(x)\\
         &=&\prod^{n-1}_{k_1=0}g'(g^{n-k_1-1}(x))\\
\end{eqnarray*}
\begin{eqnarray}           
 Dg^n(0)&=&g'(0)^n \\
        &=&g_1^n = g_1^{1/2}
\end{eqnarray}
For the second derivative 
\begin{eqnarray}            
 D^2g^n(0)&=&g_2g_1^{2n-2}+g_1 D^2g^{n-1}(0) \\
            &=&g_1^0g_2 g_1^{2n-2} \\
            &&+g_1^1g_2 g_1^{2n-4} \\
            &&+\cdots \\
            &&+g_1^{n-2}g_2 g_1^2 \\
            &&+g_1^{n-1}g_2 g_1^0 \\
            &=&g_2\sum_{k_1=0}^{n-1}g_1^{2n-k_1-2} 
\end{eqnarray}
For the third derivative 
\begin{eqnarray}
 D^3g^n(0)&=&g_3g_1^{3n-3}+3 g_2^2\sum_{k_1=0}^{n-1}g_1^{3n-k_1-5}
               +g_1 D^3g^{n-1}(0) \\
            &=&g_3\sum_{k_1=0}^{n-1}g_1^{3n-2k_1-3}
               +3g_2^2 \sum_{k_1=0}^{n-1}
                       \sum_{k_2=0}^{n-k_1-2}
                          g_1^{3n-2k_1-k_2-5}                                                                           
\end{eqnarray}
Note that the index $k_1$ from the second derivative is renamed $k_2$ in the final summation of the third derivative. A certain amount of renumbering is unavoidable in order to use a simple index scheme.
$\blacksquare$
Algebraic expansion shows that 
$$g^m(g^n(x))=g^{m+n}(x)$$.
Classification of Fixed Points
Note that $g^n(x)$ was not evaluated at $n=1/2$. What at first may appear to be a defect is a strength as it reveals the impact of the classification of fixed points. If $g'(0)$ equals a root of unity, then the geometrical progressions evaluate very differently than the standard method of evaluation. This formula works regardless of whether dealing with  hyperbolic fixed points, parabolically neutral fixed points and rationally and irrationally neutral fixed points.
See Bell polynomials of iterated functions and What are efficient ways to compute the derivatives of iterated functions? for more details.
A: My previous answer gives you a hint at analytic solution, but one may ask how to find the result in closed form.
To do so, one has to find family of flows $\phi_C(t)$ of the function $g(x)$ from the condition $\phi_C(t+1)=g(\phi_C(t))$. Then evaluate $\phi_C(1/2)$ and finally substitute $x$ instead of $C$.
For example, for the problem $f^{[2]}(x)=x^2$ the flow equation will be $\phi_C(t+1)=\phi_C(t)^2$ Taking logarithm with $C$ base one arrives at $\Delta\log_C(\phi_C(t))=\log_C\phi_C(t)$. Substituting $p(t)=\log_C\phi_C(t)$ we get a first-order linear difference equation $\Delta p(t)=p(t)$. Since $2^t$ is the only non-trivial function known to satisfy this equation, we get $\log_C\phi_C(t)=2^t$, and so $\phi_C(t)=C^{2^t}$. Evaluating at $1/2$ and substituting $x$ for $C$ we get $f(x)=x^{\sqrt{2}}$.
Of course, in many cases the equation may be non-linear and difficult to solve.
In general, any iterative equation of order $n$ (that is involving $f^{[n]}(x)$) is reducible to difference equation of order $n-1$ over flows. Such equation is expected to have no more than $n-1$-dimensional space of solutions, except maybe some trivial cases so the original equation is also in non-trivial cases expected to have no more than $n-1$ solutions. Since your problem is iterative equation of order $2$, the solution is expected unique if exists, except some trivial cases such as $f^{[2]}(x)=x$ and $f^{[2]}(x)=-x$ where any 1-periodic and 1-antiperiodic functions (respectively) will satisfy as a flow.
A: As is typically done in iterated systems, let $f(0)=0$ where $f(z)$ is smooth. Then applying Riordan's approach to solving inverse Bell polynomials produces the derivatives of $f(z)$ which are used for the power series $f(z)$ giving $$f(z)=\sum_{k=1}^{\infty}\frac{1}{k!}f^{(k)}(0) z^k.$$ There are two solutions based on whether $f'(0) = -\sqrt{g'(0)}$ or $f'(0) = \sqrt{g'(0)}$. Due to space, only the first three derivitives are shown.
Mathematica Code
f[0]=0;
max = 3;
Solve[Table[D[g[z] == f[f[z]], {z, i}] /. z -> 0 , {i, max}],
      Table[D[f[z], {z, i}] /. z -> 0, {i, max}]]

Solution 1
$$f'(0) = -\sqrt{g'(0)}$$
$$f''(0) = \frac{g'(0) g''(0)+\sqrt{g'(0)} g''(0)}{g'(0)^2-g'(0)} $$ 
$$f'''(0) = -\frac{3 g''(0)^2+g^{(3)}(0) \sqrt{g'(0)} \left(\sqrt{g'(0)}-1\right)^2}{\left(\sqrt{g'(0)}-1\right)^2 g'(0) \left(g'(0)+1\right)}$$

Solution 2
$$f'(0) = \sqrt{g'(0)} $$
$$f''(0) = \frac{g'(0) g''(0)-\sqrt{g'(0)} g''(0)}{g'(0)^2-g'(0)} $$ 
$$f'''(0) = \frac{g^{(3)}(0) \left(\sqrt{g'(0)}+1\right)^2 \sqrt{g'(0)}-3 g''(0)^2}{\left(\sqrt{g'(0)}+1\right)^2 g'(0) \left(g'(0)+1\right)}$$
Riordan, John, Combinatorial identities, New York- London-Sydney: John Wiley and Sons, Inc. 1968. XII, 256 p. (1968). ZBL0194.00502.
A: Assuming that $g(x) > x$ rather than $g(x) \geq x$ for all $x$, and also that $g$ is strictly increasing (no critical points), one can obtain a solution $f$ to $g = f \circ f$ which piecewise has roughly similar properties to $g$, in particular $f$ will also be strictly increasing with $f(x)>x$.  Things look to be more interesting if one relaxes these hypotheses.
Here's the construction.  Firstly, define $x_n$ iteratively for $n=0,1,\dots$ by $x_0 := 0$, $x_{n+1} := g(x_n)$, then the $x_n$ will be increasing to infinity (they cannot accumulate at any finite point $x_*$ as one would then have $g(x_*)=x_*$).  
Now also define $x_n$ for half-integers $n=1/2,3/2,\dots$ by picking $x_{1/2}$ arbitrarily between $x_0,x_1$ and then setting $x_{n+1} := g(x_n)$ for $n=1/2,3/2,\dots$.  The $x_n$ for half-integer $n$ interlace between the $x_n$ for integer $n$, so the $x_n$ still increase to infinity as $n$ ranges over the combined index set $0,1/2,1,3/2,\dots$.
For any $n$ in this combined index set, the function $g$ is a continuous increasing map from $[x_n,x_{n+1/2}]$ to $[x_{n+1},x_{n+3/2}]$ that maps endpoints to endpoints, and must therefore be a bijection (by the intermediate value theorem) and thus a homeomorphism (as both domain and range are compact Hausdorff).  Call $g_n: [x_n,x_{n+1/2}] \to [x_{n+1}, x_{n+3/2}]$ the restriction of $g$ to these intervals.  To finish the job it will suffice to find increasing homeomorphisms $f_n:[x_n,x_{n+1/2}] \to [x_{n+1/2},x_{n+1}]$ mapping endpoints to endpoints such that $g_n = f_{n+1/2} \circ f_n$ for all $n=0,1/2,\dots$.  This has "decoupled" the two factors of $f$ in the original equation $g = f \circ f$ and it is now easy to describe the general solution to this by writing the equation as $f_{n+1/2} = g_n \circ f_n^{-1}$ and solving recursively.  More explicitly: pick an arbitrary increasing homeomorphism $f_0: [x_0,x_{1/2}] \to [x_{1/2},x_1]$ mapping endpoints to endpoints, and set $f_n := g^n \circ f_0 \circ g^{-n}$ for integer $n$ and $f_n := g^{n+1/2} \circ f_0^{-1} \circ g^{1/2-n}$ for half-integer $n$.
The function $f$ produced here by gluing together the $f_n$ will be continuous and strictly increasing and obey the required equation $g = f \circ f$ (in fact this is the general solution to this equation with the stated properties, bearing in mind that $x_{1/2}$ is also arbitrary between $x_0$ and $x_1$).  If $g$ has no critical points then the $f_n$ will be diffeomorphisms and thence the $f_n$ will also.  So $f$ will be smooth except possibly at the transition points $x_n$.  Actually one can fix things up to be smooth at the endpoints also: if one extends $f_0$ a bit past $x_{1/2}$ to an extension $\tilde f_0$ that stays smooth near $x_{1/2}$ without critical points, and then modifies $f_0(x)$ smoothly near $x_0$ to equal $\tilde f_0^{-1}(g(x))$ for $x$ sufficiently close to $x_0$, then one can check that $f$ is now smooth everywhere.  The situation becomes more interestingly complicated when $g$ has critical points $g'(x)=0$ or fixed points $g(x)=x$, but I haven't looked into this carefully (presumably the analysis of the analytic case in other answers will indicate what the behaviour should be there).  
A: I do not know for "the same properties", but Newton series does the trick for a large class of functions. 
$$f(x)=g^{[1/2]}(x)=\sum_{m=0}^{\infty} \binom {1/2}m \sum_{k=0}^m\binom mk(-1)^{m-k}g^{[k]}(x)$$
A more extended answer you can find following the link 
A: Different conditions that may still be of interest:

If $g$ is a strictly increasing continuous function on $\mathbb R$ and $g(x)>x$ for all $x$, then there exists uncountably many families of strictly increasing continuous functions $f_\alpha$ indexed by reals satisfying $f_1=g$ and $f_{\alpha+\beta}=f_\alpha\circ f_\beta$.

The construction is not so bad. First pick a point $x_0$ and define $f_\star(\alpha)=f_\alpha(x_0)$ so that $f_\star(0)=x_0$ and $f_\star(\alpha+1)=g(f_\star(\alpha))$. Using the latter equation, we have $f_\star(1)=g(x_0)$. From here, we may define $f_\star$ on $[0,1]$ as any strictly increasing continuous function matching the provided bounds. Due to the functional equation, this furthermore defines $f_\star$ on $(a,\infty)$ for some $a\le0$ (the inverse of $g$ may not be defined on all of $\mathbb R$).
Now that we have our $f_\star$, observe that it is strictly increasing, continuous, and has range $\mathbb R$, so $f_\alpha(x)=f_\alpha(f_{f_\star^{-1}(x)}(x_0))=f_{\alpha+f_\star^{-1}(x)}(x_0)=f_\star(\alpha+f_\star^{-1}(x))$ is defined everywhere that has $\alpha+f_\star^{-1}(x)>a$. Particularly we have $f_{0.5}(x)=f_\star(0.5+f_\star^{-1}(x))$ and $f_{0.5}\circ f_{0.5}=g$.
So the only remaining variable in this construction is how we choose to define $f_\star$ on $[0,1]$, which can be manipulated to try and satisfy additional criteria.
