Here are some preliminary functions. If you just want to read the document, then you can ignore this first cell.

```
[1]:
```

```
#preliminary functions
def hyperplane_arrangement_to_fan(hyperplane_arrangement):
'''
Input: Linear hyperplane arrangement
Output: Maximal cones of the induced fan. Each maximal cone is the closure of a connected component
of R^d after removal of the hyperplane arrangement
'''
n = len(hyperplane_arrangement)
p = list(powerset(range(n)))
d = hyperplane_arrangement[0].ambient_dim()
cones = []
for entry in p:
signs = [1 if i in entry else -1 for i in range(n)]
inequalities = []
for i in range(n):
hyperplane = hyperplane_arrangement[i]
A = hyperplane.Hrepresentation()[0].A()
ieq = [0]+[signs[i]*a for a in A]
inequalities.append(ieq)
C = Polyhedron(ieqs = inequalities)
if C.dimension() == d:
cones.append(C)
return cones
def intersecting_edges(polytope, u):
'''
Input: A polytope and a normal vector u
Output: Set of edges of polytope that has nonempty intersection with the hyperplane orthogonal to u
'''
u_perp = Polyhedron(eqns= [ [0]+list(u) ])
intersecting_edges = []
for edge in polytope.faces(1):
if edge.as_polyhedron().intersection(u_perp).dimension()>-1:
intersecting_edges.append(edge)
return intersecting_edges
```

# Computing the intersection body of the cube.ΒΆ

We consider the \(3\)-dimensional cube \([-1,1]^3\) and compute the radial function of its intersection body step by step.

```
[2]:
```

```
P = polytopes.cube()
```

The combinatorial type of the intersection of \(P\) with a hyperplane through the origin is fixed for each region of the normal fan of the following zonotope Z. You can visualize this fan by looking at the dual polytope.

```
[3]:
```

```
Z = sum([Polyhedron(vertices = [-vector(v),vector(v)]) for v in P.vertices()])
dualZ = Z.polar()
dualZ.show(frame=False, point = False)
```

We now construct the associated hyperplane arrangement \(H\).

```
[4]:
```

```
H = list(set([Polyhedron(eqns = [[0]+list(vertex)]) for vertex in P.vertices()]))
Fan = hyperplane_arrangement_to_fan(H)
```

```
[5]:
```

```
C = Fan[1]
sample_vector = sum([vector(ray) for ray in C.rays()])
sample_vector
```

```
[5]:
```

```
(-2, 2, 2)
```

We now collect all facets of \(P\) that intersect the hyperplane `perp`

orthogonal to `sample_vector`

. Note that in this case \(Q = P\cap\)`perp`

is a hexagon: indeed we are intersecting \(6\) facets of \(P\).

```
[6]:
```

```
perp = Polyhedron(eqns= [ [0]+list(sample_vector) ])
intersecting_facets = []
if P.intersection(perp).dimension()==2:
for facet in P.facets():
if facet.as_polyhedron().intersection(perp).dimension() == 1:
intersecting_facets.append(facet)
intersecting_facets
```

```
[6]:
```

```
[A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices]
```

The vertices of such a hexagon arise as the intersection of `perp`

with some edges of the above facets of \(P\). We now compute these edges.

```
[7]:
```

```
intersecting_edges_of_facets = []
for facet in intersecting_facets:
edges = intersecting_edges(facet.as_polyhedron(),sample_vector)
intersecting_edges_of_facets.append(edges)
intersecting_edges_of_facets
```

```
[7]:
```

```
[[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices]]
```

Let \(x \in C\) arbitrary. By the theory we know that \(x^\perp\) intersects \(P\) in the set of edges in `intersecting_edges_of_facets`

. Therefore, for each such edge \(\hbox{conv}(a,b)\) we write the corresponding vertex \(v\) of \(P\cap x^\perp\) as:

```
[8]:
```

```
x = vector(var('x1,x2,x3'))
x1, x2, x3 = x
v = []
for edges_of_facet in intersecting_edges_of_facets:
intersecting_points_of_facets = []
for edge in edges_of_facet:
[a,b] = list(edge.vertices())
a = vector(a)
b = vector(b)
intersection_point = ((b.inner_product(x)*a - a.inner_product(x)*b )/ (b-a).inner_product(x))
intersecting_points_of_facets.append(intersection_point)
v.append(intersecting_points_of_facets)
v
```

```
[8]:
```

```
[[((x2 - x3)/x1, -1, 1), (-1, -1, (x1 + x2)/x3)],
[(-(x2 - x3)/x1, 1, -1), (-1, (x1 + x3)/x2, -1)],
[(-1, (x1 + x3)/x2, -1), (-1, -1, (x1 + x2)/x3)],
[(1, -(x1 + x3)/x2, 1), ((x2 - x3)/x1, -1, 1)],
[(1, 1, -(x1 + x2)/x3), (-(x2 - x3)/x1, 1, -1)],
[(1, -(x1 + x3)/x2, 1), (1, 1, -(x1 + x2)/x3)]]
```

We are now ready to compute the volume of the hexagon. This will be done by subdividing the polygon into simplices, and by computing the volume of each simplex as a determinant.

```
[9]:
```

```
determinants = []
for simplex in v:
M = matrix(
simplex + [x]
)
det_M = M.determinant()
if det:
sgn = sign(det_M(x1 = sample_vector[0], x2 = sample_vector[1], x3 = sample_vector[2]))
determinants.append([sgn,M.determinant()])
determinants
```

```
[9]:
```

```
[[-1, -x1*((x1 + x2)/x3 - 1) - x2 - (x2 - x3)*((x1 + x2)*x2/x3 + x3)/x1 - x3],
[1, x1*((x1 + x3)/x2 - 1) - (x2 + (x1 + x3)*x3/x2)*(x2 - x3)/x1 + x2 + x3],
[1,
x1*((x1 + x2)*(x1 + x3)/(x2*x3) - 1) + x2 + (x1 + x2)*x2/x3 + (x1 + x3)*x3/x2 + x3],
[-1, -x1*((x1 + x3)/x2 - 1) + (x2 + (x1 + x3)*x3/x2)*(x2 - x3)/x1 - x2 - x3],
[1, x1*((x1 + x2)/x3 - 1) + x2 + (x2 - x3)*((x1 + x2)*x2/x3 + x3)/x1 + x3],
[1,
x1*((x1 + x2)*(x1 + x3)/(x2*x3) - 1) + x2 + (x1 + x2)*x2/x3 + (x1 + x3)*x3/x2 + x3]]
```

Finally, by adding all the determinants with a correction factor, we obtain the radial function of \(IP\) restricted to the region \(C\).

```
[10]:
```

```
correcting_factor = 1/(factorial(P.dimension()-1) * sum([x_i^2 for x_i in x]))
rho = correcting_factor*sum([sng*det for [sng,det] in determinants])
rho.factor()
```

```
[10]:
```

```
(x1^2 + 2*x1*x2 + x2^2 + 2*x1*x3 - 2*x2*x3 + x3^2)/(x1*x2*x3)
```

From the radial function we can compute the algebraic boundary of \(IP\). In \(C\) this is the zero set of the following polynomial.

```
[11]:
```

```
p, q = rho.numerator_denominator()
alg_boundary = q-p
alg_boundary
```

```
[11]:
```

```
x1*x2*x3 - x1^2 - 2*x1*(x2 + x3) - x2^2 + 2*x2*x3 - x3^2
```