How to compute $\int_{\mathbb{C}} \frac{|dz|^2}{|z-a_1| \cdot |z-a_2| \cdot |z-a_3|}$ explicitly? I have these integrals :
$$I_1= \int_{\mathbb{C}} \frac{|dz|^2}{|z-a_1| \cdot |z-a_2| \cdot |z-a_3|},$$
and
$$I_2= \int_{\mathbb{C}} \frac{|dz|^2}{|z-a_1| \cdot |z-a_2| \cdot |z-a_3| \cdot |z-a_4|}.$$
wherer $a_i \in \mathbb{C}$, and the $a_i$ are distincts.
I can see that those integrals are well defined, but I need to evaluate them. Can they be computed explicitly ?
Note : I would like an explicit value for these integrals, but failing that I am interested in the asymptotic value when $a_2 \rightarrow a_1$.
 A: Let us consider the first integral. Make the change of variables $z=a_1+(a_2-a_1)\xi$ 
so that
$$I_1=\frac{1}{|a_3-a_1|}\iint_{\mathbb{C}}\frac{|d\xi|^2}{|\xi||1-\xi||1-t\xi|}=\frac{I(t,\bar{t})}{|a_3-a_1|},\tag{1}$$
with $\displaystyle t=\frac{a_2-a_1}{a_3-a_1}$. Note that the limit $a_2\rightarrow a_1$ corresponds to $t\rightarrow0$. But we will calculate $I(t,\bar{t})$ in (1) 
explicitly in terms of elliptic integrals without taking any limit, see formula (5) below.



*

*Since we know that our integral is well-defined, instead of integrating over 
$\mathbb{C}$, we can first integrate over $\mathbb{C}\backslash S_{\epsilon}$ (where 
$S$ consists of three disks of radius $\epsilon$ centered at $0$, $1$ and $t^{-1}$), and then consider the limit $\epsilon\rightarrow 0$. Note that our integrand is of the form 
$f(\xi,t)f(\bar{\xi},\bar{t})$ with $$f(\xi,t)=\frac{1}{\sqrt{\xi(1-\xi)(1-t\xi)}}.$$ 

*If only this function $f(\xi,t)$ were simpler, we would be able find its antiderivative $F(\xi,t)$ and then used Stokes theorem to  integrate instead of 2-form $f(\xi,t)f(\bar{\xi},\bar{t})d\xi\wedge d\bar{\xi}$ (over 
$\mathbb{C}\backslash S_{\epsilon}$) the 1-form $F(\bar{\xi},\bar{t})f(\xi,t)d\xi$ over the 
boundary of the $\mathbb{C}\backslash S_{\epsilon}$ consisting of the complex 
contours surrounding the branch cuts. Unfortunately, $f(\xi,t)$ is not that simple 
and we will need the following trick.

*Rewrite the 
integral as
$$I(t,\bar{t})=\frac{1}{2i}\iint f(\xi,t)d\xi\wedge f(\bar{\xi},\bar{t})d\bar{\xi}$$
Then notice that
$$\left\{t(1-t)\frac{\partial^2}{\partial t^2}+
(1-2t)\frac{\partial}{\partial t}-\frac14\right\}f(\xi,t)=
-\frac12\frac{\partial}{\partial\xi}\frac{\xi^2(1-\xi)^2}{\left[
\xi(1-\xi)(1-t\xi)\right]^{3/2}}.$$
Where does this come from? Well, it is a rather well-known fact that if you have an integral of the form
$$\int_{\gamma}\frac{\xi^{a}(1-\xi)^{b}d\xi}{(1-t \xi)^{c}},$$
where $\gamma$ is an arbitrary complex contour starting and ending at one of the three points $0,1,\infty$ (possibly with winding), it solves a hypergeometric equation with parameters simply depending on $a,b,c$ whenever the integral is well-defined. And the standard way to show this is what I used above: one constructs a differential operator whose action on the integrand gives a derivative of something nice. So I only had to find what are the parameters of the hypergeometric equation in your particular case.

