$f^3,f^2$ are the cube and quadratic of f respectively and both infinite differentiable on $R$,how to show so is $f$

The following papers prove this:

MR0682456 Reviewed Joris, Henri Une C∞-application non-immersive qui possède la propriété universelle des immersions. (French) [A nonimmersive C∞ mapping having the universal property of immersions] Arch. Math. (Basel) 39 (1982), no. 3, 269–277.

MR0833407 Reviewed Duncan, John; Krantz, Steven G.; Parks, Harold R. Nonlinear conditions for differentiability of functions. J. Analyse Math. 45 (1985), 46–68. (Reviewer: Wiesław Pleśniak) 26E10 (58C25)

MR2179865 Reviewed Myers, Robert An elementary proof of Joris's theorem. Amer. Math. Monthly 112 (2005), no. 9, 829–831. (Reviewer: Clifford E. Weil) 26A24

Posted here because Criterion for smooth functions got closed when I was typing all this. Look at the comment chain there if the opening phrases make no sense to you :-)

Since I couldn't recall the name of the theorem, let me at least try to recall the $\bar\partial$ proof of it :-). Also, if I remember things right (no warranty on that), the proof in Monthly is one of the flawed ones, so at the very least be extremely careful when reading it. The proof referred to by Iosif is correct, brilliant, and short too, works in much more general setting, but still has one small drawback: looking at it, it is not so easy to understand what happens when only some finite smoothness is assumed, i.e., what is the implication of $f^2,f^3\in C^N$ with some finite $N$. Note that this implication cannot be that $f\in C^N$ because $f(x)=0$ for $x<0$, $f(x)=x^{N/2+\delta}$ for $x\ge 0$ satisfies the requirements but is only $C^{N/2}$ smooth. However, it is, indeed, possible to show that $f\in C^{rN}$ for some fixed $r>0$ (though I'll not try to optimize $r$ in the argument below)

We need three pieces of standard machinery. I'll present the cheapest versions that will work for us though almost everything can be refined quite a bit.

The first one is that a function $f(x)$ on an interval $I\subset\mathbb R$ is in $C^N$ essentially if and only if there exists a uniformly bounded family of functions $f_\varepsilon(z)$ ($\varepsilon\in(0,1)$) defined and analytic in the strips $D_\varepsilon=I\times(-\varepsilon,\varepsilon)$ such that $|f-f_\varepsilon|\le\varepsilon^N$ on $I$.

To get one direction, just define $f(x+iy)=\sum_{k=0}^{N-1}\frac{f^{(k)}(x)}{k!}(iy)^k$ formally extending the Taylor series of $f$ to the complex plane from the nearest point. Then $(\bar\partial f)(x+iy)=\frac 12(\partial_x+i\partial_y)f(x+iy)=\frac 12\frac{f^{(N)}(x)}{(N-1)!}(iy)^{N-1}$, which is bounded by $C\varepsilon^{N-1}$ in $D_\varepsilon$. To get an analytic approximation, it is enough just to subtract the solution of the $\bar\partial$ problem that is obtained by the convolution of $\chi_{D_\varepsilon}\bar\partial f$ with $\frac 1{\pi z}$ as usual.

To get the other direction, let us notice that $|f_{2\varepsilon}-f_{\varepsilon}|\le C\varepsilon^N$ on $I$, so, due to the uniform boundedness and something like the Hadamard 3 lines lemma (or, if you prefer, crude estimates for the harmonic measure), we get $|f_{2\varepsilon}-f_{\varepsilon}|\le C\varepsilon^{N(1-q)}$ in $D_{q\varepsilon}$ for any $q\in(0,\frac 12)$. Thus, decomposing
$$ f=f_1-\sum_{\varepsilon\in E} (f_{2\varepsilon}-f_{\varepsilon}) $$ where $E=\{2^{-k},k=1,2,3,\dots\}$ and applying the Cauchy estimate for the $m$-th derivative, we conclude that the formally differentiated $m$ times series on the right is dominated by $$ Cm!\sum_{\varepsilon\in E}(q\varepsilon)^{-m} \varepsilon^{N(1-q)}\,, $$ which is a convergent series when $m<N(1-q)$, so if we choose $q<N^{-1}$, we shall get $f\in C^{N-1}$.

Note that we were a bit sloppy here: we went from $N$ to $N-1$ and the "three line lemma" works the way we want only if we are separated from the endpoints of $I$, but since the smoothness is a local property and we are ready to lose even some portion of $N$ in our count of derivatives, such sloppiness can be easily tolerated.

