An easier evaluation of $\det\limits_{1\leqslant i,j\leqslant n}\left\{\frac{x_i-x_j}{x_i+x_j}\right\}$

One can prove the identity with induction on the dimension $n$, using elementary transformations on the determinant.

The identity holds for $n=0$, with the empty product and the empty determinant both being equal to $1$.

For the induction step $n \to n+2$ we compute $$ R_{n+2}(x_1, \ldots, x_n, a, b) = \begin{vmatrix} q(x_1, x_1) & \cdots & q(x_1, x_n) & q(x_1, a) & q(x_1, b) \\ \vdots & \ddots & \vdots &\vdots &\vdots \\ q(x_n, x_1) & \cdots & q(x_n, x_n) & q(x_n, a) & q(x_n, b) \\ q(a, x_1) & \cdots & q(a, x_n) & 0 & q(a, b) \\ q(b, x_1) & \cdots & q(b, x_n) & q(b, a) & 0 \end{vmatrix} $$ where I have used $$ q(x, y) = \frac{x-y}{x+y} $$ to shorten the expressions a bit. We add multiples of the last two rows to the first $n$ rows to eliminate the last two entries in those rows. This gives $$ R_{n+2}(x_1, \ldots, x_n, a, b) = \begin{vmatrix} d_{1,1} & \cdots & d_{1,n} & 0 & 0 \\ \vdots & \ddots & \vdots &\vdots &\vdots \\ d_{n,1} & \cdots & d_{n,n} & 0 & 0 \\ * & \cdots & * & 0 & q(a, b) \\ * & \cdots & * & q(b, a) & 0 \end{vmatrix} $$ with $$ d_{i,j} = q(x_i, x_j) - \frac{q(x_i, a)q(b,x_j)}{q(b, a)} - \frac{q(x_i, b)q(a,x_j)}{q(a, b)} \, . $$ Now the magic (?) happens: This expression simplifies to $$ d_{i,j} = q(x_i, a)q(x_i, b)q(x_j, a)q(x_j, b)q(x_i, x_j) \, . $$

The common factors $q(x_i, a)$ and $q(x_i, b)$ can be extracted from each row and each column in $ \det \{ d_{i,j} \}$, and it follows that $$ R_{n+2}(x_1, \ldots, x_n, a, b) = \det \{ d_{i,j} \} q(a, b)^2 \\ = \det \{ q(x_i, x_j) \}\left(\prod_{i=1}^n q(x_i, a)^2 q(x_i, b)^2\right) q(a, b)^2 \\ = R_{n}(x_1, \ldots, x_n)\left(\prod_{i=1}^n q(x_i, a)^2 q(x_i, b)^2\right) q(a, b)^2 \, , $$ which is exactly the induction step.

To evaluate $D_n(\bar{x},\bar{y},\bar{z})$, we write $$\frac{x_i+z_j}{x_i+y_j}=(z_j-y_j)\left(\frac{1}{x_i+y_j}-\frac{1}{y_j-z_j}\right),$$ so that $$D_n(\bar{x},\bar{y},\bar{z})=\prod_j(z_j-y_j) \begin{vmatrix} 1&\frac{1}{y_1-z_1}&\ldots&\frac{1}{y_n-z_n} \\1&\frac{1}{x_1+y_1}&\ldots&\frac{1}{x_1+y_n} \\\cdot&\cdots&\ddots&\cdots \\1&\frac{1}{x_n+y_1}&\ldots&\frac{1}{x_n+y_n} \end{vmatrix}.$$

Expanding the determinant along the first row, we get $$D_n(\bar{x},\bar{y},\bar{z})=\prod_j(z_j-y_j)\left(C_n(\bar{x},\bar{y})+\sum_{k=1}^{n}\frac{C_n^{(k)}(\bar{x},\bar{y})}{y_k-z_k}\right),$$ where $C_n$ is the well-known Cauchy determinant: $$C_n(\bar{x},\bar{y}):=\det_{1\leqslant i,j\leqslant n}\left\{\frac{1}{x_i+y_j}\right\}=\frac{\prod_{i<j}(x_i-x_j)(y_i-y_j)}{\prod_{i,j}(x_i+y_j)},$$ and $C_n^{(k)}$ is obtained by replacing the $k$-th column of $C_n$ by a column of $\color{blue}{1}$s: $$C_n^{(k)}(\bar{x},\bar{y})=\lim_{y_k\to\infty}y_k C_n(\bar{x},\bar{y})=C_n(\bar{x},\bar{y})\frac{\prod_i(x_i+y_k)}{\prod_{i\neq k}(y_k-y_i)}.$$

Thus, we obtain Duparc's result $$D_n(\bar{x},\bar{y},\bar{z})=C_n(\bar{x},\bar{y})\prod_j(z_j-y_j)\left(1-\sum_i\frac{y_i+x_i}{y_i-z_i}\prod_{j\neq i}\frac{y_i+x_j}{y_i-y_j}\right).$$

Since $R_n(\bar{x})=D_n(\bar{x},\bar{x},-\bar{x})$, we're left to show $$\sum_i\prod_{j\neq i}\frac{x_i+x_j}{x_i-x_j}=\begin{cases}0,&n\text{ is even}\\1,&n\text{ is odd}\end{cases}.\tag{*}\label{essential}$$

For this, we do partial fraction expansion of $$F(x):=\prod_j\frac{x+x_j}{x-x_j}=A_0+\sum_i\frac{A_i}{x-x_i},$$ with $A_0=1$ and $A_i=\lim\limits_{x\to x_i}(x-x_i)F(x)=\ldots$ resulting in $$\prod_j\frac{x+x_j}{x-x_j}=1+\sum_i\frac{2x_i}{x-x_i}\prod_{j\neq i}\frac{x_i+x_j}{x_i-x_j}.$$ To obtain \eqref{essential}, it just remains to put $x=0$.

This question was linked to that question which brings up the fact that you can prove the identity above by squaring the identity $$\mathrm{Pf}\left(\frac{z_j - z_i}{z_j + z_i}\right)_{i,j} = \prod_{i<j} \frac{z_j - z_i}{z_j + z_i}$$ using the standard fact that the square of the Pfaffian is the determinant.

This identity in turn comes from the fact the LHS can be put over the common denominator on the RHS with a numerator which is of the same degree as that of the RHS. Then one reasons the LHS is antisymmetric in $z$ and hence the numerator of the RHS divides that of the LHS. Then checking one coefficient gives the identity.