Integrals of $\sqrt{x+\sqrt{\phantom|\dots+\sqrt{x+1}}}$ in elementary functions Let $f_n(x)$ be recursively defined as
$$f_0(x)=1,\ \ \ f_{n+1}(x)=\sqrt{x+f_n(x)},\tag1$$
i.e. $f_n(x)$ contains $n$ radicals and $n$ occurences of $x$:
$$f_1(x)=\sqrt{x+1},\ \ \ f_2(x)=\sqrt{x+\sqrt{x+1}},\ \ \ f_3(x)=\sqrt{x+\sqrt{x+\sqrt{x+1}}},\ \dots\tag2$$
The functions $f_0(x)$, $f_1(x)$ and $f_2(x)$ are integrable in elementary functions, e.g.:
$$\int\sqrt{x+\sqrt{x+1}}\,dx=\left(\frac{2\,x}3+\frac{\sqrt{x+1}}6-\frac14\right)\sqrt{x+\sqrt{x+1}}+\frac58\ln\left(2\,\sqrt{x+1}+2\,\sqrt{x+\sqrt{x+1}}+1\right).\tag3$$


Question: Is there an integer $n>2$ such that $f_n(x)$ is integrable in elementary functions?


Update: The question is reposted at MathOverflow as was suggested by moderator.
 A: I posted an answer to this question as an answer to the related question on MathOverflow. See my (second) answer at https://mathoverflow.net/questions/171733.  However, as Krokop pointed out, it would be better to have a copy here on Math.SE to the original question.  Thanks to Krokop for copying my answer over to Math.SE.

I am going to show that there is no elementary antiderivative of $f_n$ when  $n>2$.

Assume $n>2$ (NB: This is important, because the argument below will not work for $n\le2$; the reader may enjoy finding where it breaks down), and let $K_n = {\mathbb C}\bigl(x,f_n(x)\bigr)$ be the elementary differential field generated by $x$ and $f_n(x)$.  Then $K_n$ is the field of meromorphic functions on the normalization $\hat C_n$ of the algebraic curve $C_n$ defined by the minimal degree $y$-monic polynomial $P_n(x,y)$ that satisfies $P_n\bigl(x,f_n(x)\bigr) \equiv 0$.  This minimal degree is $2^n$; for example, $P_2(x,y) = (y^2-x)^2-x-1$ and $P_3(x,y) = \bigl((y^2-x)^2-x\bigr)^2-x-1$, etc.

Since $P_{n+1}(x,y) = (P_n(x,y)+1)^2-x-1$ for $n\ge 1$ with $P_1(x,y)=y^2-x-1$, one sees, by applying the Eisenstein Criterion to $P_n(x,y)$ regarded as an element of $D[y]$ with $D$ being the integral domain ${\mathbb C}[x]$, that $P_n(x,y)$ is irreducible for all $n\ge 1$.  Hence, $\hat C_n$ is connected. 

It will be important in what follows to observe that $K_n$ has an involution $\iota$ that fixes $x$ and sends $f_n(x)$ to $-f_n(x)$; this is because $P_n(x,y)$ is an even polynomial in $y$.  The fixed field of $\iota$ is ${\mathbb C}\bigl(x,\,f_n(x)^2\bigr)$, and the $(-1)$-eigenspace of $\iota$ is ${\mathbb C}\bigl(x,\,f_n(x)^2\bigr)f_n(x) = K_{n-1}{\cdot}f_n(x)$.

Now, the curve $C_n\subset \mathbb{CP}^2$ has only one point on the line at infinity, namely $[1,0,0]$, but the normalization $\hat C_n$ has $2^{n-1}$ points lying over this point.  They can be parametrized as follows:  First, establish the convention that $\sqrt{u}$ means the unique analytic function on the complex $u$-plane minus its negative axis and $0$ that satisfies $\sqrt1 = 1$ and $\bigl(\sqrt{u}\bigr)^2 = u$.  Let $\epsilon = (\epsilon_1,\ldots,\epsilon_{n-1})$ be any sequence with ${\epsilon_k}^2=1$ and consider the sequence of functions $g^\epsilon_k(t)$ defined by the criteria $g^\epsilon_1(t) = \sqrt{1+t^2}$ and $g^\epsilon_{k+1}(t) = \sqrt{1+\epsilon_{n-k}t g^\epsilon_k(t)}$ for $1\le k < n$.  Choose, as one may, a $\delta_n>0$ sufficiently small so that, when $t$ is complex and satisfies $|t|<\delta_n$, all of the functions $g^\epsilon_k$ are analytic when $|t|<\delta_n$.  In particular, one finds an expansion
$$
g^\epsilon_n(t) 
= 1+\tfrac12\epsilon_1\,t + \tfrac18(2\epsilon_1\epsilon_2-1)t^2 + O(t^3).
$$

