Integration of multiplication of two generalized function I know the generalized function relation $\frac
1{x-\omega+i\epsilon}=\frac 1{x-\omega}-\pi i \delta(x-\omega)$,
where the
first term on the right hand side is understood as principle value. 
However there is an integration which could give different answer when
using this relation differently. The problem is as follows, in one way
\begin{eqnarray}
&&\int d\omega \frac
{f(\omega)}{(x-\omega+i\epsilon)(y-\omega-i\epsilon)}
\\&=&\int d\omega\, f(\omega)\left( \frac
1  {(x-\omega)}-\pi i \delta(x-\omega)\right)\left( \frac
1  {(y-\omega)}+\pi i \delta(y-\omega)\right)
\\&=&\int d\omega \frac{f(\omega)}{(x-\omega)(y-\omega)}-\pi i \frac
{f(x)}{y-x}+\pi i \frac {f(y)}{x-y}+\pi^2\delta(x-y) f(x)
\end{eqnarray}
and in another way,
\begin{eqnarray}
&&\int d\omega \frac
{f(\omega)}{(x-\omega+i\epsilon)(y-\omega-i\epsilon)}
\\
&=&\frac 1
{y-x-i\epsilon}\int d\omega\,f(\omega)\left( \frac
1  {x-\omega+i\epsilon}-\frac
1  {y-\omega-i\epsilon}\right)
\\&=&
\left(\frac 1
{y-x}+\pi i\delta(y-x)\right)\int d\omega\,f(\omega)\left( \frac
1  {(x-\omega)}-\pi i \delta(x-\omega)- \frac
1  {(y-\omega)}-\pi i \delta(y-\omega)\right)
\\&=&\frac 1{y-x}\int d\omega f(\omega)\left(\frac1{x-\omega}-\frac
1{y-\omega}\right)-\pi i \frac
{f(x)}{y-x}+\pi i \frac {f(y)}{x-y}+2\pi^2\delta(x-y)
\\&=&\int d\omega \frac{f(\omega)}{(x-\omega)(y-\omega)}-\pi i \frac
{f(x)}{y-x}+\pi i \frac {f(y)}{x-y}+2\pi^2\delta(x-y)f(x)\,.
\end{eqnarray}
These two ways give different answers. Why and which is correct? In fact, I really want the
second one to be correct. 
And another question, how to do the other integration if both sign of
the $i \epsilon$ are the same in the original integral? That is,
\begin{eqnarray}
\int d\omega \frac
{f(\omega)}{(x-\omega+i\epsilon)(y-\omega+i\epsilon)}
\end{eqnarray}
Thanks.
 A: Short answer: I can only point out some of the conceptual problems. The initial expression does make sense as a distribution in the sole $\omega$ variable but already in the result of the first calculation, one has to understand it as a distribution in $x$ and $y$ also. This is even more manifest in the second line of the second calculation where $\frac{1}{y-x - i\epsilon}$ is introduced although for $x\neq y$ fixed (not variables), the $i\epsilon$ could have been dropped. Now as a distribution in $(x,y),\ i\pi \frac{f(x)}{y-x}$ and $i\pi \frac{f(y)}{x-y}$ are ambiguous: to bring this back to the single variable case $\frac{1}{x}$ is ambiguous as a distribution, it could be $\frac{1}{x\pm i\epsilon}$ or $\operatorname{pv}\left(\frac{1}{x}\right)$. (A very interesting question is "what are all the possibilities? as a distribution or a a tempered distribution". This reminds me of the discussion Example 9 p.140 in "Methods of Modern Mathematical Physics 1: Functional analysis", Michael Reed, Barry Simon). These terms are obtained by multiplying a $\delta$ and a $\operatorname{pv}$ so they should be $\operatorname{pv}$.


*

*Puting $\operatorname{pv}$ in the three first terms, the first calculation seems correct to me. (Nevertheless, we should check if $\operatorname{pv}\left(\frac{1}{x-\omega}\right) \times \operatorname{pv}\left(\frac{1}{y-\omega}\right) = "\operatorname{pv}"\left(\frac{1}{(x-\omega)(y-\omega)}\right)$, in particular what is the definition of r.h.s.).

*while for the second computation, we should check that $$\operatorname{pv}\left(\frac{1}{y-x}\right) \times \left( \operatorname{pv}\left(\frac{1}{x-\omega}\right) -\operatorname{pv}\left(\frac{1}{y-\omega}\right) \right) = "\operatorname{pv}"\left(\frac{1}{(x-\omega)(y-\omega)}\right)$$
All in all, I'm not even sure whether your two results are really different. 
Edit n°20000: it seems that as a distribution in $\omega$ and when $x\neq y$, the $\delta(x-y)$ should be replaced by $0$. Remains to see what happens when $x=y$.
To Do: We should clarify the "prescription" for the distributions $\frac{1}{(x-\omega)(y-\omega)}, \frac{1}{(x-\omega)}$ and $\frac{1}{(y-\omega)}$ and check if the different expressions coincide when "smeared" with test functions in $x$ and $y$.

