# Find closed form for quadruple integral

I am trying to find a closed form of the following integral $$\int _0^{\infty }\int _0^x\int _0^y\int _0^z \exp \left( -\frac{a x^2}{2}-\frac{b y^2}{2}-\frac{c z^2}{2}-\frac{d w^2}{2} \right) \,\mathrm{d}w\,\mathrm{d}z\,\mathrm{d}y\,\mathrm{d}x$$ where $$a,b,c,d>0$$ are some constants.

My idea is to change variable to to polar system by letting $$x=\frac{r \cos (\alpha ) \cos (\beta ) \cos (\theta )}{\sqrt{a}}, \quad y=\frac{r \cos (\alpha ) \cos (\beta ) \sin (\theta )}{\sqrt{a}}$$ $$z=\frac{r \sin (\alpha ) \cos (\beta )}{\sqrt{c}}, \quad w=\frac{r \sin (\beta )}{\sqrt{d}}$$

This reduces the original integral into $$\int_0^{\tan ^{-1}\left(\sqrt{\frac{b}{a}}\right)} \frac{\sin (\theta ) \tan ^{-1}\left(\sin (\theta ) \sqrt{\frac{d}{b+c \sin ^2(\theta )}}\right)}{\sqrt{a b d \left(b+c \sin ^2(\theta )\right)}} \, d\theta$$ But then I get stuck here.

PS: I am interested in this because I found that $$\int _0^{\infty }\int _0^x \exp \left(-\frac{a x^2}{2}-\frac{b y^2}{2}\right)\,\mathrm{d}y\,\mathrm{d}x = \frac{\tan ^{-1}\left(\sqrt{\frac{b}{a}}\right)}{\sqrt{a b}}$$ and $$\int _0^{\infty }\int _0^x\int _0^y \exp \left( -\frac{a x^2}{2}-\frac{b y^2}{2}-\frac{c z^2}{2} \right) \,\mathrm{d}z\,\mathrm{d}y\,\mathrm{d}x = \frac{\sqrt{\pi/2 } }{\sqrt{a b c}} \left(\tan ^{-1}\left(\sqrt{\frac{c}{b}}\right)-\tan ^{-1}\left(\sqrt{\frac{a c}{b (a+b+c)}}\right)\right)$$ So I am trying to generalize this. Maybe this is already known?

• At least for all constants equal to $1$ Wolfram Alpha claims that the answer is $\tfrac {\pi^2}{96}$. – Alex Ravsky Apr 28 at 15:29
• Well from the pattern I deduce that there must be a $\sf{\dfrac{\tan^{-1}\left(\sqrt{\dfrac dc}\right)}{\sqrt{abcd}}}$ term... – TheSimpliFire Apr 29 at 6:26

Define: $$\begin{eqnarray} I(a,b,c,d):=\int_0^{\tan ^{-1}\left(\sqrt{\frac{b}{a}}\right)} \frac{\sin (\theta ) \tan ^{-1}\left(\sin (\theta ) \sqrt{\frac{d}{b+c \sin ^2(\theta )}}\right)}{\sqrt{a b d \left(b+c \sin ^2(\theta )\right)}} \, d\theta \quad (i) \end{eqnarray}$$ Then also define:

$$\begin{eqnarray} {\mathfrak F}^{(A,B)}_{a,b} &:=& \int\limits_A^B \frac{\log(z+a)}{z+b} dz\\ &=& F[B,a,b] - F[A,a,b] + 1_{t^* \in (0,1)} \left( -F[A+(t^*+\epsilon)(B-A),a,b] + F[A+(t^*-\epsilon)(B-A),a,b] \right) \quad (ii) \end{eqnarray}$$ where $$\begin{eqnarray} t^*:=-\frac{Im[(A+b)(b^*-a^*)]}{Im[(B-A)(b^*-a^*)]} \end{eqnarray}$$ and $$$$F[z,a,b] := \log(z+a) \log\left( \frac{z+b}{b-a}\right) + Li_2\left( \frac{z+a}{a-b}\right)$$$$ for $$a$$,$$b$$,$$A$$,$$B$$ being complex.

