The left inequality.
We need to prove that
$$\sum_{i=1}^n\left(\frac{(n - y)^2}{y^2(x_i^2 + x_i + 1) - 3ny(x_i + 1) + 3n^2}-\frac{1}{3}\right)\geq0$$ or
$$\sum_{i=1}^n\frac{(x_i-1)(3n-2y-yx_i)}{y^2(x_i^2 + x_i + 1) - 3ny(x_i + 1) + 3n^2}\geq0$$ or $$\sum_{i=1}^n\left(\frac{(x_i-1)(3n-2y-yx_i)}{y^2(x_i^2 + x_i + 1) - 3ny(x_i + 1) + 3n^2}-\frac{x_i-1}{n-y}\right)\geq0$$ or
$$\sum_{i=1}^n\frac{(x_i-1)^2(2n-y-yx_i)}{y^2(x_i^2 + x_i + 1) - 3ny(x_i + 1) + 3n^2}\geq0,$$ which is obvious.
For $n\geq3$ we can prove a right inequality by the same way.
Also, easy to check that the right inequality is true for $n=2$.
Define $f(x, y) = \sqrt{|x + y|} - \sqrt{|x-y|}$. Then we want to prove that
\begin{align*}
\mathbb{E}[f(X, Y)] \ge 0
\end{align*}
where $X, Y$ independent and uniformly distributed over the set of values $\{x_1, \cdots, x_n\}$.
[Remark 1: I know the IMO is motivated to use techniques only available at the secondary school level, but the probabilistic notation allows me to naturally represent a bunch of summation notations in a succinct form, and can be translated easily to the summation notation in the original problem.]
First, let's observe the region in the Cartesian plane where $f(x, y) \ge a$ and $f(x, y) \le -a$, respectively:
And therefore
\begin{align*}
\mathbb{P}(f(X, Y) \ge a) &= \mathbb{P}\left(X \ge \frac{a}{2}, Y \ge \frac{a}{2}\right) + \mathbb{P}\left(X \le -\frac{a}{2}, Y \le -\frac{a}{2}\right) \\
&= \mathbb{P}\left(X \ge \frac{a}{2}\right)^2 + \mathbb{P}\left(X \le -\frac{a}{2}\right)^2
\end{align*}
and
\begin{align*}
\mathbb{P}(f(X, Y) \le a) &= \mathbb{P}\left(X \ge \frac{a}{2}, Y \le -\frac{a}{2}\right) + \mathbb{P}\left(X \le -\frac{a}{2}, Y \ge \frac{a}{2}\right) \\
&= 2\mathbb{P}\left(X \ge \frac{a}{2}\right) \left(X \le -\frac{a}{2}\right)
\end{align*}
Then, we need the lemma that
\begin{align*}
\mathbb{E}[A] = \int_0^\infty [\mathbb{P}(A \ge a) - \mathbb{P}(A \le -a)]da
\end{align*}
which can be proved by noticing that $\mathbb{E}[A] = \mathbb{E}[A^+] - \mathbb{E}[A^-]$, where $A^+ = \max(A, 0)$ and $A^- = -\min(A, 0)$ and then applying the proof for the tail sum formula. Letting $A = f(X, Y)$,
\begin{align*}
\mathbb{E}[f(X, Y)] &= \int_0^\infty[\mathbb{P}(f(X, Y) \ge a) - \mathbb{P}(f(X, Y) \le -a)]da \\
&= \int_0^\infty\left[\mathbb{P}\left(X \ge \frac{a}{2}\right) - \mathbb{P}\left(X \le -\frac{a}{2}\right)\right]^2da \\
&\ge 0
\end{align*}
Remark 2: Note this proof works for any $f(x, y) = g(|x+y|) - g(|x - y|)$ for non-decreasing $g: \mathbb{R}_{\ge 0} \rightarrow \mathbb{R}$ and general distributions for $X, Y$. For example, we also have the result
\begin{align*}
\sum_{i,j} w_i w_j \log(1 + |x_i - x_j|) \le \sum_{i,j} w_i w_j \log(1 + |x_i + x_j|)
\end{align*}
for weights $w_1, \cdots, w_n \ge 0$.
Remark 3: Completely forgot to answer the original problem. Note that if we choose $g(x) = \sqrt{x}$ and define $X, Y$ to be random variables distributed over $x_1, \cdots, x_n$ with weights $w_1, \cdots, w_n \ge 0$, then we effectively get, for $\mathbf{M} = [M_{ij}]_{ij}$,
\begin{align*}
\mathbf{w}^\intercal \mathbf{M}\mathbf{w} \ge 0
\end{align*}
for $\mathbf{w} \in \mathbb{R}^n_{\ge 0}$; we need $\mathbf{w} \in \mathbb{R}^n$ in order to get genuine positive semidefiniteness. But it appears that if we replace $\mathbb{P}$ with quasiprobability distribution (e.g. negative measures are allowed), everything in the proof still seems to go through.
Best Answer
Edit. Eventhough I had in mind the version of the ACM I needed, it all messed up when I tried to find a reference. As for the valid version, I will refer to Vasile Cirtoaje's Algebraic Inequalities. Old and new Methods, p. $267$.
We will employ a powerful technique developed by Vasile Cirtoaje back in $2006$ called the Arithmetic Compensation Method (see, for instance,
this document).Let to this end $$F(x_1, \ldots, x_n):=\sum_{1\leqslant i<j\leqslant n}\frac{x_ix_j}{(1-x_i)(1-x_j)}$$ which is clearly symmetric and continuous on $S:=\left\{(x_1, \ldots,x_n)\mid \sum_{i=1}^n x_i=\frac12, \forall i: x_i\geqslant 0\right\}$. We will now refer to
the Remark 1.1 fromthe document linked above, which basically states thatNotice that (\ref{(i)}) is equivalent to (where we will denote, for brevity, $y:=\frac{1}2(x_1+x_2)$) \begin{align*} \frac{x_1x_2}{(1-x_1)(1-x_2)}-\frac{y^2}{(1-y)^2}+\left(\frac{x_1}{1-x_1}+\frac{x_2}{1-x_2}-\frac{2y}{1-y}\right)\sum_{i=3}^n \frac{x_i}{1-x_i}>0\\ \iff \frac{y^2-x_1x_2}{(1-x_1)(1-x_2)(1-y)^2}\left[2y-1+2(1-y)\sum_{i=3}^n \frac{x_i}{1-x_i}\right]>0\\ \iff 2y-1+2(1-y)\sum_{i=3}^n \frac{x_i}{1-x_i}>0 \end{align*}
Where the last equivalence follows from $y^2-x_1x_2\geqslant 0$. We will now turn to the second inequality: \begin{align*} F(x_1,x_2,x_3,\ldots,x_n)- F(0, x_1+x_2, x_3, \ldots, x_n)\leqslant0\\ \iff \frac{x_1x_2}{(1-x_1)(1-x_2)}+\left(\frac{x_1}{1-x_1}+\frac{x_2}{1-x_2}-\frac{2y}{1-y}\right)\sum_{i=3}^n \frac{x_i}{1-x_i}\leqslant0\\ \iff \frac{x_1x_2}{(1-x_1)(1-x_2)(1-2y)}\left[1-2y+2(y-1)\sum_{i=3}^n \frac{x_i}{1-x_i}\right]\leqslant 0\\ \end{align*}
Which is easily seen to follow from (\ref{(i)}). Thus, we conclude that \begin{align*} F(x_1,x_2,x_3,\ldots,x_n)&\leqslant \underset{1\leqslant k\leqslant n}\max F\left(\frac1{2k}, \ldots, \frac1{2k},0,\ldots ,0\right)\\ &= \underset{1\leqslant k\leqslant n}\max \binom{k}{2}\frac{1}{(2k-1)^2}\\ &= \underset{1\leqslant k\leqslant n}\max \frac{k(k-1)}{2(2k-1)^2} \end{align*}
But $f: x\mapsto \frac{x(x-1)}{2(2x-1)^2}$ is increasing on $[1,\infty)$, and, hence, the result follows.
Remark. Based on the difficulty of the other contest problems, this solution seems a bit overkill to me, but I have failed in the attempt to find a simpler method, since most well-known inequalities work in the other direction...