Evaluation of an $n$-dimensional integral I asked the same question on math.se but got no answer there. Since it pertains to my current research, I decided to ask here:
Let $n\in 2\mathbb{N}$ be an even number. I want to evaluate
$$I_n
:=
\int_0^1\mathrm{d} u_1 \cdots \int_0^1 \mathrm{d} u_n \frac{\delta(1-u_1-\cdots-u_n)}{(u_1+u_2)(u_2+u_3)\cdots(u_{n-1}+u_n)(u_n+u_1)}.
$$
For small $n$, this is computable by simply parameterizing the $\delta$ function, and
$I_2 = 1$, $I_3 = \pi^2/4$, $I_4 = 2\pi^2/3$. The values of $I_5$ and $I_6$ are numerically $18.2642 \approx 3\pi^4/16$ and $51.9325\approx 8\pi^4/15$. I strongly suspect that
$$
I_{2n+2} \stackrel{?}{=} (2\pi)^{2n} \frac{(n!)^2}{(2n+1)!} =
\frac{(2\pi)^{2n}}{\binom{2n+1}{n}(n+1)}
= (2\pi)^{2n}\mathrm{B}(n+1,n+1),
$$
where $\mathrm{B}$ is the Beta function. Dividing by $(2\pi)^{2n}$, this is Sloane's A002457. For $I_6$, this conjecture is equivalent to
$$\int_0^1\mathrm{d}x \Bigl(\mathrm{Li}_2(\frac{x-1}{x})\Bigr)^2 \stackrel{?}{=} \frac{17}{180}\pi^4$$
(with $\mathrm{Li}_2$ the dilogarithm), which seems to be true numerically, but I could neither prove it nor find it in the literature.
As a last remark, it is possible to get rid of the $\delta$ function by using the identity
$$I_n = \int_{(0,\infty)^n}\mathrm{d}u \frac{f(\lvert u\rvert_1)}{(u_1+u_2)\cdots(u_n+u_1)} \Bigm/\int_0^\infty\mathrm{d}t \frac{f(t)}{t}$$
for any $f:(0,\infty)\to\mathbb{R}$ that makes both integrals finite. Using $f(t) = t 1_{[0,1]}(t)$ where $1_{[0,1]}$ is the characteristic function of the interval $[0,1]$, one can write $I_n$ as an integral over an $n$-dimensional simplex.
 A: Here is another proof, the main part of which was communicated to me by Dr. Peter Otte of Bochum University:
\begin{equation}
I_n := \int_{[0,1]^n}\mathrm{d}u\,\delta(1-\lvert u\rvert_1)
\frac{1}{\prod_{j=1}^n (u_j + u_{j+1})}
=
(2\pi)^{n-2}
\frac{[\Gamma(\frac{n}{2})]^2}{\Gamma(n)}.
\end{equation}
First, define
$$J_n(t) := \int_{[0,1]^n}\mathrm{d}u\,\delta(t-\lvert u\rvert_1)
           \frac{1}{\prod_{j=1}^{n-1}(u_j + u_{j+1})}.$$
for $t>0$. By scaling, $J_n(t) = J_n(1) =: J_n$ for all $t > 0$. Also,
    \begin{align}
      I_n &
      = \frac{1}{2}\int_{[0,1]^n}\mathrm{d}u\, \delta(1-\lvert u\rvert_1)
      \frac{2\lvert u\rvert_1}{\prod_{j=1}^n (u_j + u_{j+1})} \notag\\
      & =
      \frac{1}{2}\sum_{k=1}^n \int_{[0,1]^n}\mathrm{d}u\, \delta(1-\lvert u\rvert_1)
      \frac{u_k+u_{k+1}}{\prod_{j=1}^n (u_j + u_{j+1})} \notag\\
      & = \frac{n}{2} \int_{[0,1]^n}\mathrm{d}u\,
      \frac{\delta(1-\lvert u\rvert_1)}{(u_1+u_2)\dotsm(u_{n-1}+u_n)} =
      \frac{n}{2} J_n.
    \end{align}
Next, let $f\in L_1(0,\infty)$. Then
  \begin{equation}
    J_n
    =
    \int_{(0,\infty)^n}\mathrm{d}u\, \frac{f(\lvert u\rvert_1)}{\prod_{j=1}^{n-1}(u_j + u_{j+1})}
    \Bigm/\!
    \int_0^\infty\mathrm{d}t\, f(t).
  \end{equation}
  In particular,
  \begin{equation}
    J_n = \int_{(0,\infty)^n}\mathrm{d}u\, \frac{e^{-\lvert u\rvert_1}}{\prod_{j=1}^{n-1}(u_j + u_{j+1})},
  \end{equation}
We will need the Rosenblum-Rovnyak integral operator
$T: L_2(0,\infty)\to L_2(0,\infty)$,
see Rosenblum (1958) and Rovnyak (1970), defined via
\begin{equation}
  (Tf)(x) := \int_0^\infty \mathrm{d}y\, \frac{e^{-(x+y)/2}}{x+y} f(y)
  \quad (x\in(0,\infty)).
