A known closed form for Borchardt mean (generalization of AGM) - why doesn't it work? There is a curious four parameter iteration introduced by Borchardt:
$$a_{n+1}=\frac{a_n+b_n+c_n+d_n}{4} \\ b_{n+1}=\frac{\sqrt{a_n b_n}+\sqrt{c_n d_n}}{2} \\ c_{n+1}=\frac{\sqrt{a_n c_n}+\sqrt{b_n d_n}}{2} \\ d_{n+1}=\frac{\sqrt{a_n d_n}+\sqrt{b_n c_n}}{2} $$
Apparently, the limit of this iteration, denoted hereafter $B(a_0,b_0,c_0,d_0)$ has a closed form in terms of a certain double integral. However, this 'closed form' doesn't check out when I try implementing it in Mathematica.
For a complete description see this paper, section 2.5. For more general information see this paper.
Here is the description of the closed form from the first link above. First, define:
$$A = a_0 + b_0 + c_0 + d_0 \\
B = a_0 + b_0 - c_0 - d_0 \\
C = a_0 - b_0 + c_0 - d_0 \\
D = a_0 - b_0 - c_0 + d_0$$
(Note the problem - some of the numbers above can be negative).
Then define:
$$B_1=\frac{\sqrt{A B}+\sqrt{C D}}{2} \\ B_2=\frac{\sqrt{A B}-\sqrt{C D}}{2} \\ C_1=\frac{\sqrt{A C}+\sqrt{B D}}{2} \\ C_2=\frac{\sqrt{A C}-\sqrt{B D}}{2} \\ D_1=\frac{\sqrt{A D}+\sqrt{B C}}{2} \\ D_2=\frac{\sqrt{A D}-\sqrt{B C}}{2}
$$
(Note the problem - some of the numbers above can be complex).
Then define:
$$\Delta=\sqrt[4]{ABCDB_1C_1D_1B_2C_2D_2}$$
And finally, the main parameters:
$$\alpha_0=\frac{A C B_1}{\Delta} \\ \alpha_1=\frac{C C_1 D_1}{\Delta} \\ \alpha_2=\frac{A C_2 D_1}{\Delta} \\ \alpha_3=\frac{B_1 C_1 C_2}{\Delta}$$
And define a function:
$$R(x):=x(x-\alpha_0)(x-\alpha_1)(x-\alpha_2)(x-\alpha_3)$$
According to the paper, the limit of the iteration above is equal to:
$$\frac{\pi^2}{B(a_0,b_0,c_0,d_0)}=\int _0^{\alpha_3}\int _{\alpha_2}^{\alpha_1}\frac{x-y}{\sqrt{R(x) R(y)}}dxdy \tag{1}$$
I tried to implement this 'closed form' in Mathematica. However, for most initial conditions, even simple ones, Mathematica has trouble computing the integral and gives complex values anyway.

How can I tell if $(1)$ is correct? Is there a typo somewhere in my formulas? Or maybe there is some typo in the linked paper?

I want to know if $(1)$ is correct, and if not what is the correct form. I don't need the full proof. I tried to find Borchardt's original paper, but I couldn't. And it's not in English anyway.
 A: Here are just a few remarks.
Borchart's original article (referred to as Berl. Monatsber. , pages 611–621) is fully available (in German) at  http://bibliothek.bbaw.de/bibliothek-digital/digitalequellen/schriften/anzeige/index_html?band=09-mon/1876&seite:int=647.
However this article itself uses work already done in a longer older article (where most of the hard compatations are done), (referred to as (Mem. Acad. Berlin, 1878, pp. 33-96), and also fully available online at http://bibliothek.bbaw.de/bbaw/bibliothek-digital/digitalequellen/schriften/anzeige/index_html?band=07-abh/1878&seite:int=93.
Your square-root of negative number problem can be solved very simply by starting from $1$ (or any nonzero index) instead of $0$.  Because of Borchardt's identity 
$$
a_{n+1}+\varepsilon b_{n+1}+ \varepsilon' c_{n+1}+\varepsilon \varepsilon'd_{n+1}=\frac{(\sqrt{a_n}+\varepsilon \sqrt{b_n}+ \varepsilon'\sqrt{c_n}+\varepsilon \varepsilon'\sqrt{d_n})^2}{4} \tag{1}
$$
for any two signs $\varepsilon,\varepsilon'$ we have for $n\geq 1$,
$$a_{n+1}+\varepsilon b_{n+1}+ \varepsilon ' c_{n+1}+\varepsilon \varepsilon 'd_{n+1} \geq 0 \tag{2}$$ 
So if you use any $a_1,b_1,c_1,d_1$ instead of $a_0,b_0,c_0,d_0$ to define the $B,C,\alpha$ constants, you will get only real and positive numbers.
