Torque on Non-Uniform Sphere

Laying out steps of evaluation. Location of point mass:

$\vec{r}_c = (1,2,1)$

Location of any point on the sphere in spherical coordinates:

$\vec{r}_p = (r \sin{\theta} \cos{\phi},r \sin{\theta} \sin{\phi},r \cos{\theta} )$

Force acting on a volume mass element near the point on the sphere:

$d\vec{F} = \frac{G \rho m \ dV \left(\vec{r}_c-\vec{r}_p\right)}{\lvert \vec{r}_c-\vec{r}_p \rvert^3}$

The density is written in term of spherical coordinates and:

$dV = r^2 \sin{\theta} \ dr \ d\theta \ d\phi$

The torque due to this elementary force about the origin is:

$d\vec{\tau} = \vec{r}_p \times d\vec{F}$

Substituting all expressions and simplifying:

$d\vec{\tau} = d\tau_x \ \hat{i} + d\tau_y \ \hat{j} + d\tau_z \ \hat{k}$ $d \vec{\tau} = \left(f_x \ \hat{i} + fy \ \hat{j} + f_z \ \hat{k}\right) dr \ d\theta \ d\phi$

$\tau = \left(\int_{0}^{1} \int_{0}^{\pi} \int_{0}^{2\pi} f_x \ dr \ d\theta \ d\phi\right) \ \hat{i} + \left(\int_{0}^{1}\int_{0}^{\pi} \int_{0}^{2\pi} f_y \ dr \ d\theta \ d\phi\right) \ \hat{j} + \left(\int_{0}^{1} \int_{0}^{\pi} \int_{0}^{2\pi}f_z \ dr \ d\theta \ d\phi\right) \ \hat{k}$ $\tau = \tau_x \ \hat{i} + \tau_y \ \hat{j} + \tau_z \ \hat{k}$

The required answer is:

$\lvert \tau \rvert = \sqrt{\tau_x^2 + \tau_y^2 + \tau_z^2}$

The numerical integration is very time-consuming as it requires high numerical resolution for decent convergence. I consider myself lucky to have got the answer.