The following answer is only a workaround in a special case, *not* a general solution.
It seems that the implemented method `volume`

is buggy with a good rate, the bug is related
to the method `triangulate`

of a polyhedron!

The SHORT answer:

Use the following function, which splits the polyhedron into parts,
then computes their volume, the bug shows when trying to triangulate!
This is only important message of this answer!

```
def triangulated_volume( P ):
"""compute the triangulation of P, compute the volume as a sum of the pieces.
"""
d = P.dimension()
vertices = [ list(v) for v in P.vertices() ]
return sum( [ abs( matrix( d+1, d+1, [ [ 1, ] + vertices[ t ] for t in T ] ).det() )
for T in P.triangulate() ] ) / factorial(d)
```

Alternatively, one can try to split into pyramids based in origin, which is an interior
point in our case. But this fails relatively often!

```
def my_volume( P, interior_point = vector( [0,0,0] ) ):
"""compute the faces of P, use the interior point to split it
into faces pyramids, compute the volume of each pyramid, sum them up.
"""
return sum( [ Polyhedron( [ interior_point, ]
+ [ vector( v )
for face in P.faces(2)
for v in face.vertices() ] ).volume() ] )
```

Finally, starting with a vector `c`

with components $(x,y,z)$, $0 \leq x \leq y\leq z$, for example `( 0.1, 1.0, 1.1 )`

, and associating the polyhedron for it as in the posted question, then the volume is given by the formula:

```
4*x^2*y - 8/3*y^3 - 8*x^2*z + 4*y^2*z + 4*y*z^2 + 4/3*z^3
```

Of course, this is always working! That is all.

LOOONG answer:

This goes deeper (and also far away), but it may be relevant.

(And there was a lot of geometric fun while studying the object. Thanks for the question!)

(Moreover, in m.h.o. it is important to support computer algebra system with illustrative examples.
We are in a century of all possibilities for the software, by the "good", inspiring examples are missing.)

The plan of the longer answer is as follows.

- Confirm the bug and make a statistic of errors for some randomly generated data.
- Compute the volume starting from the vector
`vector( QQ, [ 13/17, 14/17, 19/17 ] )`

as in the posted question. No error. - Working over
`RR`

, or with a random example runs into a ZeroDivisionError. - Define an ad-hoc version of the volume, that works in our case. The one in the short answer. This ad-hoc function
`my_volume`

has an implementation based on the method `volume()`

of a special class of polyhedra, pyramids, where it seems to work well, and the fact we already know an interior point for our polyhedra. - Check that the same result is obtained for
`volume`

and `my_volume`

and `triangulated_volume`

over `QQ`

. - Find the formula for the volume of the specific polyhedra associated to the group $G$ and a starting vector.
This is a mathematical issue, but visualization and computations are supported by
`sage`

,
and sage is a pretty good mean to put the hands on
the shape of the polyhedron. (Without it, the trial and error cycles would have been too many in my case.) - Check that
`my_volume`

(over `QQ`

) and `triangulated_volume`

always gets this formula.

ONE: We start by confirming the error:

```
import random
import traceback
A = diagonal_matrix( 3, [-1,-1, 1] )
B = diagonal_matrix( 3, [ 1,-1,-1] )
s, t = SymmetricGroup(3).gens()
GROUP = MatrixGroup( [ A, B, s.matrix(), t.matrix() ] )
print "The considered group has order %s" % GROUP.order()
G = GROUP.list()
G = list( G ) # since the list method above gives me a tuple...
G.sort()
def vol( c ):
return Polyhedron( [ g * c for g in G ] ).volume()
vectors = [ vector( QQ, [ randint( 5, 15 ) for _ in range(3) ] )
for trial in range( 10 ) ]
vectors += [ vector( RR, [ random.uniform( 5, 15 ) for _ in range(3) ] )
for trial in range( 10 ) ]
for v in vectors:
try:
volv = vol(v)
print "%s has associated volume %s" % ( v, vol( v ) )
except ZeroDivisionError:
# traceback.print_exc()
print "%s [got ZeroDivisionError]" % v
```

We let this run and obtain this time:

```
The considered group has order 24
(5, 6, 14) has associated volume 22808/3
(12, 15, 15) has associated volume 13860
(5, 6, 9) has associated volume 2436
(8, 8, 6) has associated volume 6784/3
(13, 9, 12) has associated volume 28156/3
(7, 9, 7) has associated volume 5800/3
(7, 9, 9) has associated volume 3096
(12, 9, 10) has associated volume 16984/3
(7, 7, 7) has associated volume 2744/3
(14, 8, 5) has associated volume 30448/3
(12.9589728317627, 5.64395666109784, 9.16678680046540) [got ZeroDivisionError]
(7.01088675869899, 9.62861840238589, 7.41339518763128) [got ZeroDivisionError]
(7.71732224544219, 9.83842107910933, 13.1338440570290) [got ZeroDivisionError]
(13.1463788341743, 13.4704614482142, 11.4765415343236) [got ZeroDivisionError]
(5.68466999612573, 14.0075400171660, 14.7885047686280) [got ZeroDivisionError]
(9.13414931152926, 11.6321442939304, 10.2356341843752) [got ZeroDivisionError]
(8.20989707720079, 8.72810914926507, 9.28924497131092) has associated volume 6534.91775036
(10.3005535665733, 11.1392338434077, 8.80644198437511) [got ZeroDivisionError]
(9.92815378030292, 7.91391408097008, 12.3418545472610) [got ZeroDivisionError]
(5.66242280708254, 13.8289596505955, 12.1218231049418) has associated volume 44714.4988912
```

So we recover the error in most randomly generated cases.

The "group" `G`

(a list) will be used in the sequel.

TWO: Recall that for the `vector( QQ, [ 13/17, 14/17, 19/17 ] )`

we can compute the volume without error:

```
c = vector( QQ, [ 13/17, 14/17, 19/17 ] )
print Polyhedron( [ g * c for g in G ] ).volume()
```

This gives:

```
20716/4913
```

THREE: The ZeroDivisionError occurs over `RR`

...

```
c = vector( RR, [ 13/17., 14/17., 19/17. ] )
print "c =", c
try:
P = Polyhedron( [ g * c for g in G ] )
vol = P.volume()
print "The volume is ", vol
except ZeroDivisionError:
print "... sorry, ZeroDivisionError"
```

This gives:

```
c = (0.764705882352941, 0.823529411764706, 1.11764705882353)
... sorry, ZeroDivisionError
```

And further:

```
try:
P.triangulate()
except ZeroDivisionError:
print "... sorry, ZeroDivisionError during P.triangulate()"
```

and the expected echo is comming:

```
... sorry, ZeroDivisionError during P.triangulate()
```

FOUR & FIVE: Let us define our own functions `my_volume`

and `triangulated_volume`

as in the short answer,
based on the fact that
the `volume`

method seems to work better for pyramids, and always -- avoiding it via `det`

for the computed triangulation.
Then we can compare in a random case:

```
c = vector( [ 12034702141, 13274012341, 230946121123212341 ] )
P = Polyhedron( [ g * c for g in G ] )
print "The volume of P = %s" % P.volume()
print "My volume of P = %s" % my_volume(P)
print "Triangulated... = %s" % triangulated_volume(P)
```

This gives:

```
The volume of P = 16423693391061491154712244902262907888642459299750456
My volume of P = 16423693391061491154712244902262907888642459299750456
Triangulated... = 16423693391061491154712244902262907888642459299750456
```

SIX: Let us start with a general vector `c`

of the form $(x,y,z)$.
Without restricting the generality, we may and do assume $0 < x\leq y \leq z$ . Then the faces of the associated polyhedron are...

eight hexagons, four plus four, and six rectangles.

One hexagon of the first four group has the vertices $\sigma(x,y,z)$,
$\sigma$ in `S(3)`

being a permutation
that acts on the components. The six vertices live in the $(X,Y,Z)$ plane $X+Y+Z=x+y+z$,
the hexagon
$$(x,y,z)\to(x,z,y)\to(y,z,x)\to(z,y,x)\to(z,x,y)\to(y,x,z)\to(x,y,z)$$
(only $(yz)$ or $(xy)$ appear as permutation for each side)

