half iterate of $x^2+c$

This may be helpful.

Let $$ f(x) = \frac{-1 + \sqrt{1 + 4 x}}{2}, \; \; x > 0 $$ We use a technique of Ecalle to solve for the Fatou coordinate $\alpha$ that solves $$ \alpha(f(x)) = \alpha(x) + 1. $$ For any $x > 0,$ let $x_0 = x, \; x_1 = f(x), \; x_2 = f(f(x)), \; x_{n+1} = f(x_n).$ Then we get the exact $$ \alpha(x) = \lim_{n \rightarrow \infty} \frac{1}{x_n} - \log x_n + \frac{x_n}{2} - \frac{x_n^2}{3} + \frac{13 x_n^3}{36} - \frac{113 x_n^4}{ 240} + \frac{1187 x_n^5}{ 1800} - \frac{877 x_n^6}{ 945} - n. $$ The point is that this expression converges far more rapidly than one would expect, and we may stop at a fairly small $n.$ It is fast enough that we may reasonably expect to solve numerically for $\alpha^{-1}(x).$

We have $$ f^{-1}(x) = x + x^2. $$ Note $$ \alpha(x) = \alpha(f^{-1}(x)) + 1, $$ $$ \alpha(x) - 1 = \alpha(f^{-1}(x)) , $$ $$ \alpha^{-1} \left( \alpha(x) - 1 \right) = f^{-1}(x). $$ It follows that if we define $$ g(x) = \alpha^{-1} \left( \alpha(x) - \frac{1}{2} \right), $$ we get the miraculous $$ g(g(x)) = \alpha^{-1} \left( \alpha(x) - 1 \right) = f^{-1}(x) = x + x^2. $$

I put quite a number of relevant pdfs at BAKER. The host computer for this was down for about a year but has recently been repaired.

EDIT, TUESDAY:

Note that $\alpha$ is actually holomorphic in an open sector that does not include the origin, such as real part positive. That is the punchline here, $\alpha$ cannot be extended around the origin as single-valued holomorphic. So, since we are finding a power series around $0,$ not only are there a $1/z$ term, which would not be so bad, but there is also a $\log z$ term. So the $\ldots -n$ business is crucial.

I give a complete worked example at my question https://mathoverflow.net/questions/45608/formal-power-series-convergence as my answer https://mathoverflow.net/questions/45608/formal-power-series-convergence/46765#46765

The Ecalle technique is described in English in a book, see K_C_G PDF or go to BAKER and click on K_C_G_book_excerpts.pdf The Julia equation is Theorem 8.5.1 on page 346 of KCG. It would be no problem to produce, say, 50 terms of $\alpha(x)$ with some other computer algebra system that allows longer power series and enough programming that the finding of the correct coefficients, which i did one at a time, can be automated. No matter what, you always get the $\alpha = \mbox{stuff} - n$ when $f \leq x.$

As I said in comment, the way to improve this is to take a few dozen terms in the expansion of $\alpha(x)$ so as to get the desired decimal precision with a more reasonable number of evaluations of $f(x).$ So here is a brief version of the GP-PARI session that produced $\alpha(x):$

=======

    ? taylor( (-1 + sqrt(1 + 4 * x))/2  , x  )
    %1 = x - x^2 + 2*x^3 - 5*x^4 + 14*x^5 - 42*x^6 + 132*x^7 - 429*x^8 + 1430*x^9 - 4862*x^10 + 16796*x^11 - 58786*x^12 + 208012*x^13 - 742900*x^14 + 2674440*x^15 + O(x^16) 

    f = x - x^2 + 2*x^3 - 5*x^4 + 14*x^5 - 42*x^6 + 132*x^7 - 429*x^8 + 1430*x^9 - 4862*x^10 + 16796*x^11 - 58786*x^12 + 208012*x^13 - 742900*x^14 + 2674440*x^15  

    ? fp = deriv(f) 
    %3 = 40116600*x^14 - 10400600*x^13 + 2704156*x^12 - 705432*x^11 + 184756*x^10 - 48620*x^9 + 12870*x^8 - 3432*x^7 + 924*x^6 - 252*x^5 + 70*x^4 - 20*x^3 + 6*x^2 - 2*x + 1 

    L = - f^2 + a * f^3 

    R = - x^2 + a * x^3

    compare = L - fp * R 

    19129277941464384000*a*x^45 - 15941064951220320000*a*x^44 +
 8891571783902889600*a*x^43 - 4151151429711140800*a*x^42 + 
