Evaluating $\int_{0}^{1}\cdots\int_{0}^{1}\bigl\{\frac{1}{x_{1}\cdots x_{n}}\bigr\}^{2}\:\mathrm{d}x_{1}\cdots\mathrm{d}x_{n}$

For any $\alpha > 0$ and $n \in \mathbb{Z}_{+}$. Let $I_n(\alpha)$ be the $n$-dimensional integral

$$I_n(\alpha) = \int_{[0,1]^n} \left\{1\left/\prod_{i=1}^n x_i\right.\right\}^\alpha \prod_{i=1}^n dx_i$$

It is clear all these $I_n(\alpha) \in (0,1)$. As a result, the corresponding generating function

$$I(\alpha,z) = \sum_{k=0}^\infty I_{k+1}(\alpha)z^k$$ converges absolutely for any $|z| < 1$.

We will first assume $z$ is small and positive.

For any given $n$, substitute $x_i$ by $e^{-t_i}$ for $i = 1,\ldots,n$. Let $t = \sum\limits_{i=1}^n t_i$ and $x = e^t$. We have

$$I_n(\alpha) = \frac{1}{(n-1)!}\int_0^\infty \left\{ e^t \right\}^\alpha e^{-t} t^{n-1} dt$$

and hence $$ I(\alpha,z) = \int_0^\infty \left\{ e^t \right\}^\alpha e^{(z-1)t} dt = \int_1^\infty \left\{ x \right\}^\alpha x^{z-2} dx = \sum_{k=1}^\infty \int_0^1 \frac{x^\alpha}{(x+k)^{2-z}} dx $$ In the last expression, if we replace the denominator $(x + k)^{2-z}$ by its integral representation

$$\frac{1}{(x+k)^{2-z}} = \frac{1}{\Gamma(2-z)}\int_0^\infty s^{1-z} e^{-(x+k)s} ds$$

We get

$$\begin{align} I(\alpha,z) &= \frac{1}{\Gamma(2-z)} \sum_{k=1}^\infty \int_0^1 x^\alpha \left[ \int_0^\infty s^{1-z} e^{-(x+k)s} ds\right] dx\\ &= \frac{1}{\Gamma(2-z)} \int_0^1 x^\alpha \left[ \int_0^\infty s^{1-z} \frac{e^{-xs}}{e^s - 1} ds \right] dx\\ &= \frac{1}{\Gamma(2-z)} \int_0^\infty \frac{s^{1-z}}{e^s - 1}\left[\int_0^1 x^\alpha e^{-xs} dx \right] ds\\ &= \frac{1}{\Gamma(2-z)} \int_0^\infty \frac{s^{1-z}}{e^s - 1}\left[\frac{1}{s^{\alpha+1}}\int_0^s x^\alpha e^{-x} dx \right] ds\tag{*1} \end{align} $$ Up to this step, we are using the assumption $z$ is small and positive. A consequence of this assumption is all the terms involved in above steps are non-negative numbers. The replacement of the denominator by its integral representation, the switching order of summation and integrations are automatically valid.

To proceed further, we need to split what's inside the square bracket in last expression of $(*1)$. If we do that, we will notice there are terms that are no longer non-negative. Furthermore, if $z$ remains to be small and positive, some of the terms simply diverge.

To bypass this obstacle, we will use the fact as long as $\Re z < 1$, the last expression of $(*1)$ defines an analytic function in $z$. Instead of sticking with $z$ is small and positive,

we will switch our assumption to $z$ is real and sufficiently negative

After we work out what the last expression of $(*1)$ really are, we will analytic continue the result back to small and positive $z$.

For integer $\alpha > 0$, we have

$$\begin{align} I(\alpha,z) &= \frac{1}{\Gamma(2-z)}\int_0^\infty \frac{s^{-(\alpha+z)}}{e^s - 1}\left[ -e^{-x} \sum_{k=0}^\alpha \frac{\alpha!}{(\alpha-k)!}x^{\alpha-k}\right]_0^s ds\\ &= \frac{1}{\Gamma(2-z)}\int_0^\infty \frac{s^{-(\alpha+z)}}{e^s - 1} \left[ \alpha! - e^{-s} \sum_{k=0}^\alpha \frac{\alpha!}{(\alpha-k)!}s^{\alpha-k} \right] ds\\ &= \frac{\alpha!}{\Gamma(2-z)}\int_0^\infty \left[s^{-(\alpha+z)} - \frac{1}{e^s-1}\sum_{k=0}^{\alpha-1}\frac{s^{-(k+z)}}{(\alpha-k)!} \right] e^{-s} ds \end{align} $$ Since $\displaystyle\; \int_0^\infty \frac{s^{\beta-1}}{e^s-1} e^{-s} ds = \Gamma(\beta)( \zeta(\beta) - 1) \;$, we find $$I(\alpha,z) = \frac{\alpha!}{\Gamma(2-z)} \left[\Gamma(1-\alpha-z) - \sum_{k=0}^{\alpha-1}\Gamma(1-k-z)\frac{\zeta(1-k-z)-1}{(\alpha-k)!} \right] $$