is centered in
$$G=(g,g,g)\ , \ g=\frac 13(x+y+z)\ ,$$ $OG$ is perpendicular to the plane and of length
$ g\sqrt 3 $,
the sides have length either $a=(z-y)\sqrt2$ or $b=(y-x)\sqrt2$, they alternate $a,b,a,b,a,b$,
and we can draw the full lines supporting the segments $a,a,a $ inside to get an
equilateral triangle with the side $a+2b$. The area of the hexagon is thus
$$\frac{\sqrt3}4(\ (a+2b)^2-3b^2\ ) = \frac{\sqrt3}4(\ a^2+4ab+b^2\ ) = \frac{\sqrt3}2(\ (z-y)^2+4(z-y)(y-x)+(y-x)^2\ ) \ .$$ We know $OG$, and thus the volume:
$$ \frac 16(x+y+z)\left[\ (z-y)^2+4(z-y)(y-x)+(y-x)^2\ \right] \ .$$

The eight planes have equations $X+Y+Z=\pm x\pm y\pm z$, and four out of the eight hexagons are of the
same shape as computed above.

One of the other four hexagons has the vertices
$$(-x,y,-z)\to(-x,z,-y)\to(y,z,x)\to(z,y,x)\to(z,-x,-y)\to(y,-x,-z)\dots $$
(twisted permutations),
and is centered in
$$H=(h,h,-h)\ , \ h=\frac 13(y+z-x)\ ,$$ $OH$ is perpendicular to the plane and of length
$ h\sqrt 3 $,
the sides have length either $a=(z-y)\sqrt2$ or $b=(y+x)\sqrt2$, they alternate $a,b,a,b,a,b$,
and we can draw the full lines supporting the segments $a,a,a $ inside to get an
equilateral triangle with the side $a+2b$. The area of the hexagon is thus
$$\frac{\sqrt3}4(\ (a+2b)^2-3b^2\ ) = \frac{\sqrt3}4(\ a^2+4ab+b^2\ ) = \frac{\sqrt3}2(\ (z-y)^2+4(z-y)(y+x)+(y+x)^2\ ) \ .$$
We know $OH$, and thus the volume:
$$ \frac 16(-x+y+z)\left[\ (z-y)^2+4(z-y)(y+x)+(y+x)^2\ \right] \ .$$

The six rectangles are of the same shape now,
one of them has vertices ${z}\times{\ (x,y),(y,x),(-x,-y),(-y,-x)\ }$,
the other being obtained by (the action of $G$, so by) changing the component where $z$ is placed in,
and some signs.
The corresponding six pyramids of the same shape have each height $z$,
and base area $2(y^2-x^2)$. Each volume is thus $$\frac 23z(y+x)(y-x)\ .$$

This gives a formula for the total volume.

For this we define:

```
def V1(x,y,z): return 2/3*z*(y-x)*(y+x)
def V2(x,y,z): return 1/6*( x+y+z)*((z-y)^2+4*(z-y)*(y-x)+(y-x)^2)
def V3(x,y,z): return 1/6*(-x+y+z)*((z-y)^2+4*(z-y)*(y+x)+(y+x)^2)
def V (x,y,z): return 6*V1(x,y,z) + 4*V2(x,y,z) + 4*V3(x,y,z)
def Vinfo(x,y,z): print V1(x,y,z), V2(x,y,z), V3(x,y,z), V(x,y,z)
```

And the formula is:

```
var( 'x,y,z' );
V(x,y,z).factor()
```

This gives:

```
4*x^2*y - 8/3*y^3 - 8*x^2*z + 4*y^2*z + 4*y*z^2 + 4/3*z^3
```

Let us verify this formula.
It may seem useful to get an explicit example an require from `sage`

all information on all computed volumes.
It is simple to rewrite `my_volume`

to get also the information
during the computations, for instance:

```
def info_my_volume( P, interior_point = vector( [0,0,0] ) ):
"""compute the faces of P, use the interior point to split it
into faces pyramids, compute the volume of each pyramid, sum them up.
"""
vol = 0
for face in P.faces(2):
Q = Polyhedron( [ interior_point, ]
+ [ vector( v )for v in face.vertices() ] )
print "Face: %s %s vol = %s" % ( face
, [ list(v) for v in face.vertices() ]
, Q.volume() )
vol += Q.volume()
return vol
```

Then we have in a concrete example:

```
x,y,z = 1, 6, 9
c = vector( [ x,y,z ] )
P = Polyhedron( [ g*c for g in G ] )
vol = info_my_volume( P )
print "P.volume() =", P.volume()
print "Our volume =", vol
```

And this gives:

```
Face: <0,1,2,3> [[-9, -6, 1], [-9, -1, 6], [-9, 1, -6], [-9, 6, -1]] vol = 210
Face: <0,1,4,5,8,9> [[-9, -6, 1], [-9, -1, 6], [-6, -9, 1], [-6, -1, 9], [-1, -9, 6], [-1, -6, 9]] vol = 752/3
Face: <2,3,6,7,10,11> [[-9, 1, -6], [-9, 6, -1], [-6, 1, -9], [-6, 9, -1], [-1, 6, -9], [-1, 9, -6]] vol = 752/3
Face: <0,2,4,6,12,13> [[-9, -6, 1], [-9, 1, -6], [-6, -9, 1], [-6, 1, -9], [1, -9, -6], [1, -6, -9]] vol = 994/3
Face: <1,3,5,7,14,15> [[-9, -1, 6], [-9, 6, -1], [-6, -1, 9], [-6, 9, -1], [1, 6, 9], [1, 9, 6]] vol = 994/3
Face: <4,8,12,16> [[-6, -9, 1], [-1, -9, 6], [1, -9, -6], [6, -9, -1]] vol = 210
Face: <6,10,13,17> [[-6, 1, -9], [-1, 6, -9], [1, -6, -9], [6, -1, -9]] vol = 210
Face: <5,9,14,18> [[-6, -1, 9], [-1, -6, 9], [1, 6, 9], [6, 1, 9]] vol = 210
Face: <7,11,15,19> [[-6, 9, -1], [-1, 9, -6], [1, 9, 6], [6, 9, 1]] vol = 210
Face: <12,13,16,17,20,21> [[1, -9, -6], [1, -6, -9], [6, -9, -1], [6, -1, -9], [9, -6, -1], [9, -1, -6]] vol = 752/3
Face: <8,9,16,18,20,22> [[-1, -9, 6], [-1, -6, 9], [6, -9, -1], [6, 1, 9], [9, -6, -1], [9, 1, 6]] vol = 994/3
Face: <10,11,17,19,21,23> [[-1, 6, -9], [-1, 9, -6], [6, -1, -9], [6, 9, 1], [9, -1, -6], [9, 6, 1]] vol = 994/3
Face: <14,15,18,19,22,23> [[1, 6, 9], [1, 9, 6], [6, 1, 9], [6, 9, 1], [9, 1, 6], [9, 6, 1]] vol = 752/3
Face: <20,21,22,23> [[9, -6, -1], [9, -1, -6], [9, 1, 6], [9, 6, 1]] vol = 210
P.volume() = 3588
Our volume = 3588
```

The information agrees with:

```
sage: Vinfo(x,y,z)
210 752/3 994/3 3588
```

SEVEN: Even now after performing all computations and after some dozen `P.show()`

calls for different shapes,
i am not quite sure that the splitted computation always works, and the formula above is recovering the volume.
Let us make some 100 tests over `QQ`

.

```
okCounter = 0 # so far
def my_formula( x,y,z ):
return 4*x^2*y - 8/3*y^3 - 8*x^2*z + 4*y^2*z + 4*y*z^2 + 4/3*z^3
NTRIALS = 100
for test in xrange(NTRIALS):
c = [ randint( 1, 10000 ) for _ in range(3) ]
c . sort()
x, y, z = c
c = vector( QQ, c )
P = Polyhedron( [ g*c for g in G ] )
vol1 = P.volume()
vol2 = my_volume( P )
vol3 = triangulated_volume( P )
vol4 = triangulated_volume( P )
# print "[%s] %s %s %s -> %s" % ( test, x, y, z, vol1 )
if len( { vol1, vol2, vol3, vol4 } ) == 1:
okCounter += 1 # ok, all answers coincide!
print "%s trials, okCounter = %s" % ( NTRIALS, okCounter )
```

And i've got this time:

```
100 trials, okCounter = 100
```

For some reason this site won't let me format the first code block

see https://trac.sagemath.org/ticket/18214