1752764158206050880*a*x^41 - 694541260905326880*a*x^40 + 
263750697873178528*a*x^39 - 97281246609064752*a*x^38 + 35183136631942128*a*x^37 
- 12571609170862072*a*x^36 + 4469001402841488*a*x^35 - 1592851713897816*a*x^34 + 
575848308018344*a*x^33 - 216669955210116*a*x^32 + 96991182256584*a*x^31 + 
(-37103739145436*a - 7152629313600)*x^30 + (13153650384828*a + 
3973682952000)*x^29 + (-4464728141142*a - 1664531636560)*x^28 + (1475471500748*a 
+ 623503489280)*x^27 + (-479514623058*a - 220453019424)*x^26 + (154294360974*a + 
75418138224)*x^25 + (-49409606805*a - 25316190900)*x^24 + (15816469500*a + 
8416811520)*x^23 + (-5083280370*a - 2792115360)*x^22 + (1648523850*a + 
930705120)*x^21 + (-543121425*a - 314317080)*x^20 + (183751830*a + 
108854400)*x^19 + (-65202585*a - 39539760)*x^18 + (-14453775*a + 15967980)*x^17 
+ (3380195*a + 30421755)*x^16 + (-772616*a - 7726160)*x^15 + (170544*a + 
1961256)*x^14 + (-35530*a - 497420)*x^13 + (6630*a + 125970)*x^12 + (-936*a - 
31824)*x^11 + 8008*x^10 + (77*a - 2002)*x^9 + (-45*a + 495)*x^8 + (20*a - 
120)*x^7 + (-8*a + 28)*x^6 + (3*a - 6)*x^5 + (-a + 1)*x^4 

    Therefore a = 1  !!! 

    ? 
    L = - f^2 +  f^3 + a * f^4

    R = - x^2 +  x^3 + a * x^4 

    compare = L - fp * R 
     ....+ (1078*a + 8008)*x^10 + (-320*a - 1925)*x^9 + (95*a + 450)*x^8 + (-28*a - 100)*x^7 + (8*a + 20)*x^6 + (-2*a - 3)*x^5 

    This time a = -3/2  !


    L = - f^2 +  f^3  - 3 * f^4 / 2  + c * f^5 

    R = - x^2 +  x^3 - 3 * x^4 / 2  + c * x^5  

     compare = L - fp * R
    ...+ (2716*c - 27300)*x^11 + (-749*c + 6391)*x^10 + (205*c - 1445)*x^9 + (-55*c + 615/2)*x^8 + (14*c - 58)*x^7 + (-3*c + 8)*x^6 


    So c = 8/3 . 

    The printouts began to get too long, so I said no using semicolons, and requested coefficients one at a time..

    L = - f^2 +  f^3  - 3 * f^4 / 2  + 8 * f^5 / 3 + a * f^6; 

    R = - x^2 +  x^3 - 3 * x^4 / 2  + 8 * x^5 / 3  + a * x^6; 

       compare = L - fp * R;

    ? polcoeff(compare,5)
    %22 = 0
    ? 
    ?  polcoeff(compare,6)
    %23 = 0
    ? 
    ?  polcoeff(compare,7)
    %24 = -4*a - 62/3

    So this a = -31/6 


    I ran out of energy about here:
      L = - f^2 +  f^3  - 3 * f^4 / 2  + 8 * f^5 / 3 - 31 * f^6 / 6 + 157 * f^7 / 15 - 649 * f^8 / 30 + 9427 * f^9 / 210 + b * f^10 ; 

      R = - x^2 +  x^3 - 3 * x^4 / 2  + 8 * x^5 / 3  - 31 * x^6 / 6 + 157 * x^7 / 15 - 649 * x^8 / 30 + 9427 * x^9 / 210  + b * x^10;

       compare = L - fp * R; 
    ? 
    ?  polcoeff(compare, 10 )
    %56 = 0
    ? 
    ? 
    ?  polcoeff(compare, 11 ) 
    %57 = -8*b - 77692/105
    ? 
    ? 
      L = - f^2 +  f^3  - 3 * f^4 / 2  + 8 * f^5 / 3 - 31 * f^6 / 6 + 157 * f^7 / 15 - 649 * f^8 / 30 + 9427 * f^9 / 210 - 19423 * f^10 / 210 ; 

      R = - x^2 +  x^3 - 3 * x^4 / 2  + 8 * x^5 / 3  - 31 * x^6 / 6 + 157 * x^7 / 15 - 649 * x^8 / 30 + 9427 * x^9 / 210 - 19423 * x^10 / 210;

       compare = L - fp * R; 
    ?  polcoeff(compare, 10 )
    %61 = 0
    ? 
    ?  polcoeff(compare, 11 ) 
    %62 = 0
    ? 
    ?  polcoeff(compare, 12) 
    %63 = 59184/35
    ? 

    So R = 1 / alpha' solves the Julia equation   R(f(x)) = f'(x) R(x).

    Reciprocal is alpha'

    ? S =   taylor( 1 / R, x)
    %65 = -x^-2 - x^-1 + 1/2 - 2/3*x + 13/12*x^2 - 113/60*x^3 + 1187/360*x^4 - 1754/315*x^5 + 14569/1680*x^6 + 532963/3024*x^7 + 1819157/151200*x^8 - 70379/4725*x^9 + 10093847/129600*x^10 - 222131137/907200*x^11 + 8110731527/12700800*x^12 - 8882574457/5953500*x^13 + 24791394983/7776000*x^14 - 113022877691/18144000*x^15 + O(x^16) 

    The bad news is that Pari refuses to integrate 1/x, 
