# How to exactly calculate the volume?

Let us consider the intersection of four cylinders of the unit radius along the big diagonals of the cube $$[-10,10]^3$$ and the cylinder of the unit radius with the $$z$$-axis as its axis. More exactly,

reg = ImplicitRegion[(2 x/3 + y/3 + z/3)^2 + (2 y/3 + x/3 - z/3)^2 + (2 z/3 + x/3 - y/3)^2 <= 1 &&
(2 x/3 - y/3 - z/3)^2 + (2 y/3 - x/3 - z/3)^2 + (2 z/3 - x/3 - y/3)^2 <=  1 &&
(2 x/3 + y/3 - z/3)^2 + (2 y/3 + x/3 + z/3)^2 + (2 z/3 - x/3 + y/3)^2 <= 1 &&
(2 x/3 - y/3 + z/3)^2 + (2 y/3 - x/3 + z/3)^2 + (2 z/3 +  x/3 + y/3)^2 <= 1 && x^2 + y^2 <= 1, {x, y, z}];

RegionPlot3D[(2 x/3 + y/3 + z/3)^2 + (2 y/3 + x/3 - z/3)^2 + (2 z/3 +
x/3 - y/3)^2 <=   1 && (2 x/3 - y/3 - z/3)^2 + (2 y/3 - x/3 - z/3)^2 + (2 z/3 - x/3 - y/3)^2 <=  1 &&
(2 x/3 + y/3 - z/3)^2 + (2 y/3 + x/3 + z/3)^2 + (2 z/3 - x/3 +
y/3)^2 <=  1 && (2 x/3 - y/3 + z/3)^2 + (2 y/3 - x/3 + z/3)^2 + (2 z/3 + x/3 +
y/3)^2 <= 1 && x^2 + y^2 <= 1, {x, -3/2, 3/2}, {y, -3/2, 3/2}, {z, -3/2, 3/2}, PlotPoints -> 50]


Its volume

Volume[reg]


4.40045

It is very probable that the exact result equals 22/5. Just for sportive interest, how can I prove or disprove it with Mathematica? I don't find the answer here.

• I would be really surprised if the true value were a rational number. Yhinking of the Steinmetz solid, a rational times $\pi$ is more likely. Feb 9 at 7:28
• One possibility to make this problem easier for Mathematica is to attempt splitting the solid into a polyhedron and "caps" on top of it and calculating these volumes separately. Feb 9 at 11:58
• @HenrikSchumacher Surprising fact about Steinmetz solids is that their volume does not contain $\pi$. This was already known to Archimedes. The same is true for the OP (see my answer). Feb 9 at 21:42
• @yarchik Ha! So I got the Steinmetz solid wrong! You see me really surprised! =D Feb 10 at 6:41
• @HenrikSchumacher Although the hypothesis in this case didn't hold, it wasn't entirely unfounded - after all, the original Steinmetz solid (intersection of two cylinders at right angles) doesn't just have an algebraic volume, but it's actually rational ($\frac{16}{3}$)! Feb 11 at 7:09

No numerics hacks here; this really computes the volume symbolically. It is a bit tedious and demands some tricks which may appear more obvious in this answer than they would really be on the first try.

Don't expect this code to magically be useful for other such problems; it definitely has its fragile parts.

The idea of the code below is to split the Steinmetz(-like?) solid into easier pieces for volume calculation; specifically involving just one Cylinder and a handful of HalfSpaces each. As we later see, this is actually a necessity to work around orientation problems in Mathematica. Splitting could be visualised with an exploded view drawing:

In order to accomplish this, we have to find an accurate definition for each subregion. Cylinder-cylinder surface intersections on the surface of the solid (parts of ellipses around the origin with various orientations), and points connecting those different parts (points where different pairs of intersections meet) help. These ellipses naturally lie on planes, and define HalfSpaces which can be used to constrain a subregion, and graph cycles of points on the surface of the same Cylinder define which cylinder is used for each of them as a "cap".

