As Mathematica does, let's define the complete elliptic integral of the first kind in terms of the parameter $m=z^2$ so that we can write the integrand as $x^{n} K(x)K(1-x)$.
For nonnegative integers $n$, the result can also be expressed as $$\int_{0}^{1} x^{n} K(x) K(1-x) \, \mathrm dx = \frac{\pi}{2} \frac{(-1)^{n}}{n!} \lim_{\alpha \to 0} \frac{\mathrm d^{n}}{\mathrm d \alpha^{n}} \left( K(- \alpha)\right)^{2}. $$
I will use contour integration to show this.
The principal branch of $K(z)$ has a branch cut on $[1, \infty)$.
On the lower side of the branch cut, $$K(z) = \frac{1}{\sqrt{z}} \, K \left(\frac{1}{z}\right) - i K(1-z), \quad z >1.$$ While on the upper side of the branch cut, $$K(z) = \frac{1}{\sqrt{z}} \, K \left(\frac{1}{z}\right) + i K(1-z), \quad z >1.$$
(See this answer for an explanation.)
At $z=1$, $K(z)$ has the series expansion $$K(z) = - \frac{1}{2} \, \ln \left(1-z \right) +O(1).$$
And as $|z| \to \infty$, $K(z) $ has the asymptotic form $$K(z) \sim \frac{\ln(-z)}{2\sqrt{-z}}. $$
So by integrating the function $$f(z) = \frac{\left(K(z)\right)^{2}}{\alpha+z}, \quad \alpha >-1,$$ around a keyhole contour that is deformed around the branch cut, we get
$$\begin{align} \oint f(z) \, \mathrm dz &= \int_{1}^{\infty} \frac{\left(\frac{1}{\sqrt{x}} \, K \left(\frac{1}{x}\right) + i K(1-x) \right)^{2}}{\alpha+x} \, \mathrm dx + \int_{\infty}^{1} \frac{\left(\frac{1}{\sqrt{x}} \, K \left(\frac{1}{x}\right) - i K(1-x) \right)^{2}}{\alpha +x} \, \mathrm dx \\ &= 4i \int_{1}^{
\infty} \frac{K \left(\frac{1}{x} \right)K(1-x)}{\sqrt{x} (\alpha+x)} \, \mathrm dx \\ &= 2 \pi i \operatorname*{Res}_{z=-\alpha} \frac{\left(K(z)\right)^{2}}{\alpha +z} \\ &= 2 \pi i \, \left(K(- \alpha)\right)^{2}. \end{align}$$
Therefore, $$ \int_{1}^{\infty} \frac{K \left(\frac{1}{x} \right)K(1-x)}{\sqrt{x} (\alpha+x)} \, \mathrm dx = \frac{\pi}{2} \, \left(K(-\alpha)\right)^{2}. $$
And by making the substitution $ u = \frac{1}{x}$, we have $$\int_{0}^{1} \frac{K(u) K \left(1- \frac{1}{u} \right)}{(1+\alpha u) \sqrt{u}} \, \mathrm du \overset{\spadesuit}{=} \int_{0}^{1} \frac{K(u) K(1-u)}{1+\alpha u} \, \mathrm du = \frac{\pi}{2} \left(K(-\alpha)\right)^{2}. $$
If we differentiate both sides of the above equation $n$ times and then let $\alpha \to 0$, we get $$(-1)^n \, n! \int_{0}^{1} u^{n} K(u) K(1-u) \, \mathrm du = \frac{\pi}{2} \lim_{\alpha \to 0} \frac{\mathrm d^{n}}{\mathrm d \alpha^{n}} \left(K(-\alpha) \right)^{2}. $$
For $n=0$, we have $$\int_{0}^{1} K(x) K(1-x) \, \mathrm dx = \frac{\pi}{2} \left(K(0) \right)^{2} = \frac{\pi}{2} \left(\frac{\pi}{2} \right)^{2} = \frac{\pi^{3}}{8}. $$
For $n=1$, we have $$\int_{0}^{1} x K(x) K(1-x) \, \mathrm dx = - \frac{\pi}{2} \lim_{\alpha \to 0} 2 K(-\alpha) \frac{\mathrm d}{\mathrm d \alpha} K(-\alpha) = -\pi \left(\frac{\pi}{2} \right) \left(-\frac{\pi}{8} \right) = \frac{\pi^{3}}{16}.$$
For $n=2$, we have $$ \begin{align} \int_{0}^{1} x^{2} K(x) K(1-x) \, \mathrm dx &= \frac{\pi}{2} \frac{1}{2!} \lim_{\alpha \to 0} 2\left(\left(\frac{\mathrm d}{\mathrm d \alpha }K(-\alpha) \right)^{2} +K(-\alpha) \frac{\mathrm d^{2}}{\mathrm d \alpha^{2}} K(-\alpha) \right) \\ &= \frac{\pi}{2} \left( \left(-\frac{\pi}{8} \right)^{2}+\left(\frac{\pi}{2} \right) \left(\frac{9 \pi}{64} \right) \right) \\ &= \frac{11 \pi^{3}}{256}. \end{align}$$
And so on.
The derivatives of $K(-\alpha)$ at $\alpha=0$ can be extracted from coefficients of the Maclaurin series $$K(-\alpha) = \frac{\pi}{2} \, _1F_{2} \left(\frac{1}{2}, \frac{1}{2}; 1; -\alpha \right) = \frac{\pi}{2} \sum_{n=0}^{\infty} \left(\frac{(2n)!}{2^{2n} (n!)^{2}} \right)^{2} (-\alpha)^{n}.$$
$\spadesuit$ For $u >0$, we have
$$ \begin{align} \frac{1}{\sqrt{u}} \, K \left(1- \tfrac{1}{u} \right) &= \frac{1}{\sqrt{u}} \int_{0}^{\pi/2} \frac{d\theta}{\sqrt{1-\left(1-\tfrac{1}{u}\right)\sin^{2}(\theta})} \\ &= \int_{0}^{\pi/2} \frac{\mathrm d \theta}{\sqrt{u-(u-1)\sin^{2}(\theta)}} \\ &= - \int_{\pi/2}^{0} \frac{\mathrm d \phi}{\sqrt{u-(u-1) \sin^{2}(\pi/2- \phi)}} \\ &= \int_{0}^{\pi/2} \frac{\mathrm d \phi}{\sqrt{u-(u-1) \cos^{2}(\phi)}} \\ &= \int_{0}^{\pi/2} \frac{\mathrm d \phi}{\sqrt{u-(u-1) (1-\sin^{2}(\phi))}} \\ &= \int_{0}^{\pi/2} \frac{\mathrm d \phi }{\sqrt{1-(1-u) \sin^{2}(\phi)}} \\ &=K(1-u). \end{align}$$