Long answer: In a first part, I explain what could be a multiplication of distributions. Then I try to no avail to get some explicit formulae (which are not given by theorems).
There are several attempts at defining multiplication of distributions (at the same point!). One by Jean-François Colombeau who has published a book with exactly that title (about which I know nothing except the title...) and (at least) a second approach by Hörmander. Of course the product cannot be defined for any pair of distributions but he gives a sufficient condition.
Remark: (cf. Cours d'analyse, Théorie des distributions et analyse de Fourier, Jean-Michel Bony, Remarque 9.5.14 p.183. He talks about extending a bilinear map from $\mathcal{A} \times \mathcal{B}$. Contrary to linear maps e.g. Fourier transform that is extended from a space to a bigger one in which the previous set is dense, one cannot just look for the "biggest" space in which the bilinear could be defined. Indeed if one chooses a "big" space for $\mathcal{A}$ then one may have to take a more restrictive space for $\mathcal{B}$, and there are hence different choices of pairs $\mathcal{A} , \mathcal{B}$.)
Here is the idea of Hörmander (I suppose...): for $f,g \in L^2 \big(\mathbb{R}^n\big)$ or $\mathcal{S}\big(\mathbb{R}^n\big)$ or the Wiener algebra, or any space on which the Fourier transform is defined and where one can take the inverse, one has
$$ \mathcal{F}(f g) = \mathcal{F}(f)\ast \mathcal{F}(g) \quad \text{so}\quad fg = \mathcal{F}^{-1} \Big(\mathcal{F}(f)\ast \mathcal{F}(g) \Big) \tag{0}\label{Idea}$$
Of course one already knows how to multiply $fg$ in these spaces, but the idea is to see if the r.h.s. still makes sense for $f,g$ distributions. 
Now comes two subtleties: 


*

*Fourier transform for distributions are only defined for tempered ones (they are defined by $\langle \mathcal{F}(T),\varphi \rangle :=\langle T, \mathcal{F}\varphi \rangle$ but then one should check that $\mathcal{F}\varphi$ is still a test function. It turns out that the Fourier transform of a function with compact support is analytic (Paley-Wiener), it cannot be of compact support unless it is constant at $0$. Therefore one takes Schwartz' space as the space of test functions which is stable under $\mathcal{F}$).

*One will in fact define the product $TS$ of two distributions $T$ and $S$ in $\mathcal{D}'(\mathbb{R}^n)$ (the usual space of distributions, not tempered) "locally". Indeed for these, if one knows the restriction $T|_{\mathcal{U}}: \mathcal{C}^{\infty}_c(\mathcal{U}) \longrightarrow \mathbb{C}$ for a family of open subsets $\mathcal{U}$ covering $\mathbb{R}^n$ then $T:\mathcal{C}^{\infty}_c(\mathbb{R}^n) \longrightarrow \mathbb{C}$ is uniquely defined. This local definition does not work for tempered distributions...


Here is the theorem (sufficient condition), from "Methods of Modern Mathematical Physics 2: Fourier Analysis, Self-adjointness", Michael Reed, Barry Simon, Theorem IX.45 p.95: (first a def.)
Def: Let $T\in \mathcal{D}'(\mathbb{R}^n)$ be a distribution. A point $(x_0,k_0)\in \mathbb{R}^n \times \mathbb{R}^n\backslash\{0\}$ is regular directed if there exists a neighborhood $U\times V$ of it and a "plateau" function ($\chi\in \mathcal{C}^{\infty}_c(\mathbb{R}^n)$, identically $1$ on $U$) such that for any $m\in \mathbb{N}^*,\ \exists\, C_m > 0$,
 $$ \lvert \mathcal{F}(\chi T)(\lambda k) \rvert \leq \frac{C_m}{(1+ \lambda)^m},\quad \forall\, k\in V,\ \lambda > 0 $$
The wave front set $\operatorname{WF}(T)$ is the complementary in $\mathbb{R}^n \times \mathbb{R}^n\backslash\{0\}$, i.e. the set of points which are not regular.
Remark: If a function $f$ is smooth then its Fourier transform $\mathcal{F}(f)(k)$ decreases faster than any $(1+\lVert k\rVert)^{-m},\ m\in \mathbb{N}$. In fact if for a fixed $x_0\in \mathbb{R}^n,\ (x_0,k)$ is regular directed in all direction $k\in \mathbb{R}^n\backslash \{0\}$ then the "localized" distribution $\chi T$ is a smooth function.
Thm: Let $S,T \in \mathcal{D}'(\mathbb{R}^n)$. If
$$ \operatorname{WF} (S) \oplus \operatorname{WF}(T) := \left\lbrace (x,k_1 + k_2) \in \mathbb{R}^n \times \mathbb{R}^d,\ (x,k_1) \in  \operatorname{WF} (S) , (x,k_2) \in \operatorname{WF} (T) \right\rbrace $$
does not contain any element of the form $(x,0)$ (i.e. one never has $k_2 = - k_1$) then the product is well-defined as the unique distribution $ST\in \mathcal{D}'(\mathbb{R}^n)$ such that for any $x\in \mathbb{R}^n,\ \exists\ U$ neighborhood of $x$ and $ \chi \in \mathcal{C}^{\infty}_c(\mathbb{R}^n)$ identically 1 on $U$ such that (p.90)
$$ \mathcal{F}(\chi^2 ST)(k) = \mathcal{F}(\chi S)\ast \mathcal{F}(\chi T)(k) , \quad \forall k\in \mathbb{R}\tag{1}\label{Hoermander}$$
and $\ \operatorname{WF} (ST) \subseteq \operatorname{WF} (S) \cup \operatorname{WF} (T) \cup \big(\operatorname{WF} (S) \oplus \operatorname{WF}(T) \big)$.
Remark: Writing $\pi_1: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n, (x,k) \mapsto x$ the projection on the first component, one has $\pi_1 \big(\operatorname{WF} (S) \oplus \operatorname{WF}(T) \big) \subseteq \pi_1 \big(\operatorname{WF} (S) \big) \cap \pi_1 \big(\operatorname{WF}(T) \big)  $. The conditions is only about points where both $S$ and $T$ are singular. Indeed outside, there are regular or more precisely, their Fourier transform do not "decrease slowly" in the same direction, so in these points the definition (\ref{Hoermander}) of $ST$ is not problematic.
Actually when one looks in the book "The Analysis of Linear Partial Differential Operators 1 - Distribution Theory and Fourier Analysis" (1998), Lars Hörmander, Theorem 8.2.10 p. 267, the product is defined as the restriction of the tensor product $T\otimes S="T(x)\, S(y)"$ (which would take test functions $\varphi(x,y)$) to the "diagonal" $(x,x) \in \mathbb{R}^n\times \mathbb{R}^n$ (i.e. the couples $(x,y)\in\mathbb{R}^n\times \mathbb{R}^n$ such that $ y=x$ and thus the test functions one should take are functions of $x$ only.). I think the conditions are the same.
Two brief introductory texts are A smooth intro to WFS, Brouder, Dang, Hélein and Fronts d'onde et op. sur les distrib., J.M. Bony

Here are the wave front sets of the different distributions:


*

*$\mathcal{F}\left(\frac{1}{x \pm i \epsilon}\right)(k) = \mp 2i\pi \theta(\pm k)$ with the convention $\mathcal{F}(f)(k):= \int e^{-ikx}f(x)\, dx$. Such a function of $k$ is either $0$ or $\pm 2 i \pi$ depending on the sign of $k$, i.e. decreases faster than any polynomial  or does not. Hence
$$\operatorname{WF}\left( \frac{1}{x \pm i \epsilon} \right)= \Big\lbrace (0,k)\in \mathbb{R}\times \mathbb{R}^*,\ \pm k > 0 \Big\rbrace \tag{2}\label{WFS1}$$
(cf. e.g. Brouder, Dang, Hélein, Example 5 & 16 p.7 & 13, where unfortunately they take a very strange convention, so mind the signs...)

*Similarly, $\mathcal{F}(\delta)(k) =1$, no decrease in either positive or negative direction:
$$\operatorname{WF}\left( \delta \right)= \Big\lbrace (0,k)\in \mathbb{R}\times \mathbb{R}^* \Big\rbrace \tag{3}\label{WFS2}$$ 

*By the so-called Sokhotski-Plemelj formula mentionned by the OP, $ pv\left(\frac{1}{x} \right) = \frac{1}{2}\left(\frac{1}{x + i \epsilon} + \frac{1}{x - i \epsilon} \right)$
$$ \mathcal{F}\left( pv\left(\frac{1}{x} \right) \right)(k) =\frac{1}{2}\left(\mathcal{F}\left(\frac{1}{x + i \epsilon}\right)(k) + \mathcal{F}\left(\frac{1}{x - i \epsilon}\right)(k) \right) = \frac{1}{2}\left(- 2i\pi\, \theta(k) + 2 i \pi\, \theta(-k)\right)= -i \pi \operatorname{sign}(k)$$
where $\theta$ denotes the Heaviside function. One could also first show that $\frac{x}{x^2 + \epsilon} \overset{\text{weak}}{\underset{\epsilon \to 0}{\longrightarrow}} pv\left(\frac{1}{x}\right)$ and use the Fourier transform of the so-called Cauchy-Lorentz distribution (or actually its first moment), cf. also this answer. So as for the $\delta$ distribution, one has
$$\operatorname{WF}\left( pv\left(\frac{1}{x} \right) \right)= \Big\lbrace (0,k)\in \mathbb{R}\times \mathbb{R}^* \Big\rbrace \tag{4}\label{WFS3}$$



*

*When $x\neq y$, the distribution in the $\omega$ variable $\frac{1}{(x-\omega+i\epsilon)(y-\omega-i\epsilon)}$ is well-defined: indeed $\frac{1}{x-\omega+i\epsilon}=- \frac{1}{\omega -x -i\epsilon}$ has wave front set
$$\operatorname{WF}\left(\frac{1}{x-\omega+i\epsilon}\right) = \big\lbrace (x,k),\ k<0 \big\rbrace $$
while $\operatorname{WF}\left(\frac{1}{y-\omega-i\epsilon}\right)=\big\lbrace (y,k),\ k>0 \big\rbrace $ so that (here, assume $x\neq y$, fixed, not variables)
$$\operatorname{WF}\left(\frac{1}{x-\omega+i\epsilon}\right) \oplus \operatorname{WF}\left(\frac{1}{y-\omega-i\epsilon}\right)= \emptyset $$
So in principle, the original expression should make sense as a distribution in $\omega$ alone. Nevertheless, the first calculation (which perfectly makes sense, understanding $\frac{1}{(x-\omega)(y-\omega)}$ as a product of pv, principal value, whose singularities are at $x$ and $y \neq x$) produces a $\pi^2 \delta(x-\omega)\, \delta (y-\omega)$. Again, the condition given by the theorem holds $\operatorname{WF}\big(\delta(x-\omega)\big)\oplus \operatorname{WF}\big(\delta(y-\omega)\big) = \emptyset$ does not contain elements of the form $(\omega,0)$ so the product is defined, nevertheless writing $\pi^2 \delta (x-y)$ implicitly means that we understand it as a distribution in $x$ or $y$ or both, while we have just said that it should make sense as a distribution in $\omega$ alone. Let us try to get more explicit expression, e.g. using (\ref{Idea}))
$$ \delta(x-\omega)\, \delta (y-\omega) := \mathcal{F}^{-1}\left( \int e^{-ixl} e^{-iy(k-l)}\, dl  \right)(\omega)= \mathcal{F}^{-1}\left( e^{-iyk} 2 \pi \, \delta(x-y)\right)(\omega) = \delta(\omega-y)\,2\pi\,\delta(x-y)$$
This does not help and I get a mysterious $2\pi$ factor... Let us try approximations of the identity $\frac{e^{-x^2/\epsilon}}{\sqrt{\epsilon}}\underset{\epsilon \to 0}{\overset{w}{\longrightarrow}} \delta(x)$. Does the following make sense as a distribution in $\omega$ only:
$$ \delta(x-\omega)\, \delta (y-\omega) = \lim_{\epsilon\to 0} \frac{1}{\epsilon} \exp\left(- \frac{(\omega-x)^2 + (\omega-y)^2}{\epsilon}\right) = \lim_{\epsilon\to 0} \frac{1}{\epsilon} \exp\left(- \frac{2\left[\left(\omega-\frac{(x+y)}{2}\right)^2  -\left(\frac{(x+y)}{2}\right)^2\right] + x^2 + y^2}{\epsilon}\right) = \lim_{\epsilon\to 0} \frac{1}{\epsilon} \exp\left(- \frac{2\left(\omega-\frac{(x+y)}{2}\right)^2  +\left(\frac{(x-y)}{2}\right)^2}{\epsilon}\right) =0$$
(recall $x\neq y$). So in fact, one should probably ignore the $\delta (x-y)$ in that case.

*When $x=y$, the initial quantity is not defined (cf. Brouder, Dang, Hélein, Ex. 6 p.8. In fact you see that $-\pi i \frac{f(x)}{y-x}$ and $\pi i \frac {f(y)}{x-y}$ aren't defined either in that case). Therefore one has to understand everything as a distribution in $(\omega,x,y)$ or at least in two variables, e.g. $(\omega,x)$.

*For those who had the patience to read up to here, the last expression of the OP does make sense even when $x=y$, cf. Brouder, Dang, Hélein, Ex. 5 p.7.