In particular, when $\alpha = 2$, this leads to $$ I(2,z) = \frac{2!}{\Gamma(2-z)} \left[ \Gamma(-1-z) -\frac12\Gamma(1-z)(\zeta(1-z)-1) -\Gamma(-z)(\zeta(-z)-1) \right] $$ Since $\;\Gamma(1+w) = w\Gamma(w)\;$, this can be simplified as $$\begin{align} I(2,z) &= \frac{2}{z(1-z^2)} - \frac{1}{1-z}(\zeta(1-z)-1) + \frac{2}{z(1-z)}(\zeta(-z)-1)\\ &= -\frac{1}{z+1} - \frac{1}{1-z}\zeta(1-z) + \frac{2}{z(1-z)}\zeta(-z)\tag{*2} \end{align} $$ We know $\zeta(w)$ is analytic over the whole complex plane except at $w = 1$. It has a simple pole with residue $1$ there. Using this info, we know RHS$(*2)$ may have poles at $z = 0$ and $\pm 1$.

• We don't care what happens at $z = 1$.
• Near $z = -1$, RHS$(*2)$ $\sim -\frac{1}{z+1} + O(1) + (-1 + O(z+1))(\frac{1}{-z-1} + O(1)) = O(1)$.
• Near $z = 0$, RHS$(*2)$ $\sim O(1) - (1 + O(z))(\frac{1}{-z} + O(1)) + \frac{1}{z}(2\zeta(0) + O(z)) = O(1)$ again.

Combine these, we know RHS$(*2)$ is analytic at $z = -1$ and $0$. There is nothing from extending its validity to the unit disk $|z| < 1$ where $I(\alpha,z)$ coincides with its definition as a generating function. As a result,

$$\bbox[4pt,border:1px solid blue;]{ I_n(2) = (n-1)!\frac{d^{n-1}}{dz^{n-1}}\left[-\frac{1}{z+1} - \frac{1}{1-z}\zeta(1-z) + \frac{2}{z(1-z)}\zeta(-z) \right]_{z=0} }\tag{*3}$$

If we throw RHS$(*2)$ to WA and ask for its Taylor expansion, WA returns

$$I(2,z) = (-1-\gamma + \log(2\pi)) + z \left(1 - \gamma+ \frac{\gamma^2}{2}-\frac{\pi^2}{24} + \log(2\pi) - \frac{\log^2(2\pi)}{2}\right) + O(z^2)$$

This implies $$ \begin{align} I_1(2) &= -1-\gamma + \log(2\pi)\\ I_2(2) &= 1 - \gamma+ \frac{\gamma^2}{2}-\frac{\pi^2}{24} + \log(2\pi) - \frac{\log^2(2\pi)}{2} \end{align}$$ as expected.

Here is an approach.

Theorem 1. Let $s$ be a complex number such that $-1<\Re{s}<1$, $s\neq 0$.

Then $$ \int_{0}^{1} x^{s}\left\{1/x\right\}^{2}\:\mathrm{d}x = -\frac{2\zeta(s)}{s(1+s)}-\frac{\zeta(1+s)}{1+s}-\frac{1}{1-s} \tag1 $$ where $\left\{x\right\}=x-\lfloor x\rfloor$ denotes the fractional part of $x$ and where $\zeta$ denotes the Riemann zeta function.

Proof. Let us assume first that $\Re{s}>1$. We have

\begin{align} \displaystyle \int_{0}^{1} x^{s}\left\{1/x\right\}^{2}\:\mathrm{d}x & = \sum_{k=1}^{\infty} \displaystyle \int_{1/(k+1)}^{1/k} x^{s}\left\{1/x\right\}^{2} \:\mathrm{d}x \\ & = \sum_{k=1}^{\infty} \displaystyle \int_{k}^{k+1} (x-k)^{2} \frac{\mathrm{d}x}{x^{s+2}} \\ & = \sum_{k=1}^{\infty} \displaystyle \int_{0}^{1}\frac{x^{2}}{(x+k)^{s+2}}\mathrm{d}x \\ & = \sum_{k=1}^{\infty} \displaystyle \int_{0}^{1}\left(\frac{1}{(x+k)^{s}}-\frac{2k}{(x+k)^{s+1}}+ \frac{k^{2}}{(x+k)^{s+2}}\right)\mathrm{d}x \\ & = \sum_{k=1}^{\infty} \displaystyle \left. \left( -\frac{1}{(s-1)} \frac{1}{(x+k)^{s-1}} +\frac{2}{s}\frac{k}{(x+k)^{s}} -\frac{k^{2}}{(s+1)(x+k)^{s+1}} \right) \right|_{0}^{1}\\ & = -\frac{1}{1-s} -\frac{2\zeta(s)}{s(1+s)}-\frac{\zeta(1+s)}{1+s}. \end{align} Since the integral on the left hand side is clearly analytic for $\Re s >-1$ and since all singularities on the right hand side (including $s=0$) are removable, we may extend the preceding identity by analytic continuation to obtain $(1)$.