even when I took out that term it put it all on a common denominator,
 so i integrated one term at a time to get

alpha = integral(S)

and i had to type in the terms myself, especially the log(x)

    ? alpha = 1 / x - log(x) + x / 2 - x^2 / 3 + 13 * x^3 / 36 - 113 * x^4 / 240 + 1187 * x^5 / 1800 - 877 * x^6 / 945 + 14569 * x^7 / 11760 + 532963 * x^8 / 24192 

======

Arguably an off-subject remark:

If only you relented to allow c < 0, there is the celebrated ("chaotic " logistic map) closed form example (p302) of Ernst Schroeder himself (1870); namely, for
$$ h(x)= x^2-2, $$ it follows directly that for $$ y=\frac{x\pm \sqrt{x^2-4}}{2} $$ that is $$ x=y+y^{-1}, $$ one has $$ h(x)=y^2+y^{-2}\equiv h_1(x). $$ Whence, subscripting the iteration index, $$ h_n(x)= y^{2^n}+ y^{-2^n}. $$

This, then, specifies the whole iteration group: so your functional square root is just $$ h_{\sqrt2} (x)=y^{\sqrt 2} +y^{-\sqrt 2}. $$

Pardon if the point has been made, explicitly, or implicitly, in the outstanding answers above. If not, it might well offer guidance or continuation ideas.

More formally, in E.S.'s language of conjugacy, $\psi(x)=\frac{x\pm \sqrt{x^2-4}}{2}$, $~f(y)=y^2$, $~f_n(y)=y^{2^n}$; so that $h(x)= \psi^{-1} \circ f \circ \psi (x)$, and $$h_n= \psi^{-1} \circ f_n \circ \psi ~.$$

A conjugacy iteration approximation method is available in our 2011 paper: Approximate solutions of Functional equations. Apologies if this late lark answer is only proffering coals to Newcastle, but, in my experience, this is the canonical gambit of chaos discussions--naturally, domains and ranges are chosen suitably for the answer to make sense.

a plug

For some material on fractional iterates of $x^2+c$ see the last section of...
"Fractional Iteration of Series and Transseries" by G. A. Edgar ... LINK
To appear in Trans. Amer. Math. Soc.