How comes that integrating factor = $1/(x\,y)$
You may have to guess, but in this case you can make it an educated guess.
If you multiply by $dx$ and thus rewrite the differential as
$M dx + N dy=(x+y)dx - (x^2/y)dy$
you see that if $x$ and $y$ were multiplied by a constant $\lambda$, you would get just the same differential multiplied by a constant:
$(\lambda x+\lambda y)d(\lambda x) - [(\lambda x)^2/\lambda y]d(\lambda y) = (\lambda^2)[(x+y)dx - (x^2/y)dy]$
A differential with this property is called homogeneous.
Homogeneous differentials tend to lend themselves to integrating factors having the form $x^ay^b$, which is itself homogeneous. With that in mind let us try this type of integrating factor. We then have the differential
$P dx +Q dy=x^ay^b(x+y)dx - x^ay^b(x^2/y)dy=(x^{a+1}y^b+x^ay^{b+1})dx - x^{a+2}y^{b-1}dy$
Applying the exactness criterion we have
$\partial P/\partial y = \partial Q/\partial x$
$bx^{a+1}y^{b-1}+(b+1)x^ay^b=-(a+2)x^{a+1}y^{b-1}$
Combining like terms reduces the exactness criterion to
$(a+b+2)x^{a+1}y^{b-1}+(b+1)x^ay^b=0$
So we must have $a+b+2=0,b+1=0$. From this we get $b=-1$ and then $a=-1$, so we find the integrating factor $x^ay^b=1/(xy)$.