Consider a function of a complex variable $\omega$ which is given by $f(\omega) = e^{-\omega^2/2}$.
This function is symmetric, holomorphic everywhere, and vanishes as $|\omega| \rightarrow \infty$. Thus the Kramers-Kronig relations tell us that, given the real part we can find the imaginary part:
$$ \Im[f(\omega)] = \frac{2\omega}{\pi} \mathcal{P} \int_{0}^{\infty} d\omega' \frac{\Re[f(\omega')]}{\omega'^2 - \omega^2} $$
Here $\mathcal{P}\int$ represents the Cauchy principal value.
We expect of course to find that $\Im[f(\omega)] = 0$ for real $\omega$, but Mathematica seems to disagree. Attempting to solve this numerically, we have
f = Exp[-x^2/2]
g[y_] := (2 y)/Pi * NIntegrate[-(Re[f]/(x^2 - y^2)), {x, 0, y, Infinity}, Method -> "PrincipalValue"]
Plot[{g[y], Im[f /. x -> y]}, {y, 0, 5}]
we obtain an output
which clearly is not zero everywhere. I'm inclined to think that this a problem with how I've written my Mathematica code, but I really can't work my head around it. Mathematica doesn't complain about convergence, the function satisfies all the constraints (I think!) and I'm fairly certain I've specified the integral correctly (i.e. giving the excluded singular point in the integration region).
What's going on?