\end{equation}
for $f\in L_2(0,\infty)$. This is the special case $T = \mathcal{H}_0$ in
Rosenblum (1958), Formula (2.3).
The operator $T$ is unitary equivalent to the Hilbert matrix
$H:\ell_2(\mathbb{N})\to\ell_2(\mathbb{N})$,
\begin{equation}
  (H x)_j = \sum_{k=1}^\infty \frac{x_k}{j+k-1} \quad(j\in\mathbb{N}, x\in\ell_2(\mathbb{N}))
\end{equation}
and can be explicitly diagonalized: Following
Yafaev (2010), Sec. 4.2, we define the unitary operator
$U: L_2(0,\infty)\to L_2(0,\infty)$ via
\begin{equation}
  (Uf)(k) = \pi^{-1}\sqrt{k\sinh 2\pi k} \, \lvert \Gamma(1/2 - ik)\rvert
            \int_0^\infty\mathrm{d}x\, x^{-1} W_{0,ik}(x)f(x)
\end{equation}
for $f\in L_2(0,\infty)$ and $k\in(0,\infty)$, where the Whittaker
functions are given by
\begin{equation}
  W_{0,\nu}(x) = \sqrt{x/\pi} K_\nu(x/2) \quad (\nu, x\in(0,\infty)),
\end{equation}
with $K_\nu$ as the modified Bessel function of the second kind, see
DLMF.
In order to compute $J_n$, we will employ the following result due to
Rosenblum, see Yafaev, Prop. 4.1:
\begin{equation}
 (UTf)(k) = \frac{\pi}{\cosh(k\pi)}(Uf)(k)  \quad (k\in(0,\infty),
  f\in L_2(0,\infty).
\end{equation}
Proof of $I_n = (2\pi)^{n-2} \frac{[\Gamma(\frac{n}{2})]^2}{\Gamma(n)}$.
  Let $n\in\mathbb{N}_{\ge 2}$. From the definition of $T$ and the
  identity of $J_n$ above, we see that
  \begin{equation}
    J_n = \langle f_0, T^{n-1}f_0\rangle
  \end{equation}
  with $f_0(x) := e^{-x/2}$.
  From this and the identity of $UT$ above, we obtain
  \begin{equation}
    J_n =
    \langle Uf_0, UT^{n-1}f_0\rangle
    =
    \int_0^\infty\mathrm{d}k\, \lvert \hat{f}_0(k)\rvert^2 \Bigl(\frac{\pi}{\cosh(k\pi)}\Bigr)^{n-1},
  \end{equation}
  where $\hat{f}_0 := Uf_0$. In order to compute $\hat{f}_0$, we
  employ the classical formula
  \begin{equation}
    \lvert\Gamma(1/2 - ik)\rvert^2 = \frac{\pi}{\cosh(k\pi)}  \quad (k\in\mathbb{R}),
  \end{equation}
  which is a consequence of the reflection formula for the Gamma
  function, and
  \begin{equation}
    \int_0^\infty\mathrm{d}x\, x^{-1} W_{0,ik}(x)e^{-x/2} =
    \frac{\pi}{\cosh(k\pi)}
    \quad(k > 0),
  \end{equation}
  which follows from the special case $z=1/2$ and $\nu = \kappa = 0$ in
  DLMF. From the definition of $U$ above and the last two equations,
  we deduce
  \begin{equation}
    \lvert\hat{f}_0(k)\rvert^2 = 2\pi k\frac{\sinh(k\pi)}{\cosh(k\pi)^2}
    \quad (k > 0).
  \end{equation}
  This yields
  \begin{equation}
    J_n = 2\pi^{n-2}\int_0^\infty\mathrm{d}k\, k
    \frac{\sinh(k)}{\cosh(k)^{n+1}}
    =
    \frac{2\pi^{n-2}}{n}\int_0^\infty\mathrm{d}k\,\frac{1}{\cosh(k)^n}
  \end{equation}
  where we applied the substitution $\tilde{k} = k\pi$ and integrated
  by parts. This integral can be evaluated using the substitutions
  $y = \cosh(k)^{-1}$ and $x = y^2$, one after the other:
  \begin{align}
    J_n
    =
    \frac{2\pi^{n-2}}{n}
    \int_0^1\mathrm{d}y\, \frac{y^{n-1}}{\sqrt{1-y^2}}
    & =
    \frac{\pi^{n-2}}{n}
    \int_0^1\mathrm{d}x\, x^{n/2-1}(1-x)^{-1/2} \\
    & =
    \frac{\pi^{n-2}}{n} \mathrm{B}(n/2, 1/2),
  \end{align}
  since $k'(y) = - y^{-1}(1-y^2)^{-1/2}$.
  The claim then follows by expressing the Beta function via the Gamma
  function and then applying the classical duplication formula.