Then we have: $$\begin{eqnarray} I(a,b,c,d)&=& \frac{1}{\sqrt{a}} \int\limits_0^{\sqrt{\frac{d}{a+b+c}}} \frac{u \tan ^{-1}(u)}{\left(d-c u^2\right) \sqrt{d-u^2 (b+c)}} du\\ &=& \frac{\sqrt{d}}{\sqrt{a} (b+c)}\int\limits_0^{\sin ^{-1}\left(\sqrt{\frac{b+c}{a+b+c}}\right)} \frac{\sin (\phi ) \tan ^{-1}\left(\sin (\phi ) \sqrt{\frac{d}{b+c}}\right)}{d-\frac{c d \sin ^2(\phi )}{b+c}} d\phi\\ &=&-\frac{2 i}{\sqrt{a} \sqrt{d}} \int\limits_0^{\frac{\sqrt{\frac{b+c}{a+b+c}}}{\sqrt{\frac{a}{a+b+c}}+1}} \frac{t}{b \left(t^2+1\right)^2+c \left(t^2-1\right)^2} \log \left(\frac{2 i t \sqrt{\frac{d}{b+c}}+t^2+1}{-2 i t \sqrt{\frac{d}{b+c}}+t^2+1}\right) dt\\ &=&\frac{1}{4} \frac{1}{\sqrt{a b c d}} \sum\limits_{\xi=1}^4 \sum\limits_{\eta=1}^4 (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\left\lfloor \frac{\xi -1}{2}\right\rfloor } \int\limits_0^{\frac{\sqrt{\frac{b+c}{a+b+c}}}{\sqrt{\frac{a}{a+b+c}}+1}} \frac{\log \left(i (-1)^{\left\lfloor \frac{\xi -1}{2}\right\rfloor } \sqrt{\frac{d}{b+c}}+i (-1)^{\xi -1} \sqrt{\frac{b+c+d}{b+c}}+t\right)}{t-i (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\eta +1} e^{i (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor } \tan ^{-1}\left(\frac{\sqrt{c}}{\sqrt{b}}\right)}} dt \\ &=&\frac{1}{4} \frac{1}{\sqrt{a b c d}} \sum\limits_{\xi=1}^4 \sum\limits_{\eta=1}^4 (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\left\lfloor \frac{\xi -1}{2}\right\rfloor } % {\mathfrak F}^{(0,\frac{\sqrt{\frac{b+c}{a+b+c}}}{\sqrt{\frac{a}{a+b+c}}+1})}_{i (-1)^{\left\lfloor \frac{\xi -1}{2}\right\rfloor } \sqrt{\frac{d}{b+c}}+i (-1)^{\xi -1} \sqrt{\frac{b+c+d}{b+c}},-i (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\eta +1} e^{i (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor } \tan ^{-1}\left(\frac{\sqrt{c}}{\sqrt{b}}\right)}} \end{eqnarray}$$

In the top line we substituted for $$u=\sin(\theta) \sqrt{d/(b+c \sin(\theta)^2)}$$. In the second line we substituted $$u = \sqrt{d/(c+b)} \sin(\phi)$$. In the third line we substituted $$t=\tan(\phi/2)$$. In the forth line we used partial fraction decomposition and properties of the logarithm. Finally in the fifth line we used the anti-derivative defined in $$(ii)$$.

Clear[F]; Clear[FF];
F[z_, a_, b_] :=
Log[a + z] Log[(b + z)/(-a + b)] + PolyLog[2, (a + z)/(a - b)];
FF[A_, B_, a_, b_] :=
Module[{result, ts, zs, zsp, zsm, eps = 10^(-15)},
(*This is Integrate[Log[z+a]/(z+b),{z,A,B}] where all a,b,A,
and B are complex. *)
result = F[B, a, b] - F[A, a, b];

ts = - (Im[(A + b) (Conjugate[b] - Conjugate[a])]/
Im[(B - A) (Conjugate[b] - Conjugate[a])]);
If[0 <= ts <= 1,
zsp = A + (ts + eps) (B - A);
zsm = A + (ts - eps) (B - A);
result += -F[zsp, a, b] + F[zsm, a, b];
];

result
];

