<s>[INST] Combinatorial identity: $\sum_{i,j \ge 0} \binom{i+j}{i}^2 \binom{(a-i)+(b-j)}{a-i}^2=\frac{1}{2} \binom{(2a+1)+(2b+1)}{2a+1}$ In my research, I found this identity and as I experienced, it's surely right. But I can't give a proof for it.
Could someone help me?
This is the identity:
let $a$ and $b$ be two positive integers; then:
$\sum_{i,j \ge 0} \binom{i+j}{i}^2 \binom{(a-i)+(b-j)}{a-i}^2=\frac{1}{2} \binom{(2a+1)+(2b+1)}{2a+1}$.
 [/INST] Let us denote
$$S=\sum_{i,j \ge 0} \binom{i+j}{i}^2 \binom{(a-i)+(b-j)}{a-i}^2.$$
First, let $s=i+j$ so that
$$S = \sum_{s\geq 0}\sum_{i=0}^s \binom{s}{i}^2 \binom{a+b-s}{a-i}^2.$$
Consider the generating function
$$F(x,y) = \sum_{s,i} \binom{s}{i}^2 x^i y^s = (1-2y+y^2-2xy-2xy^2+x^2y^2)^{-1/2}.$$
Then $S$ is nothing else but the coefficient of $x^a y^{a+b}$ in
$$F(x,y)^2 = (1-2y+y^2-2xy-2xy^2+x^2y^2)^{-1}$$
$$ = \frac{1}{4y\sqrt{x}}\left(\frac{1}{1-y(1+x+2\sqrt{x})} - \frac{1}{1-y(1+x-2\sqrt{x})}\right)$$
$$ = \frac{1}{4y\sqrt{x}}\left(\frac{1}{1-y(1+\sqrt{x})^2} - \frac{1}{1-y(1-\sqrt{x})^2}\right).$$
(derivation simplified)
The coefficient of $y^{a+b}$ equals
$$[y^{a+b}]\ F(x,y)^2
=\frac{1}{4} \frac{(1+\sqrt{x})^{2(a+b+1)} - (1-\sqrt{x})^{2(a+b+1)}}{\sqrt{x}}.$$
Now we trivially conclude that
$$S = [x^ay^{a+b}]\ F(x,y)^2 = \frac{1}{2}\binom{2(a+b+1)}{2a+1}.$$

UPDATE. Alternatively to computing the coefficient of $x^ay^{a+b}$, one can follow the venue of Fedor Petrov's proof. This way one needs to consider the generating function
$$G(x,y) = \sum_{m,n}\binom{m}{n} x^ny^m = \frac{1}{1-y-xy}$$
and verify that
$$8xy^2F(x^2,y^2)^2 = G(x,y) + G(x,-y) - G(-x,y) - G(-x,-y).$$
</s>