*Therefore, applying hypergeometric operator from the left side of the last relation 
to $I(t,\bar{t})$, we can realize the previous program (Stokes theorem) and obtain zero 
(+ similarly for $\bar{t}$).
This implies that
$$I(t,\bar{t})=A\cdot G_a(t)G_a(\bar{t})+B\cdot G_b(t)G_b(\bar{t})+C\cdot
\left[G_a(t)G_b(\bar{t})+G_b(t)G_a(\bar{t})\right],\tag{2}$$
where $G_{a,b}(t)$ are two independent solutions of the hypergeometric equation
$$\left\{t(1-t)\frac{d^2}{d t^2}+
(1-2t)\frac{d}{dt}-\frac14\right\}G(t)=0,\tag{3}$$
and $A,B,C$ are so far unknown constants.

*So now the problem is almost solved. It remains to take two independent solutions of 
(3) and determine the coefficients $A,B,C$ in (2). The hypergeometric functions 
in this case reduce to elliptic integrals and we can choose 
$$G_a(t)= K(t),\qquad G_b(t)=\cdot K(1-t),$$
where $K(t)$ denotes complete elliptic integral of the first kind
$$K(t)=\frac{\pi}{2}{}_2F_1\left(\frac12,\frac12;1;t\right)=\int_0^1\frac{dx}{\sqrt{(1-x^2)(1-tx^2)}}.$$
There are different conventions for the argument of $K(t)$. I use the one used in Mathematica, but
in the notation of the referenced page my $K(t)$ would be $K(\sqrt{t})$.

*Another simple observation is that $I_1$ is symmetric with respect to permutations of $a_{1,2,3}$. But, for example, the exchange $a_1\leftrightarrow a_3$ translates into $t\leftrightarrow 1-t$, and therefore we can conclude that $A=B$. This means that
\begin{align}I\left(t,\bar{t}\right)=A\left[K(t)K\left(\bar{t}\right)+K(1-t)K\left(1-\bar{t}\right)\right]+\\ +C\left[K(t)K\left(1-\bar{t}\right)+K(1-t)K\left(\bar{t}\right)\right].\tag{4}
\end{align}

*My simplest numerical experiments (it would be nice if somebody could confirm them) suggest that $A=0$ and $C=4$, so that the final formula for $I_1$ is given by
$$\boxed{\displaystyle I_1=4\cdot\frac{K(t)K\left(1-\bar{t}\right)+K(1-t)K\left(\bar{t}\right)}{|a_3-a_1|}}\tag{5}$$
with $\displaystyle t=\frac{a_2-a_1}{a_3-a_1}$.
In the end let me  remark that, as $t\rightarrow0$,
$$I(t,\bar{t})=2\pi\left(4\ln2-\ln|t|\right)+o(1).$$