{a, b, c, d} = RandomReal[{0, 3}, 4, WorkingPrecision -> 50];
NIntegrate[
Exp[-a/2 x^2 - b/2 y^2 - c/2 z^2 - d/2 w^2], {x, 0, Infinity}, {y, 0,
x}, {z, 0, y}, {w, 0, z}]
NIntegrate[
Sin[th]/Sqrt[a b d (b + c Sin[th]^2)] ArcTan[
Sin[th] Sqrt[d/(b + c Sin[th]^2)]], {th, 0, ArcTan[Sqrt[b/a]]}]
1/Sqrt[a ] NIntegrate[
u ArcTan[u] 1/((d - c u^2) Sqrt[d - (c + b) u^2]), {u, 0, Sqrt[
d/ (a + b + c)]}]
Sqrt[d]/(Sqrt[a ] (b + c))
NIntegrate[
Sin[phi] ArcTan[
Sqrt[d/(c + b)] Sin[phi]] 1/(d - c (d/(c + b) Sin[phi]^2)) , {phi,
0, ArcSin[ Sqrt[( c + b)/ (a + b + c)]]}]
- I 2/(Sqrt[a ] Sqrt[d])
NIntegrate[
t /(c (-1 + t^2)^2 + b (1 + t^2)^2) Log[(
1 + t^2 + 2 I Sqrt[d/(b + c)] t)/(
1 + t^2 - 2 I Sqrt[d/(b + c)] t)], {t, 0, Sqrt[(b + c)/(
a + b + c)]/(1 + Sqrt[a/(a + b + c)])}]
- I 2/(Sqrt[a ] Sqrt[d])
NIntegrate[
t /(c (-1 + t^2)^2 +
b (1 + t^2)^2) Log[((1/
2 (2 I Sqrt[d/(b + c)] - Sqrt[-4 - (4 d)/(b + c)]) +
t) (1/2 (2 I Sqrt[d/(b + c)] + Sqrt[-4 - (4 d)/(b + c)]) +
t))/((1/2 (-2 I Sqrt[d/(b + c)] - Sqrt[-4 - (4 d)/(b + c)]) +
t) (1/2 (-2 I Sqrt[d/(b + c)] + Sqrt[-4 - (4 d)/(b + c)]) +
t))], {t, 0, Sqrt[(b + c)/(a + b + c)]/(
1 + Sqrt[a/(a + b + c)])}]
1/Sqrt[a b c d]  1/4 NIntegrate[
Sum[(-1)^(Floor[(eta - 1)/2]) (-1)^
Floor[(xi - 1)/2] Log[
t + (-1)^Floor[(xi - 1)/2] I Sqrt[d/(b + c)] + (-1)^(xi - 1)
I Sqrt[( b + c + d)/(b + c)]]/(
t - (-1)^(1 + eta +
Floor[(eta - 1)/2]) I Exp[(-1)^(Floor[(eta - 1)/2]) I ArcTan[
Sqrt[c]/Sqrt[b]]]), {xi, 1, 4}, {eta, 1, 4}], {t, 0, Sqrt[(
b + c)/(a + b + c)]/(1 + Sqrt[a/(a + b + c)])}]
1/Sqrt[a b c d]  1/4 Sum[(-1)^(Floor[(eta - 1)/2]) (-1)^
Floor[(xi - 1)/2] FF[0, Sqrt[(b + c)/(a + b + c)]/(
1 + Sqrt[a/(
a + b + c)]), (-1)^Floor[(xi - 1)/2] I Sqrt[d/(b + c)] + (-1)^(
xi - 1) I Sqrt[( b + c + d)/(
b + c)], -(-1)^(1 + eta +
Floor[(eta - 1)/2]) I Exp[(-1)^(Floor[(eta - 1)/2]) I ArcTan[
Sqrt[c]/Sqrt[b]]]], {xi, 1, 4}, {eta, 1, 4}]


