This question is more complicated than it seems

If $f^{-1}(x)=f'(x)$ then $f(f^{-1}(x))=f(f'(x))$ or $x=f(f'(x))$

First derivative is giving

$\displaystyle 1=f'(f'(x))f''(x)$

Now replace $y=f'(x)$

$\displaystyle y(y(x))y'(x)=1$

Instead of the above we will solve a more general:

$\displaystyle g(y(x))y'(x)=1$

This one has a solution expressed as

$\displaystyle \int_{1}^{g(x)}f(r) \mathrm{d} r = x+c$

Now since we want $f(x)=g(x)$ we have

$\displaystyle \int_{1}^{f(x)}f(r) \mathrm{d} r = x+c$

$\int_{1}^{f(x)}f(r) \mathrm{d} r$ is actually a curious operator which we will represent as $\mathfrak{R}(f(x))_{1}$

Now this operator has its inverse and this is what we are looking for. For a function $h(x)$ we want to know $s(x)$ so that

$\displaystyle \mathfrak{R}(s(x))=h(x)$

We will call $\mathfrak{R}^{-1}()$ the integral root.

As a new operator we need to examine it a little bit. We can create a table for some known functions first

$\begin{matrix} h(x) & \mathfrak{R}^{-1}(h(x))_{1} \\ & \\ \frac{1}{2}(x^2-1) & x\\ \frac{1}{3}(x^6-1) & x^2\\ 1 - \ln(x)(1 - \ln(\ln(x))) & \ln(x)\\ e^{e^x}-e & e^x\\ -\ln(x) & \frac{1}{x}\\ \cos(1)-\cos(\sin(x)) & \sin(x)\\ \frac{2}{3}(x^{\frac{3}{4}}-1) & \sqrt{x} \end{matrix}$

We do not ignore constant term since integral root is very sensitive operator and we define in general integral root with base $b$ as:

$\mathfrak{R}(h(x))_{b}=\int_{b}^{f(x)}f(r) \mathrm{d} r$

We will return to the importance of base later but for now we say that the function does not have to be defined everywhere and base can help about it. So, all we need to find the integral root of $h(x)=x$.

The operator resembles a normal derivative/integration except that it is extending it all, it is giving sort of a faster result.

It is a deep question if this operator has unique values. (We will hint something into that direction but overall this is a very good subject for some semester work.)

Notice that the operator is very sensitive to the constant value and this one cannot be ignored.

Since powers of $x$ appears from the powers of $x$, we could try the solution in the most general form as $g(x)=ax^b+f(x)$ however this will convince us very quickly that for any solution defined everywhere $f(x)=0$, since such solutions are very restrictive and it is not possible to drive constants the way we would like to, even adding $ax^b+c$ would create a solution that is useless. The way of dealing with constant term is by changing the base in general.

So, if we restrict ourselves to global world we have

$\displaystyle \mathfrak{R}(ax^b)_{1}=\frac{a^{b+2}x^{b(b+1)}}{b+1}-\frac{a}{b+1}$

Since we want

$\displaystyle \mathfrak{R}(ax^b)_{1}= x+c$

that requires $b(b+1)=1$ making the result

$ba^{b+2}x-ab$

This requires $ba^{b+2}=1$ as well which makes $a=(b+1)^{\frac{1}{b+2}}$

(Notice that these are not two different integral roots of the same function.)

Now our solution is $y(x)=f'(x)$ so we have it back as

$\displaystyle f(x)=\frac{a}{b+1}x^{b+1}=(b+1)^{-\frac{1}{b+1}}x^{b+1}$

where as we have mentioned $b(b+1)=1$

The reason we opt for this operator is that the function in question has derivative thus it is assumed that it is a nicely behaving function.

Let us prove that these two are the solutions

Derivative

$\displaystyle f'(x)=(b+1)^{-\frac{1}{b+1}}(b+1)x^{b}$

$\displaystyle f'(x)=(b+1)^{-\frac{1}{b+1}+1}x^{b}$

$\displaystyle f'(x)=(b+1)^{\frac{b+1-1}{b+1}}x^{b}$

$\displaystyle f'(x)=(b+1)^{\frac{b}{b+1}}x^{b}$

Inverse

$\displaystyle x=(b+1)^{-\frac{1}{b+1}}h(x)^{b+1}$

$\displaystyle (b+1)^{\frac{1}{b+1}}x=h(x)^{b+1}$

$\displaystyle (b+1)^{\frac{1}{(b+1)(b+1)}}x^{\frac{1}{b+1}}=h(x)$

$\displaystyle f^{-1}=h(x)=(b+1)^{\frac{b}{b+1}}x^{b}$

Finally let us use golden ratio to find what the solutions we are talking about

$\displaystyle b_{1}=-\phi$,$ b_{2}=\phi-1$

making the first solution actually complex and the second real:

$\displaystyle (-\phi+1)^{-\frac{1}{-\phi+1}}x^{-\phi+1}=(-\phi)^{-\phi}x^{-\frac{1}{\phi}}$

$\displaystyle \phi^{-\frac{1}{\phi}}x^{\phi}$

The two solutions we have are driven by differing a constant only, otherwise the integral root has the unique solution just like normal root has. It is possible that there are solutions that are not analytical or defined everywhere, otherwise the two given solutions are the only nice global solutions.

Apart from these global solutions that do not change if we change the integral root base, there is an option of having a local solution. For this we extend the function analytically at some point $x_{0}$ and ask that the function behaves as $x$ regarding the integral root at that point.

Then near $x_{0}$ that function will behave as required.

$\displaystyle f(x)=\sum\limits_{n=0}^{+\infty}c_{n}(x-x_{0})^n $

Basically we are asking

$\displaystyle \mathfrak{R}(\sum\limits_{n=0}^{+\infty}c_{n}(x-x_{0})^n)_{x_{0}} = x+c$

Notice that the change of base is very important, albeit purely technical for this particular task, as we have assumed that we know this function only around $x_{0}$

Then we need to have:

$\displaystyle \sum\limits_{n=0}^{+\infty}\frac{c_{n}}{n+1}(\sum\limits_{m=0}^{+\infty}c_{m}(x-x_{0})^m-x_{0})^{n+1}=x+f(x_{0})+c$

Obviously we can match all coefficients since the final coefficient by $x$ should be $1$ and all higher powers $0$. Solving all coefficient for $c_{k}$ will give local solutions for each $x_{0}$. We already have one such non trivial solution, but there could be more than that global one.

However, notice that the integral root base does not affect the answer, since the base deals with constant term only, and constant term is arbitrary. Instead of $1$ for our function we could take any constant.

This is to say that our local solution again has the same resolution which we have found when we started from $1$.

Unless there is some curious and undiscovered way the integral root can have more than two given solutions, one real and another complex, each local solution is equal to either of the two found global solutions.