Prove a Vandermonde Determinant Formula I seem to have stumbled upon a strange formula, and I'm not quite sure that I understand it. Consider the determinant standard Vandermonde matrix
$${}^nD = \left|\begin{matrix}
    1 & x_{0} & x^{2}_{0} & \cdots & x^{n}_{0} \\
    1 & x_{1} & x^{2}_{1} & \cdots & x^{n}_{1} \\
\vdots & \vdots & \vdots & \ddots &  \vdots \\
    1 & x_{n} & x^{2}_{n} & \cdots & x^{n}_{n} \\
    \end{matrix}\right|$$
In numerical analysis, we can often use this matrix (along with smaller, similar matrices) to generate the coefficients of the $n$-th degree interpolating polynomial for a general function $f$ passing through the set of points $\{(x_i,f(x_i))\}_{i = 0}^n$. Indeed, if we define the notation
$${}^nD_i(f) = \left|\begin{matrix}
    1 & x_{0} & x^{2}_{0} & \cdots & x^{i-1}_0 & f(x_0) & x^{i + 1}_0 & \cdots & x^{n}_{0} \\
    1 & x_{1} & x^{2}_{1} & \cdots & x^{i-1}_1 & f(x_1) & x^{i + 1}_1 & \cdots & x^{n}_{1} \\
    \vdots & \vdots & \vdots & \ddots &  \vdots & \vdots & \vdots & \ddots &\vdots \\
    1 & x_{n} & x^{2}_{n} & \cdots & x^{i-1}_n & f(x_n) & x^{i + 1}_n & \cdots & x^{n}_{n} \\
    \end{matrix}\right|$$
then the $n$-th order Lagrange approximation for the function $f$ can be obtained from the formula
$$p(x) = f(x_0) + \frac{\det\left({}^1D_1(f)\right)}{\det(D)}(x - x_0) + \cdots +  \frac{\det\left({}^nD_n(f)\right)}{\det(D)}(x - x_0)\cdots (x-x_{n-1}).$$
Now we can define an analogous notation wherein we replace not just the column corresponding to $x^i$ but also some other $x^k$ by
$${}^nD_{i,k}(f,g) = \left|\begin{matrix}
    1 & x_{0} & x^{2}_{0} & \cdots & x^{i-1}_0 & f(x_0) & x^{i + 1}_0 & \cdots & x^{k-1}_0 & g(x_0) & x^{k + 1}_0 & \cdots & x^{n}_{0} \\
    1 & x_{1} & x^{2}_{1} & \cdots & x^{i-1}_1 & f(x_1) & x^{i + 1}_1 & \cdots & x^{k-1}_0 & g(x_1) & x^{k + 1}_0 & \cdots & x^{n}_{1} \\
    \vdots & \vdots & \vdots & \ddots &  \vdots & \vdots & \vdots & \ddots &\vdots & \vdots & \vdots & \ddots &\vdots\\
    1 & x_{n} & x^{2}_{n} & \cdots & x^{i-1}_n & f(x_n) & x^{i + 1}_n & \cdots & x^{k-1}_0 & g(x_n) & x^{k + 1}_0 & \cdots & x^{n}_{n} \\
    \end{matrix}\right|$$
and it would seem that the following identity holds:
$${}^nD_{i,k}(f,g)\:{}^nD = {}^nD_i(f)\:{}^nD_k(g) - {}^nD_i(g)\:{}^nD_k(f).$$
The first few cases are fairly easy to verify (some naive Gaussian elimination / symmetric polynomial tricks can handle up to $n = 5$ nicely, and I verified up to $n = 12$ with a computer), but I'm having trouble coming up with a general proof of this fact. If anyone has an idea as to how to show this, I would appreciate the insight. Thank you!
 A: Hint
Say that the polynomial functions are given below:
$$
f(x) = \sum_{j=0}^{n} f_{j}x^{j}
$$
Let $M$ be an identity matrix of size $n+1$ but with its $i$-th column replaced with
$
\begin{bmatrix}
f_{0} & f_{1} & \dots & f_{n}
\end{bmatrix}^{T}
$. Then we have the following relation:
$$
\begin{align}
D_{i}(f)=DM
\implies
\det{\left(D_{i}(f)\right)} &=
\det{\left(D\right)}
\cdot 
\det{\left(M\right)}
\\
&=
\det{\left(D\right)}
\cdot
f_{i}
\end{align}
$$
--
Say that the second polynomial functions are given below:
$$
g(x) = \sum_{j=0}^{n} g_{j}x^{j}
$$
Now let $M$ be an identity matrix of size $n+1$ but with its $i$-th and $k$-th columns replaced with
$
\begin{bmatrix}
f_{0} & f_{1} & \dots & f_{n}
\end{bmatrix}^{T}
$ and $
\begin{bmatrix}
g_{0} & g_{1} & \dots & g_{n}
\end{bmatrix}^{T}
$ respectively. Then we have the following relation:
$$
\begin{align}
D_{i,k}(f,g)=DM
\implies
\det{\left(D_{i,k}(f,g)\right)} &=
\det{\left(D\right)}
\cdot
\det{\left(M\right)}
\\
&=
\det{\left(D\right)}
\cdot
\det{\left(\begin{bmatrix}f_{i} & f_{k} \\ g_{i} & g_{k}\end{bmatrix}\right)}
\end{align}
$$
--
Hope this can help you see why you have such identity. Next, what about $D_{i,j,k,...}(f,g,h,...)$?
