Here is the integral:
May you please suggest some beautiful idea on using Riemann surface, or some Gauss-Ostrogradsky at the beginning. Also, the initial integral looks really symmetric, so maybe there is some way to simplify form or switch to spherical coordinates...
And here are steps which I performed by myself:
By integrating over φ, this can be reduced to one variable integral:
Now I switch to complex plane:
Now poles (there are 4, all lying on real axis, 2 of them are inside contour |z|=1:
And suddenly, here comes sensation that there is root of polynomial, so there is branching at this points.
I've put this equation into Wolfram, but it returns something alienish: This integral on integrals.wolfram.com. Maybe it is possible to calculate all that limits at points θ=0 and θ=π. But there will be no way to prove that Mathematica solution, which is necessary for me :(
Thank you for your attention.