The second piece of machinery is the computation of the Laplacian $$ \Delta \log(|h|^2+\rho^2)=4\frac{\rho^2|h'|^2}{(|h|^2+\rho^2)^2} $$ and the Green formula $$ \int_{\mathbb D}\Delta u(z)\log(1/|z|)\,dA(z)=\int_{\mathbb T} u(z)dm(z)-2\pi u(0) $$ where $\mathbb D=\{z:|z|<1\}$, $\mathbb T=\{z:|z|=1\}$, $dA$ is the area measure and $dm$ is the length measure.

It follows immediately that if $h$ is a bounded analytic function in the unit disk, then $$ \int_{\frac 12\mathbb D}|h'|^2\chi_{\{|h|\le\rho\}}\le C\log\frac{\|h\|_\infty^2+\rho^2}{\rho^2}\rho^2\,. $$ Rescaling to $\varepsilon\mathbb D$ and covering $D_{\varepsilon/3}$ by $\approx \varepsilon^{-1}$ disks of radius $\varepsilon/2$, we conclude that for any analytic function $h:D_\varepsilon\to\mathbb C$ that is bounded by a constant, we have $$ \int_{D_{\varepsilon/3}}|h'|^2\chi_{\{|h|\le\rho\}}\le C\varepsilon^{-1}\log\frac{1+\rho^2}{\rho^2}\rho^2\,. $$

The last piece of machinery is the bound for the uniform norm of the convolution $w*\frac 1{\pi z}$ in a bounded domain. Split $\frac 1{\pi z}=K_r+K^r$ where $K_r$ is the part in the disk $r\mathbb D$ and $K^r$ is the part outside the disk (we assume that the function is truncated beyond the diameter of the domain). Then $\|K_r\|_{L^1}\le Cr$ and $\|K^r\|_{L^2}\le C\sqrt{\log(1/r)}$ for small $r$. Hence, we get $$ \left\|w*\frac 1{\pi z}\right\|\le Cr\|w\|_{L^\infty}+C\sqrt{\log(1/r)}\|w\|_{L^2}. $$ for every $r\in(0,\frac 12)$, say.

Now everything is ready for the "punchline". Assume that $f^2$ and $f^3$ are $C^{N}$ on $I$ with derivatives up to $N$ bounded by $1$, say. Fix $\varepsilon>0$. Construct their bounded analytic approximations in $D_{\varepsilon}$ with error $\varepsilon^N$ and call them $h$ and $g$ respectively. Then $|g^2-h^3|\le C\varepsilon^N$ on the real line, so if $q<1/N$, this estimate (with $N(1-q)>N-1$ instead of $N$) spreads to $D_{q\varepsilon}$. Thus, replacing $\varepsilon$ by $q\varepsilon$, we can safely assume that $|g^2-h^3|\le\delta=C\varepsilon^{N-1}$ in $D_{\varepsilon}$.

How would a normal person recover $f$ grom $g\approx f^3$ and $h\approx f^2$? One possibility is just to say $f\approx g/h$. This would be perfect if we didn't have the small (or even $0$) denominator problem, so let us take $F=\frac{g\bar h}{\max(|h|^2,\rho^2)}$ instead with some $\rho\gg \delta$ to be chosen later and put $f_\varepsilon=F-(\bar\partial F)*\frac 1{\pi z}$. We just need to do two things now: to estimate $|F-f|$ on $I$ and to bound the convolution in the uniform norm.

Estimate of $|F-f|$

Case 1: $|h|<\rho$. In this case $|f|\le \sqrt(|h|+\delta)\le \sqrt{\rho+\delta}$ and $\frac{|g||h|}{\rho^2}\le \frac{\sqrt{\rho^3+\delta}}{\rho}\le C\sqrt{\rho}$, provided that $\rho^3\ge \delta$. Thus $|F-f|\le C\sqrt{\rho}$.

Case 2: $|h|\ge\rho$. In this case $h=f^2+O(\delta)$, $g=f^3+O(\delta)$, so $|f|\ge \sqrt{\rho-O(\delta)}\approx \sqrt\rho$, so $|f|^3\ge \rho^{3/2}\gg \delta$ and $$ F=\frac gh=\frac{f^3+O(\delta)}{f^2+O(\delta)}=f+O\left( |f|\left(\frac{\delta}{|f|^3}+\frac{\delta}{|f|^2}\right) \right)=f+O(\delta/|f|^2)=f+O(\delta/\rho)\,, $$ whence $|F-f|\le \frac\delta\rho$ in this case.

Thus, we get $|F-f|\le C(\sqrt{\rho}+\frac{\delta}{\rho})$ whenever $\rho^3\ge\delta$.

Estimate of the convolution

Notice that $w=\bar\partial F=\rho^{-2}g\bar{h'}\chi_{\{|h|\le\rho\}}$. Also, we have already seen that $|h|\le\rho$ implies $|g|\le C\rho^{3/2}$ if $\rho^3\ge\delta$, so we get $$ |w|\le C\rho^{-1/2}|h'|\chi_{\{|h|\le\rho\}}\,. $$ Thus, $\|w\|_{L^2}^2\le C\varepsilon^{-1}\log\frac{1+\rho^2}{\rho^2}\rho$ in $D_{\varepsilon/3}$. Also, since $h$ is bounded in $D_\varepsilon$, we have $|h'|\le C/\varepsilon$ and $\|w\|_{L^\infty}\le C\rho^{-1/2}\varepsilon^{-1}$ in $D_{\varepsilon/3}$.

Hence, by the last standard estimate, we get the bound $$ Cr\rho^{-1/2}\varepsilon^{-1}+C\sqrt{\varepsilon^{-1}\rho\log\frac{1+\rho^2}{\rho^2}}\sqrt{\log(1/r)} $$ for the $L^\infty$ norm of $w*\frac 1{\pi z}$. If we choose $r=\rho$, this becomes $C\varepsilon^{-1}\sqrt\rho\log(1/\rho)$.

Now it is easy to see that with $\rho^3=\delta=C\varepsilon^{N-1}$, all our bounds are at most $C\varepsilon^{\frac N6-2}\log(1/\varepsilon)$, so we have shown that we have at least $\frac N6-3$ continuous derivatives (and there is still some room for improvement).