UPD: The comment of achille hui allows to rigorously show that the values of constants $A$ and $C$ are indeed $0$ and $4$, which finishes the derivation of (5).
A: This is not an independent nor complete answer, consider this as some thoughts inspired by O.L's answer.
To integrate the $2^{nd}$ integral $I_{2}$, we ultimately need to evaluate integrals of the form:
$$I(s,t;\bar{s},\bar{t}) = \frac{1}{2i}\int_{\mathbb{C}} \frac{d\bar{z} \wedge dz}{|z(1-z)(1-sz)(1-tz)|} = \frac{1}{2i}\int_{\mathbb{C}} \bar{f}(\bar{z}) d\bar{z} \wedge f(z) dz$$
where 
$$f(z) = \frac{1}{\sqrt{z(1-z)(1-sz)(1-tz)}}$$
Let $\alpha = \sqrt{1-s}, \beta = \sqrt{1-t}, \gamma = \frac{\beta-\alpha}{\beta+\alpha}$ and consider following Möbius transform:
$$\omega = \omega(z) := \frac{(1+\alpha\beta) z - 1}{(\alpha\beta-1) z + 1}
\quad\longleftrightarrow\quad
z = z(\omega) := \frac{\omega+1}{(1-\alpha\beta)\omega + (1+\alpha\beta)}
$$
Under the forward transform $z \mapsto \omega$, the singularities of $f(z)$ get mapped to:
$$\begin{align}
0 & \quad\mapsto\quad \omega(0) = -1\\
1 & \quad\mapsto\quad \omega(1) =  1\\
s^{-1} & \quad\mapsto\quad \omega(s^{-1}) = \omega(\frac{1}{1-\alpha^2}) = \frac{\beta+\alpha}{\beta-\alpha} = \gamma^{-1}\\
t^{-1} & \quad\mapsto\quad \omega(t^{-1}) = \omega(\frac{1}{1-\beta^2}) = \frac{\beta + \alpha}{\alpha-\beta} = -\gamma^{-1}
\end{align}$$
Since $\omega(0) = -\omega(1)$ and $\omega(s^{-1}) = -\omega(t^{-1})$, the integral
will become more symmetric when everything is expressed in terms of $\omega$.
In particular, we have:
$$\begin{align}
dz &= \frac{2\alpha\beta}{((1-\alpha\beta)\omega + (1+\alpha\beta))^2} d\omega\\
z(1-z) &= \frac{\alpha\beta}{((1-\alpha\beta)\omega + (1+\alpha\beta))^2}(1-\omega^2)\\
(1-sz)(1-tz) &= \frac{\alpha\beta}{((1-\alpha\beta)\omega + (1+\alpha\beta))^2}( (\alpha+\beta)^2 - (\beta-\alpha)^2\omega^2)
\end{align}$$
This implies
$$ f(z) dz = \frac{2}{\alpha+\beta}\frac{d\omega}{\sqrt{(1-\omega^2)(1-\gamma^2 \omega^2)}}\tag{*1}$$
Let $\text{sn}(u), \text{cn}(u), \text{dn}(u)$ be the Jacobi elliptic functions corresponds to modulus $\gamma$.
Consider an addition layer of parametrization $\omega = \text{sn}(u)$. The R.H.S of $(*1)$ becomes:
$$\frac{2}{\alpha+\beta}\frac{d\;\text{sn}(u)}{\sqrt{(1-\text{sn}(u)^2)(1-\gamma^2 \text{sn}(u)^2)}} = \frac{2}{\alpha+\beta} \frac{\text{cn}(u)\text{dn}(u)du}{\sqrt{\text{cn}(u)^2 \text{dn}(u)^2}} = \frac{2}{\alpha+\beta} du$$
As a result, we get:
$$
I(s,t;\bar{s}\bar{t}) 
= \frac{4}{|\alpha+\beta|^2}\int_{\mathscr{V}} \frac{d\bar{u} \wedge du}{2i} 
= \frac{4}{|\alpha+\beta|^2} \text{Area}(\mathscr{V})
= \frac{4}{|\alpha+\beta|^2} \frac{\text{Area}(\mathscr{D})}{\mathscr{N}}
$$
where $\mathscr{V}$ is a domain in the parametrization space $u$ such that the map
$\mathscr{V} \ni u \mapsto \omega = \text{sn}(u) \in \mathbb{C}$ covered $\mathbb{C}$ once.
$\mathscr{D}$ is the fundamental domain of $\text{sn}(u)$ and $\mathscr{N}$ is the number of times the map $\mathscr{D} \ni u \mapsto \text{sn}(u)$ covered $\mathbb{C}$.
The fundamental domain $\mathscr{D}$ is spanned by two complex vectors:
$$4 K(\gamma)\quad\text{ and }\quad 2 i K'(\gamma) = 2 i K(\sqrt{1 - \gamma^2})$$
where $K(\gamma) = \int_{0}^{1}\frac{dv}{\sqrt{(1-v^2)(1-\gamma^2 v^2)}}$ is the complete elliptic integral of the $1^{st}$ kind (Please note that I'm using a different convention of K from O.L's answer, this definition is matching the one used in the Wiki's page). We have:
$$\text{Area}(\mathscr{D}) = 8 |K(\gamma) |^2 \left|\Re \frac{K'(\gamma)}{K(\gamma)} \right| = 4 ( K(\gamma)K'(\bar{\gamma}) + K(\bar{\gamma})K'(\gamma) )$$
It is well known the map $u \mapsto \text{sn}(u)$ is a double cover of $\mathbb{C}$. This implies $\mathscr{N} = 2$ and we get:
$$I(s,t;\bar{s},\bar{t}) = \frac{8}{|\alpha+\beta|^2}\left( K(\gamma)K'(\bar{\gamma}) + K(\bar{\gamma})K'(\gamma) \right)\tag{*2}$$
The elliptic integrals of different modulus are related by various transformations. 
In particular, the Landen's transformation allow us to rewrite $(*2)$ into a form more suitable for comparison with O.L's answer. Basically, we have the identities:
$$\begin{align}
K(k)  &= \frac{\pi}{2\text{AGM}(1+k,1-k)} = \frac{\pi}{2\text{AGM}(1,\sqrt{1-k^2})}\\
K'(k) = K(\sqrt{1-k^2}) &= \frac{\pi}{2\text{AMG}(1,k)}
\end{align}$$
where $\text{AGM}(a,b)$ is the arithmetic-geometric mean of $a$ and $b$ which satisfies following functional equation:
$$\text{AGM}(a,b) = \text{AGM}(\frac{a+b}{2},\sqrt{ab})$$
Using this, we have:
$$
\begin{align}
K(\gamma)
=\frac{\pi}{2\text{AGM}(\frac{2\beta}{\beta+\alpha},\frac{2\alpha}{\beta+\alpha})}
&=\left(\frac{\beta+\alpha}{2\beta}\right)\frac{\pi}{2\text{AGM}(1,\sqrt{\frac{1-s}{1-t}})}
=\left(\frac{\beta+\alpha}{2\beta}\right)K(\sqrt{\frac{s-t}{1-t}})\\
K'(\gamma) 
= \frac{\pi}{2\text{AGM}(1,\frac{\beta-\alpha}{\beta+\alpha})}
&= \left(\frac{\beta+\alpha}{\beta}\right)\frac{\pi}{2\text{AGM}(1,\sqrt{\frac{s-t}{1-t}})}
= \left(\frac{\beta+\alpha}{\beta}\right)K'(\sqrt{\frac{s-t}{1-t}})
\end{align}$$
Substitute these into $(*2)$, we get:
$$I(s,t;\bar{s},\bar{t}) = \frac{4}{|1-t|}\left(
K(\sqrt{\frac{s-t}{1-t}}) K'(\sqrt{\frac{\bar{s}-\bar{t}}{1-\bar{t}}}) +
K(\sqrt{\frac{\bar{s}-\bar{t}}{1-\bar{t}}}) K'(\sqrt{\frac{s-t}{1-t}})
\right)$$
If we fix $s$ and send $t$ to $0$, we obtain:
$$\begin{align}\lim_{t\to 0} I(s,t;\bar{s},\bar{t}) 
&= 4\left( K(\sqrt{s})K'(\sqrt{\bar{s}}) + K(\sqrt{\bar{s}})K'(\sqrt{s})\right)\\
&= 4\left( K(\sqrt{s})K(\sqrt{1-\bar{s}})+ K(\sqrt{\bar{s}})K(\sqrt{1-s})\right)\end{align}
$$
Reproducing the result for $I(s,\bar{s})$ in O.L's answer (up to convention used for K).
