How to compute this triple integral?

Questions
1

Who know how to count

Integrate[
 Integrate[
  Integrate[
   Sin[x] Cos[y] Cos[z]/(1 + Cos[x] Cos[y] Cos[z]), {x, 0, Pi/2}], {y,
    0, Pi/2}], {z, 0, Pi/2}]
2

This integral is little tricky to handle directly. But with some manual help (double angle formula + Feynman integration technique), Mathematica can carry out the computation symbolically.

Numeric check

In[]:= NIntegrate[Sin[x] Cos[y] Cos[z]/(1 + Cos[x] Cos[y] Cos[z]), {x, 0, Pi/2}, {y, 0,Pi/2}, {z, 0, Pi/2}]
Out = 0.78665

Integrate about x and let g(z) = Log[1+Cos[y] Cos[z]] 

In[]:= Integrate[Sin[x] Cos[y] Cos[z]/(1+Cos[x] Cos[y] Cos[z]),{x,0,Pi/2},Assumptions->0<y<Pi/2&&0<z<Pi/2]
Out= Log[1+Cos[y] Cos[z]]

Use Feynman’s differentiation to take derivative of g[z] about z and then integrate about y. Call it g’[z]:

In[]:= D[Log[1 + Cos[y] Cos[z]], z]
Out= -((Cos[y] Sin[z])/(1 + Cos[y] Cos[z]))

In[]:= Integrate[%,{y,0,\[Pi]/2},Assumptions->0<z<\[Pi]/2]
Out= -(1/2) (\[Pi]+4 ArcTan[Cot[z]-Csc[z]] Csc[z]) Tan[z]

Then use integral by parts to bring it back: \int g(z) dz= z \cdot g - \int z \cdot g' , you can check that the first term is 0 at z=0 and z = \pi /2:

In[]:= -Integrate[%*z,{z,0,Pi/2}]
Out= -1/2 (Pi^2 Log[2]-7 Zeta[3])

The numeric value is

In[]:= N[%]
Out= 0.78665

I think the place where Wolfram can improve is to automatically further simplify this part: ArcTan[Cot[z]-Csc[z]] . If we replace z with z=2u, we have:

In[]:= FullSimplify[ArcTan[Cot[2u]-Csc[2u]],Assumptions->0<u<Pi/4]
Out= -u

This step is correct because cos(2 \cdot u) = 1 - 2 \cdot sin(u)^2 and sin (2 \cdot u) = 2 \cdot sin (u) \cdot cos (u) . Then the expression inside ArcTan is simply -Tan[u]. Recall ArcTan is an odd function so this compound part with ArcTan[-Tan[u]] reduces to -u = -z/2.

3

Thank You very much! will be possible in the future to do general alghorhitm to count multi integrals and adopt them in Wolfram programmes (Mathematica, Cloud,Alpha)?

4

Here is another way, doing iterated integration with a formula simplification half way through:

f[x_, y_, z_] := Sin[x] Cos[y] Cos[z]/(1 + Cos[x] Cos[y] Cos[z]);
fx[y_, z_] = Integrate[f[x, y, z], {x, 0, Pi/2},
  Assumptions -> 0 < y < Pi/2 && 0 < z < Pi/2]
fxy[z_] = Integrate[fx[y, z], {y, 0, Pi/2},
  Assumptions -> 0 < z < Pi/2]
simpRule = Arg[pt_] :> FullSimplify[VectorAngle[
     FullSimplify[ReIm[pt], 0 < z < Pi/2], {1, 0}],
    0 < z < Pi/2];
fxySimp[z_] = FullSimplify[fxy[z] /. simpRule,
  0 < z < Pi/2]
fxyz[z_] = Integrate[fxySimp[z], z]
Limit[fxyz[z], z -> Pi/2,
   Direction -> "FromBelow"] - fxyz[0] //
 FullSimplify
% // N
5

How to use this procedure when end of integration is infinity e.g.
Integrate[
1/Sqrt[(4 + x^2) (4 + y^2) (4 + (x + y)^2 (4 + (x - y)^2))], {x, 0,
Infinity}, {y, 0, Infinity}]