If we consider that
$$\nabla^2\left(\frac{1}{r}\right) = -4\pi\delta(\vec{r})$$
we can explain the dirac-delta function here via the results of Gauss' law $$\int_V \nabla^2\left(\frac{1}{r}\right) = \int_S\nabla\left(\frac{1}{r}\right)\cdot d\vec{a} = \int_S \frac{-\hat{r}}{r^2}\cdot d\vec{a} = -4\pi\frac{R^2}{r^2} $$
where the surface integral has been taken over a sphere of radius $R$. Now, if we say $r>0$ and $R\to 0$, we get $$\int_V \nabla^2\left(\frac{1}{r}\right) = 0. $$ Furthermore, if we set $r=R$ and send $R\to0$ we get $$\int_V \nabla^2\left(\frac{1}{r}\right) = -4\pi $$
Thus we have demonstrated the validity of the first equation. However, if we were to directly calculate the Laplacian of $\frac 1r$ we obtain $$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial}{\partial r}\left(\frac{1}{r}\right)\right) = 0$$
So it seems like the -$4\pi\delta(\vec{r})$ has somehow been neglected by the direct calculation. That being said, what about the case where we have, for some $\epsilon\in\mathbb{R} > 0$
$$\nabla^2\left(\frac{1}{r^\epsilon}\right) = \frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial}{\partial r}\left(\frac{1}{r^\epsilon}\right)\right) = \frac{(\epsilon-1)\epsilon}{r^{\epsilon+2}}$$
Has there also been a dirac delta function neglected here? In general, when does one need to consider contributions due to the dirac-delta function?