Also, it is easy to verify that the disk in $\mathbb{CP}^2$ defined by
$$
[x,y,1] = [1,\ t g^\epsilon_n(t),\ t^2]\qquad\text{for}\quad |t|<\delta_n
$$
is a nonsingular parametrization of a branch of $C_n$ in a neighborhood of the point $[1,0,0]$.  In the normalization $\hat C_n$, this is then a local parametrization of a neighborhood of a point $p_\epsilon\in \hat C_n$.
Obviously, this describes $2^{n-1}$ distinct points on $\hat C_n$. 

When $x$ and $f_n$ are regarded as meromorphic functions on $\hat C_n$, 
it follows that there is a unique local coordinate chart $t_\epsilon:D_\epsilon\to D(0,\delta_n)\subset \mathbb{C}$ of an open disk $D_\epsilon\subset \hat C_n$ about $p_\epsilon$ such that $t_\epsilon(p_\epsilon)=0$ and on which one
has formulae
$$
x = \frac1{{t_\epsilon}^2}
\quad\text{and}\quad
f_n(x) = \frac{g^\epsilon_n(t_\epsilon)}{t_\epsilon} 
= \frac{1+\tfrac12\epsilon_1\ t_\epsilon 
         +\tfrac18(2\epsilon_1\epsilon_2-1)\ {t_\epsilon}^2}
   {t_\epsilon} 
+ O({t_\epsilon}^2).
$$ 
In particular, it follows that $f_n(x)$, as a meromorphic function on $\hat C_n$,
has polar divisor equal to the sum of the $p_\epsilon$ and hence has degree $2^{n-1}$. Of course, this implies that the zero divisor of $f_n(x)$ on $\hat C_n$ must be of degree $2^{n-1}$ as well. 

Note that the functions $g^\epsilon_k$ satisfy $g^{-\epsilon}_k(-t) = g^{\epsilon}_k(t)$, where $-\epsilon = (-\epsilon_1,\ldots,-\epsilon_{n-1})$.
This implies that $\iota(p_\epsilon) = p_{-\epsilon}$ and that
$t_\epsilon\circ\iota = -t_{-\epsilon}$.

