Different possible solutions for the wave equation? The Wave equation is:
$$\nabla^2\psi(\mathbf{x},t)-\frac{1}{c}\frac{\partial^2 \psi(\mathbf{x},t)}{\partial t^2}=f(\mathbf{x},t)$$
The Green function is then
$$\nabla^2G(\mathbf{x},t)-\frac{1}{c}\frac{\partial^2 G(\mathbf{x},t)}{\partial t^2}=\delta(\mathbf{x}-\mathbf{x}')\delta(t-t')$$
Using the Fourier transform
$$(-k^2+\frac{\omega^2}{c})G(\mathbf{k},\omega)=e^{i\mathbf{k}\cdot\mathbf{x}'}e^{-i\omega t'}$$
Then
$$G(\mathbf{x},t)=\int_{\mathbb{R}^4}\frac{e^{i\mathbf{k}\cdot(\mathbf{x}'-\mathbf{x})}e^{-i\omega(t'-t)}}{(-k^2+\frac{\omega^2}{c})}d^3k\,d\omega$$
Let's consider only the integral in the frequency space
$$I_{\mathbb{R}}=\int_{\mathbb{R}}\frac{e^{-i\omega(t-t')}}{(\omega^2-c^2k^2)}d\omega$$
The integrand has two simple poles in $\omega=\pm ck$ and has to be solved using the Residue Theorem. Let's choose the retarded solution ($t'<t$).
If my path ($\Gamma$) is just a closed semi-circle the poles are not inside the path, so I have two options:


*

*Move the poles: $\pm ck = \displaystyle\lim_{\varepsilon \to 0}\pm ck-i\varepsilon$ with $\varepsilon\in \mathbb{R}^+$

*Add two small semicircles to the path


Let's follow method (1) first, choosing the path with the semicircle ($SC$) closing in the lower complex plane, its contribution is zero.
$$I_{\Gamma}=I_{\mathbb{R}}+I_{SC}=I_{\mathbb{R}}=-2\pi i\left(\mathrm{Res}\left(\frac{e^{-i\omega(t-t')}}{(\omega^2-c^2k^2)},+ck\right)+\mathrm{Res}\left(\frac{e^{-i\omega(t-t')}}{(\omega^2-c^2k^2)},-ck\right)\right) \tag{a}$$
But if I move the poles in the upper plane (i.e. $\pm ck = \displaystyle\lim_{\varepsilon \to 0}\pm ck+i\varepsilon$ with $\varepsilon\in \mathbb{R}^+$) there are no poles in the path so
$$I_{\Gamma}=I_{\mathbb{R}}+I_{SC}=I_{\mathbb{R}}=0 \tag{b}$$
Following instead method (2):
If I choose the small semicircles ($sc_\pm$) to be in the upper plane, then the poles are inside $\Gamma$ and then 
$$I_{\Gamma}=I_{\mathbb{R}}+I_{SC}+I_{sc_-}+I_{sc_+}=I_{\mathbb{R}}+I_{sc_-}+I_{sc_+}=-2\pi i\left(\mathrm{Res}\left(\frac{e^{-i\omega(t-t')}}{(\omega^2-c^2k^2)},+ck\right)+\mathrm{Res}\left(\frac{e^{-i\omega(t-t')}}{(\omega^2-c^2k^2)},-ck\right)\right)$$
so
$$I_{\mathbb{R}}=I_{\Gamma}-I_{sc_-}-I_{sc_+}=-\pi i\left(\mathrm{Res}\left(\frac{e^{-i\omega(t-t')}}{(\omega^2-c^2k^2)},+ck\right)+\mathrm{Res}\left(\frac{e^{-i\omega(t-t')}}{(\omega^2-c^2k^2)},-ck\right)\right) \tag{c}$$
because 
$$ I_{sc_\pm} = -\pi i\mathrm{Res}\left(\frac{e^{-i\omega(t-t')}}{(\omega^2-c^2k^2)},\pm ck\right)$$
If I choose the small semicircles ($sc_\pm$) to be in the lower plane, there are no poles inside $\Gamma$ and then I'm getting again
$$I_{\mathbb{R}}=I_{\Gamma}-I_{sc_-}-I_{sc_+}=-\pi i\left(\mathrm{Res}\left(\frac{e^{-i\omega(t-t')}}{(\omega^2-c^2k^2)},+ck\right)+\mathrm{Res}\left(\frac{e^{-i\omega(t-t')}}{(\omega^2-c^2k^2)},-ck\right)\right) \tag{d}$$
because in this case
$$ I_{sc_\pm} = \pi i\mathrm{Res}\left(\frac{e^{-i\omega(t-t')}}{(\omega^2-c^2k^2)},\pm ck\right)$$
I've then found three different solutions, the first one should be the correct one, but why? How do I choose it? Am I doing something wrong?

Maybe this analogy shouldn't be made, but wikipedia claims that in the resolution of the Klein-Gordon equation (very similar to the wave equation) is the same thing to add a small semicircle around the pole or modify the integrand adding a small term $\varepsilon \to 0$ to the poles. Why my case should be different? 

 A: The first one, as you note, is the expected answer. The second and third answers are in fact the same, and there you've computed the Cauchy principal value of the integral. The difference lies primarily in how you're taking the limit of some $\epsilon \rightarrow 0$. 
In the first case, you're changing the integrand slightly, and in the second/third case you're changing the domain of integration slightly. There's no reason to expect these integrals to be the same.
In fact, one can view the first and third approaches to be the same so long as one includes the semi-circle portions in the third answer. Similarly, if one does this for the second answer one gets the advanced propagator.
At the end of the day, what we have really learned is that extracting sensible answers from divergent integrals is a bit arbitrary. One should think of these various prescriptions either as tricks or mnemonics to get the correct answer but not as mathematically rigorous solutions to the problem. In other words, there are many possible ways to regularize a divergent integral to extract a convergent answer, but in order to get the correct "physical" answer corresponding to the original problem only particular prescriptions are allowed. 
There are ways to mathematically formalize the problem such that one doesn't have to resort to such tricks to get the desired answer. As noted by DanielSank, one can encode the correct prescription/get the right result rigorously by enforcing the proper boundary conditions. I am not familiar with how to do this, unfortunately.