Update: As a sanity check look at the case $$a=b=c=d=1$$. Then define: $$\begin{eqnarray} M1&:=&\left( \begin{array}{cccc} -1+\sqrt{3} & \sqrt{2} & \sqrt{2} & -1+\sqrt{3} \\ 1 & \sqrt{2-\sqrt{3}} & \sqrt{2-\sqrt{3}} & 1 \\ \sqrt{2-\sqrt{3}} & 1 & 1 & \sqrt{2-\sqrt{3}} \\ \sqrt{2} & -1+\sqrt{3} & -1+\sqrt{3} & \sqrt{2} \\ \end{array} \right)\\ M2&:=&\left( \begin{array}{cccc} \frac{1}{\sqrt{2}} & \frac{1}{2} \left(1+\sqrt{3}\right) & \frac{1}{2} \left(1+\sqrt{3}\right) & \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & \frac{1}{2} \left(-1+\sqrt{3}\right) & \frac{1}{2} \left(-1+\sqrt{3}\right) & \frac{1}{\sqrt{2}} \\ \frac{1}{2} \left(-1+\sqrt{3}\right) & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & \frac{1}{2} \left(-1+\sqrt{3}\right) \\ \frac{1}{2} \left(1+\sqrt{3}\right) & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & \frac{1}{2} \left(1+\sqrt{3}\right) \\ \end{array} \right)\\ A1&:=&\left( \begin{array}{cccc} -\frac{\pi }{6} & \frac{\pi }{12} & -\frac{\pi }{4} & 0 \\ \frac{5 \pi }{6} & \frac{\pi }{12} & \frac{5 \pi }{12} & -\frac{\pi }{3} \\ -\frac{5 \pi }{12} & \frac{\pi }{3} & -\frac{5 \pi }{6} & -\frac{\pi }{12} \\ \frac{\pi }{4} & 0 & \frac{\pi }{6} & -\frac{\pi }{12} \\ \end{array} \right)\\ A2&:=&\left( \begin{array}{cccc} -\frac{\pi }{12} & \frac{\pi }{6} & -\frac{\pi }{6} & \frac{\pi }{12} \\ \frac{7 \pi }{12} & -\frac{\pi }{6} & \frac{\pi }{6} & -\frac{7 \pi }{12} \\ -\frac{\pi }{6} & \frac{7 \pi }{12} & -\frac{7 \pi }{12} & \frac{\pi }{6} \\ \frac{\pi }{6} & -\frac{\pi }{12} & \frac{\pi }{12} & -\frac{\pi }{6} \\ \end{array} \right) \end{eqnarray}$$ and we have $$\begin{eqnarray} I(1,1,1,1)=\frac{1}{4} \sum\limits_{\xi=1}^4 \sum\limits_{\eta=1}^4 (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\left\lfloor \frac{\xi -1}{2}\right\rfloor } \left( Li_2(M1_{\xi,\eta}\exp(\imath A1_{\xi,\eta}))- Li_2(M2_{\xi,\eta}\exp(\imath A2_{\xi,\eta})) \right) \end{eqnarray}$$ We have checked numerically that this quantity above coincides with $$\pi^2/96$$ to one hundred digits. It would be interesting to prove this analytically.

• This is an impressive amount of algebra, but looking at Your Mathematica code at the end we're still left with one integral to be computed. It's unclear to me how is this more efficient than the original integral. What were we going for here? – Radost Apr 30 at 16:03
• @Radost Sorry, the code was outdated. Now the code uses no numerical integrals. You can easily get the result with fifty digits of accuracy as I do above. But going back to your point I have played around with integrals like that many times. They simply reduce to poly-logarithms. there is no way around that. Yet polylogarithms can be computed with arbitrary precision whereas your original integral cannot. – Przemo Apr 30 at 16:22

This integral is up to normalization constant integral of multivariate gaussian distribution. Due too lack of cross terms we get that there are 4 independent zero mean gaussian variables involved.

We can begin by rewriting this a probability of an event.

Let $$X,Y,Z,W$$ independent normally distributed with zero mean and possibly different variances.

Then integral reduces to: $$\mathbb{P} (0 < W < Z < Y < X)$$

Which in turn is the same as:

$$\mathbb{P}(W > 0 \land Z-W >0 \land Y-Z >0 \land X-Y >0)$$

and one can observe that $$W,Z-W,Y-Z,X-Y$$ are correlated jointly normal variables. Probability that all components of a jointly normal vector are positive is called orthant probability and in general doesn't have closed form expression. I suppose this is quite special case with covariance matrix almost diagonal so maybe there are some articles about how to approach this special case.

For the case with 3 or less variables formulas are known (cf this question for example) and I guess they would coincide with what you've found.

