Pairwise $L_p$ distance matrix between rows of x1 and x2.

Specialised paths exist for p in $\{1, 2, \infty, 0\}$; other p values fall back to the general $(\sum_k |a_k - b_k|^p)^{1/p}$ formula. For p=2 the numerically stable expansion $\|a-b\|^2 = \|a\|^2 + \|b\|^2 - 2 a \cdot b^\top$ is used so large coordinates don't lose precision.