Gradient of a multi-dimensional multi-variable function (First of all, feel free to suggest a better title for the question, I might just be totally missing the naming, hence not finding my answer because of that :) )
I understand how to compute the partial derivative of some function $f(x, y)$, with respect its different variables, and how to get the gradient of the function from that.
Now, if I have a function $f$ that takes, let say, two 2D vectors $p1$ and $p2$ as inputs and I want to find the gradient of this function with respect to each point $\nabla_{p1}f(p1, p2)$ and $\nabla_{p2}f(p1, p2)$. This is where I'm totally lost. How is this computed ?
For example, if $f(p1,p2) = |p1 - p2| - d$, $|p1 - p2|$ being the distance between the two points (or the norm of the vector defined by those points) and $d$ being a constant, how does one compute $\nabla_{p1}f(p1, p2)$ and $\nabla_{p2}f(p1, p2)$ ?
In that case, the results I need to find are
$$\nabla_{p1}f(p1, p2) = \frac{p1-p2}{|p1 - p2|}$$
and
$$\nabla_{p2}f(p1, p2) = -\frac{p1-p2}{|p1 - p2|}$$
but I do not understand how to find this result.
Edit:
To add a bit more context, I want to understand how to compute those formulas to be able to put them in a computer graphics physics simulation loop (namely using Position-Based Dynamics).
The function $f$ is in fact a constraint between the inputs (in that case, we want the two points to keep a certain distance from each other).
 A: Expanding my comment. You can use the calculus of differentials for this kind of calculations. My favorite tutorial on this is Practical Guide to Matrix Calculus for Deep Learning by Andrew Delong. It is about the use of differentials for matrices, but of course it also can be used for vectors. Rule (17) in the paper is exactly about your question.
$$dy=\sum_kA^k\cdot dX^k \rightarrow \nabla_ky=A^k$$
so applying the rule to your example
$$df=d(\sqrt{(p_1-p_2)\cdot(p_1-p_2)}-d)=\text{by (15)}=\frac{1}{2|p_1-p_2|}d[(p_1-p_2)\cdot(p_1-p_2)]=$$
$$=\text{by (13)}=\frac{1}{2|p_1-p_2|}[d(p_1-p_2)\cdot(p_1-p_2)+(p_1-p_2)\cdot d(p_1-p_2)]=$$
$$=\text{by (11)}=\frac{1}{2|p_1-p_2|}[(dp_1-dp_2)\cdot(p_1-p_2)+(p_1-p_2)\cdot (dp_1-dp_2)]=$$
$$=\frac{p_1-p_2}{|p_1-p_2|}\cdot dp_1+\frac{-(p_1-p_2)}{|p_1-p_2|}\cdot dp_2$$
So by the rule (17):
$$\nabla_{p_1}f=\frac{p_1-p_2}{|p_1-p_2|},\ \nabla_{p_2}f=\frac{-(p_1-p_2)}{|p_1-p_2|}$$
In order to learn the theory behind the differentials, you can read, for example, the book  "Matrix Differential Calculus with Applications in Statistics and Econometrics" by J. Magnus et al.
A: We are given
$$ f(x,y) := |x-y|-d$$

*

*$f$ is not actually differentiable everywhere. It is not differentiable on the set $x=y$.

*Note that $$f(x,y)=f(y,x).\tag{$*$}$$ So in fact, $(\nabla_y f)(a,b) = (\nabla_x f)(b,a)$. Easy to verify by checking the directional derivatives:
$$ (\partial_{y_i} f)(a,b) = \lim_{t\downarrow 0} \frac{f(a,b+te_i)-f(a,b)}t\overset{ (*)}= \lim_{t\downarrow 0} \frac{f(b+te_i,a)-f(b,a)}t= (\partial_{x_i} f)(b,a).$$
Once we know this, we only need to compute one of them, say $\nabla_x f(x,y)$.

*Lets consider $y$ as a fixed vector. Then $f$ is the composition of two functions $f=g\circ h$
\begin{align} g&:\mathbb R^n \to \mathbb R &g(z) &:= |z|-d\\
h&:\mathbb R^n \to \mathbb R^n &h(x) &:= x-y
\end{align}
and chain rule gives
$$ \nabla_x f(x,y)^T = \nabla g(h(x))^T \nabla h(x)$$

*The derivatives are
$$ \nabla g(z) = \frac z{|z|}, \ (z\neq 0), \quad \nabla h(x) = I$$

*putting it together
$$  \nabla_x f(x,y)^T = \nabla g(h(x))^T \nabla h(x) = \frac{h(x)^TI}{|h(x)|} = \frac{(x-y)^T}{|x-y|}$$
taking transposes gives the answer.

