(Fundamental) Solution of the Helmholtz equation The fundamental solution of the Helmholtz equation in $\mathbb{R}^3$
$$(\Delta+k^2)u=-\delta \tag{1}$$
is well known: 
$$u(x)=\frac{e^{\pm ik|x|}}{4\pi |x|}$$
solves the Helmholtz equation in distributional sense. The usual ansatz to obtain fundamental solutions is to Fourier transform both sides. Then $(1)$ becomes
$$(-|x|^2+k^2)\hat{u}(x)=-1 \implies \hat{u}(x)=\frac{1}{k^2-|x|^2}.$$
The problem now arises is that we want to inverse Fourier transform $\hat{u}$ to obtain the solution. But $\hat{u}$ has singularities on the sphere of radius $k$.
How do we proceed in a rigorous distributional way to obtain the above fundamental solution? What do people usually do and is there any book about such problems?
 A: 
There were some errors in the OP.  First, the Fourier Transform of the Helmholtz Equation (Equation $(1)$ in the OP) is given by
$$(-k'^2+k^2)\hat u(\vec k')=-1$$
where $\vec k'$ is the $3$-D transform variable with $k'=|\vec k'|$.
Second, it then follows that $$\hat u(\vec k')=\frac{1}{k'^2-k^2}\tag 1$$  In what follows, we  propose a way forward to take the inverse Fourier Transform of $\hat u(\vec k')$ as given by $(1)$.


METHODOLGY $(1)$:
One way to handle the singularities at $k'=\pm k$ is to let $k$ have a "small," non-zero imaginary part.  We will assume that $\text{Im}(k)<0$.
Then, we can write
$$\begin{align}
u(\vec r)&=\frac{1}{(2\pi)^3}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \frac{e^{i\vec k'\cdot \vec r}}{k'^2-k^2}\,dk'_x\,dk'_y\,dk'_z\\\\
&=\frac{1}{(2\pi)^3}\int_0^{2\pi}\int_0^\pi\int_0^\infty \frac{e^{ik'r\cos(\theta)}}{k'^2-k^2}\,k'^2\sin(\theta)\,dk'\,d\theta\,d\phi\\\\
&=\frac{1}{(2\pi)^2}\int_0^\pi\int_0^\infty \frac{e^{ik'r\cos(\theta)}}{k'^2-k^2}\,k'^2\sin(\theta)\,dk'\,d\theta\\\\
&=\frac{1}{(2\pi)^2}\int_0^\infty \frac{k'}{k'^2-k^2}\frac{2\sin(k'r)}{r}\,dk'\\\\
&=\frac{1}{(2\pi)^2}\int_{-\infty}^\infty \frac{k'}{k'^2-k^2}\frac{\sin(k'r)}{r}\,dk'\\\\
&=\frac{e^{-ikr}}{4\pi r}
\end{align}$$
Had we assumed that $\text{Im}(k)>0$, we would have recovered the solution $\frac{e^{ikr}}{4\pi r}$.
Finally, we can let the imaginary part of $k$ approach $0$ and recover the result for real $k$.


Recommended references include (i) "Field Theory of Guided Waves," Robert Collin, IEEE Press, (ii) "Waves and Fields in Inhomogeneous Media," W.C. Chew, Van Nostrand Reinhold, and (iii) Radiation and Scattering of Waves, Felsen and Marcuvitz, IEEE Press.


METHODOLGY $(2)$:
We could solve Equation $(1)$ in the OP without the use of integral transformation.  Instead, we write
$$\nabla^2 u(\vec r)+k^2u(\vec r)=0$$
for $\vec r\ne 0$.  Exploiting the spherical symmetry of the problem, we have
$$\nabla^2 u+k^2 u=\frac1{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial u}{\partial r}\right)+k^2u=0$$
which can be rearranged as
$$\frac{\partial^2}{\partial r^2}(r u)+k^2 (ru)=0 \tag 2$$
Solutions to $(2)$ are trivially seen to be
$$u=C^{\pm}\frac{e^{\pm ikr}}{r}$$
We find the constant $C^\pm$ by enforcing the condition $\int_{|\vec r|\le \epsilon}\nabla^2 u(\vec r)\,dV=-1$.
Applying the Divergence Theorem in the sense of distributions reveals
$$\begin{align}
\oint_{|\vec r|=\epsilon}\left.\left(\frac{\partial u(\vec r)}{\partial r}\right)\right|_{|\vec r|=\epsilon}\,\epsilon^2 \sin(\theta)\,d\theta\,d\phi&=-4\pi C^{\pm}\\\\
&=-1\end{align}$$
whereupon solving for $C^{\pm}$ yields $C^\pm=\frac{1}{4\pi}$.
