# Square root of a function (in the sense of composition)

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?)

-
–  Samuel Aug 30 '10 at 11:30
Thanks. The linked posting is very useful. –  user1551 Aug 30 '10 at 11:51
This has also been asked on stackoverflow, but of course, those solutions are more program-oriented. –  KennyTM Sep 4 '10 at 15:16
–  Qiaochu Yuan Dec 3 '10 at 20:50

http://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.

-
I don't see, how this can effectively be used. I tried it with $x=0.5$ for some small n; "small" because you cannot iterate the function much without getting an overflow. See my answer for a short list of results... –  Gottfried Helms Aug 4 '13 at 20:26

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$.

-

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.

-
It seems to me that you can simply show the existence of a solution of a difference equation. The question here,it seems, would be of uniqueness and continuity with respect to the initial condition. –  Baby Dragon Aug 4 '13 at 18:45

(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 solution can be generated based on the Abel function, $\alpha(z)$ of $f(z)=z^2+1$. If we have the Abel function, then the half iterate of z can be generated as $h(z)=\alpha^{-1}(\alpha(z)+0.5)$ Though it is not the only solution (nor in my opinion, the best), the most accessible Abel function solution is based on a Boettcher function for the fixed point at infinity; I wrote a program last year that does just that, from which the results I posted earlier were quickly generated. It is easiest to work with the inverse Boettcher function, from which the inverse Abel function can easily be generated. I'm using use the symbol $\beta$ for the Boettcher function. The problem is that f(x) actually has a super attracting fixed point at infinity, not a fixed point at zero. So, we work with the reciprocal of the $\beta^{-1}$ function. We defined the formal $\beta^{-1}$ function via the following relationship.

$\beta^{-1}(z^2)=\frac{1}{f( \, 1 \, / \, {\beta^{-1}(z) \, })}$

First generate the formal power series for the reciprocal function, 1/f(1/z), which allows one to generate the formal $\beta^{-1}$ series.

$fi(z)=\frac{1}{f(\frac{1}{z})} = z^2 - z^4 + z^6 - z^8 + z^{10} - z^{12} ...$

${\beta^{-1}(z^2)}=fi({\beta^{-1}(z)})$

Now, all you need is the formal power series for the $\beta^{-1}(z)$, along with the equation for the inverse Abel function, in terms of the Boettcher function, and the equation to generate the half iterate in terms of the Abel function, $h(z)=\alpha^{-1}(\alpha(z)+0.5)$.

$\alpha^{-1}(z)=\frac{1}{\beta^{-1}(\exp(-2^{z}))}$

$\beta^{-1}(z)=$

z +
z^ 3*  1/2 +
z^ 5*  5/8 +
z^ 7*  11/16 +
z^ 9*  131/128 +
z^11*  335/256 +
z^13*  1921/1024 +
z^15*  5203/2048 +
z^17*  122531/32768 +
z^19*  342491/65536 +
z^21*  1992139/262144 +
z^23*  5666653/524288 +
z^25*  66211495/4194304 +
z^27*  190532251/8388608 +
z^29*  1112640185/33554432 +
z^31*  3225372323/67108864 +
z^33*  151170463523/2147483648 +
z^35*  440562661907/4294967296 +
z^37*  2583809849479/17179869184 +
z^39*  7558966177753/34359738368 +
z^41*  88836407031661/274877906944...


So, this yields an approximate solution for the superfunction or $\alpha^{-1}(z)$ of $\exp(2^z)$, which is the superfunction for x^2. This approximation is modified by the Boettcher function to become exactly, $\frac{1}{\beta^{-1}(\exp(-2^z)}$. Notice that as z increases, $\exp(-2^z)$ rapidly goes to zero, as long as $|\Im(z)|<\frac{\pi}{2\log(2)}$, and the approximation for the superfunction becomes more and more accurate. This is the Taylor series centered so that $\alpha^{-1}(0)=2$. $\alpha^{-1}(z)=$

        2.00000000000000000000000
+x^ 1*  1.47042970479728200070736
+x^ 2*  0.762480577752927164660093
+x^ 3*  0.424267970164226197579471
+x^ 4*  0.195424007045383357908720
+x^ 5*  0.0885363745236815506982063
+x^ 6*  0.0359598551892287716903761
+x^ 7*  0.0144792452984198575961554
+x^ 8*  0.00535551113121023140421654
+x^ 9*  0.00201219850895305456107215
+x^10*  0.000694227259952985754369526
+x^11*  0.000244367434796641079018478
+x^12*  0.0000769214480826208320220663
+x^13*  0.0000269925934667063689310974
+x^14*  0.00000813609797954979262652707
+x^15*  0.00000283560192079757251765790
+x^16*  0.000000705532363923839906429084
+x^17*  0.000000277796704124709172266365
+x^18*  0.0000000569382653822375560531824
+x^19*  0.0000000291321329124127631158831
+x^20*  0.00000000199960494407016834679507
+x^21*  0.00000000353966190200798175752179
+x^22* -1.84359576880995872838519 E-10
+x^23*  0.000000000489582426965793452585949
+x^24* -1.69340677715894785103962 E-10
+x^25*  9.49659586691303353973779 E-11
+x^26* -2.92631386240628006146382 E-11
+x^27*  1.88357410017244782298422 E-11
+x^28* -9.69806059398720144574851 E-12
+x^29*  4.16913890865704504495135 E-12
+x^30* -1.73667913416272696484187 E-12
+x^31*  9.23380420463300741831335 E-13
+x^32* -4.74750042625944938044382 E-13
+x^33*  2.15998350305014866568442 E-13
+x^34* -9.49184477375019128289258 E-14
+x^35*  4.68362987742936825161002 E-14
+x^36* -2.44077226697947882519346 E-14
+x^37*  1.15019105307459064415620 E-14
+x^38* -5.06531638741476544065356 E-15
+x^39*  2.48236503924756119664440 E-15


The Abel function, and its inverse the superfunction=$\alpha^{-1}(z)$, combine to yield a valid solution for the half iterate using numerical methods to get a Taylor series for $h(z)=\alpha^{-1}(\alpha(z)+0.5)$. I prefer the Cauchy integral, to generate each coefficient of the Taylor series for the half iterate. So below this paragraph is the half iterate, generated by putting iterations of $x^2+1$ into correspondence with iterations of the $x^2$ via the Boettcher super attracting fixed point of infinity/zero. My main reason for preferring the Kneser type solution is that the superfunction generated from the Kneser type solution has no singularities in the upper half of the complex plane, where as the Bottcher function solution is not nearly so well behaved, with an infinite number of singularities as $|\Im(z)|$ approaches $\frac{\pi}{2\log(2)}$. But the Kneser solution requires a Riemann mapping so it is not as accessible as this Boettcher function solution. At the real axis, both functions are very close in values to each other. I haven't studied the half iterates of either in much detail; although the nearest singularity defines the radius of convergence, $\sqrt{1-a_0}\approx 0.598252i$, as noted in my comments above. Here is the half iterate, $h(z)$, for $f(z)=z^2+1$. Notice that the radius of convergence is a little bit too small, so that $h(h(z))$ doesn't converge to $z^2+1$.

         0.642094504390828381495363 +
x^ 2 *  1.00690224867415593906994 +
x^ 4 * -0.514130215253435435599237 +
x^ 6 *  0.733302763753911249332061 +
x^ 8 * -1.32058357980755641903265 +
x^10 *  2.63883845336820960564369 +
x^12 * -5.60443025341316030005301 +
x^14 * 12.4064479200198191890023 +
x^16 *-28.3152137792182421744708 +
x^18 * 66.1663983446023842196175 +
x^20 *-157.550867142011717456622

-

(This is not an answer but a table of data after I've tried Anixx's first/accepted solution)
I get for the first few approximations using Pari/GP

\\ this is the function to be iterated to height "h"
f(x,h) = for(k=1,h,x=x^2+1);x

x = 0.5;
n=4;                                              \\ assume some n
su1=sum(k=0,n,(-1)^k * f(x,k)/k!/(n-k)!/(1/2-k)); \\ compute numerator
su2=sum(k=0,n,(-1)^k         /k!/(n-k)!/(1/2-k)); \\ compute denominator
print([n,su1,su2,  su1/su2]);         \\ the last entry (ratio) is the approximation

%315 = [4, -0.157781292143, 0.304761904762,                   -0.517719864845]
%319 = [5, 5.81430361558, 0.0677248677249,                    85.8518268237]
%323 = [6, -2903.09654248, 0.0123136123136,              -235763.191868]
%327 = [7, 4051027927.89, 0.00189440189440,        2.13842054311 E12]
%331 = [8, -5.82421095518 E22, 0.000252586919254, -2.30582445535 E26]


I don't see how this could be made convergent to some value..

-
around z=0, the functional square root is also an even function, with a leading constant coefficient, and x^2 as the first coefficient. I was able to generate a taylor series, but it has a limited radius of convergence. –  Sheldon L Aug 10 '13 at 1:46
@Sheldon: I would really like to see it. And how you did arrive at it. –  Gottfried Helms Aug 10 '13 at 5:12
Hmm, using Newton iteration to arrive at a squareroot of the carleman-matrix for $f(x)=1+x^2 \qquad (= g(g(x)))$ I could do at best with $g(x) \approx 0.642124541629 + 1.00681024077 x^2 - 0.515947721599 x^4 + 0.750971743189 x^6 + O(x^7)$ (I have a handful more coefficients, but only to at most 14 digits, only at even exponents , alternating in signs and also slowly diverging). Does it approximately match your result? –  Gottfried Helms Aug 10 '13 at 6:09
Hmm, instead of using 64x64 size for the Carleman matrices and 800 digits precision I took 128x128 and 1200 digits precision. This gives $g(x) \approx 0.64209540 + 1.0068539 x^2 - 0.51369998 x^4 + 0.73159216 x^6 + O(x^7)$. It gives $g(0)$ being the constant and $g(g(0)) \approx 1$ to 10 digits precision. However - the change caused by bigger matrix-size occurs even in the fifth digit of the constant... So this left still much room for improvement... –  Gottfried Helms Aug 10 '13 at 8:00
I'm getting a different solution then you. This solution is based on the Kneser method for the superfunction of x^2+1 with two complex fixed points, $0.5+/-\sqrt{-0.75}$. $ht(z)=\alpha^{-1}(\alpha(z)+0.5)$ {ht= 0.642094752506609371698073 +x^ 2* 1.00690493722501474594184 +x^ 4* -0.514140931459684938879986 +x^ 6* 0.733283234297571957782760 +x^ 8* -1.32052177688948603556692 +x^10* 2.63888706127476651531197 +x^12* -5.60464285450601626636908 +x^14* 12.4064962313713263598877 +x^16* -28.3149176972867198028224 +x^18* 66.1664084652015917411194 +x^20* -157.552176646482538089039 } –  Sheldon L Aug 10 '13 at 20:41

This is an additional example for my last comment @Sheldon made as an "answer" because of number of bytes required. The range of convergence of the (first) Kneser-type "ht(x)"-function in Sheldon's comment can be extended by Eulersummation; this can be implemented by simply introducing suitable cofactors for the power series. Here is Pari-code:

{fshel(x,eo=0)=local(evec);
evec=ESumVec(eo,11);   \\ generates the vector of coefficients for E.summation
0.642094752506609371698073      *evec[1]
+x^2* 1.00690493722501474594184 *evec[2]
+x^4* -0.514140931459684938879986  *evec[3]
...
+x^18* 66.1664084652015917411194  *evec[10]
+x^20* -157.552176646482538089039  *evec[11] }


and here the comparision based on the 11 terms only:

fshel(fshel(0,0.0),0.0)     \\ Euler-order 0 = no Euler-summation
%2156 = 0.98887772          \\ should be 1
fshel(fshel(0.0,0.0),0.45)   \\ first call Euler-order 0, second call 0.45
%2157 = 0.99999781          \\ much better

fshel(fshel(0.1,0),0.0)      \\ Euler-order 0 = no Euler-summation
%2152 = 0.99460603          \\ should be   1.01
fshel(fshel(0.1,0.1),0.45)   \\ first call Euler-order 0.1, second call 0.45
%2153 = 1.0099971


The best order of the Euler-summation depends on the value of the x-argument; for instance to sum the alternating series $1-2+4-8+...-...$ we need order 2 and to sum $1-3+9-27+...-...$ we need order 3. The procedure for the computation of the vector is relatively simple; I can put it here if wanted.

-
Time is limited this morning. The solution I gave has a radius of converegence limited by a singularity at $z=\sqrt{a_0-1} \approx 0.598252i$ This is due to a simple square root branch, since $z^2+1=a_0$, and $h(z)=0$. $f(0)=z^2+1$, which seems to cause a simple square root branch singularity, which seen in the half iterate of z, since $h(z)=\sqrt{h(z^2+1)-1}$. So the radius of convergence really is limited to $\approx 0.598252$, and there's nothing you can do to improve the convergence since the function becomes multiple valued. But $\alpha^{-1}(\alpha(z)+0.5)$ solution is well defined. –  Sheldon L Aug 11 '13 at 12:12
I'll post the Bottcher function formal power series solution for the fixed point at infinity later. –  Sheldon L Aug 11 '13 at 12:18