$\,\Box$

A first application of $(1)$.

Example. $$ I_{2}= \int_{0}^{1}\!\!\int_{0}^{1} \left\{\frac{1}{x y}\right\}^{2} \:\mathrm{d}x \mathrm{d}y = 1-\gamma+\dfrac{\gamma^2}{2}-\dfrac{\pi^2}{24}+\ln(2\pi)-\dfrac{\ln^2(2\pi)}{2} \tag2 $$ where $\left\{x\right\}=x-\lfloor x\rfloor$ denotes the fractional part of $x$ and $\gamma$ is the Euler constant.

Proof. We use a change of variable in the given integral: $$ \begin{align} \displaystyle u & = xy \\ v & = x \\ \end{align} $$ then the absolute value of the Jacobian determinant and the new domain of integration are $$ \left| \frac{D(u,v)}{D(x,y)}\right|=v, \quad 0<u<v<1, $$ leading to $$ \begin{align} I_{2} & = \int_{0}^{1}\!\!\int_{0}^{1} \left\{\frac{1}{x y}\right\}^{2} \:\mathrm{d}x \mathrm{d}y \\\\ & = \int_{0}^{1} \! \int_{u}^{1} \left\{1/u\right\}^2\frac{\mathrm{d}v}{v} \mathrm{d}u \\\\ & =-\int_{0}^{1} \left\{1/u\right\}^{2} \ln u\: \mathrm{d}u \\\\ & =-\left.\frac{d}{ds}\left(\int_{0}^{1} u^s \left\{1/u\right\}^{2} \mathrm{d}u\right) \right|_{s=0} \\\\ & =-\left.\frac{d}{ds}\left(-\frac{1}{1-s} -\frac{2\zeta(s)}{s(1+s)}-\frac{\zeta(1+s)}{1+s}\right) \right|_{s=0} \\\\ & = 1-\gamma+\dfrac{\gamma^2}{2}-\dfrac{\pi^2}{24}+\ln(2\pi)-\dfrac{\ln^2(2\pi)}{2} \end{align} $$ where, for $s$ near $0$, we have used $$ \begin{align} & \zeta(1+s) =\frac1s+\gamma-\gamma_1 s +\mathcal{O}(s^2), \\& \zeta(s) =-\frac12-\dfrac{\ln(2\pi)}{2} s +\left(\dfrac{\gamma^2}{4}-\dfrac{\pi^2}{48}+\ln(2\pi)-\dfrac{\ln^2(2\pi)}{4}+\dfrac{\gamma_1}{2}\right)s^2+\mathcal{O}(s^3). \end{align} $$ $\Box$

Now one may observe that $(1)$ has the following consequence.

Theorem 2. Let $n=1,2,3, \cdots.$

Then $$ \frac{(-1)^{n-1}}{(n-1)!} \int_{0}^{1} \left\{1/x\right\}^{2} \ln ^{n-1} x\: \mathrm{d}x = (-1)^{n}+1 -\sum_{k=0}^{n-1}\frac{\gamma_{k}}{k!}+2\sum_{k=0}^{n} \frac{(-1)^k}{k!}\zeta^{(k)}(0) \tag3 $$ where $\left\{x\right\}=x-\lfloor x\rfloor$ denotes the fractional part of $x$ and where $\gamma_{k}$ denote the Laurent-Stieltjes constants which may be defined by $$ \gamma_{k} = \displaystyle \lim_{N \rightarrow \infty}\left(\sum_{m=1}^{N}\frac{\ln^{k}m}{m}-\frac{\ln^{k+1}N}{k+1} \right). $$