Here we give an answer using a different method. Assume that $$a\ge0$$, $$b\ge 0$$, $$c\ge 0$$ and $$d\ge 0$$. Define: $$$$I(a,b,c,d):=\int _0^{\infty }\int _0^x\int _0^y\int _0^z \exp \left( -\frac{a x^2}{2}-\frac{b y^2}{2}-\frac{c z^2}{2}-\frac{d w^2}{2} \right) \,\mathrm{d}w\,\mathrm{d}z\,\mathrm{d}y\,\mathrm{d}x$$$$ Then we have: $$\begin{eqnarray} &&I(a,b,c,d)=\\ &&\sqrt{\frac{\pi}{2 c d}} \int\limits_0^\infty \int\limits_0^x \int\limits_0^{\sqrt{c} y}\exp(-\frac{a x^2}{2} - \frac{b y^2}{2}-\frac{ z^2}{2})\cdot \text{erf}(\sqrt{\frac{d}{2 c}} z)\; dz \; dy \;dx=\\ &&\frac{2 \pi}{\sqrt{c d}} \int\limits_0^\infty \int\limits_0^x \exp(-\frac{a x^2}{2}-\frac{b y^2}{2}) \left(\frac{1}{2\pi} \arctan[\sqrt{\frac{d}{c}}] - T(\sqrt{c} y,\sqrt{\frac{d}{c}}\right)\; dy \; dx=\\ && \frac{2 \pi^2}{\sqrt{a b c d}}\int\limits_0^\infty \frac{\exp(-\frac{y^2}{2})}{\sqrt{2\pi}} \text{erfc}(\sqrt{\frac{a}{2 b}y})\left(\frac{1}{2\pi} \arctan[\sqrt{\frac{d}{c}}] - T(\sqrt{\frac{c}{b}} y,\sqrt{\frac{d}{c}}\right) \; dy=\\ && \frac{2 \pi^2}{\sqrt{a b c d}} \left( \frac{\left(\pi -2 \tan ^{-1}\left(\sqrt{\frac{a}{b}}\right)\right) \tan ^{-1}\left(\sqrt{\frac{d}{c}}\right)-\pi \sin ^{-1}\left(\frac{\sqrt{b d}}{\sqrt{(b+c) (c+d)}}\right)}{4 \pi ^2}+ \int\limits_0^\infty \frac{\exp(-\frac{y^2}{2})}{\sqrt{2 \pi}} \text{erf}(\sqrt{\frac{a}{2 b}} y) T(\sqrt{\frac{c}{b}} y,\sqrt{\frac{d}{c}} )\; dy \right)=\\ && \frac{2 \pi^2}{\sqrt{a b c d}} \left( \frac{\tan ^{-1}\left(\sqrt{\frac{d}{c}}\right)-\tan ^{-1}\left(\frac{\sqrt{b d}}{\sqrt{c (b+c+d)}}\right)}{4 \pi } +\right.\\ &&\left. \frac{1}{8 \pi^2} \sum\limits_{i=1}^4 \sum\limits_{j=1}^4 (-1)^{\left\lfloor \frac{i-1}{2}\right\rfloor +j} {\mathfrak F}^{(1,\sqrt{\frac{a+b+c}{b+c}}-\sqrt{\frac{a}{b+c}})}_{i \left((-1)^{\left\lfloor \frac{j-1}{2}\right\rfloor } \sqrt{\frac{b+c+d}{b+c}}+(-1)^j \sqrt{\frac{d}{b+c}}\right),(-1)^{i+1} \sqrt{\frac{c}{b+c}}+i \sqrt{\frac{b}{b+c}} (-1)^{\left\lceil \frac{i-1}{2}\right\rceil +1}} \right) \end{eqnarray}$$ In the first line we integrated over $$w$$ using the definition of the error function. In the second line we integrated over $$z$$ using the definition of the Owen's T function. In the third line we swapped the order of integration and integrated over $$x$$ using the definition of the complementary error function. In the forth line we split the integral into doable integrals and a more complicated one and finally in the fifth line we evaluated the remaining integral using An integral involving a Gaussian, error functions and the Owen's T function. .

Clear[F]; Clear[FF];
F[z_, a_, b_] :=
Log[a + z] Log[(b + z)/(-a + b)] + PolyLog[2, (a + z)/(a - b)];
FF[A_, B_, a_, b_] :=
Module[{result, ts, zs, zsp, zsm, eps = 10^(-50)},
(*This is Integrate[Log[z+a]/(z+b),{z,A,B}] where all a,b,A,
and B are complex. *)
result = F[B, a, b] - F[A, a, b];

ts = - (Im[(A + b) (Conjugate[b] - Conjugate[a])]/
Im[(B - A) (Conjugate[b] - Conjugate[a])]);
If[0 <= ts <= 1,
zsp = A + (ts + eps) (B - A);
zsm = A + (ts - eps) (B - A);
result += -F[zsp, a, b] + F[zsm, a, b];
];

result
];
J[a_, b_, c_] :=
1/ Pi^2 (ArcTan[Sqrt[2] a]/2 ArcTan[ c] +
1/8  Sum[
FF[1, ( Sqrt[1 + 2 a^2 + b^2] - Sqrt[2] a)/Sqrt[
1 + b^2], ((-1)^j I b c + (-1)^Floor[(j - 1)/2] I Sqrt[
1 + b^2 + b^2 c^2])/Sqrt[
1 + b^2], -(((-1)^Ceiling[(i - 1)/2] I + (-1)^i b)/Sqrt[
1 + b^2])] (-1)^(j + Floor[(i - 1)/2]), {i, 1, 4}, {j, 1,
4}] );

{a, b, c, d} = RandomReal[{0, 10}, 4, WorkingPrecision -> 50];
NIntegrate[
Exp[-a/2 x^2 - b/2 y^2 - c/2 z^2 - d/2 w^2], {x, 0, Infinity}, {y, 0,
x}, {z, 0, y}, {w, 0, z}]
Sqrt[\[Pi]/2]/(Sqrt[c] Sqrt[d])
NIntegrate[
Exp[-a/2 x^2 - b/2 y^2 - 1/2 z^2] Erf[Sqrt[d/(2 c)] z], {x, 0,
Infinity}, {y, 0, x}, {z, 0, Sqrt[c] y}]
(2 \[Pi])/Sqrt[c d]
NIntegrate[
Exp[-a/2 x^2 - b/2 y^2] (ArcTan[Sqrt[d/c]]/(2 \[Pi]) -
OwenT[Sqrt[c] y, Sqrt[d/c]]), {x, 0, Infinity}, {y, 0, x}]
(2  \[Pi]^2)/Sqrt[a b c d]
NIntegrate[
Exp[ -1/2 y^2]/Sqrt[2 Pi]
Erfc[Sqrt[a/(2 b)] y] (ArcTan[Sqrt[d/c]]/(2 \[Pi]) -
OwenT[Sqrt[c/b] y, Sqrt[d/c]]), {y, 0, Infinity}]
(2  \[Pi]^2)/Sqrt[
a b c d] ((-\[Pi] ArcSin[Sqrt[b d]/
Sqrt[(b + c) (c + d)]] + (\[Pi] - 2 ArcTan[Sqrt[a/b]]) ArcTan[
Sqrt[d/c]])/(4 \[Pi]^2) +
NIntegrate[
Exp[ -1/2 y^2]/Sqrt[2 Pi]
Erf[Sqrt[a/(2 b)] y] OwenT[Sqrt[c/b] y, Sqrt[d/c]], {y, 0,
Infinity}])
(2  \[Pi]^2)/Sqrt[
a b c d] ((-ArcTan[Sqrt[b d]/ Sqrt[c (b + c + d)]] +
ArcTan[Sqrt[d/c]])/(4 \[Pi]) +
1/( 8 Pi^2)
Sum[FF[1, Sqrt[(a + b + c)/(b + c)] - Sqrt[a/(b + c)],
I  ((-1)^j Sqrt[d/(b + c)] + (-1)^Floor[1/2 (-1 + j)] Sqrt[(
b + c + d)/(b + c)]), (-1)^(i + 1) Sqrt[c/(b + c)] +
Sqrt[b/(b + c)] I (-1)^(1 + Ceiling[1/2 (-1 + i)])] (-1)^(
j + Floor[(i - 1)/2]), {i, 1, 4}, {j, 1, 4}] )