I'm looking for the above commutator where I suppose $p=|\vec{p}|$. Using $$[A,BC]=[A,B]C+B[A,C] \, ,$$ I came up with $$[x_i,p_j p^{-1}]=[x_i,p_j]p^{-1} + p_j[x_i,p^{-1}] = i \delta_{ij} p^{-1} + p_j[x_i,p^{-1}] \tag{0}$$ where I used the canonical commutator $[x_i,p_j]=i\delta_{ij}$ (pre-factor is the imaginary unit). For the second commutator I have $$0=[x_i,pp^{-1}]=[x_i,p]p^{-1}+p[x_i,p^{-1}]$$ $$\Rightarrow \qquad [x_i,p^{-1}]=-p^{-1}[x_i,p]p^{-1} \, , \tag{1}$$ which leaves me with $$2ip_i=[x_i,p^2]=[x_i,p]p + p[x_i,p]$$ and now I'm stuck. I see that if $[x_i,p]=i p_i/p$, both sides work out, but it's not really a formal proof.
Is it really necessary to write for each component (for instance $i=1$) $$p=\sqrt{p_1^2+p_2^2+p_3^2}=\sqrt{p_2^2+p_3^2} \sqrt{1+\frac{p_1^2}{p_2^2+p_3^2}}=\sqrt{p_2^2+p_3^2} \sum_{k=0}^\infty \binom{1/2}{k} \left( \frac{p_1^2}{p_2^2+p_3^2} \right)^k$$ and use the fact that $$[x_1,p_2]=[x_1,p_3]=0 \\ [x_1,p_1^{2k}]=i2k\,p_1^{2k-1}=i\frac{{\rm d}}{{\rm d}p_1}p_1^{2k} \, ?$$
This gives $$[x_i,p]=i\frac{{\rm d}}{{\rm d}p_i} p=i p_i/p \qquad \stackrel{(1)}{\Rightarrow} \qquad [x_i,p^{-1}]=-i p_i/p^3 \, .$$ Surely, I could have used the binomial expansion in (1) directly, but that's besides the point.
Are there other approaches, without relying on the series expansion?
The canonical commutation relations $$[x_i, p_j]=i\delta_{ij} \tag{1} \label{eq1}$$ (valid on a suitable common domain for the components of the position and the momentum operator) imply $$[x_i, f(\vec{p})]= i \frac{\partial}{\partial p_j} f(\vec{p}), \tag{2} \label{eq2}$$ for a suitable differentiable function $f$. (Remember that the position operator generates translations in momentum space and vice versa.)
The validity of \eqref{eq2} can easily seen by realizing \eqref{eq1} in the momentum representation with $x_i \to i \partial /\partial p_i$ and $p_j$ as a multiplication operator acting on momentum-space wave functions $\psi(\vec{p})$. Acting on $\psi(\vec{p})$, \eqref{eq2} becomes indeed $$[x_i, f(\vec{p})] \psi(\vec{p})=i \frac{\partial}{\partial p_i} \left(f(\vec{p}) \psi(\vec{p})\right)- i f(\vec{p}) \frac{\partial}{\partial p_i} \psi(\vec{p}) = i \frac{\partial f(\vec{p})}{\partial p_i} \psi(\vec{p}). \tag{3} \label{eq3}$$ Note that the position-space representation (where $x_i$ acts as a multiplication operator on wave functions $\phi(\vec{x})$ and $p_i \to -i \partial/\partial x_i$) and the momentum-space representation are equivalent as they are connected a unitary mapping, the Fourier tranformation $\mathcal{F}: L^2(\mathbb{R}^3) \to L^2(\mathbb{R}^3)$ with $$\mathcal{F}^{-1} x_i \mathcal{F} = i \partial/\partial p_i, \quad \mathcal{F}^{-1} (-i \partial/ \partial x_i) \mathcal{F} = p_i.$$ For showing \eqref{eq2}, the momentum-space representation turns out to be simply more convenient.
In your case, $f(\vec{p}) = p_j/|\vec{p}|$, one finds thus the result $$\left[x_i, \frac{p_j}{|\vec{p}|}\right]=i\left( \frac{\delta_{ij}}{|\vec{p}|}-\frac{p_i p_j}{|\vec{p}|^3} \right). \tag{4} \label{eq4}$$