According to Mathematica, the sum is
$$ \frac{3}{\Gamma(\frac13)\Gamma(\frac23)}\left( -1 + {}_3F_2\left(\frac14,\frac12,\frac34; \frac23,\frac43; -1\right) \right). $$
This form is actually quite straightforward if you write out $(4n)!$ as
$$ 4^{4n}n!(1/4)_n (1/2)_n (3/4)_n $$
using rising powers ("Pochhammer symbols") and then use the definition of a hypergeometric function.
The hypergeometric function there can be handled with equation 25 here: http://mathworld.wolfram.com/HypergeometricFunction.html:
$$ {}_3F_2\left(\frac14,\frac12,\frac34; \frac23,\frac43; y\right)=\frac{1}{1-x^k},$$
where $k=3$, $0\leq x\leq (1+k)^{-1/k}$ and
$$ y = \left(\frac{x(1-x^k)}{f_k}\right)^k, \qquad f_k = \frac{k}{(1+k)^{(1+1/k)}}. $$
Now setting $y=-1$, we get the polynomial equation in $x$
$$ \frac{256}{27} x^3 \left(1-x^3\right)^3 = -1,$$
which has two real roots, neither of them in the necessary interval $[0,(1+k)^{-1/k}=4^{-1/3}]$, since one is $-0.43\ldots$ and the other $1.124\ldots$. However, one of those roots, $x_1=-0.436250\ldots$ just happens to give the (numerically at least) right answer, so never mind that.
Also, note that
$$ \Gamma(1/3)\Gamma(2/3) = \frac{2\pi}{\sqrt{3}}. $$
The polynomial equation above is in terms of $x^3$, so we can simplify that too a little,
so the answer is that the sum equals
$$ \frac{3^{3/2}}{2\pi} \left(-1+(1-z_1)^{-1}\right), $$
where $z_1$ is a root of the polynomial equation
$$ 256z(1-z)^3+27=0, \qquad z_1=-0.0830249175076244\ldots $$
(The other real root is $\approx 1.42$.)
How did you find the conjectured closed form?
In this answer, whereas the conjectured formula will not be confirmed so far, a simpler non-recursive, explicit formula for the Fabius function will be offered, which is expressed in terms similar to, but simpler than, those in the conjectured formula.
As noted in the Wikipedia article on the Fabius function, on the interval $I:=[0,1]$ the Fabius function coincides with the cumulative distribution function (cdf) of
$$\sum_{j=1}^\infty 2^{-j}U_j,$$
where the $U_j$'s are independent random variables uniformly distributed on $I$. So, for each $x\in I$
$$F(x)=\lim_{n\to\infty} F_n(x),\tag{1}$$
where $F_n$ is the cdf of $\sum_{j=1}^n 2^{-j}U_j$.
Next (see e.g. formula (2.2)), for any real $x$
$\newcommand\vp{\varepsilon}$
\begin{equation}
F_n(x)=\text{vol}_n(I^n\cap H_{n;c^{(n)},x})
=\frac1{n!\prod_1^n c_i}\,\sum_{\vp\in\{0,1\}^n}(-1)^{|\vp|}\,
\big(x-c^{(n)}\cdot\vp\big)_+^n,
\end{equation}
where $\text{vol}_n$ is the Lebesgue measure on $\mathbb R^n$, $H_{n;b,x}:=\{v\in\mathbb R^n\colon b\cdot v\le x\}$, $c^{(n)}:=(c_1,\dots,c_n)$, $c_j:=2^{-j}$, $|\vp|:=\vp_1+\dots+\vp_n$, $\cdot$ denotes the dot product, and $t_+^n:=\max(0,t)^n$. So, for $x\in I$
\begin{equation}
F_n(x)=\frac{2^{n(n+1)/2}}{n!}\,\sum_{y\in D_{n,x}}(-1)^{s(y)}\,
\big(x-y\big)^n, \tag{2}
\end{equation}
where $D_{n,x}$ is the set of all dyadic numbers in $[0,x]$ of the form $m2^{-n}$ for integers $m$, and $s(y)$ is the sum of the binary digits of $y$.
So, this answer consists of formulas (1) and (2).
Some of the differences between (1)--(2) and the conjectured formula are as follows:
- The main conjectured formula for $F(x)$ is stated only for dyadic numbers $x$, and then extended to all values of $x$ by the continuity of $F$.
- The expression in the OP for $F(x)$ for dyadic $x$ contains a double summation and a number of $q$-Pochhammer symbols (I have had no experience with those symbols).
N.B.: this answer is cross-posted, with some edits, from MathOverflow.
Best Answer
As @Mathphile pointed out in the comments, the series can be written as $$S=1+\sum_{n=1}^{\infty} \frac{(2n-1)!!}{(2n-2)!!\cdot(2n)^n} \left( \frac{1}{n} - \frac{1}{(n+1)^2} \right)=1+\sum_{n=1}^{\infty} \frac{a_nb_n}{2^{3n-2}}$$ where $$a_n=\frac{(2n-1)!}{n^n\cdot(n-1)!^2}\quad\text{and}\quad b_n=\frac{1}{n} - \frac{1}{(n+1)^2},$$ since $(2n-1)!!=(2n - 1)! / [2^{n - 1} \cdot (n - 1)!]$ and $(2n-2)!!=2^{n - 1} \cdot (n - 1)!$. Clearly the terms of $b_n$ form a decreasing sequence, and this also holds for $a_n$ since $$\lim_{n\to\infty}\frac{a_{n+1}}{a_n}=\lim_{n\to\infty}\frac{2\left(\frac2n+\frac1{n^2}\right)}{1+\frac1n}\left(1-\frac1{n+1}\right)^n<1.$$ Thus $$S=1+\sum_{n=1}^{10}\frac{a_nb_n}{2^{3n-2}}+\sum_{m>10}\frac{a_{10}b_{10}}{2^{3m-2}}<1+0.414+10^{-14}<\sqrt2.$$