I found the following result on this webpage:
$$\int_0^{\infty } \left(\prod _{k=0}^7 \frac{\sin \left(\frac{x}{2 k+1}\right)}{\frac{x}{2 k+1}}\right) \, \mathbb{d}x= \frac{\pi}{2} - \frac{6879714958723010531}{935615849440640907310521750000} \pi $$
However, I can't determine how to prove it.
I have a visual explanation, why its below $\frac{\pi}{2}$
In 2001, the Borwein Brothers surprised their readers with the trick you referred to:
The pattern suggests that the definite integral of the sinc function ($=\pi/2$), does not alter when this sinc is multiplied with other sincs having lower frequencies. But unexpectedly, they show that the eighth integral is about one billionth of a percent smaller. How is that possible?
French mathematician Fourier allows us to look behind the scenes. From his perspective, we see a rectangular pedestal that gradually erodes: although in the beginning only its corners are rounded, inevitably the middle of the plateau will be affected as well, thereby lowering its actual height. To understand this, let's first recollect some Fourier theory:
Five familiar Fourier facts
The Fourier transform of a sinc, is a rectangular function (SINC-RECT pair). For the Borwein Integrals, it is instructive to see each sinc-function as the frequency spectrum of a rectangular signal.
The 'Uncertainty Principle' tells us: 'the wider the sinc, the narrower the block'. Thus the decreasing frequency of each new sinc factor, translates in a decreasing width of its corresponding rectangular pulse.
The Convolution Theorem tells us: a multiplication of spectra is equivalent to the convolution of two signals. Therefore, we can understand the multiplication of these sincs, as a convolution of their corresponding block pulses.
The inverse Fourier Transform tells us, that the integral of a spectrum is the central value in the time domain (similar as the integral of a signal, is the central component in the Frequency domain). Therefore these Borwein Integrals are proportional to the central value of the convolution outcome of the corresponding block pulses.
The convolution of two equivalent block pulses (equal width), results in a triangular function (an isosceles): the left and right corner of the rectangle are completely eroded (see animation below), but the center of the plateau is preserved (which btw explains why $\int_0^{\infty} \operatorname{sinc}(x)^2 dx = \frac{\pi}{2}$ (see: integralOfSquaredSincs))
Source: Convolution
Convolving a wide rectangle with a narrow one
With this theory in mind, let's look at the second Borwein Integral. The product of the first two sincs, is equivalent to convolving a rectangle having unit width (image below, left lower graph) with a (1/3)-kernel (upper graph in the middle). The result is a trapezoid, in which only two third of the initial unit plateau width is preserved, and one third is eroded:
Why? Because in the interval [-1/3, 1/3], the (1/3)-kernel was completely embedded within the unit block. Just outside this domain, however, an increasing part of the kernel is not 'matched', resulting in two downward slopes moving out of the center. However, since the center of the initial plateau is not (yet) affected, the second Borwein integral equals still $\frac{\pi}{2}$.
The result after some more convolutions:
In the third Borwein integral, the resulting trapezoid is convolved with an (1/5)-kernel, which takes away another 20% of the original plateau width. Therefore the new residual plateau is reduced to 1 - 1/3 - 1/5 = 7/15 part of the initial plateau width:
With each new convolution, the residual plateau gets smaller and smaller:
*(which is not that well visible, because each convolution also 'smoothens' the outcome, exactly how the probability distribution of a coin flipping experiment gets smoother when more coins are flipped)*
Not enough plateau left to completely match the (1/15)-kernel
Once the outcome of the (1/11)-convolution in turn is convolved with the (1/13)-kernel, the residual plateau equals: $1 - \frac{1}{3}-\frac{1}{5}-...-\frac{1}{13}= \frac{2021}{45045}<\frac{1}{15}$ which is for the first time smaller than the next convolution kernel.
So for the eighth Borwein integral, the remaining plateau is too small to completely match the (1/15)-kernel in the center, resulting in a slightly lower integration outcome for the main frequency component.
And since this main frequency component is proportional to this eighth integral:
$\int_0^\infty \frac{\sin(x)}{x}\frac{\sin(x/3)}{x/3}\cdots\frac{\sin(x/15)}{x/15}dx$
we now see why it's for the first time slightly below $\frac{\pi}{2}$.
Job Bouwman's answer in terms of convolutions can be framed nicely in terms of probability.
In that answer, $\mathop{\rm sinc}(t) = \frac{\sin(t)}{t}$ plays the role of the Fourier transform of the indicator function of an interval. This indicator function can be interpreted as the probability density function of a uniformly distributed random variable. Then its Fourier transform is the characteristic function of that random variable. Convolution of probability density functions corresponds to addition of independent random variables, which in turn corresponds to multiplication of their characteristic functions.
Specifically, let $X_0, X_1, X_2, \dots$ be independent and uniformly distributed in $[-1,1]$. As in Job Bouwman's answer, we care about the value at $x=0$ of the convolution of the density functions of $X_0$, $X_1/3$, $X_2/5$, and so on. So we're interested in the value at $x=0$ of the probability density function of
$$X_0 + \frac13 X_1 + \frac15 X_2 + \dots + \frac1{2N+1} X_N.$$
That is, how likely is that sum to land in a small interval centered at $0$, compared to the size of that small interval? Up to a factor of 2, this is the same as asking how likely the sum
$$Y = \frac13 X_1 + \frac15 X_2 + \dots + \frac1{2N+1} X_N$$
is to land in the interval $[-1, 1].$ Since each $X_i$ has absolute value at most $1$, this happens with probability $1$ when $N<7$, but with probability slightly less than $1$ when $N \geq 7$.
Explicitly, let $f$ denote the probability density function of $Y$ and $\varphi$ denote the characteristic function of $Y$. We have
\begin{align*} \mathop{\mathbb P}(Y \in [-1,1]) &= \int_{-1}^1 f(x)\,dx\\ &= \int_{\mathbb R} \chi_{[-1,1]}(x)\, f(x)\,dx\\ &= \frac1{2\pi} \int_{\mathbb R} 2 \mathop{\rm sinc}(t)\, \varphi(t)\,dt, \end{align*} where the last equality is due to the fact that Fourier transforms preserve the inner product between functions, up to a factor of $2\pi$.
We also have
$$\varphi(t) = \mathop{\mathbb E}(e^{itY}) = \prod_{k=1}^N \mathop{\mathbb E}(e^{i\,t/(2k+1)\,X_k}) = \prod_{k=1}^N \mathop{\rm sinc}\left(\frac{t}{2k+1}\right).$$
Hence
\begin{align*} \mathop{\mathbb P}(Y \in [-1,1]) &= \frac1{2\pi} \int_{\mathbb R} 2 \mathop{\rm sinc}(t)\, \prod_{k=1}^N \mathop{\rm sinc}\left(\frac{t}{2k+1}\right)\,dt\\ &= \frac2{\pi} \int_0^\infty \prod_{k=0}^N \mathop{\rm sinc}\left(\frac{t}{2k+1}\right)\,dt\\ &= \frac2{\pi} B_N, \end{align*}
where $B_N$ is the integral we're interested in. So $B_N = \pi/2$ when $N<7$ and $B_n$ is just a little bit less than $\pi/2$ when $N \geq 7$. Note that even if you take $N$ to infinity, the probability that
$$Y = \frac13 X_1 + \frac15 X_2 + \frac17 X_3 + \dots$$
is outside $[-1, 1]$ is still quite small, and
$$B_\infty = \int_0^\infty \prod_{k=0}^\infty \mathop{\rm sinc}\left(\frac{t}{2k+1}\right)\,dt = \int_0^\infty \prod_{j=1}^\infty \cos\left(\frac{t/2}{j}\right)\,dt = 2\int_0^\infty \prod_{j=1}^\infty \cos\left(\frac{u}{j}\right)\,du $$
is still quite close to $\pi/2$, correct to four digits I believe.
You can observe an even more striking version of this behavior by evaluating the probability density function of
$$X_0 + \frac13 X_1 + \frac15 X_2 + \dots + \frac1{2N+1} X_N$$
at $x = 1$ rather than $x = 0$. Up to a factor of $2$, this is the probability that
$$Y = \frac13 X_1 + \frac15 X_2 + \dots + \frac1{2N+1} X_N$$
lands in the interval $[0,2]$. By symmetry around the origin, when $N < 56$, this probability is exactly $1/2$, but when $N \geq 56$, it is very slightly smaller than $1/2$. Translating this into facts about integrals, we get that
$$C_N = \frac12 \int_{\mathbb R} e^{it} \prod_{k=0}^N \mathop{\rm sinc}\left(\frac{t}{2k+1}\right)\,dt = \int_0^\infty \cos(t)\, \prod_{k=0}^N \mathop{\rm sinc}\left(\frac{t}{2k+1}\right)\,dt$$
is exactly $\pi/4$ when $N < 56$ and very slightly smaller when $N \geq 56$.
And because the probability that
$$\frac13 X_1 + \frac15 X_2 + \frac17 X_3 + \dots$$
falls outside $[-2, 2]$ is miniscule but nonzero, we have that
$$C_\infty = \int_0^\infty \cos(t)\,\prod_{k=0}^\infty \mathop{\rm sinc}\left(\frac{t}{2k+1}\right)\,dt = 2\int_0^\infty \cos(2u) \prod_{j=1}^\infty \cos\left(\frac{u}{j}\right)\,du $$
is incredibly close to $\pi/4$, correct to more than forty digits!
As for the original question of calculating $B_7$ exactly, it boils down to calculating the exact probability $p$ that
$$\frac13 X_1 + \frac15 X_2 + \frac17 X_3 + \frac19 X_4 + \frac1{11} X_5 + \frac1{13} X_6 + \frac1{15} X_7 > 1.$$
Up to a factor of $2^7$, this is the volume of a tetrahedral section of the 7-dimensional cube $[-1, 1]^7$, sliced off by the plane $x_1 / 3 + \dots + x_7/15 = 1$. This tetrahedron has orthogonal sides of length $3s, 5s, \dots, 15s$, where $$s = \left(\frac13 + \frac15 + \dots + \frac1{15}\right) - 1 = \frac{44190}{15!!}.$$ Hence,
$$p = \frac{1}{2^7} \frac{1}{7!} 3s \dots 15s = \frac{15!! s^7}{2^7 7!} = \frac{44190^7}{2^7 7! (15!!)^6}$$
and
$$B_7 = \frac\pi2 (1-2p) = \frac\pi2\left(1 - \frac{44190^7}{2^6 7! (15!!)^6} \right).$$
The great reference is:
Borwein, David; Borwein, Jonathan M. "Some remarkable properties of sinc and related integrals." Ramanujan J. 5 (2001) 73–89.
DOI 10.1023/A:1011497229317
$\newcommand{\S}{\operatorname{sinc}}$
A 'brute force' solution to a beautiful problem. I don't claim this answer is as insightful as the others and the question is somewhat old, but I feel it is relevant and unique enough to merit posting.
As usual, we will let $\S(x) = \sin(x)/x$ with $\S(0)=1$.
The big idea: use angle-addition, Taylor series, and integration by parts to 'revert' the integral to a bunch- $128$, to be precise- of weighted integrals of the form $\int_0^{\infty}\sin(a x)/x\,dx$, $a>0$, each of which is $\pi/2$. Then we just add up the weights to produce the answer.
Recall that \begin{align} 2\sin(\theta)\sin(\phi)=\cos(\theta-\phi)-\cos(\theta+\phi)\\ 2\sin(\theta)\cos(\phi)=\sin(\theta-\phi)-\sin(\theta+\phi)\\ 2\cos(\theta)\cos(\phi)=\cos(\theta-\phi)+\cos(\theta+\phi)\\ \end{align}
For instance, this tells us $\S(x)\S(x/3) = \frac{3}{2}x^{-2}(\cos(2x/3)-\cos(4x/3))$. You probably see what lies ahead, even if it's somewhat unappetizing. We have $$ \S(x)\S(x/3)\cdots \S(x/15) = $$ $$ \frac{1\cdot 3\cdots 15}{2^7 x^8} \sum_{e_i\in\{\pm 1\}} (-1)^{\#(e_i=-1)} \cos\left( x(1+e_1 \cdot 1/3 + e_2\cdot 1/5 +\cdots +e_7 \cdot 1/15)\right) $$Thus, if we denote $W=(1\cdot 3\cdot \cdots \cdot 15)/2^7$, we have to evaluate $128$ integrals that look like $\displaystyle{W \int_0^{\infty}\frac{\cos(a x)}{x^8}\,dx}$, for some constant $a$. The trouble is that none of them are improperly integrable near $x=0$, which brings us to the second step in the process.
The idea here is to add and subtract the same number of lower-order terms of the Taylor series of each $\cos(a x)$ to make the integrals converge. For example, $\displaystyle{\int_0^{\infty}\frac{\cos(x)}{x^8}\,dx}$ does not converge but $\displaystyle{\int_0^{\infty}\frac{\cos(x)-(1-x^2/2+x^4/24-x^6/720)}{x^8}\,dx}$ does. Since the sum of the arguments of the cosines is zero (check! A short explanation is that the $+$ and $-$ signs are equal in number) we will add and subtract up to the degree six terms to make the integrals convergent. Therefore, we have to evaluate $$ \int_0^{\infty}\frac{\cos(a x)-(1-(ax)^2/2+(ax)^4/24-(ax)^6/720)}{x^8}\,dx $$
We have $$ \int_0^{\infty}\frac{\cos(a x)-(1-(ax)^2/2+(ax)^4/24-(ax)^6/720)}{x^8}\,dx $$ \begin{align} &=\left.-\frac{\cos(a x)-(1-(ax)^2/2+(ax)^4/24-(ax)^6/720)}{7x^7}\right|_0^{\infty}\\ &+\frac{1}{7}\int _0^{\infty} \frac{-a\sin(ax)+a^6 x^5/120 - a^4 x^3/6 + a^2 x }{x^7}\,dx \end{align}The boundary terms vanish: at $\infty$, the denominator is a greater power and at $0$ the cosine series is $O(x^8)$. Continue this procedure until the denominator is just $x$: $$ \cdots = \frac{a^7}{7!}\int_0^{\infty}\frac{\sin(ax)}{x}\,dx = \frac{a^7}{7!}\cdot\frac{\pi}{2} $$
All that's left is to add up these values for all $128$ values of $a$ corresponding to the $2^7$ terms of the sum. This becomes $$ W \cdot \frac{\pi}{2}\cdot \frac{1}{7!} \cdot \sum_{e_i\in\{\pm 1\}} (-1)^{\#(e_i=-1)} (1+e_1 \cdot 1/3 + e_2\cdot 1/5 +\cdots +e_7 \cdot 1/15)^7 $$A somewhat tedious but routine calculation gives the final answer of $$ \frac{467807924713440738696537864469 \pi}{935615849440640907310521750000} $$
For clarity, I'll do a smaller example (the three sinc case) that will evaluate to $\pi/2$. $$ \int _0^{\infty} \S(x)\S(x/3)\S(x/5)\,dx = 15 \int_0^{\infty} \frac{\sin(x)\sin(x/3)\sin(x/5)}{x^3}\,dx $$ $$ =\frac{15}{4} \int_0^{\infty} \frac{-\sin(7/15 x)+\sin(13/15 x)+\sin(17/15 x)-\sin(23/15 x)}{x^3}\,dx $$ $$ =\frac{15}{4} \int_0^{\infty} \frac{-\sin(7/15 x)+\sin(13/15 x)+\sin(17/15 x)-\sin(23/15 x) }{x^3}+\frac{7/15 x-13/15 x-17/15 x+23/15 x}{x^3}\,dx $$I'll just pick one of these, say the first one; the others are exactly the same. $$ =\int_0^{\infty} \frac{-\sin(7/15 x)+7/15 x}{x^3}\,dx $$ $$ = \left.\frac{-\sin(7/15 x)+7/15 x}{-2x^2}\right|_0^{\infty} + \int_0^{\infty} \frac{-7/15 \cos(7/15 x)+7/15}{2x^2}\,dx $$ $$ =0+0 + \frac{(7/15)^2}{2} \int _0^{\infty}\frac{\sin(7/15 x)}{x}\,dx =\frac{(7/15)^2}{2} \cdot \frac{\pi}{2} $$So the total is $$ \frac{15}{4}\cdot \frac{\pi}{2}\left(\frac{(7/15)^2}{2}-\frac{(13/15)^2}{2}-\frac{(17/15)^2}{2}+\frac{(23/15)^2}{2}\right)= \frac{\pi}{2} $$
Of course, this method is merely computational: it has none of the theoretical trappings of the other answers. Nevertheless, the method could easily be generalized to other products and is quite amenable to reproduction on software. For instance, Mathematica thought long and hard about the original integral, but was able to perform the $128$ sub-integrals much quicker.
As a final remark, I'm reminded of the idea of quadrature from numerical analysis. It seems these integral-sinc-products, and more importantly their corresponding series, are 'exact' for the first seven until an error is introduced in the eighth. Perhaps I will investigate this further.
