I believe this can also be solved using double integrals.
It is possible (if I remember correctly) to justify switching the order of integration to give the equality:
$$\int_{0}^{\infty} \Bigg(\int_{0}^{\infty} e^{-xy} \sin x \,dy \Bigg)\, dx = \int_{0}^{\infty} \Bigg(\int_{0}^{\infty} e^{-xy} \sin x \,dx \Bigg)\,dy$$
Notice that
$$\int_{0}^{\infty} e^{-xy} \sin x\,dy = \frac{\sin x}{x}$$
This leads us to
$$\int_{0}^{\infty} \Big(\frac{\sin x}{x} \Big) \,dx = \int_{0}^{\infty} \Bigg(\int_{0}^{\infty} e^{-xy} \sin x \,dx \Bigg)\,dy$$
Now the right hand side can be found easily, using integration by parts.
$$\begin{align*}
I &= \int e^{-xy} \sin x \,dx = -e^{-xy}{\cos x} - y \int e^{-xy} \cos x \, dx\\
&= -e^{-xy}{\cos x} - y \Big(e^{-xy}\sin x + y \int e^{-xy} \sin x \,dx \Big)\\
&= \frac{-ye^{-xy}\sin x - e^{-xy}\cos x}{1+y^2}.
\end{align*}$$
Thus $$\int_{0}^{\infty} e^{-xy} \sin x \,dx = \frac{1}{1+y^2}$$
Thus $$\int_{0}^{\infty} \Big(\frac{\sin x}{x} \Big) \,dx = \int_{0}^{\infty}\frac{1}{1+y^2}\,dy = \frac{\pi}{2}.$$
It turns out that this integral takes on a very simple form amenable to analysis via residues. Let $u = \sin{x}/\sin{(x+\pi/3)}$. We may then find that (+)
$$\tan{x} = \frac{(\sqrt{3}/2)u}{1-(u/2)}$$
A little bit of algebra reveals a very nice form for the differential:
$$dx = \frac{\sqrt{3}}{2} \frac{du}{1-u+u^2}$$
so the original integral takes on a much simpler-looking form:
$$\frac{\sqrt{3}}{2} \int_0^1 du \frac{\log^2{u}}{1-u+u^2}$$
This is not ready for contour integration yet. We may transform this into such an integral by substituting $u=1/v$ and observing that
$$\int_0^1 du \frac{\log^2{u}}{1-u+u^2} = \int_1^\infty du \frac{\log^2{u}}{1-u+u^2} = \frac{1}{2} \int_0^\infty du \frac{\log^2{u}}{1-u+u^2}$$
We may now analyze that last integral via the residue theorem. Consider the integral
$$\oint_C dz \frac{\log^3{z}}{1-z+z^2}$$
where $C$ is a keyhole contour that passes up and back along the positive real axis. It may be shown that the integral along the large and small circular arcs vanish as the radii of the arcs goes to $\infty$ and $0$, respectively. We may then write the integral in terms of positive contributions just above the real axis and negative contributions just below. The result is
$$\oint_C dz \frac{\log^3{z}}{1-z+z^2} = \begin{array}\\ i \left ( - 6 \pi \int_0^\infty du \frac{\log^2{u}}{1-u+u^2} + 8 \pi^3 \int_0^\infty du \frac{1}{1-u+u^2} \right ) \\ + 12 \pi^2 \int_0^\infty du \frac{\log{u}}{1-u+u^2} \end{array}$$
We set this equal to $i 2 \pi$ times the sum of the residues of the poles of the integrand within $C$. The poles are $z \in \{e^{i \pi/3},e^{i 5\pi/3}\}$. The residues are
$$\mathrm{Res}_{z=e^{i \pi/3}} = -\frac{\pi^3}{27 \sqrt{3}}$$
$$\mathrm{Res}_{z=e^{i 5\pi/3}} = \frac{125 \pi^3}{27 \sqrt{3}}$$
$i 2 \pi$ times the sum of these residues is then
$$i \frac{248 \pi^4}{27 \sqrt{3}}$$
Equating imaginary parts of the integral to the above quantity, we see that
$$ - 6 \pi \int_0^\infty du \frac{\log^2{u}}{1-u+u^2} + 8 \pi^3 \int_0^\infty du \frac{1}{1-u+u^2} = \frac{248 \pi^4}{27 \sqrt{3}}$$
Now, I will state without proof for now that (++)
$$\int_0^\infty du \frac{1}{1-u+u^2} = \frac{4 \pi}{3 \sqrt{3}}$$
Then with a little arithmetic, we find that
$$\int_0^\infty du \frac{\log^2{u}}{1-u+u^2} = \frac{20 \pi^3}{81 \sqrt{3}}$$
The integral we want is $\sqrt{3}/4$ times this value; therefore
$$\int_0^{\pi/3} dx \log^2{\left [ \frac{\sin{x}}{\sin{(x + \pi/3)}}
\right ]} = \frac{5\pi^3}{81}$$
Proof of (++)
Now, to prove (++), we go right back to the observation (+) that
$$ x = \int \frac{du}{1-u+u^2} \implies \tan{\left ( \frac{\sqrt{3}}{2} x \right )} = \frac{(\sqrt{3}/2)u}{1-(u/2)}$$
Therefore
$$\int_0^1 \frac{du}{1-u+u^2} = \frac{2}{\sqrt{3}} \left [ \arctan \left ( \frac{(\sqrt{3}/2) u}{1-(u/2)} \right ) \right ]_0^1 = \frac{2}{\sqrt{3}} \frac{\pi}{3}$$
and we showed that this is $1/2$ the integral over $[0,\infty)$, and
$$\int_0^\infty \frac{du}{1-u+u^2} = \frac{4 \pi}{3 \sqrt{3}}$$
QED
Best Answer
Here's another approach.
We have $$\begin{eqnarray*} \int_0^\infty dx\, \left(\frac{\sin x}{x}\right)^n &=& \lim_{\epsilon\to 0^+} \frac{1}{2} \int_{-\infty}^\infty dx\, \left(\frac{\sin x}{x-i\epsilon}\right)^n \\ &=& \lim_{\epsilon\to 0^+} \frac{1}{2} \int_{-\infty}^\infty dx\, \frac{1}{(x-i\epsilon)^n} \left(\frac{e^{i x}-e^{-i x}}{2i}\right)^n \\ &=& \lim_{\epsilon\to 0^+} \frac{1}{2} \frac{1}{(2i)^n} \int_{-\infty}^\infty dx\, \frac{1}{(x-i\epsilon)^n} \sum_{k=0}^n (-1)^k {n \choose k} e^{i x(n-2k)} \\ &=& \lim_{\epsilon\to 0^+} \frac{1}{2} \frac{1}{(2i)^n} \sum_{k=0}^n (-1)^k {n \choose k} \int_{-\infty}^\infty dx\, \frac{e^{i x(n-2k)}}{(x-i\epsilon)^n}. \end{eqnarray*}$$ If $n-2k \ge 0$ we close the contour in the upper half-plane and pick up the residue at $x=i\epsilon$. Otherwise we close the contour in the lower half-plane and pick up no residues. The upper limit of the sum is thus $\lfloor n/2\rfloor$. Therefore, using the Cauchy differentiation formula, we find $$\begin{eqnarray*} \int_0^\infty dx\, \left(\frac{\sin x}{x}\right)^n &=& \frac{1}{2} \frac{1}{(2i)^n} \sum_{k=0}^{\lfloor n/2\rfloor} (-1)^k {n \choose k} \frac{2\pi i}{(n-1)!} \left.\frac{d^{n-1}}{d x^{n-1}} e^{i x(n-2k)}\right|_{x=0} \\ &=& \frac{1}{2} \frac{1}{(2i)^n} \sum_{k=0}^{\lfloor n/2\rfloor} (-1)^k {n \choose k} \frac{2\pi i}{(n-1)!} (i(n-2k))^{n-1} \\ &=& \frac{\pi}{2^n (n-1)!} \sum_{k=0}^{\lfloor n/2\rfloor} (-1)^k {n \choose k} (n-2k)^{n-1}. \end{eqnarray*}$$ The sum can be written in terms of the hypergeometric function but the result is not particularly enlightening.