Proof. Let $s$ be a complex number such that $0<\Re{s}<1$. Using Theorem 1, one may just write $$ \int_{0}^{1} x^{s}\left\{1/x\right\}^2 \:\mathrm{d}x = \sum_{n=0}^{\infty}\left(\int_{0}^{1} \left\{1/x\right\}^2 \ln^{n}x \:\mathrm{d}x \right) \frac{s^{n}}{n!}, $$ and one may take into account the standard Laurent series $$\begin{align} \displaystyle \zeta(1+s) =\frac{1}{s}+\sum_{k=0}^{\infty} (-1)^{k} \frac{\gamma_{k}}{k!} s^{k} \end{align}$$ to obtain $(3)$ by identifying the Taylor coefficients on both sides of $(1)$.

$\,\Box$

From $(3)$ we deduce the following result.

Theorem 3. Let $n=1,2,3, \cdots.$ Set $$ I_{n}:= \int_{0}^{1} \! \cdots \! \int_{0}^{1} \left\{\frac{1}{x_{1}x_{2} \cdots x_{n}}\right\}^{2} \:\mathrm{d}x_{1}\mathrm{d}x_{2} \cdots \mathrm{d}x_{n}. $$ Then $$ I_{n}= (-1)^{n}+1 -\sum_{k=0}^{n-1}\frac{\gamma_{k}}{k!}+2\sum_{k=0}^{n} \frac{(-1)^k}{k!}\zeta^{(k)}(0) \tag4 $$ where $\left\{x\right\}=x-\lfloor x\rfloor$ denotes the fractional part of $x$ and where $\gamma_{k}$ denote the Laurent-Stieltjes constants defined by $$ \gamma_{k} = \displaystyle \lim_{N \rightarrow \infty}\left(\sum_{m=1}^{N}\frac{\ln^{k}m}{m}-\frac{\ln^{k+1}N}{k+1} \right). $$

Proof. We reduce $ I_{n} $ to a one dimensional integral by making the change of variables $$ \begin{align} \displaystyle u_{1} & = x_{1} \\ u_{2} & = x_{1}x_{2} \\ & \vdots \\ u \, & = x_{1}x_{2} \cdots x_{n} . \end{align} $$ The initial domain of integration is transformed into $ \displaystyle 0<u< \cdots <u_{2}< u_{1}<1 $ with a Jacobian equal to $\displaystyle u_{1} \cdots u_{n-1}$. Hence $$ \begin{align} I_{n} = \displaystyle \int_{0}^{1} \! \int_{u}^{1} \cdots \! \int_{u_{2}}^{1} \left\{1/u\right\}^{2} \: \frac{\mathrm{d}u_{1} \cdots \mathrm{d}u_{n-1}}{u_{1} \cdots u_{n-1}} \mathrm{d}u \end{align} $$ then integrating out successively $\displaystyle u_{1}, u_{2}, \dots ,u_{n-1}$ gives $$ I_{n} = \frac{(-1)^{n-1}}{(n-1)!}\displaystyle \int_{0}^{1} \left\{1/u\right\}^{2} \ln ^{n-1} u\: \mathrm{d}u $$ and we are directly led to $(4)$ by $(3)$.

$\Box$

Here is a general formula answering an open problem from Ovidiu Furdui (2006) on this class of multiple integrals involving the fractional part.

Theorem (O. Oloa, general case, 18th Aug 2014). Let $n=1,2,3...$ and let $p=1,2,3....$

Set $$ \begin{align} I_{n,p} :=\int_{0}^{1} \! \cdots \! \int_{0}^{1} \left\{\frac{1}{x_{1} \cdots x_{n}}\right\}^{p} \:\mathrm{d}x_{1} \cdots \mathrm{d}x_{n} \end{align} $$ where $\left\{x\right\}=x-\lfloor x\rfloor$ denotes the fractional part of $x$.

Then $$ \begin{align} \displaystyle I_{n,1} :=\int_{0}^{1} \! \cdots \! \int_{0}^{1} \left\{\frac{1}{x_{1} \cdots x_{n}}\right\} \:\mathrm{d}x_{1} \cdots \mathrm{d}x_{n} = 1-\sum_{k=0}^{n-1}\frac{\gamma_{k}}{k!}, \end{align} $$ and, for $p=2,3,4,...,$ $$ \begin{align} I_{n,p} = \frac{(-1)^{n}}{(p-1)^{n}}- 1-\sum_{k=0}^{n-1}\frac{\gamma_{k}}{k!} - \frac{p!}{n!}\sum_{\ell=0}^{p-2}\frac{(-1)^{\ell}}{(p-\ell-1)!}\left.\left(\! \frac{\zeta(s-\ell)}{(1+s)(1-s)_{\ell}}\! \right)^{(n)}\right|_{s=0} \end{align} $$ where $\gamma_{k}$ are the Laurent Stieltjes constants and where $\zeta$ denotes the Riemann zeta function.