(* Code for calculating the probability distribution function of the least of the nine acute angles between pairs of edges of two randomly disoriented cubes, its mean, and standard deviation, based on J. K. Mackenzie, Second Paper on Statistics Associated with the Random Disorientation of Cubes, Biometrika, Vol. 45, No. 1-2 (1958), pp. 229-240. *) SetOptions[NIntegrate, WorkingPrecision -> 25, Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 5000}]; gamma[a_] := ArcSin[Sqrt[((Sqrt[2] + 1)^2 * Tan[a/2]^2 - 1)/(4 * Sqrt[2] * Tan[a/2]^2)]]; mean[k_] := (288/Pi^2) * (Pi^2/32 * NIntegrate[a^k * Sin[a], {a, 0, ArcCos[2/3]}] - NIntegrate[a^k * Sin[a] * ArcSin[Tan[a/2] * Cos[x]], {a, 0, ArcCos[2/3]}, {x, 0, Pi/4}] - (Pi/4) * NIntegrate[a^k * Sin[a] * gamma[a], {a, Pi/4, ArcCos[2/3]}] + NIntegrate[a^k * Sin[a] * ArcSin[Tan[a/2] * Cos[x]], {a, Pi/4, ArcCos[2/3]}, {x, Pi/8 - gamma[a], Pi/8 + gamma[a]}]); RealDigits[mean[1], 10, 18][[1]] RealDigits[Sqrt[mean[2] - mean[1]^2], 10, 18][[1]]