For $\alpha, \beta, \gamma \in (0,1)$ satisfying $\alpha+\beta+\gamma = 1$ and
$\mu \in \mathbb{C} \setminus [1,\infty)$, define
$$
F_{\alpha\beta}(\mu) = \int_0^1\frac{dx}{x^\alpha(1-x)^\beta(1-\mu x)^\gamma}
\quad\text{ and }\quad
\Delta = \frac{\Gamma(1-\alpha)\Gamma(1-\beta)}{\Gamma(1+\gamma)}
$$
When $|\mu| < 1$, we can rewrite the integral $F_{\alpha\beta}(\mu)$ as
$$\begin{align}
F_{\alpha\beta}(\mu)
= & \int_0^1 \frac{1}{x^\alpha(1-x)^{\beta}}\left(\sum_{n=0}^{\infty}\frac{(\gamma)_n}{n!}\mu^n x^n\right) dx
= \sum_{n=0}^{\infty}\frac{(\gamma)_n}{n!}\frac{\Gamma(n+1-\alpha)\Gamma(1-\beta)}{\Gamma(n+1+\gamma)}\mu^n\\
= & \Delta\sum_{n=0}^{\infty}\frac{(\gamma)_n (1-\alpha)_n}{n!(\gamma+1)_n}\mu^n
= \Delta\gamma \sum_{n=0}^{\infty}\frac{(1-\alpha)_n}{n!(\gamma+n)}\mu^n
\end{align}$$
This implies
$$
\mu^{-\gamma} \left(\mu\frac{\partial}{\partial \mu}\right) \mu^{\gamma} F_{\alpha\beta}(\mu) =
\Delta\gamma \sum_{n=0}^{\infty}\frac{(1-\alpha)_n}{n!}\mu^n
= \Delta\gamma\frac{1}{(1-\mu)^{1-\alpha}}
$$
and hence
$$F_{\alpha\beta}(\mu)
= \Delta\gamma \mu^{-\gamma} \int_0^\mu \frac{\nu^{\gamma-1}d\nu}{(1-\nu)^{1-\alpha}}
= \Delta\gamma \int_0^1 \frac{t^{\gamma-1} dt}{(1-\mu t)^{1-\alpha}}
= \Delta \int_0^1 \frac{dt}{(1 - \mu t^{1/\gamma})^{1-\alpha}}$$
Notice if we substitute $x$ by $y = 1-x$, we have
$$F_{\alpha\beta}(\mu) = \int_0^1 \frac{dy}{y^\beta(1-y)^\alpha(1-\mu - \mu y)^{\gamma}}
= \frac{1}{(1-\mu)^\gamma} F_{\beta\alpha}(-\frac{\mu}{1-\mu})$$
Combine these two representations of $F_{\alpha\beta}(\mu)$ and let $\omega = \left(\frac{\mu}{1-\mu}\right)^{\gamma}$, we obtain
$$F_{\alpha\beta}(\mu) = \frac{\Delta}{(1-\mu)^{\gamma}}\int_0^1 \frac{dt}{( 1 + \omega^{1/\gamma} t^{1/\gamma})^{1-\beta}} = \frac{\Delta}{\mu^\gamma}\int_0^\omega \frac{dt}{(1 + t^{1/\gamma})^{1-\beta}}$$
Let $(\alpha,\beta,\gamma) = (\frac14,\frac12,\frac14)$ and $\mu = \frac{\sqrt{3}}{2}$, the identity we want to check becomes
$$\frac{\Gamma(\frac34)\Gamma(\frac12)}{\Gamma(\frac54) (\sqrt{3})^{1/4}}\int_0^\omega \frac{dt}{\sqrt{1+t^4}} \stackrel{?}{=} \frac{2\sqrt{2}}{3\sqrt[8]{3}} \pi\tag{*1}$$
Let $K(m)$ be the complete elliptic integral of the first kind associated with modulus $m$. i.e.
$$K(m) = \int_0^1 \frac{dx}{\sqrt{(1-x^2)(1-mx^2)}}$$
It is known that $\displaystyle K(\frac12) = \frac{8\pi^{3/2}}{\Gamma(-\frac14)^2}$. In term of $K(\frac12)$, it is easy to check $(*1)$ is equivalent to
$$\int_0^\omega \frac{dt}{\sqrt{1+t^4}} \stackrel{?}{=} \frac23 K(\frac12)\tag{*2}$$
To see whether this is the case, let $\varphi(u)$ be the inverse function of above integral.
More precisely, define $\varphi(u)$ by following relation:
$$u = \int_0^{\varphi(u)} \frac{dt}{\sqrt{1+t^4}}$$
Let $\psi(u)$ be $\frac{1}{\sqrt{2}}(\varphi(u) + \varphi(u)^{-1})$. It is easy to check/verify
$$
\varphi'(u)^2 = 1 + \varphi(u)^4
\implies
\psi'(u)^2 = 4 (1 - \psi(u)^2)(1 - \frac12 \psi(u)^2)
$$
Compare the ODE of $\psi(u)$ with that of a Jacobi elliptic functions with modulus $m = \frac12$, we find
$$\psi(u) = \text{sn}(2u + \text{constant} | \frac12 )\tag{*3}$$
Since we are going to deal with elliptic functions/integrals with $m = \frac12$ only,
we will simplify our notations and drop all reference to modulus, i.e
$\text{sn}(u)$ now means $\text{sn}(u|m=\frac12)$ and $K$ means $K(m = \frac12)$.
Over the complex plane, it is known that $\text{sn}(u)$ is doubly periodic with
fundamental period $4 K$ and $2i K$. It has two poles at $i K$ and $(2 + i)K$ in the fundamental domain.
When $u = 0$, we want $\varphi(u) = 0$ and hence $\psi(u) = \infty$. So the constant
in $(*3)$ has to be one of the pole. For small and positive $u$, we want $\varphi(u)$ and hence $\psi(u)$ to be positive. This fixes the constant to $i K$. i.e.
$$\psi(u) = \text{sn}(2u + iK )$$
and the condition $(*2)$ becomes whether following equality is true or not.
$$\frac{1}{\sqrt{2}} (\omega + \omega^{-1}) \stackrel{?}{=} \text{sn}( \frac43 K + i K)\tag{*4}$$
Notice $ 3( \frac43 K + i K) = 4 K + 3 i K $ is a pole of $\text{sn}(u)$. if one repeat
apply the addition formula for $\text{sn}(u+v)$
$$\text{sn}(u+v) = \frac{\text{sn}(u)\text{cn}(v)\text{dn}(v)+\text{sn}(v)\text{cn}(u)\text{dn}(u)}{1-m\,\text{sn}(u)^2 \text{sn}(v)^2}$$
One find in order for $\text{sn}(3u)$ to blow up, $\text{sn}(u)$ will be a root of
following polynomial equation:
$$3 m^2 s^8-4 m^2 s^6-4 m s^6+6 m s^4-1 = 0$$
Substitute $m = \frac12$ and $s = \frac{1}{\sqrt{2}}(t+\frac{1}{t})$ into this, the equation $\omega$ need to satisfy is given by:
$$(t^8 - 6 t^4 - 3)(3 t^8 + 6 t^4 - 1 ) = 0$$
One can check that $\omega = \sqrt[4]{\frac{\sqrt{3}}{2-\sqrt{3}}}$ is indeed a root of this polynomial. As a result, the original equality is valid.
Best Answer
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$.