You need to be more clear about your double integral. Say you have
$$ \int_c^d \left(\int_a^b f(x,y)g(x,y)dx\right) dy$$
And you need to know the antiderivative of $g(x,y)$ with respect to $x$. So the information $\int_X g(x,y)dx=w(y)$ is not enough. Because this is not an antiderivative of $g$ with respect to the $x$ direction. Instead, you need to have
$$ \int_a^x g(s,y)dy=h(x,y)$$
That is,
$$ \int_a^b g(x,y)dy=h(b,y)-h(a,y)=h(b,y)=w(y)$$
Under this assumption, you can mimic the way how we do integration by parts in 1 dimensional.
$$ \int_c^d\left(\int_a^b f(x,y)g(x,y)dx\right)dy=\int_c^d\left(f(x,y)h(x,y)\Big|_a^b-\int_a^b h(x,y)\frac{\partial}{\partial x}f(x,y)dx\right)dy$$
$$ =\int_c^d f(b,y)w(y)dy-\int_c^d\left(\int_a^b h(x,y)\frac{\partial}{\partial x}f(x,y)dx\right)dy$$
If you know the antiderivative of $h(x,y)$ with respect to $y$, i.e.
$$ \int_c^y h(x,t)dt=\int_c^y\int_a^x g(s,t)ds\, dt=l(x,y)$$
you can keep going:
$$ \int_c^d\left(\int_a^b f(x,y)g(x,y)dx\right)dy=\int_c^d f(b,y)w(y)dy-\int_c^d\left(\int_a^b h(x,y)\frac{\partial}{\partial x}f(x,y)dx\right)dy=\int_c^d f(b,y)w(y)dy-\int_a^b\left(l(x,d)\frac{\partial}{\partial x}f(x,d)-\int_c^d l(x,y)\frac{\partial^2}{\partial x\partial y}f(x,y)dy\right)dx=...$$
It all depends on what information you have for $f(x,y)$ and $g(x,y)$.
Part $a$ diverges for the reason you found: the integral blows up logarithmically as $\rho\rightarrow 0$. Since the integrand is strictly nonnegative, there's no need to worry about cancellations between positive and negative.
Part $b$ is more interesting. You can set it up as an iterated integral in cartesian or polar, integrating $x$ first, $y$ first, $\rho$ first, or $\varphi$ first. Try it and see what you get.
If you're having trouble with the cartesian integral, here's a handy integral table:
$$
\int\frac{x-y}{(x^2+y^2)^{3/2}}dx = -\frac{x + y}{y \sqrt{x^2 + y^2}}
$$
$$
\int \frac{\sqrt{1+y^2}-1-y}{y\sqrt{1+y^2}}dy = \ln\left(\frac{1+\sqrt{1+y^2}}{y+\sqrt{1+y^2}}\right)
$$
Best Answer
To get rid of the absolute value signs in $$\frac{x^2-y^2}{(x^2+y^2)^2}$$ we need to rule out where in $[0,1]\times [0,1]$ this is positive and negative. But $x^2-y^2=(x-y)(x+y)$. The term $x+y$ is always nonnegative, while $x-y$ is positive or negative according as $x>y$ or $x<y$. This is why in the article, they split the integral along the two triangles the line $x=y$ divides $[0,1]\times [0,1]$ into, and leave the appropriate sign in the integrand.