Now, the $2^{n-1}$ zeroes of $f_n(x)$ on $\hat C_n$ are distinct, for they are the zeros of the polynomial $q_n(x) = P_n(x,0) = (q_{n-1}+1)^2-x-1$, and the discriminant of $q_n$, being the resultant of $q_n$ and $q_n'$, is clearly an odd integer, and hence is not zero.  Thus, $C_n$ is a branched double cover of $C_{n-1}$, branched exactly where $f_{n}$ has its zeros. This induces a branched cover $\pi_n:\hat C_n\to \hat C_{n-1}$ that is exactly the quotient of $\hat C_n$ by the involution $\iota$ (whose fixed points are where $f_n$ has its zeros).  Since one then has the Riemann-Hurwitz formula
$$
\chi(\hat C_n) = 2\chi(\hat C_{n-1}) - B_n = 2\chi(\hat C_{n-1}) - 2^{n-1},
$$
and $\chi(\hat C_1) = \chi(\hat C_2) = 2$, induction gives $\chi(\hat C_n) = 
(3{-}n)2^{n-1}$, so the genus of $\hat C_n$ is $(n{-}3) 2^{n-2} + 1$.  (This won't actually be needed below, but it is interesting.)

The only poles of $x$ and $f_n(x)$ on $\hat C_n$ are the points $p_\epsilon$, 
and computation using the above expansions shows that, 
in a neighborhood of $p_\epsilon$, one has an expansion of the form
$$
f_n(x)\,\mathrm{d} x 
- \mathrm{d}\left(f_n(x)\bigl(\tfrac12\ x + \tfrac16\ f_n(x)^2\bigr) \right)
= \left(\frac{ (1-\epsilon_1\epsilon_2) }
         {4{t_\epsilon}^2}
          + O({t_\epsilon}^{-1})\right)\ \mathrm{d} t_\epsilon\ .
$$
Thus, the meromorphic differential $\eta$ on $\hat C_n$ 
defined by the left hand side of this equation has, at worst, double poles 
at the points $p_\epsilon$ and no other poles.

Now, by Liouville's Theorem, $f_n$ has an elementary antiderivative if and only if $f_n(x)\ \mathrm{d} x$ and, hence, the form $\eta$ are expressible as finite linear combinations of exact differentials and log-exact differentials.
Thus, $f_n(x)$ has an elementary antiderivative if and only if $\eta$ is expressible in the form
$$
\eta = \mathrm{d} h + \sum_{i=1}^m c_i\,\frac{\mathrm{d} g_i}{g_i}
$$
for some $h,g_1,\cdots g_m\in K_n$ and some constants $c_1,\ldots,c_m$.  Suppose that these exist.  Since $\eta$ has, at worst, double poles at the $p_\epsilon$ and no other poles, it follows that $h$ must have, at worst, simple poles at the points $p_\epsilon$ and no other poles; in fact, $h$ is uniquely determined up to an additive constant because its expansion at $p_\epsilon$ in terms of $t_\epsilon$ must be of the form
$$
h = \frac{\epsilon_1\epsilon_2-1}{4t_\epsilon} + O(1).
$$
Moreover, because $\eta$ is odd with respect to $\iota$, it follows that $h$ (after adding a suitable constant if necessary) must also be odd with respect to $\iota$.  This implies, in particular, that $h$ vanishes at each of the zeros of $f_n$ (which, by the argument above, are simple zeros).  This implies that $h = r\,f_n$ for some $r\in K_{n-1}$ that has no poles and satisfies $r(p_\epsilon) = (\epsilon_1\epsilon_2-1)/4$ for each $\epsilon$.  However, since $r$ has no poles and $\hat C_n$ is connected, it follows that $r$ is constant.  Thus, it cannot take the two distinct values $0$ and $-1/2$, as the equation $r(p_\epsilon) = (\epsilon_1\epsilon_2-1)/4$ implies.

Thus, the desired $h$ does not exist, and $f_n$ cannot be integrated in elementary terms for any $n>2$.
A: If you pick $f_0 = 1, K_0 = \Bbb C(X)$ and $f_{n+1} = \sqrt{X+f_n}, K_{n+1} = K_n(f_{n+1})$, the first three fields have genus $0$, so integrating $f_{n}$ is as easy as integrating a rational fraction.
For example, if we put $w = \frac {f_1-1}{f_2-1}$ then $K_2 = \Bbb C(w)$ and we have formulas
$$ x = \frac{3w^4-8w^3+10w^2-4w}{(w^2-1)^2},
f_1 = \frac{-2w^2+2w-1}{w^2-1},
f_2 = \frac{w^2-3w+1}{w^2-1}
$$
Then integrating $f_2dx$ is integrating $f_2(dx/dw) dw$, and we obtain a result in terms of $w$ and $\log(w^2-1)$.
However, $K_3 = K_2(\sqrt{x+f_2}) = K_2(v = \sqrt{4w^4-11w^3+10w^2-w-1})$ is obtained by adjoining a square root of an element of degree $4$, so $K_3$ is an elliptic curve (the following fields' genus then grow exponentially). And you want to integrate the $1$-form $f_3 dx$ on this curve.
If I am not mistaken, the poles of the $1$-form are at the $4$ points at $x = \infty$ (or rather $w^2 = 1$). A computation shows that the residues are $+1/4$ around the points (in the coordinates $(v,w)$) $A =(+5,-1), B = (+1,+1)$ and $-1/4$ around the points $C = (-5,-1), D = (-1,+1)$.
Sadly, $[A]+[B]-[C]-[D] \neq 0$. I think this is a point of infinite order on the elliptic curve so there is no rational function we can plug into a $\log$ to cancel the residues of the $1$-form.
So unless I'm mistaken, a primitive of $f_3dx$ can't even be expressed using logarithms and elliptic functions.
A: I think there is highly have no chance to have an integer $n>2$ such that $f_n(x)$ is integrable in elementary functions.
The reasons are mainly because of the processes of eliminating the radicals.
For $f_0(x)$ and $f_1(x)$ , they are trivially integrable in elementary functions.
For $\int f_2(x)~dx$ ,
Let $u=\sqrt{x+1}$ ,
Then $x=u^2-1$
$dx=2u~du$
$\therefore\int\sqrt{x+\sqrt{x+1}}~dx$
$=\int2u\sqrt{u^2+u-1}~du$ , which can express in elementary functions.
Starting from $\int f_3(x)~dx$ ,
Let $u=\sqrt{x+1}$ ,
Then $x=u^2-1$
$dx=2u~du$
$\therefore\int\sqrt{x+\sqrt{x+\sqrt{x+1}}}~dx$
$=\int2u\sqrt{u^2-1+\sqrt{u^2+u-1}}~du$
Introduce the Euler substitution:
Let $v=u+\sqrt{u^2+u-1}$ ,
Then $u=\dfrac{v^2+1}{2v+1}$
$du=\dfrac{2v(2v+1)-(v^2+1)2}{(2v+1)^2}dv=\dfrac{2v^2+2v-2}{(2v+1)^2}dv$
$\therefore\int2u\sqrt{u^2-1+\sqrt{u^2+u-1}}~du$
$=\int2\dfrac{v^2+1}{2v+1}\sqrt{\left(\dfrac{v^2+1}{2v+1}\right)^2-1+v-\dfrac{v^2+1}{2v+1}}\dfrac{2v^2+2v-2}{(2v+1)^2}dv$
$=\int4\dfrac{(v^2+1)(v^2+v-1)}{(2v+1)^3}\sqrt{\dfrac{(v^2+1)^2-(v^2+1)(2v+1)+(v-1)(2v+1)^2}{(2v+1)^2}}~dv$ , which is highly have no chance to express in elementary functions.
