(Argument transcribed from my work on AoPS here). The original problem of that thread, from the 2006 Colombian olympiad, was exactly the "original version" you referenced here.
We will start all sums at zero; the rows and columns of our matrices are indexed by $\{0,1,\cdots,n-1\}$.
Using the inner product notation $\langle f,g\rangle = \int_0^1 fg$, we have that the $ik$ entry $\frac1{i+k+1}$ of the Hilbert matrix $H$ is $\langle x^{i},x^{k}\rangle.$ If $f(x)=a_{0}+a_{1}x+\cdots+a_{n-1}x^{n-1}$ and $g(x)=b_{0}+b_{1}x+\cdots+b_{n-1}x^{n-1},$ this gives $\langle f,g\rangle=\sum_{i,k}a_{i}b_{k}h_{ik}=a^{T}Hb,$ where $a$ and $b$ are the column vectors of coefficients. Suppose $q$ is a polynomial of degree at most $n$ satisfying $\langle q,x^{k}\rangle=1$ for each $k,$ and let $a=(a_{0},a_{1},\dots,a_{n-1})^{T}$ be its coefficient vector in the standard basis. Our inner product relation becomes
$(0,\dots,0,1,0,\dots,0)\cdot H\cdot (a_{0},\dots,a_{k-1},a_{k},a_{k+1},\dots,a_{n-1})= 1$, or equivalently $I\cdot H\cdot\begin{pmatrix}a_{0}\\ a_{1}\\ \vdots\\ a_{n-1}\end{pmatrix}=\begin{pmatrix}1\\ 1\\ \vdots\\ 1\end{pmatrix}=j.$ (In all following computations, $j$ is this vector.)
Multiplying by $H^{-1}$, $a=H^{-1}j,$ and $\langle q,q\rangle=(H^{-1}j)^{T}H(H^{-1}j)=j^{T}H^{-1}j.$ For any matrix $M,$ $j^{T}Mj$ is the sum of the entries of $M$, and we have that the two versions of this problem are equivalent.
Now, we need to compute some details of $H^{-1}.$ To do this, we work with an orthonormal set of polynomials $p_i$, which are scaled Legendre polynomials. For each $k,$ $p_{k}(x)=\frac{\sqrt{2k+1}}{k!}\frac{d^{k}}{dx^{k}}(x-x^{2})^{k}.$ Orthogonality can be proven by integrating by parts repeatedly, differentiating the polynomial of lower degree and integrating the polynomial of higher degree; the boundary terms vanish at each step. The normalization constant also comes from an integration by parts: $\int_{0}^{1}\left(\frac{d^{k}}{dx^{k}}(x-x^{2})^{k}\right)^{2}\,dx=(2k)!\int_{0}^{1}x^{j}(1-x)^{j}\,dx,$ and this latter integral can be computed by integrating by parts another $j$ times.
Let $P$ be the upper triangular matrix with $ik$ entry the $x^{i}$ coefficient of $p_{k}.$ We have $P^{T}HP=\begin{pmatrix}\langle p_{0},p_{0}\rangle&\langle p_{0},p_{1}\rangle&\cdots&\langle p_{0},p_{n-1}\rangle\\ \langle p_{1},p_{0}\rangle&\langle p_{1},p_{1}\rangle&\cdots&\langle p_{1},p_{n-1}\rangle\\ \vdots&\vdots&\ddots&\vdots\\ \langle p_{n-1},p_{0}\rangle&\langle p_{n-1},p_{1}\rangle&\cdots&\langle p_{n-1},p_{n-1}\rangle\end{pmatrix}=I.$
Therefore, $H=(P^{T})^{-1}IP^{-1}=(PP^{T})^{-1},$ and $H^{-1}=PP^{T}.$ We get $j^{T}H^{-1}j=(j^{T}P)(P^{T}j)=\|P^{T}j\|^{2}.$ The $k$th entry of $P^{T}j$ is the sum of the entries in the $k$th column of $P,$ which is $p_{k}(1).$ By the symmetry of our definition, $p_{k}(1)=(-1)^{k}p_{k}(0)=(-1)^{k}\sqrt{2k+1}.$ Therefore, $j^{T}H^{-1}j=\sum_{k=0}^{n-1}2k+1=n^{2}$ and we have solved the problem.
That's not all, of course: we can get plenty of information on $H^{-1}$ from $P$. The upper left entry of $H^{-1}$ is the sum of the squares of the entries in the first row of $P;$ the $k$th entry in the first row is $p_{k}(0)=\sqrt{2k+1},$ and the upper left entry is also equal to $n^{2}.$
We can show that the entries of $H^{-1}$ are all integers: since $\frac1{k!}\frac{d^{k}}{dx^{k}}$ maps polynomials with integer coefficients to other polynomials with integer coefficients, $P=\begin{pmatrix}q_{0}&q_{1}&\cdots&q_{n-1}\end{pmatrix}\cdot \begin{pmatrix}1&0&\cdots&0\\ 0&\sqrt{3}&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&\sqrt{2n-1}\end{pmatrix},$ where each $q_{k}$ is an integer vector (which we will also interpret as the integer polynomial $\frac1{\sqrt{2k+1}}p_{k}$). This makes $PP^{T}=\begin{pmatrix}q_{0}&q_{1}&\cdots&q_{n-1}\end{pmatrix}\cdot \begin{pmatrix}1&0&\cdots&0\\ 0&3&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&2n-1\end{pmatrix}\cdot\begin{pmatrix}q_{0}^{T}\\ q_{1}^{T}\\ \cdots\\ q_{n-1}^{T}\end{pmatrix},$ a product of integer matrices.
We can go a step further and write the exact entries of $H^{-1}$ as sums; $(x-x^{2})^{k}=\sum_{i=0}^{k}(-1)^{i}\binom{k}{i}x^{i+k},$ and $q_{k}(x)=\frac1{k!}\sum_{i=0}^{k}(-1)^{i}\binom{k}{i}\frac{(k+i)!}{i!}x^{i}= \sum_{i=0}^{k}(-1)^{i}\binom{k}{i}\cdot\binom{k+i}{i}x^{i}.$ Let $(q_{k})_{i}$ be the coefficient of $x^{i}$ in $q_{k}.$
The $ik$ entry of $H^{-1}$ is then equal to $\sum_{l=0}^{n-1}(q_\ell)_{i}\cdot (2\ell+1)\cdot (q_\ell)_{k}=(-1)^{i+k}\sum_{\ell=0}^{n-1}(2\ell+1)\binom{\ell}{i}\cdot\binom{\ell}{k}\cdot\binom{\ell+i}{i}\cdot\binom{\ell+k}{k}.$ For many entries, most terms in this sum are zero. In an extreme case, the lower right entry of $H^{-1}$ has just one nonzero term in the sum; that entry is $(2n-1)\binom{2n-2}{n-1}^{2}.$
Since $P$ is triangular, the determinant of $H^{-1}$ is relatively easy to calculate; it's $\prod_{i=0}^{n-1}(2i+1)(q_{i})_{i}^{2}=\prod_{i=0}^{n-1}(2i+1)\binom{2i}{i}^{2}.$
If we write $\binom{2i}{i}=2^{2i}\cdot\frac{1\cdot 3\cdot 5\cdots (2i-1)}{2\cdot 4\cdot 6\cdots (2i)},$ the determinant becomes $\det(H^{-1})$
$=2^{4+8+\cdots+(4n-4)}\cdot 1^{2n-1}\cdot 2^{-(2n-2)}\cdot 3^{2n-3}\cdot 4^{-(2n-4)}\cdots (2n-1)^{1}\cdot (2n)^{0}$
$=2^{2n^{2}-2n}\prod_{k=1}^{2n}k^{(-1)^{k+1}(2n-k)}.$