CylindricalDecomposition with "Components" is used to guarantee that each individual curve expression is a connected component (there may be several such components on each ellipse), and thus can be used to figure out which points lie on which continuous curve.

(* Unit-length cylinder direction vectors. *)
directions = Append[
Table[RotationTransform[a, {0, 0, 1}][{1, 1, 1}/Sqrt[3]],
{a, 0, 3 Pi/2, Pi/2}],
{0, 0, 1}];

(* Used often. *)
simplify = FullSimplify[#, Element[x | y | z, Reals]] &;

(* For extra warnings if things don't go as planned on comparisons. *)
eqSimplify = (If[! BooleanQ[#],
Echo[#, "Didn't simplify to a Boolean:"]]; #) &@*FullSimplify;

(* Helper function which returns the two ellipses on
surface intersection of two cylinders passing through
the origin, based on their unit-length direction. *)
ClearAll[cylinderRings];
cylinderRings[v1_List, v2_List] :=
With[{transform =
{x, y, z} ->
{Normalize@Cross[v1, v2],
Sqrt[(1 - v1 . v2)/2] Normalize[v1 + v2],
Sqrt[(1 + v1 . v2)/2] Normalize@Cross[Cross[v1, v2], v1 + v2]} .
{x, y, z}]},
ImplicitRegion[simplify[# /. transform], {x, y, z}] & /@
{x^2 + y^2 == 1 && z == 0, x^2 + z^2 == 1 && y == 0}]

(* All expressions for separate dimensional components of
intersections of two cylinders on the surface of the solid. *)
curves =
With[
{solid = RegionIntersection @@ (Cylinder[2 {-#, #}] & /@ directions)},
(* Find topologically separate components of each ring on the surface. *)
(CylindricalDecomposition[
RegionMember[RegionIntersection[solid, #], {x, y, z}],
{x, y, z}, "Components"] & /@
cylinderRings @@ #) & /@
Subsets[directions, {2}] //
simplify // Flatten //
(* Filter out dimensionless solutions. *)
Select[RegionDimension[ImplicitRegion[#, {x, y, z}]] != 0 &]];

(* All points which lie on two separate curves. *)
points = SolveValues[#, {x, y, z}] & /@ Subsets[curves, {2}] //
FullSimplify // Flatten[#, 1] & // DeleteDuplicates;

(* The surface graph; vertices are points above, edges the curves. *)
graph =
(* Create a list of points on each curve. *)
Select[points,
eqSimplify@*RegionMember[ImplicitRegion[#, {x, y, z}]]] & /@
curves //
(* Find the shortest geometric ordering of these points as a
line per each curve. This is a bit of a hack... *)
Map[
First@TakeSmallestBy[Permutations[#],
RegionMeasure[Line[#], 1] &, 1] &] //
(* Create graph edges on basis of these line segments. *)
Map[UndirectedEdge @@@ Partition[#, 2, 1] &] // Flatten // Graph


This graph corresponds to the following (observe the matching number of edges connected to different vertices):

With this graph we can find out subgraphs of points for each Cylinder:

(* Compute per-cylinder subgraphs on the solid. *)
subgraphs =
Subgraph[graph,
(* Select subgraphs with vertices on each cylinder surface. *)
Select[VertexList[graph],
eqSimplify@*RegionMember[
(* Hack around deficiency in RegionBoundary;
implicit cylinder surface (boundary) regions. *)
RotationTransform[{{0, 0, 1}, #}]@
ImplicitRegion[x^2 + y^2 == 1, {x, y, z}]]]] & /@
directions


Now each subregion can be obtained with (region) intersections of aforementioned primitives. This code has a trick up its sleeve: technically one should care about handedness of cycles which form the definition of HalfSpaces constraining the region, but since these Cylinders are in symmetric regarding the origin it's not really necessary - wrong handedness results a symmetric mirror image of the subregion and has the same volume.

In order to pacify Volume computation we reorient the subregion so that the Cylinder is always oriented towards the $$z$$ axis by rotating HalfSpaces; otherwise Mathematica seems to have trouble succeeding in this task.

(* Compute the solid volume by summing up subregion volumes for
intersection of half-spaces and one cylinder each. *)
Function[{dir, subgraph},
Volume[
RegionIntersection[
(* Important hack: avoid arbitrarily oriented cylinders.
Use cylinder oriented towards the z axis instead
and rotate half-spaces to match. Without this,
Mathematica stumbles and fails to compute exact volumes;
this is probably an interaction between internal
CylindricalDecomposition sub-region result and
Integrate over regions. *)
Append[
(* Half-spaces' normal vectors are computed from
vertices as spanning vectors on each graph edge. *)
RotationTransform[{dir, {0, 0, 1}}]@
HalfSpace[Cross @@ #, {0, 0, 0}] & /@ #,
Cylinder[{{0, 0, -2}, {0, 0, 2}}]]]] & /@
(* The volume computation is performed for each 4-cycle on
their corresponding cylinder surfaces. *)
FindCycle[subgraph, {4}, All]],
{directions, subgraphs}] //
Flatten // Total // FullSimplify


$$-\frac{4}{3} \left(-2-10 \sqrt{2}+\sqrt{6 \left(19+6 \sqrt{2}\right)}\right)$$

or

ResourceFunction["RadicalDenest"][%]


$$-\frac{4}{3} \left(-2-10 \sqrt{2}+6 \sqrt{3}+\sqrt{6}\right)$$

The volume computation takes a while.

This matches numerically acquired results:

N[%, 50]

(* 4.4004547140460115048732334911411402560985863336674 *)


The above volume computation can be sped up in this case (but not necessarily even on rather similar cases) by implementing it as an integration in cylindrical coordinates over the region in the cylinder oriented on the $$z$$ axis:

(* Create rotated intersections of HalfSpaces like before. *)
RotationTransform[{dir, {0, 0, 1}}]@
RegionIntersection[
HalfSpace[Cross @@ #, {0, 0, 0}] & /@ #] & /@
FindCycle[subgraph, {4}, All]],
{directions, subgraphs}] //
Flatten //
ParallelMap[
(* Convert Cartesian coordinates to cylindrical before integration. *)
Activate@IntegrateChangeVariables[
Inactive[Integrate][
(* Integrate over the region inside the z-oriented cylinder. *)
Boole[RegionMember[#, {x, y, z}]],
Element[{x, y, z}, Cylinder[2 {{0, 0, -1}, {0, 0, 1}}, 1]]],
{r, \[Theta], zz}, "Cartesian" -> "Cylindrical"] &, #] & //
Total // FullSimplify


No matter what the coordinate system, I haven't had success on getting integration results from cylinders on their alignment as they are actually present in the solid. If that would be easily achievable all the work with piecing up the solid would be unnecessary, but it would appear to be less trivial an endeavour.

Code for Steinmetz solid visualisations featured above:

Show[
With[{solid =
RegionIntersection @@ (Cylinder[2 {-#, #}] & /@ directions)},
BoundaryDiscretizeRegion[solid, RegionBounds[solid],
MaxCellMeasure -> {1 -> 0.02}]],
Graphics3D[
{Thick, Black,
ScalingTransform[{1.001, 1.001, 1.001}] /@
Flatten@Parallelize[
MeshPrimitives[
DiscretizeRegion[ImplicitRegion[#, {x, y, z}],
Method -> "Semialgebraic"], {1}] & /@
curves],
Sphere[#, 0.03] & /@ points}]]


Note Method -> "Semialgebraic" - which is used to discretize implicit single-dimensional regions (that's lines and other curves) in 3D embedding. It's slow but reliable in this task unlike other methods.

MapThread[
Function[{dir, subgraph},
RegionIntersection[
BooleanRegion[
BooleanCountingFunction[{{0, Length@#}}, Length@#],
HalfSpace[Cross @@ #, {0, 0, 0}] & /@ #],
Cylinder[{-dir, dir}],
HalfSpace[-Total[VertexList@#], {0, 0, 0}]] & /@
FindCycle[subgraph, {4}, All]],
{directions, subgraphs}] // Flatten //
Map[BoundaryDiscretizeRegion[#, 1.3 {{-1, 1}, {-1, 1}, {-1, 1}},
MaxCellMeasure -> {1 -> 0.0075}] &] //
Map[TranslationTransform[
Normalize[RegionCentroid[#]]/2][#] &] // Show


Unlike the volume computing code this visualisation needs to care about handedness of subregion boundaries (clockwise/counter-clockwise). This is handled by BooleanRegion to first include both, and then accepting only the side where surface points actually lie with another half-space intersection. Since discretisation is just numerics the more convoluted implicit regions don't really matter much here.

This "exploded view" is important for visually inspecting that all pieces are actually accounted for in the volume computation. Recognising if complicated equations actually are exactly zero is a complicated problem, and can lead to Selects on basis of RegionMember results to drop matches on this indecisiveness. (This is also the reason for couple eqSimplify@*RegionMember constructs in the main code.)

Some additional comments on CylindricalDecomposition, more commonly know as cylindrical algebraic decomposition:

It has nothing to do with Cylinders in the sense of this question; it's just a coincidental name of a method. In general it's a very useful tool in real algebraic geometry, easing mechanised symbolic computation on semialgebraic sets (which correspond with what people often consider in "constructive solid geometry"), at least in low enough dimensions, polynomial orders and coefficient complexities.

Regarding my commentary on Mathematica failing to compute Volumes when Cylinders are not oriented "nicely", or in this case on a coordinate axis: CylindricalDecomposition is used, also to my knowledge internally, to split regions into subregions which are easier to handle by other functions. This splitting produces sub-region splits on coordinate axis directions. I tried to figure out what Mathematica finds hard in this case, and got hints that it fails to compute some volume integrals with some cylindrical decomposition generated subregion variations. How these regions are split depend on the orientation of the region, although it's only a difference in rotation...

• +1 for your masterpiece. Can you kindly explain what CylindricalDecomposition does? Also "Important hack: avoid arbitrarily oriented cylinders. Use cylinder oriented towards the z axis instead and rotate half-spaces to match. Without this, Mathematica stumbles and fails to compute exact volumes; this is probably an interaction between internal CylindricalDecomposition sub-region result and Integrate over regions" is unclear to me. I have to think about your answer before accepting it. Feb 10 at 8:51
• @user64494 Ah, nothing to do with cylinders per se, the naming is a coincidence. :) CylindricalDecomposition is a generally useful tool when working with semialgebraic sets, that is Boolean expressions consisting of comparisons of polynomials, take a look at the documentation! In this case its primary use is to guarantee that each curve is continuous and separate from others from the same "cylinder ring." Topological arguments such as "Components" that one can nowadays give to CylindricalDecomposition are quite useful! Feb 10 at 8:57
• @user64494 Removing single points or theoretically even empty regions is mostly a convenience; they would definitely be both useless and possibly also problematic on the graph. The graph, in this case, shows points at intersection points of the curves ("corners" in the visual sense of the solid) of the solid as vertices, and cylinder-cylinder intersection curves ("ridges") as edges between them. These edges also correspond with planes intersecting with end vertices of edges and the origin, and can be used to construct these simplified, rotatable subregions for volume computation. Feb 10 at 9:53
• +1 Excellent solution! This should be the accepted answer. Feb 10 at 13:11
• @kirma: My colleague presented only the idea and the exact result. He didn't present any calculations. Feb 12 at 10:22

Surprisingly, the volume of this complicated shape can be computed analytically. Some steps exceed MA capabilities, however, using arbitrary precision (I check up to 200 digits) arithmetics one can arrive at numerically exact answer.

r = CylindricalDecomposition[(2 x/3 + y/3 + z/3)^2 + (2 y/3 + x/3 -
z/3)^2 + (2 z/3 + x/3 - y/3)^2 <= 1 &&
(2 x/3 - y/3 - z/3)^2 + (2 y/3 - x/3 - z/3)^2 + (2 z/3 - x/3 - y/3)^2 <= 1 &&
(2 x/3 + y/3 - z/3)^2 + (2 y/3 + x/3 + z/3)^2 + (2 z/3 - x/3 + y/3)^2 <= 1 &&
(2 x/3 - y/3 + z/3)^2 + (2 y/3 - x/3 + z/3)^2 + (2 z/3 + x/3 + y/3)^2 <= 1 &&
x^2 + y^2 <= 1, {x, y, z}];


The result has 10 terms of which only 3 are full-dimensional.

R7 = NIntegrate[Boole[r[[7]]], {x, 0, 2}, {y, 0, 2}, {z, 0, 2}, WorkingPrecision -> 20]
Out[2]= 0.44089520247881330042

R9 =
NIntegrate[Boole[r[[9]]], {x, 0, 2}, {y, 0, 2}, {z, 0, 2}, WorkingPrecision -> 20]
Out[2]= 0.044043057059126566059

R11 =
NIntegrate[Boole[r[[11]]], {x, 0, 2}, {y, 0, 2}, {z, 0, 2}, WorkingPrecision -> 20]
Out[3]= 0.065118579717663229092


Combining together

(R7 + R9 + R11)*8
(*4.4004547140448247646*)


The final result can also be reduced to 1D integral, which is too complicated to present here and which cannot be computed in closed form, but can be evaluated to arbitrary numerical precision.

4.40045471404601150487323349114114025609858633366736911333235058125968\
8407614606852636185052837712051303716016786758083926093099371567866549\
2606639295944504935130060271113944852033606388454773344568577


This numerically exact number can be recognized with RootApproximant yielding the following simple form:

Root[2065408 - 76800 #1 - 86976 #1^2 - 864 #1^3 + 81 #1^4 &, 3]


or

$$\frac{4}{3} \left(10 \sqrt{2}-\sqrt{6 \left(6 \sqrt{2}+19\right)}+2\right).$$

EDIT It was asked how I integrated numerically. Here is the answer.

r11 = FullSimplify[r[[11]]];
a11 = Integrate[Boole[r11], {z, 0, 2}];
b11 = FullSimplify[a11 && 0 < y < 1];
c11 = Integrate[b11, {y, 0, 1}];

r7 = FullSimplify[r[[7]]];
a7 = Integrate[Boole[r7], {z, 0, 2}];
b7 = FullSimplify[a7 && 0 < y < 1];
c7 = Integrate[b7, {y, 0, 1}];

r9 = FullSimplify[r[[9]]];
a9 = Integrate[Boole[r9], {z, 0, 2}];
b9 = FullSimplify[a9 && 0 < y < 1];
c9 = Integrate[b9, {y, 0, 1}];


After these preparations, we can finally do the last integral to arbitrary precision numerically:

vol = NIntegrate[8 (c7 + c9 + c11), {x, 0, 1}, WorkingPrecision -> 200]
RootApproximant[vol] // FullSimplify

• Using ResourceFunction["RadicalDenest"] the expression can be simplified to (-(4/3))*(-2 - 10*Sqrt[2] + 6*Sqrt[3] + Sqrt[6]) Feb 9 at 22:42
• +1. Thank you. You wrote "The final result can also be reduced to 1D integral...". Can you kindly present that integral in your answer? TIA. What does CylindrigalDecomposition do? Feb 10 at 8:47
• Heh, clearly I didn't pay much attention as I only now noticed that you have used CylindricalDecomposition in association with numeric integration... yes, it comes handy on these sort of problems in various forms. :) Feb 10 at 17:17
• I have strong doubts about " but can be evaluated to arbitrary numerical precision. 4.40045471404601150487323349114114025609858633366736911333235058125968\ 8407614606852636185052837712051303716016786758083926093099371567866549\ 2606639295944504935130060271113944852033606388454773344568577". Feb 11 at 19:31
• @yarchik: I changed my mind after your edit. Feb 12 at 9:44

Not an answer, but too long for a comment.

I don't know exactly about the details of the meshing algorithm that Mathematica uses. But I expect that it generates a region whose boundary vertices all lie on the true boundary surface (up to a very small error due to solving the boundary equations iteratively). As the body at hand is convex (it is the intersection of convex cylinders), and as all triangles/tetrahedra have their vertices inside the body (or on its boundary), this would imply that the discretized body lies completely in the inside of the true body. Hence Mathematica's result for the volume should be lower than the true volume (plus some rounding errors). Because 4.40045 is greater than 22/5 (and since I deem 0.00045 as too much for a rounding error), I would say that the 22/5 conjecture is not true.

Edit

When you use DiscretizeRegion and BoundaryDiscretizeRegion one should prescibe a maximal edge length of the generated tets and triangles. For example like this:

M = BoundaryDiscretizeRegion[reg, MaxCellMeasure -> (1 -> 0.01)];
Volume[M]


4.40038

If plot M and zoom-in, then you will see that the sharp edges of the body are cut off. Thus, the actual volume should be larger than 4.40038.

Edit 2

Here is a function that can be used to check whether a vertex lies within the body. It returns 1 if the 3-vector X lies inside and 0 otherwise.

cf = Compile[{{X, _Real, 1}},
Block[{x, y, z},
x = CompileGetElement[X, 1];
y = CompileGetElement[X, 2];
z = CompileGetElement[X, 3];

Boole[
(0.3333333333333333 x + 0.6666666666666666 y -
0.3333333333333333 z)^2 + (0.6666666666666666 x +
0.3333333333333333 y +
0.3333333333333333 z)^2 + (0.3333333333333333 x -
0.3333333333333333 y + 0.6666666666666666 z)^2 <=
1. && (0.6666666666666666 x - 0.3333333333333333 y -
0.3333333333333333 z)^2 + (-0.3333333333333333 x +
0.6666666666666666 y -
0.3333333333333333 z)^2 + (-0.3333333333333333 x -
0.3333333333333333 y + 0.6666666666666666 z)^2 <=
1. && (0.6666666666666666 x + 0.3333333333333333 y -
0.3333333333333333 z)^2 + (0.3333333333333333 x +
0.6666666666666666 y +
0.3333333333333333 z)^2 + (-0.3333333333333333 x +
0.3333333333333333 y + 0.6666666666666666 z)^2 <=
1. && (0.6666666666666666 x - 0.3333333333333333 y +
0.3333333333333333 z)^2 + (-0.3333333333333333 x +
0.6666666666666666 y +
0.3333333333333333 z)^2 + (0.3333333333333333 x +
0.3333333333333333 y + 0.6666666666666666 z)^2 <= 1. &&
x^2 + y^2 <= 1.
]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
];


If I apply this to the mesh M from above then I get

b = cf[MeshCoordinates[M]];
Count[b, 0]
Count[b, 1]


109

399393

So the vast majority of points lies inside. But few of them loe outside. I guess they deviate so little that this would not matter. Maybe one can move them just a little inside? For simplicity, I just move all the points of the discrete region a little bit towards the origin:

Mscaled = TransformedRegion[M, {x, y, z} |-> 0.99999999 {x, y, z}];
Count[cf[MeshCoordinates[Mscaled]], 0]
Volume[Mscaled]


0

4.40038

So now all of the discrete regions vertices lies within the convex body, but the volume enclosed by the mesh is still greater than 22/5.

Edit 3

By the way, in order to have at least on higher resolution picture of the body:

• Thank you for your work. Unfortunately, DiscretizeRegion[reg, MaxCellMeasure -> Infinity] produces "DiscretizeRegion::drf: DiscretizeRegion was unable to discretize the region ImplicitRegion[<<2>>]." and the result of DiscretizeRegion[reg] is not good. There is a possibility Mathematica uses numeric integration to calculate Volume[reg]. Feb 9 at 8:45
• Volume[DiscretizeRegion[reg]] results in 0.323708. Feb 9 at 8:52
• Thank you for your edit. I can't reproduce your M on my weak comp. I can execute M1 = BoundaryDiscretizeRegion[reg, MaxCellMeasure -> (1 -> 0.2)]. Then P = Region[M1, PlotRange -> {{-3/2, 3/2}, {-3/2, 3/2}, {-3/2, 3/2}}];. Feb 9 at 10:18
• Show[{P, RegionPlot3D[(2 x/3 + y/3 + z/3)^2 + (2 y/3 + x/3 - z/3)^2 + (2 z/3 + x/3 - y/3)^2 <= 1 && (2 x/3 - y/3 - z/3)^2 + (2 y/3 - x/3 - z/3)^2 + (2 z/3 - x/3 - y/3)^2 <= 1 && (2 x/3 + y/3 - z/3)^2 + (2 y/3 + x/3 + z/3)^2 + (2 z/3 - x/3 + y/3)^2 <= 1 && (2 x/3 - y/3 + z/3)^2 + (2 y/3 - x/3 + z/3)^2 + (2 z/3 + x/3 + y/3)^2 <= 1 && x^2 + y^2 <= 1, {x, -3/2, 3/2}, {y, -3/2, 3/2}, {z, -3/2, 3/2}, PlotPoints -> 50]}, Opacity -> 0.5] demonstrates that some pieces of M1 are outside of reg. Feb 9 at 10:20
• +1. Thank you. Your answer is solid, as usually. Let us wait some time for an exact answer. Feb 9 at 16:55
Clear["Global*"];

reg = ImplicitRegion[(2 x/3 + y/3 + z/3)^2 + (2 y/3 + x/3 -
z/3)^2 + (2 z/3 + x/3 - y/3)^2 <=
1 && (2 x/3 - y/3 - z/3)^2 + (2 y/3 - x/3 - z/3)^2 + (2 z/3 -
x/3 - y/3)^2 <=
1 && (2 x/3 + y/3 - z/3)^2 + (2 y/3 + x/3 + z/3)^2 + (2 z/3 -
x/3 + y/3)^2 <=
1 && (2 x/3 - y/3 + z/3)^2 + (2 y/3 - x/3 + z/3)^2 + (2 z/3 +
x/3 + y/3)^2 <= 1 && x^2 + y^2 <= 1, {x, y, z}];

tab = Table[{wp, vol = Volume[reg, WorkingPrecision -> wp], 22/5 - vol},
{wp, Prepend[Range[10, 30, 5], MachinePrecision]}];

Prepend[tab, {"Prec", "volume", StringForm[" - volume", 22/5]}] /.
MachinePrecision -> "MP" //
Grid[#, Alignment -> {{Center, {"."}}, Automatic}, Frame -> All] &


It appears unlikely that the exact value is 22/5

• +1. Thank you. I executed Table[{wp, vol = Volume[reg, WorkingPrecision -> wp], 22/5 - vol}, {wp, 10, 20, 5}]`. Your answer looks convincing. Let us wait some time for an exact answer. Feb 9 at 16:52
• Dividing your answers by π didn't seem to make it rational either. Feb 10 at 13:22